Ase Chemist

ASE Simulation Skill (v1.4)

Always do this first

Before any non-trivial task, run:

python scripts/check_env.py

It prints which calculators and analysis tools are actually installed, and ends with a one-line "what you can run right now" summary. Recommend a method that the environment supports — do not ask the user to install xTB if EMT or LJ already covers the question.

If a backend is missing and the user wants it, prefer the conda install on HPC / conda systems:

conda install -c conda-forge ase tblite-python mdanalysis matplotlib

…and pip only when conda isn't available:

pip install ase tblite mdanalysis matplotlib

For MACE foundation-model support (v1.2+), install separately:

pip install mace-torch

tblite ships GFN1-xTB and GFN2-xTB and is the supported successor to the deprecated xtb-python. If check_env.py reports [BROKEN] tblite ... C extension unloadable, the pip wheel is libgfortran-incompatible — switch to conda install -c conda-forge tblite-python. The standalone xtb binary (Grimme group) adds GFN0 and GFN-FF if it's on PATH.

mace-torch provides the MACE-MP-0 and MACE-OFF calculators (element-set routing and the ~1k-atom size cliff are covered in Step-2 rule 6). MACE requires torch; CUDA is strongly recommended (CPU mode is ~10× slower). check_env.py reports CUDA status and a size-cliff warning.

Method selection

Walk these three steps in order. Each rule names what to do and why; if the user's case doesn't fit the "because", the rule probably doesn't apply and you should keep walking.

Step 1 — what task is this?

Step 2 — pick the calculator

Apply the first rule that fits the system, in this order:

1. If the user explicitly named a calculator (xTB, EMT, GFN2, TIP3P, …), use that one. Why: don't second-guess an explicit choice.

2. If the system contains only EMT-supported metals (Al, Cu, Ag, Au, Ni, Pd, Pt, plus H/C/N/O as adsorbates), prefer EMT. Why: it's free, instant, and the user probably wants a quick metallic-system answer.

3. If the system is pure water (H₂O molecules only), the choice depends on the task:

   • Production MD → TIP3P. Why: parameterized for exactly this case. Requires ase.constraints.FixBondLengths (rigid-body model; bare ASE lets the O–H/H–H bonds vibrate and the run blows up). The bundled scripts do not auto-attach it — add it inline. See references/ase_core.md §Water (TIP3P + FixBondLengths).
   • One-off relaxations or quick energy / single-point checks on small water systems → GFN2-xTB. Why: simpler — no constraints to set up — and small water clusters are well within xTB's accuracy range. scripts/optimize.py and scripts/single_point.py work as expected with --calculator xtb.

4. If the system has organic / main-group chemistry (heteroatoms, non-EMT elements, organic functional groups, ionic bonding), use tblite GFN2-xTB. Why: EMT will silently give nonsense for non-metals; GFN2-xTB is the cheapest method that knows real chemistry. Sub-rule: if the user wants production-length explicit-solvent MD (≥ 100 ps in a TIP3P/OPC box) on a single small organic, switch to GAFF2 + AM1-BCC via scripts/parameterize_gaff2.py → scripts/run_amber.py. xTB MD with explicit solvent past ~100 ps is impractical (the box pushes well past 1k atoms once water is added); GAFF2 is the right tool for that task and the v1.3 scripts handle the antechamber → tleap → pmemd pipeline. See references/amber.md for force-field and water-model details.

   ⚠️ Architecture note (v1.3 Amber). Amber is the only engine in the skill that does not run through ASE — parameterize_gaff2.py and run_amber.py shell out to AmberTools and pmemd, and the MD loop runs natively in pmemd. This is a performance choice, not forced: ASE exposes a CPU-only in-process path (ase.calculators.amber.SANDER), but pmemd.cuda is ~10–50× faster on production-sized systems. The trade is under review (four options open). See references/amber.md §1 and PLAN.md §"Phase 3", and surface the carve-out when recommending GAFF2 so the user can decide whether they want it.

5. If the system is a transition-metal complex and GFN2 fails to converge, fall back to GFN1-xTB. Why: GFN1 is more robust on d-block elements at the cost of some accuracy.

