diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index 12670dcfab..925cb717d3 100644 --- a/.github/workflows/testing.yml +++ b/.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: diff --git a/docs/user/codes/vasp.md b/docs/user/codes/vasp.md index 06e4372bc2..1d8258cc6f 100644 --- a/docs/user/codes/vasp.md +++ b/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). diff --git a/src/atomate2/ase/md.py b/src/atomate2/ase/md.py index 4a52c29e2d..14ea89d751 100644 --- a/src/atomate2/ase/md.py +++ b/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 diff --git a/src/atomate2/common/flows/finite_temperature_phonons.py b/src/atomate2/common/flows/finite_temperature_phonons.py new file mode 100644 index 0000000000..f03a6442a1 --- /dev/null +++ b/src/atomate2/common/flows/finite_temperature_phonons.py @@ -0,0 +1,500 @@ +"""Flow for effective harmonic phonons at a finite temperature.""" + +from __future__ import annotations + +from abc import ABC, abstractmethod +from dataclasses import dataclass +from enum import Enum +from typing import TYPE_CHECKING, Literal + +import numpy as np +from jobflow import Flow, Maker, OutputReference +from pymatgen.util.due import Doi, due + +from atomate2 import SETTINGS +from atomate2.common.flows.md import ChainedMDMaker +from atomate2.common.jobs.finite_temperature_phonons import ( + fit_finite_temperature_phonons, + get_md_supercell, + get_npt_structure, + select_md_snapshots, +) +from atomate2.common.jobs.pheasy import get_supercell_size +from atomate2.common.jobs.phonons import run_phonon_displacements +from atomate2.common.utils import check_class_name + +if TYPE_CHECKING: + from pathlib import Path + + from emmet.core.math import Matrix3D + from jobflow import Job + from pymatgen.core import Structure + +SUPPORTED_CODES = frozenset(("vasp", "forcefields")) + + +def _get_force_field(maker: Maker, code: str) -> tuple[str | None, dict | None]: + """Get the calculator name and keyword arguments of a force field maker.""" + if code != "forcefields": + return None, None + name = maker.calculator_meta + return str(name.value if isinstance(name, Enum) else name), dict( + maker.calculator_kwargs + ) + + +@due.dcite( + Doi("10.48550/arXiv.2508.01020"), + description="Pheasy code for (an)harmonic force constants.", +) +@due.dcite( + Doi("10.1103/PhysRevB.84.180301"), + description="Effective harmonic force constants fitted to MD, as in TDEP.", +) +@due.dcite( + Doi("10.1103/PhysRevB.87.104111"), + description="Temperature dependent effective potential method.", +) +@due.dcite( + Doi("10.1088/1361-648X/acd831"), + description="Implementation strategies in phonopy and phono3py.", +) +@due.dcite( + Doi("10.7566/JPSJ.92.012001"), + description="Phonopy and phono3py.", +) +@dataclass +class BaseFiniteTemperaturePhononMaker(Maker, ABC): + """ + Maker for effective harmonic phonons at a finite temperature. + + The structure is relaxed first. An NVT MD run at the temperature then starts + from the undisplaced supercell. The frames of the first equilibration_time + are left out, and n_snapshots snapshots are picked evenly spread over the + rest of the trajectory. The phonon displacement maker computes the forces on + these snapshots and on the undisplaced supercell. The displacements of the + snapshots are measured from the relaxed positions. pheasy fits one set of + second-order force constants to the displacements and forces with LASSO. + Fitting effective harmonic force constants to MD is the idea of the + temperature-dependent effective potential method. The force constants + include the effect of the anharmonic forces on the frequencies, at the + temperature and at the volume of the relaxed structure. They give no + phonon lifetimes. The MD is classical. Thermal expansion is not included. + The force constants give the phonon band structure and density of states. + Imaginary modes are reported, not removed. + + If npt_maker is set, thermal expansion is included. An NPT MD run at the + temperature and pressure first starts from the undisplaced supercell of the + relaxed structure. Its cell, averaged over the frames after + npt_equilibration_time, gives the unit cell at the temperature, see + :obj:`.get_npt_structure`. If fixed_cell_relax_maker is set, it relaxes the + atoms in this cell. The NVT MD, the phonon displacement calculations and the + fit then use this cell. + + The MD maker and the phonon displacement maker are independent. Either can + be a VASP or a force field maker, so the forces can come from a different + level of theory than the trajectory. The residual forces of the undisplaced + supercell always come from the phonon displacement maker. + + The trajectory is also checked. 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 the output + document. A warning is raised if the check fails, but the fit is still done. + + This workflow is new and has not been tested widely. It might still change + in future versions. + + .. Note:: + The atoms of the relaxed structure are sorted by electronegativity + before the supercell is built, as the VASP input sets do. The atoms are + then in the same order in every calculation. The magnetic moments of the + relaxed structure, if any, are carried over to the supercell, the + snapshots and every MD job. phonopy and pheasy find the symmetry + without the magnetic moments, as in the atomate2 phonon workflows. The + band structure uses the seekpath k-path. + + Parameters + ---------- + name: str + Name of the flows produced by this maker. + temperature: float + MD temperature in K. + md_time: float + Length of the MD trajectory in ps. + md_time_step: float + MD time step in fs. + equilibration_time: float + Time at the start of the trajectory that is left out, in ps. + n_snapshots: int + Number of MD snapshots in the fit. + thermostat: Literal["nose-hoover", "langevin"] + Thermostat of the NVT MD. Langevin, the default, samples the canonical + ensemble also for nearly harmonic modes (Bussi and Parrinello, Phys. Rev. + E 75, 056707 (2007)). For a harmonic oscillator weakly coupled to a + Nose-Hoover thermostat, the dynamics is not ergodic (Legoll et al., + Arch. Ration. Mech. Anal. 184, 449 (2007)). The NPT MD does not use + this setting. + md_runs: int + Number of consecutive MD jobs that make up the trajectory. Each job + continues from the positions and velocities of the previous one. The + thermostat variables start again from zero in each job. + npt_time: float + Length of the NPT MD in ps, if there is one. + npt_equilibration_time: float + Time at the start of the NPT trajectory that is left out of the average + cell, in ps. + pressure: float + Pressure of the NPT MD in kbar. + min_length: float + Each lattice vector of the diagonal supercell is at least min_length + long, in Angstrom. + symprec: float + Symmetry precision for phonopy and pheasy. + rotational_sum_rule: Literal["BH", "H", "BHH"] | None + Rotational sum rule passed to pheasy with --rasr. None imposes none. + alpha_min: int + Base-10 exponent of the smallest LASSO penalty in the cross-validation, + passed to pheasy with --alpha_min. The default of -6 is pheasy's own + default. A warning is raised if the chosen penalty is on either bound of + the search. + random_seed: int | None + Seed of the LASSO fit. It also seeds the initial velocities and the + Langevin random forces of a force field MD. MD job i, counted from 0, + uses random_seed + i and the NPT MD uses random_seed + md_runs. + tol_imaginary_modes: float + Frequencies below -tol_imaginary_modes in THz are imaginary. + store_force_constants: bool + Whether to store the force constants in the output document. + code: str + Code of the phonon displacement calculations, 'vasp' or 'forcefields'. + md_code: str + Code of the MD, 'vasp' or 'forcefields'. + socket: bool + If True, the phonon displacement calculations run as a single batched + calculation. This is not supported by VASP. + bulk_relax_maker: Maker | None + Maker for the relaxation of the unit cell. None skips the relaxation. + born_maker: Maker | None + Maker for the Born effective charges and the dielectric tensor, used + for the non-analytical correction, for example a + ForceFieldDielectricMaker or a VASP DielectricMaker. None skips it. + npt_maker: Maker | None + Maker for the NPT MD. The flow sets its temperature, pressure, time + step and number of steps. None skips the NPT MD. + fixed_cell_relax_maker: Maker | None + Maker for the relaxation of the atoms in the cell from the NPT MD. It + must keep the cell. None keeps the fractional coordinates of the relaxed + structure. Only used with npt_maker. + md_maker: Maker + Maker for the MD. The flow sets its temperature, time step, number of + steps and thermostat. + phonon_displacement_maker: Maker + Maker for the static calculations on the snapshots and the undisplaced + supercell. + """ + + name: str = "finite temperature phonon" + temperature: float = 300.0 + md_time: float = 8.0 + md_time_step: float = 1.0 + equilibration_time: float = 1.0 + n_snapshots: int = 50 + thermostat: Literal["nose-hoover", "langevin"] = "langevin" + md_runs: int = 1 + npt_time: float = 8.0 + npt_equilibration_time: float = 2.0 + pressure: float = 0.0 + min_length: float = 12.0 + symprec: float = SETTINGS.PHONON_SYMPREC + rotational_sum_rule: Literal["BH", "H", "BHH"] | None = "BHH" + alpha_min: int = -6 + random_seed: int | None = 103 + tol_imaginary_modes: float = 0.1 + store_force_constants: bool = True + code: str | None = None + md_code: str | None = None + socket: bool = False + bulk_relax_maker: Maker | None = None + born_maker: Maker | None = None + npt_maker: Maker | None = None + fixed_cell_relax_maker: Maker | None = None + md_maker: Maker | None = None + phonon_displacement_maker: Maker | None = None + + def __post_init__(self) -> None: + """Check the settings before any calculation is run.""" + if self.md_maker is None or self.phonon_displacement_maker is None: + raise ValueError("md_maker and phonon_displacement_maker must be set.") + for key in ("code", "md_code"): + if getattr(self, key) not in SUPPORTED_CODES: + raise ValueError( + f"{key} must be one of {sorted(SUPPORTED_CODES)}, not " + f"{getattr(self, key)!r}." + ) + if self.alpha_min >= -2: + raise ValueError(f"alpha_min must be below -2, not {self.alpha_min}.") + if min(self.temperature, self.md_time, self.md_time_step) <= 0: + raise ValueError("temperature, md_time and md_time_step must be positive.") + n_steps = round(self.md_time * 1000 / self.md_time_step) + n_equil = round(self.equilibration_time * 1000 / self.md_time_step) + if self.equilibration_time < 0 or n_steps - n_equil < self.n_snapshots: + raise ValueError( + f"The MD has {n_steps} steps. After leaving out equilibration_time, " + f"fewer than the {self.n_snapshots} snapshots are left." + ) + if self.npt_maker is not None and not ( + 0 <= self.npt_equilibration_time < self.npt_time + ): + raise ValueError( + "npt_equilibration_time must be at least zero and shorter than " + "npt_time." + ) + + def get_md_steps(self) -> list[int]: + """Get the number of MD steps of each MD job.""" + n_steps = round(self.md_time * 1000 / self.md_time_step) + base, extra = divmod(n_steps, self.md_runs) + return [base + (idx < extra) for idx in range(self.md_runs)] + + def make( + self, + structure: Structure, + prev_dir: str | Path | None = None, + born: list[Matrix3D] | None = None, + epsilon_static: Matrix3D | None = None, + supercell_matrix: Matrix3D | None = None, + ) -> Flow: + """ + Make a flow to calculate effective harmonic phonons at a temperature. + + Parameters + ---------- + structure: Structure + The unit cell. + prev_dir: str | Path | None + A previous calculation directory. It is passed to the relaxation. + All later jobs get the relaxation directory, or this directory if + there is no relaxation. + born: list[Matrix3D] | None + Born effective charges of the structure the NVT MD runs in, with its + atoms sorted by electronegativity. This is the relaxed structure, or + the structure at the temperature with npt_maker. If given with + epsilon_static, the born_maker is not run. + epsilon_static: Matrix3D | None + High-frequency dielectric tensor. + supercell_matrix: Matrix3D | None + Diagonal supercell matrix. If None, it is chosen from min_length. + + Returns + ------- + Flow + """ + if supercell_matrix is not None and not isinstance( + supercell_matrix, OutputReference + ): + matrix = np.array(supercell_matrix) + if not np.allclose(matrix, np.diag(np.diag(matrix))): + raise ValueError("pheasy needs a diagonal supercell matrix.") + + # The VASP input sets sort the atoms by electronegativity. Sorting here as + # well keeps the Born charges of a force field born_maker in the order of + # the fit. + if not isinstance(structure, OutputReference): + structure = structure.get_sorted_structure() + + jobs: list[Job | Flow] = [] + optimization_run_job_dir = optimization_run_uuid = None + born_run_job_dir = born_run_uuid = None + npt_kwargs: dict = {} + + if self.bulk_relax_maker is not None: + relax = self.bulk_relax_maker.make(structure, prev_dir=prev_dir) + jobs.append(relax) + structure = relax.output.structure + prev_dir = relax.output.dir_name + optimization_run_job_dir = relax.output.dir_name + optimization_run_uuid = relax.output.uuid + + if supercell_matrix is None: + # pheasy needs a diagonal supercell matrix. With force_diagonal, + # the size only depends on min_length. + supercell_job = get_supercell_size( + structure, + self.min_length, + None, + force_90_degrees=False, + force_diagonal=True, + ) + jobs.append(supercell_job) + supercell_matrix = supercell_job.output + + if self.npt_maker is not None: + start_job = get_md_supercell( + structure, supercell_matrix, self.symprec, self.code + ) + n_steps = round(self.npt_time * 1000 / self.md_time_step) + npt_job = self.get_npt_maker(n_steps).make( + start_job.output, prev_dir=prev_dir + ) + npt_job.append_name(" NPT") + npt_cell_job = get_npt_structure( + npt_job.output.dir_name, + self.md_code, + structure, + supercell_matrix, + start_job.output, + self.md_time_step, + self.npt_equilibration_time, + self.symprec, + ) + jobs.extend([start_job, npt_job, npt_cell_job]) + npt_kwargs = { + "npt_input_structure": structure, + "pressure": self.pressure, + "npt_time": self.npt_time, + "npt_equilibration_time": self.npt_equilibration_time, + "npt_trajectory_health": npt_cell_job.output["trajectory_health"], + "npt_uuid": npt_job.uuid, + "npt_job_dir": npt_job.output.dir_name, + } + structure = npt_cell_job.output["structure"] + if self.fixed_cell_relax_maker is not None: + relax = self.fixed_cell_relax_maker.make(structure, prev_dir=prev_dir) + jobs.append(relax) + npt_kwargs["fixed_cell_relax_uuid"] = relax.uuid + npt_kwargs["fixed_cell_relax_job_dir"] = relax.output.dir_name + structure = relax.output.structure + + if self.born_maker is not None and (born is None or epsilon_static is None): + # as in the phonon workflow, so that a VASP born_maker can run after a + # force field relaxation, whose folder has no VASP files + born_kwargs = {} + if self.prev_calc_dir_argname is not None: + born_kwargs[self.prev_calc_dir_argname] = prev_dir + born_job = self.born_maker.make(structure, **born_kwargs) + jobs.append(born_job) + if check_class_name(self.born_maker, "ForceFieldDielectricMaker"): + born = born_job.output.born + epsilon_static = born_job.output.epsilon_static + else: + born = born_job.output.calcs_reversed[0].output.outcar["born"] + epsilon_static = born_job.output.calcs_reversed[0].output.epsilon_static + born_run_job_dir = born_job.output.dir_name + born_run_uuid = born_job.output.uuid + + reference_job = get_md_supercell( + structure, supercell_matrix, self.symprec, self.code + ) + jobs.append(reference_job) + reference = reference_job.output + + md_flow = ChainedMDMaker( + md_makers=[ + self.get_md_maker(n_steps, idx) + for idx, n_steps in enumerate(self.get_md_steps()) + ], + ).make(reference, prev_dir=prev_dir) + jobs.append(md_flow) + md_dirs = md_flow.output["dir_names"] + + snapshot_job = select_md_snapshots( + md_dirs, + self.md_code, + reference, + self.md_time_step, + self.equilibration_time, + self.n_snapshots, + ) + jobs.append(snapshot_job) + + displacement_calcs = run_phonon_displacements( + displacements=snapshot_job.output["structures"], + structure=structure, + supercell_matrix=supercell_matrix, + phonon_maker=self.phonon_displacement_maker, + socket=self.socket, + prev_dir_argname=self.prev_calc_dir_argname, + prev_dir=prev_dir, + store_displaced_structures=True, + ) + jobs.append(displacement_calcs) + + force_field_name, force_field_kwargs = _get_force_field( + self.phonon_displacement_maker, self.code + ) + md_force_field_name, md_force_field_kwargs = _get_force_field( + self.md_maker, self.md_code + ) + fit_job = fit_finite_temperature_phonons( + structure=structure, + supercell_matrix=supercell_matrix, + snapshot_data=snapshot_job.output, + displacement_data=displacement_calcs.output, + temperature=self.temperature, + thermostat=self.thermostat, + md_time_step=self.md_time_step, + equilibration_time=self.equilibration_time, + code=self.code, + md_code=self.md_code, + symprec=self.symprec, + force_field_name=force_field_name, + force_field_kwargs=force_field_kwargs, + md_force_field_name=md_force_field_name, + md_force_field_kwargs=md_force_field_kwargs, + md_uuids=md_flow.output["uuids"], + md_job_dirs=md_dirs, + optimization_run_uuid=optimization_run_uuid, + optimization_run_job_dir=optimization_run_job_dir, + born_run_uuid=born_run_uuid, + born_run_job_dir=born_run_job_dir, + rotational_sum_rule=self.rotational_sum_rule, + alpha_min=self.alpha_min, + random_seed=self.random_seed, + tol_imaginary_modes=self.tol_imaginary_modes, + born=born, + epsilon_static=epsilon_static, + store_force_constants=self.store_force_constants, + **npt_kwargs, + ) + jobs.append(fit_job) + return Flow(jobs, fit_job.output, name=self.name) + + @property + @abstractmethod + def prev_calc_dir_argname(self) -> str | None: + """Name of the prev_dir argument of the phonon displacement and Born makers. + + As this differs between codes, it is implemented by the inheriting class. + """ + + @abstractmethod + def get_npt_maker(self, n_steps: int) -> Maker: + """ + Get the maker of the NPT MD job, with the settings of this flow. + + Parameters + ---------- + n_steps: int + Number of MD steps. + + Returns + ------- + Maker + """ + + @abstractmethod + def get_md_maker(self, n_steps: int, index: int = 0) -> Maker: + """ + Get the maker of one MD job, with the settings of this flow. + + Parameters + ---------- + n_steps: int + Number of MD steps of the job. + index: int + Position of the job in the MD, counted from 0. + + Returns + ------- + Maker + """ diff --git a/src/atomate2/common/flows/md.py b/src/atomate2/common/flows/md.py new file mode 100644 index 0000000000..1739cdebb4 --- /dev/null +++ b/src/atomate2/common/flows/md.py @@ -0,0 +1,88 @@ +"""Flows for molecular dynamics that are shared between codes.""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import TYPE_CHECKING + +from jobflow import Flow, Maker + +from atomate2.ase.md import AseMDMaker +from atomate2.common.jobs.md import get_md_restart_structure + +if TYPE_CHECKING: + from pathlib import Path + + from jobflow import Job + from pymatgen.core import Structure + + +@dataclass +class ChainedMDMaker(Maker): + """ + Maker to run one MD as several consecutive MD jobs, with VASP or a force field. + + The first MD job starts from the input structure. Each later one starts from + the final positions and velocities of the previous one. The thermostat + variables are not carried over, so they start again from zero in each job. + + Parameters + ---------- + name: str + Name of the flows produced by this maker. + md_makers: list[Maker] + Makers of the MD jobs, in the order of the trajectory. VASP makers and + ASE makers, such as the force field ones, are supported. An ASE maker + that is followed by another job must write its trajectory in the ASE + format, with traj_file set. + """ + + name: str = "chained MD" + md_makers: list[Maker] = field(default_factory=list) + + def make(self, structure: Structure, prev_dir: str | Path | None = None) -> Flow: + """ + Make a flow of consecutive MD jobs. + + Parameters + ---------- + structure: Structure + The structure the first MD job starts from. Its magnetic moments, if + any, are also given to the later MD jobs. + prev_dir: str | Path | None + A previous calculation directory, passed to each MD job. + + Returns + ------- + Flow + Its output holds the directories and the uuids of the MD jobs, in + the order of the trajectory. + """ + jobs: list[Job] = [] + md_jobs: list[Job] = [] + md_structure = structure + for idx, maker in enumerate(self.md_makers, start=1): + md_job = maker.make(md_structure, prev_dir=prev_dir) + if len(self.md_makers) > 1: + md_job.append_name(f" {idx}/{len(self.md_makers)}") + jobs.append(md_job) + md_jobs.append(md_job) + if idx < len(self.md_makers): + traj_file = None + if isinstance(maker, AseMDMaker): + if maker.traj_file is None or maker.traj_file_fmt != "ase": + raise ValueError( + f"{maker.name} must write its trajectory in the ASE " + "format, so that the next MD job can continue from it." + ) + traj_file = str(maker.traj_file) + restart = get_md_restart_structure( + md_job.output.dir_name, structure, traj_file + ) + jobs.append(restart) + md_structure = restart.output + output = { + "dir_names": [md_job.output.dir_name for md_job in md_jobs], + "uuids": [md_job.uuid for md_job in md_jobs], + } + return Flow(jobs, output, name=self.name) diff --git a/src/atomate2/common/flows/pheasy.py b/src/atomate2/common/flows/pheasy.py index 910b606866..5e3b215254 100644 --- a/src/atomate2/common/flows/pheasy.py +++ b/src/atomate2/common/flows/pheasy.py @@ -142,9 +142,10 @@ class BasePhononMaker(PurePhonopyMaker, ABC): cutoff distance in Bohr for each FC order, starting at second order. The default value is [-1, 12, 10]. The first entry is not used, since pheasy fits the second-order FCs without a cutoff. The second and third entries - are the cutoffs for third- and fourth-order FCs. The cutoff of each fitted - order, up to anhar_max_order, must be positive. Longer cutoffs increase - the number of free FCs, and with it the number of displaced supercells. + are the cutoffs for third- and fourth-order FCs. The defaults of 12 and 10 + Bohr are 6.35 and 5.29 Å. The cutoff of each fitted order, up to + anhar_max_order, must be positive. Longer cutoffs increase the number of + free FCs, and with it the number of displaced supercells. min_length: float minimum length of lattice constants will be used to create the supercell, the default value is 8.0 A. It can be increased for larger supercells. diff --git a/src/atomate2/common/jobs/finite_temperature_phonons.py b/src/atomate2/common/jobs/finite_temperature_phonons.py new file mode 100644 index 0000000000..ad5483c948 --- /dev/null +++ b/src/atomate2/common/jobs/finite_temperature_phonons.py @@ -0,0 +1,871 @@ +"""Jobs for effective harmonic phonons from MD snapshots and pheasy.""" + +from __future__ import annotations + +import logging +import warnings +from pathlib import Path +from typing import TYPE_CHECKING + +import numpy as np +from ase.io import read as ase_read +from ase.units import kB +from jobflow import job +from monty.io import zopen +from monty.os.path import zpath +from phonopy.file_IO import parse_FORCE_CONSTANTS +from phonopy.harmonic.dynmat_to_fc import get_commensurate_points +from phonopy.interface.vasp import write_vasp +from pymatgen.core import Structure +from pymatgen.io.phonopy import get_phonopy_structure, get_pmg_structure +from pymatgen.io.vasp import Incar, Kpoints, Xdatcar +from pymatgen.phonon.bandstructure import PhononBandStructureSymmLine +from pymatgen.phonon.dos import PhononDos +from pymatgen.symmetry.analyzer import SpacegroupAnalyzer + +from atomate2.common.jobs.md import _get_site_properties +from atomate2.common.jobs.pheasy import ( + _DEFAULT_FILE_PATHS, + _check_lasso_alpha, + _run_harmonic_fit, +) +from atomate2.common.jobs.phonons import ( + _generate_phonon_object, + _get_kpath, + _run_band_structure_and_plot, + _run_total_dos_and_plot, +) +from atomate2.common.schemas.finite_temperature_phonons import ( + FiniteTemperaturePhononDoc, + TrajectoryHealth, +) +from atomate2.common.schemas.phonons import ( + ForceConstants, + PhononComputationalSettings, + PhononJobDirs, + PhononUUIDs, + _set_nac_params, +) +from atomate2.utils.path import strip_hostname + +if TYPE_CHECKING: + from collections.abc import Sequence + + from emmet.core.math import Matrix3D + from phonopy import Phonopy + +logger = logging.getLogger(__name__) + +# file of the force field MD trajectory, written in the ASE format +ASE_TRAJECTORY_FILE = "md_trajectory.traj" +_KPATH_SCHEME = "seekpath" + +# Limits of the trajectory check +# Rule of thumb: in the first 10 steps thermal motion moves an atom by about a +# tenth of an Angstrom or less, while a wrong atom order moves it by a bond length. +_N_START_FRAMES = 10 +_START_LIMIT = 0.5 # Angstrom +# Lindemann ratio at melting of an fcc solid. It is about 0.18 for a bcc solid. +# Saija et al., J. Chem. Phys. 124, 244504 (2006) +_LINDEMANN_LIMIT = 0.15 +# Rule of thumb: u_ref / u_vib = 1.5 means the mean positions moved by about as +# much as the atoms vibrate (u_shift = 1.1 u_vib). +_SHIFT_RATIO_LIMIT = 1.5 +# Rule of thumb: the mean positions are only checked if the second half of the +# trajectory lasts at least one period of a 1 THz vibration. +_MIN_SHIFT_WINDOW = 1.0 # ps +# Rule of thumb: three standard deviations of the energy of single frames +_DRIFT_SIGMA = 3.0 +# Rule of thumb: a third of a typical nearest-neighbor distance of 3 Angstrom, +# about twice the vibration at the Lindemann limit for that distance. +_RMS_LIMIT = 1.0 # Angstrom + + +def _get_phonopy( + structure: Structure, supercell_matrix: Matrix3D, symprec: float, code: str +) -> Phonopy: + """Get the phonopy object of the sorted structure and the supercell.""" + return _generate_phonon_object( + structure.get_sorted_structure(), + supercell_matrix, + displacement=0.01, + sym_reduce=True, + symprec=symprec, + use_symmetrized_structure=None, + kpath_scheme=_KPATH_SCHEME, + code=code, + ) + + +@job +def get_md_supercell( + structure: Structure, supercell_matrix: Matrix3D, symprec: float, code: str +) -> Structure: + """ + Build the supercell used as the reference structure of the MD and the fit. + + The atoms of the structure are sorted by electronegativity, as the VASP + input sets do. The supercell is then the one phonopy builds, so that its + atoms are in the order of the force constant fit. Each atom of the + supercell gets the magnetic moment of its atom in the unit cell, if the + structure has magnetic moments. + + Parameters + ---------- + structure: Structure + Relaxed unit cell. + supercell_matrix: Matrix3D + Supercell matrix. + symprec: float + Symmetry precision for phonopy. + code: str + Code of the phonon displacement calculations. + + Returns + ------- + Structure + The undisplaced supercell. + """ + structure = structure.get_sorted_structure() + supercell = _get_phonopy(structure, supercell_matrix, symprec, code).supercell + site_properties = {} + if "magmom" in structure.site_properties: + unit_index = [supercell.u2u_map[idx] for idx in supercell.s2u_map] + magmoms = structure.site_properties["magmom"] + site_properties["magmom"] = [magmoms[idx] for idx in unit_index] + return Structure( + supercell.cell, + supercell.symbols, + supercell.scaled_positions, + site_properties=site_properties, + ) + + +def _read_vasp_md( + directory: Path, +) -> tuple[np.ndarray, np.ndarray, np.ndarray, list[str], float]: + """Read the frames, cells, energies and time step of a VASP MD run.""" + incar = Incar.from_file(zpath(str(directory / "INCAR"))) + xdatcar = Xdatcar(zpath(str(directory / "XDATCAR"))) + frac_coords = np.array([frame.frac_coords for frame in xdatcar.structures]) + cells = np.array([frame.lattice.matrix for frame in xdatcar.structures]) + species = [str(site.specie) for site in xdatcar.structures[0]] + + # the free energy F of each ionic step + with zopen(zpath(str(directory / "OSZICAR")), mode="rt") as file: + energies = [ + float(line.split(" F= ")[1].split()[0]) + for line in file + if " T= " in line and " F= " in line + ] + if len(energies) != len(frac_coords): + raise ValueError( + f"{directory} has {len(frac_coords)} frames in XDATCAR but " + f"{len(energies)} MD steps in OSZICAR. XDATCAR must have a frame at " + "every step, with NBLOCK = 1." + ) + return frac_coords, cells, np.array(energies), species, float(incar["POTIM"]) + + +def _read_ase_md( + directory: Path, time_step: float +) -> tuple[np.ndarray, np.ndarray, np.ndarray, list[str], float]: + """Read the frames, cells and energies of a force field MD run.""" + # the first frame is the starting structure, before the first step, which + # XDATCAR leaves out as well + frames = ase_read(directory / ASE_TRAJECTORY_FILE, index=":")[1:] + frac_coords = np.array([atoms.get_scaled_positions() for atoms in frames]) + cells = np.array([atoms.cell[:] for atoms in frames]) + energies = np.array([atoms.get_potential_energy() for atoms in frames]) + return frac_coords, cells, energies, frames[0].get_chemical_symbols(), time_step + + +def _read_md( + md_dirs: Sequence[str], md_code: str, reference: Structure, md_time_step: float +) -> tuple[np.ndarray, np.ndarray, np.ndarray, float]: + """ + Read and join the trajectories of MD runs, in order. + + For VASP, XDATCAR, OSZICAR and the POTIM of INCAR are read from each run + directory. For force fields, the ASE trajectory file is read, and + md_time_step is the time step. Returns the fractional coordinates, the cell + and the potential energy of each frame, and the time step in fs. + """ + frac_coords, cells, energies = [], [], [] + for md_dir in md_dirs: + directory = Path(strip_hostname(md_dir)) + if md_code == "vasp": + coords, cell, energy, species, time_step = _read_vasp_md(directory) + else: + coords, cell, energy, species, time_step = _read_ase_md( + directory, md_time_step + ) + if species != [str(site.specie) for site in reference]: + raise ValueError( + f"The atoms of the MD in {directory} are not in the order of the " + "reference supercell." + ) + frac_coords.append(coords) + cells.append(cell) + energies.append(energy) + return ( + np.concatenate(frac_coords), + np.concatenate(cells), + np.concatenate(energies), + time_step, + ) + + +def _get_displacements(frac_coords: np.ndarray, reference: Structure) -> np.ndarray: + """Cartesian displacements from the reference with the minimum image.""" + diff = frac_coords - reference.frac_coords + diff -= np.round(diff) + return diff @ reference.lattice.matrix + + +def _remove_center_of_mass(disps: np.ndarray, reference: Structure) -> np.ndarray: + """Remove the displacement of the center of mass from each frame.""" + masses = np.array([site.specie.atomic_mass for site in reference]) + center = np.einsum("i,fia->fa", masses, disps) / masses.sum() + return disps - center[:, None, :] + + +def _assess_trajectory( + frac_coords: np.ndarray, + energies: np.ndarray, + reference: Structure, + time_step: float, +) -> TrajectoryHealth: + """ + Check whether the MD trajectory stayed at the reference structure. + + The displacement of the center of mass of each frame is removed from its + displacements. The checks are applied in this order. The limits are rules + of thumb, except where a source is given next to them in this module. + + - A root mean square displacement above 0.5 Angstrom in the first 10 + frames means the atoms are not in the order of the reference. + - A root mean square vibration u_vib above 0.15 of the nearest-neighbor + distance means the structure melted. u_vib and u_ref are measured in + the second half of the frames. + - A total displacement u_ref above 1.5 times u_vib means the mean + positions moved away from the reference. This check is skipped if the + second half of the trajectory is shorter than 1 ps. + - A change of the mean potential energy from the second to the last + quarter of the frames above 3 times the standard deviation of the + energy in the last fifth means the structure transformed (falling energy) + or is disordering (rising energy). This check is skipped below 20 + frames. + - A root mean square displacement above 1 Angstrom in the last fifth of the + frames points at diffusion or very soft motion. + + The energy drift is measured late in the trajectory, because a run that + starts from the reference must gain 3/2 k_B T of potential energy to reach + equipartition. + + Parameters + ---------- + frac_coords: np.ndarray + Fractional coordinates of each frame, with shape (n_frames, n_atoms, 3). + energies: np.ndarray + Potential energy of each frame in eV. + reference: Structure + The undisplaced supercell. + time_step: float + MD time step in fs. + + Returns + ------- + TrajectoryHealth + """ + n_frames, n_atoms = frac_coords.shape[:2] + disps = _remove_center_of_mass( + np.array([_get_displacements(frame, reference) for frame in frac_coords]), + reference, + ) + rms = np.sqrt(np.mean(np.sum(disps**2, axis=2), axis=1)) + rms_start = float(rms[:_N_START_FRAMES].mean()) + rms_end = float(rms[-max(1, n_frames // 5) :].mean()) + + # split the displacement in the second half into a static and a vibrational + # part, u_ref**2 = u_shift**2 + u_vib**2 + tail = disps[n_frames // 2 :] + msd_ref = float(np.mean(np.sum(tail**2, axis=2))) + msd_shift = float(np.mean(np.sum(tail.mean(axis=0) ** 2, axis=1))) + u_ref = np.sqrt(msd_ref) + u_shift = np.sqrt(msd_shift) + u_vib = np.sqrt(max(msd_ref - msd_shift, 0.0)) + + distances = reference.distance_matrix.copy() + np.fill_diagonal(distances, np.inf) + d_nn = float(distances.min(axis=1).mean()) + lindemann_ratio = u_vib / d_nn + shift_ratio = u_ref / max(u_vib, 1e-9) + shifted = ( + shift_ratio > _SHIFT_RATIO_LIMIT + and len(tail) * time_step / 1000 >= _MIN_SHIFT_WINDOW + ) + + energy = energies / n_atoms + quarter = len(energy) // 4 + drift = ( + float(energy[3 * quarter :].mean() - energy[quarter : 2 * quarter].mean()) + if quarter >= 5 + else None + ) + fluctuation = float(energy[-max(1, len(energy) // 5) :].std()) + big_drift = drift is not None and abs(drift) > _DRIFT_SIGMA * max(fluctuation, 1e-9) + + if rms_start > _START_LIMIT: + verdict = "reference_mismatch" + elif lindemann_ratio > _LINDEMANN_LIMIT: + verdict = "melted" + elif shifted: + verdict = "transformed" if big_drift and drift < 0 else "shifted_or_diffusing" + elif big_drift: + verdict = "transformed" if drift < 0 else "disordering" + elif rms_end > _RMS_LIMIT: + verdict = "diffusing_or_soft" + else: + verdict = "stable" + + return TrajectoryHealth( + verdict=verdict, + n_frames=n_frames, + rms_displacement_start=rms_start, + rms_displacement_end=rms_end, + u_ref=u_ref, + u_shift=u_shift, + u_vib=u_vib, + nearest_neighbor_distance=d_nn, + lindemann_ratio=lindemann_ratio, + shift_ratio=shift_ratio, + energy_drift=drift, + energy_fluctuation=fluctuation, + ) + + +def average_npt_structure( + cells: np.ndarray, + structure: Structure, + supercell_matrix: Matrix3D, + symprec: float, +) -> Structure: + """ + Average the supercell lattices of an NPT MD into a unit cell. + + The metric tensor L L^T of the supercell, with the lattice vectors as the + rows of L, is averaged over the cells. Unlike the lattice vectors, it does + not change when the cell rotates. It is converted to the metric tensor of + the unit cell with the supercell matrix and averaged over the point group of + the structure, so that the cell keeps its symmetry. The new lattice is the + stretch of the lattice of the structure, without a rotation, that has this + metric tensor. The atoms keep their fractional coordinates. + + Parameters + ---------- + cells: np.ndarray + Supercell lattices of the frames to average, with the lattice vectors + as rows, in Angstrom. + structure: Structure + Unit cell whose supercell the NPT MD started from. + supercell_matrix: Matrix3D + Supercell matrix. + symprec: float + Symmetry precision for the point group of the structure. + + Returns + ------- + Structure + The unit cell with the averaged lattice. + """ + metric = np.mean([cell @ cell.T for cell in cells], axis=0) + inv_matrix = np.linalg.inv(np.array(supercell_matrix, dtype=float)) + metric = inv_matrix @ metric @ inv_matrix.T + # a rotation W of fractional coordinates leaves the metric tensor of the + # symmetric cell unchanged, W^T G W = G + analyzer = SpacegroupAnalyzer(structure, symprec=symprec) + rotations = [op.rotation_matrix for op in analyzer.get_symmetry_operations()] + metric = np.mean([rot.T @ metric @ rot for rot in rotations], axis=0) + lattice = structure.lattice.matrix + inv_lattice = np.linalg.inv(lattice) + eigvals, eigvecs = np.linalg.eigh(inv_lattice @ metric @ inv_lattice.T) + stretch = eigvecs @ np.diag(np.sqrt(eigvals)) @ eigvecs.T + return Structure( + lattice @ stretch, + structure.species, + structure.frac_coords, + site_properties=structure.site_properties, + ) + + +@job +def get_npt_structure( + npt_dir: str, + md_code: str, + structure: Structure, + supercell_matrix: Matrix3D, + reference: Structure, + md_time_step: float, + equilibration_time: float, + symprec: float, +) -> dict: + """ + Get the unit cell at the temperature from an NPT MD run. + + The cells after equilibration_time are averaged with + :obj:`average_npt_structure`. The NPT trajectory is checked as in + :obj:`select_md_snapshots`. + + Parameters + ---------- + npt_dir: str + Directory of the NPT MD run. + md_code: str + Code of the MD, "vasp" or "forcefields". + structure: Structure + Unit cell whose supercell the NPT MD started from. + supercell_matrix: Matrix3D + Supercell matrix. + reference: Structure + The undisplaced supercell the NPT MD started from. + md_time_step: float + MD time step in fs. Only used for force field trajectories. + equilibration_time: float + Time at the start of the trajectory that is left out, in ps. + symprec: float + Symmetry precision for the point group of the structure. + + Returns + ------- + dict + The unit cell at the temperature and the check of the NPT trajectory. + """ + frac_coords, cells, energies, time_step = _read_md( + [npt_dir], md_code, reference, md_time_step + ) + n_equil = round(equilibration_time * 1000 / time_step) + if n_equil >= len(cells): + raise ValueError( + f"The NPT trajectory has {len(cells)} frames of {time_step} fs, none " + f"of them after the {equilibration_time} ps that are left out." + ) + npt_structure = average_npt_structure( + cells[n_equil:], structure, supercell_matrix, symprec + ) + + health = _assess_trajectory(frac_coords, energies, reference, time_step) + if health.verdict != "stable": + warnings.warn( + f"The NPT trajectory did not stay at the reference structure (verdict: " + f"{health.verdict}). The averaged cell may not describe it.", + stacklevel=2, + ) + return {"structure": npt_structure, "trajectory_health": health.model_dump()} + + +@job(data=["structures"]) +def select_md_snapshots( + md_dirs: Sequence[str], + md_code: str, + reference: Structure, + md_time_step: float, + equilibration_time: float, + n_snapshots: int, +) -> dict: + """ + Pick snapshots from the MD trajectory and check the trajectory. + + The trajectories of all MD runs are joined in order, with a frame after + every time step. The time of a frame is its step number times the time + step. The frames of the first equilibration_time are left out, and + n_snapshots frames are picked evenly spread over the rest, from the first + to the last. For VASP, XDATCAR, OSZICAR and the POTIM of INCAR are + read from each run directory. For force fields, the ASE trajectory file is + read. The run directories must be readable from where this job runs. + + Parameters + ---------- + md_dirs: Sequence[str] + Directories of the MD runs, in order. + md_code: str + Code of the MD, "vasp" or "forcefields". + reference: Structure + The undisplaced supercell the MD started from. + md_time_step: float + MD time step in fs. Only used for force field trajectories. + equilibration_time: float + Time at the start of the trajectory that is left out, in ps. + n_snapshots: int + Number of snapshots. + + Returns + ------- + dict + The structures for the phonon displacement calculations, the + snapshots followed by the undisplaced supercell. The MD time of each + snapshot in ps. The number of frames times the time step in ps. The + root mean square displacement of the snapshots in Angstrom, with the + displacement of the center of mass removed. The trajectory check. + """ + all_coords, _, all_energies, time_step = _read_md( + md_dirs, md_code, reference, md_time_step + ) + n_frames = len(all_coords) + + n_equil = round(equilibration_time * 1000 / time_step) + n_left = n_frames - n_equil + if n_left < n_snapshots: + raise ValueError( + f"The trajectory has {n_frames} frames of {time_step} fs. After " + f"leaving out {equilibration_time} ps, {max(n_left, 0)} frames are " + f"left, fewer than the {n_snapshots} snapshots." + ) + indices = np.linspace(n_equil, n_frames - 1, n_snapshots).round().astype(int) + + site_properties = _get_site_properties(reference) + snapshots = [ + Structure( + reference.lattice, + reference.species, + all_coords[idx], + site_properties=site_properties, + ) + for idx in indices + ] + disps = _remove_center_of_mass( + np.array([_get_displacements(all_coords[idx], reference) for idx in indices]), + reference, + ) + health = _assess_trajectory(all_coords, all_energies, reference, time_step) + if health.verdict != "stable": + warnings.warn( + f"The MD trajectory did not stay at the reference structure (verdict: " + f"{health.verdict}). The fitted force constants may not describe it.", + stacklevel=2, + ) + return { + "structures": [*snapshots, reference], + "snapshot_times": [(idx + 1) * time_step / 1000 for idx in indices], + "md_time": n_frames * time_step / 1000, + "rms_displacement": float(np.sqrt(np.mean(np.sum(disps**2, axis=2)))), + "trajectory_health": health.model_dump(), + } + + +def _get_rms_displacement(phonon: Phonopy, temperature: float, tol: float) -> float: + """ + Get the classical root mean square displacement from the force constants. + + The modes are those of the supercell at the Gamma point. A mode with a + frequency above tol in THz adds k_B T divided by its eigenvalue of the + dynamical matrix, times the sum over the atoms and Cartesian components of + its squared eigenvector component divided by the atomic mass. The sum over + the modes is divided by the number of atoms. The modes at or below tol are + left out. They include the three acoustic modes at the Gamma point, which + move the center of mass. + """ + masses = np.repeat(np.array(phonon.supercell.masses), 3) + n_dof = len(masses) + force_constants = phonon.force_constants.transpose(0, 2, 1, 3).reshape(n_dof, n_dof) + inv_sqrt_mass = 1 / np.sqrt(masses) + dynmat = force_constants * np.outer(inv_sqrt_mass, inv_sqrt_mass) + eigvals, eigvecs = np.linalg.eigh((dynmat + dynmat.T) / 2) + frequencies = ( + np.sign(eigvals) * np.sqrt(np.abs(eigvals)) * phonon.unit_conversion_factor + ) + keep = frequencies > tol + weights = np.sum(eigvecs[:, keep] ** 2 / masses[:, None], axis=0) + msd = kB * temperature * np.sum(weights / eigvals[keep]) + return float(np.sqrt(msd / len(phonon.supercell.masses))) + + +@job( + output_schema=FiniteTemperaturePhononDoc, + data=[PhononDos, PhononBandStructureSymmLine, ForceConstants], +) +def fit_finite_temperature_phonons( + structure: Structure, + supercell_matrix: Matrix3D, + snapshot_data: dict, + displacement_data: dict, + temperature: float, + thermostat: str, + md_time_step: float, + equilibration_time: float, + code: str, + md_code: str, + symprec: float, + force_field_name: str | None = None, + force_field_kwargs: dict | None = None, + md_force_field_name: str | None = None, + md_force_field_kwargs: dict | None = None, + md_uuids: list[str] | None = None, + md_job_dirs: list[str] | None = None, + optimization_run_uuid: str | None = None, + optimization_run_job_dir: str | None = None, + born_run_uuid: str | None = None, + born_run_job_dir: str | None = None, + npt_input_structure: Structure | None = None, + pressure: float | None = None, + npt_time: float | None = None, + npt_equilibration_time: float | None = None, + npt_trajectory_health: dict | None = None, + npt_uuid: str | None = None, + npt_job_dir: str | None = None, + fixed_cell_relax_uuid: str | None = None, + fixed_cell_relax_job_dir: str | None = None, + rotational_sum_rule: str | None = "BHH", + alpha_min: int = -6, + random_seed: int | None = 103, + tol_imaginary_modes: float = 0.1, + born: list[Matrix3D] | None = None, + epsilon_static: Matrix3D | None = None, + store_force_constants: bool = True, + npoints_band: int = 101, + kpoint_density_dos: int = 7_000, +) -> FiniteTemperaturePhononDoc: + """ + Fit effective second-order force constants to MD snapshots with pheasy. + + The forces on the undisplaced supercell, the last phonon displacement + calculation, are subtracted from the forces on the snapshots. pheasy then + fits the second-order force constants with LASSO, with standardized data + and without a cutoff. phonopy imposes translational and permutation + symmetry on them. They give the phonon band structure along the seekpath + k-path and the density of states. Imaginary modes are reported, not + removed. + + Parameters + ---------- + structure: Structure + Unit cell of the NVT MD, the relaxed unit cell or the cell from the NPT + MD. This job sorts its atoms by electronegativity. + supercell_matrix: Matrix3D + Diagonal supercell matrix. + snapshot_data: dict + Output of :obj:`select_md_snapshots`. + displacement_data: dict + Output of the phonon displacement calculations on the snapshots and on + the undisplaced supercell, which comes last. + temperature: float + MD temperature in K. + thermostat: str + Thermostat of the MD. + md_time_step: float + MD time step in fs. + equilibration_time: float + Time at the start of the trajectory that was left out, in ps. + code: str + Code of the phonon displacement calculations. + md_code: str + Code of the MD. + symprec: float + Symmetry precision for phonopy and pheasy. + force_field_name: str | None + Force field of the phonon displacement calculations, if any. + force_field_kwargs: dict | None + Keyword arguments of the force field calculator of the phonon + displacement calculations. + md_force_field_name: str | None + Force field of the MD, if any. + md_force_field_kwargs: dict | None + Keyword arguments of the force field calculator of the MD. + md_uuids: list[str] | None + UUIDs of the MD jobs. + md_job_dirs: list[str] | None + Directories of the MD jobs. + optimization_run_uuid: str | None + UUID of the relaxation. + optimization_run_job_dir: str | None + Directory of the relaxation. + born_run_uuid: str | None + UUID of the Born charge calculation. + born_run_job_dir: str | None + Directory of the Born charge calculation. + npt_input_structure: Structure | None + Unit cell whose supercell the NPT MD started from, if there was one. + pressure: float | None + Pressure of the NPT MD in kbar. + npt_time: float | None + Length of the NPT MD in ps. + npt_equilibration_time: float | None + Time at the start of the NPT trajectory that was left out, in ps. + npt_trajectory_health: dict | None + Check of the NPT trajectory. + npt_uuid: str | None + UUID of the NPT MD job. + npt_job_dir: str | None + Directory of the NPT MD job. + fixed_cell_relax_uuid: str | None + UUID of the relaxation of the atoms in the cell from the NPT MD. + fixed_cell_relax_job_dir: str | None + Directory of the relaxation of the atoms in the cell from the NPT MD. + rotational_sum_rule: str | None + Rotational sum rule passed to pheasy with --rasr, or None for none. + alpha_min: int + Base-10 exponent of the smallest LASSO penalty in the cross-validation. + random_seed: int | None + Seed of the LASSO fit in pheasy. + tol_imaginary_modes: float + Frequencies below -tol_imaginary_modes in THz are imaginary. + born: list[Matrix3D] | None + Born effective charges of the sorted unit cell, for the non-analytical + correction of the dynamical matrix. + epsilon_static: Matrix3D | None + High-frequency dielectric tensor, needed together with born. + store_force_constants: bool + Whether to store the force constants in the output document. + npoints_band: int + Number of q-points per band structure segment. + kpoint_density_dos: int + Number of q-points per reciprocal atom of the density of states mesh. + + Returns + ------- + FiniteTemperaturePhononDoc + """ + supercell_matrix = np.array(supercell_matrix) + structure = structure.get_sorted_structure() + if born is not None and len(born) != len(structure): + raise ValueError("The number of Born charges is not the number of atoms.") + phonon = _get_phonopy(structure, supercell_matrix, symprec, code) + supercell = get_pmg_structure(phonon.supercell) + + forces = np.array(displacement_data["forces"]) + structures = displacement_data["displaced_structures"] + residual_forces = forces[-1] + fit_forces = forces[:-1] - residual_forces + disps = np.array( + [_get_displacements(s.frac_coords, supercell) for s in structures[:-1]] + ) + n_data = len(disps) + + write_vasp("POSCAR", get_phonopy_structure(structure)) + write_vasp("SPOSCAR", phonon.supercell) + np.save(_DEFAULT_FILE_PATHS["harmonic_displacements"], disps) + np.save(_DEFAULT_FILE_PATHS["harmonic_force_matrix"], fit_forces) + + log_file = Path(_DEFAULT_FILE_PATHS["harmonic_fit_log"]) + fc_file = Path(_DEFAULT_FILE_PATHS["force_constants"]) + _run_harmonic_fit( + supercell_matrix, + symprec, + n_data, + rotational_sum_rule=rotational_sum_rule, + alpha_min=alpha_min, + random_seed=random_seed, + log_file=str(log_file), + ) + + if not fc_file.exists(): + raise RuntimeError(f"pheasy did not write {fc_file}.") + lasso_alpha = _check_lasso_alpha(log_file, alpha_min, alpha_min_name="alpha_min") + phonon.force_constants = parse_FORCE_CONSTANTS(filename=str(fc_file)) + phonon.symmetrize_force_constants() + + # the fitted forces are F = -Phi u + predicted = -np.einsum("ijab,mjb->mia", phonon.force_constants, disps) + force_rmse = float(np.sqrt(np.mean((fit_forces - predicted) ** 2))) + + borns, epsilon = _set_nac_params(phonon, born, epsilon_static, symprec, code) + + # frequencies at the q-points commensurate with the supercell + matrix = np.linalg.inv(phonon.primitive_matrix) @ phonon.supercell_matrix + phonon.run_qpoints(get_commensurate_points(np.rint(matrix).astype(int))) + frequencies = phonon.qpoints.frequencies + + kpath_dict, kpath_concrete = _get_kpath( + structure=get_pmg_structure(phonon.primitive), + kpath_scheme=_KPATH_SCHEME, + symprec=symprec, + ) + bs_symm_line, has_imaginary_modes = _run_band_structure_and_plot( + phonon, + kpath_dict, + kpath_concrete, + _DEFAULT_FILE_PATHS["band_structure"], + has_nac=phonon.nac_params is not None, + npoints_band=npoints_band, + filename_bs=_DEFAULT_FILE_PATHS["band_structure_plot"], + tol_imaginary_modes=tol_imaginary_modes, + ) + kpoint = Kpoints.automatic_density( + structure=get_pmg_structure(phonon.primitive), + kppa=kpoint_density_dos, + force_gamma=True, + ) + dos = _run_total_dos_and_plot( + phonon, + kpoint, + _DEFAULT_FILE_PATHS["dos"], + filename_dos=_DEFAULT_FILE_PATHS["dos_plot"], + ) + phonon.save(_DEFAULT_FILE_PATHS["phonopy"]) + + return FiniteTemperaturePhononDoc.from_structure( + meta_structure=structure, + structure=structure, + temperature=temperature, + code=code, + md_code=md_code, + force_field_name=force_field_name, + force_field_kwargs=force_field_kwargs, + md_force_field_name=md_force_field_name, + md_force_field_kwargs=md_force_field_kwargs, + thermostat=thermostat, + md_time_step=md_time_step, + md_time=snapshot_data["md_time"], + equilibration_time=equilibration_time, + n_snapshots=n_data, + snapshot_times=snapshot_data["snapshot_times"], + supercell_matrix=phonon.supercell_matrix.tolist(), + primitive_matrix=phonon.primitive_matrix.tolist(), + rotational_sum_rule=rotational_sum_rule, + lasso_alpha=lasso_alpha, + force_rmse=force_rmse, + max_residual_force=float(np.abs(residual_forces).max()), + rms_displacement=snapshot_data["rms_displacement"], + rms_displacement_from_force_constants=_get_rms_displacement( + phonon, temperature, tol_imaginary_modes + ), + trajectory_health=snapshot_data["trajectory_health"], + tol_imaginary_modes=tol_imaginary_modes, + has_imaginary_modes=has_imaginary_modes, + n_imaginary_modes=int(np.sum(frequencies < -tol_imaginary_modes)), + lowest_frequency=float(frequencies.min()), + phonon_bandstructure=bs_symm_line, + phonon_dos=dos, + phonopy_settings=PhononComputationalSettings( + npoints_band=npoints_band, + kpath_scheme=_KPATH_SCHEME, + kpoint_density_dos=kpoint_density_dos, + ), + force_constants=ForceConstants(phonon.force_constants.tolist()) + if store_force_constants + else None, + born=borns.tolist() if borns is not None else None, + epsilon_static=epsilon.tolist() if epsilon is not None else None, + md_uuids=md_uuids, + md_job_dirs=md_job_dirs, + npt_input_structure=npt_input_structure, + pressure=pressure, + npt_time=npt_time, + npt_equilibration_time=npt_equilibration_time, + npt_trajectory_health=npt_trajectory_health, + npt_uuid=npt_uuid, + npt_job_dir=npt_job_dir, + fixed_cell_relax_uuid=fixed_cell_relax_uuid, + fixed_cell_relax_job_dir=fixed_cell_relax_job_dir, + uuids=PhononUUIDs( + optimization_run_uuid=optimization_run_uuid, + displacements_uuids=displacement_data["uuids"], + born_run_uuid=born_run_uuid, + ), + jobdirs=PhononJobDirs( + displacements_job_dirs=displacement_data["dirs"], + born_run_job_dir=born_run_job_dir, + optimization_run_job_dir=optimization_run_job_dir, + taskdoc_run_job_dir=str(Path.cwd()), + ), + ) diff --git a/src/atomate2/common/jobs/md.py b/src/atomate2/common/jobs/md.py new file mode 100644 index 0000000000..dd8a520393 --- /dev/null +++ b/src/atomate2/common/jobs/md.py @@ -0,0 +1,66 @@ +"""Jobs for molecular dynamics that are shared between codes.""" + +from __future__ import annotations + +from pathlib import Path + +from ase.io import read as ase_read +from jobflow import job +from monty.os.path import zpath +from pymatgen.core import Structure +from pymatgen.io.vasp import Poscar + +from atomate2.utils.path import strip_hostname + + +def _get_site_properties(reference: Structure) -> dict: + """Get the magnetic moments of the reference, the only site property kept.""" + if "magmom" in reference.site_properties: + return {"magmom": reference.site_properties["magmom"]} + return {} + + +@job +def get_md_restart_structure( + md_dir: str, reference: Structure, traj_file: str | None = None +) -> Structure: + """ + Get the final positions and velocities of an MD run. + + For VASP, they are read from CONTCAR, with the velocities in Angstrom/fs. + For an ASE MD, they are read from the last frame of its trajectory file, + with the velocities in ASE units. In both cases, the velocities are a + site property that the next MD run starts from. The thermostat variables + are not carried over, so they start again from zero in the next run. The + magnetic moments of the reference, if any, are added as a site property. + + Parameters + ---------- + md_dir: str + Directory of the MD run. + reference: Structure + The structure the MD started from. + traj_file: str | None + Name of the ASE trajectory file of an ASE MD, written in the ASE format. + None for a VASP MD. + + Returns + ------- + Structure + The final structure, with a "velocities" site property. + """ + directory = Path(strip_hostname(md_dir)) + if traj_file is None: + structure = Poscar.from_file(zpath(str(directory / "CONTCAR"))).structure + else: + atoms = ase_read(directory / traj_file, index=-1) + structure = Structure( + atoms.cell[:], + atoms.get_chemical_symbols(), + atoms.get_positions(), + coords_are_cartesian=True, + site_properties={"velocities": atoms.get_velocities().tolist()}, + ) + for key, values in _get_site_properties(reference).items(): + structure.add_site_property(key, values) + return structure diff --git a/src/atomate2/common/jobs/pheasy.py b/src/atomate2/common/jobs/pheasy.py index 3389d8f689..b6984f372b 100644 --- a/src/atomate2/common/jobs/pheasy.py +++ b/src/atomate2/common/jobs/pheasy.py @@ -37,6 +37,7 @@ CubicSupercellTransformation, ) +from atomate2 import SETTINGS from atomate2.common.jobs.phonons import ( ANGSTROM_TO_BOHR, _generate_phonon_object, @@ -74,6 +75,7 @@ "one_shot_dir": "one_shot", "refit_dir": "short_cutoff_refit", "anharmonic_fit_log": "pheasy_anharmonic_fit.log", + "harmonic_fit_log": "pheasy_harmonic_fit.log", } # The anharmonic training set is sized so that the fit has this many force @@ -220,7 +222,9 @@ def _check_anharmonic_settings( ) -def _check_lasso_alpha(log_file: Path, alpha_min: int) -> None: +def _check_lasso_alpha( + log_file: Path, alpha_min: int, alpha_min_name: str = "anhar_alpha_min" +) -> float: """ Warn if the cross-validated LASSO penalty is on either bound of the search. @@ -235,11 +239,18 @@ def _check_lasso_alpha(log_file: Path, alpha_min: int) -> None: Log file written by pheasy during the fit. alpha_min: int Base-10 exponent of the smallest penalty in the search. + alpha_min_name: str + Name of the setting that gives alpha_min, used in the warning. + + Returns + ------- + float + The penalty chosen by cross-validation. """ if not log_file.exists(): raise FileNotFoundError( f"pheasy exited without writing {log_file}, so the result of the " - "anharmonic fit is unknown." + "fit is unknown." ) match = re.search(r"alpha_opt:\s*([-+0-9.eE]+)", log_file.read_text()) if match is None: @@ -258,11 +269,90 @@ def _check_lasso_alpha(log_file: Path, alpha_min: int) -> None: if alpha_opt <= 10.0**alpha_min * (1 + 1e-6): warnings.warn( f"The LASSO penalty chosen by cross-validation ({alpha_opt:e}) is on " - f"the lower bound 1e{alpha_min} in {log_file}. The anharmonic force " - "constants depend on this bound and may be wrong. Lower " - "anhar_alpha_min and refit.", + f"the lower bound 1e{alpha_min} in {log_file}. The force constants " + f"depend on this bound and may be wrong. Lower {alpha_min_name} and " + "refit.", stacklevel=2, ) + return alpha_opt + + +def _run_harmonic_fit( + supercell_matrix: np.ndarray, + symprec: float, + num_har: int, + use_lasso: bool = True, + rotational_sum_rule: str | None = "BHH", + alpha_min: int | None = None, + random_seed: int | None = 103, + log_file: str | None = None, +) -> None: + """ + Fit the second-order force constants with pheasy in the current folder. + + pheasy reads POSCAR, SPOSCAR and the harmonic displacement and force + matrices, and writes FORCE_CONSTANTS. The caller checks the result, since + the pheasy phonon workflow and the finite-temperature workflow check + different things. + + Parameters + ---------- + supercell_matrix: np.ndarray + Diagonal supercell matrix. + symprec: float + Symmetry precision. + num_har: int + Number of displaced supercells in the fit. + use_lasso: bool + If True, fit with LASSO on standardized data, converged to --tol 1e-8. + If False, fit with pheasy's default least squares. + rotational_sum_rule: str | None + Rotational sum rule passed to pheasy with --rasr, or None for none. + alpha_min: int | None + Base-10 exponent of the smallest LASSO penalty in the search. None + keeps pheasy's default. Ignored without LASSO. + random_seed: int | None + Seed for the LASSO fit. Ignored without LASSO. + log_file: str | None + Log file of the fit. None keeps pheasy's default. + """ + dim = " ".join(str(int(supercell_matrix[i][i])) for i in range(3)) + base = ( + f"{SETTINGS.PHEASY_CMD} --scell SPOSCAR --dim {dim} -w 2 " + f"--symprec {float(symprec)}" + ) + fit = f"{base} -f --full_ifc" + if use_lasso: + # pheasy's default --tol of 1e-4 leaves the fit unconverged, so the force + # constants change between machines and package versions + fit += " -l LASSO --std --tol 1e-8" + if alpha_min is not None: + fit += f" --alpha_min {int(alpha_min)}" + if random_seed is not None: + fit += f" --seed {int(random_seed)}" + if rotational_sum_rule is not None: + fit += f" --rasr {rotational_sum_rule}" + fit += ( + f" --ndata {int(num_har)} " + f"--force_matrix_file {_DEFAULT_FILE_PATHS['harmonic_force_matrix']}" + ) + if log_file is not None: + fit += f" -o {log_file}" + + commands = [ + # clusters and orbits of the second-order force constants + f"{base} -s --nbody 2", + # null space + f"{base} -c", + # sensing matrix from the displacement matrix + ( + f"{base} -d --ndata {int(num_har)} --disp_file " + f"--disp_matrix_file {_DEFAULT_FILE_PATHS['harmonic_displacements']}" + ), + fit, + ] + for cmd in commands: + subprocess.run(shlex.split(cmd), check=True) def _run_anharmonic_fit( @@ -315,7 +405,7 @@ def _run_anharmonic_fit( """ dim = " ".join(str(int(supercell_matrix[i][i])) for i in range(3)) base = ( - f"pheasy --scell SPOSCAR --dim {dim} -w {anhar_max_order} " + f"{SETTINGS.PHEASY_CMD} --scell SPOSCAR --dim {dim} -w {anhar_max_order} " f"--symprec {float(symprec)}" ) cutoffs = f"--c3 {float(fcs_cutoff_radius[1] / ANGSTROM_TO_BOHR)}" @@ -728,71 +818,17 @@ def generate_frequencies_eigenvectors( prim = ase_read("POSCAR") supercell = ase_read("SPOSCAR") - # Create the clusters and orbitals for second order force constants. - # The harmonic fit is always second order (-w 2, --nbody 2). - pheasy_cmd_1 = ( - f"pheasy --scell SPOSCAR --dim {int(supercell_matrix[0][0])} " - f"{int(supercell_matrix[1][1])} " - f"{int(supercell_matrix[2][2])} " - f"-s -w 2 --symprec {float(symprec)} --nbody 2" - ) - - # Create the null space to further reduce the free parameters for - # specific force constants and make them physically correct. - pheasy_cmd_2 = ( - f"pheasy --scell SPOSCAR --dim {int(supercell_matrix[0][0])} " - f"{int(supercell_matrix[1][1])} " - f"{int(supercell_matrix[2][2])} -c --symprec " - f"{float(symprec)} -w 2" - ) - - # Generate the Compressive Sensing matrix,i.e., displacement matrix - # for the input of machine leaning method.i.e., LASSO, - pheasy_cmd_3 = ( - f"pheasy --scell SPOSCAR --dim {int(supercell_matrix[0][0])} " - f"{int(supercell_matrix[1][1])} " - f"{int(supercell_matrix[2][2])} -w 2 -d " - f"--symprec {float(symprec)} " - f"--ndata {int(num_har)} --disp_file " - f"--disp_matrix_file {_DEFAULT_FILE_PATHS['harmonic_displacements']}" - ) - - # Here we set a criteria to determine which method to use to generate the - # force constants. If the number of displacements is larger than 3, we - # will use the LASSO method to generate the force constants. Otherwise, - # we will use the least-squred method to generate the force constants. - if len(phonon.displacements) > 3: - # Calculate the force constants using the LASSO method due to the - # random-displacement method Obviously, the rotaional invariance - # constraint, i.e., tag: --rasr BHH, is enforced during the - # fitting process. - pheasy_cmd_4 = ( - f"pheasy --scell SPOSCAR --dim {int(supercell_matrix[0][0])} " - f"{int(supercell_matrix[1][1])} " - f"{int(supercell_matrix[2][2])} -f --full_ifc " - f"-w 2 --symprec {float(symprec)} " - f"-l LASSO --std --rasr BHH --ndata {int(num_har)} " - f"--force_matrix_file {_DEFAULT_FILE_PATHS['harmonic_force_matrix']}" - + (f" --seed {int(random_seed)}" if random_seed is not None else "") - ) - - else: - # Calculate the force constants using the least-squred method - pheasy_cmd_4 = ( - f"pheasy --scell SPOSCAR --dim {int(supercell_matrix[0][0])} " - f"{int(supercell_matrix[1][1])} " - f"{int(supercell_matrix[2][2])} -f --full_ifc " - f"-w 2 --symprec {float(symprec)} " - f"--rasr BHH --ndata {int(num_har)} " - f"--force_matrix_file {_DEFAULT_FILE_PATHS['harmonic_force_matrix']}" - ) - + # With more than 3 displacements, the random-displacement data are fitted + # with LASSO. Otherwise, least squares is used. The rotational sum rule + # BHH is enforced in both. logger.info("Start running pheasy in cluster") - - subprocess.call(shlex.split(pheasy_cmd_1)) - subprocess.call(shlex.split(pheasy_cmd_2)) - subprocess.call(shlex.split(pheasy_cmd_3)) - subprocess.call(shlex.split(pheasy_cmd_4)) + _run_harmonic_fit( + supercell_matrix, + symprec, + num_har, + use_lasso=len(phonon.displacements) > 3, + random_seed=random_seed, + ) fc_file = Path(_DEFAULT_FILE_PATHS["force_constants"]) if cal_anhar_fcs and not fc_file.exists(): @@ -946,7 +982,8 @@ def generate_frequencies_eigenvectors( shutil.copy(filename, refit_dir / filename) pheasy_cmd_11 = ( - f"pheasy --scell SPOSCAR --dim {int(supercell_matrix[0][0])} " + f"{SETTINGS.PHEASY_CMD} --scell SPOSCAR " + f"--dim {int(supercell_matrix[0][0])} " f"{int(supercell_matrix[1][1])} " f"{int(supercell_matrix[2][2])} -s -w 2 --c2 " f"10.0 --symprec {float(symprec)} " @@ -954,14 +991,16 @@ def generate_frequencies_eigenvectors( ) pheasy_cmd_12 = ( - f"pheasy --scell SPOSCAR --dim {int(supercell_matrix[0][0])} " + f"{SETTINGS.PHEASY_CMD} --scell SPOSCAR " + f"--dim {int(supercell_matrix[0][0])} " f"{int(supercell_matrix[1][1])} " f"{int(supercell_matrix[2][2])} -c --symprec " f"{float(symprec)} --c2 10.0 -w 2" ) pheasy_cmd_13 = ( - f"pheasy --scell SPOSCAR --dim {int(supercell_matrix[0][0])} " + f"{SETTINGS.PHEASY_CMD} --scell SPOSCAR " + f"--dim {int(supercell_matrix[0][0])} " f"{int(supercell_matrix[1][1])} " f"{int(supercell_matrix[2][2])} -w 2 -d --symprec " f"{float(symprec)} --c2 10.0 " @@ -973,17 +1012,19 @@ def generate_frequencies_eigenvectors( if len(phonon.displacements) > 3: pheasy_cmd_14 = ( - f"pheasy --scell SPOSCAR --dim {int(supercell_matrix[0][0])} " + f"{SETTINGS.PHEASY_CMD} --scell SPOSCAR " + f"--dim {int(supercell_matrix[0][0])} " f"{int(supercell_matrix[1][1])} " f"{int(supercell_matrix[2][2])} -f --c2 10.0 " f"--full_ifc -w 2 --symprec {float(symprec)} " - f"-l LASSO --std --rasr BHH --ndata {int(num_har)} " + f"-l LASSO --std --tol 1e-8 --rasr BHH --ndata {int(num_har)} " f"--force_matrix_file {_DEFAULT_FILE_PATHS['harmonic_force_matrix']}" ) else: pheasy_cmd_14 = ( - f"pheasy --scell SPOSCAR --dim {int(supercell_matrix[0][0])} " + f"{SETTINGS.PHEASY_CMD} --scell SPOSCAR " + f"--dim {int(supercell_matrix[0][0])} " f"{int(supercell_matrix[1][1])} " f"{int(supercell_matrix[2][2])} -f --full_ifc " f"--c2 10.0 -w 2 --symprec {float(symprec)} " diff --git a/src/atomate2/common/schemas/finite_temperature_phonons.py b/src/atomate2/common/schemas/finite_temperature_phonons.py new file mode 100644 index 0000000000..85482155fc --- /dev/null +++ b/src/atomate2/common/schemas/finite_temperature_phonons.py @@ -0,0 +1,214 @@ +"""Schemas for the finite-temperature phonon workflow.""" + +from typing import Literal + +from pydantic import BaseModel, Field +from pymatgen.core import Structure + +from atomate2.common.schemas.phonons import PhononBSDOSDoc + + +class TrajectoryHealth(BaseModel): + """Check whether the MD trajectory stayed at the reference structure. + + A second-order fit describes small vibrations around the reference + structure. If the atoms melt, move to another structure or diffuse, the + fitted force constants may not describe small vibrations around the + reference. The fit error does not show this, so it is measured on the + trajectory. The displacement of the center of mass of each frame is + removed from its displacements. + """ + + verdict: ( + Literal[ + "stable", + "melted", + "transformed", + "shifted_or_diffusing", + "disordering", + "diffusing_or_soft", + "reference_mismatch", + ] + | None + ) = Field( + None, + description="Result of the check. 'reference_mismatch' means the atoms of " + "the MD are not in the order of the reference. Every other verdict except " + "'stable' means the trajectory may have left the reference structure.", + ) + n_frames: int | None = Field(None, description="Number of MD frames analysed.") + rms_displacement_start: float | None = Field( + None, + description="Root mean square displacement from the reference structure in " + "the first 10 frames, in Angstrom. The MD starts from the reference, so a " + "large value points at a mismatch in atom order.", + ) + rms_displacement_end: float | None = Field( + None, + description="Root mean square displacement from the reference structure in " + "the last fifth of the frames, in Angstrom.", + ) + u_ref: float | None = Field( + None, + description="Root mean square displacement from the reference structure in " + "the second half of the frames, in Angstrom.", + ) + u_shift: float | None = Field( + None, + description="Root mean square displacement of the mean atomic positions from " + "the reference structure in the second half of the frames, in Angstrom.", + ) + u_vib: float | None = Field( + None, + description="Root mean square vibration around the mean atomic positions in " + "the second half of the frames, in Angstrom. u_ref**2 = u_shift**2 + " + "u_vib**2.", + ) + nearest_neighbor_distance: float | None = Field( + None, + description="Nearest-neighbor distance of each atom in the reference " + "structure, averaged over the atoms, in Angstrom.", + ) + lindemann_ratio: float | None = Field( + None, description="u_vib divided by the nearest-neighbor distance." + ) + shift_ratio: float | None = Field(None, description="u_ref divided by u_vib.") + energy_drift: float | None = Field( + None, + description="Mean potential energy in the last quarter of the frames minus " + "that in the second quarter, in eV/atom. None for fewer than 20 frames.", + ) + energy_fluctuation: float | None = Field( + None, + description="Standard deviation of the potential energy in the last fifth of " + "the frames, in eV/atom.", + ) + + +class FiniteTemperaturePhononDoc(PhononBSDOSDoc): + """Effective harmonic phonons at a finite temperature. + + structure is the unit cell of the NVT MD and the fit, with its atoms sorted + by electronegativity. It is the relaxed unit cell, or the cell from the NPT + MD if there was one. code is the code of the phonon displacement + calculations. born and epsilon_static are the symmetrized values used for + the non-analytical correction. has_imaginary_modes is True if the band + structure has imaginary modes. The undisplaced supercell is the last + phonon displacement calculation in uuids and jobdirs. The thermodynamic + properties are not computed. + """ + + temperature: float | None = Field(None, description="MD temperature in K.") + md_code: str | None = Field( + None, description="Code of the MD, 'vasp' or 'forcefields'." + ) + force_field_name: str | None = Field( + None, + description="Force field of the phonon displacement calculations, if a " + "force field was used.", + ) + force_field_kwargs: dict | None = Field( + None, + description="Keyword arguments of the force field calculator of the phonon " + "displacement calculations.", + ) + md_force_field_name: str | None = Field( + None, description="Force field of the MD, if a force field was used." + ) + md_force_field_kwargs: dict | None = Field( + None, description="Keyword arguments of the force field calculator of the MD." + ) + thermostat: str | None = Field(None, description="Thermostat of the NVT MD.") + md_time_step: float | None = Field(None, description="MD time step in fs.") + md_time: float | None = Field( + None, + description="Length of the joined MD trajectory, the number of frames times " + "the time step, in ps.", + ) + equilibration_time: float | None = Field( + None, description="Time at the start of the trajectory left out, in ps." + ) + n_snapshots: int | None = Field( + None, description="Number of MD snapshots in the force constant fit." + ) + snapshot_times: list[float] | None = Field( + None, + description="MD time of each snapshot in ps, its step number times the time " + "step.", + ) + rotational_sum_rule: str | None = Field( + None, description="Rotational sum rule passed to pheasy with --rasr." + ) + lasso_alpha: float | None = Field( + None, description="LASSO penalty chosen by cross-validation in pheasy." + ) + force_rmse: float | None = Field( + None, + description="Root mean square error of the symmetrized force constants on " + "the snapshot forces used in the fit, in eV/Angstrom.", + ) + max_residual_force: float | None = Field( + None, + description="Largest force component on the undisplaced supercell, in " + "eV/Angstrom. The residual forces of this supercell are subtracted from the " + "snapshot forces before the fit.", + ) + rms_displacement: float | None = Field( + None, + description="Root mean square displacement of the snapshots from the " + "reference structure, with the displacement of the center of mass of each " + "snapshot removed, in Angstrom.", + ) + rms_displacement_from_force_constants: float | None = Field( + None, + description="Classical root mean square displacement at the temperature from " + "the fitted force constants, in Angstrom. It sums the modes of the " + "supercell at the Gamma point with a frequency above tol_imaginary_modes.", + ) + trajectory_health: TrajectoryHealth | None = Field( + None, description="Check of the MD trajectory." + ) + tol_imaginary_modes: float | None = Field( + None, description="Frequencies below -tol_imaginary_modes in THz are imaginary." + ) + n_imaginary_modes: int | None = Field( + None, + description="Number of frequencies below -tol_imaginary_modes, summed over " + "all q-points commensurate with the supercell. q-points that are equivalent " + "by symmetry are each counted.", + ) + lowest_frequency: float | None = Field( + None, + description="Lowest frequency at the q-points commensurate with the " + "supercell, in THz. Imaginary frequencies are negative.", + ) + md_uuids: list[str] | None = Field(None, description="UUIDs of the MD jobs.") + md_job_dirs: list[str] | None = Field( + None, description="Directories of the MD jobs." + ) + npt_input_structure: Structure | None = Field( + None, + description="Unit cell whose supercell the NPT MD started from. None without " + "an NPT MD.", + ) + pressure: float | None = Field(None, description="Pressure of the NPT MD in kbar.") + npt_time: float | None = Field(None, description="Length of the NPT MD in ps.") + npt_equilibration_time: float | None = Field( + None, + description="Time at the start of the NPT trajectory left out of the " + "average cell, in ps.", + ) + npt_trajectory_health: TrajectoryHealth | None = Field( + None, description="Check of the NPT trajectory." + ) + npt_uuid: str | None = Field(None, description="UUID of the NPT MD job.") + npt_job_dir: str | None = Field(None, description="Directory of the NPT MD job.") + fixed_cell_relax_uuid: str | None = Field( + None, + description="UUID of the relaxation of the atoms in the cell from the NPT MD.", + ) + fixed_cell_relax_job_dir: str | None = Field( + None, + description="Directory of the relaxation of the atoms in the cell from the " + "NPT MD.", + ) diff --git a/src/atomate2/forcefields/flows/finite_temperature_phonons.py b/src/atomate2/forcefields/flows/finite_temperature_phonons.py new file mode 100644 index 0000000000..8709dc67af --- /dev/null +++ b/src/atomate2/forcefields/flows/finite_temperature_phonons.py @@ -0,0 +1,519 @@ +"""Define the force field makers for finite-temperature phonons.""" + +from __future__ import annotations + +import math +from dataclasses import dataclass, field, replace +from typing import TYPE_CHECKING + +from ase import units + +from atomate2.ase.md import MDEnsemble +from atomate2.common.flows.finite_temperature_phonons import ( + BaseFiniteTemperaturePhononMaker, +) +from atomate2.common.jobs.finite_temperature_phonons import ASE_TRAJECTORY_FILE +from atomate2.forcefields.jobs import ForceFieldRelaxMaker, ForceFieldStaticMaker +from atomate2.forcefields.md import ForceFieldMDMaker +from atomate2.vasp.flows.finite_temperature_phonons import FiniteTemperaturePhononMaker + +if TYPE_CHECKING: + from jobflow import Maker + from typing_extensions import Self + + from atomate2.forcefields import MLFF + +_DEFAULT_FORCE_FIELD = "MACE-MP-0" + +# VASP sets the Nose mass for SMASS = 0 so that the temperature oscillates with a +# period of about 40 time steps. With ASE's thermostat mass Q = 3 N k_B T tdamp**2, the +# linearised Nose-Hoover period is close to pi * sqrt(2) * tdamp. +_NOSE_HOOVER_TDAMP_STEPS = 40 / (math.pi * math.sqrt(2)) +# time constant of the barostat of the NPT MD +_BAROSTAT_PDAMP_STEPS = 1000 + + +def _get_md_maker( + force_field_name: str | MLFF | dict, calculator_kwargs: dict | None = None +) -> ForceFieldMDMaker: + """Get the force field MD maker, whose settings the flow sets.""" + return ForceFieldMDMaker( + force_field_name=force_field_name, + calculator_kwargs=calculator_kwargs or {}, + store_trajectory="no", + ionic_step_data=("energy",), + ) + + +def _get_force_field_md_maker( + maker: BaseFiniteTemperaturePhononMaker, + n_steps: int, + index: int = 0, + npt: bool = False, +) -> ForceFieldMDMaker: + """ + Get the force field MD maker of one MD job of a finite-temperature flow. + + See :obj:`ForceFieldFiniteTemperaturePhononMaker` for the MD settings. + + Parameters + ---------- + maker: .BaseFiniteTemperaturePhononMaker + The finite-temperature phonon maker. + n_steps: int + Number of MD steps of the job. + index: int + Position of the job in the MD, counted from 0. It is added to + maker.random_seed. The NPT MD uses maker.md_runs instead. + npt: bool + If True, get the NPT MD maker from maker.npt_maker. If False, get the + NVT MD maker from maker.md_maker. + + Returns + ------- + ForceFieldMDMaker + """ + md_maker = maker.npt_maker if npt else maker.md_maker + name = "NPT maker" if npt else "MD maker" + if not isinstance(md_maker, ForceFieldMDMaker): + raise TypeError( + f"The {name} must be a ForceFieldMDMaker, not {type(md_maker).__name__}." + ) + if md_maker.dynamics is not None or md_maker.ase_md_kwargs: + raise ValueError( + f"The flow sets the thermostat. The {name} must not set dynamics or " + "ase_md_kwargs." + ) + if npt or maker.thermostat == "nose-hoover": + dynamics = "nose-hoover-chain" + md_kwargs = { + "tdamp": _NOSE_HOOVER_TDAMP_STEPS * maker.md_time_step * units.fs, + "tchain": 1, + } + if npt: + md_kwargs["pdamp"] = _BAROSTAT_PDAMP_STEPS * maker.md_time_step * units.fs + else: + # AseMDMaker sets its default friction of 10 ps^-1 + dynamics = "langevin" + md_kwargs = {} + return replace( + md_maker, + ensemble=MDEnsemble.npt if npt else MDEnsemble.nvt, + temperature=maker.temperature, + pressure=maker.pressure if npt else md_maker.pressure, + n_steps=n_steps, + time_step=maker.md_time_step, + dynamics=dynamics, + ase_md_kwargs=md_kwargs, + traj_file=ASE_TRAJECTORY_FILE, + traj_file_fmt="ase", + traj_interval=1, + store_trajectory="no", + mb_velocity_seed=( + None + if maker.random_seed is None + else maker.random_seed + (maker.md_runs if npt else index) + ), + zero_linear_momentum=True, + ) + + +@dataclass +class ForceFieldFiniteTemperaturePhononMaker(BaseFiniteTemperaturePhononMaker): + """ + Maker for effective harmonic phonons at a finite temperature with a force field. + + The relaxation, the MD and the phonon displacement calculations all use one + force field, MACE-MP-0 by default. Use :obj:`from_force_field_name` to use + another one. The MD writes an ASE trajectory file with a frame at every + time step. The snapshots are read from this file, and the trajectory is not + stored in the output document of the MD job. The Langevin thermostat has + the default friction of :obj:`.AseMDMaker`, 10 ps^-1. The Nose-Hoover + thermostat is ASE's NoseHooverChainNVT with one thermostat variable, like + MDALGO = 2 in VASP. Its time constant gives a period of about 40 time + steps, like SMASS = 0 in VASP. The initial velocities of the first MD job + follow the Maxwell-Boltzmann distribution, with zero total momentum. The + NPT MD, if any, uses ASE's MTKNPT, a Nose-Hoover thermostat and barostat + that change the whole cell, whatever the thermostat setting (Martyna et + al., J. Chem. Phys. 101, 4177 (1994)). Its thermostat time constant is + that of the Nose-Hoover thermostat above, and its barostat time constant + is 1000 time steps. MD job i, counted from 0, seeds its initial velocities + and its Langevin random forces with random_seed + i. The NPT MD uses + random_seed + md_runs. + + See :obj:`.BaseFiniteTemperaturePhononMaker` for the workflow. + + This workflow is new and has not been tested widely. It might still change + in future versions. + + Parameters + ---------- + name: str + Name of the flows produced by this maker. + bulk_relax_maker: .ForceFieldRelaxMaker | None + Maker for the relaxation of the unit cell. It keeps the symmetry of + the structure. None skips the relaxation. + born_maker: .Maker | None + Maker for the Born effective charges and the dielectric tensor, for + example a ForceFieldDielectricMaker with MACE-Field or a VASP + DielectricMaker. It is None by default, as in the force field pheasy + phonon workflow. + npt_maker: .ForceFieldMDMaker | None + Maker for the NPT MD. The flow sets the same settings as for md_maker, + and the pressure. It must not set dynamics or ase_md_kwargs. None skips + the NPT MD. + fixed_cell_relax_maker: .ForceFieldRelaxMaker | None + Maker for the relaxation of the atoms in the cell from the NPT MD. It + must keep the cell. + md_maker: .ForceFieldMDMaker + Maker for the MD. The flow sets its ensemble, temperature, n_steps, + time_step, dynamics, ase_md_kwargs, traj_file, traj_file_fmt, + traj_interval, store_trajectory, mb_velocity_seed and + zero_linear_momentum. It must not set dynamics or ase_md_kwargs. + phonon_displacement_maker: .ForceFieldStaticMaker + Maker for the static calculations on the snapshots and the undisplaced + supercell. + code: str + Code of the phonon displacement calculations. + md_code: str + Code of the MD. + """ + + name: str = "force field finite temperature phonon" + bulk_relax_maker: ForceFieldRelaxMaker | None = field( + default_factory=lambda: ForceFieldRelaxMaker( + force_field_name=_DEFAULT_FORCE_FIELD, + relax_kwargs={"fmax": 1e-5}, + fix_symmetry=True, + ) + ) + md_maker: ForceFieldMDMaker = field( + default_factory=lambda: _get_md_maker(_DEFAULT_FORCE_FIELD) + ) + phonon_displacement_maker: ForceFieldStaticMaker = field( + default_factory=lambda: ForceFieldStaticMaker( + force_field_name=_DEFAULT_FORCE_FIELD + ) + ) + code: str = "forcefields" + md_code: str = "forcefields" + + @property + def prev_calc_dir_argname(self) -> None: + """Name of the prev_dir argument of the phonon displacement and Born makers. + + Force field makers take no previous directory, so it is None. A VASP + born_maker then runs without a previous directory. + """ + return + + def get_md_maker(self, n_steps: int, index: int = 0) -> ForceFieldMDMaker: + """ + Get the force field MD maker of one MD job. + + Parameters + ---------- + n_steps: int + Number of MD steps of the job. + index: int + Position of the job in the MD, counted from 0. + + Returns + ------- + ForceFieldMDMaker + """ + return _get_force_field_md_maker(self, n_steps, index) + + def get_npt_maker(self, n_steps: int) -> ForceFieldMDMaker: + """ + Get the force field maker of the NPT MD job. + + Parameters + ---------- + n_steps: int + Number of MD steps. + + Returns + ------- + ForceFieldMDMaker + """ + return _get_force_field_md_maker(self, n_steps, npt=True) + + @classmethod + def from_force_field_name( + cls, + force_field_name: str | MLFF | dict, + calculator_kwargs: dict | None = None, + relax_initial_structure: bool = True, + run_npt: bool = False, + **kwargs, + ) -> Self: + """ + Create a finite-temperature phonon flow from a force field name. + + The relaxation, the MD and the phonon displacement calculations use the + force field, and born_maker is None. These makers replace any given in + kwargs. + + Parameters + ---------- + force_field_name : str or .MLFF or dict + The name of the force field. + calculator_kwargs : dict | None + The keyword arguments to pass to the calculator. + relax_initial_structure: bool = True + Whether to relax the initial structure. + run_npt: bool = False + Whether to run an NPT MD with the force field first. The NVT MD then + uses its average cell, with the atoms relaxed in that cell. + **kwargs + Additional kwargs to pass to ForceFieldFiniteTemperaturePhononMaker. + + Returns + ------- + ForceFieldFiniteTemperaturePhononMaker + """ + calculator_kwargs = calculator_kwargs or {} + phonon_displacement_maker = ForceFieldStaticMaker( + force_field_name=force_field_name, + calculator_kwargs=calculator_kwargs, + ) + npt_maker: ForceFieldMDMaker | None = None + fixed_cell_relax_maker: ForceFieldRelaxMaker | None = None + if run_npt: + npt_maker = _get_md_maker(force_field_name, calculator_kwargs) + fixed_cell_relax_maker = ForceFieldRelaxMaker( + force_field_name=force_field_name, + calculator_kwargs=calculator_kwargs, + relax_cell=False, + relax_kwargs={"fmax": 1e-5}, + fix_symmetry=True, + ) + kwargs.update( + npt_maker=npt_maker, + fixed_cell_relax_maker=fixed_cell_relax_maker, + bulk_relax_maker=( + ForceFieldRelaxMaker( + force_field_name=force_field_name, + calculator_kwargs=calculator_kwargs, + relax_kwargs={"fmax": 1e-5}, + fix_symmetry=True, + ) + if relax_initial_structure + else None + ), + md_maker=_get_md_maker(force_field_name, calculator_kwargs), + phonon_displacement_maker=phonon_displacement_maker, + born_maker=None, + ) + return cls( + name=( + f"{phonon_displacement_maker.mlff.name} Finite Temperature Phonon Maker" + ), + **kwargs, + ) + + +@dataclass +class VaspMDMLFFStaticFiniteTemperaturePhononMaker(FiniteTemperaturePhononMaker): + """ + Maker for finite-temperature phonons from VASP MD and force field statics. + + The relaxation and the MD are those of + :obj:`.FiniteTemperaturePhononMaker`. The forces on the snapshots and on + the undisplaced supercell come from a force field, MACE-MP-0 by default. + Use :obj:`from_force_field_name` to use another one. + + This workflow is new and has not been tested widely. It might still change + in future versions. + + Parameters + ---------- + name: str + Name of the flows produced by this maker. + bulk_relax_maker: .Maker | None + Maker for the relaxation of the unit cell. None skips the relaxation. + born_maker: .Maker | None + Maker for the Born effective charges and the dielectric tensor. It is + None by default, as in the force field pheasy phonon workflow, since + the force constants come from the force field. + md_maker: .MDMaker + Maker for the MD, see :obj:`.FiniteTemperaturePhononMaker`. + phonon_displacement_maker: .ForceFieldStaticMaker + Maker for the static calculations on the snapshots and the undisplaced + supercell. + code: str + Code of the phonon displacement calculations. + """ + + name: str = "vasp md mlff static finite temperature phonon" + born_maker: Maker | None = None + phonon_displacement_maker: Maker = field( + default_factory=lambda: ForceFieldStaticMaker( + force_field_name=_DEFAULT_FORCE_FIELD + ) + ) + code: str = "forcefields" + + @property + def prev_calc_dir_argname(self) -> None: + """Name of the prev_dir argument of the phonon displacement and Born makers. + + Force field makers take no previous directory, so it is None. A VASP + born_maker then runs without a previous directory. + """ + return + + @classmethod + def from_force_field_name( + cls, + force_field_name: str | MLFF | dict, + calculator_kwargs: dict | None = None, + **kwargs, + ) -> Self: + """ + Create a finite-temperature phonon flow whose statics use a force field. + + The phonon displacement calculations use the force field, and born_maker + is None. These makers replace any given in kwargs. + + Parameters + ---------- + force_field_name : str or .MLFF or dict + The name of the force field. + calculator_kwargs : dict | None + The keyword arguments to pass to the calculator. + **kwargs + Additional kwargs to pass to + VaspMDMLFFStaticFiniteTemperaturePhononMaker. + + Returns + ------- + VaspMDMLFFStaticFiniteTemperaturePhononMaker + """ + phonon_displacement_maker = ForceFieldStaticMaker( + force_field_name=force_field_name, + calculator_kwargs=calculator_kwargs or {}, + ) + kwargs.update( + phonon_displacement_maker=phonon_displacement_maker, born_maker=None + ) + return cls( + name=( + f"VASP MD {phonon_displacement_maker.mlff.name} Static Finite " + "Temperature Phonon Maker" + ), + **kwargs, + ) + + +@dataclass +class MLFFMDVaspStaticFiniteTemperaturePhononMaker(FiniteTemperaturePhononMaker): + """ + Maker for finite-temperature phonons from force field MD and VASP statics. + + The relaxation, the Born charges and the phonon displacement calculations + are those of :obj:`.FiniteTemperaturePhononMaker`. The MD uses a force + field, MACE-MP-0 by default, and starts from the supercell of the + VASP-relaxed structure. Use :obj:`from_force_field_name` to use another + force field. The MD settings are those of + :obj:`ForceFieldFiniteTemperaturePhononMaker`. + + This workflow is new and has not been tested widely. It might still change + in future versions. + + Parameters + ---------- + name: str + Name of the flows produced by this maker. + bulk_relax_maker: .Maker | None + Maker for the relaxation of the unit cell. None skips the relaxation. + born_maker: .BaseVaspMaker | None + Maker for the Born effective charges and the dielectric tensor, used + for the non-analytical correction. It is a DielectricMaker by default, + as in the VASP pheasy phonon workflow. None skips it. + md_maker: .ForceFieldMDMaker + Maker for the MD. The flow sets its ensemble, temperature, n_steps, + time_step, dynamics, ase_md_kwargs, traj_file, traj_file_fmt, + traj_interval, store_trajectory, mb_velocity_seed and + zero_linear_momentum. It must not set dynamics or ase_md_kwargs. + npt_maker: .ForceFieldMDMaker | None + Maker for the NPT MD, under the same conditions as md_maker. None skips + the NPT MD. + phonon_displacement_maker: .BaseVaspMaker + Maker for the static calculations on the snapshots and the undisplaced + supercell. + md_code: str + Code of the MD. + """ + + name: str = "mlff md vasp static finite temperature phonon" + md_maker: Maker = field(default_factory=lambda: _get_md_maker(_DEFAULT_FORCE_FIELD)) + md_code: str = "forcefields" + + def get_md_maker(self, n_steps: int, index: int = 0) -> ForceFieldMDMaker: + """ + Get the force field MD maker of one MD job. + + Parameters + ---------- + n_steps: int + Number of MD steps of the job. + index: int + Position of the job in the MD, counted from 0. + + Returns + ------- + ForceFieldMDMaker + """ + return _get_force_field_md_maker(self, n_steps, index) + + def get_npt_maker(self, n_steps: int) -> ForceFieldMDMaker: + """ + Get the force field maker of the NPT MD job. + + Parameters + ---------- + n_steps: int + Number of MD steps. + + Returns + ------- + ForceFieldMDMaker + """ + return _get_force_field_md_maker(self, n_steps, npt=True) + + @classmethod + def from_force_field_name( + cls, + force_field_name: str | MLFF | dict, + calculator_kwargs: dict | None = None, + **kwargs, + ) -> Self: + """ + Create a finite-temperature phonon flow whose MD uses a force field. + + The MD maker built for the force field replaces any given in kwargs. + + Parameters + ---------- + force_field_name : str or .MLFF or dict + The name of the force field. + calculator_kwargs : dict | None + The keyword arguments to pass to the calculator. + **kwargs + Additional kwargs to pass to + MLFFMDVaspStaticFiniteTemperaturePhononMaker. + + Returns + ------- + MLFFMDVaspStaticFiniteTemperaturePhononMaker + """ + md_maker = _get_md_maker(force_field_name, calculator_kwargs) + kwargs.update(md_maker=md_maker) + return cls( + name=( + f"{md_maker.mlff.name} MD VASP Static Finite Temperature Phonon Maker" + ), + **kwargs, + ) diff --git a/src/atomate2/settings.py b/src/atomate2/settings.py index 32bcbf2da2..99c21f13c1 100644 --- a/src/atomate2/settings.py +++ b/src/atomate2/settings.py @@ -258,6 +258,8 @@ class Atomate2Settings(BaseSettings): "GBRV_v1.5", description="location of JDFTX pseudopotentials." ) + PHEASY_CMD: str = Field("pheasy", description="The command to run pheasy.") + LAMMPS_CMD: str = Field("lmp", description="The command to run LAMMPS.") LAMMPS_MPICMD: str | None = Field( diff --git a/src/atomate2/vasp/flows/finite_temperature_phonons.py b/src/atomate2/vasp/flows/finite_temperature_phonons.py new file mode 100644 index 0000000000..8bf78746a8 --- /dev/null +++ b/src/atomate2/vasp/flows/finite_temperature_phonons.py @@ -0,0 +1,267 @@ +"""Define the VASP maker for finite-temperature phonons.""" + +from __future__ import annotations + +from dataclasses import dataclass, field, fields, replace +from typing import TYPE_CHECKING + +from atomate2.common.flows.finite_temperature_phonons import ( + BaseFiniteTemperaturePhononMaker, +) +from atomate2.vasp.flows.core import DoubleRelaxMaker +from atomate2.vasp.jobs.core import DielectricMaker, TightRelaxMaker +from atomate2.vasp.jobs.md import MDMaker +from atomate2.vasp.jobs.phonons import PhononDisplacementMaker +from atomate2.vasp.sets.core import ( + LangevinMDSetGenerator, + MDSetGenerator, + StaticSetGenerator, + TightRelaxSetGenerator, +) + +if TYPE_CHECKING: + from jobflow import Maker + + from atomate2.vasp.jobs.base import BaseVaspMaker + +# INCAR tags of the MD that the flow sets +_MD_FLOW_TAGS = ( + "IBRION", + "ISIF", + "LANGEVIN_GAMMA", + "MDALGO", + "NSW", + "POTIM", + "SMASS", + "TEBEG", + "TEEND", +) + + +def _check_md_incar(user_incar: dict, name: str) -> None: + """Check that an MD maker leaves the INCAR tags that the flow sets alone.""" + if tags := sorted(set(_MD_FLOW_TAGS) & set(user_incar)): + raise ValueError(f"The flow sets {tags}. Remove them from the {name}.") + if int(user_incar.get("NBLOCK", 1)) != 1: + raise ValueError( + "The MD must write every step to XDATCAR, so NBLOCK must be 1." + ) + + +@dataclass +class FiniteTemperaturePhononMaker(BaseFiniteTemperaturePhononMaker): + """ + Maker for effective harmonic phonons at a finite temperature with VASP. + + The relaxation, the MD and the phonon displacement calculations use VASP. + The functional, the spin settings and +U follow the atomate2 VASP defaults. + All three use ISMEAR = 0 with SIGMA = 0.05 and the same reciprocal_density + setting of 100. The relaxation is a double tight relaxation with ENCUT = + 600 eV, ENAUG = 1360 eV and PREC = Accurate, as in the phonon displacement + calculations. These otherwise keep the settings of PhononDisplacementMaker, + such as EDIFF = 1e-7. The MD uses ENCUT = 500 eV, EDIFF = 1e-5 and PREC = + Normal, and writes every step to XDATCAR with NBLOCK = 1. The Langevin + thermostat uses MDALGO = 3 with LANGEVIN_GAMMA = 10 ps^-1 for each species. + The Nose-Hoover thermostat uses MDALGO = 2 and SMASS = 0. VASP draws the + initial velocities of the first MD job. The MD and the phonon displacement + calculations start from the magnetic moments of the relaxed structure, if + it has any. They all get the same prev_dir, the relaxation directory by + default. auto_ispin of the MD and of the phonon displacement calculations + sets ISPIN from it. + + See :obj:`.BaseFiniteTemperaturePhononMaker` for the workflow. + + This workflow is new and has not been tested widely. It might still change + in future versions. + + Parameters + ---------- + name: str + Name of the flows produced by this maker. + bulk_relax_maker: .Maker | None + Maker for the relaxation of the unit cell. None skips the relaxation. + born_maker: .BaseVaspMaker | None + Maker for the Born effective charges and the dielectric tensor, used + for the non-analytical correction. It is a DielectricMaker by default, + as in the VASP pheasy phonon workflow. None skips it. + md_maker: .MDMaker + Maker for the MD, with an :obj:`.MDSetGenerator`. Its + user_incar_settings must not contain IBRION, ISIF, LANGEVIN_GAMMA, + MDALGO, NSW, POTIM, SMASS, TEBEG or TEEND, since the flow sets them. + NBLOCK must be 1 if it is set. + npt_maker: .MDMaker | None + Maker for the NPT MD, with an :obj:`.MDSetGenerator`, under the same + conditions as md_maker. It must not set PSTRESS either. Its ENCUT should + be larger than that of the NVT MD, since the cell changes. None skips + the NPT MD. + fixed_cell_relax_maker: .Maker | None + Maker for the relaxation of the atoms in the cell from the NPT MD, for + example with ISIF = 2. + phonon_displacement_maker: .BaseVaspMaker + Maker for the static calculations on the snapshots and the undisplaced + supercell. + code: str + Code of the phonon displacement calculations. + md_code: str + Code of the MD. + """ + + bulk_relax_maker: Maker | None = field( + default_factory=lambda: DoubleRelaxMaker.from_relax_maker( + TightRelaxMaker( + input_set_generator=TightRelaxSetGenerator( + user_incar_settings={ + "ENCUT": 600, + "ENAUG": 1360, + "PREC": "Accurate", + "LASPH": True, + "ISMEAR": 0, + "SIGMA": 0.05, + }, + user_kpoints_settings={"reciprocal_density": 100}, + ) + ) + ) + ) + born_maker: BaseVaspMaker | None = field(default_factory=DielectricMaker) + md_maker: Maker = field( + default_factory=lambda: MDMaker( + input_set_generator=MDSetGenerator( + user_incar_settings={ + "ENCUT": 500, + "EDIFF": 1e-5, + "ALGO": "Normal", + "LREAL": "Auto", + "NBLOCK": 1, + "LWAVE": False, + "ISMEAR": 0, + "SIGMA": 0.05, + }, + user_kpoints_settings={"reciprocal_density": 100}, + ) + ) + ) + phonon_displacement_maker: Maker = field( + default_factory=lambda: PhononDisplacementMaker( + input_set_generator=StaticSetGenerator( + user_incar_settings={ + "IBRION": 2, + "ISIF": 3, + "NSW": 0, + "ENCUT": 600, + "ENAUG": 1360, + "EDIFF": 1e-7, + "PREC": "Accurate", + "ALGO": "Normal", + "LASPH": True, + "NELM": 200, + "LREAL": False, + "LAECHG": False, + "LCHARG": False, + "LWAVE": False, + "ISMEAR": 0, + "SIGMA": 0.05, + }, + user_kpoints_settings={"reciprocal_density": 100}, + auto_ispin=True, + ) + ) + ) + code: str = "vasp" + md_code: str = "vasp" + + @property + def prev_calc_dir_argname(self) -> str | None: + """Name of the prev_dir argument of the phonon displacement and Born makers. + + Returns + ------- + str + """ + return "prev_dir" + + def get_md_maker(self, n_steps: int, index: int = 0) -> MDMaker: + """ + Get the VASP MD maker of one MD job. + + Parameters + ---------- + n_steps: int + Number of MD steps of the job. + index: int + Position of the job in the MD, counted from 0. Unused. + + Returns + ------- + MDMaker + """ + md_maker = self.md_maker + if not isinstance(md_maker, MDMaker) or not isinstance( + md_maker.input_set_generator, MDSetGenerator + ): + raise TypeError( + "The MD maker must be a VASP MDMaker with an MDSetGenerator." + ) + generator = md_maker.input_set_generator + _check_md_incar(generator.user_incar_settings, "MD maker") + if isinstance(generator, LangevinMDSetGenerator): + if self.thermostat != "langevin": + raise ValueError( + "The MD maker has a Langevin thermostat. Set thermostat to " + "'langevin'." + ) + elif self.thermostat == "langevin": + init_fields = [f.name for f in fields(generator) if f.init] + generator = LangevinMDSetGenerator( + **{name: getattr(generator, name) for name in init_fields} + ) + generator = replace( + generator, + ensemble="nvt", + start_temp=self.temperature, + end_temp=self.temperature, + nsteps=n_steps, + time_step=self.md_time_step, + ) + return replace(md_maker, input_set_generator=generator) + + def get_npt_maker(self, n_steps: int) -> MDMaker: + """ + Get the VASP maker of the NPT MD job. + + The NPT MD uses the Langevin thermostat with the Parrinello-Rahman + barostat of VASP, MDALGO = 3 with ISIF = 3, as MDSetGenerator sets them + for ensemble="npt". PSTRESS is the pressure of the flow. + + Parameters + ---------- + n_steps: int + Number of MD steps. + + Returns + ------- + MDMaker + """ + md_maker = self.npt_maker + if ( + not isinstance(md_maker, MDMaker) + or type(md_maker.input_set_generator) is not MDSetGenerator + ): + raise TypeError( + "The NPT maker must be a VASP MDMaker with an MDSetGenerator." + ) + generator = md_maker.input_set_generator + user_incar = generator.user_incar_settings + _check_md_incar(user_incar, "NPT maker") + if "PSTRESS" in user_incar: + raise ValueError("The flow sets PSTRESS. Remove it from the NPT maker.") + generator = replace( + generator, + ensemble="npt", + start_temp=self.temperature, + end_temp=self.temperature, + nsteps=n_steps, + time_step=self.md_time_step, + user_incar_settings={**user_incar, "PSTRESS": self.pressure}, + ) + return replace(md_maker, input_set_generator=generator) diff --git a/src/atomate2/vasp/flows/pheasy.py b/src/atomate2/vasp/flows/pheasy.py index d934c538dd..1b7bd79a5b 100644 --- a/src/atomate2/vasp/flows/pheasy.py +++ b/src/atomate2/vasp/flows/pheasy.py @@ -105,9 +105,10 @@ class PhononMaker(BasePhononMaker): cutoff distance in Bohr for each FC order, starting at second order. The default value is [-1, 12, 10]. The first entry is not used, since pheasy fits the second-order FCs without a cutoff. The second and third entries - are the cutoffs for third- and fourth-order FCs. The cutoff of each fitted - order, up to anhar_max_order, must be positive. Longer cutoffs increase - the number of free FCs, and with it the number of displaced supercells. + are the cutoffs for third- and fourth-order FCs. The defaults of 12 and 10 + Bohr are 6.35 and 5.29 Å. The cutoff of each fitted order, up to + anhar_max_order, must be positive. Longer cutoffs increase the number of + free FCs, and with it the number of displaced supercells. min_length: float minimum length of lattice constants will be used to create the supercell, the default value is 8.0 A. It can be increased for larger supercells. diff --git a/src/atomate2/vasp/sets/core.py b/src/atomate2/vasp/sets/core.py index a56e943e10..5594d35a89 100644 --- a/src/atomate2/vasp/sets/core.py +++ b/src/atomate2/vasp/sets/core.py @@ -723,6 +723,42 @@ def _get_ensemble_defaults(structure: Structure, ensemble: str) -> dict[str, Any raise ValueError(f"Expect {ensemble=} to be one of {supported}") from err +@dataclass +class LangevinMDSetGenerator(MDSetGenerator): + """ + Class to generate VASP input sets for NVT MD with a Langevin thermostat. + + Parameters + ---------- + langevin_gamma + Friction coefficient in ps^-1, the same for each species. The VASP + ``LANGEVIN_GAMMA`` parameter. + **kwargs + Other keyword arguments that will be passed to :obj:`MDSetGenerator`. + """ + + langevin_gamma: float = 10.0 + + @property + def incar_updates(self) -> dict: + """Get updates to the INCAR for a Langevin NVT MD job. + + Returns + ------- + dict + A dictionary of updates to apply. + """ + if self.ensemble.lower() != "nvt": + raise ValueError(f"Expect ensemble='nvt', not {self.ensemble!r}.") + updates = super().incar_updates + # the Nose mass is only used by the Nose-Hoover thermostat + updates.pop("SMASS", None) + # one value for each element, as for ensemble="npt" + n_elements = len(self.structure.composition) + updates.update(MDALGO=3, LANGEVIN_GAMMA=[self.langevin_gamma] * n_elements) + return updates + + @dataclass class LobsterTightStaticSetGenerator(LobsterSet): """ diff --git a/tests/ase/test_md.py b/tests/ase/test_md.py index f764908555..010e7415fc 100644 --- a/tests/ase/test_md.py +++ b/tests/ase/test_md.py @@ -74,6 +74,47 @@ def test_npt_init_kwargs(si_structure, clean_dir, caplog): assert "The `NPT` module in ASE is no longer recommended" in caplog.text +def test_nvt_nose_hoover_chain(si_structure, clean_dir): + """NoseHooverChainNVT runs at a constant temperature only.""" + maker = LennardJonesMDMaker( + ensemble="nvt", + dynamics="nose-hoover-chain", + temperature=300, + n_steps=5, + ase_md_kwargs={"tdamp": 10, "tchain": 1}, + ) + result = maker.run_ase(si_structure) + assert len(result.trajectory) == 6 + + maker = LennardJonesMDMaker( + ensemble="nvt", + dynamics="nose-hoover-chain", + temperature=[300, 600], + n_steps=5, + ase_md_kwargs={"tdamp": 10}, + ) + with pytest.raises(ValueError, match="cannot follow a temperature schedule"): + maker.run_ase(si_structure) + + +def test_langevin_seed(si_structure, clean_dir): + """mb_velocity_seed also seeds the random forces of the Langevin thermostat.""" + si_structure.add_site_property("velocities", [[0.0, 0.0, 0.0]] * len(si_structure)) + + def positions(seed): + maker = LennardJonesMDMaker( + ensemble="nvt", + dynamics="langevin", + temperature=300, + n_steps=5, + mb_velocity_seed=seed, + ) + return maker.run_ase(si_structure).final_mol_or_struct.cart_coords + + assert positions(1) == pytest.approx(positions(1)) + assert positions(1) != pytest.approx(positions(2)) + + @pytest.mark.parametrize("calculator_name", list(name_to_maker)) def test_ase_nvt_maker(calculator_name, lj_fcc_ne_pars, fcc_ne_structure, clean_dir): # Langevin thermostat no longer works with single atom structures in ase>3.24.x diff --git a/tests/common/jobs/test_finite_temperature_phonons.py b/tests/common/jobs/test_finite_temperature_phonons.py new file mode 100644 index 0000000000..0991b229d8 --- /dev/null +++ b/tests/common/jobs/test_finite_temperature_phonons.py @@ -0,0 +1,486 @@ +"""Tests of the jobs of the finite-temperature phonon workflow. + +The tests of the whole force field workflow are in tests/forcefields/flows. +""" + +import gzip +from pathlib import Path + +import numpy as np +import pytest +from ase import units +from ase.build import bulk +from ase.calculators import emt +from ase.calculators.singlepoint import SinglePointCalculator +from ase.io import write as ase_write +from phonopy import Phonopy +from pymatgen.core import Lattice, Structure +from pymatgen.io.ase import AseAtomsAdaptor +from pymatgen.io.phonopy import get_phonopy_structure, get_pmg_structure + +from atomate2.common.jobs.finite_temperature_phonons import ( + ASE_TRAJECTORY_FILE, + _assess_trajectory, + _get_rms_displacement, + _remove_center_of_mass, + average_npt_structure, + fit_finite_temperature_phonons, + get_md_supercell, + get_npt_structure, + select_md_snapshots, +) + +TEMPERATURE = 300.0 +# a long time step, so that the second half of the synthetic trajectories of 400 +# frames lasts 2 ps and the mean positions are checked +TIME_STEP = 10.0 +TEST_DIR = Path(__file__).resolve().parents[2] / "test_data" + + +@pytest.fixture +def cu_supercell(): + """A 2x2x2 supercell of the cubic cell of fcc Cu, 32 atoms.""" + atoms = bulk("Cu", "fcc", a=3.61, cubic=True) * (2, 2, 2) + return AseAtomsAdaptor.get_structure(atoms) + + +def _frames(reference, n_frames, sigma, rng, shift=None): + """Fractional coordinates of frames with Gaussian displacements in Angstrom.""" + disps = rng.normal(0, sigma, size=(n_frames, len(reference), 3)) + if shift is not None: + disps += shift + cart = reference.cart_coords + disps + return np.array([reference.lattice.get_fractional_coords(c) for c in cart]) + + +def _trajectory(reference, rng, sigma=0.05, shift=None, energies=None): + """A trajectory of 400 frames that starts at the reference, and its energies.""" + frac = _frames(reference, 400, sigma, rng) + frac[:10] = _frames(reference, 10, 0.01, rng) + if shift is not None: + frac[200:] = _frames(reference, 200, sigma, rng, shift=shift) + if energies is None: + energies = rng.normal(0, 1e-3, 400) * len(reference) + return frac, energies + + +def _verdict(reference, frac, energies, time_step=TIME_STEP): + return _assess_trajectory(frac, energies, reference, time_step).verdict + + +def test_assess_trajectory(cu_supercell): + rng = np.random.default_rng(1) + n_atoms = len(cu_supercell) + + frac, energies = _trajectory(cu_supercell, rng) + health = _assess_trajectory(frac, energies, cu_supercell, TIME_STEP) + assert health.verdict == "stable" + assert health.n_frames == 400 + # u_vib is sqrt(3) sigma for Gaussian displacements + assert health.u_vib == pytest.approx(np.sqrt(3) * 0.05, rel=0.05) + assert health.u_ref**2 == pytest.approx(health.u_shift**2 + health.u_vib**2) + assert health.nearest_neighbor_distance == pytest.approx(3.61 / np.sqrt(2)) + + # a translation of the whole cell is not a displacement + moved = frac + np.array([0.1, 0.05, 0]) + health_moved = _assess_trajectory(moved, energies, cu_supercell, TIME_STEP) + assert health_moved.verdict == "stable" + assert health_moved.u_ref == pytest.approx(health.u_ref) + + # a vibration of 0.25 A per direction is above the Lindemann limit + frac, energies = _trajectory(cu_supercell, rng, sigma=0.25) + assert _verdict(cu_supercell, frac, energies) == "melted" + + # the mean positions move in the second half, with a flat energy + shift = rng.normal(0, 0.4, size=(n_atoms, 3)) + frac, energies = _trajectory(cu_supercell, rng, shift=shift) + health = _assess_trajectory(frac, energies, cu_supercell, TIME_STEP) + assert health.verdict == "shifted_or_diffusing" + assert health.shift_ratio > 1.5 + # in a second half shorter than 1 ps the mean positions are not checked + assert _verdict(cu_supercell, frac, energies, time_step=1.0) == "stable" + + # the same move with a falling energy is a transformation + falling = np.concatenate([np.zeros(200), np.linspace(0, -0.1, 200)]) * n_atoms + frac, energies = _trajectory(cu_supercell, rng, shift=shift, energies=falling) + assert _verdict(cu_supercell, frac, energies) == "transformed" + + # a falling energy without a move is a transformation too + frac, energies = _trajectory(cu_supercell, rng, energies=falling) + assert _verdict(cu_supercell, frac, energies) == "transformed" + + # a rising energy at the reference is disordering + rising = np.linspace(0, 0.1, 400) * n_atoms + frac, energies = _trajectory(cu_supercell, rng, energies=rising) + health = _assess_trajectory(frac, energies, cu_supercell, TIME_STEP) + assert health.verdict == "disordering" + assert health.energy_drift == pytest.approx(0.05, rel=0.05) + + # a large move at the very end, in a structure with long bonds, stays below + # the other limits + sparse = Structure(Lattice.cubic(8.0), ["Cu"], [[0, 0, 0]]) * (3, 3, 3) + frac, energies = _trajectory(sparse, rng) + move = rng.choice([-2.2, 2.2], size=(len(sparse), 3)) / np.sqrt(3) + frac[360:] = _frames(sparse, 40, 0.05, rng, shift=move) + health = _assess_trajectory(frac, energies, sparse, TIME_STEP) + assert health.verdict == "diffusing_or_soft" + assert health.rms_displacement_end > 1.0 + + # large displacements in the first frames mean the atom order is wrong + frac, energies = _trajectory(cu_supercell, rng) + frac[:10] = frac[:10][:, rng.permutation(n_atoms)] + assert _verdict(cu_supercell, frac, energies) == "reference_mismatch" + + # too few frames to compare the quarters + health = _assess_trajectory(frac[:15], energies[:15], cu_supercell, TIME_STEP) + assert health.energy_drift is None + + +def _write_vasp_md(directory, reference, frac, energies, gz=False): + """Write the INCAR, XDATCAR and OSZICAR of a VASP MD run with POTIM = 2.""" + directory.mkdir() + lat = reference.lattice.matrix + lines = ["Cu", "1.0", *[" ".join(f"{x:.10f}" for x in row) for row in lat]] + lines += ["Cu", str(len(reference))] + oszicar = [] + for idx, (coords, energy) in enumerate(zip(frac, energies, strict=True), 1): + lines.append(f"Direct configuration= {idx:5d}") + lines += [" ".join(f"{x:.10f}" for x in row) for row in coords] + oszicar += [ + "DAV: 1 -0.1E+02 -0.1E+02 -0.1E+02 100 0.1E+00", + ( + f"{idx:6d} T= 300. E= -.1E+02 F= {energy:.8E} " + f"E0= {energy:.8E} EK= 0.1E+00 SP= 0.0E+00 SK= 0.0E+00" + ), + ] + files = { + "INCAR": "IBRION = 0\nPOTIM = 2.0\n", + "XDATCAR": "\n".join(lines) + "\n", + "OSZICAR": "\n".join(oszicar) + "\n", + } + for name, text in files.items(): + if gz: + with gzip.open(directory / f"{name}.gz", "wt") as file: + file.write(text) + else: + (directory / name).write_text(text) + + +def test_select_md_snapshots_vasp(cu_supercell): + rng = np.random.default_rng(2) + frac = _frames(cu_supercell, 300, 0.05, rng) + energies = rng.normal(-100, 0.01, 300) + # two runs with POTIM = 2 fs, the second one gzipped + _write_vasp_md(Path("md1"), cu_supercell, frac[:150], energies[:150]) + _write_vasp_md(Path("md2"), cu_supercell, frac[150:], energies[150:], gz=True) + md_dirs = ["host:" + str(Path("md1").resolve()), str(Path("md2").resolve())] + + # 0.1 ps is 50 frames of 2 fs, so 250 frames are left for 10 snapshots + output = select_md_snapshots.original(md_dirs, "vasp", cu_supercell, 1.0, 0.1, 10) + structures = output["structures"] + assert len(structures) == 11 + assert structures[-1] == cu_supercell + # from the first to the last frame after 0.1 ps + indices = [50, 78, 105, 133, 161, 188, 216, 244, 271, 299] + for structure, idx in zip(structures[:-1], indices, strict=True): + diff = structure.frac_coords - frac[idx] + assert np.allclose(diff - np.round(diff), 0, atol=1e-8) + # XDATCAR has the frames after steps 1 to 300 + assert output["snapshot_times"] == pytest.approx( + [(idx + 1) * 0.002 for idx in indices] + ) + assert output["md_time"] == pytest.approx(0.6) + assert output["rms_displacement"] == pytest.approx(np.sqrt(3) * 0.05, rel=0.1) + assert output["trajectory_health"]["n_frames"] == 300 + + with pytest.raises(ValueError, match="fewer than the 300 snapshots"): + select_md_snapshots.original(md_dirs[:1], "vasp", cu_supercell, 1, 0.1, 300) + + # atoms of another species at the same positions + other = cu_supercell.copy() + other.replace_species({"Cu": "Ag"}) + with pytest.raises(ValueError, match="not in the order of the reference"): + select_md_snapshots.original(md_dirs[:1], "vasp", other, 1, 0.1, 10) + + Path("md4").mkdir() + for name in ("INCAR", "XDATCAR"): + (Path("md4") / name).write_text((Path("md1") / name).read_text()) + (Path("md4") / "OSZICAR").write_text("") + with pytest.raises(ValueError, match=r"0 MD steps in OSZICAR\..*NBLOCK = 1"): + select_md_snapshots.original( + [str(Path("md4").resolve())], "vasp", cu_supercell, 1, 0.1, 10 + ) + + +def _write_ase_md(directory, reference, frac, energies, cells=None): + frames = [] + for idx, (coords, energy) in enumerate(zip(frac, energies, strict=True)): + atoms = AseAtomsAdaptor.get_atoms(reference) + if cells is not None: + atoms.set_cell(cells[idx]) + atoms.set_scaled_positions(coords) + atoms.calc = SinglePointCalculator(atoms, energy=energy) + frames.append(atoms) + directory.mkdir() + ase_write(directory / ASE_TRAJECTORY_FILE, frames) + + +def test_select_md_snapshots_forcefields(cu_supercell): + rng = np.random.default_rng(3) + frac = _frames(cu_supercell, 151, 0.05, rng) + energies = rng.normal(0, 0.01, 151) + # the second run starts from the last frame of the first one + _write_ase_md(Path("md1"), cu_supercell, frac[:101], energies[:101]) + _write_ase_md(Path("md2"), cu_supercell, frac[100:], energies[100:]) + + output = select_md_snapshots.original( + [str(Path(name).resolve()) for name in ("md1", "md2")], + "forcefields", + cu_supercell, + 2.0, + 0.02, + 14, + ) + # the starting structure of each run is left out, so frame 0 is step 1 + assert output["trajectory_health"]["n_frames"] == 150 + # 10 frames of 2 fs are left out, and 140 frames are left for 14 snapshots + indices = np.linspace(10, 149, 14).round().astype(int) + assert output["snapshot_times"] == pytest.approx((indices + 1) * 0.002) + for structure, idx in zip(output["structures"][:-1], indices, strict=True): + diff = structure.frac_coords - frac[idx + 1] + assert np.allclose(diff - np.round(diff), 0, atol=1e-8) + + # the snapshots keep the magnetic moments of the reference + magmoms = [1.0] * len(cu_supercell) + magnetic = cu_supercell.copy(site_properties={"magmom": magmoms}) + output = select_md_snapshots.original( + [str(Path("md1").resolve())], "forcefields", magnetic, 2.0, 0.02, 14 + ) + for structure in output["structures"]: + assert structure.site_properties == {"magmom": magmoms} + + # a trajectory that leaves the reference gives a warning + melted = _frames(cu_supercell, 101, 0.3, rng) + melted[:10] = cu_supercell.frac_coords + _write_ase_md(Path("md3"), cu_supercell, melted, energies[:101]) + with pytest.warns(UserWarning, match="verdict: melted"): + output = select_md_snapshots.original( + [str(Path("md3").resolve())], + "forcefields", + cu_supercell, + 2.0, + 0.02, + 14, + ) + assert output["trajectory_health"]["verdict"] == "melted" + + +def test_average_npt_structure(): + """Rotated cells of an expanded tetragonal supercell average to its unit cell.""" + structure = Structure( + Lattice.tetragonal(3.0, 5.0), ["Si", "Si"], [[0, 0, 0], [0.5, 0.5, 0.5]] + ) + matrix = np.diag([2, 2, 1]) + supercell = np.diag([1.01, 1.01, 1.03]) @ matrix @ structure.lattice.matrix + rng = np.random.default_rng(7) + cells = [] + for _ in range(5): + rotation, _ = np.linalg.qr(rng.normal(size=(3, 3))) + cells.append(supercell @ rotation.T) + averaged = average_npt_structure(np.array(cells), structure, matrix, 1e-3) + assert np.allclose(averaged.lattice.matrix, np.diag([3.03, 3.03, 5.15])) + assert np.allclose(averaged.frac_coords, structure.frac_coords) + + +def test_get_npt_structure(cu_supercell): + unit_cell = AseAtomsAdaptor.get_structure(bulk("Cu", "fcc", a=3.61, cubic=True)) + rng = np.random.default_rng(4) + n_frames = 201 + frac = _frames(cu_supercell, n_frames, 0.05, rng) + # the cell is 2% longer in each direction, with random strains, and it + # rotates about the z axis + cells = [] + for angle in rng.uniform(0, 2 * np.pi, n_frames): + strain = np.eye(3) + rng.normal(0, 0.005, (3, 3)) + cos, sin = np.cos(angle), np.sin(angle) + rotation = np.array([[cos, -sin, 0], [sin, cos, 0], [0, 0, 1]]) + cells.append(1.02 * cu_supercell.lattice.matrix @ strain @ rotation.T) + energies = rng.normal(0, 0.01, n_frames) + _write_ase_md(Path("npt"), cu_supercell, frac, energies, cells=cells) + + args = ("forcefields", unit_cell, (2 * np.eye(3)).tolist(), cu_supercell) + output = get_npt_structure.original( + str(Path("npt").resolve()), *args, 2.0, 0.02, 1e-4 + ) + structure = output["structure"] + # the cell stays cubic and keeps the orientation of the unit cell + assert np.allclose( + structure.lattice.matrix, 1.02 * unit_cell.lattice.matrix, atol=5e-3 + ) + assert structure.lattice.abc == pytest.approx([structure.lattice.a] * 3) + assert np.allclose(structure.frac_coords, unit_cell.frac_coords) + assert output["trajectory_health"]["verdict"] == "stable" + + with pytest.raises(ValueError, match="none of them after"): + get_npt_structure.original(str(Path("npt").resolve()), *args, 2.0, 1.0, 1e-4) + + +def test_get_md_supercell(): + """The supercell has the atom order of phonopy and the magnetic moments.""" + structure = Structure( + Lattice.cubic(4.17), + ["O", "Ni", "O", "Ni"], + [[0.5, 0, 0], [0, 0, 0], [0, 0.5, 0], [0.5, 0.5, 0]], + site_properties={"magmom": [0.0, 2.0, 0.0, -2.0]}, + ) + matrix = np.diag([2, 2, 1]) + supercell = get_md_supercell.original(structure, matrix, 1e-3, "vasp") + assert len(supercell) == 16 + assert [str(site.specie) for site in supercell] == ["Ni"] * 8 + ["O"] * 8 + # each atom of the sorted unit cell has four images in a row + assert supercell.site_properties["magmom"] == [2.0] * 4 + [-2.0] * 4 + [0.0] * 8 + + supercell = get_md_supercell.original( + structure.copy().remove_site_property("magmom"), + np.diag([2, 2, 1]), + 1e-3, + "vasp", + ) + assert supercell.site_properties == {} + + +def test_get_rms_displacement(): + """An Einstein crystal has = 3 k_B T / k per atom.""" + structure = AseAtomsAdaptor.get_structure(bulk("Cu", "fcc", a=3.61, cubic=True)) + phonon = Phonopy(get_phonopy_structure(structure), np.diag([2, 2, 2])) + n_atoms = len(phonon.supercell) + spring = 2.0 # eV/A^2 + force_constants = np.zeros((n_atoms, n_atoms, 3, 3)) + force_constants[np.arange(n_atoms), np.arange(n_atoms)] = spring * np.eye(3) + phonon.force_constants = force_constants + rms = _get_rms_displacement(phonon, TEMPERATURE, 0.1) + assert rms == pytest.approx(np.sqrt(3 * units.kB * TEMPERATURE / spring)) + + +def test_rms_displacement_matches_mode_sampling(): + """Canonical mode sampling of Cu3Au gives the rms of the force constants. + + Each sample gets a random translation. Removing the displacement of the + center of mass takes it out again. The unweighted mean of the displacements + is not zero for two species, so it would not. + """ + structure = Structure( + Lattice.cubic(3.75), + ["Cu", "Cu", "Cu", "Au"], + [[0, 0.5, 0.5], [0.5, 0, 0.5], [0.5, 0.5, 0], [0, 0, 0]], + ) + phonon = _emt_phonon(structure, np.diag([2, 2, 2])) + masses = np.repeat(phonon.supercell.masses, 3) + n_atoms = len(phonon.supercell) + dynmat = phonon.force_constants.transpose(0, 2, 1, 3).reshape( + 3 * n_atoms, 3 * n_atoms + ) / np.sqrt(np.outer(masses, masses)) + eigvals, eigvecs = np.linalg.eigh((dynmat + dynmat.T) / 2) + keep = eigvals > 1e-6 + + rng = np.random.default_rng(6) + n_samples = 4000 + amplitudes = rng.normal(size=(n_samples, keep.sum())) * np.sqrt( + units.kB * TEMPERATURE / eigvals[keep] + ) + disps = (amplitudes @ eigvecs[:, keep].T / np.sqrt(masses)).reshape( + n_samples, n_atoms, 3 + ) + disps += rng.normal(0, 0.1, size=(n_samples, 1, 3)) + + reference = get_pmg_structure(phonon.supercell) + removed = _remove_center_of_mass(disps, reference) + rms = np.sqrt(np.mean(np.sum(removed**2, axis=2))) + assert rms == pytest.approx( + _get_rms_displacement(phonon, TEMPERATURE, 0.1), rel=0.02 + ) + unweighted = disps - disps.mean(axis=1, keepdims=True) + assert not np.allclose(unweighted, removed, atol=1e-3) + + +def _emt_phonon(structure, supercell_matrix): + """0 K phonopy object of EMT, with its force constants.""" + phonon = Phonopy( + get_phonopy_structure(structure), supercell_matrix, primitive_matrix="auto" + ) + phonon.generate_displacements(distance=0.01) + forces = [] + for cell in phonon.supercells_with_displacements: + atoms = AseAtomsAdaptor.get_atoms(get_pmg_structure(cell)) + atoms.calc = emt.EMT() + forces.append(atoms.get_forces()) + phonon.forces = forces + phonon.produce_force_constants() + return phonon + + +def test_fit_recovers_harmonic_force_constants(): + """Forces from known force constants give them back, with the NAC applied.""" + structure = Structure( + Lattice.cubic(3.75), + ["Cu", "Cu", "Cu", "Au"], + [[0, 0.5, 0.5], [0.5, 0, 0.5], [0.5, 0.5, 0], [0, 0, 0]], + ) + supercell_matrix = np.diag([2, 2, 2]).tolist() + phonon = _emt_phonon(structure, supercell_matrix) + supercell = get_pmg_structure(phonon.supercell) + disps = np.random.default_rng(0).normal(0, 0.03, (20, len(supercell), 3)) + forces = -np.einsum("ijab,mjb->mia", phonon.force_constants, disps) + snapshots = [ + Structure( + supercell.lattice, + supercell.species, + supercell.cart_coords + disp, + coords_are_cartesian=True, + ) + for disp in disps + ] + displacement_data = { + "forces": [*forces, np.zeros((len(supercell), 3))], + "displaced_structures": [*snapshots, supercell], + "uuids": ["uuid"] * 21, + "dirs": ["dir"] * 21, + } + snapshot_data = { + "snapshot_times": [0.0] * 20, + "md_time": 1.0, + "rms_displacement": 0.05, + "trajectory_health": {}, + } + born = [np.eye(3).tolist()] * 3 + [(-3 * np.eye(3)).tolist()] + fit_kwargs = { + "structure": structure, + "supercell_matrix": supercell_matrix, + "snapshot_data": snapshot_data, + "displacement_data": displacement_data, + "temperature": TEMPERATURE, + "thermostat": "langevin", + "md_time_step": 1.0, + "equilibration_time": 0.5, + "code": "forcefields", + "md_code": "forcefields", + "symprec": 1e-3, + "rotational_sum_rule": None, + "epsilon_static": (10 * np.eye(3)).tolist(), + "npoints_band": 11, + } + doc = fit_finite_temperature_phonons.original(born=born, **fit_kwargs) + assert doc.force_rmse < 0.01 * np.sqrt(np.mean(forces**2)) + force_constants = np.array(doc.force_constants.force_constants) + assert np.abs(force_constants - phonon.force_constants).max() < 0.05 + assert np.array(doc.born).shape == (4, 3, 3) + assert doc.phonon_bandstructure.has_nac + assert not doc.has_imaginary_modes + assert doc.uuids.displacements_uuids == ["uuid"] * 21 + assert doc.jobdirs.displacements_job_dirs == ["dir"] * 21 + + # the Born charges are checked before pheasy runs + fc_time = Path("FORCE_CONSTANTS").stat().st_mtime_ns + with pytest.raises(ValueError, match="number of Born charges"): + fit_finite_temperature_phonons.original(born=born[:3], **fit_kwargs) + assert Path("FORCE_CONSTANTS").stat().st_mtime_ns == fc_time diff --git a/tests/common/jobs/test_md.py b/tests/common/jobs/test_md.py new file mode 100644 index 0000000000..34629d3987 --- /dev/null +++ b/tests/common/jobs/test_md.py @@ -0,0 +1,62 @@ +"""Tests of the shared MD jobs and flows.""" + +from pathlib import Path + +import numpy as np +import pytest +from ase.build import bulk +from ase.io import write as ase_write +from pymatgen.core import Structure +from pymatgen.io.ase import AseAtomsAdaptor + +from atomate2.common.flows.md import ChainedMDMaker +from atomate2.common.jobs.md import get_md_restart_structure +from atomate2.forcefields.md import ForceFieldMDMaker + +TEST_DIR = Path(__file__).resolve().parents[2] / "test_data" + + +def test_get_md_restart_structure(): + cu_supercell = AseAtomsAdaptor.get_structure( + bulk("Cu", "fcc", a=3.61, cubic=True) * (2, 2, 2) + ) + rng = np.random.default_rng(4) + velocities = rng.normal(0, 0.01, size=(len(cu_supercell), 3)) + reference = cu_supercell.copy(site_properties={"magmom": [1.0] * 32}) + + # a CONTCAR of a VASP MD, with the predictor-corrector block after the + # velocities + md_dir = TEST_DIR / "vasp/Si_multi_md/molecular_dynamics_1/outputs" + si_reference = Structure( + np.eye(3) * 3.87, ["Si", "Si"], [[0, 0, 0], [0.25, 0.25, 0.25]] + ) + restart = get_md_restart_structure.original(str(md_dir), si_reference) + assert restart.site_properties["velocities"][0] == pytest.approx( + [0.40252568e-03, -0.31480439e-02, -0.16045943e-02] + ) + assert "magmom" not in restart.site_properties + + Path("vasp").mkdir() + structure = cu_supercell.copy(site_properties={"velocities": velocities}) + structure.to(filename="vasp/CONTCAR", fmt="poscar") + restart = get_md_restart_structure.original(str(Path("vasp").resolve()), reference) + assert np.allclose(restart.site_properties["velocities"], velocities, atol=1e-8) + assert restart.site_properties["magmom"] == [1.0] * 32 + + Path("ase").mkdir() + atoms = AseAtomsAdaptor.get_atoms(cu_supercell) + atoms.set_velocities(velocities) + ase_write(Path("ase") / "md.traj", [atoms.copy(), atoms]) + restart = get_md_restart_structure.original( + str(Path("ase").resolve()), reference, "md.traj" + ) + assert np.allclose(restart.site_properties["velocities"], velocities) + assert np.allclose(restart.cart_coords, cu_supercell.cart_coords) + assert restart.site_properties["magmom"] == [1.0] * 32 + + +def test_chained_md_maker_needs_ase_trajectory(): + structure = Structure(np.eye(3) * 3.61, ["Cu"], [[0, 0, 0]]) + maker = ChainedMDMaker(md_makers=[ForceFieldMDMaker(), ForceFieldMDMaker()]) + with pytest.raises(ValueError, match="trajectory in the ASE format"): + maker.make(structure) diff --git a/tests/common/jobs/test_pheasy.py b/tests/common/jobs/test_pheasy.py index fc1d816d29..056b1e47a7 100644 --- a/tests/common/jobs/test_pheasy.py +++ b/tests/common/jobs/test_pheasy.py @@ -21,6 +21,7 @@ from atomate2.common.jobs.pheasy import ( _check_lasso_alpha, _get_num_anharmonic_supercells, + _run_harmonic_fit, generate_frequencies_eigenvectors, generate_phonon_displacements, ) @@ -260,17 +261,59 @@ def test_get_num_anharmonic_supercells(monkeypatch): _get_num_anharmonic_supercells(num_disp_anhar=0, **kwargs) +def test_run_harmonic_fit(monkeypatch): + """The fit flags of the pheasy workflow and the finite-temperature fit.""" + calls = [] + + def fake_run(args, check): + calls.append((args, check)) + + monkeypatch.setattr(pheasy_jobs.subprocess, "run", fake_run) + matrix = np.diag([2, 3, 4]) + + _run_harmonic_fit(matrix, 1e-3, 5) + assert [args[:6] for args, _ in calls] == [ + ["pheasy", "--scell", "SPOSCAR", "--dim", "2", "3"] + ] * 4 + assert all(check for _, check in calls) + assert [args[11] for args, _ in calls] == ["-s", "-c", "-d", "-f"] + fit = " ".join(calls[-1][0]) + assert "-l LASSO --std --tol 1e-8 --seed 103 --rasr BHH --ndata 5" in fit + assert "--alpha_min" not in fit + assert " -o " not in fit + + calls.clear() + _run_harmonic_fit(matrix, 1e-3, 3, use_lasso=False) + fit = " ".join(calls[-1][0]) + assert "-f --full_ifc --rasr BHH --ndata 3" in fit + for flag in ("-l", "--std", "--tol", "--seed"): + assert flag not in calls[-1][0] + + calls.clear() + _run_harmonic_fit( + matrix, 1e-3, 5, rotational_sum_rule=None, alpha_min=-8, log_file="x.log" + ) + fit = " ".join(calls[-1][0]) + assert "--alpha_min -8 --seed 103 --ndata 5" in fit + assert fit.endswith("-o x.log") + assert "--rasr" not in fit + + def test_check_lasso_alpha(tmp_dir): log_file = Path("pheasy_anharmonic_fit.log") log_file.write_text("- alpha_min: 1e-12\n- alpha_opt: 2.947052e-09\n") with warnings.catch_warnings(): warnings.simplefilter("error") - _check_lasso_alpha(log_file, alpha_min=-12) + assert _check_lasso_alpha(log_file, alpha_min=-12) == pytest.approx( + 2.947052e-09 + ) log_file.write_text("- alpha_min: 1e-12\n- alpha_opt: 1.000000e-12\n") - with pytest.warns(UserWarning, match="on the lower bound"): + with pytest.warns(UserWarning, match="on the lower bound.*anhar_alpha_min"): _check_lasso_alpha(log_file, alpha_min=-12) + with pytest.warns(UserWarning, match="Lower alpha_min and refit"): + _check_lasso_alpha(log_file, alpha_min=-12, alpha_min_name="alpha_min") log_file.write_text("- alpha_max: 1e-2\n- alpha_opt: 1.000000e-02\n") with pytest.warns(UserWarning, match="on the upper bound"): diff --git a/tests/forcefields/flows/test_finite_temperature_phonons.py b/tests/forcefields/flows/test_finite_temperature_phonons.py new file mode 100644 index 0000000000..fa5f3841ec --- /dev/null +++ b/tests/forcefields/flows/test_finite_temperature_phonons.py @@ -0,0 +1,268 @@ +"""Tests for the force field finite-temperature phonon workflow. + +The tests run the workflow with ASE's EMT potential, so they need no reference +data and no machine-learned force field. pheasy is only installed in the +numpy-limited forcefield CI job, so the tests are skipped in the other +forcefield jobs. +""" + +import pytest + +pytest.importorskip("pheasy") + +import numpy as np +from ase import units +from jobflow import run_locally +from pymatgen.core import Lattice, Structure + +from atomate2.ase.md import MDEnsemble +from atomate2.forcefields.flows.finite_temperature_phonons import ( + ForceFieldFiniteTemperaturePhononMaker, + _get_force_field_md_maker, +) +from atomate2.forcefields.jobs import ForceFieldDielectricMaker, ForceFieldStaticMaker +from atomate2.forcefields.md import ForceFieldMDMaker +from atomate2.vasp.jobs.core import DielectricMaker + +EMT = {"@module": "ase.calculators.emt", "@callable": "EMT"} + + +@pytest.fixture +def cu3au(): + """L1_2 Cu3Au with Au first, so that sorting by electronegativity reorders it.""" + return Structure( + Lattice.cubic(3.75), + ["Au", "Cu", "Cu", "Cu"], + [[0, 0, 0], [0, 0.5, 0.5], [0.5, 0, 0.5], [0.5, 0.5, 0]], + ) + + +def run_emt_flow(structure, **kwargs): + """Run the force field workflow with EMT and return its output document.""" + maker = ForceFieldFiniteTemperaturePhononMaker.from_force_field_name( + EMT, min_length=7.0, md_time_step=2.0, **kwargs + ) + flow = maker.make(structure) + responses = run_locally(flow, create_folders=True, ensure_success=True) + return responses[flow.output.uuid][1].output + + +@pytest.mark.parametrize( + ("kwargs", "match"), + [ + ({"alpha_min": -2}, "alpha_min must be below -2"), + ({"md_time_step": 0}, "must be positive"), + ({"equilibration_time": 8.0}, "fewer than the 50 snapshots"), + ({"equilibration_time": -1.0}, "fewer than the 50 snapshots"), + ({"md_maker": None}, "md_maker and phonon_displacement_maker must be set"), + ({"code": "aims"}, "code must be one of"), + ({"md_code": None}, "md_code must be one of"), + ( + {"npt_maker": ForceFieldMDMaker(), "npt_equilibration_time": 8.0}, + "npt_equilibration_time must be", + ), + ], +) +def test_maker_checks(kwargs, match): + with pytest.raises(ValueError, match=match): + ForceFieldFiniteTemperaturePhononMaker(**kwargs) + + +def test_force_field_md_maker(cu3au): + maker = ForceFieldFiniteTemperaturePhononMaker.from_force_field_name( + EMT, + md_time=1.0, + md_time_step=2.0, + equilibration_time=0.5, + md_runs=3, + random_seed=5, + ) + assert maker.get_md_steps() == [167, 167, 166] + md_maker = maker.get_md_maker(167) + assert md_maker.n_steps == 167 + assert md_maker.time_step == 2.0 + assert md_maker.temperature == 300 + assert md_maker.mb_velocity_seed == 5 + assert maker.get_md_maker(166, 2).mb_velocity_seed == 7 + assert md_maker.store_trajectory == "no" + # Langevin by default, with the default friction of AseMDMaker + assert md_maker.dynamics == "langevin" + assert md_maker.ase_md_kwargs == {} + # the maker of the flow is not changed + assert maker.md_maker.n_steps == 1000 + assert maker.md_maker.dynamics is None + # the force field relaxation keeps the symmetry + assert maker.bulk_relax_maker.fix_symmetry + assert ForceFieldFiniteTemperaturePhononMaker().bulk_relax_maker.fix_symmetry + + # the MD jobs have a number only when there are several + flow = maker.make(cu3au, supercell_matrix=np.eye(3).tolist()) + md_flow = next(job for job in flow.jobs if job.name == "chained MD") + md_names = [job.name for job in md_flow.jobs if "MD" in job.name] + assert md_names == ["ASE MD 1/3", "ASE MD 2/3", "ASE MD 3/3"] + maker = ForceFieldFiniteTemperaturePhononMaker(md_time=2.0) + flow = maker.make(cu3au, supercell_matrix=np.eye(3).tolist()) + md_flow = next(job for job in flow.jobs if job.name == "chained MD") + assert [job.name for job in md_flow.jobs] == ["ASE MD"] + + nose_hoover = ForceFieldFiniteTemperaturePhononMaker( + thermostat="nose-hoover", md_time_step=2.0 + ).get_md_maker(10) + assert nose_hoover.dynamics == "nose-hoover-chain" + assert nose_hoover.ase_md_kwargs["tchain"] == 1 + # a thermostat period of 40 steps, 80 fs + period = np.pi * np.sqrt(2) * nose_hoover.ase_md_kwargs["tdamp"] / units.fs + assert period == pytest.approx(80.0) + + maker = ForceFieldFiniteTemperaturePhononMaker( + md_maker=ForceFieldMDMaker(ase_md_kwargs={"tdamp": 10}) + ) + with pytest.raises(ValueError, match="must not set dynamics or ase_md_kwargs"): + maker.get_md_maker(10) + with pytest.raises(TypeError, match="must be a ForceFieldMDMaker"): + _get_force_field_md_maker( + ForceFieldFiniteTemperaturePhononMaker(md_maker=ForceFieldStaticMaker()), 10 + ) + + with pytest.raises(ValueError, match="diagonal supercell matrix"): + maker.make(cu3au, supercell_matrix=[[1, 1, 0], [0, 1, 0], [0, 0, 1]]) + + +def test_force_field_npt_maker(cu3au): + maker = ForceFieldFiniteTemperaturePhononMaker.from_force_field_name( + EMT, run_npt=True, md_time_step=2.0, pressure=10.0 + ) + npt_maker = maker.get_npt_maker(100) + assert npt_maker.ensemble == MDEnsemble.npt + assert npt_maker.dynamics == "nose-hoover-chain" + assert npt_maker.n_steps == 100 + assert npt_maker.pressure == 10.0 + # a barostat time constant of 1000 steps + assert npt_maker.ase_md_kwargs["pdamp"] / units.fs == pytest.approx(2000) + assert maker.get_md_maker(100).ensemble == MDEnsemble.nvt + # the NPT MD has its own seed, after those of the md_runs NVT MD jobs + assert npt_maker.mb_velocity_seed == 104 + # the atoms are relaxed in the cell from the NPT MD + assert not maker.fixed_cell_relax_maker.relax_cell + assert maker.fixed_cell_relax_maker.fix_symmetry + + flow = maker.make(cu3au, supercell_matrix=np.eye(3).tolist()) + names = [job.name for job in flow.jobs] + assert "ASE MD NPT" in names + assert "get_npt_structure" in names + + # no NPT MD by default + maker = ForceFieldFiniteTemperaturePhononMaker.from_force_field_name(EMT) + assert maker.npt_maker is None + assert maker.fixed_cell_relax_maker is None + flow = maker.make(cu3au, supercell_matrix=np.eye(3).tolist()) + assert not any("NPT" in job.name for job in flow.jobs) + + +def test_force_field_with_vasp_born_charges(cu3au): + """A VASP born_maker after a force field relaxation gets no prev_dir.""" + maker = ForceFieldFiniteTemperaturePhononMaker.from_force_field_name(EMT) + maker.born_maker = DielectricMaker() + flow = maker.make(cu3au, supercell_matrix=np.eye(3).tolist()) + born_job = next(job for job in flow.jobs if job.name == "dielectric") + assert born_job.function_kwargs.get("prev_dir") is None + + +def test_finite_temperature_phonon_maker_emt(clean_dir, cu3au): + """Run the whole workflow with EMT on Cu3Au at 300 K.""" + # a 2x2x2 supercell of the cubic cell, 32 atoms, and 2 ps of MD at 2 fs in + # two MD jobs + doc = run_emt_flow( + cu3au, md_time=2.0, md_runs=2, equilibration_time=0.5, n_snapshots=12 + ) + + assert [str(site.specie) for site in doc.structure] == ["Cu", "Cu", "Cu", "Au"] + assert doc.temperature == 300 + assert doc.thermostat == "langevin" + assert doc.md_code == doc.code == "forcefields" + assert doc.md_force_field_name == doc.force_field_name == "ase.calculators.emt.EMT" + assert doc.force_field_kwargs == {} + assert doc.md_time == pytest.approx(2.0) + assert doc.n_snapshots == 12 + # 250 steps of 2 fs are left out + assert doc.snapshot_times[0] == pytest.approx(0.502) + assert doc.snapshot_times[-1] == pytest.approx(2.0) + assert len(doc.md_uuids) == len(doc.md_job_dirs) == 2 + assert len(doc.uuids.displacements_uuids) == 13 + assert doc.uuids.optimization_run_uuid is not None + assert doc.uuids.born_run_uuid is None + assert doc.supercell_matrix == ((2, 0, 0), (0, 2, 0), (0, 0, 2)) + assert doc.trajectory_health.verdict == "stable" + assert doc.max_residual_force < 1e-8 + assert np.array(doc.force_constants.force_constants).shape == (32, 32, 3, 3) + # Cu3Au is stable at 300 K + assert not doc.has_imaginary_modes + assert doc.n_imaginary_modes == 0 + # the snapshot amplitude agrees with the one of the fitted force constants + assert doc.rms_displacement == pytest.approx( + doc.rms_displacement_from_force_constants, rel=0.25 + ) + # values of this run, to catch changes of the results + assert np.max(doc.phonon_bandstructure.bands) == pytest.approx(6.815, rel=0.05) + # the acoustic modes at Gamma + assert doc.lowest_frequency == pytest.approx(0, abs=1e-3) + assert doc.rms_displacement == pytest.approx(0.1235, rel=0.05) + assert doc.force_rmse == pytest.approx(0.1045, rel=0.15) + + +def test_finite_temperature_phonon_maker_emt_npt(clean_dir, cu3au): + """Run the workflow with an NPT MD first, with EMT on Cu3Au at 300 K.""" + doc = run_emt_flow( + cu3au, + run_npt=True, + npt_time=2.0, + npt_equilibration_time=0.5, + md_time=1.0, + equilibration_time=0.2, + n_snapshots=10, + ) + + assert doc.pressure == 0.0 + assert doc.npt_time == 2.0 + assert doc.npt_equilibration_time == 0.5 + assert doc.npt_trajectory_health.verdict == "stable" + assert doc.npt_uuid is not None + assert doc.npt_job_dir is not None + assert doc.fixed_cell_relax_uuid is not None + # the atoms are relaxed in the cell from the NPT MD + assert doc.max_residual_force < 1e-3 + # the NVT MD and the fit use the cubic cell from the NPT MD, which expanded + assert doc.structure.lattice.abc == pytest.approx([doc.structure.lattice.a] * 3) + assert doc.structure.lattice.angles == pytest.approx((90, 90, 90)) + assert doc.structure.volume > doc.npt_input_structure.volume + assert doc.structure.lattice.a == pytest.approx(3.7328, rel=2e-3) + assert doc.trajectory_health.verdict == "stable" + assert not doc.has_imaginary_modes + assert np.max(doc.phonon_bandstructure.bands) == pytest.approx(6.554, rel=0.05) + + +def test_finite_temperature_phonon_maker_emt_born_charges( + clean_dir, cu3au, fake_dielectric_calculator +): + """Born charges from a force field dielectric job enter the NAC.""" + maker = ForceFieldFiniteTemperaturePhononMaker.from_force_field_name( + EMT, + min_length=7.0, + md_time_step=2.0, + md_time=1.0, + equilibration_time=0.2, + n_snapshots=10, + ) + maker.born_maker = ForceFieldDielectricMaker() + flow = maker.make(cu3au) + responses = run_locally(flow, create_folders=True, ensure_success=True) + doc = responses[flow.output.uuid][1].output + + # the fake calculator gives +2 to the first species of the sorted cell, Cu + born = np.array(doc.born) + assert born.shape == (4, 3, 3) + assert np.all(np.diagonal(born[:3], axis1=1, axis2=2) > 0) + assert np.all(np.diagonal(born[3:], axis1=1, axis2=2) < 0) + assert np.array(doc.epsilon_static) == pytest.approx(4 * np.eye(3)) + assert doc.uuids.born_run_uuid is not None + assert doc.phonon_bandstructure.has_nac diff --git a/tests/vasp/flows/test_finite_temperature_phonons.py b/tests/vasp/flows/test_finite_temperature_phonons.py new file mode 100644 index 0000000000..87dcd5f959 --- /dev/null +++ b/tests/vasp/flows/test_finite_temperature_phonons.py @@ -0,0 +1,296 @@ +import pytest +from jobflow import Flow, OutputReference +from pymatgen.core import Lattice, Structure + +from atomate2.common.jobs.finite_temperature_phonons import ASE_TRAJECTORY_FILE +from atomate2.forcefields.flows.finite_temperature_phonons import ( + MLFFMDVaspStaticFiniteTemperaturePhononMaker, + VaspMDMLFFStaticFiniteTemperaturePhononMaker, +) +from atomate2.forcefields.jobs import ForceFieldStaticMaker +from atomate2.forcefields.md import ForceFieldMDMaker +from atomate2.vasp.flows.finite_temperature_phonons import FiniteTemperaturePhononMaker +from atomate2.vasp.jobs.core import StaticMaker +from atomate2.vasp.jobs.md import MDMaker +from atomate2.vasp.sets.core import LangevinMDSetGenerator, MDSetGenerator + + +@pytest.fixture +def nio_supercell(): + """A 3x3x3 supercell of rock-salt NiO, which gets +U in VASP.""" + structure = Structure( + Lattice.cubic(4.17), ["Ni", "O"], [[0, 0, 0], [0.5, 0.5, 0.5]] + ) + return structure * (3, 3, 3) + + +@pytest.mark.parametrize( + ("maker_cls", "md_code", "code"), + [ + (FiniteTemperaturePhononMaker, "vasp", "vasp"), + (VaspMDMLFFStaticFiniteTemperaturePhononMaker, "vasp", "forcefields"), + (MLFFMDVaspStaticFiniteTemperaturePhononMaker, "forcefields", "vasp"), + ], +) +def test_vasp_flows(si_structure, maker_cls, md_code, code): + maker = maker_cls(md_runs=2) + assert (maker.md_code, maker.code) == (md_code, code) + flow = maker.make(si_structure) + names = [job.name for job in flow.jobs] + # the Born charges are computed when the statics use VASP, as in the + # pheasy phonon workflows + born_names = ["dielectric"] if code == "vasp" else [] + md_names = ["molecular dynamics" if md_code == "vasp" else "ASE MD"] * 2 + assert names == [ + "double relax", + "get_supercell_size", + *born_names, + "get_md_supercell", + "chained MD", + "select_md_snapshots", + "run_phonon_displacements", + "fit_finite_temperature_phonons", + ] + relax = flow.jobs[0] + reference, md_flow, snapshots, statics, fit = flow.jobs[-5:] + md_1, restart, md_2 = md_flow.jobs + assert [job.name for job in md_flow.jobs] == [ + f"{md_names[0]} 1/2", + "get_md_restart_structure", + f"{md_names[1]} 2/2", + ] + assert isinstance(relax, Flow) + + # a diagonal supercell + supercell = flow.jobs[1] + assert supercell.function_args[1:] == (12.0, None) + assert supercell.function_kwargs["force_diagonal"] + + # the first MD starts from the undisplaced supercell, the second from the + # end of the first, and both get the relaxation directory + assert md_1.function_args[0].uuid == reference.uuid + assert md_2.function_args[0].uuid == restart.uuid + traj_file = ASE_TRAJECTORY_FILE if md_code == "forcefields" else None + assert restart.function_args == (md_1.output.dir_name, reference.output, traj_file) + for md_job in (md_1, md_2): + assert md_job.function_kwargs["prev_dir"].uuid == relax.output.uuid + assert maker.get_md_steps() == [4000, 4000] + + args = snapshots.function_args + assert [ref.uuid for ref in args[0]] == [md_1.uuid, md_2.uuid] + assert args[1:] == (md_code, reference.output, 1.0, 1.0, 50) + + # the statics get the snapshots followed by the undisplaced supercell, and + # only VASP statics get the relaxation directory + static_kwargs = statics.function_kwargs + assert static_kwargs["displacements"].uuid == snapshots.uuid + assert static_kwargs["prev_dir"].uuid == relax.output.uuid + assert static_kwargs["prev_dir_argname"] == ("prev_dir" if code == "vasp" else None) + assert static_kwargs["phonon_maker"] == maker.phonon_displacement_maker + + fit_kwargs = fit.function_kwargs + assert fit_kwargs["md_code"] == md_code + assert fit_kwargs["code"] == code + assert fit_kwargs["md_uuids"] == [md_1.uuid, md_2.uuid] + assert fit_kwargs["displacement_data"].uuid == statics.uuid + assert fit_kwargs["optimization_run_uuid"] == relax.output.uuid + assert fit_kwargs["rotational_sum_rule"] == "BHH" + assert fit_kwargs["alpha_min"] == -6 + if code == "forcefields": + assert fit_kwargs["force_field_name"] == "MACE-MP-0" + assert fit_kwargs["force_field_kwargs"] == {"model": "medium"} + assert fit_kwargs["md_force_field_name"] is None + else: + assert fit_kwargs["force_field_name"] is None + if md_code == "forcefields": + assert fit_kwargs["md_force_field_name"] == "MACE-MP-0" + if code == "vasp": + born = flow.jobs[2] + assert fit_kwargs["born"].uuid == born.uuid + assert fit_kwargs["epsilon_static"].uuid == born.uuid + assert fit_kwargs["born_run_uuid"] == born.uuid + assert born.function_kwargs["prev_dir"].uuid == relax.output.uuid + else: + assert fit_kwargs["born"] is None + assert fit_kwargs["born_run_uuid"] is None + assert fit.output.uuid == flow.output.uuid + + +def test_flow_from_output_reference(): + """The structure can be the output of an earlier job.""" + structure = OutputReference("1234", attributes=(("a", "structure"),)) + flow = FiniteTemperaturePhononMaker().make(structure) + assert flow.jobs[0].name == "double relax" + assert [job.name for job in flow.jobs][4] == "chained MD" + + +def test_non_diagonal_supercell_matrix(si_structure): + with pytest.raises(ValueError, match="diagonal supercell matrix"): + FiniteTemperaturePhononMaker().make( + si_structure, supercell_matrix=[[1, 1, 0], [0, 1, 0], [0, 0, 1]] + ) + + +def test_relax_md_and_static_settings_match(nio_supercell, test_dir): + """The relaxation, the MD and the statics must see the same Hamiltonian.""" + maker = FiniteTemperaturePhononMaker( + md_time=4.0, md_time_step=0.5, thermostat="nose-hoover" + ) + md_generator = maker.get_md_maker(maker.get_md_steps()[0]).input_set_generator + static_generator = maker.phonon_displacement_maker.input_set_generator + relax_generator = maker.bulk_relax_maker.relax_maker1.input_set_generator + # alternating magnetic moments on Ni, as carried over from a relaxation + magmoms = [2.0, -2.0] * 13 + [2.0] + [0.0] * 27 + magnetic = nio_supercell.copy(site_properties={"magmom": magmoms}) + sets = { + "md": md_generator.get_input_set(magnetic, potcar_spec=True), + "static": static_generator.get_input_set(magnetic, potcar_spec=True), + "relax": relax_generator.get_input_set(magnetic, potcar_spec=True), + } + for key in ("GGA", "ISPIN", "MAGMOM", "LDAU", "LDAUU", "LDAUL", "ISMEAR", "SIGMA"): + values = {name: input_set.incar.get(key) for name, input_set in sets.items()} + assert values["md"] == values["static"] == values["relax"], key + assert sets["md"].incar["GGA"] == "Ps" + assert sets["md"].incar["LDAUU"] == [6.2, 0] + assert sets["md"].incar["MAGMOM"] == magmoms + assert sets["md"].kpoints.kpts == sets["static"].kpoints.kpts + assert sets["relax"].kpoints.kpts == sets["static"].kpoints.kpts + for key in ("ENCUT", "ENAUG", "PREC"): + assert sets["relax"].incar[key] == sets["static"].incar[key], key + + # 4 ps at 0.5 fs, every step written to XDATCAR + md_incar, static_incar = sets["md"].incar, sets["static"].incar + assert md_incar["NSW"] == 8000 + assert md_incar["POTIM"] == 0.5 + assert md_incar["NBLOCK"] == 1 + assert md_incar["TEBEG"] == md_incar["TEEND"] == 300 + assert (md_incar["IBRION"], md_incar["ISIF"]) == (0, 2) + assert (md_incar["MDALGO"], md_incar["SMASS"]) == (2, 0) + assert md_incar["ENCUT"] == 500 + assert static_incar["ENCUT"] == 600 + assert static_incar["EDIFF"] == 1e-7 + assert static_incar["NSW"] == 0 + + # the same holds after a previous calculation sets the spin + si = Structure.from_file(test_dir / "structures" / "Si.cif") * (3, 3, 3) + prev_dir = test_dir / "vasp" / "Si_pheasy" / "tight_relax_2" / "outputs" + prev_sets = { + name: generator.get_input_set(si, prev_dir=prev_dir, potcar_spec=True) + for name, generator in ( + ("md", md_generator), + ("static", static_generator), + ("relax", relax_generator), + ) + } + for key in ("ISPIN", "MAGMOM", "ISMEAR", "SIGMA"): + assert prev_sets["md"].incar.get(key) == prev_sets["static"].incar.get(key) + assert prev_sets["md"].kpoints.kpts == prev_sets["static"].kpoints.kpts + assert prev_sets["relax"].kpoints.kpts == prev_sets["static"].kpoints.kpts + + +def test_langevin_thermostat(nio_supercell): + # Langevin by default + maker = FiniteTemperaturePhononMaker(temperature=600) + md_maker = maker.get_md_maker(100) + assert isinstance(md_maker.input_set_generator, LangevinMDSetGenerator) + incar = md_maker.input_set_generator.get_input_set( + nio_supercell, potcar_spec=True + ).incar + assert incar["MDALGO"] == 3 + assert incar["LANGEVIN_GAMMA"] == [10.0, 10.0] + assert incar["TEBEG"] == incar["TEEND"] == 600 + assert incar["ENCUT"] == 500 + # the maker of the flow is not changed + assert type(maker.md_maker.input_set_generator) is MDSetGenerator + + # a Langevin MD maker needs the Langevin thermostat + maker = FiniteTemperaturePhononMaker(md_maker=md_maker, thermostat="nose-hoover") + with pytest.raises(ValueError, match="Set thermostat to 'langevin'"): + maker.get_md_maker(100) + maker = FiniteTemperaturePhononMaker(md_maker=md_maker, thermostat="langevin") + assert maker.get_md_maker(100).input_set_generator.nsteps == 100 + + +def test_md_maker_checks(): + maker = FiniteTemperaturePhononMaker( + md_maker=MDMaker( + input_set_generator=MDSetGenerator(user_incar_settings={"POTIM": 2}) + ) + ) + with pytest.raises(ValueError, match=r"The flow sets \['POTIM'\]"): + maker.get_md_maker(10) + + maker = FiniteTemperaturePhononMaker( + md_maker=MDMaker( + input_set_generator=MDSetGenerator(user_incar_settings={"NBLOCK": 10}) + ) + ) + with pytest.raises(ValueError, match="NBLOCK must be 1"): + maker.get_md_maker(10) + + with pytest.raises(TypeError, match="VASP MDMaker with an MDSetGenerator"): + FiniteTemperaturePhononMaker(md_maker=StaticMaker()).get_md_maker(10) + + +def test_npt_maker(si_structure, nio_supercell): + npt_maker = MDMaker( + input_set_generator=MDSetGenerator(user_incar_settings={"ENCUT": 700}) + ) + maker = FiniteTemperaturePhononMaker( + temperature=600, pressure=5.0, npt_maker=npt_maker + ) + incar = ( + maker.get_npt_maker(100) + .input_set_generator.get_input_set(nio_supercell, potcar_spec=True) + .incar + ) + assert incar["ISIF"] == 3 + assert incar["MDALGO"] == 3 + assert incar["PSTRESS"] == 5.0 + assert incar["TEBEG"] == incar["TEEND"] == 600 + assert incar["NSW"] == 100 + assert incar["ENCUT"] == 700 + flow = maker.make(si_structure) + assert any(job.name.endswith(" NPT") for job in flow.jobs) + + maker = FiniteTemperaturePhononMaker( + npt_maker=MDMaker( + input_set_generator=MDSetGenerator(user_incar_settings={"PSTRESS": 1}) + ) + ) + with pytest.raises(ValueError, match="The flow sets PSTRESS"): + maker.get_npt_maker(10) + maker = FiniteTemperaturePhononMaker( + npt_maker=MDMaker(input_set_generator=LangevinMDSetGenerator()) + ) + with pytest.raises(TypeError, match="VASP MDMaker with an MDSetGenerator"): + maker.get_npt_maker(10) + + # a force field MD gets a force field NPT MD + maker = MLFFMDVaspStaticFiniteTemperaturePhononMaker(npt_maker=ForceFieldMDMaker()) + assert maker.get_npt_maker(10).pressure == 0.0 + + +def test_from_force_field_name(): + maker = VaspMDMLFFStaticFiniteTemperaturePhononMaker.from_force_field_name( + "MACE-MP-0", + calculator_kwargs={"model": "medium-omat-0"}, + temperature=500.0, + phonon_displacement_maker=StaticMaker(), + ) + assert maker.name == ("VASP MD MACE_MP_0 Static Finite Temperature Phonon Maker") + assert maker.temperature == 500.0 + assert isinstance(maker.md_maker, MDMaker) + # the maker built for the force field replaces the one given + displacement_maker = maker.phonon_displacement_maker + assert isinstance(displacement_maker, ForceFieldStaticMaker) + assert displacement_maker.calculator_kwargs["model"] == "medium-omat-0" + assert maker.born_maker is None + + maker = MLFFMDVaspStaticFiniteTemperaturePhononMaker.from_force_field_name( + "MACE-MP-0" + ) + assert maker.name == "MACE_MP_0 MD VASP Static Finite Temperature Phonon Maker" + assert isinstance(maker.md_maker, ForceFieldMDMaker) + assert (maker.md_code, maker.code) == ("forcefields", "vasp") + assert type(maker.born_maker).__name__ == "DielectricMaker" diff --git a/tests/vasp/test_sets.py b/tests/vasp/test_sets.py index 483ee027c6..bce3610021 100644 --- a/tests/vasp/test_sets.py +++ b/tests/vasp/test_sets.py @@ -9,6 +9,7 @@ HSERelaxSetGenerator, HSEStaticSetGenerator, HSETightRelaxSetGenerator, + LangevinMDSetGenerator, LobsterTightStaticSetGenerator, MDSetGenerator, NonSCFSetGenerator, @@ -337,3 +338,24 @@ def test_md_set_generator_sorts_structure(): n_types_in_poscar = len(vasp_input["POSCAR"].natoms) n_langevin_gamma = len(vasp_input["INCAR"]["LANGEVIN_GAMMA"]) assert n_types_in_poscar == n_langevin_gamma + + +def test_langevin_md_set_generator(): + """LANGEVIN_GAMMA has one value for each species block of the POSCAR.""" + structure = Structure( + lattice=Lattice.cubic(10), + species=["Al", "Cl", "Al", "Cl", "O", "Li"], + coords=[[0.1 * i, 0.1 * i, 0.1 * i] for i in range(6)], + ) + input_gen = LangevinMDSetGenerator(start_temp=600, end_temp=600, langevin_gamma=5) + vasp_input = input_gen.get_input_set(structure, potcar_spec=True) + incar = vasp_input["INCAR"] + assert incar["MDALGO"] == 3 + assert incar["LANGEVIN_GAMMA"] == [5] * len(vasp_input["POSCAR"].natoms) + assert incar["TEBEG"] == incar["TEEND"] == 600 + assert "SMASS" not in incar + + with pytest.raises(ValueError, match="Expect ensemble='nvt'"): + LangevinMDSetGenerator(ensemble="npt").get_input_set( + structure, potcar_spec=True + ) diff --git a/tutorials/finite_temperature_phonons.ipynb b/tutorials/finite_temperature_phonons.ipynb new file mode 100644 index 0000000000..5726ffe022 --- /dev/null +++ b/tutorials/finite_temperature_phonons.ipynb @@ -0,0 +1,2277 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Finite-temperature phonons with pheasy and machine-learned potentials\n", + "\n", + "This notebook runs `ForceFieldFiniteTemperaturePhononMaker` with three MACE\n", + "potentials. The workflow fits effective harmonic force constants to MD\n", + "snapshots at a given temperature. The workflow is new and less tested than the\n", + "other phonon workflows, so its defaults may still change.\n", + "\n", + "## Contents\n", + "\n", + "1. Five materials from 100 to 900 K\n", + "2. Thermal expansion with an NPT MD\n", + "3. MD length and number of snapshots\n", + "4. Li3PS4, a crystal that is unstable at 0 K\n", + "5. Li3PS4 with hopping Li: diffusion, DOS and free energies from the MD\n", + "6. Known limitations\n", + "\n", + "The long runs are switched off. Our results are shown as figures and tables." + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "## How the workflow works\n", + "\n", + "1. Relax the structure (optional).\n", + "2. Build a diagonal supercell with lattice vectors of at least `min_length`,\n", + " 12 Å by default.\n", + "3. Optional: run an NPT MD to get the cell at the temperature (Section 2).\n", + "4. Run an NVT MD, 8 ps with a 1 fs time step by default.\n", + "5. Leave out the first 1 ps and pick 50 evenly spaced snapshots.\n", + "6. Compute the forces on the snapshots and on the undisplaced supercell.\n", + "7. Fit the second-order force constants with pheasy (LASSO).\n", + "8. Compute the band structure and the DOS. Imaginary modes are reported, not\n", + " removed.\n", + "\n", + "A check of the trajectory flags melting, a shift away from the starting\n", + "structure and a drift of the energy. It raises a warning, and the fit is still\n", + "done. The idea is that of TDEP (Hellman et al., Phys. Rev. B 84, 180301\n", + "(2011)). The force constants include the effect of the anharmonic forces on the\n", + "frequencies. They give no lifetimes, and the MD is classical. The VASP version\n", + "is `atomate2.vasp.flows.finite_temperature_phonons.FiniteTemperaturePhononMaker`." + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "## Thermostat\n", + "\n", + "The NVT MD uses a Langevin thermostat with a friction of 10 ps⁻¹. It samples\n", + "the canonical ensemble for every mode (Bussi and Parrinello, Phys. Rev. E 75,\n", + "056707 (2007)). `thermostat=\"nose-hoover\"` is also available, but a weakly\n", + "coupled Nose-Hoover thermostat is not ergodic for a harmonic oscillator (Legoll\n", + "et al., Arch. Ration. Mech. Anal. 184, 449 (2007)). The NPT MD uses ASE's\n", + "`MTKNPT`. `random_seed` seeds the initial velocities and the Langevin forces." + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "## Installation\n", + "\n", + "```\n", + "pip install 'atomate2[pheasy,ase]'\n", + "pip install 'mace-torch>=0.3.16' mp-api\n", + "```\n", + "\n", + "The `pheasy` extra builds ALM from source. The hiPhive tutorial lists what the\n", + "build needs." + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "## Potentials\n", + "\n", + "- MACE-OMAT-0-medium. `mace_mp(model=\"medium-omat-0\")` downloads it.\n", + "- MACE-MATPES-PBE-0 and MACE-MATPES-r2SCAN-0. Download\n", + " [MACE-matpes-pbe-omat-ft.model](https://github.com/ACEsuit/mace-foundations/releases/download/mace_matpes_0/MACE-matpes-pbe-omat-ft.model)\n", + " and\n", + " [MACE-matpes-r2scan-omat-ft.model](https://github.com/ACEsuit/mace-foundations/releases/download/mace_matpes_0/MACE-matpes-r2scan-omat-ft.model)\n", + " and set their paths in the next cell.\n", + "\n", + "All three use the Academic Software License. Section 1 only needs\n", + "MACE-OMAT-0-medium. We ran in `float64` on one NVIDIA A100 GPU." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import warnings\n", + "\n", + "# macOS: conda's llvm-openmp and torch's bundled libomp both load and the\n", + "# duplicate aborts the process. This must be set before torch is imported.\n", + "os.environ.setdefault(\"KMP_DUPLICATE_LIB_OK\", \"TRUE\")\n", + "# hide the warnings of torch, e3nn and MACE when the models are loaded\n", + "warnings.filterwarnings(\"ignore\", module=r\"(torch|e3nn|mace)(\\.|$)\")\n", + "\n", + "import torch # noqa: E402\n", + "from mace.calculators import mace_mp # noqa: E402\n", + "\n", + "# Replace the two paths with the files you downloaded.\n", + "MODEL_FILES = {\n", + " \"OMAT-0-medium\": \"medium-omat-0\",\n", + " \"MATPES-PBE-0\": \"/path/to/MACE-matpes-pbe-omat-ft.model\",\n", + " \"MATPES-r2SCAN-0\": \"/path/to/MACE-matpes-r2scan-omat-ft.model\",\n", + "}\n", + "DEVICE = \"cuda\" if torch.cuda.is_available() else \"cpu\"\n", + "\n", + "\n", + "def calc_kwargs(model: str) -> dict:\n", + " \"\"\"Return the calculator settings for one of the models above.\"\"\"\n", + " return {\"model\": MODEL_FILES[model], \"device\": DEVICE, \"default_dtype\": \"float64\"}\n", + "\n", + "\n", + "_ = mace_mp(**calc_kwargs(\"OMAT-0-medium\")) # downloads and caches on first use" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "## One calculator for all jobs\n", + "\n", + "`run_locally` runs all jobs in one process, and each force field job loads its\n", + "own calculator. The next cell makes them share one, as in the thermal expansion\n", + "tutorial. A workflow manager does not need this." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "from ase.calculators.calculator import Calculator\n", + "\n", + "import atomate2.forcefields.utils as ff_utils\n", + "\n", + "_ase_calculator = ff_utils.ase_calculator\n", + "_calculators: dict[str, Calculator] = {}\n", + "\n", + "\n", + "def shared_ase_calculator(calculator_meta: object, **kwargs: object) -> Calculator:\n", + " \"\"\"Build each calculator once and reuse it for every job.\"\"\"\n", + " key = repr((calculator_meta, sorted(kwargs.items())))\n", + " if key not in _calculators:\n", + " _calculators[key] = _ase_calculator(calculator_meta, **kwargs)\n", + " _calculators[key].reset()\n", + " return _calculators[key]\n", + "\n", + "\n", + "ff_utils.ase_calculator = shared_ase_calculator" + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "## 1. Five materials from 100 to 900 K\n", + "\n", + "We first run Ca3Ir4Sn13 (mp-1200211) at 300 K with MACE-OMAT-0-medium. It has\n", + "40 atoms in the primitive cell and imaginary modes at 0 K. We take its PBEsol\n", + "cell from the Materials Project's Harmonic Phonon Database. Set\n", + "`MP_API_KEY` first." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "from mp_api.client import MPRester\n", + "from pymatgen.core import Structure\n", + "\n", + "if not os.environ.get(\"MP_API_KEY\"):\n", + " raise OSError(\"export MP_API_KEY before running this cell\")\n", + "\n", + "\n", + "def get_pheasy_structure(mp_id: str) -> Structure:\n", + " \"\"\"Get the PBEsol cell of a material in the Harmonic Phonon Database.\"\"\"\n", + " with MPRester() as mpr:\n", + " (doc,) = mpr.materials.phonon.search(\n", + " material_ids=[mp_id],\n", + " phonon_method=\"pheasy\",\n", + " fields=[\"structure\"],\n", + " num_chunks=1,\n", + " chunk_size=10,\n", + " )\n", + " return doc.structure\n", + "\n", + "\n", + "structure = get_pheasy_structure(\"mp-1200211\")\n", + "structure.lattice" + ] + }, + { + "cell_type": "markdown", + "id": "10", + "metadata": {}, + "source": [ + "## Building the workflow\n", + "\n", + "`from_force_field_name` uses one force field for every job. \"MACE-MP-0\" picks\n", + "the MACE calculator, and `calculator_kwargs` picks the model. The cell is\n", + "already relaxed with PBEsol, so we skip the relaxation. We use the 2x2x2\n", + "supercell (320 atoms) of our AIMD reference." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "from atomate2.forcefields.flows.finite_temperature_phonons import (\n", + " ForceFieldFiniteTemperaturePhononMaker,\n", + ")\n", + "\n", + "maker = ForceFieldFiniteTemperaturePhononMaker.from_force_field_name(\n", + " \"MACE-MP-0\",\n", + " calculator_kwargs=calc_kwargs(\"OMAT-0-medium\"),\n", + " relax_initial_structure=False,\n", + " temperature=300,\n", + ")\n", + "flow = maker.make(structure, supercell_matrix=[[2, 0, 0], [0, 2, 0], [0, 0, 2]])\n", + "flow.draw_graph().show()" + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "## Running the workflow\n", + "\n", + "The flow has one MD job and 51 force calculations. It took 12 minutes on one\n", + "A100 GPU. `create_folders=True` is needed, since the snapshots are read from the\n", + "trajectory file of the MD job." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "from jobflow import JobStore, run_locally\n", + "from maggma.stores import MemoryStore\n", + "\n", + "job_store = JobStore(MemoryStore(), additional_stores={\"data\": MemoryStore()})\n", + "responses = run_locally(flow, store=job_store, create_folders=True, ensure_success=True)\n", + "doc = responses[flow.output.uuid][1].output" + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": {}, + "source": [ + "The output is a `FiniteTemperaturePhononDoc`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "health = doc.trajectory_health\n", + "{\n", + " \"highest frequency in THz\": float(np.max(doc.phonon_bandstructure.bands)),\n", + " \"imaginary modes on the band path\": doc.has_imaginary_modes,\n", + " \"lowest frequency at the commensurate q-points in THz\": doc.lowest_frequency,\n", + " \"LASSO penalty\": doc.lasso_alpha,\n", + " \"force RMSE in eV/A\": doc.force_rmse,\n", + " \"trajectory check\": health.verdict,\n", + " \"Lindemann ratio\": health.lindemann_ratio,\n", + "}" + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "Our run gave a highest frequency of 5.97 THz, no imaginary modes and a stable\n", + "trajectory. `has_imaginary_modes` checks the band path. `lowest_frequency` uses\n", + "only the q-points of the supercell, where the fit sets the frequencies directly\n", + "(Section 3). The check reports melting above a Lindemann ratio of 0.15, the\n", + "value at melting of fcc solids (Saija et al., J. Chem. Phys. 124, 244504\n", + "(2006))." + ] + }, + { + "cell_type": "markdown", + "id": "17", + "metadata": {}, + "source": [ + "## Comparing with 0 K\n", + "\n", + "The next cell plots our 300 K bands next to the 0 K PBEsol bands of the\n", + "Materials Project's Harmonic Phonon Database (Sahasrabuddhe et al., ChemRxiv\n", + "(2026), doi:10.26434/chemrxiv.15004632/v1)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18", + "metadata": {}, + "outputs": [], + "source": [ + "import re\n", + "\n", + "import matplotlib.pyplot as plt\n", + "from pymatgen.phonon.bandstructure import PhononBandStructureSymmLine\n", + "from pymatgen.phonon.dos import PhononDos\n", + "\n", + "\n", + "def plot_bands(\n", + " ax: plt.Axes,\n", + " bs: PhononBandStructureSymmLine,\n", + " label: str | None = None,\n", + " distance: list[float] | None = None,\n", + " **style,\n", + ") -> None:\n", + " \"\"\"Plot a phonon band structure, one segment of the path at a time.\n", + "\n", + " distance replaces the distances along the path of bs, for a band structure\n", + " computed on the same q-points as another one.\n", + " \"\"\"\n", + " distance = np.array(bs.distance if distance is None else distance)\n", + " positions, names = [], []\n", + " for idx, branch in enumerate(bs.branches):\n", + " part = slice(branch[\"start_index\"], branch[\"end_index\"] + 1)\n", + " lines = ax.plot(distance[part], bs.bands[:, part].T, **style)\n", + " if idx == 0 and label:\n", + " lines[0].set_label(label)\n", + " start, end = branch[\"name\"].split(\"-\")\n", + " if positions and np.isclose(positions[-1], distance[part.start]):\n", + " if names[-1] != start:\n", + " names[-1] += \"|\" + start\n", + " else:\n", + " positions.append(distance[part.start])\n", + " names.append(start)\n", + " positions.append(distance[branch[\"end_index\"]])\n", + " names.append(end)\n", + " names = [name.replace(\"\\\\Gamma\", \"Γ\").replace(\"GAMMA\", \"Γ\") for name in names]\n", + " names = [re.sub(r\"_(\\w+)\", r\"$_{\\1}$\", name) for name in names]\n", + " ax.set_xticks(positions, names)\n", + " ax.set_xlim(distance[0], distance[-1])\n", + " ax.axhline(0, color=\"gray\", lw=0.5)\n", + "\n", + "\n", + "def plot_dos(ax: plt.Axes, dos: PhononDos, label: str | None = None, **style) -> None:\n", + " \"\"\"Plot a phonon DOS sideways, to sit next to a band structure.\"\"\"\n", + " ax.plot(dos.densities, dos.frequencies, label=label, **style)\n", + " ax.set_xlabel(\"DOS\")\n", + "\n", + "\n", + "with MPRester() as mpr:\n", + " bs_0k = mpr.materials.phonon.get_bandstructure_from_material_id(\n", + " \"mp-1200211\", phonon_method=\"pheasy\"\n", + " ).to_pmg\n", + " dos_0k = mpr.materials.phonon.get_dos_from_material_id(\n", + " \"mp-1200211\", phonon_method=\"pheasy\"\n", + " ).to_pmg\n", + "\n", + "fig, axes = plt.subplots(\n", + " 1, 3, figsize=(13, 4), sharey=True, gridspec_kw={\"width_ratios\": [3, 3, 1]}\n", + ")\n", + "plot_bands(axes[0], bs_0k, color=\"gray\", lw=0.8)\n", + "plot_bands(axes[1], doc.phonon_bandstructure, color=\"C0\", lw=0.8)\n", + "plot_dos(axes[2], dos_0k, label=\"0 K\", color=\"gray\")\n", + "plot_dos(axes[2], doc.phonon_dos, label=\"300 K\", color=\"C0\")\n", + "axes[0].set_title(\"0 K, PBEsol\")\n", + "axes[1].set_title(\"300 K, MACE-OMAT-0-medium\")\n", + "axes[0].set_ylabel(\"Frequency (THz)\")\n", + "axes[2].legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "19", + "metadata": {}, + "source": [ + "![Band structure and DOS of Ca3Ir4Sn13 at 0 K and at 300 K](finite_temperature_phonons_figures/s1_ca3ir4sn13_300k.png)\n", + "\n", + "The imaginary modes of 0 K, down to -0.75 THz, are gone at 300 K. AIMD with\n", + "PBEsol in the same supercell gives a highest frequency of 6.03 THz. Over the\n", + "band path, MACE-OMAT-0-medium differs from it by 0.07 THz on average." + ] + }, + { + "cell_type": "markdown", + "id": "20", + "metadata": {}, + "source": [ + "## All temperatures, potentials and materials\n", + "\n", + "The loop runs five materials with three potentials at five temperatures. AlAsPt5,\n", + "Zn3P2 and Er5Tl3 also have imaginary modes at 0 K, and Bi4S3N2 is polar. The\n", + "loop is off by default. Each run took 6 to 22 minutes." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21", + "metadata": {}, + "outputs": [], + "source": [ + "from atomate2.common.schemas.finite_temperature_phonons import (\n", + " FiniteTemperaturePhononDoc,\n", + ")\n", + "\n", + "# mp-id and the diagonal of the supercell matrix of our AIMD reference\n", + "MATERIALS = {\n", + " \"Ca3Ir4Sn13\": (\"mp-1200211\", (2, 2, 2)),\n", + " \"AlAsPt5\": (\"mp-1025306\", (4, 4, 2)),\n", + " \"Zn3P2\": (\"mp-2071\", (2, 2, 2)),\n", + " \"Bi4S3N2\": (\"mp-1245549\", (2, 2, 2)),\n", + " \"Er5Tl3\": (\"mp-1105965\", (2, 2, 2)),\n", + " \"KNaICl\": (\"mp-1002081\", (3, 3, 2)),\n", + "}\n", + "TEMPERATURES = (100, 300, 500, 700, 900)\n", + "\n", + "\n", + "def run_ft(\n", + " structure: Structure,\n", + " supercell: tuple[int, int, int],\n", + " model: str,\n", + " root_dir: str,\n", + " temperature: float,\n", + " md_time: float = 8.0,\n", + " n_snapshots: int = 50,\n", + " relax: bool = False,\n", + " run_npt: bool = False,\n", + ") -> FiniteTemperaturePhononDoc:\n", + " \"\"\"Run the finite-temperature phonon workflow and return its output.\"\"\"\n", + " maker = ForceFieldFiniteTemperaturePhononMaker.from_force_field_name(\n", + " \"MACE-MP-0\",\n", + " calculator_kwargs=calc_kwargs(model),\n", + " relax_initial_structure=relax,\n", + " run_npt=run_npt,\n", + " temperature=temperature,\n", + " md_time=md_time,\n", + " n_snapshots=n_snapshots,\n", + " )\n", + " flow = maker.make(structure, supercell_matrix=np.diag(supercell).tolist())\n", + " os.makedirs(root_dir, exist_ok=True)\n", + " responses = run_locally(\n", + " flow,\n", + " store=JobStore(MemoryStore(), additional_stores={\"data\": MemoryStore()}),\n", + " create_folders=True,\n", + " root_dir=root_dir,\n", + " ensure_success=True,\n", + " )\n", + " return responses[flow.output.uuid][1].output\n", + "\n", + "\n", + "RUN_ALL = False\n", + "docs = {}\n", + "if RUN_ALL:\n", + " for name in (\"Ca3Ir4Sn13\", \"AlAsPt5\", \"Zn3P2\", \"Bi4S3N2\", \"Er5Tl3\"):\n", + " mp_id, supercell = MATERIALS[name]\n", + " unit_cell = get_pheasy_structure(mp_id)\n", + " for model in MODEL_FILES:\n", + " for temperature in TEMPERATURES:\n", + " docs[name, model, temperature] = run_ft(\n", + " unit_cell,\n", + " supercell,\n", + " model,\n", + " f\"runs/{name}_{model}_{temperature}K\",\n", + " temperature,\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "22", + "metadata": {}, + "source": [ + "The next cell plots the bands at 100, 500 and 900 K. For the two materials with\n", + "40 atoms in the primitive cell it also plots the DOS." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23", + "metadata": {}, + "outputs": [], + "source": [ + "def plot_temperatures(name: str, with_dos: bool, label: str | None = None) -> None:\n", + " \"\"\"Plot the band structures of one material at 100, 500 and 900 K.\n", + "\n", + " label picks the runs of Section 2, \"relaxed\" or \"npt\". None picks Section 1.\n", + " \"\"\"\n", + " suffix = () if label is None else (label,)\n", + " ratios = [3, 1] * len(MODEL_FILES) if with_dos else [3] * len(MODEL_FILES)\n", + " fig, axes = plt.subplots(\n", + " 1,\n", + " len(ratios),\n", + " figsize=(4.5 * len(MODEL_FILES), 4),\n", + " sharey=True,\n", + " gridspec_kw={\"width_ratios\": ratios},\n", + " )\n", + " colors = dict(\n", + " zip(TEMPERATURES, plt.cm.plasma(np.linspace(0, 0.85, 5)), strict=True)\n", + " )\n", + " step = 2 if with_dos else 1\n", + " for idx, model in enumerate(MODEL_FILES):\n", + " ax = axes[step * idx]\n", + " for temperature in (100, 500, 900):\n", + " plot_bands(\n", + " ax,\n", + " docs[(name, model, temperature, *suffix)].phonon_bandstructure,\n", + " label=f\"{temperature} K\",\n", + " color=colors[temperature],\n", + " lw=0.6,\n", + " )\n", + " ax.set_title(model)\n", + " if with_dos:\n", + " for temperature in TEMPERATURES:\n", + " plot_dos(\n", + " axes[step * idx + 1],\n", + " docs[(name, model, temperature, *suffix)].phonon_dos,\n", + " label=f\"{temperature} K\",\n", + " color=colors[temperature],\n", + " )\n", + " axes[0].set_ylabel(\"Frequency (THz)\")\n", + " axes[-1].legend(fontsize=8)\n", + " fig.suptitle(name if label is None else f\"{name}, {label} cell\", y=1.02)\n", + " plt.show()\n", + "\n", + "\n", + "if RUN_ALL:\n", + " for name in (\"Ca3Ir4Sn13\", \"AlAsPt5\", \"Zn3P2\", \"Bi4S3N2\", \"Er5Tl3\"):\n", + " plot_temperatures(name, with_dos=name in (\"Ca3Ir4Sn13\", \"Zn3P2\"))" + ] + }, + { + "cell_type": "markdown", + "id": "24", + "metadata": {}, + "source": [ + "The figures below show our runs.\n", + "\n", + "![Band structures of Ca3Ir4Sn13 at 100, 500 and 900 K](finite_temperature_phonons_figures/s1_bands_Ca3Ir4Sn13.png)\n", + "\n", + "![Band structures of AlAsPt5 at 100, 500 and 900 K](finite_temperature_phonons_figures/s1_bands_AlAsPt5.png)\n", + "\n", + "![Band structures of Zn3P2 at 100, 500 and 900 K](finite_temperature_phonons_figures/s1_bands_Zn3P2.png)\n", + "\n", + "![Band structures of Bi4S3N2 at 100, 500 and 900 K](finite_temperature_phonons_figures/s1_bands_Bi4S3N2.png)\n", + "\n", + "![Band structures of Er5Tl3 at 100, 500 and 900 K](finite_temperature_phonons_figures/s1_bands_Er5Tl3.png)" + ] + }, + { + "cell_type": "markdown", + "id": "25", + "metadata": {}, + "source": [ + "## Our results\n", + "\n", + "The dashed line is the highest frequency of AIMD at 300 K." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "26", + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "\n", + "# highest frequency on the band path in THz, at 100, 300, 500, 700 and 900 K\n", + "MAX_FREQUENCY = {\n", + " \"Ca3Ir4Sn13\": {\n", + " \"OMAT-0-medium\": (6.03, 5.97, 5.94, 5.89, 5.78),\n", + " \"MATPES-PBE-0\": (5.80, 5.75, 5.74, 5.66, 5.64),\n", + " \"MATPES-r2SCAN-0\": (6.03, 5.95, 5.91, 5.89, 5.78),\n", + " },\n", + " \"AlAsPt5\": {\n", + " \"OMAT-0-medium\": (7.80, 7.76, 7.70, 7.77, 7.78),\n", + " \"MATPES-PBE-0\": (7.58, 7.66, 7.68, 7.67, 7.80),\n", + " \"MATPES-r2SCAN-0\": (7.65, 7.58, 7.54, 7.68, 7.58),\n", + " },\n", + " \"Zn3P2\": {\n", + " \"OMAT-0-medium\": (10.73, 10.44, 10.05, 9.80, 9.42),\n", + " \"MATPES-PBE-0\": (10.42, 10.08, 9.91, 9.63, 9.27),\n", + " \"MATPES-r2SCAN-0\": (11.26, 10.91, 10.69, 10.33, 10.14),\n", + " },\n", + " \"Bi4S3N2\": {\n", + " \"OMAT-0-medium\": (17.35, 16.93, 15.98, 15.85, 15.12),\n", + " \"MATPES-PBE-0\": (15.64, 15.21, 15.28, 14.56, 14.53),\n", + " \"MATPES-r2SCAN-0\": (17.60, 17.23, 16.52, 16.00, 15.73),\n", + " },\n", + " \"Er5Tl3\": {\n", + " \"OMAT-0-medium\": (4.36, 4.25, 4.28, 4.24, 4.19),\n", + " \"MATPES-PBE-0\": (4.15, 4.12, 4.05, 3.91, 3.86),\n", + " \"MATPES-r2SCAN-0\": (4.18, 4.14, 4.15, 4.11, 3.98),\n", + " },\n", + "}\n", + "# AIMD reference at 300 K: highest frequency in THz, and the mean absolute\n", + "# difference of each potential from it over all branches of the band path\n", + "AIMD_MAX_300K = {\n", + " \"Ca3Ir4Sn13\": 6.03,\n", + " \"AlAsPt5\": 7.79,\n", + " \"Zn3P2\": 10.09,\n", + " \"Bi4S3N2\": 17.03,\n", + " \"Er5Tl3\": 3.97,\n", + "}\n", + "MAE_300K = {\n", + " \"Ca3Ir4Sn13\": {\n", + " \"OMAT-0-medium\": 0.07,\n", + " \"MATPES-PBE-0\": 0.17,\n", + " \"MATPES-r2SCAN-0\": 0.10,\n", + " },\n", + " \"AlAsPt5\": {\"OMAT-0-medium\": 0.11, \"MATPES-PBE-0\": 0.19, \"MATPES-r2SCAN-0\": 0.17},\n", + " \"Zn3P2\": {\"OMAT-0-medium\": 0.21, \"MATPES-PBE-0\": 0.10, \"MATPES-r2SCAN-0\": 0.44},\n", + " \"Bi4S3N2\": {\"OMAT-0-medium\": 0.15, \"MATPES-PBE-0\": 0.46, \"MATPES-r2SCAN-0\": 0.16},\n", + " \"Er5Tl3\": {\"OMAT-0-medium\": 0.18, \"MATPES-PBE-0\": 0.17, \"MATPES-r2SCAN-0\": 0.13},\n", + "}\n", + "# runs with imaginary modes on the band path, and runs whose trajectory check\n", + "# was not \"stable\"\n", + "IMAGINARY = {\n", + " (\"AlAsPt5\", \"MATPES-PBE-0\"): (100,),\n", + " (\"Bi4S3N2\", \"OMAT-0-medium\"): (100, 300, 500),\n", + " (\"Bi4S3N2\", \"MATPES-PBE-0\"): (300, 500, 700, 900),\n", + " (\"Bi4S3N2\", \"MATPES-r2SCAN-0\"): (100, 300, 500, 700, 900),\n", + "}\n", + "NOT_STABLE = {\n", + " (\"Ca3Ir4Sn13\", \"OMAT-0-medium\", 100): \"shifted_or_diffusing\",\n", + " (\"Ca3Ir4Sn13\", \"MATPES-r2SCAN-0\", 100): \"shifted_or_diffusing\",\n", + " (\"Er5Tl3\", \"MATPES-r2SCAN-0\", 100): \"shifted_or_diffusing\",\n", + " (\"Zn3P2\", \"OMAT-0-medium\", 900): \"melted\",\n", + " (\"Zn3P2\", \"MATPES-PBE-0\", 900): \"melted\",\n", + "}\n", + "\n", + "fig, axes = plt.subplots(1, 5, figsize=(16, 3.4))\n", + "for ax, (name, by_model) in zip(axes, MAX_FREQUENCY.items(), strict=True):\n", + " for model, values in by_model.items():\n", + " ax.plot(TEMPERATURES, values, \"o-\", label=model)\n", + " ax.axhline(AIMD_MAX_300K[name], color=\"black\", ls=\"--\", lw=0.8)\n", + " ax.set_title(name)\n", + " ax.set_xlabel(\"Temperature (K)\")\n", + "axes[0].set_ylabel(\"Highest frequency (THz)\")\n", + "fig.tight_layout()\n", + "fig.legend(\n", + " *axes[0].get_legend_handles_labels(),\n", + " loc=\"upper center\",\n", + " ncol=3,\n", + " bbox_to_anchor=(0.5, 0),\n", + ")\n", + "plt.show()\n", + "\n", + "pd.DataFrame(MAE_300K).T.assign(**{\"AIMD highest frequency\": pd.Series(AIMD_MAX_300K)})" + ] + }, + { + "cell_type": "markdown", + "id": "27", + "metadata": {}, + "source": [ + "![Highest frequency against temperature](finite_temperature_phonons_figures/s1_max_frequency.png)\n", + "\n", + "| Material | OMAT-0-medium | MATPES-PBE-0 | MATPES-r2SCAN-0 | AIMD highest frequency (THz) |\n", + "|---|---|---|---|---|\n", + "| Ca3Ir4Sn13 | 0.07 | 0.17 | 0.10 | 6.03 |\n", + "| AlAsPt5 | 0.11 | 0.19 | 0.17 | 7.79 |\n", + "| Zn3P2 | 0.21 | 0.10 | 0.44 | 10.09 |\n", + "| Bi4S3N2 | 0.15 | 0.46 | 0.16 | 17.03 |\n", + "| Er5Tl3 | 0.18 | 0.17 | 0.13 | 3.97 |\n", + "\n", + "The table gives the mean absolute difference from AIMD at 300 K over the band\n", + "path, in THz.\n", + "\n", + "- No potential is best for all materials. MACE-OMAT-0-medium is closest for\n", + " three of the five.\n", + "- Most frequencies soften with temperature. The highest frequency of Bi4S3N2\n", + " drops by 2.2 THz from 100 to 900 K. AlAsPt5 barely changes.\n", + "- The imaginary modes of 0 K are gone for Ca3Ir4Sn13, Zn3P2 and Er5Tl3. Bi4S3N2\n", + " keeps some in 12 of 15 runs, and so does its AIMD reference.\n", + "- Zn3P2 melts at 900 K with two of the potentials. At 100 K three runs are\n", + " reported as shifted. Their atoms likely sit in a slightly distorted\n", + " structure." + ] + }, + { + "cell_type": "markdown", + "id": "28", + "metadata": {}, + "source": [ + "## 2. Thermal expansion with an NPT MD\n", + "\n", + "Each force field has its own 0 K cell, and the crystal expands as it heats up.\n", + "With `run_npt=True`, the workflow relaxes the structure and runs an NPT MD at\n", + "the temperature, 8 ps by default. It averages the cell after the first 2 ps,\n", + "keeping the symmetry, and relaxes the atoms in that cell. The NVT MD and the fit\n", + "then use this cell. We run Ca3Ir4Sn13 at 900 K." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "29", + "metadata": {}, + "outputs": [], + "source": [ + "maker = ForceFieldFiniteTemperaturePhononMaker.from_force_field_name(\n", + " \"MACE-MP-0\",\n", + " calculator_kwargs=calc_kwargs(\"OMAT-0-medium\"),\n", + " run_npt=True,\n", + " temperature=900,\n", + ")\n", + "flow = maker.make(structure, supercell_matrix=[[2, 0, 0], [0, 2, 0], [0, 0, 2]])\n", + "responses = run_locally(\n", + " flow,\n", + " store=JobStore(MemoryStore(), additional_stores={\"data\": MemoryStore()}),\n", + " create_folders=True,\n", + " ensure_success=True,\n", + ")\n", + "doc_npt = responses[flow.output.uuid][1].output\n", + "{\n", + " \"lattice at 0 K in A\": doc_npt.npt_input_structure.lattice.abc,\n", + " \"lattice at 900 K in A\": doc_npt.structure.lattice.abc,\n", + " \"volume ratio\": doc_npt.structure.volume / doc_npt.npt_input_structure.volume,\n", + " \"NPT trajectory check\": doc_npt.npt_trajectory_health.verdict,\n", + " \"highest frequency in THz\": float(np.max(doc_npt.phonon_bandstructure.bands)),\n", + " \"imaginary modes on the band path\": doc_npt.has_imaginary_modes,\n", + " \"trajectory check\": doc_npt.trajectory_health.verdict,\n", + "}" + ] + }, + { + "cell_type": "markdown", + "id": "30", + "metadata": {}, + "source": [ + "Our run took about 20 minutes. The 0 K cell of MACE-OMAT-0-medium is 3.8% larger\n", + "than the PBEsol cell, and the 900 K cell 4.3% larger again. Both trajectories\n", + "were stable. The highest frequency is 4.91 THz, against 5.78 THz in the PBEsol\n", + "cell." + ] + }, + { + "cell_type": "markdown", + "id": "31", + "metadata": {}, + "source": [ + "## All runs\n", + "\n", + "The loop runs the 75 settings of Section 1 again, in the force field's 0 K cell\n", + "and with the NPT MD. It is off by default." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "32", + "metadata": {}, + "outputs": [], + "source": [ + "RUN_NPT = False\n", + "if RUN_NPT:\n", + " for name in (\"Ca3Ir4Sn13\", \"AlAsPt5\", \"Zn3P2\", \"Bi4S3N2\", \"Er5Tl3\"):\n", + " mp_id, supercell = MATERIALS[name]\n", + " unit_cell = get_pheasy_structure(mp_id)\n", + " for model in MODEL_FILES:\n", + " for temperature in TEMPERATURES:\n", + " for label, run_npt in ((\"relaxed\", False), (\"npt\", True)):\n", + " docs[name, model, temperature, label] = run_ft(\n", + " unit_cell,\n", + " supercell,\n", + " model,\n", + " f\"runs/{name}_{model}_{temperature}K_{label}\",\n", + " temperature,\n", + " relax=True,\n", + " run_npt=run_npt,\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "33", + "metadata": {}, + "source": [ + "The next cell plots the bands in the NPT cell at 100, 500 and 900 K, as in\n", + "Section 1." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "34", + "metadata": {}, + "outputs": [], + "source": [ + "if RUN_NPT:\n", + " for name in (\"Ca3Ir4Sn13\", \"AlAsPt5\", \"Zn3P2\", \"Bi4S3N2\", \"Er5Tl3\"):\n", + " plot_temperatures(name, with_dos=name in (\"Ca3Ir4Sn13\", \"Zn3P2\"), label=\"npt\")" + ] + }, + { + "cell_type": "markdown", + "id": "35", + "metadata": {}, + "source": [ + "The figures below show our runs.\n", + "\n", + "![Band structures of Ca3Ir4Sn13 in the NPT cell at 100, 500 and 900 K](finite_temperature_phonons_figures/s2_bands_Ca3Ir4Sn13.png)\n", + "\n", + "![Band structures of AlAsPt5 in the NPT cell at 100, 500 and 900 K](finite_temperature_phonons_figures/s2_bands_AlAsPt5.png)\n", + "\n", + "![Band structures of Zn3P2 in the NPT cell at 100, 500 and 900 K](finite_temperature_phonons_figures/s2_bands_Zn3P2.png)\n", + "\n", + "![Band structures of Bi4S3N2 in the NPT cell at 100, 500 and 900 K](finite_temperature_phonons_figures/s2_bands_Bi4S3N2.png)\n", + "\n", + "![Band structures of Er5Tl3 in the NPT cell at 100, 500 and 900 K](finite_temperature_phonons_figures/s2_bands_Er5Tl3.png)\n", + "\n", + "- From 100 to 900 K the highest frequency drops by 0.3 to 3.2 THz in the NPT\n", + " cell, against at most 2.3 THz in the PBEsol cell of Section 1.\n", + "- AlAsPt5 barely changes in the PBEsol cell, but softens by 0.3 to 0.6 THz in\n", + " the NPT cell." + ] + }, + { + "cell_type": "markdown", + "id": "36", + "metadata": {}, + "source": [ + "## Our results\n", + "\n", + "The plots give the volume relative to the PBEsol cell, and the highest frequency\n", + "in three cells: the PBEsol cell, the force field's 0 K cell and the NPT cell." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "37", + "metadata": {}, + "outputs": [], + "source": [ + "# volume relative to the PBEsol cell: the force field's relaxed 0 K cell, then\n", + "# the NPT cell at 100, 300, 500, 700 and 900 K\n", + "VOLUME = {\n", + " \"Ca3Ir4Sn13\": {\n", + " \"OMAT-0-medium\": (1.0382, 1.0420, 1.0517, 1.0612, 1.0716, 1.0828),\n", + " \"MATPES-PBE-0\": (1.0432, 1.0475, 1.0568, 1.0665, 1.0774, 1.0890),\n", + " \"MATPES-r2SCAN-0\": (1.0192, 1.0227, 1.0311, 1.0398, 1.0491, 1.0591),\n", + " },\n", + " \"AlAsPt5\": {\n", + " \"OMAT-0-medium\": (1.0341, 1.0382, 1.0463, 1.0551, 1.0642, 1.0744),\n", + " \"MATPES-PBE-0\": (1.0175, 1.0229, 1.0313, 1.0400, 1.0501, 1.0501),\n", + " \"MATPES-r2SCAN-0\": (1.0021, 1.0061, 1.0133, 1.0210, 1.0302, 1.0392),\n", + " },\n", + " \"Zn3P2\": {\n", + " \"OMAT-0-medium\": (1.0533, 1.0585, 1.0681, 1.0780, 1.0884, 1.1012),\n", + " \"MATPES-PBE-0\": (1.0510, 1.0551, 1.0636, 1.0723, 1.0813, 1.0934),\n", + " \"MATPES-r2SCAN-0\": (1.0087, 1.0125, 1.0211, 1.0301, 1.0387, 1.0491),\n", + " },\n", + " \"Bi4S3N2\": {\n", + " \"OMAT-0-medium\": (1.1075, 1.1153, 1.1269, 1.1424, 1.1515, 1.1741),\n", + " \"MATPES-PBE-0\": (1.0457, 1.0594, 1.0764, 1.0967, 1.1515, 1.1625),\n", + " \"MATPES-r2SCAN-0\": (1.0292, 1.0381, 1.0502, 1.0656, 1.0807, 1.0861),\n", + " },\n", + " \"Er5Tl3\": {\n", + " \"OMAT-0-medium\": (1.0291, 1.0356, 1.0465, 1.0573, 1.0684, 1.0801),\n", + " \"MATPES-PBE-0\": (1.1224, 1.1283, 1.1431, 1.1589, 1.1744, 1.1939),\n", + " \"MATPES-r2SCAN-0\": (1.0646, 1.0697, 1.0806, 1.0925, 1.1034, 1.1164),\n", + " },\n", + "}\n", + "# highest frequency on the band path in THz, at 100, 300, 500, 700 and 900 K,\n", + "# in the force field's relaxed 0 K cell\n", + "MAX_RELAXED = {\n", + " \"Ca3Ir4Sn13\": {\n", + " \"OMAT-0-medium\": (5.58, 5.58, 5.51, 5.47, 5.40),\n", + " \"MATPES-PBE-0\": (5.37, 5.31, 5.26, 5.19, 5.17),\n", + " \"MATPES-r2SCAN-0\": (5.78, 5.76, 5.69, 5.59, 5.54),\n", + " },\n", + " \"AlAsPt5\": {\n", + " \"OMAT-0-medium\": (7.19, 7.07, 7.09, 7.20, 7.17),\n", + " \"MATPES-PBE-0\": (7.23, 7.24, 7.11, 7.09, 7.30),\n", + " \"MATPES-r2SCAN-0\": (7.52, 7.43, 7.42, 7.49, 7.45),\n", + " },\n", + " \"Zn3P2\": {\n", + " \"OMAT-0-medium\": (10.25, 9.82, 9.67, 9.39, 8.91),\n", + " \"MATPES-PBE-0\": (9.99, 9.73, 9.35, 9.06, 8.70),\n", + " \"MATPES-r2SCAN-0\": (11.19, 10.92, 10.46, 10.40, 9.84),\n", + " },\n", + " \"Bi4S3N2\": {\n", + " \"OMAT-0-medium\": (16.74, 15.80, 15.47, 14.80, 14.41),\n", + " \"MATPES-PBE-0\": (15.40, 15.02, 14.32, 14.55, 14.42),\n", + " \"MATPES-r2SCAN-0\": (17.27, 16.95, 16.22, 15.74, 15.42),\n", + " },\n", + " \"Er5Tl3\": {\n", + " \"OMAT-0-medium\": (4.10, 4.07, 4.06, 4.00, 4.01),\n", + " \"MATPES-PBE-0\": (3.70, 3.62, 3.56, 3.57, 3.40),\n", + " \"MATPES-r2SCAN-0\": (3.83, 3.83, 3.77, 3.78, 3.72),\n", + " },\n", + "}\n", + "# the same in the NPT cell\n", + "MAX_NPT = {\n", + " \"Ca3Ir4Sn13\": {\n", + " \"OMAT-0-medium\": (5.56, 5.39, 5.26, 5.11, 4.91),\n", + " \"MATPES-PBE-0\": (5.33, 5.19, 5.03, 4.89, 4.71),\n", + " \"MATPES-r2SCAN-0\": (5.72, 5.59, 5.44, 5.26, 5.06),\n", + " },\n", + " \"AlAsPt5\": {\n", + " \"OMAT-0-medium\": (7.18, 6.90, 6.90, 6.51, 6.88),\n", + " \"MATPES-PBE-0\": (7.15, 6.94, 6.90, 6.63, 6.64),\n", + " \"MATPES-r2SCAN-0\": (7.39, 7.27, 7.26, 7.11, 6.77),\n", + " },\n", + " \"Zn3P2\": {\n", + " \"OMAT-0-medium\": (10.23, 9.82, 9.31, 8.91, 8.42),\n", + " \"MATPES-PBE-0\": (9.97, 9.68, 9.19, 8.86, 8.19),\n", + " \"MATPES-r2SCAN-0\": (11.16, 10.72, 10.39, 9.87, 9.66),\n", + " },\n", + " \"Bi4S3N2\": {\n", + " \"OMAT-0-medium\": (16.39, 15.58, 15.15, 13.75, 13.63),\n", + " \"MATPES-PBE-0\": (15.50, 15.04, 14.74, 13.39, 12.34),\n", + " \"MATPES-r2SCAN-0\": (17.38, 16.63, 16.28, 15.33, 14.67),\n", + " },\n", + " \"Er5Tl3\": {\n", + " \"OMAT-0-medium\": (4.04, 3.99, 3.79, 3.77, 3.66),\n", + " \"MATPES-PBE-0\": (3.67, 3.56, 3.46, 3.32, 3.25),\n", + " \"MATPES-r2SCAN-0\": (3.78, 3.72, 3.64, 3.56, 3.44),\n", + " },\n", + "}\n", + "\n", + "fig, axes = plt.subplots(1, 5, figsize=(16, 3.4))\n", + "for ax, (name, by_model) in zip(axes, VOLUME.items(), strict=True):\n", + " for model, values in by_model.items():\n", + " ax.plot((0, *TEMPERATURES), values, \"o-\", label=model)\n", + " ax.axhline(1, color=\"black\", ls=\"--\", lw=0.8)\n", + " ax.set_title(name)\n", + " ax.set_xlabel(\"Temperature (K)\")\n", + "axes[0].set_ylabel(\"Volume / PBEsol volume\")\n", + "fig.tight_layout()\n", + "fig.legend(\n", + " *axes[0].get_legend_handles_labels(),\n", + " loc=\"upper center\",\n", + " ncol=3,\n", + " bbox_to_anchor=(0.5, 0),\n", + ")\n", + "plt.show()\n", + "\n", + "fig, axes = plt.subplots(3, 5, figsize=(16, 8.5), sharex=True)\n", + "for col, name in enumerate(VOLUME):\n", + " for row, model in enumerate(MODEL_FILES):\n", + " ax = axes[row, col]\n", + " for label, values, style in (\n", + " (\"PBEsol cell\", MAX_FREQUENCY, \"o-\"),\n", + " (\"relaxed 0 K cell\", MAX_RELAXED, \"s--\"),\n", + " (\"NPT cell\", MAX_NPT, \"^:\"),\n", + " ):\n", + " ax.plot(TEMPERATURES, values[name][model], style, label=label)\n", + " if row == 0:\n", + " ax.set_title(name)\n", + " if col == 0:\n", + " ax.set_ylabel(f\"{model}\\nHighest frequency (THz)\")\n", + " if row == 2:\n", + " ax.set_xlabel(\"Temperature (K)\")\n", + "fig.tight_layout()\n", + "fig.legend(\n", + " *axes[0, 0].get_legend_handles_labels(),\n", + " loc=\"upper center\",\n", + " ncol=3,\n", + " bbox_to_anchor=(0.5, 0),\n", + ")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "38", + "metadata": {}, + "source": [ + "![Volume against temperature](finite_temperature_phonons_figures/s2_volume.png)\n", + "\n", + "![Highest frequency in three cells](finite_temperature_phonons_figures/s2_max_frequency.png)\n", + "\n", + "- The 0 K cells of the force fields are up to 12% larger than the PBEsol cell.\n", + " MACE-MATPES-r2SCAN-0 is closest to PBEsol for four of the five materials.\n", + "- From 0 to 900 K the volume grows by 3 to 6.5%, and by 11% for Bi4S3N2 with\n", + " MACE-MATPES-PBE-0.\n", + "- The larger 0 K cell alone lowers the frequencies, by up to 0.5 THz at 300 K.\n", + "- Thermal expansion lowers them by 0.15 to 0.8 THz at 900 K, and by 2.1 THz for\n", + " Bi4S3N2 with MACE-MATPES-PBE-0.\n", + "- Zn3P2 and Bi4S3N2 melt at 900 K in the NPT runs. AlAsPt5 with\n", + " MACE-MATPES-PBE-0 is reported as shifted at 900 K, and its cell does not grow\n", + " from 700 to 900 K.\n", + "\n", + "To compare with DFT in a given cell, run in that cell as in Section 1. To\n", + "predict the phonons at a temperature with the force field alone, use\n", + "`run_npt=True`." + ] + }, + { + "cell_type": "markdown", + "id": "39", + "metadata": {}, + "source": [ + "## 3. MD length and number of snapshots\n", + "\n", + "We run KNaICl (mp-1002081, 72 atoms in the supercell) at 300 K with 4, 8 and 16 ps\n", + "of MD and 25, 50 and 100 snapshots, and compare with AIMD. KNaICl is polar, so\n", + "the comparison is without the non-analytical correction. The grid is off by\n", + "default." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "40", + "metadata": {}, + "outputs": [], + "source": [ + "RUN_GRID = False\n", + "if RUN_GRID:\n", + " mp_id, supercell = MATERIALS[\"KNaICl\"]\n", + " knaicl = get_pheasy_structure(mp_id)\n", + " for model in MODEL_FILES:\n", + " for md_time in (4, 8, 16):\n", + " for n_snapshots in (25, 50, 100):\n", + " docs[\"KNaICl\", model, md_time, n_snapshots] = run_ft(\n", + " knaicl,\n", + " supercell,\n", + " model,\n", + " f\"runs/KNaICl_{model}_{md_time}ps_{n_snapshots}\",\n", + " 300,\n", + " md_time,\n", + " n_snapshots,\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "41", + "metadata": {}, + "outputs": [], + "source": [ + "MD_TIMES = (4, 8, 16)\n", + "N_SNAPSHOTS = (25, 50, 100)\n", + "# by potential, rows 4, 8 and 16 ps, columns 25, 50 and 100 snapshots\n", + "MAE_GRID = { # mean absolute difference from AIMD at 300 K in THz\n", + " \"OMAT-0-medium\": (\n", + " (0.179, 0.181, 0.147),\n", + " (0.167, 0.170, 0.148),\n", + " (0.162, 0.171, 0.158),\n", + " ),\n", + " \"MATPES-PBE-0\": (\n", + " (0.158, 0.137, 0.133),\n", + " (0.129, 0.108, 0.120),\n", + " (0.142, 0.089, 0.105),\n", + " ),\n", + " \"MATPES-r2SCAN-0\": (\n", + " (0.188, 0.161, 0.162),\n", + " (0.191, 0.184, 0.136),\n", + " (0.152, 0.154, 0.169),\n", + " ),\n", + "}\n", + "MAX_GRID = { # highest frequency on the band path in THz\n", + " \"OMAT-0-medium\": ((5.48, 5.52, 5.66), (5.35, 5.56, 5.63), (5.57, 5.62, 5.61)),\n", + " \"MATPES-PBE-0\": ((5.20, 5.18, 5.13), (5.07, 5.22, 5.33), (5.05, 5.31, 5.24)),\n", + " \"MATPES-r2SCAN-0\": ((5.64, 5.51, 5.57), (5.32, 5.58, 5.42), (5.42, 5.53, 5.53)),\n", + "}\n", + "IMAGINARY_GRID = { # True if the band path has imaginary modes\n", + " \"OMAT-0-medium\": ((False, True, True), (True, True, False), (True, True, False)),\n", + " \"MATPES-PBE-0\": ((True, True, True), (False, True, False), (False, False, True)),\n", + " \"MATPES-r2SCAN-0\": ((True, True, True), (True, True, False), (True, True, True)),\n", + "}\n", + "\n", + "fig, axes = plt.subplots(1, 3, figsize=(13, 3.6))\n", + "for ax, model in zip(axes, MODEL_FILES, strict=True):\n", + " grid = np.array(MAE_GRID[model])\n", + " image = ax.imshow(grid, cmap=\"viridis\", vmin=0.08, vmax=0.20)\n", + " for i in range(3):\n", + " for j in range(3):\n", + " mark = \" *\" if IMAGINARY_GRID[model][i][j] else \"\"\n", + " ax.text(\n", + " j, i, f\"{grid[i, j]:.3f}{mark}\", ha=\"center\", va=\"center\", color=\"w\"\n", + " )\n", + " ax.set_xticks(range(3), N_SNAPSHOTS)\n", + " ax.set_yticks(range(3), [f\"{t} ps\" for t in MD_TIMES])\n", + " ax.set_xlabel(\"Number of snapshots\")\n", + " ax.set_title(model)\n", + "fig.colorbar(image, ax=axes, label=\"Mean absolute difference from AIMD (THz)\")\n", + "plt.show()\n", + "\n", + "pd.concat(\n", + " {\n", + " model: pd.DataFrame(\n", + " MAX_GRID[model], index=[f\"{t} ps\" for t in MD_TIMES], columns=N_SNAPSHOTS\n", + " )\n", + " for model in MODEL_FILES\n", + " },\n", + " axis=1,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "42", + "metadata": {}, + "source": [ + "![Difference from AIMD against MD length and snapshots](finite_temperature_phonons_figures/s3_knaicl_grid.png)\n", + "\n", + "| Potential | MD length | 25 snapshots | 50 snapshots | 100 snapshots |\n", + "|---|---|---|---|---|\n", + "| OMAT-0-medium | 4 ps | 5.48 | 5.52 | 5.66 |\n", + "| OMAT-0-medium | 8 ps | 5.35 | 5.56 | 5.63 |\n", + "| OMAT-0-medium | 16 ps | 5.57 | 5.62 | 5.61 |\n", + "| MATPES-PBE-0 | 4 ps | 5.20 | 5.18 | 5.13 |\n", + "| MATPES-PBE-0 | 8 ps | 5.07 | 5.22 | 5.33 |\n", + "| MATPES-PBE-0 | 16 ps | 5.05 | 5.31 | 5.24 |\n", + "| MATPES-r2SCAN-0 | 4 ps | 5.64 | 5.51 | 5.57 |\n", + "| MATPES-r2SCAN-0 | 8 ps | 5.32 | 5.58 | 5.42 |\n", + "| MATPES-r2SCAN-0 | 16 ps | 5.42 | 5.53 | 5.53 |\n", + "\n", + "A star marks a run with imaginary modes on the band path. The table gives the\n", + "highest frequency in THz.\n", + "\n", + "- The difference from AIMD is 0.09 to 0.19 THz, with no trend. It is the noise\n", + " of a single trajectory.\n", + "- 19 of the 27 runs have imaginary modes, mostly at q-points between those of\n", + " the supercell, where the frequencies are interpolated. AIMD has none. An\n", + " earlier round with a Nose-Hoover thermostat gave 12 of 27.\n", + "- More snapshots cost little with a force field, but they did not remove these\n", + " modes. The defaults stay at 8 ps and 50 snapshots, since with VASP each\n", + " snapshot is a DFT calculation." + ] + }, + { + "cell_type": "markdown", + "id": "43", + "metadata": {}, + "source": [ + "## 4. Li3PS4, a crystal that is unstable at 0 K\n", + "\n", + "In experiment Li3PS4 is γ up to 523 to 573 K, β up to 723 to 748 K and α above\n", + "(Homma et al., Solid State Ionics 182, 53 (2011); Kaup et al., J. Mater. Chem. A\n", + "8, 12446 (2020)). mp-985583 is one ordering of the partly occupied Li sites of\n", + "β, and it has imaginary modes at 0 K. We run it from 100 to 1000 K with each\n", + "potential, in the relaxed 0 K cell (NVT) and in the NPT cell. The supercell has\n", + "128 atoms." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "44", + "metadata": {}, + "outputs": [], + "source": [ + "li3ps4 = get_pheasy_structure(\"mp-985583\")\n", + "LI3PS4_SUPERCELL = (2, 2, 1) # 128 atoms\n", + "LI3PS4_TEMPERATURES = tuple(range(100, 1001, 100))\n", + "li3ps4.composition, li3ps4.get_space_group_info()" + ] + }, + { + "cell_type": "markdown", + "id": "45", + "metadata": {}, + "source": [ + "## Phonons at 0 K\n", + "\n", + "We use the phonopy workflow here, since the pheasy workflow refits the force\n", + "constants to remove imaginary modes. The relaxation keeps the symmetry." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "46", + "metadata": {}, + "outputs": [], + "source": [ + "from atomate2.common.schemas.phonons import PhononBSDOSDoc\n", + "from atomate2.forcefields.flows.phonons import PhononMaker\n", + "\n", + "\n", + "def li3ps4_harmonic(model: str) -> PhononBSDOSDoc:\n", + " \"\"\"Run the phonopy workflow at 0 K, in the cell relaxed by the potential.\"\"\"\n", + " maker = PhononMaker.from_force_field_name(\n", + " \"MACE-MP-0\", calculator_kwargs=calc_kwargs(model)\n", + " )\n", + " maker.bulk_relax_maker.fix_symmetry = True\n", + " flow = maker.make(li3ps4, supercell_matrix=np.diag(LI3PS4_SUPERCELL).tolist())\n", + " root_dir = f\"runs/Li3PS4_{model}_0K\"\n", + " os.makedirs(root_dir, exist_ok=True)\n", + " responses = run_locally(\n", + " flow,\n", + " store=JobStore(MemoryStore(), additional_stores={\"data\": MemoryStore()}),\n", + " create_folders=True,\n", + " root_dir=root_dir,\n", + " ensure_success=True,\n", + " )\n", + " return responses[flow.output.uuid][1].output\n", + "\n", + "\n", + "doc_0k = li3ps4_harmonic(\"OMAT-0-medium\")\n", + "float(np.min(doc_0k.phonon_bandstructure.bands))" + ] + }, + { + "cell_type": "markdown", + "id": "47", + "metadata": {}, + "source": [ + "Our run took 2 minutes. The lowest frequency is -2.59 THz. With\n", + "MACE-MATPES-PBE-0 it is -1.97 THz and with MACE-MATPES-r2SCAN-0 -3.48 THz." + ] + }, + { + "cell_type": "markdown", + "id": "48", + "metadata": {}, + "source": [ + "## One temperature\n", + "\n", + "We run 600 K with MACE-OMAT-0-medium and the NPT MD." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "49", + "metadata": {}, + "outputs": [], + "source": [ + "doc_li3ps4 = run_ft(\n", + " li3ps4,\n", + " LI3PS4_SUPERCELL,\n", + " \"OMAT-0-medium\",\n", + " \"runs/Li3PS4_OMAT-0-medium_600K_npt\",\n", + " 600,\n", + " relax=True,\n", + " run_npt=True,\n", + ")\n", + "{\n", + " \"volume at 0 K in A^3\": doc_li3ps4.npt_input_structure.volume,\n", + " \"volume at 600 K in A^3\": doc_li3ps4.structure.volume,\n", + " \"lowest frequency in THz\": float(np.min(doc_li3ps4.phonon_bandstructure.bands)),\n", + " \"trajectory check\": doc_li3ps4.trajectory_health.verdict,\n", + " \"Lindemann ratio\": doc_li3ps4.trajectory_health.lindemann_ratio,\n", + "}" + ] + }, + { + "cell_type": "markdown", + "id": "50", + "metadata": {}, + "source": [ + "Our run took 31 minutes. The volume grows from 664.8 to 698.5 ų. The lowest\n", + "frequency is -0.57 THz, and the check reports \"melted\" with a Lindemann ratio of\n", + "0.24. Section 5 shows that only the Li atoms move." + ] + }, + { + "cell_type": "markdown", + "id": "51", + "metadata": {}, + "source": [ + "## All runs\n", + "\n", + "The loop runs the 60 settings and the three 0 K runs. It is off by default. Each\n", + "run took 15 to 62 minutes." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "52", + "metadata": {}, + "outputs": [], + "source": [ + "RUN_LI3PS4 = False\n", + "li3ps4_docs = {}\n", + "if RUN_LI3PS4:\n", + " for model in MODEL_FILES:\n", + " li3ps4_docs[model, 0, \"harmonic\"] = li3ps4_harmonic(model)\n", + " for temperature in LI3PS4_TEMPERATURES:\n", + " for label, run_npt in ((\"nvt\", False), (\"npt\", True)):\n", + " li3ps4_docs[model, temperature, label] = run_ft(\n", + " li3ps4,\n", + " LI3PS4_SUPERCELL,\n", + " model,\n", + " f\"runs/Li3PS4_{model}_{temperature}K_{label}\",\n", + " temperature,\n", + " relax=True,\n", + " run_npt=run_npt,\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "53", + "metadata": {}, + "source": [ + "The next cell plots the bands at 0, 300 and 700 K. Every run keeps the Pnma\n", + "symmetry, so all runs of a potential share the band path of its 0 K run." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "54", + "metadata": {}, + "outputs": [], + "source": [ + "def plot_li3ps4_bands(model: str) -> None:\n", + " \"\"\"Plot the bands of one potential at 0, 300 and 700 K, in both cells.\"\"\"\n", + " reference = li3ps4_docs[model, 0, \"harmonic\"].phonon_bandstructure\n", + " _, axes = plt.subplots(1, 2, figsize=(14, 4), sharey=True)\n", + " for ax, label in zip(axes, (\"nvt\", \"npt\"), strict=True):\n", + " plot_bands(ax, reference, label=\"0 K\", color=\"gray\", lw=0.5)\n", + " for temperature, color in ((300, \"C0\"), (700, \"C3\")):\n", + " plot_bands(\n", + " ax,\n", + " li3ps4_docs[model, temperature, label].phonon_bandstructure,\n", + " label=f\"{temperature} K\",\n", + " distance=reference.distance,\n", + " color=color,\n", + " lw=0.5,\n", + " )\n", + " ax.set_title(f\"{model}, {label.upper()}\")\n", + " axes[0].set_ylabel(\"Frequency (THz)\")\n", + " axes[1].legend(fontsize=8)\n", + " plt.show()\n", + "\n", + "\n", + "if RUN_LI3PS4:\n", + " for model in MODEL_FILES:\n", + " plot_li3ps4_bands(model)" + ] + }, + { + "cell_type": "markdown", + "id": "55", + "metadata": {}, + "source": [ + "The figures below show our runs.\n", + "\n", + "![Band structures of Li3PS4 with OMAT-0-medium at 0, 300 and 700 K](finite_temperature_phonons_figures/s4_bands_OMAT-0-medium.png)\n", + "\n", + "![Band structures of Li3PS4 with MATPES-PBE-0 at 0, 300 and 700 K](finite_temperature_phonons_figures/s4_bands_MATPES-PBE-0.png)\n", + "\n", + "![Band structures of Li3PS4 with MATPES-r2SCAN-0 at 0, 300 and 700 K](finite_temperature_phonons_figures/s4_bands_MATPES-r2SCAN-0.png)" + ] + }, + { + "cell_type": "markdown", + "id": "56", + "metadata": {}, + "source": [ + "## Our results\n", + "\n", + "The lowest frequency, the Lindemann ratio and the NPT volume of each run." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "57", + "metadata": {}, + "outputs": [], + "source": [ + "# lowest frequency on the band path in THz and Lindemann ratio of the NVT MD,\n", + "# at 100 to 1000 K, in the relaxed 0 K cell (nvt) and in the NPT cell (npt)\n", + "LOWEST_FREQUENCY = {\n", + " \"nvt\": {\n", + " \"OMAT-0-medium\": (-0.8, -0.4, -0.9, -0.8, -0.4, -1.6, -1.1, -1.1, -1.5, -1.5),\n", + " \"MATPES-PBE-0\": (-0.2, -0.9, -1.0, -2.2, -1.4, -1.3, -1.7, -1.6, -1.5, -1.6),\n", + " \"MATPES-r2SCAN-0\": (-0.7, -0.3, -0.6, -0.3, -0.7, -1.1, -1.6, -1.5, -1.5, -1.4),\n", + " },\n", + " \"npt\": {\n", + " \"OMAT-0-medium\": (-0.7, -0.4, -0.7, -0.7, -0.5, -0.6, -1.4, -1.5, -1.3, -1.4),\n", + " \"MATPES-PBE-0\": (-0.4, -0.5, -1.4, -0.5, -1.3, -1.4, -1.3, -1.7, -1.6, -1.6),\n", + " \"MATPES-r2SCAN-0\": (-0.5, -0.6, -1.4, 0.0, -0.7, -1.0, -1.1, -1.3, -1.8, -1.6),\n", + " },\n", + "}\n", + "LINDEMANN = {\n", + " \"nvt\": {\n", + " \"OMAT-0-medium\": (0.06, 0.1, 0.12, 0.17, 0.18, 0.26, 0.28, 0.33, 0.51, 0.54),\n", + " \"MATPES-PBE-0\": (0.08, 0.13, 0.16, 0.2, 0.24, 0.36, 0.42, 0.41, 0.55, 0.84),\n", + " \"MATPES-r2SCAN-0\": (0.06, 0.1, 0.13, 0.17, 0.21, 0.26, 0.28, 0.36, 0.58, 0.45),\n", + " },\n", + " \"npt\": {\n", + " \"OMAT-0-medium\": (0.08, 0.11, 0.15, 0.17, 0.19, 0.24, 0.39, 0.36, 0.58, 0.77),\n", + " \"MATPES-PBE-0\": (0.08, 0.13, 0.18, 0.19, 0.23, 0.32, 0.33, 0.41, 0.61, 0.74),\n", + " \"MATPES-r2SCAN-0\": (0.07, 0.1, 0.13, 0.18, 0.21, 0.24, 0.27, 0.4, 0.52, 0.61),\n", + " },\n", + "}\n", + "# volume of the unit cell in A^3, relaxed at 0 K and from the NPT MD\n", + "VOLUME_0K = {\"OMAT-0-medium\": 664.7, \"MATPES-PBE-0\": 651.6, \"MATPES-r2SCAN-0\": 649.6}\n", + "NPT_VOLUME = {\n", + " \"OMAT-0-medium\": (670, 675, 680, 687, 693, 698, 702, 713, 776, 767),\n", + " \"MATPES-PBE-0\": (656, 658, 661, 660, 663, 662, 667, 672, 678, 692),\n", + " \"MATPES-r2SCAN-0\": (656, 661, 666, 669, 673, 673, 675, 683, 691, 699),\n", + "}\n", + "\n", + "fig, axes = plt.subplots(1, 3, figsize=(14, 3.6))\n", + "for idx, model in enumerate(MODEL_FILES):\n", + " for label, style in ((\"nvt\", \"-o\"), (\"npt\", \"--s\")):\n", + " kwargs = {\"color\": f\"C{idx}\", \"ms\": 3, \"label\": f\"{model}, {label}\"}\n", + " axes[0].plot(\n", + " LI3PS4_TEMPERATURES, LOWEST_FREQUENCY[label][model], style, **kwargs\n", + " )\n", + " axes[1].plot(LI3PS4_TEMPERATURES, LINDEMANN[label][model], style, **kwargs)\n", + " ratio = np.array(NPT_VOLUME[model]) / VOLUME_0K[model]\n", + " axes[2].plot(LI3PS4_TEMPERATURES, ratio, \"--s\", color=f\"C{idx}\", ms=3)\n", + "axes[1].axhline(0.15, color=\"gray\", ls=\":\")\n", + "axes[0].set_ylabel(\"Lowest frequency (THz)\")\n", + "axes[1].set_ylabel(\"Lindemann ratio\")\n", + "axes[2].set_ylabel(\"V(NPT) / V(0 K)\")\n", + "for ax in axes:\n", + " ax.set_xlabel(\"Temperature (K)\")\n", + "axes[0].legend(fontsize=7)\n", + "fig.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "58", + "metadata": {}, + "source": [ + "![Lowest frequency, Lindemann ratio and volume of Li3PS4](finite_temperature_phonons_figures/s4_runs.png)\n", + "\n", + "- The imaginary modes shrink with temperature but do not go away. The one\n", + " exception is MACE-MATPES-r2SCAN-0 at 400 K in the NPT cell.\n", + "- The Lindemann ratio passes 0.15 at 300 to 400 K, and the check reports\n", + " melting. Li3PS4 does not melt there in experiment.\n", + "- The NPT cell grows by 3 to 7% up to 800 K.\n", + "- MACE-MATPES-PBE-0 gives the softest bands, with 15.4 THz at 0 K against 17.8\n", + " and 17.9 THz." + ] + }, + { + "cell_type": "markdown", + "id": "59", + "metadata": {}, + "source": [ + "## 5. Li3PS4 with hopping Li\n", + "\n", + "At 0 K the imaginary modes come from the structure, since this ordering is not\n", + "at an energy minimum. Above about 500 K the Li atoms also start to hop, while\n", + "the PS4 framework stays solid. A fit around fixed sites reads this as very soft\n", + "springs. The MD itself needs no fixed sites. We get the diffusion, the\n", + "vibrational DOS and the free energy from the saved trajectories,\n", + "`md_trajectory.traj`, starting with the 600 K run above." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "60", + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "\n", + "from ase import Atoms\n", + "from ase.io import Trajectory\n", + "\n", + "SKIP = 1000 # frames left out at the start, as in the fit\n", + "TIME_STEP = 1e-3 # ps per frame\n", + "\n", + "\n", + "def nvt_trajectory(folder: str) -> list[Atoms]:\n", + " \"\"\"Get the frames of the NVT MD in a run folder, after the first SKIP.\"\"\"\n", + " for path in Path(folder).rglob(\"md_trajectory.traj\"):\n", + " traj = Trajectory(str(path))\n", + " if np.allclose(traj[0].cell.array, traj[-1].cell.array):\n", + " return [traj[i] for i in range(SKIP, len(traj))]\n", + " raise FileNotFoundError(f\"no NVT trajectory in {folder}\")\n", + "\n", + "\n", + "def unwrapped_positions(frames: list[Atoms]) -> np.ndarray:\n", + " \"\"\"Get the Cartesian positions of every frame, without jumps across the cell.\"\"\"\n", + " cell = frames[0].cell.array\n", + " frac = np.array([atoms.get_scaled_positions(wrap=False) for atoms in frames])\n", + " steps = np.diff(frac, axis=0)\n", + " steps -= np.round(steps)\n", + " frac = np.concatenate([frac[:1], frac[:1] + np.cumsum(steps, axis=0)])\n", + " return frac @ cell\n", + "\n", + "\n", + "frames_600k = nvt_trajectory(\"runs/Li3PS4_OMAT-0-medium_600K_npt\")" + ] + }, + { + "cell_type": "markdown", + "id": "61", + "metadata": {}, + "source": [ + "## Diffusion\n", + "\n", + "$D$ is the slope of the mean-square displacement divided by 6, fitted from 1 to\n", + "3 ps." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "62", + "metadata": {}, + "outputs": [], + "source": [ + "def diffusion(\n", + " frames: list[Atoms], max_lag: int = 3000, fit: tuple[int, int] = (1000, 3000)\n", + ") -> dict[str, tuple[np.ndarray, np.ndarray, float]]:\n", + " \"\"\"Get the MSD (A^2) against time (ps) and D (cm^2/s) of each element.\n", + "\n", + " The motion of the center of mass is removed. fit gives the range of lags in\n", + " frames where the slope of the MSD is fitted.\n", + " \"\"\"\n", + " masses = frames[0].get_masses()\n", + " symbols = np.array(frames[0].get_chemical_symbols())\n", + " pos = unwrapped_positions(frames)\n", + " pos -= ((pos * masses[None, :, None]).sum(1) / masses.sum())[:, None]\n", + " lags = np.arange(0, max_lag + 1, 10)\n", + " origins = np.arange(0, len(pos) - max_lag, 100)\n", + " msd = np.array(\n", + " [((pos[origins + lag] - pos[origins]) ** 2).sum(2).mean(0) for lag in lags]\n", + " )\n", + " keep = (lags >= fit[0]) & (lags <= fit[1])\n", + " out = {}\n", + " for element in (\"Li\", \"P\", \"S\"):\n", + " curve = msd[:, symbols == element].mean(1)\n", + " slope = np.polyfit(lags[keep] * TIME_STEP, curve[keep], 1)[0]\n", + " out[element] = (lags * TIME_STEP, curve, slope / 6 * 1e-4)\n", + " return out\n", + "\n", + "\n", + "{element: f\"{d:.1e}\" for element, (_, _, d) in diffusion(frames_600k).items()}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "63", + "metadata": {}, + "outputs": [], + "source": [ + "# diffusion coefficient in cm^2/s of each element in the NVT MD in the relaxed\n", + "# 0 K cell, at 300 to 1000 K\n", + "D_TEMPERATURES = (300, 400, 500, 600, 700, 800, 900, 1000)\n", + "DIFFUSION = {\n", + " \"OMAT-0-medium\": {\n", + " \"Li\": (-4.4e-07, 2e-06, -4.9e-08, 7.6e-06, 8e-06, 5e-06, 3.1e-05, 4.7e-05),\n", + " \"P\": (-1.3e-08, -3.3e-08, 1.2e-08, 1.4e-07, 1.3e-07, 6.9e-08, 9.6e-07, 2.6e-07),\n", + " \"S\": (-3.3e-08, -1.9e-08, 2.1e-08, 2.4e-07, 2.1e-07, 4.7e-07, 1.9e-06, 1.3e-06),\n", + " },\n", + " \"MATPES-PBE-0\": {\n", + " \"Li\": (1.9e-06, 3e-06, 6.3e-06, 1.5e-05, 2.6e-05, 2.9e-05, 4.6e-05, 6.4e-05),\n", + " \"P\": (4.6e-08, 2e-08, -1.9e-08, 1.4e-07, -2.6e-07, 2.9e-07, 1.2e-07, 3.5e-07),\n", + " \"S\": (1.2e-07, 1.8e-08, 3.9e-08, 2.7e-07, -3.7e-07, 5.4e-07, -2.7e-07, 8.1e-07),\n", + " },\n", + " \"MATPES-r2SCAN-0\": {\n", + " \"Li\": (1.6e-06, 7.9e-07, 2.6e-06, 4.8e-06, 6.8e-06, 2.6e-05, 3.8e-05, 3.4e-05),\n", + " \"P\": (\n", + " 3.2e-09,\n", + " -4.1e-09,\n", + " -6.6e-08,\n", + " 8.4e-08,\n", + " 1.5e-07,\n", + " 1.4e-07,\n", + " 4.2e-07,\n", + " -2.2e-07,\n", + " ),\n", + " \"S\": (1e-07, 9.1e-09, 7e-08, 1.8e-07, 7e-08, 3e-07, 9.3e-07, -3.7e-09),\n", + " },\n", + "}\n", + "\n", + "fig, axes = plt.subplots(1, 3, figsize=(14, 3.6), sharey=True)\n", + "inverse_t = 1000 / np.array(D_TEMPERATURES)\n", + "for ax, model in zip(axes, MODEL_FILES, strict=True):\n", + " ax.axvspan(1000 / 748, 1000 / 523, color=\"0.9\")\n", + " for element, color in ((\"Li\", \"C0\"), (\"P\", \"C2\"), (\"S\", \"C3\")):\n", + " d = np.array(DIFFUSION[model][element])\n", + " keep = d > 1e-8\n", + " ax.semilogy(inverse_t[keep], d[keep], \"-o\", color=color, ms=3, label=element)\n", + " ax.axhline(3e-7, color=\"gray\", ls=\":\")\n", + " ax.set_title(model)\n", + " ax.set_xlabel(\"1000 / T (1/K)\")\n", + "axes[0].set_ylabel(\"D (cm$^2$/s)\")\n", + "axes[0].legend()\n", + "fig.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "64", + "metadata": {}, + "source": [ + "![Diffusion coefficients of Li3PS4](finite_temperature_phonons_figures/s5_diffusion.png)\n", + "\n", + "D of Li in 10⁻⁶ cm²/s, in the relaxed 0 K cell:\n", + "\n", + "| Potential | 300 K | 400 K | 500 K | 600 K | 700 K | 800 K | 900 K | 1000 K |\n", + "|---|---|---|---|---|---|---|---|---|\n", + "| OMAT-0-medium | -0.4 | 2.0 | 0.0 | 7.6 | 8.0 | 5.0 | 31.0 | 47.0 |\n", + "| MATPES-PBE-0 | 1.9 | 3.0 | 6.3 | 15.0 | 26.0 | 29.0 | 46.0 | 64.0 |\n", + "| MATPES-r2SCAN-0 | 1.6 | 0.8 | 2.6 | 4.8 | 6.8 | 26.0 | 38.0 | 34.0 |\n", + "\n", + "A $D$ below 3×10⁻⁷ cm²/s is noise in these 7 ps runs.\n", + "\n", + "- Li diffuses from about 600 K, and from about 500 K with MACE-MATPES-PBE-0.\n", + "- P and S start to diffuse at 900 to 1000 K.\n", + "- $D$ depends on the Langevin friction. For transport, run an NVE MD." + ] + }, + { + "cell_type": "markdown", + "id": "65", + "metadata": {}, + "source": [ + "## Vibrational DOS\n", + "\n", + "The DOS is the mass-weighted power spectrum of the velocities,\n", + "\n", + "$$g(\\nu) \\propto \\sum_i m_i \\left| \\int v_i(t)\\, e^{2\\pi i \\nu t}\\, dt \\right|^2$$\n", + "\n", + "It includes all the anharmonicity and has no imaginary frequencies. The\n", + "Langevin friction $\\gamma$ broadens every peak by $\\gamma/2\\pi$ = 1.6 THz." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "66", + "metadata": {}, + "outputs": [], + "source": [ + "def md_dos(\n", + " frames: list[Atoms], smearing: float = 0.15, max_frequency: float = 25.0\n", + ") -> tuple[np.ndarray, dict[str, np.ndarray]]:\n", + " \"\"\"Get the DOS of each element from the velocities, in states/THz per Li3PS4.\n", + "\n", + " The DOS is smeared with a Gaussian of width smearing in THz, mirrored at 0.\n", + " It integrates to 3 states per atom.\n", + " \"\"\"\n", + " velocities = np.array([atoms.get_velocities() for atoms in frames])\n", + " masses = frames[0].get_masses()\n", + " symbols = np.array(frames[0].get_chemical_symbols())\n", + " nu = np.fft.rfftfreq(len(velocities), TIME_STEP)\n", + " keep = nu <= max_frequency\n", + " power = (np.abs(np.fft.rfft(velocities, axis=0)) ** 2).sum(2)[keep] * masses\n", + " step = nu[1]\n", + " half = int(4 * smearing / step)\n", + " kernel = np.exp(-0.5 * (np.arange(-half, half + 1) * step / smearing) ** 2)\n", + " n_formula_units = len(symbols) / 8\n", + " out = {}\n", + " for element in (\"Li\", \"P\", \"S\"):\n", + " g = power[:, symbols == element].sum(1)\n", + " g = np.convolve(np.concatenate([g[half:0:-1], g]), kernel, \"same\")[half:]\n", + " n_states = 3 * (symbols == element).sum() / n_formula_units\n", + " out[element] = g * n_states / (g.sum() * step)\n", + " return nu[keep], out\n", + "\n", + "\n", + "nu_600k, dos_600k = md_dos(frames_600k)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "67", + "metadata": {}, + "outputs": [], + "source": [ + "# vibrational DOS of the NVT MD in the relaxed 0 K cell, against the DOS of\n", + "# the force constants fitted to the same MD, at 300, 500 and 700 K: the fraction\n", + "# of the Li DOS below 5 THz, and the frequency in THz of the P-S stretch, the\n", + "# highest peak of the P DOS\n", + "VDOS_TEMPERATURES = (300, 500, 700)\n", + "LI_BELOW_5_THZ = {\n", + " \"OMAT-0-medium\": {\"MD\": (0.15, 0.18, 0.23), \"fit\": (0.18, 0.27, 0.9)},\n", + " \"MATPES-PBE-0\": {\"MD\": (0.21, 0.23, 0.25), \"fit\": (0.53, 0.98, 1.0)},\n", + " \"MATPES-r2SCAN-0\": {\"MD\": (0.14, 0.17, 0.21), \"fit\": (0.19, 0.41, 0.99)},\n", + "}\n", + "PS_STRETCH = {\n", + " \"OMAT-0-medium\": {\"MD\": (16.4, 16.0, 16.4), \"fit\": (15.9, 15.2, 12.9)},\n", + " \"MATPES-PBE-0\": {\"MD\": (14.6, 15.0, 15.0), \"fit\": (13.8, 12.4, 12.0)},\n", + " \"MATPES-r2SCAN-0\": {\"MD\": (16.3, 16.3, 16.3), \"fit\": (16.1, 15.2, 13.2)},\n", + "}\n", + "\n", + "fig, axes = plt.subplots(1, 2, figsize=(12, 3.6))\n", + "for idx, model in enumerate(MODEL_FILES):\n", + " for source, style in ((\"MD\", \"-o\"), (\"fit\", \"--s\")):\n", + " kwargs = {\"color\": f\"C{idx}\", \"ms\": 3, \"label\": f\"{model}, {source}\"}\n", + " axes[0].plot(VDOS_TEMPERATURES, LI_BELOW_5_THZ[model][source], style, **kwargs)\n", + " axes[1].plot(VDOS_TEMPERATURES, PS_STRETCH[model][source], style, **kwargs)\n", + "axes[0].set_ylabel(\"Fraction of the Li DOS below 5 THz\")\n", + "axes[1].set_ylabel(\"P-S stretch (THz)\")\n", + "for ax in axes:\n", + " ax.set_xlabel(\"Temperature (K)\")\n", + "axes[0].legend(fontsize=7)\n", + "fig.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "68", + "metadata": {}, + "source": [ + "![Li DOS below 5 THz and P-S stretch, MD against fit](finite_temperature_phonons_figures/s5_vdos.png)\n", + "\n", + "Li DOS below 5 THz and P-S stretch, in the relaxed 0 K cell:\n", + "\n", + "| Potential | T (K) | Li DOS below 5 THz, MD | fit | P-S stretch (THz), MD | fit |\n", + "|---|---|---|---|---|---|\n", + "| OMAT-0-medium | 300 | 15% | 18% | 16.4 | 15.9 |\n", + "| OMAT-0-medium | 500 | 18% | 27% | 16.0 | 15.2 |\n", + "| OMAT-0-medium | 700 | 23% | 90% | 16.4 | 12.9 |\n", + "| MATPES-PBE-0 | 300 | 21% | 53% | 14.6 | 13.8 |\n", + "| MATPES-PBE-0 | 500 | 23% | 98% | 15.0 | 12.4 |\n", + "| MATPES-PBE-0 | 700 | 25% | 100% | 15.0 | 12.0 |\n", + "| MATPES-r2SCAN-0 | 300 | 14% | 19% | 16.3 | 16.1 |\n", + "| MATPES-r2SCAN-0 | 500 | 17% | 41% | 16.3 | 15.2 |\n", + "| MATPES-r2SCAN-0 | 700 | 21% | 99% | 16.3 | 13.2 |\n", + "\n", + "- At 300 K the MD and the fit agree, except with MACE-MATPES-PBE-0.\n", + "- At 700 K the fit puts 90 to 100% of the Li DOS below 5 THz. The MD puts 21 to\n", + " 25% there.\n", + "- The P-S stretch stays at 15 to 16.4 THz in the MD, but drops to 12 to 13 THz in\n", + " the fit. The framework does not soften. The fit around fixed sites makes it\n", + " look soft." + ] + }, + { + "cell_type": "markdown", + "id": "69", + "metadata": {}, + "source": [ + "## Free energy from the MD\n", + "\n", + "1. **F at 100 K by Frenkel-Ladd.** We switch from an Einstein crystal, with\n", + " springs at the mean positions of the 100 K MD, to the MACE potential at 8\n", + " values of $\\lambda$. No Li hops at 100 K.\n", + " $$F = F_E + \\int_0^1 \\langle U_\\mathrm{MLIP} - U_E \\rangle_\\lambda\\, d\\lambda$$\n", + "2. **F at higher temperatures by Gibbs-Helmholtz,** with $U$ the mean energy of\n", + " the MD, or the mean enthalpy for NPT. We stop at 800 K, where P and S start to\n", + " diffuse.\n", + " $$\\frac{d(F/T)}{dT} = -\\frac{U}{T^2}$$\n", + "3. **A quantum correction** from the MD DOS, since the MD is classical." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "70", + "metadata": {}, + "outputs": [], + "source": [ + "from ase import units\n", + "from ase.calculators.calculator import Calculator, all_changes\n", + "from ase.constraints import FixCom\n", + "from ase.md.langevin import Langevin\n", + "from ase.md.velocitydistribution import MaxwellBoltzmannDistribution, Stationary\n", + "from scipy import constants\n", + "\n", + "\n", + "class Switched(Calculator):\n", + " \"\"\"Energy and forces of lam U_MLIP + (1 - lam) U_E.\n", + "\n", + " delta holds U_MLIP - U_E of the last call.\n", + " \"\"\"\n", + "\n", + " implemented_properties = (\"energy\", \"forces\")\n", + "\n", + " def __init__(\n", + " self,\n", + " mlip: Calculator,\n", + " positions: np.ndarray,\n", + " springs: np.ndarray,\n", + " lam: float,\n", + " ) -> None:\n", + " super().__init__()\n", + " self.mlip, self.r0, self.k, self.lam = mlip, positions, springs, lam\n", + "\n", + " def calculate(\n", + " self,\n", + " atoms: Atoms | None = None,\n", + " properties: tuple[str, ...] = (\"energy\",),\n", + " system_changes: list[str] = all_changes,\n", + " ) -> None:\n", + " \"\"\"Compute the switched energy and forces.\"\"\"\n", + " super().calculate(atoms, properties, system_changes)\n", + " self.mlip.calculate(self.atoms, [\"energy\", \"forces\"], system_changes)\n", + " d = self.atoms.positions - self.r0\n", + " e_einstein = 0.5 * (self.k * (d**2).sum(1)).sum()\n", + " self.delta = self.mlip.results[\"energy\"] - e_einstein\n", + " self.results = {\n", + " \"energy\": self.lam * self.mlip.results[\"energy\"]\n", + " + (1 - self.lam) * e_einstein,\n", + " \"forces\": self.lam * self.mlip.results[\"forces\"]\n", + " - (1 - self.lam) * self.k[:, None] * d,\n", + " }\n", + "\n", + "\n", + "def einstein_free_energy(\n", + " masses: np.ndarray, springs: np.ndarray, volume: float, temperature: float\n", + ") -> float:\n", + " \"\"\"Get the classical F (eV) of the Einstein crystal with a fixed center of mass.\n", + "\n", + " It includes the term for a center of mass that is free in the volume of the\n", + " cell.\n", + " \"\"\"\n", + " kt = units.kB * temperature\n", + " hbar = constants.hbar / constants.e * units.second # eV in ASE time units\n", + " free = 3 * kt * np.log(hbar * np.sqrt(springs / masses) / kt).sum()\n", + " variance = ((masses / masses.sum()) ** 2 * kt / springs).sum()\n", + " return free + kt * np.log((2 * np.pi * variance) ** 1.5 / volume)\n", + "\n", + "\n", + "def frenkel_ladd(\n", + " frames: list[Atoms],\n", + " mlip: Calculator,\n", + " temperature: float = 100.0,\n", + " n_lambda: int = 8,\n", + " steps: tuple[int, int] = (2000, 8000),\n", + " seed: int = 103,\n", + ") -> float:\n", + " \"\"\"Get F (eV per atom) of the potential at the temperature of the frames.\n", + "\n", + " steps gives the equilibration and production steps of the MD at each\n", + " lambda.\n", + " \"\"\"\n", + " rng = np.random.default_rng(seed)\n", + " atoms = frames[0].copy()\n", + " masses = atoms.get_masses()\n", + " pos = unwrapped_positions(frames)\n", + " com = (pos * masses[None, :, None]).sum(1) / masses.sum()\n", + " pos = pos - com[:, None] + com[0]\n", + " atoms.positions = pos.mean(0)\n", + " msd = ((pos - atoms.positions) ** 2).sum(2).mean(0)\n", + " symbols = np.array(atoms.get_chemical_symbols())\n", + " springs = np.empty(len(atoms))\n", + " for element in set(symbols):\n", + " springs[symbols == element] = (\n", + " 3 * units.kB * temperature / msd[symbols == element].mean()\n", + " )\n", + " nodes, weights = np.polynomial.legendre.leggauss(n_lambda)\n", + " integral = 0.0\n", + " for node, weight in zip(nodes, weights, strict=True):\n", + " run = atoms.copy()\n", + " run.set_constraint(FixCom())\n", + " run.calc = Switched(mlip, atoms.positions.copy(), springs, (node + 1) / 2)\n", + " MaxwellBoltzmannDistribution(run, temperature_K=temperature, rng=rng)\n", + " Stationary(run)\n", + " dyn = Langevin(\n", + " run,\n", + " units.fs,\n", + " temperature_K=temperature,\n", + " friction=0.01 / units.fs,\n", + " fixcm=False,\n", + " rng=rng,\n", + " )\n", + " dyn.run(steps[0])\n", + " deltas = []\n", + " for _ in range(steps[1] // 10):\n", + " dyn.run(10)\n", + " deltas.append(run.calc.delta)\n", + " integral += weight / 2 * np.mean(deltas)\n", + " f_e = einstein_free_energy(masses, springs, atoms.get_volume(), temperature)\n", + " return (f_e + integral) / len(atoms)" + ] + }, + { + "cell_type": "markdown", + "id": "71", + "metadata": {}, + "source": [ + "On a Cu-Ag crystal with EMT at 10 K, these functions gave 0.247542 ± 0.000087 eV\n", + "against the exact 0.247599 eV. The next cell runs them on the 100 K run of\n", + "Section 4. It is off by default." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "72", + "metadata": {}, + "outputs": [], + "source": [ + "RUN_FL = False\n", + "if RUN_FL:\n", + " mlip = mace_mp(**calc_kwargs(\"OMAT-0-medium\"))\n", + " frames = nvt_trajectory(\"runs/Li3PS4_OMAT-0-medium_100K_nvt\")\n", + " f_100 = frenkel_ladd(frames, mlip) # eV per atom" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "73", + "metadata": {}, + "outputs": [], + "source": [ + "PLANCK = constants.h / constants.e * 1e12 # eV/THz\n", + "\n", + "\n", + "def gibbs_helmholtz(\n", + " t0: float, f0: float, temperatures: list[float], energies: list[float]\n", + ") -> np.ndarray:\n", + " \"\"\"Get F at each temperature from F at t0 and the mean energy U.\n", + "\n", + " F and U are in eV per atom. temperatures starts at t0, and U is linear in T\n", + " between the temperatures.\n", + " \"\"\"\n", + " f_over_t = [f0 / t0]\n", + " for i in range(1, len(temperatures)):\n", + " t1, t2 = temperatures[i - 1], temperatures[i]\n", + " b = (energies[i] - energies[i - 1]) / (t2 - t1)\n", + " a = energies[i - 1] - b * t1\n", + " f_over_t.append(f_over_t[-1] - (a * (1 / t1 - 1 / t2) + b * np.log(t2 / t1)))\n", + " return np.array(f_over_t) * np.array(temperatures)\n", + "\n", + "\n", + "def quantum_correction(\n", + " nu: np.ndarray, dos: np.ndarray, temperature: float\n", + ") -> tuple[float, float]:\n", + " \"\"\"Get the quantum minus the classical harmonic F (eV) and S (eV/K) of a DOS.\n", + "\n", + " nu is in THz. The result is per formula unit for a DOS per formula unit, as\n", + " md_dos gives it. Divide by 8 for a value per atom.\n", + " \"\"\"\n", + " keep = nu > 0\n", + " nu, g, step = nu[keep], dos[keep], nu[1] - nu[0]\n", + " kt = units.kB * temperature\n", + " x = PLANCK * nu / kt\n", + " df = PLANCK * nu / 2 + kt * (np.log1p(-np.exp(-x)) - np.log(x))\n", + " ds = units.kB * (x / np.expm1(x) - np.log1p(-np.exp(-x)) - 1 + np.log(x))\n", + " return (df * g).sum() * step, (ds * g).sum() * step" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "74", + "metadata": {}, + "outputs": [], + "source": [ + "# per formula unit, at 100 to 800 K: F in kJ/mol relative to the relaxed 0 K\n", + "# energy and S in J/(K mol). From the MD in the relaxed 0 K cell (nvt) and in the\n", + "# NPT cell (npt), from the force constants fitted in the 0 K cell with U0 added,\n", + "# and from the 0 K force constants.\n", + "THERMO_T = (100, 200, 300, 400, 500, 600, 700, 800)\n", + "MD_F = {\n", + " \"nvt\": {\n", + " \"OMAT-0-medium\": (31.5, 21.2, 2.7, -22.9, -52.3, -87.4, -125.5, -167.4),\n", + " \"MATPES-PBE-0\": (27.2, 13.7, -7.6, -34.6, -67.0, -103.4, -143.5, -186.4),\n", + " \"MATPES-r2SCAN-0\": (31.2, 19.2, 0.3, -25.0, -55.6, -90.5, -129.0, -170.8),\n", + " },\n", + " \"npt\": {\n", + " \"OMAT-0-medium\": (31.1, 19.6, -0.4, -26.5, -58.7, -94.0, -134.6, -178.8),\n", + " \"MATPES-PBE-0\": (27.6, 13.6, -9.0, -36.5, -69.8, -106.9, -148.2, -192.3),\n", + " \"MATPES-r2SCAN-0\": (30.9, 18.5, -1.4, -27.5, -58.6, -94.8, -134.6, -178.4),\n", + " },\n", + "}\n", + "MD_S = {\n", + " \"nvt\": {\n", + " \"OMAT-0-medium\": (63, 150, 218, 273, 317, 363, 401, 435),\n", + " \"MATPES-PBE-0\": (79, 169, 242, 296, 343, 384, 415, 441),\n", + " \"MATPES-r2SCAN-0\": (68, 154, 221, 278, 326, 361, 400, 433),\n", + " },\n", + " \"npt\": {\n", + " \"OMAT-0-medium\": (68, 158, 230, 286, 334, 380, 419, 461),\n", + " \"MATPES-PBE-0\": (87, 178, 246, 306, 351, 392, 426, 456),\n", + " \"MATPES-r2SCAN-0\": (76, 159, 229, 284, 334, 379, 417, 449),\n", + " },\n", + "}\n", + "FIT_F = {\n", + " \"OMAT-0-medium\": (32.2, 19.7, 0.3, -36.5, -60.3, -118.8, -166.6, -219.7),\n", + " \"MATPES-PBE-0\": (27.1, 12.4, -17.5, -49.0, -95.1, -139.5, -184.4, -240.1),\n", + " \"MATPES-r2SCAN-0\": (31.3, 18.9, -0.7, -29.4, -68.5, -115.1, -168.9, -225.2),\n", + "}\n", + "FIT_S = {\n", + " \"OMAT-0-medium\": (68, 155, 223, 305, 333, 410, 452, 493),\n", + " \"MATPES-PBE-0\": (82, 179, 273, 326, 391, 433, 457, 495),\n", + " \"MATPES-r2SCAN-0\": (65, 155, 224, 288, 350, 401, 449, 489),\n", + "}\n", + "HARMONIC_F = {\n", + " \"OMAT-0-medium\": (35.4, 25.5, 8.8, -13.4, -40.1, -70.3, -103.7, -139.7),\n", + " \"MATPES-PBE-0\": (30.0, 18.0, -1.4, -26.7, -56.5, -90.1, -126.8, -166.2),\n", + " \"MATPES-r2SCAN-0\": (35.6, 26.0, 9.6, -12.3, -38.4, -68.2, -101.0, -136.5),\n", + "}\n", + "HARMONIC_S = {\n", + " \"OMAT-0-medium\": (59, 135, 197, 246, 286, 319, 348, 373),\n", + " \"MATPES-PBE-0\": (75, 160, 226, 277, 318, 352, 381, 407),\n", + " \"MATPES-r2SCAN-0\": (58, 132, 193, 241, 281, 314, 342, 367),\n", + "}\n", + "\n", + "fig, axes = plt.subplots(2, 3, figsize=(14, 7), sharex=True)\n", + "for col, model in enumerate(MODEL_FILES):\n", + " for row, (md, fit, harmonic) in enumerate(\n", + " ((MD_F, FIT_F, HARMONIC_F), (MD_S, FIT_S, HARMONIC_S))\n", + " ):\n", + " ax = axes[row, col]\n", + " ax.plot(THERMO_T, md[\"nvt\"][model], \"-o\", ms=3, label=\"MD, 0 K cell\")\n", + " ax.plot(THERMO_T, md[\"npt\"][model], \"-s\", ms=3, label=\"MD, NPT cell\")\n", + " ax.plot(THERMO_T, fit[model], \"--\", label=\"fit\")\n", + " ax.plot(\n", + " THERMO_T, harmonic[model], \":\", color=\"gray\", label=\"0 K force constants\"\n", + " )\n", + " axes[0, col].set_title(model)\n", + " axes[1, col].set_xlabel(\"Temperature (K)\")\n", + "axes[0, 0].set_ylabel(\"F (kJ/mol)\")\n", + "axes[1, 0].set_ylabel(\"S (J/(K mol))\")\n", + "axes[1, 0].legend(fontsize=8)\n", + "fig.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "75", + "metadata": {}, + "source": [ + "![Free energy and entropy of Li3PS4](finite_temperature_phonons_figures/s5_thermo.png)\n", + "\n", + "S at 700 K in J/(K mol) per formula unit:\n", + "\n", + "| Potential | MD, 0 K cell | MD, NPT cell | fit | 0 K force constants |\n", + "|---|---|---|---|---|\n", + "| OMAT-0-medium | 401 | 419 | 452 | 348 |\n", + "| MATPES-PBE-0 | 415 | 426 | 457 | 381 |\n", + "| MATPES-r2SCAN-0 | 400 | 417 | 449 | 342 |\n", + "\n", + "The fitted F includes U₀, the mean potential energy of its MD minus 3kT/2 per\n", + "atom.\n", + "\n", + "- Up to 300 K the fit and the MD agree within 1 to 2.5 kJ/mol, except with\n", + " MACE-MATPES-PBE-0.\n", + "- Where Li hops, the fit overestimates S at 700 K by 40 to 50 J/(K mol).\n", + "- The 0 K force constants underestimate it by 35 to 60 J/(K mol). So does the\n", + " entropy in the Materials Project's Harmonic Phonon Database." + ] + }, + { + "cell_type": "markdown", + "id": "76", + "metadata": {}, + "source": [ + "## γ against β\n", + "\n", + "The phase with the lowest G is stable, so a crossing of G(β) and G(γ) is a\n", + "transition. The Materials Project has no γ entry and no α framework. We take γ\n", + "from Kaup et al. (Table S1), which matches Homma et al. (COD 1570309)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "77", + "metadata": {}, + "outputs": [], + "source": [ + "from pymatgen.core import Lattice\n", + "\n", + "gamma = Structure.from_spacegroup(\n", + " \"Pmn2_1\",\n", + " Lattice.orthorhombic(7.7557, 6.5627, 6.1362),\n", + " [\"Li\", \"Li\", \"P\", \"S\", \"S\", \"S\"],\n", + " [\n", + " [0.2441, 0.3122, -0.0029],\n", + " [0, 0.1493, 0.475],\n", + " [0, 0.8174, 0.9955],\n", + " [0.2177, 0.6717, 0.8864],\n", + " [0, 0.1115, 0.8935],\n", + " [0, 0.8094, 0.3274],\n", + " ],\n", + ")\n", + "GAMMA_SUPERCELL = (2, 2, 2) # 128 atoms\n", + "gamma.composition, gamma.get_space_group_info()" + ] + }, + { + "cell_type": "markdown", + "id": "78", + "metadata": {}, + "source": [ + "At 0 K β lies above γ by 4.7 meV/atom with MACE-OMAT-0-medium, 11.9 with\n", + "MACE-MATPES-PBE-0 and 10.5 with MACE-MATPES-r2SCAN-0. We run γ like β, with the\n", + "NPT MD and Frenkel-Ladd. The loop is off by default." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "79", + "metadata": {}, + "outputs": [], + "source": [ + "RUN_GAMMA = False\n", + "gamma_docs = {}\n", + "if RUN_GAMMA:\n", + " for model in MODEL_FILES:\n", + " for temperature in LI3PS4_TEMPERATURES:\n", + " gamma_docs[model, temperature] = run_ft(\n", + " gamma,\n", + " GAMMA_SUPERCELL,\n", + " model,\n", + " f\"runs/gamma_{model}_{temperature}K_npt\",\n", + " temperature,\n", + " relax=True,\n", + " run_npt=True,\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "80", + "metadata": {}, + "outputs": [], + "source": [ + "# G(beta) - G(gamma) in meV/atom and S of each phase in J/(K mol) per formula\n", + "# unit, at 100 to 800 K and zero pressure, with the quantum correction\n", + "PHASE_T = (100, 200, 300, 400, 500, 600, 700, 800)\n", + "DELTA_G = {\n", + " \"OMAT-0-medium\": (-0.3, -1.8, -3.6, -5.1, -6.6, -7.3, -6.9, -4.9),\n", + " \"MATPES-PBE-0\": (7.1, 5.5, 3.4, 1.1, -1.1, -2.2, -1.9, -0.4),\n", + " \"MATPES-r2SCAN-0\": (3.0, 1.6, -0.5, -2.4, -4.6, -7.2, -9.5, -10.5),\n", + "}\n", + "PHASE_S = {\n", + " \"gamma\": {\n", + " \"OMAT-0-medium\": (80, 156, 222, 278, 327, 380, 430, 484),\n", + " \"MATPES-PBE-0\": (92, 172, 238, 292, 340, 391, 436, 474),\n", + " \"MATPES-r2SCAN-0\": (79, 154, 219, 273, 319, 360, 403, 450),\n", + " },\n", + " \"beta\": {\n", + " \"OMAT-0-medium\": (90, 168, 235, 290, 338, 383, 422, 463),\n", + " \"MATPES-PBE-0\": (103, 188, 255, 311, 355, 395, 429, 458),\n", + " \"MATPES-r2SCAN-0\": (91, 169, 235, 288, 337, 382, 419, 451),\n", + " },\n", + "}\n", + "\n", + "fig, axes = plt.subplots(1, 2, figsize=(12, 3.8))\n", + "for ax in axes:\n", + " ax.axvspan(523, 573, color=\"0.9\")\n", + "axes[0].axhline(0, color=\"gray\", lw=0.5)\n", + "for idx, model in enumerate(MODEL_FILES):\n", + " axes[0].plot(PHASE_T, DELTA_G[model], \"-o\", color=f\"C{idx}\", ms=3, label=model)\n", + " for phase, style in ((\"gamma\", \"--\"), (\"beta\", \"-\")):\n", + " axes[1].plot(PHASE_T, PHASE_S[phase][model], style, color=f\"C{idx}\")\n", + "axes[0].set_ylabel(\"G(beta) - G(gamma) (meV/atom)\")\n", + "axes[1].set_ylabel(\"S (J/(K mol)), beta solid, gamma dashed\")\n", + "for ax in axes:\n", + " ax.set_xlabel(\"Temperature (K)\")\n", + "axes[0].legend(fontsize=8)\n", + "fig.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "81", + "metadata": {}, + "source": [ + "![G(beta) - G(gamma) and entropies](finite_temperature_phonons_figures/s5_phases.png)\n", + "\n", + "G(β) - G(γ) in meV/atom:\n", + "\n", + "| Potential | 100 K | 200 K | 300 K | 400 K | 500 K | 600 K | 700 K | 800 K |\n", + "|---|---|---|---|---|---|---|---|---|\n", + "| OMAT-0-medium | -0.3 | -1.8 | -3.6 | -5.1 | -6.6 | -7.3 | -6.9 | -4.9 |\n", + "| MATPES-PBE-0 | 7.1 | 5.5 | 3.4 | 1.1 | -1.1 | -2.2 | -1.9 | -0.4 |\n", + "| MATPES-r2SCAN-0 | 3.0 | 1.6 | -0.5 | -2.4 | -4.6 | -7.2 | -9.5 | -10.5 |\n", + "\n", + "β becomes stable above 278 K with MACE-MATPES-r2SCAN-0 and above 450 K with\n", + "MACE-MATPES-PBE-0, against 523 to 573 K in experiment. With MACE-OMAT-0-medium β\n", + "is already lower at 100 K. The second crossing of MACE-MATPES-PBE-0 at 822 K\n", + "does not count, since γ flows there.\n", + "\n", + "- β is softer. At 100 K its classical F lies only 0.2 to 7.9 meV/atom above γ,\n", + " much less than the gap at 0 K.\n", + "- The quantum correction is about 5 meV/atom for each phase at 500 K, but it\n", + " changes G(β) - G(γ) by 1.1 meV/atom or less.\n", + "- The error of each Frenkel-Ladd F is about 0.1 meV/atom." + ] + }, + { + "cell_type": "markdown", + "id": "82", + "metadata": {}, + "source": [ + "## What this approach cannot do\n", + "\n", + "- **α.** Its ordered versions relax into distorted structures 18 to 22 meV/atom\n", + " above γ, so Frenkel-Ladd has no stable start. α needs a treatment of its Li\n", + " disorder, for example a cluster expansion.\n", + "- **Melting.** It needs a reference for the liquid.\n", + "- **Kinetics.** The free energies give the equilibrium temperature only, not the\n", + " hysteresis seen in experiment.\n", + "\n", + "## Why the γ to β temperature is off\n", + "\n", + "- **The potentials.** 1 meV/atom in G(β) - G(γ) moves the crossing by about\n", + " 50 K, and the 0 K gaps differ by 7 meV/atom between the potentials. Kam et al.\n", + " found r²SCAN DFT 200 to 300 K too low (arXiv:2307.00878), and\n", + " MACE-MATPES-r2SCAN-0 falls in the same range.\n", + "- **One ordering of β.** Below about 500 K the MD misses part of the Li disorder\n", + " of β.\n", + "- **Size and length.** 128 atoms and 8 ps of MD." + ] + }, + { + "cell_type": "markdown", + "id": "83", + "metadata": {}, + "source": [ + "## Known limitations\n", + "\n", + "- Without `run_npt=True` the MD keeps the volume of the starting structure. The\n", + " NPT cell comes from one MD run and carries its noise.\n", + "- MACE gives no Born charges, so no non-analytical correction is applied. A\n", + " `born_maker` adds it, either `ForceFieldDielectricMaker` with MACE-Field or a\n", + " VASP `DielectricMaker`.\n", + "- The trajectory check uses fixed limits. The Lindemann limit of 0.15 is that of\n", + " fcc solids, and light atoms that hop can pass it while the crystal stays solid.\n", + "- Each result comes from one trajectory. Runs differ by up to 0.1 THz on average\n", + " and 0.3 THz in the highest frequency (Section 3)." + ] + } + ], + "metadata": { + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/tutorials/finite_temperature_phonons_figures/s1_bands_AlAsPt5.png b/tutorials/finite_temperature_phonons_figures/s1_bands_AlAsPt5.png new file mode 100644 index 0000000000..c8ad751da7 Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s1_bands_AlAsPt5.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s1_bands_Bi4S3N2.png b/tutorials/finite_temperature_phonons_figures/s1_bands_Bi4S3N2.png new file mode 100644 index 0000000000..d176a4dec5 Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s1_bands_Bi4S3N2.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s1_bands_Ca3Ir4Sn13.png b/tutorials/finite_temperature_phonons_figures/s1_bands_Ca3Ir4Sn13.png new file mode 100644 index 0000000000..9ff7bd7362 Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s1_bands_Ca3Ir4Sn13.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s1_bands_Er5Tl3.png b/tutorials/finite_temperature_phonons_figures/s1_bands_Er5Tl3.png new file mode 100644 index 0000000000..99b519990a Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s1_bands_Er5Tl3.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s1_bands_Zn3P2.png b/tutorials/finite_temperature_phonons_figures/s1_bands_Zn3P2.png new file mode 100644 index 0000000000..681f1488c5 Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s1_bands_Zn3P2.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s1_ca3ir4sn13_300k.png b/tutorials/finite_temperature_phonons_figures/s1_ca3ir4sn13_300k.png new file mode 100644 index 0000000000..10af076ff0 Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s1_ca3ir4sn13_300k.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s1_max_frequency.png b/tutorials/finite_temperature_phonons_figures/s1_max_frequency.png new file mode 100644 index 0000000000..a5aaa72c4c Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s1_max_frequency.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s2_bands_AlAsPt5.png b/tutorials/finite_temperature_phonons_figures/s2_bands_AlAsPt5.png new file mode 100644 index 0000000000..52408670de Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s2_bands_AlAsPt5.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s2_bands_Bi4S3N2.png b/tutorials/finite_temperature_phonons_figures/s2_bands_Bi4S3N2.png new file mode 100644 index 0000000000..36a15ceba3 Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s2_bands_Bi4S3N2.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s2_bands_Ca3Ir4Sn13.png b/tutorials/finite_temperature_phonons_figures/s2_bands_Ca3Ir4Sn13.png new file mode 100644 index 0000000000..9d7d822cd0 Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s2_bands_Ca3Ir4Sn13.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s2_bands_Er5Tl3.png b/tutorials/finite_temperature_phonons_figures/s2_bands_Er5Tl3.png new file mode 100644 index 0000000000..69bc97bb05 Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s2_bands_Er5Tl3.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s2_bands_Zn3P2.png b/tutorials/finite_temperature_phonons_figures/s2_bands_Zn3P2.png new file mode 100644 index 0000000000..990147fb62 Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s2_bands_Zn3P2.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s2_max_frequency.png b/tutorials/finite_temperature_phonons_figures/s2_max_frequency.png new file mode 100644 index 0000000000..2dec04186e Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s2_max_frequency.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s2_volume.png b/tutorials/finite_temperature_phonons_figures/s2_volume.png new file mode 100644 index 0000000000..98dda08837 Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s2_volume.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s3_knaicl_grid.png b/tutorials/finite_temperature_phonons_figures/s3_knaicl_grid.png new file mode 100644 index 0000000000..eeb92b02f6 Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s3_knaicl_grid.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s4_bands_MATPES-PBE-0.png b/tutorials/finite_temperature_phonons_figures/s4_bands_MATPES-PBE-0.png new file mode 100644 index 0000000000..e1890c21a5 Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s4_bands_MATPES-PBE-0.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s4_bands_MATPES-r2SCAN-0.png b/tutorials/finite_temperature_phonons_figures/s4_bands_MATPES-r2SCAN-0.png new file mode 100644 index 0000000000..19611ce979 Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s4_bands_MATPES-r2SCAN-0.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s4_bands_OMAT-0-medium.png b/tutorials/finite_temperature_phonons_figures/s4_bands_OMAT-0-medium.png new file mode 100644 index 0000000000..a74f4c7eff Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s4_bands_OMAT-0-medium.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s4_runs.png b/tutorials/finite_temperature_phonons_figures/s4_runs.png new file mode 100644 index 0000000000..3a720f5192 Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s4_runs.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s5_diffusion.png b/tutorials/finite_temperature_phonons_figures/s5_diffusion.png new file mode 100644 index 0000000000..5dc96b9d9f Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s5_diffusion.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s5_phases.png b/tutorials/finite_temperature_phonons_figures/s5_phases.png new file mode 100644 index 0000000000..311cb27e49 Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s5_phases.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s5_thermo.png b/tutorials/finite_temperature_phonons_figures/s5_thermo.png new file mode 100644 index 0000000000..6da0ddf7ff Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s5_thermo.png differ diff --git a/tutorials/finite_temperature_phonons_figures/s5_vdos.png b/tutorials/finite_temperature_phonons_figures/s5_vdos.png new file mode 100644 index 0000000000..81bdde0fd6 Binary files /dev/null and b/tutorials/finite_temperature_phonons_figures/s5_vdos.png differ diff --git a/tutorials/tutorials.md b/tutorials/tutorials.md index 487c6da5e9..c530f277bc 100644 --- a/tutorials/tutorials.md +++ b/tutorials/tutorials.md @@ -19,6 +19,7 @@ hiphive_workflow force_fields/phonon_workflow grueneisen_workflow cte_workflow +finite_temperature_phonons qha_workflow torchsim_tutorial ```