Add a finite-temperature phonon workflow

[hrushikesh-s]

Summary

This PR adds a workflow for effective harmonic phonons at a finite temperature. pheasy fits second-order force constants with LASSO to snapshots of an MD run at the temperature, as in the temperature-dependent effective potential (TDEP) method. The force constants include the anharmonic effects at that temperature.

What the workflow does

1. Relaxes the structure. This can be switched off.
2. Builds a diagonal supercell, at least 12 Å along each lattice vector by default.
3. Optional, off by default: runs an NPT MD to get the cell at the temperature (see below).
4. Runs an NVT MD at the temperature. The defaults are 8 ps, a 1 fs time step and a Nose-Hoover thermostat. A Langevin thermostat is also available.
5. Leaves out the first 1 ps and picks 50 snapshots evenly from the rest.
6. Computes the forces on the snapshots and on the undisplaced supercell.
7. Fits the force constants with pheasy (LASSO). The forces of the undisplaced supercell are subtracted from those of each snapshot.
8. Computes the band structure and DOS. Imaginary modes are reported, not removed.

The MD trajectory is checked for melting, for a move away from the starting structure and for a drift of the potential energy. The result is stored in the output. A warning is raised if the check fails, but the fit is still done. A run is reported as melted when its Lindemann ratio is above 0.15. This is the ratio at melting of fcc solids (Saija et al., J. Chem. Phys. 124, 244504 (2006)).

Makers

• FiniteTemperaturePhononMaker in atomate2.vasp.flows: VASP MD and VASP statics.
• In atomate2.forcefields.flows.finite_temperature_phonons:
  • ForceFieldFiniteTemperaturePhononMaker: force field MD and force field statics.
  • VaspMDMLFFStaticFiniteTemperaturePhononMaker: VASP MD and force field statics.
  • MLFFMDVaspStaticFiniteTemperaturePhononMaker: force field MD and VASP statics.

Optional NPT MD

Without it, the MD runs at the volume of the starting structure, so thermal expansion is missing. With npt_maker set, an NPT MD runs first, by default for 8 ps at 0 kbar. Its cell is averaged after the first 2 ps and symmetrized to the point group of the structure. The atoms keep their fractional coordinates and can be relaxed in the new cell with fixed_cell_relax_maker. The NVT MD and the fit then use this cell. For force fields, from_force_field_name(..., run_npt=True) sets both makers. For VASP, get_npt_maker gives an NPT MD maker (MDALGO = 3, ISIF = 3, PSTRESS from pressure). No VASP NPT maker is set by default.

Changes to existing code

• common/jobs/pheasy.py: the harmonic fit of the pheasy phonon workflow and of this workflow now share one function. Two things change for the existing pheasy phonon workflow:
  • The pheasy commands run with subprocess.run(check=True). Before, subprocess.call ignored a failed pheasy run, and the workflow went on with missing or old force constants.
  • The LASSO fit runs with --tol 1e-8. With pheasy's default of 1e-4 the fit stops before it converges, and the force constants change between machines and package versions.
• ase/md.py: adds the nvt_nose_hoover_chain preset (ASE's NoseHooverChainNVT). This thermostat cannot follow a temperature schedule, so asking it to follow one now raises an error.
• vasp/sets/core.py: adds LangevinMDSetGenerator for NVT MD with a Langevin thermostat.

Tutorial

tutorials/finite_temperature_phonons.ipynb runs the workflow with three MACE potentials.

• Section 1 runs five materials from 100 to 900 K in the PBEsol cell of the Materials Project's Harmonic Phonon Database. At 300 K the potentials agree with effective harmonic force constants from AIMD to 0.07 to 0.43 THz on average over the band path.
• Section 2 runs the same materials with the NPT MD. It separates the effect of the force field's own 0 K cell from thermal expansion.
• Section 3 tests the MD length and the number of snapshots for KNaICl.

The loops are off by default and our results are hardcoded, as in the CTE tutorial. The notebook is left out of the nbmake job, since it needs a GPU and the MACE models.

Limitations

• A force field gives no Born charges, so the non-analytical correction needs VASP (born_maker).
• Symmetry is found without the magnetic moments, as in the other phonon workflows.
• There is no end-to-end mock_vasp test yet.

Checklist

• NumPy docstring format
• Type annotations
• Tests for the new code
• pre-commit and tests pass locally

[hrushikesh-s]

@JaGeo , I will send in the benchmarking tests on this workflow in the next few days.

[hrushikesh-s]

@JaGeo, given that both #1558 and #1559 are merged into main, I will get this PR ready for review.

[JaGeo]

This will now take a bit of time @hrushikesh-s . Maybe end of this week.

[hrushikesh-s]

@JaGeo , sounds good. Thanks!