6. If the system is past the xTB size cliff (~1k atoms; xTB MD on a system that size is impractical past a few ps), reach for a MACE foundation model. Why: GFN2-xTB MD stops being practical at 1k atoms; MACE foundation models (MACE-OFF for organics, MACE-MP-0 for crystals/materials) deliver roughly DFT-quality energies and forces in that 1k–2k atom range on a 40 GB GPU. Use --calculator mace in optimize.py / run_md.py; routing is automatic by element set. Cross-validation against GFN2-xTB is on by default for MD — every 1 ps the script recomputes E and F on the latest frame through xTB and aborts the run when MAE_F > 100 meV/Å. This is the contract under which MACE is recommended at all; do not turn it off (--no-validate) without a specific reason. Read references/ml_potentials.md for the full method-selection rules and known failure modes (liquid mixtures, OOD geometries).

7. If the system is past the MACE ceiling too (>2k atoms on a 40 GB GPU, >~1k on CPU, or anything past ~50k atoms in v1), say so out loud: "v1.2 caps at MACE-medium on a single GPU; v2.2 is slated to add larger ML potentials (CHGNet, Orb) and v2.3 adds Amber for biomolecular MD beyond GAFF2 small molecules." See references/ase_core.md §Appendix for the full size table.

8. If the user explicitly wants DFT (B3LYP, ωB97X-D, M06-2X, PBE0, post-HF, "publication thermochem", "transition-metal barriers within 1 kcal/mol", or "compute G298") and Gaussian is available — use gaussian_sp.py / gaussian_opt.py / gaussian_freq.py. Why: xTB tops out at ~few-kcal/mol error on relative energies and is unreliable on transition metals; DFT is the right tool. No method/basis defaults — surface a recommendation (ωB97X-D/def2-TZVP for organics; PBE0-D3(BJ)/def2- TZVP for transition metals; see references/gaussian_method_selection.md) and confirm before running. The scripts also require explicit --charge, --multiplicity, --mem, --nproc. Solvent → SMD by default. Thermochem parsing is in-house (_gaussian_log.py), no cclib dependency.

Step 3 — confirm the calculator is installed

Read the [OK] / [MISSING] lines from scripts/check_env.py. If your chosen calculator is [MISSING], ask the user to install it; do not silently substitute a wrong-physics fallback. EMT on an organic is the classic failure mode — it will return numbers that look fine and are meaningless.

Verification & clarification

Don't ask what's already named — and frame what you do ask

The two failure modes to avoid: silently picking wrong physics (EMT on an organic returns plausible nonsense; 2 fs timestep without hydrogen constraints; MACE-OFF on a system with metals; Gaussian DFT without explicit method/basis), and re-asking the user something the prompt or structure file already names.

When the answer is genuinely underdetermined, frame the question with the option you'd pick and the reason — e.g., "GFN2-xTB looks right here because the system has heteroatoms; want me to fall back to GFN1 for d-block robustness instead?" — beats a blank "which method?".

Ask the user to verify before recommending execution

After choosing parameters, restate them in a short block and ask the user to confirm before suggesting they run anything. The minimum to surface:

• Calculator (and GFN level if applicable)
• Optimizer / ensemble
• For MD: temperature, friction, timestep, n_steps
• For optimization: fmax, max_steps
• Output paths that will be written

Keep it tight — a 4-6 line summary, not a paragraph. If the user has already approved the plan, don't re-ask.

Scripts — when to invoke each

All scripts live in scripts/ and are parameterized via argparse. Run with --help to see options.

Default: call the bundled CLI. Write inline only for a specific, named capability the script lacks (e.g., "run_md.py has no --barostat, user asked for NPT") — confirmable via --help. Not justifications: user phrasing ("write me a script" — the CLI one-liner IS a script), readability or pedagogy (point at scripts/<name>.py instead), or tweaks already covered by flags (fmax, timestep, ensemble, calculator, seed). If the gap is a missing flag, prefer adding the flag over a one-off rewrite. Name the gap in one sentence before any inline code — no sentence, no carve-out.

Per-script use:

