VARNEB — Variable-Cell NEB

This repository now contains the calculator-agnostic VARNEB (VC-NEB) toolkit in vcneb/. It is meant for crystal phase-transition barriers where the cell changes along the path.

The project direction is a pure-Python OpenVCNEB toolkit: no MATLAB or USPEX runtime is required. VASP, ABACUS, and future calculators are external backends behind an ASE-compatible energy/force/stress contract.

ASE status

ASE's built-in ase.mep.NEB does not optimize variable cells in periodic directions. Its interpolation helper can interpolate cells, and ASE's cell filters show the stress-to-generalized-force pattern, but the NEB object itself rejects periodic images with different cells. This prototype fills that gap by optimizing an extended coordinate vector:

• atomic block: fractional coordinates mapped through the reference cell;
• cell block: deformation gradient F, where cell = cell0 @ F.T;
• force block: Cartesian forces and stress transformed into those coordinates.

The implementation works with normal ASE optimizers (FIRE, BFGS, LBFGS) and with any ASE calculator that provides energy, forces, and stress.

The package is installable without MATLAB or USPEX:

python -m pip install .
varneb --version
# ``vcneb`` remains a compatible legacy console alias.

Before an optimization, run_vcneb() checks every image calculator for this contract. Missing stress is a hard error because a variable-cell calculation must not silently replace the cell force by zero. The same preflight report is available through inspect_calculator() and validate_image_calculators().

After a run, chain.saddle_diagnostics() reports the highest interior image, its relative enthalpy, residual generalized force, tangent/normal force components, a local finite-difference tangent curvature, and whether an interior peak rises above both endpoints. A negative curvature alone is not enough: a monotonic band can make the highest interior image look like a CI candidate even though no transition-state barrier is present. The curvature and peak flags are path-local diagnostics, not a full Hessian or proof of first-order saddle character.

Files

• vcneb/core.py: VC-NEB algorithm and optimizer-compatible object.
• vcneb/modes.py: mode-guided initial paths and modal path projections.
• vcneb/vasp.py: VASP input parsing and per-image calculator setup.
• vcneb/abacus.py: ABACUS calculator factory adapter.
• vcneb/calculator.py: capability preflight and image-aware calculator diagnostics.
• vcneb/executor.py: optional image-level concurrent calculator executor; the controller remains single-process and each external calculator job step must use isolated directories (and srun --exclusive on Slurm).
• examples/run_toy_vcneb.py: analytic smoke test with a known 0.25 eV barrier.
• examples/run_hfo2_t_po_model_vcneb.py: mapped 12-atom HfO2 T -> PO geometry smoke test with a synthetic endpoint double-well calculator.
• examples/compare_initial_cell_paths.py: calculator-free comparison of linear and logarithmic-strain initial paths for any ASE-readable endpoint pair (defaults to HfO2); supports --mapping auto.
• examples/run_vcneb_vasp.py: VASP driver based on the existing endpoint layout.
• examples/run_fixed_cell_ase_comparison.py: ASE CINEB versus fixed-cell VCNEB comparison.
• examples/run_vasp_single_image_smoke.py: real VASP energy/force/stress smoke driver.
• examples/run_vcneb_abacus.py: ABACUS driver skeleton.
• examples/relax_abacus_native.py: ABACUS-native atomic/cell endpoint relaxation (calculation cell-relax, relax_method bfgs); ASE is used only for structure conversion and final CONTCAR export.
• scripts/audit_native_endpoint.py: read-only endpoint gate for native summaries (return code, convergence flag, composition, force and stress).
• scripts/promote_hfo2_endpoints.py: validates both endpoint summaries and atomically publishes only passing CONTCAR files to the production relaxed_T/ and relaxed_PO/ directories.
• examples/relax_abacus_endpoint.py: ASE-driven endpoint adapter. It is a numerically equivalent fallback to native cell-relax when the latter has step-control trouble; pass --stress-kbar to require a force-and-stress gate.
• scripts/setup_hfo2_t_po_validation.py: builds the HfO2 T -> PO validation fixture from local source structures or portable copies.
• scripts/validate_vcneb_inputs.py: static dry-run validator for VASP/ABACUS VC-NEB image directories.
• scripts/audit_vcneb_result.py: calculator-free audit of a completed summary, including generalized-force, barrier, volume, geometry and interior-barrier gates. Use --max-stress-kbar when the production endpoint/path policy requires a common stress threshold in addition to the generalized-force gate.
• scripts/compare_vcneb_images.py: calculator-free comparison of completed 5/7/9-image summaries, including shared calculator settings and explicit handling of consistent barrierless paths.
• docs/batio3_validation_protocol.md: fixed BTO settings, image-count convergence gates, CI staging and recovery/archive requirements.
• tests/check_vcneb_forces.py: finite-difference checks for force/stress transforms.

