Add a finite-temperature phonon workflow

.github/workflows/testing.yml

@@ -311,7 +311,7 @@ jobs:
           MP_API_KEY: ${{ secrets.MP_API_KEY }}
         run: |
           micromamba activate a2
-          pytest --durations=5 -n auto --nbmake ./tutorials --ignore=./tutorials/openmm_tutorial.ipynb --ignore=./tutorials/force_fields --ignore=./tutorials/torchsim_tutorial.ipynb --ignore=./tutorials/lammps_workflow.ipynb --ignore=./tutorials/pheasy_workflow.ipynb --ignore=./tutorials/hiphive_workflow.ipynb --ignore=./tutorials/cte_workflow.ipynb 
+          pytest --durations=5 -n auto --nbmake ./tutorials --ignore=./tutorials/openmm_tutorial.ipynb --ignore=./tutorials/force_fields --ignore=./tutorials/torchsim_tutorial.ipynb --ignore=./tutorials/lammps_workflow.ipynb --ignore=./tutorials/pheasy_workflow.ipynb --ignore=./tutorials/hiphive_workflow.ipynb --ignore=./tutorials/cte_workflow.ipynb --ignore=./tutorials/finite_temperature_phonons.ipynb 
 
       - name: Test ASE
         env:

docs/user/codes/vasp.md

