2026-10-03: the batch path for the MBE step (MBE phase 3)#

Part of the MBE campaign: the design is ~/Work/SEAMM/MBE_correction_step_design.rst, and the step is mbe_step. vasp_step gains the Model Chemistry batch contract (get_model_chemistry_options, get_task, analyze_task, can_run_task) with the “fragment registered to a parent cell” mode the MBE step’s periodic low level needs. Reviewed by the “design” session; approved by Paul.

Step 1: the input builder, with no behaviour change#

vasp_step/inputs.py builds INCAR, POTCAR, KPOINTS and POSCAR from explicit settings, and the Energy substep now calls it. Values that need the flowchart (the initial WAVECAR, an ENCUT expression, extra keywords with variables) are resolved by the substep and passed in.

Identity test.

  • tests/data/inputs holds the substep’s four input files for four representative settings, captured from the unchanged code:

    • water at Gamma with stress;

    • Si by k-spacing (odd, Monkhorst-Pack, Methfessel-Paxton, pressure only);

    • spin-polarized LiF (Li_sv, explicit grid, tetrahedron + Blöchl + Fermi-Dirac, extra keywords);

    • a water Optimization.

  • The test requires the files to stay byte-identical. The POTCAR library is fake, since licensed files can’t be shipped.

  • To regenerate them from the code before a change:

    • git worktree add <tmp> <old commit>, then copy the installed _version.py into it;

    • run tests/input_harness.py’s capture() with that worktree first on PYTHONPATH.

    This is how the fixtures were made and re-verified (vasp_step 1690fe2).

Headless fixes.

  • The default potentials moved from tk_energy.py, which imports tkinter, to potentials.py. A calculation set up without the dialog has an empty “potentials” parameter: it now gets the set’s defaults instead of a KeyError.

  • An ENCUT expression of ENMAX takes ENMAX from the potentials when the dialog has not set it.

Step 2: the batch path#

  • Model chemistries: VASP:DFT@<functional>/<PAW|PAW-hard|PAW-LDA>@<ENCUT eV> in the existing grammar.

    • Declared: periodic, not MDI, prefers_batch, stress_convention = "stress".

    • vasp_step already turns VASP’s kB pressure into a GPa stress (−0.1 × kB), so nothing is flipped twice.

  • Grid (grid.py) reproduces the prototype’s gen_registered.py:

    • NG is the smallest even 2·3·5·7-smooth n with h ≤ max_spacing (150 for the 12.4297 Å pilot cell);

    • the box is that n·h with n·h ≥ max(12 Å, extent + padding);

    • the shift is a whole number of grid steps, k = round((L/2 − centre)/h).

    • The tests compare NG, box and coordinates with the prototype’s own POSCARs and INCARs (cell, m00, d00_18) to 1e-9 Å. The outer pair d00_17, at image (0, 0, −1), is checked against a port of register.

    • Only orthorhombic parents are supported.

  • get_task:

    • INCAR, POTCAR, KPOINTS through the same builder as the substep. A test checks that the batch path and the substep write identical INCAR, KPOINTS and POTCAR for the same settings. That test caught string values (“no”) being truthy, so the batch settings now go through EnergyParameters exactly like the substep’s.

    • POSCAR in Cartesian with 10 decimals.

    • For a charged structure, NELECT from the POTCARs’ ZVAL; for an open shell, NUPDOWN.

    • Fragments get the dipole correction (IDIPOL = 4, DIPOL = 0.5 0.5 0.5: a registered fragment is centred to within h/2).

    • return_files: INCAR, KPOINTS, POSCAR, OUTCAR, OSZICAR, vasprun.xml, vasp.out, and the dftd4 files.

    • success_text: OUTCAR “General timing”, and the presence of dftd4.json.

    • estimated_seconds is calibrated on the prototype: 330 s for a water fragment on a 150³ grid at 8 ranks, scaled by grid points, √atoms and ranks.

  • D4:

    • Neither TinkerCliffs VASP build (6.6.1-intel2025b, 6.5.1-intel2023a) is compiled with DFTD4, so IVDW = 13 is out. The dftd4 CLI runs inside the task after VASP: --func <f> --grad --json dftd4.json --charge Q.

    • It runs on POSCAR for a cell (periodic) and on fragment.xyz for a fragment (isolated).

    • The rule and its evidence: on monomer m00, dftd4 on the fragment’s POSCAR (in its 12.43 Å box) gives −1.3399e−3 eV, while the xyz gives −1.3103e−3 eV. The prototype’s fragment D4 (−1.31026e−3, the library without a lattice) matches the xyz exactly. Image dispersion would break the regression and the consistency with the molecular level.

    • TinkerCliffs: conda-forge dftd4 4.3.0 at /projects/seamm/conda-envs/dftd4/bin/dftd4 (installed with Paul’s OK). It reproduces the prototype’s cell D4 (−2.655539152468714 eV, P_D4 −2322.500241 atm) and m00.

  • analyze_task:

    • Reads vasprun.xml, taking the calculation’s own <energy> block: the first e_0_energy inside the last <calculation> is an SCF step’s (+7052.9 eV for the pilot cell, against the real −1007.455).

    • Forces go back to the structure’s atom order; the stress is −0.1 × the kB tensor (vasprun.xml’s stress has OUTCAR’s “in kB” sign); dftd4’s energy, gradient and virial are added (σ_D4 = virial/V).

    • An SCF without “aborting loop because EDIFF is reached” in OUTCAR raises.

    • The vaspout.h5 requirement is the Energy substep’s, untouched; the batch path doesn’t need HDF5.

    • Tested on the prototype’s own pilot runs. Energies match to 1e-6 kJ/mol; forces to 1e-4 kJ/mol/Å, because vasprun.xml has 8 decimals and OUTCAR 6.

    • The stress of the cell is σ = −(P_VASP + P_D4), and gives the prototype’s −61,519.55 − 2,322.50 atm. vasprun.xml’s tensor is asymmetric by 2e-6 GPa.

  • resolver (entry point vasp): puts vasp.ini’s gamma_code (or code) and dftd4 into the command, before the executor fills in {NTASKS}.