Cluster execution policy

For long NEB or DFT runs, use the shared cluster only. The current production entry point is Hefei through ssh hf and Slurm. Before every submission, inspect sinfo, squeue, and the account's own jobs, then select a suitable partition, node allocation, and task count from the actual calculator cost. Do not mechanically request 40 cores: small cells and cheap settings may use fewer tasks, while a large k-point/plane-wave calculation may justify more. Record the partition, node, task count, elapsed time, and input version with each result. Do not run long jobs on the gateway or in local WSL, and do not use qsub/PBS for this project.

Mode-guided paths and component constraints

Before attaching an expensive calculator, an interpolated path can be checked for cell and atom geometry:

from vcneb import interpolate_vcneb, path_geometry_diagnostics

images = interpolate_vcneb(
    initial, final, n_images=9, mic=True,
    minimum_distance=1.0,
    maximum_deformation=0.8,
)
report = path_geometry_diagnostics(images)

The optional thresholds on interpolate_vcneb() raise a ValueError with the image index when a path violates them. The report itself is calculator-free and contains the per-image volume, shortest periodic interatomic distance, and deformation norm.

The cell path can use the historical linear deformation interpolation, a logarithmic strain path for aligned symmetric-positive deformation gradients, or a project-specific callback:

images = interpolate_vcneb(
    initial, final, n_images=9, mic=True,
    cell_interpolation="log_strain",
)

def interpolate_cell(lam, deform0, deform1):
    return (1.0 - lam) * deform0 + lam * deform1

images = interpolate_vcneb(
    initial, final, n_images=9,
    cell_interpolation=interpolate_cell,
)

log_strain removes the relative rigid rotation when align_cells=True and rejects unsuitable non-positive or non-symmetric deformation gradients. A custom callback must return a finite 3x3 deformation matrix; every resulting cell is still checked for a positive determinant.

The endpoint atom order can be kept explicitly or inferred by element and periodic geometry. The inferred permutation is returned by validate_atom_mapping() and should be saved with the run manifest:

from vcneb import validate_atom_mapping

mapping_report = validate_atom_mapping(initial, final, "auto", mic=True)
images = interpolate_vcneb(
    initial, final, n_images=9, mapping=mapping_report["mapping"], mic=True,
)

Automatic mapping is a geometry heuristic, not a chemical identity proof. For large reconstructive transitions, inspect mapping_report or provide an explicit final-atom permutation.

For periodic crystals whose endpoint coordinate origins differ by a lattice translation, enable the joint translation/mapping alignment before any calculator is attached:

images = interpolate_vcneb(
    initial, final, n_images=9, mapping="auto", mic=True,
    align_translation=True, cell_interpolation="log_strain",
    minimum_distance=1.0,
)

This changes only the periodic coordinate gauge of the final endpoint. The translation is chosen together with the element-grouped assignment, because independent MIC choices can otherwise create artificial short bonds in the interior path. Always inspect the resulting mapping and geometry report for reconstructive transitions.

The mode feature is intentionally split into an initial-path generator and a diagnostic projection. The former bends the interior images with an endpoint-zero envelope, then a normal unconstrained NEB calculation can relax to the full-space MEP:

from vcneb import Mode, mode_guided_path, project_path_onto_modes

mode = Mode.from_file("soft_mode.dat", n_atoms=len(initial))
images = mode_guided_path(
    initial,
    final,
    n_images=7,
    mode=mode,
    amplitude=0.20,
    envelope="sin",
    mic=True,
)
modal_coordinates = project_path_onto_modes(images, initial, mode)

mode_guided_path() accepts the same cell_interpolation, mapping, align_translation, minimum_distance, and maximum_deformation controls as interpolate_vcneb(), so a mode-guided path cannot silently bypass endpoint mapping or calculator-free geometry checks.

For an explicit preflight before creating an optimizer:

from vcneb import validate_image_calculators

reports = validate_image_calculators(images)

Each VASP or ABACUS image should have its own calculator directory. A custom launcher can additionally call validate_image_calculators(images, require_directory=True, require_unique_directories=True).

Text mode files contain n_atoms rows of three numbers or a flattened 3N vector. JSON and NPZ files can also carry an optional cell 3x3 deformation mode. A normalized atomic mode has an amplitude in Angstrom; a cell mode is dimensionless in deformation-gradient space.

For ABINIT-like fixed components, pass an (n_atoms, 3) mask. 1 keeps a component active and 0 keeps it fixed:

chain = VCNEB(images, atom_mask=mask)

This is a constrained calculation when used during optimization. A mode-guided initial path by itself is not constrained and should be preferred when the final result is intended to be a full MEP.

For strict mode-subspace dynamics, first build a basis in the same extended coordinate space as VCNEB:

from vcneb import Mode, build_mode_basis