@@ -349,6 +349,8 @@ See the [notes on MACE-Field](forcefields.md#mace-field-notes).
 
 Alternatively, users can accelerate the calculation of interatomic force constants using the machine-learning-based [Pheasy code](https://doi.org/10.48550/arXiv.2508.01020).
 The `pheasy` extra, `pip install "atomate2[pheasy]"`, installs the pheasy version this workflow needs, together with phonopy and ALM.
+The workflows call the `pheasy` command.
+A different command can be set with `PHEASY_CMD` in the atomate2 settings.
 ALM is compiled from source. If that build fails, see the ALM instructions below.
 By design, these workflows have the same basic structure as the harmonic forcefield workflows and use [Phonopy](https://doi.org/10.7566/JPSJ.92.012001) in part to compute the phonon spectrum.
 To use `Pheasy` in the previous example, we would replace the import string to `from atomate2.vasp.flows.pheasy import PhononMaker`.
@@ -411,6 +413,121 @@ The notebook `tutorials/hiphive_workflow.ipynb` reproduces it.
 Results on low-symmetry cells should be checked against the phonopy workflow before being trusted.
 The cluster space grows faster there than the number of displacements the workflow generates, and the fit can return spurious soft modes.
 
+#### Finite-temperature phonons
+
+`FiniteTemperaturePhononMaker` fits effective harmonic force constants at a finite temperature, as in the temperature-dependent effective potential ([TDEP](https://doi.org/10.1103/PhysRevB.84.180301)) method.
+They include the effect of the anharmonic forces on the frequencies at that temperature.
+They give no phonon lifetimes.
+By default, the MD runs at the volume of the relaxed structure, so thermal expansion is not included.
+It can be switched on with an NPT MD before the NVT MD, see below.
+The structure is relaxed first.
+An NVT MD run then samples the displacements at the temperature.
+By default, it runs for 8 ps at 300 K with a time step of 1 fs and a Langevin thermostat.
+The first 1 ps is left out, and 50 snapshots are picked evenly spread over the rest.
+Static calculations give the forces on the snapshots and on the undisplaced supercell.
+The displacements are measured from the undisplaced supercell.
+Pheasy then fits the second-order force constants to them with LASSO.
+As in the pheasy phonon workflow, a `DielectricMaker` computes the Born charges and the dielectric tensor for the non-analytical correction.
+It needs the `pheasy` extra, like the pheasy phonon workflow above.
+The notebook `tutorials/finite_temperature_phonons.ipynb` runs it for seven materials with three MACE potentials.
+
+```{warning}
+This workflow is new and has not been tested widely.
+It might still change in future versions.
+```
+
+```python
+from atomate2.vasp.flows.finite_temperature_phonons import (
+    FiniteTemperaturePhononMaker,
+)
+
+flow = FiniteTemperaturePhononMaker(temperature=300).make(structure)
+```
+
+The fit only uses the positions of the MD, not its dynamics.
+The MD has to sample the canonical ensemble.
+A Langevin thermostat does this for every mode ([Bussi and Parrinello](https://doi.org/10.1103/PhysRevE.75.056707)).
+For a harmonic oscillator weakly coupled to a Nosé-Hoover thermostat, the dynamics is not ergodic ([Legoll et al.](https://doi.org/10.1007/s00205-006-0029-1)).
+`thermostat="nose-hoover"` switches the NVT MD to a Nosé-Hoover thermostat with one thermostat variable.
+In the KNaICl test of the notebook, 12 of 27 runs with it had imaginary modes, against 19 of 27 with Langevin.
+
+The MD uses looser settings than the statics, for example ENCUT = 500 eV instead of 600 eV.
+Only the forces of the statics enter the fit.
+The relaxation, the MD and the statics use the same functional, +U, smearing and k-point settings.
+The MD and the statics start from the magnetic moments of the relaxed structure, if it has any.
+Both set ISPIN from the relaxation directory with `auto_ispin`.
+`md_runs` splits the MD into consecutive jobs, for example to stay within the walltime of a queue.
+Each job continues from the positions and velocities of the previous one.
+The thermostat variables start again from zero in each job.
+`ChainedMDMaker` in `atomate2.common.flows.md` makes these jobs, for VASP and for force fields.
+The MD files are read by the job that picks the snapshots, so the MD run directories must be on a file system that this job can read.
+
+An NPT MD can be run first to include thermal expansion.
+Set `npt_maker` to a VASP `MDMaker` with an `MDSetGenerator`, and use a larger ENCUT than for the NVT MD, since the cell changes.
+The flow runs it at the temperature and at `pressure`, 0 kbar by default, with MDALGO = 3 and ISIF = 3.
+VASP only runs NPT MD with its Langevin thermostat.
+The cell is averaged over the frames after `npt_equilibration_time`, given the symmetry of the relaxed structure, and used for the NVT MD, the statics and the fit.
+`fixed_cell_relax_maker` can relax the atoms in this cell, for example with ISIF = 2.
+
+The trajectory is also checked before the fit.
+The check looks for melting and for a move away from the reference structure.
+It also looks at the drift of the potential energy late in the run.
+The result is stored in `trajectory_health` of the output `FiniteTemperaturePhononDoc`, and a warning is raised when the check fails.
+The Lindemann ratio of 0.15 at melting is that of an fcc solid ([Saija et al.](https://doi.org/10.1063/1.2208357)).
+The other limits of the check are rules of thumb.
+`has_imaginary_modes` checks the band path.
+`n_imaginary_modes` and `lowest_frequency` use the q-points commensurate with the supercell.
+Imaginary modes are reported and not removed.
+
+The MD and the statics can also come from different codes.
+Three more makers are in `atomate2.forcefields.flows.finite_temperature_phonons`.
+`VaspMDMLFFStaticFiniteTemperaturePhononMaker` runs the MD with VASP and the statics with a force field.
+`MLFFMDVaspStaticFiniteTemperaturePhononMaker` runs the MD with a force field and the statics with VASP.
+`ForceFieldFiniteTemperaturePhononMaker` uses a force field for all steps.
+Each has a `from_force_field_name` method to choose the force field.
+They need the package of the force field, for example `mace-torch` for MACE.
+All of them use the same fit.
+A force field NPT MD uses ASE's `MTKNPT`, a Nosé-Hoover thermostat with one thermostat variable and the barostat of [Martyna et al.](https://doi.org/10.1063/1.467468).
+
+```python
+from atomate2.forcefields.flows.finite_temperature_phonons import (
+    ForceFieldFiniteTemperaturePhononMaker,
+)
+
+maker = ForceFieldFiniteTemperaturePhononMaker.from_force_field_name(
+    "MACE-MP-0", temperature=600, run_npt=True
+)
+flow = maker.make(structure)
+```
+
+`run_npt=True` sets the NPT MD and the relaxation of the atoms in its cell.
+
+Most force fields give no Born charges.
+`ForceFieldFiniteTemperaturePhononMaker` and `VaspMDMLFFStaticFiniteTemperaturePhononMaker` therefore skip the non-analytical correction by default.
+`MLFFMDVaspStaticFiniteTemperaturePhononMaker` keeps the VASP `DielectricMaker`.
+The Born charges and the dielectric tensor can come from MACE-Field, see the [force field notes](forcefields.md#mace-field-notes), or from VASP.
+Either pass `born` and `epsilon_static` to `make`, or set a `born_maker`:
+
+```python
+from atomate2.forcefields.jobs import ForceFieldDielectricMaker
+from atomate2.vasp.jobs.core import DielectricMaker
+
+# with MACE-Field
+maker.born_maker = ForceFieldDielectricMaker(
+    calculator_kwargs={"model": "/path/to/MACEField-MH-0-omat-dielectric.model"}
+)
+# or with VASP
+maker.born_maker = DielectricMaker()
+```
+
+Known limitations:
+
+- pheasy needs a diagonal supercell matrix.
+- The MD is classical, so there is no zero-point motion.
+- Each result comes from one MD trajectory, so it carries the noise of that trajectory.
+- The Lindemann ratio is averaged over all atoms. It can pass 0.15 when only the light atoms move far, for example Li in a superionic conductor, while the rest of the crystal stays solid.
+- With a force field MD and VASP statics, the trajectory and the forces come from different potential energy surfaces.
+
 ### Grüneisen parameter workflow
 
 Calculates mode-dependent Grüneisen parameters with the help of [Phonopy](https://doi.org/10.7566/JPSJ.92.012001).

src/atomate2/ase/md.py

@@ -3,6 +3,7 @@
 from __future__ import annotations
 
 import contextlib
+import inspect
 import io
 import logging
 import os
@@ -62,6 +63,7 @@ class DynamicsPresets(Enum):
     nvt_berendsen = "ase.md.nvtberendsen.NVTBerendsen"
     nvt_langevin = "ase.md.langevin.Langevin"
     nvt_nose_hoover = "ase.md.npt.NPT"
+    nvt_nose_hoover_chain = "ase.md.nose_hoover_chain.NoseHooverChainNVT"
     npt_berendsen = "ase.md.nptberendsen.NPTBerendsen"
     npt_nose_hoover = "ase.md.npt.NPT"  # noqa: PIE796
     npt_nose_hoover_chain = "ase.md.nose_hoover_chain.MTKNPT"
@@ -155,7 +157,8 @@ class AseMDMaker(AseMaker, ABC):
         The step interval for saving the trajectories.
     mb_velocity_seed : int or None
         If an int, a random number seed for generating initial velocities
-        from a Maxwell-Boltzmann distribution.
+        from a Maxwell-Boltzmann distribution. It also seeds the random forces
+        of dynamics that take an rng argument, such as Langevin.
     zero_linear_momentum : bool = False
         Whether to initialize the atomic velocities with zero linear momentum
     zero_angular_momentum : bool = False
@@ -380,13 +383,12 @@ def run_ase(
             # ASE NPT implementation requires upper triangular cell
             atoms.set_cell(atoms.cell.standard_form(form="upper")[0])
 
+        rng = np.random.default_rng(seed=self.mb_velocity_seed)
         if initial_velocities:
             atoms.set_velocities(initial_velocities)
         elif not np.isnan(self.t_schedule).any():
             MaxwellBoltzmannDistribution(
-                atoms=atoms,
-                temperature_K=self.t_schedule[0],
-                rng=np.random.default_rng(seed=self.mb_velocity_seed),
+                atoms=atoms, temperature_K=self.t_schedule[0], rng=rng
             )
             if self.zero_linear_momentum:
                 Stationary(atoms)
@@ -397,19 +399,34 @@ def run_ase(
 
         md_observer = TrajectoryObserver(atoms, store_md_outputs=True)
 
+        md_kwargs = dict(self.ase_md_kwargs)
+        if (
+            self.mb_velocity_seed is not None
+            and "rng" in inspect.signature(dynamics).parameters
+        ):
+            md_kwargs.setdefault("rng", rng)
         md_runner = dynamics(
-            atoms=atoms, timestep=self.time_step * units.fs, **self.ase_md_kwargs
+            atoms=atoms, timestep=self.time_step * units.fs, **md_kwargs
         )
 
         md_runner.attach(md_observer, interval=self.traj_interval)
 
+        can_set_temperature = hasattr(md_runner, "_temperature_K") or hasattr(
+            md_runner, "set_temperature"
+        )
+        if not can_set_temperature and np.ptp(self.t_schedule) > 0:
+            raise ValueError(
+                f"{type(md_runner).__name__} cannot follow a temperature schedule."
+            )
+
         def _callback(dyn: MolecularDynamics = md_runner) -> None:
             if self.ensemble == MDEnsemble.nve:
                 return
             if hasattr(dyn, "_temperature_K"):
                 dyn._temperature_K = self.t_schedule[dyn.nsteps]  # noqa: SLF001
-            else:
+            elif hasattr(dyn, "set_temperature"):
                 dyn.set_temperature(temperature_K=self.t_schedule[dyn.nsteps])
+            # NoseHooverChainNVT has neither and keeps its initial temperature
             if self.ensemble == MDEnsemble.nvt:
                 return