• scripts/check_env.py — Reports installed backends and a one-line capability summary. Run first on any non-trivial task so you recommend a method the environment actually supports.
• scripts/optimize.py — Geometry optimization with BFGS / FIRE / LBFGS, calculator EMT / LJ / TIP3P / xTB / MACE. Real-gas LJ via --epsilon/--sigma/--rc. MACE via --calculator mace (auto-routed to MACE-OFF for pure organics, MACE-MP-0 otherwise).
• scripts/run_md.py — NVE / NVT-Langevin / NVT-Nose-Hoover MD with EMT / LJ / TIP3P / xTB / MACE. Defaults tuned for organic molecules (1 fs, 300 K, Langevin friction 0.01/fs, log every 100 steps). With --calculator mace, GFN2-xTB cross-validation runs by default (the contract in Step-2 rule 6; flags --validate-every, --abort-mae-f, --no-validate).
• scripts/ml_calculator.py — Helper exposing make_ml_calc(atoms, system_class=, device=, model_size=). Imported by optimize.py and run_md.py when --calculator mace. Run as python scripts/ml_calculator.py --structure mol.xyz to print routing without loading weights.
• scripts/validate_ml_md.py — Post-hoc cross-validation of a saved MACE trajectory against GFN2-xTB. Same MAE_F threshold as run_md.py runtime validation; writes validation.csv. For trajectories produced with --no-validate or to re-validate with a different reference / stride.
• scripts/parameterize_gaff2.py — Drives antechamber -c bcc (AM1-BCC) → parmchk2 (frcmod) → tleap (solvate in TIP3P/OPC, neutralize with Na+/Cl-) for a small organic. Output: .prmtop / .rst7 pair for run_amber.py. Mandatory: --net-charge matches the formal charge — getting it wrong silently shifts every partial charge.
• scripts/run_amber.py — Runs Amber MD on a .prmtop (--prmtop) / .rst7 (--rst, not --rst7) pair. --protocol standard runs min → heat (50 ps NVT, 0→300 K) → density (100 ps NPT) → prod (default 500 ps NPT). Engine auto-picks pmemd.cuda > pmemd > sander (--engine). Outputs NetCDF .nc. v1.3 mdin defaults are GAFF2-tuned (protein/NA prmtops run but may want different cutoffs/restraints). For Amber-deep workflows (REMD, MMPBSA, restart/extend, implicit-GB), use amber-chemist instead.
• scripts/gaussian_sp.py — DFT single-point E/F/dipole via ase.calculators.gaussian.Gaussian (g09 fallback). Mulliken charges and HOMO/LUMO parsed by _gaussian_log.py. Required flags + SMD default per Step-2 rule 8 — no silent defaults.
• scripts/gaussian_opt.py — DFT geometry optimization via GaussianOptimizer (Gaussian's L103, one g16/g09 invocation). --convergence is a string (loose/default/tight/verytight), not a numeric eV/Å. Use tight or verytight if the optimized geometry feeds into a Freq job.
• scripts/gaussian_freq.py — DFT frequency + thermochemistry (vib_freqs / ZPE / enthalpy / Gibbs G) via _gaussian_log.py. Reports imaginary modes (warning + nonzero exit). Freq method/basis must match the optimization — not enforced; surface it.
• scripts/_gaussian_log.py — Helper module: regex parsers for Gaussian .log fields ASE doesn't cover (vib_freqs, thermochem, Mulliken charges, MO eigenvalues). Imported by gaussian_sp.py and gaussian_freq.py. Stdlib-only.
• scripts/single_point.py — Single-point energy plus xTB electronic observables (dipole, Mulliken charges, Wiberg bond orders, HOMO-LUMO raw eigenvalue gap). Tagged key=value output. Optimize first — single-point observables on a strained geometry are nonsense. For binding-energy decomposition, run three times.
• scripts/analyze_traj.py — RMSD, RMSF, energy drift, optional RDF from a trajectory. Saves PNG plots and CSV data alongside the input. These analyses ARE the script's primary purpose — do not write a substitute for any of them inline. Handles edge cases (Kabsch alignment, missing-calculator fallback for energy drift, periodic unwrapping for RDF) that an inline rewrite will get wrong.

Growing the skill: when to offer to bundle new scripts

When inline code looks like recurring work, offer to promote it to a bundled script. Offer only when all hold: the code is substantial (>~30 lines or a parametric workflow), no existing scripts/ entry covers it, and the request reads as recurring ("for each molecule", "every time I get a new structure"). Don't offer for trivial one-shots, already-covered tasks, or exploratory/definitional questions.

If the user says yes: refactor into scripts/<verb>.py (naming like optimize.py / run_md.py) with argparse + a top-of-file docstring saying when to reach for it, match output conventions (banner, tagged [OK] / [INFO] lines, plots/CSVs alongside input, meaningful exit codes), add a one-line SKILL.md §Scripts bullet, and verify with --help + one example. If no, leave it and don't re-ask this session.

References — read these on demand

Each file is short and topic-scoped; read the one whose topic comes up.

• references/ase_core.md — structure I/O, ase.build, optimizers, MD integrators, units, Trajectory format, NEB scaffolding.
• references/xtb.md — tblite install, GFN1 vs GFN2, the standalone xtb binary (GFN0/GFN-FF), xTB observables, limitations.
• references/analysis.md — ASE readers vs MDAnalysis, recipes for the analyze_traj.py analyses, pitfalls.
• references/ml_potentials.md — MACE vs xTB, the cross-validation contract, MACE failure modes, the GPU ceiling, troubleshooting.
• references/amber.md — when GAFF2 wins, the antechamber→parmchk2→tleap→pmemd pipeline, force-field/water choices, engine selection, failure modes. Protein/NA (ff19SB/OL21) deferred to v2.3.
• references/gaussian.md — when Gaussian beats xTB, the no-defaults policy + recommended method/basis, SMD vs PCM, g16/g09, the _gaussian_log.py parser, failure modes. Opt=TS/IRC/NBO/TDDFT/post-HF deferred to v3+.

Smell test — don't fabricate technical semantics. If you are about to write "I think <keyword> defaults to ..." or "the standard value for <flag> is roughly ..." — for a Gaussian route line, an xTB GFN convention, a MACE element-set rule, an Amber mdin keyword, or any other domain-specific knob — stop and check the right reference file (or its upstream manual). Hallucinated semantics is a high-cost, hard-to-detect failure mode because the calculation often runs with the wrong value and produces plausible-looking output.

Defaults and conventions

• Units: ASE uses eV, Å, ASE-time-units. Use ase.units.fs / ase.units.kB rather than raw numbers. Temperature kwarg is temperature_K= (canonical since ASE 3.21.0).
• Timestep: 1 fs is safe for organic molecules with all-atom dynamics. Bump to 2 fs only if you constrain hydrogen bonds (ASE doesn't do RATTLE/SHAKE elegantly, so 1 fs is the safer default).
• Friction (Langevin): 0.01 / fs is a reasonable thermostat coupling for production. Higher values (0.1 / fs) for fast equilibration.
• Nose-Hoover coupling: --tdamp 100 fs is the default characteristic timescale for the deterministic nvt-nose-hoover thermostat in run_md.py.
• Optimization tolerance: fmax=0.05 eV/Å for production geometries; 0.01 eV/Å for vibrational analysis input.
• Trajectory format: prefer .traj (ASE binary, includes calculator results) over .xyz (positions only) when energies/forces matter downstream.
• Random seed: Set seed in MD integrators if reproducibility matters to the user.

Reporting results

When you finish a task, report:

1. The method used (calculator + integrator/optimizer) and why it was chosen given system size, available backends, and accuracy needed.
2. Final numbers (energy, fmax, temperature, etc.) with units.
3. Where outputs were written (trajectory, plots, CSVs).
4. Any caveats (e.g., "GFN2-xTB; transition-metal accuracy is limited", "NVE energy drift was 0.3 meV/atom over 1 ps — reasonable").

What v1 does NOT support

Be honest about scope. Deferrals:

• Biopolymer Amber MD (ff19SB+OPC / OL21) → v2.3. BYO-prmtop runs work with run_amber.py but mdin defaults are GAFF2-tuned — flag the mismatch.
• Gaussian Opt=TS / QST / IRC, anharmonic Freq, NBO/NPA, post-HF (CCSD/MP2/CASSCF), excited states (TDDFT/CIS/EOM-CCSD) → v3+; see references/gaussian_failure_modes.md §"Out of scope". (Method strings pass to Gaussian verbatim, so a post-HF route runs, but gets no method-specific parsing or validation — prefer DFT.)
• ML potentials beyond MACE (CHGNet, Orb-v3, M3GNet, SevenNet) → v2.2+; MACE-MP-0 covers most of the same scope today (see references/ml_potentials.md).
• VASP, Quantum ESPRESSO → no v2 plan; CP2K / FHI-aims bridges may land in v3.
• Free-energy (TI/FEP/MBAR), enhanced sampling (REMD, metadynamics, umbrella), QM/MM, constant-pH.
• RESP charges via Gaussian — AM1-BCC only in v1.3.
• SLURM submission scripts; web GUI / visualization servers.