basis = build_mode_basis(Mode([[1.0, 0.0, 0.0]]), initial)
chain = VCNEB(images, mode_basis=basis, constraint_mode="subspace")

constraint_mode="subspace" requires the endpoint displacement to lie in the supplied basis and projects the initial interior images and every optimizer update into that affine subspace. Use constraint_mode="projected" when the non-mode part of an existing initial path should remain fixed while only the optimization update is projected. build_direction_basis() provides the same interface for one allowed direction per atom or arbitrary atomic direction combinations. Call direction_basis_conflicts() before combining a direction basis with component masks; it reports partially clipped columns, fully inactive columns, and rank loss. These are constrained transition paths; a constrained saddle is not automatically a first-order saddle in the full configuration space.

The three semantics are compared reproducibly by examples/compare_mode_path_variants.py: an unconstrained path, a mode-guided initial path followed by unconstrained VCNEB, and a strict atomic-plus-cell mode subspace. On the coupled analytic toy surface, using seven images and fmax=0.002 eV/A, all three recover the same 0.25 eV barrier within 2.6e-5 eV; the strict case is exactly on the known MEP. The HF record is outputs/mode_path_variants_hf.json. This comparison is an algorithm test: a strict constrained saddle must still be interpreted as a saddle in the constrained space, not automatically as a full-space saddle.

For a perturbed or noisy initial path, use a two-stage CI protocol: first call run_vcneb(..., climb=False) to relax the path, then call it again on the same images with climb=True for saddle refinement. Starting CI immediately can select a wrong image and produce a folded, projection-stationary path. The robustness example tests four deterministic perturbation seeds with both linear and log_strain cell interpolation; the staged protocol recovers the analytic 0.25 eV barrier in all eight cases. Always inspect path_diagnostics() in addition to the optimizer residual: the latter is the projected NEB force and is not, by itself, a proof that the physical path is a valid MEP. The record is outputs/vcneb_robustness_staged_hf.json.

The same protocol is available directly through run_vcneb(..., climb_after=N): the first N completed optimizer steps use ordinary NEB and CI is then enabled automatically. path_geometry_diagnostics() reports adjacent extended-coordinate segment cosines; pass fold_cosine_threshold=0.0 to validate_path_geometry() when a folded path must be rejected. The default remains diagnostic-only so existing workflows are not silently changed.

For a physically interpretable CI result, require saddle_diagnostics()["has_interior_barrier"] and inspect path_diagnostics()["ci_warning"]. If the band is monotonic or its interior peak is below an endpoint, report a barrierless/unresolved path rather than calling the highest interior image a transition state.

Local checks

python tests/check_vcneb_forces.py
python examples/run_toy_vcneb.py
python scripts/setup_hfo2_t_po_validation.py
python examples/run_hfo2_t_po_model_vcneb.py
python examples/run_abacus_single_image_smoke.py
python examples/relax_abacus_endpoint.py --help
python examples/relax_abacus_native.py --help
python examples/run_vasp_single_image_smoke.py --help
python examples/run_fixed_cell_ase_comparison.py
python examples/run_vcneb_convergence.py
python examples/run_vcneb_robustness.py
python examples/run_release_and_refine.py
python examples/run_finite_difference_report.py
python scripts/audit_vcneb_result.py path/to/completed_vcneb_workdir

Expected toy output:

barrier_eV=0.250004
delta_eV=0.000000

The HfO2 fixture validates a real 12-atom variable-cell path:

• initial endpoint: tetragonal P4_2/nmc (137);
• final endpoint: polar orthorhombic Pca2_1 (29);
• model VC-NEB smoke barrier: 0.800023 eV.

These checks were copied to 235:/home/zhuxd/abacus/agent-runs/20260618-vcneb and passed on cu05.

VASP run

From this repository root:

python examples/run_vcneb_vasp.py \
  --initial initial_state/relax \
  --final final_state/relax \
  --workdir run_VCNEB/run \
  --n-images 7 \
  --fmax 0.05 \
  --steps 300 \
  --vasp-bin /home/zhuxd/Software/src/vasp/6.3.2/bin/vasp_std \
  --ncores 8

The script reads CONTCAR/POSCAR endpoints, interpolates fractional coordinates and cell deformation, copies POTCAR, and uses static single-point VASP settings with IBRION=-1, NSW=0, ISIF=2, ISYM=0.

To continue a stopped run from the latest complete chain snapshot:

python examples/run_vcneb_vasp.py \
  --initial initial_state/relax \
  --final final_state/relax \
  --workdir run_VCNEB/run \
  --n-images 7 \
  --resume

The trajectory stores one frame per image at each optimizer step. Resume reads only complete n_images frame groups and ignores an interrupted partial tail. When appending to an existing run, snapshot files under snapshots/ continue from the next available step_#### / chain_step_####.traj index.