Review (design, 2026-10-03)#

No wrong-physics findings. Fixed, with tests:

  1. Restart identity. The tasks’ fingerprints leave out the INCAR’s parallelization keywords (NCORE, KPAR, NPAR, NSIM), so a rerun with another number of ranks reuses the finished calculations. They used to be recomputed.

  2. Small cells. A cell narrower than 10 Å is refused at the Gamma point alone; options["k_spacing"] (1/Å) gives a Gamma-centred mesh and the standard build ({code}).

  3. Tests. The tests skip the level-grammar check without model_chemistry_step.

  4. Pin. seamm-exec>=2026.10.3 is pinned.

  5. Hard potentials. PAW-hard includes P_h, S_h and Cl_h, for PF₆⁻ and other electrolytes.

  6. ENCUT. An ENCUT below the largest ENMAX is refused.

  7. Nits.

    • The resolver falls back to vasp_std when there is no Gamma build.

    • The -D4BJ (IVDW = 13) levels are hidden.

    • LDA potentials go only with the LDA functionals.

    • The estimate for a cell uses its grid.

    • The POTCAR is deleted after a successful run.

Also: a -D4 level strips IVDW and VDW_* from the functional, so D4 is never added twice.

The cost estimate (after C.1)#

C.1’s first cell job hit its 1 h limit: the old estimate (330 s × grid/150³ × √(atoms/6) × 8/ranks, from one fragment) gave 930 s for the 64-water cell, which took 4,822 s.

  • The data.

    • The VASP step’s timing file on TinkerCliffs (~/.seamm.d/timing/vasp.csv, 1.7 GB) holds 434,729 runs. 396,528 of them are Gamma-point single points on 8 ranks with known electrons: r2SCAN(-D3BJ), 1-432 atoms, 2025-12 to 2026-04.

    • The prototype’s 64,656 VASP runs at the MBE settings, from the frames’ vasp.tar.gz.

    The per-run tables are /projects/seamm/psaxe/mbe_c1/vasp_timing_features.csv.gz and mbe_vasp_timings.csv.

  • The fit. log t = −7.209 + 0.636 log Ne + 1.342 log(V (ENCUT/500)^1.5) (s on 8 ranks).

    • R² = 0.67 in log t; 68% of runs within ×1.3, 95% within ×2.

    • Against the prototype on 8 ranks the median real/fit is 1.09 for monomers and 0.86 for pairs, so no settings factor is needed.

    • Rank scaling: the prototype’s 16-rank runs were 1.4-1.6× faster, so speed ∝ ranks^0.5 above 8 and ∝ ranks below.

  • A cell’s time limit. A cell gets Resources.walltime = 3 × the estimate, at least 1 h. The prototype’s slowest cell, a dense frame, took 2.8× the median. A fragment keeps the bundle’s walltime.

  • The timing file. It is large because every row carries the full POSCAR/INCAR/KPOINTS, and one line had NUL bytes (concurrent appends without a lock). Node medians range 0.68-1.40× the fit; tc103 is the slowest.

Noticed, not changed#

  • metadata.py gives plain revPBE IVDW = 12, the same as revPBE-D3BJ: vasp-step#13.

  • From the review, both pre-existing:

    • an ENCUT expression of ENMAX never reaches the substep (Parameters refuses the string): vasp-step#15;

    • initial wavefunction = "random guess" tries to copy a file named “random guess”: vasp-step#16.

  • TinkerCliffs’ /projects/seamm/SEAMM/vasp.ini uses gamma_code = mpiexec -np {NTASKS} vasp_std. vasp_gam would be faster for Gamma-only runs; the prototype used it.