Before launching VASP, run a dry static check:

python scripts/validate_vcneb_inputs.py \
  --mode vasp \
  --template-dir initial_state/relax \
  --image-root validation/hfo2_t_to_po

ABACUS run

examples/run_vcneb_abacus.py expects an ASE ABACUS calculator. The material settings can be supplied on the command line, so the same VC-NEB driver can be used with a different calculator configuration:

• ecutwfc
• kpts
• pp
• basis
• pseudo_dir
• basis_dir
• spin, smearing, van der Waals, and convergence settings

For example, the HfO2 Dojo-FR setup on the shared cluster uses:

python examples/run_vcneb_abacus.py \
  --initial validation/hfo2_t_to_po/image_00/POSCAR \
  --final validation/hfo2_t_to_po/image_06/POSCAR \
  --workdir validation/hfo2_t_to_po/abacus_vcneb_smoke \
  --command "srun -n ${SLURM_NTASKS:-8} /home/zhuxd/Software/abacus/INSTALL/3.10.0-LTS/bin/abacus" \
  --pseudo-dir /home/zhuxd/abacus/PSEUDO/ABACUS-orbitals/Dojo-NC-FR/Pseudopotential \
  --basis-dir /home/zhuxd/abacus/PSEUDO/ABACUS-orbitals/Dojo-NC-FR/selected_Orbs \
  --pp Hf=Hf.upf --pp O=O.upf \
  --basis Hf=Hf_gga_7au_100Ry_4s2p2d1f.orb \
  --basis O=O_gga_7au_100Ry_2s2p1d.orb \
  --ecutwfc 60 --kpts 1 1 1

The driver enforces cal_force=1, cal_stress=1, and out_stru=1, which are required for VC-NEB.

ABACUS runs support the same --resume and --resume-trajectory options as the VASP driver. Before a costly run, --validate-only creates the image directories, checks energy/force/stress capability and unique per-image directories, and writes vcneb_preflight.json without launching ABACUS. Completed runs include the calculator report, Slurm metadata and Git revision in an atomically updated vcneb_summary.json; resuming does not overwrite the original initial-vcneb.traj.

For image-level concurrency, pass --image-workers N. The controller remains one Python process and evaluates up to N independent calculators concurrently; on Slurm the calculator command must use srun --exclusive, and the allocation must provide N times the MPI width requested by one calculator. The default is 0 (serial images), so ordinary runs retain the reference execution path. --image-retries K retries only a failed image evaluation (default 0); the calculator must be restart-safe in its per-image directory when this is enabled. Parallel runs append image_worker_manifest.jsonl (or the path supplied by --image-manifest) after each controller evaluation batch, including status, elapsed time and per-image attempt counts. The bundled ABACUS adapter also disables ASE's process-global ase_sort.dat write when the input atoms are already grouped by species (as in the BTO and HfO₂ fixtures). This avoids a thread-level parser race; non-identity atom orders should use the serial image backend or an adapter that provides a directory-local sort file. For reproducible recovery from a known complete snapshot, combine --resume --resume-step N with --resume-trajectory; negative N counts from the end and -1 means the latest complete chain. The Hefei BTO template cluster/hf_batio3_vcneb_parallel.slurm demonstrates four 32-MPI workers in a 128-task allocation.

The HfO₂ production template cluster/hf_hfo2_vcneb_parallel.slurm uses the same controller/worker layout (four isolated image workers × 32 MPI by default), with the 100-Ry and Orb-DZP-10au Hf/O inputs. It remains a template until both native cell-relax endpoints pass their force and stress gates.

For an ABACUS template directory, the dry static check is:

python scripts/validate_vcneb_inputs.py \
  --mode abacus \
  --template-dir path/to/abacus/template \
  --image-root path/to/vcneb/images

Current limitations

• Endpoints must have the same atom count and composition. Identity order is supported for audited structures; mapping="auto" provides an element-grouped periodic geometry heuristic, and explicit permutations are recommended for reconstructive transitions.
• The code removes a global cell rotation during interpolation, but production phase-transition work still needs careful endpoint matching.
• Stress can only drive physical strain components; arbitrary cell rotations are a gauge, not a real force degree of freedom.
• Strict mode-subspace and projected-update constraints are available through mode_basis; direction-basis conflict diagnostics are available through direction_basis_conflicts(). VCNEB.path_diagnostics() reports per-image cell metrics, raw atomic force/stress, cell force, and the true/spring/NEB force decomposition needed to audit a variable-cell path. Generalized sparse projectors and release-then-refine workflows are still under development. examples/run_release_and_refine.py demonstrates the supported two-stage workflow on an analytic coupled potential.
• This is a working prototype, not yet a published SSNEB implementation. Treat DFT results as research data: compare against fixed-cell NEB and endpoint cell-relax results before trusting barriers.