From 471a07b07567d6ad88597fbe8b4354798c681c65 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Thu, 24 Sep 2026 01:05:09 -0700 Subject: [PATCH 01/28] Fix pheasy anharmonic fitting and add fit options --- docs/user/codes/vasp.md | 18 +- pyproject.toml | 5 +- src/atomate2/common/flows/hiphive.py | 12 +- src/atomate2/common/flows/pheasy.py | 141 ++++-- src/atomate2/common/jobs/hiphive.py | 76 +-- src/atomate2/common/jobs/pheasy.py | 665 +++++++++++++++++---------- src/atomate2/common/jobs/phonons.py | 111 +++++ src/atomate2/vasp/flows/pheasy.py | 93 ++-- tests/common/jobs/test_pheasy.py | 304 ++++++++++++ tests/vasp/flows/test_pheasy.py | 65 +++ tutorials/pheasy_workflow.ipynb | 10 +- 11 files changed, 1131 insertions(+), 369 deletions(-) create mode 100644 tests/common/jobs/test_pheasy.py diff --git a/docs/user/codes/vasp.md b/docs/user/codes/vasp.md index edd4f46cc3..ee9e9ba3a7 100644 --- a/docs/user/codes/vasp.md +++ b/docs/user/codes/vasp.md @@ -345,12 +345,14 @@ phonon_flow = PhononMaker(min_length=15.0, store_force_constants=False).make( #### Pheasy 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). -`Pheasy` can be installed with `pip install pheasy`. +The `pheasy` extra, `pip install "atomate2[pheasy]"`, installs the pheasy version this workflow needs, together with phonopy and ALM. +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`. This workflow was used to build the Materials Project's Harmonic Phonon Database, described in [this preprint](https://chemrxiv.org/doi/full/10.26434/chemrxiv.15004632/v1). -By default, this workflow does not compute anharmonic force constants, but can be extended to using the `cal_anhar_fcs` kwarg and the `ALAMODE` code. +By default, this workflow does not compute anharmonic force constants, but can be extended to using the `cal_anhar_fcs` kwarg. +ALM, from the ALAMODE package, counts the free force constants used to size the random displacement sets. To install ALAMODE, see their [installation guidelines](https://alamode.readthedocs.io/en/latest/install.html#). Linux and MacOS x86-64 users can try to install using conda forge: @@ -366,6 +368,7 @@ cd ALM/python python setup.py build pip install -e . ``` +The `pheasy` extra pins ALM to commit `f1d668f`. When building ALM by hand, check out that commit. NB: MacOS users will need to ensure that `gcc` and `g++` are used rather than `clang` - both can be installed with `homebrew`. Note also that `boost` and `eigen` can be installed via `homebrew`. For example, using `gcc-15` from `homebrew`, one might set: @@ -373,6 +376,17 @@ For example, using `gcc-15` from `homebrew`, one might set: export CC=gcc-15 ; CXX=g++-15 ; CXX_FLAGS=-DOPENMP ``` +With `cal_anhar_fcs=True`, the anharmonic force constants are fitted with LASSO to a second set of randomly displaced supercells. +The number of these supercells is set from the number of free force constants, so that the fit has 100 force equations per free force constant. +Unless `num_disp_anhar` is set, at least 20 supercells are used, and above 600 the job stops and asks for a shorter cutoff, a larger supercell or an explicit `num_disp_anhar`. An explicit `num_disp_anhar` is used as given. +`anhar_max_order` selects third-order force constants (3, the default) or third- and fourth-order force constants (4). +`anhar_fit_methods` selects how the second-order force constants are treated. +`"cocktail"` (the default) keeps them fixed to the harmonic fit, and `"one-shot"` fits them together with the higher orders and writes the results to a `one_shot` folder. +Both can be requested in one run. +If the cross-validated LASSO penalty lands on either end of the search, `10**anhar_alpha_min` or pheasy's `1e-2`, a warning is raised. +The harmonic and anharmonic LASSO fits are seeded, so that repeated runs give the same force constants. +The anharmonic force constants are written to files in the job folder and are not stored in the output document. + #### hiPhive The same force constants can instead be fitted with [hiPhive](https://hiphive.materialsmodeling.org/), which builds a cluster expansion of the force constant potential and fits it by regression. diff --git a/pyproject.toml b/pyproject.toml index db6701c909..75d5a4cbcc 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -52,8 +52,9 @@ mp = ["mp-api>=0.37.5"] # phonon.save -> phonopy.load round-trip used by the Grüneisen workflow # ("Force constants shape disagrees with crystal structure setting"). phonons = ["phonopy>=2.43,<5", "seekpath>=2.0.0"] -pheasy = ["hiphive==1.3.1", "numpy<=2.2", "pheasy @ git+https://gitlab.com/hpsahasrabuddhe/pheasy.git@support-disp-force-matrix-input"] -alamode = ["alm @ git+https://github.com/ttadano/ALM.git@develop#subdirectory=python"] +# the pheasy fork adds the --disp_matrix_file and --force_matrix_file options +pheasy = ["atomate2[phonons,alamode]", "hiphive==1.3.1", "numpy<=2.2", "pheasy @ git+https://gitlab.com/hpsahasrabuddhe/pheasy.git@e67a777ecfde5ab368d71167837b099df75cfd4c"] +alamode = ["alm @ git+https://github.com/ttadano/ALM.git@f1d668fdee66e7e7218a04c88daf19d0e14fce0c#subdirectory=python"] hiphive = ["hiphive==1.3.1", "trainstation>=1.0", "atomate2[phonons,alamode]"] lobster = ["ijson>=3.2.2", "lobsterpy>=0.6.0"] defects = [ diff --git a/src/atomate2/common/flows/hiphive.py b/src/atomate2/common/flows/hiphive.py index 7fb63172fb..37bbd9a543 100644 --- a/src/atomate2/common/flows/hiphive.py +++ b/src/atomate2/common/flows/hiphive.py @@ -1,4 +1,4 @@ -"""Flow for calculating (an)harmonic FCs and phonon renormalisation with hiPhive.""" +"""Flow for calculating harmonic FCs with hiPhive.""" from __future__ import annotations @@ -6,6 +6,8 @@ from dataclasses import dataclass, field from typing import TYPE_CHECKING, Literal +from pymatgen.util.due import Doi, due + from atomate2.common.flows.phonons import BasePhononMaker as PurePhonopyMaker from atomate2.common.jobs.hiphive import ( generate_frequencies_eigenvectors, @@ -28,6 +30,14 @@ SUPPORTED_CODES = frozenset(("vasp", "aims", "forcefields")) +@due.dcite( + Doi("10.1002/adts.201800184"), + description="hiPhive, force constant potentials by regression.", +) +@due.dcite( + Doi("10.1088/0953-8984/26/22/225402"), + description="ALM, used to count the free force constants.", +) @dataclass class BasePhononMaker(PurePhonopyMaker, ABC): """Maker to calculate harmonic phonons with the cluster-expansion code hiPhive. diff --git a/src/atomate2/common/flows/pheasy.py b/src/atomate2/common/flows/pheasy.py index 12ed008b71..158c90efab 100644 --- a/src/atomate2/common/flows/pheasy.py +++ b/src/atomate2/common/flows/pheasy.py @@ -1,4 +1,4 @@ -"""Flow for calculating (an)harmonic FCs and phonon renormalisation with pheasy.""" +"""Flow for calculating harmonic and anharmonic FCs with pheasy.""" from __future__ import annotations @@ -10,6 +10,7 @@ from atomate2.common.flows.phonons import BasePhononMaker as PurePhonopyMaker from atomate2.common.jobs.pheasy import ( + _check_anharmonic_settings, generate_frequencies_eigenvectors, generate_phonon_displacements, get_supercell_size, @@ -17,6 +18,7 @@ from atomate2.common.jobs.phonons import run_phonon_displacements if TYPE_CHECKING: + from collections.abc import Sequence from pathlib import Path from emmet.core.math import Matrix3D @@ -34,24 +36,34 @@ Doi("10.26434/chemrxiv.15004632/v1"), description="Materials Project's Harmonic Phonon Database.", ) +@due.dcite( + Doi("10.48550/arXiv.2508.01020"), + description="Pheasy code for (an)harmonic force constants.", +) +@due.dcite( + Doi("10.1088/0953-8984/26/22/225402"), + description="ALM, used to count the free force constants.", +) @dataclass class BasePhononMaker(PurePhonopyMaker, ABC): - """Maker to calculate harmonic phonons with LASSO-based ML code Pheasy. - - Calculate the zero-K harmonic phonons of a material and higher-order FCs. - Initially, a tight structural relaxation is performed to obtain a structure - without forces on the atoms. Subsequently, supercells with all atoms displaced - by a small amplitude (generally using 0.01 A) are generated and accurate forces - are computed for these structures for the second order force constants. With the - help of pheasy (LASSO technique), these forces are then converted into a dynamical - matrix. In this Workflow, we separate the harmonic phonon calculations and - anharmonic force constants calculations. To correct for polarization effects, a + """Maker to calculate harmonic phonons and anharmonic FCs with Pheasy. + + Calculate the zero-K harmonic phonons of a material and, optionally, its + anharmonic FCs. Initially, a tight structural relaxation is performed to obtain + a structure without forces on the atoms. Subsequently, displaced supercells + (generally using 0.01 A) are generated and accurate forces are computed for + them. If phonopy needs more than three finite displacements, all atoms are + displaced randomly and pheasy fits the second order force constants with LASSO. + Otherwise, the finite displacements are used with a least-squares fit. These + force constants are then converted into a dynamical matrix. In this Workflow, + we separate the harmonic phonon calculations and anharmonic force constants + calculations. To correct for polarization effects, a correction of the dynamical matrix based on BORN charges can be performed. Finally, phonon densities of states, phonon band structures and thermodynamic properties are computed. For the anharmonic force constants, the supercells with all atoms displaced by a larger amplitude (generally using 0.08 A) are generated and accurate forces are computed for these structures. With the help of pheasy (LASSO technique), - the third- and fourth-order force constants are extracted at once. + the third-order (and optionally fourth-order) force constants are extracted. .. Note:: It is heavily recommended to symmetrize the structure before passing it to @@ -59,10 +71,8 @@ class BasePhononMaker(PurePhonopyMaker, ABC): displacement calculations will be required for pheasy phonon calculation. It is recommended to check the convergence parameters here and adjust them if necessary. The default might not be strict enough for your specific case. - Additionally, for high-throughoput calculations, it is recommended to calculate - the residual forces on the atoms in the supercell after the relaxation. Then the - forces on displaced supercells can deduct the residual forces to reduce the - error in the dynamical matrix. + The residual forces of the undisplaced supercell are always subtracted from + the forces of the displaced supercells. Parameters ---------- @@ -83,35 +93,62 @@ class BasePhononMaker(PurePhonopyMaker, ABC): number of displacements to be generated using a random-displacement approach for harmonic phonon calculations. The default value is 0 and the number of displacements is automatically determined by the number of atoms in the - supercell and its space group. + supercell and its space group. Not used when phonopy needs at most three + finite displacements. cal_anhar_fcs: bool - if set to True, anharmonic force constants(FCs) up to fourth-order FCs will - be calculated. The default value is False, and only harmonic phonons will - be calculated. + if set to True, anharmonic force constants(FCs) up to order + anhar_max_order will be calculated. The default value is False, and only + harmonic phonons will be calculated. displacement_anhar: float - displacement distance for anharmonic force constants(FCs) up to fourth-order - FCs, for most cases 0.08 A is a good choice, but it can be increased to 0.1 A. + displacement distance for the anharmonic force constants(FCs). The default + is 0.08 A. num_disp_anhar: int number of displacements to be generated using a random-displacement approach - for anharmonic phonon calculations. The default value is 0 and the number of - displacements is automatically determined by the number of atoms in the - supercell, cutoff distance for anharmonic FCs its space group. generally, - 50 large-distance displacements are enough for most cases. + for anharmonic phonon calculations. A non-zero value is used as given. The + default value is 0, and then the number is set so that the fit has 100 force + equations per free force constant, with at least 20 displaced supercells. + Above 600 the job stops. The free force constants are counted with ALM, which + uses its own symmetry search, for this supercell and fcs_cutoff_radius. The + second-order FCs are included when "one-shot" is requested. + anhar_max_order: int + highest order of the anharmonic FCs, 3 or 4. The default of 3 fits the + third-order FCs. 4 also fits fourth-order FCs. + anhar_fit_methods: Sequence[Literal["cocktail", "one-shot"]] + how the anharmonic FCs are fitted to the large-distance dataset. + "cocktail" keeps the second-order FCs fixed to the harmonic fit and fits + only the higher orders. "one-shot" fits the second- and higher-order FCs + together from the same dataset and writes them to the "one_shot" folder. + Both can be requested. The default is ("cocktail",). The cocktail FCs are + written to the job folder, and the one-shot FCs, including fc2.hdf5, to its + one_shot subfolder: FORCE_CONSTANTS_3RD and fc3.hdf5, plus + FORCE_CONSTANTS_4TH and fc4.hdf5 at fourth order. They are not stored in the + output document. The LASSO fits are seeded, so a run is reproducible. With + few displaced supercells, the one-shot result depends on that seed. + anhar_alpha_min: int + base-10 exponent of the smallest LASSO penalty tried by cross-validation + in the anharmonic fits. It must be an integer below -2. The default of -12 + is below pheasy's own default of -6. If the chosen penalty lands on either + end of the search, 10**anhar_alpha_min or pheasy's 1e-2, a warning is + raised. fcs_cutoff_radius: list - cutoff distance for anharmonic force constants(FCs) up to fourth-order FCs. - The default value is [-1, 12, 10], which means that the cutoff distance for - second-order FCs is the Wigner-Seitz cell boundary and the cutoff distance - for third-order FCs is 12 Borh, and the cutoff distance for fourth-order FCs - is 10 Bohr. Generally, the default value is good enough. + 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. min_length: float minimum length of lattice constants will be used to create the supercell, - the default value is 14.0 A. In most cases, the default value is good - enough, but it can be increased for larger supercells. - prefer_90_degrees: bool - if set to True, supercell algorithm will first try to find a supercell - with 3 90 degree angles. + the default value is 8.0 A. It can be increased for larger supercells. + max_atoms: float | None + maximum number of atoms in the supercell. + force_90_degrees: bool + if set to True, only supercells with three 90 degree angles are allowed. + force_diagonal: bool + if set to True, only diagonal supercell matrices are allowed. The pheasy + commands take the diagonal of the supercell matrix. get_supercell_size_kwargs: dict - kwargs that will be passed to get_supercell_size to determine supercell size + not used by this workflow. use_symmetrized_structure: str allowed strings: "primitive", "conventional", None @@ -174,14 +211,12 @@ class BasePhononMaker(PurePhonopyMaker, ABC): cal_anhar_fcs: bool = False displacement_anhar: float = 0.08 num_disp_anhar: int = 0 + anhar_max_order: int = 3 + anhar_fit_methods: Sequence[Literal["cocktail", "one-shot"]] = ("cocktail",) + anhar_alpha_min: int = -12 fcs_cutoff_radius: list = field( default_factory=lambda: [-1, 12, 10] ) # units in Bohr - renorm_phonon: bool = False - renorm_temp: list = field(default_factory=lambda: [100, 700, 100]) - cal_ther_cond: bool = False - ther_cond_mesh: list = field(default_factory=lambda: [20, 20, 20]) - ther_cond_temp: list = field(default_factory=lambda: [100, 700, 100]) min_length: float | None = 8.0 max_atoms: float | None = 200 force_90_degrees: bool = True @@ -203,6 +238,17 @@ class BasePhononMaker(PurePhonopyMaker, ABC): store_force_constants: bool = True socket: bool = False + def __post_init__(self) -> None: + """Check the anharmonic settings before any calculation is run.""" + _check_anharmonic_settings( + self.anhar_max_order, + self.anhar_fit_methods, + self.cal_anhar_fcs, + self.fcs_cutoff_radius, + self.anhar_alpha_min, + self.num_disp_anhar, + ) + def get_displacements( self, structure: Structure, supercell_matrix: Matrix3D ) -> Job | Flow: @@ -232,6 +278,8 @@ def get_displacements( use_symmetrized_structure=self.use_symmetrized_structure, kpath_scheme=self.kpath_scheme, code=self.code, + anhar_max_order=self.anhar_max_order, + anhar_fit_methods=self.anhar_fit_methods, ) def run_displacements( @@ -309,11 +357,6 @@ def get_results( displacement=self.displacement, cal_anhar_fcs=self.cal_anhar_fcs, fcs_cutoff_radius=self.fcs_cutoff_radius, - renorm_phonon=self.renorm_phonon, - renorm_temp=self.renorm_temp, - cal_ther_cond=self.cal_ther_cond, - ther_cond_mesh=self.ther_cond_mesh, - ther_cond_temp=self.ther_cond_temp, sym_reduce=self.sym_reduce, symprec=self.symprec, use_symmetrized_structure=self.use_symmetrized_structure, @@ -332,6 +375,10 @@ def get_results( optimization_run_uuid=optimization_run_uuid, create_thermal_displacements=self.create_thermal_displacements, store_force_constants=self.store_force_constants, + num_displaced_supercells=self.num_displaced_supercells, + anhar_max_order=self.anhar_max_order, + anhar_fit_methods=self.anhar_fit_methods, + anhar_alpha_min=self.anhar_alpha_min, **self.generate_frequencies_eigenvectors_kwargs, ) diff --git a/src/atomate2/common/jobs/hiphive.py b/src/atomate2/common/jobs/hiphive.py index da27fb44c3..f1b073707e 100644 --- a/src/atomate2/common/jobs/hiphive.py +++ b/src/atomate2/common/jobs/hiphive.py @@ -38,6 +38,8 @@ from atomate2.common.jobs.phonons import ( _generate_phonon_object, _get_kpath, + _get_num_harmonic_supercells, + _get_num_irreducible_fcs, _run_band_structure_and_plot, _run_total_dos_and_plot, ) @@ -48,20 +50,14 @@ logger = logging.getLogger(__name__) -try: - from alm import ALM -except ImportError: - ALM = None - -# Safety margin below hiPhive's maximum allowed cutoff, in Angstrom. -# Configurations per suggested displacement set. ALM sizes num_disp_sc so the -# equation count reaches its own free-parameter tally, and this scales from -# there. pheasy uses 1.8. hiPhive's cluster space has a different, smaller -# parameter count, and the ratio between the two varies by a factor of three -# across materials, so 1.8 leaves a low-symmetry cell with too little data for -# a LASSO fit and the modes come out soft. +# Configurations per suggested displacement set. The minimum set is sized from +# ALM's free-parameter tally so the equation count reaches it (see +# _get_num_harmonic_supercells), and this scales from there. pheasy uses the +# same 1.8. hiPhive's own cluster space has a smaller parameter count than +# ALM's tally, and the ratio between the two differs between materials. _N_CONFIG_MULTIPLIER = 1.8 +# Safety margin below hiPhive's maximum allowed cutoff, in Angstrom. _CUTOFF_MARGIN = 0.1 # hiPhive refuses a cutoff that clears a neighbour shell by less than its @@ -290,11 +286,12 @@ def generate_phonon_displacements( random_seed: int | None = 103, verbose: bool = False, ) -> list[Structure]: - """Generate small-distance perturbed structures with phonopy based on two ways. + """Generate the displaced supercells with phonopy. - 1. finite-displacment method (one displaced atom) when the displacement number - is less than 3. 2. random-displacement method (all-displaced atoms) when the - displacement number is more than 3. + The finite-displacement method (one displaced atom) is used when phonopy + needs at most three displacements, and the random-displacement method (all + atoms displaced) otherwise. The undisplaced supercell is added last, for + the residual forces. Parameters ---------- @@ -322,11 +319,6 @@ def generate_phonon_displacements( Whether to log warnings. """ - # TODO: remove ALMODE dependence for 2nd order force constants - if not ALM: - raise ImportError( - "Error importing ALM. Please ensure the 'alm' library is installed." - ) phonon = _generate_phonon_object( structure, supercell_matrix, @@ -353,46 +345,28 @@ def generate_phonon_displacements( # of the matrix can not always guarantee accurate results, you # may need to displace more random configurations. Use at least one or # two more configurations based on the suggested number of displacements. - supercell_ph = phonon.supercell - lattice = supercell_ph.cell - positions = supercell_ph.scaled_positions - numbers = supercell_ph.numbers - natom = len(numbers) - - # get the number of free parameters of 2ND FCs from ALM, labeled as n_fp - with ALM(lattice, positions, numbers) as alm: - alm.define(1) - alm.suggest() - n_fp = alm._get_number_of_irred_fc_elements(1) # noqa: SLF001 - - # get the number of displaced supercells based on the number of free parameters - num_disp_sc = int(np.ceil(n_fp / (3.0 * natom))) + num_har = _get_num_harmonic_supercells( + phonon, num_displaced_supercells, multiplier=_N_CONFIG_MULTIPLIER + ) if verbose: + (n_fp,) = _get_num_irreducible_fcs(phonon.supercell, 2) logger.info( - f"There are {n_fp} free parameters for the second-order " - "force constants (FCs)." - f"There are {3 * natom * num_disp_sc} equations used to " - "obtain the second-order FCs." - "CAUTION: you may need to increase the number of " - "displacements in some cases." - "If the number of atoms in the supercell are less than 100 and " + f"There are {n_fp} free parameters for the second-order force " + f"constants (FCs), and {num_har} displaced supercells are used to " + "fit them. CAUTION: you may need to increase the number of " + "displacements in some cases. " + "If the number of atoms in the supercell is less than 100 and " "all lattice constants are less than 10 Å, the user is advised " "to use 1-2 more randomly-displaced configurations." ) - # get the number of displaced supercells from phonopy to compared with the number - # of 3, if the number of displaced supercells is less than 3, we will use the finite - # displacement method to generate the supercells. Otherwise, we will use the random - # displacement method to generate the supercells. + # if phonopy needs more than three finite displacements, we use the random + # displacement method instead if len(phonon.displacements) > 3: phonon.generate_displacements( distance=displacement, - number_of_snapshots=( - num_displaced_supercells - if num_displaced_supercells != 0 - else int(np.ceil(num_disp_sc * _N_CONFIG_MULTIPLIER)) + 1 - ), + number_of_snapshots=num_har, random_seed=random_seed, ) diff --git a/src/atomate2/common/jobs/pheasy.py b/src/atomate2/common/jobs/pheasy.py index 29c80cfb58..eb2648b817 100644 --- a/src/atomate2/common/jobs/pheasy.py +++ b/src/atomate2/common/jobs/pheasy.py @@ -3,8 +3,12 @@ from __future__ import annotations import logging +import numbers +import re import shlex +import shutil import subprocess +import warnings from pathlib import Path from typing import TYPE_CHECKING @@ -18,7 +22,7 @@ from hiphive.utilities import extract_parameters from jobflow import job from packaging.version import parse as parse_version -from phonopy.file_IO import parse_FORCE_CONSTANTS, write_force_constants_to_hdf5 +from phonopy.file_IO import parse_FORCE_CONSTANTS from phonopy.interface.vasp import write_vasp from phonopy.structure.symmetry import symmetrize_borns_and_epsilon from pymatgen.core import Structure @@ -31,25 +35,23 @@ ) from atomate2.common.jobs.phonons import ( + ANGSTROM_TO_BOHR, _generate_phonon_object, _get_kpath, + _get_num_harmonic_supercells, + _get_num_irreducible_fcs, _run_band_structure_and_plot, _run_total_dos_and_plot, ) if TYPE_CHECKING: + from collections.abc import Sequence + from emmet.core.math import Matrix3D + from phonopy.structure.atoms import PhonopyAtoms logger = logging.getLogger(__name__) -try: - from alm import ALM -except ImportError: - ALM = None - -# CODATA 2018: 1 Angstrom = 1 / 0.529177210903 Bohr -ANGSTROM_TO_BOHR = 1.8897261246257702 - _DEFAULT_FILE_PATHS = { "force_displacements": "dataset_forces.npy", "displacements": "dataset_disps.npy", @@ -65,8 +67,283 @@ "harmonic_force_matrix": "force_matrix.npy", "anharmonic_force_matrix": "force_matrix_anhar.npy", "website": "phonon_website.json", + "one_shot_dir": "one_shot", + "anharmonic_fit_log": "pheasy_anharmonic_fit.log", } +# The anharmonic training set is sized so that the fit has this many force +# equations per free force constant. The floor keeps the set from becoming +# very small. Above the ceiling the job stops instead of requesting more +# displaced supercells. +_EQUATIONS_PER_FREE_FC = 100 +_MIN_NUM_DISP_ANHAR = 20 +_MAX_NUM_DISP_ANHAR = 600 + +_ANHARMONIC_FIT_METHODS = ("cocktail", "one-shot") + +# many-body terms kept in the anharmonic fit for each maximum order, as passed +# to pheasy with --nbody and to ALM when counting the free force constants +_NBODY = {3: [2, 3], 4: [2, 3, 3]} + + +def _get_num_anharmonic_supercells( + supercell: PhonopyAtoms, + num_disp_anhar: int, + anhar_max_order: int, + fcs_cutoff_radius: Sequence[float], + anhar_fit_methods: Sequence[str], +) -> int: + """ + Get the number of randomly displaced supercells for the anharmonic fit. + + A non-zero num_disp_anhar is returned as given. Otherwise, each supercell + gives 3 * natom force equations, and the number of supercells is set so + that there are _EQUATIONS_PER_FREE_FC equations per free force constant, + with at least _MIN_NUM_DISP_ANHAR supercells. Above _MAX_NUM_DISP_ANHAR a + ValueError is raised. The free force constants are those of third (and + fourth) order, plus those of second order when the one-shot fit is + requested, since it fits all orders together. Second-order force constants + are counted without a cutoff, as pheasy fits them. ALM counts the + irreducible force constants before the acoustic sum rules are applied, with + its own symmetry search, so the count is expected to be at least the number + pheasy fits. + + Parameters + ---------- + supercell: PhonopyAtoms + Supercell used for the force constant fit. + num_disp_anhar: int + Number of displaced supercells requested by the user, 0 for automatic. + anhar_max_order: int + Highest force constant order in the anharmonic fit, 3 or 4. + fcs_cutoff_radius: Sequence[float] + Cutoff radius in Bohr for each order, starting at second order. + anhar_fit_methods: Sequence[str] + Anharmonic fit methods that will use this dataset. + + Returns + ------- + int + Number of anharmonic displaced supercells. + """ + if num_disp_anhar != 0: + return num_disp_anhar + # pheasy fits the second-order force constants without a cutoff + cutoffs = [-1, *fcs_cutoff_radius[1:]] + n_irred = _get_num_irreducible_fcs( + supercell, anhar_max_order, cutoffs, nbody=_NBODY[anhar_max_order] + ) + n_free = sum(n_irred[1:]) + if "one-shot" in anhar_fit_methods: + n_free += n_irred[0] + natom = len(supercell.numbers) + num = int(np.ceil(_EQUATIONS_PER_FREE_FC * n_free / (3.0 * natom))) + num = max(num, _MIN_NUM_DISP_ANHAR) + if num > _MAX_NUM_DISP_ANHAR: + raise ValueError( + f"{num} displaced supercells are needed for {n_free} free force " + f"constants, more than the limit of {_MAX_NUM_DISP_ANHAR}. Reduce " + "fcs_cutoff_radius, use a larger supercell, or set num_disp_anhar " + "to accept the cost." + ) + return num + + +def _check_anharmonic_settings( + anhar_max_order: int, + anhar_fit_methods: Sequence[str], + cal_anhar_fcs: bool = True, + fcs_cutoff_radius: Sequence[float] | None = None, + anhar_alpha_min: int | None = None, + num_disp_anhar: int = 0, +) -> None: + """ + Check the settings of the anharmonic force constant fit. + + Parameters + ---------- + anhar_max_order: int + Highest force constant order in the anharmonic fit. + anhar_fit_methods: Sequence[str] + Anharmonic fit methods. + cal_anhar_fcs: bool + Whether the anharmonic force constants are calculated. + fcs_cutoff_radius: Sequence[float] | None + Cutoff radius in Bohr for each order, starting at second order. + anhar_alpha_min: int | None + Base-10 exponent of the smallest LASSO penalty. pheasy only accepts an + integer, below its largest penalty of 1e-2. + num_disp_anhar: int + Number of anharmonic displaced supercells requested by the user. + """ + if anhar_max_order not in (3, 4): + raise ValueError(f"anhar_max_order must be 3 or 4, not {anhar_max_order}.") + unknown = set(anhar_fit_methods) - set(_ANHARMONIC_FIT_METHODS) + if ( + unknown + or len(anhar_fit_methods) == 0 + or len(set(anhar_fit_methods)) != len(anhar_fit_methods) + ): + raise ValueError( + f"anhar_fit_methods must be a non-empty subset of " + f"{_ANHARMONIC_FIT_METHODS} without repeats, not {list(anhar_fit_methods)}." + ) + if ( + isinstance(num_disp_anhar, bool) + or not isinstance(num_disp_anhar, numbers.Integral) + or num_disp_anhar < 0 + ): + raise ValueError( + f"num_disp_anhar must be a non-negative integer, not {num_disp_anhar!r}." + ) + if cal_anhar_fcs and fcs_cutoff_radius is not None: + # one cutoff per order from 3 to anhar_max_order, after the fc2 entry + anhar_radii = list(fcs_cutoff_radius[1 : anhar_max_order - 1]) + # pheasy reads a negative cutoff as a neighbour shell, ALM as no cutoff + if len(anhar_radii) != anhar_max_order - 2 or min(anhar_radii) <= 0: + raise ValueError( + f"fcs_cutoff_radius needs a positive cutoff in Bohr for each order " + f"from 3 to {anhar_max_order}, not {list(fcs_cutoff_radius)}." + ) + if anhar_alpha_min is not None and ( + isinstance(anhar_alpha_min, bool) + or not isinstance(anhar_alpha_min, numbers.Integral) + or anhar_alpha_min >= -2 + ): + raise ValueError( + f"anhar_alpha_min must be an integer below -2, not {anhar_alpha_min!r}." + ) + + +def _check_lasso_alpha(log_file: Path, alpha_min: int) -> None: + """ + Warn if the cross-validated LASSO penalty is on either bound of the search. + + On the lower bound, the search wanted an even weaker penalty than it was + allowed to try, so the fitted force constants depend on that bound. On the + upper bound, the penalty is the strongest one tried. Raises an error if + pheasy wrote no log or no penalty. + + Parameters + ---------- + log_file: Path + Log file written by pheasy during the fit. + alpha_min: int + Base-10 exponent of the smallest penalty in the search. + """ + if not log_file.exists(): + raise FileNotFoundError( + f"pheasy exited without writing {log_file}, so the result of the " + "anharmonic fit is unknown." + ) + match = re.search(r"alpha_opt:\s*([-+0-9.eE]+)", log_file.read_text()) + if match is None: + raise RuntimeError(f"No LASSO alpha found in {log_file}, the fit failed.") + alpha_opt = float(match.group(1)) + match_max = re.search(r"alpha_max:\s*1e([-+0-9]+)", log_file.read_text()) + if match_max is not None and alpha_opt >= 10.0 ** int(match_max.group(1)) * ( + 1 - 1e-6 + ): + warnings.warn( + f"The LASSO penalty chosen by cross-validation ({alpha_opt:e}) is on " + f"the upper bound in {log_file}. The fitted force constants may be " + "heavily penalized. Check them.", + stacklevel=2, + ) + 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.", + stacklevel=2, + ) + + +def _run_anharmonic_fit( + method: str, + supercell_matrix: np.ndarray, + symprec: float, + anhar_max_order: int, + fcs_cutoff_radius: Sequence[float], + num_anhar: int, + anhar_alpha_min: int, + work_dir: Path, + random_seed: int | None = 103, +) -> None: + """ + Fit the anharmonic force constants with pheasy using LASSO. + + Both methods use the same randomly displaced supercells. + + - "cocktail" keeps the second-order force constants fixed to the harmonic + fit (FORCE_CONSTANTS in work_dir) and fits only the higher orders. + - "one-shot" fits the second- and higher-order force constants together. + + pheasy writes FORCE_CONSTANTS_3RD and fc3.hdf5 to work_dir, and + FORCE_CONSTANTS_4TH and fc4.hdf5 for fourth order. The one-shot fit also + writes its second-order force constants, fc2.hdf5. + + Parameters + ---------- + method: str + "cocktail" or "one-shot". + supercell_matrix: np.ndarray + Diagonal supercell matrix. + symprec: float + Symmetry precision, the same one used for the displacements. + anhar_max_order: int + Highest force constant order, 3 or 4. + fcs_cutoff_radius: Sequence[float] + Cutoff radius in Bohr for each order, starting at second order. + num_anhar: int + Number of anharmonic displaced supercells. + anhar_alpha_min: int + Base-10 exponent of the smallest LASSO penalty in the search. + work_dir: Path + Folder holding POSCAR and the anharmonic displacement and force + matrices. + random_seed: int | None + Seed for the random coordinate descent in pheasy's LASSO fit. Without + it, the cross-validated penalty and the fitted force constants change + from run to run on the same forces. + """ + dim = " ".join(str(int(supercell_matrix[i][i])) for i in range(3)) + base = f"pheasy --dim {dim} -w {anhar_max_order} --symprec {float(symprec)}" + cutoffs = f"--c3 {float(fcs_cutoff_radius[1] / ANGSTROM_TO_BOHR)}" + if anhar_max_order == 4: + cutoffs += f" --c4 {float(fcs_cutoff_radius[2] / ANGSTROM_TO_BOHR)}" + nbody = "--nbody " + " ".join(str(n) for n in _NBODY[anhar_max_order]) + fix_fc2 = "--fix_fc2 --fc2_fmt PHONOPY " if method == "cocktail" else "" + log_file = _DEFAULT_FILE_PATHS["anharmonic_fit_log"] + seed = f"--seed {int(random_seed)} " if random_seed is not None else "" + + commands = [ + # clusters and orbits + f"{base} -s {nbody} {cutoffs}", + # null space + f"{base} -c", + # sensing matrix from the displacement matrix + ( + f"{base} -d --ndata {int(num_anhar)} --disp_file " + f"--disp_matrix_file {_DEFAULT_FILE_PATHS['anharmonic_displacements']}" + ), + # LASSO fit. OLS, pheasy's default, gives dense force constants. + # --std and --rasr are not passed to the anharmonic fit. + ( + f"{base} -f {fix_fc2}-l LASSO --alpha_min {anhar_alpha_min} {seed}" + f"--ndata {int(num_anhar)} --hdf5 " + f"--force_matrix_file {_DEFAULT_FILE_PATHS['anharmonic_force_matrix']} " + f"-o {log_file}" + ), + ] + # pheasy appends to its log, so remove the log of an earlier fit + (work_dir / log_file).unlink(missing_ok=True) + for cmd in commands: + subprocess.run(shlex.split(cmd), cwd=work_dir, check=True) + + _check_lasso_alpha(work_dir / log_file, anhar_alpha_min) + @job def get_supercell_size( @@ -77,7 +354,7 @@ def get_supercell_size( force_diagonal: bool, ) -> list[list[float]]: """ - Determine supercell size with given min_length and max_length. + Determine the supercell matrix with pymatgen's CubicSupercellTransformation. Parameters ---------- @@ -85,14 +362,12 @@ def get_supercell_size( Input structure that will be used to determine supercell min_length: float minimum length of cell in Angstrom - max_length: float - maximum length of cell in Angstrom - prefer_90_degrees: bool - if True, the algorithm will try to find a cell with 90 degree angles first - allow_orthorhombic: bool - if True, orthorhombic supercells are allowed - **kwargs: - Additional parameters that can be set. + max_atoms: int + maximum number of atoms in the supercell + force_90_degrees: bool + if True, only supercells with three 90 degree angles are allowed + force_diagonal: bool + if True, only diagonal supercell matrices are allowed """ transformation = CubicSupercellTransformation( min_length=min_length, @@ -115,21 +390,25 @@ def generate_phonon_displacements( cal_anhar_fcs: bool, displacement_anhar: float, num_disp_anhar: int, - fcs_cutoff_radius: list[int], + fcs_cutoff_radius: list[float], sym_reduce: bool, symprec: float, use_symmetrized_structure: str | None, kpath_scheme: str, code: str, + anhar_max_order: int = 3, + anhar_fit_methods: Sequence[str] = ("cocktail",), random_seed: int | None = 103, verbose: bool = False, ) -> list[Structure]: - """Generate small-distance perturbed structures with phonopy based on two ways. + """Generate the displaced supercells with phonopy. - (we will directly use the pheasy to generate the supercell in the near future) - 1. finite-displacment method (one displaced atom) when the displacement number - is less than 3. 2. random-displacement method (all-displaced atoms) when the - displacement number is more than 3. + The small-distance set for the harmonic force constants uses the + finite-displacement method (one displaced atom) when phonopy needs at most + three displacements, and the random-displacement method (all atoms + displaced) otherwise. With cal_anhar_fcs, a large-distance set of randomly + displaced supercells is added for the anharmonic force constants. The + undisplaced supercell is added last, for the residual forces. Parameters ---------- @@ -138,13 +417,20 @@ def generate_phonon_displacements( supercell_matrix: np.array array to describe supercell matrix displacement: float - displacement in Angstrom (default: 0.01) + displacement in Angstrom num_displaced_supercells: int - number of displaced supercells defined by users + number of harmonic random displacements, 0 for automatic. Not used when + phonopy needs at most three finite displacements. cal_anhar_fcs: bool - TODO : docstr + if True, also generate the large-distance displacements used for the + anharmonic force constants displacement_anhar: float - TODO : docstr + displacement in Angstrom for the anharmonic force constants + num_disp_anhar: int + number of anharmonic displaced supercells, 0 to set it from the number + of free force constants + fcs_cutoff_radius: list[float] + cutoff radius in Bohr for each force constant order, from second order sym_reduce: bool if True, symmetry will be used to generate displacements symprec: float @@ -155,17 +441,17 @@ def generate_phonon_displacements( scheme to generate kpath code: str code to perform the computations + anhar_max_order: int + highest anharmonic force constant order, 3 or 4 + anhar_fit_methods: Sequence[str] + anharmonic fit methods, used to size the anharmonic dataset random_seed : int | None = 103 - Random seed to use in generating randomly-displaced structures. + Random seed for the harmonic random displacements. The anharmonic set + uses random_seed + 1. verbose : bool = False - Whether to log warnings. + Whether to log warnings and the numbers of displaced supercells. """ - # TODO: remove ALMODE dependence for 2nd order force constants - if not ALM: - raise ImportError( - "Error importing ALM. Please ensure the 'alm' library is installed." - ) phonon = _generate_phonon_object( structure, supercell_matrix, @@ -178,91 +464,58 @@ def generate_phonon_displacements( verbose=verbose, ) - # 1. the ALM module is used to determine the number of free parameters - # (irreducible force constants) corresponding to the second order - # force constants (FCs) given a supercell. - # 2. Based on the number of free parameters, we can determine how many - # displaced supercells we need to use to extract the second order force - # constants. Generally, the number of free parameters should be less than - # 3 * natom(supercell) * num_displaced_supercells. However, the full rank - # of the matrix can not always guarantee accurate results, you - # may need to displace more random configurations. Use at least one or - # two more configurations based on the suggested number of displacements. supercell_ph = phonon.supercell - lattice = supercell_ph.cell - positions = supercell_ph.scaled_positions - numbers = supercell_ph.numbers - natom = len(numbers) - - # get the number of free parameters of 2ND FCs from ALM, labeled as n_fp - with ALM(lattice, positions, numbers) as alm: - alm.define(1) - alm.suggest() - n_fp = alm._get_number_of_irred_fc_elements(1) # noqa: SLF001 - - # get the number of displaced supercells based on the number of free parameters - num_disp_sc = int(np.ceil(n_fp / (3.0 * natom))) + num_har = _get_num_harmonic_supercells(phonon, num_displaced_supercells) if verbose: + (n_fp,) = _get_num_irreducible_fcs(supercell_ph, 2) logger.info( - f"There are {n_fp} free parameters for the second-order " - "force constants (FCs)." - f"There are {3 * natom * num_disp_sc} equations used to " - "obtain the second-order FCs." - "CAUTION: you may need to increase the number of " - "displacements in some cases." - "If the number of atoms in the supercell are less than 100 and " - "all lattice constants are less than 10 Å, the user is advised " - "to use 1-2 more randomly-displaced configurations." + f"There are {n_fp} free parameters for the second-order force " + f"constants (FCs), and {num_har} displaced supercells are used to " + "fit them. CAUTION: you may need to increase the number " + "of displacements in some cases. If the number of atoms in the " + "supercell is less than 100 and all lattice constants are less than " + "10 Å, the user is advised to use 1-2 more randomly-displaced " + "configurations." ) - # get the number of displaced supercells from phonopy to compared with the number - # of 3, if the number of displaced supercells is less than 3, we will use the finite - # displacement method to generate the supercells. Otherwise, we will use the random - # displacement method to generate the supercells. + # if phonopy needs more than three finite displacements, we use the random + # displacement method instead if len(phonon.displacements) > 3: phonon.generate_displacements( distance=displacement, - number_of_snapshots=( - num_displaced_supercells - if num_displaced_supercells != 0 - else int(np.ceil(num_disp_sc * 1.8)) + 1 - ), + number_of_snapshots=num_har, random_seed=random_seed, ) supercells = phonon.supercells_with_displacements displacements = [get_pmg_structure(cell) for cell in supercells] - # Here, the ALAMODE code is used to determine the number of - # third and fourth-order FCs are needed for the supercell if cal_anhar_fcs: - # Due to the cutoff radius of the force constants use the unit of Bohr in ALM, - # we need to convert the cutoff radius from Angstrom to Bohr. - with ALM(lattice * ANGSTROM_TO_BOHR, positions, numbers) as alm: - # Define the force constants up to fourth order with a list of - # cutoff radius - alm.define(3, fcs_cutoff_radius) - # Perform symmetry analysis and suggest irreducible force constants. - alm.suggest() - # Get the number of irreducible elements for both 3RD- and 4TH-order - # force constants - n_rd_anh = alm._get_number_of_irred_fc_elements( # noqa: SLF001 - 2 - ) + alm._get_number_of_irred_fc_elements(3) # noqa: SLF001 - # we can determine how many displaced supercells we need to use to extract - # the 3rd and 4th order force constants, and we can add a scaling factor - # to reduce the number of displaced supercells due to we use the lasso - # technique. - num_d_anh = int(np.ceil(n_rd_anh / (3.0 * natom))) - num_dis_cells_anhar = num_disp_anhar if num_disp_anhar != 0 else num_d_anh - - num_dis_cells_anhar = 20 - # generate the supercells for anharmonic force constants + _check_anharmonic_settings( + anhar_max_order, + anhar_fit_methods, + fcs_cutoff_radius=fcs_cutoff_radius, + num_disp_anhar=num_disp_anhar, + ) + num_dis_cells_anhar = _get_num_anharmonic_supercells( + supercell_ph, + num_disp_anhar, + anhar_max_order, + fcs_cutoff_radius, + anhar_fit_methods, + ) + if verbose: + logger.info( + f"{num_dis_cells_anhar} displaced supercells are used for the " + "anharmonic force constants." + ) + # generate the supercells for anharmonic force constants. A different + # seed keeps them from repeating the directions of a random harmonic set. phonon.generate_displacements( distance=displacement_anhar, number_of_snapshots=num_dis_cells_anhar, - random_seed=random_seed, + random_seed=None if random_seed is None else random_seed + 1, ) supercells = phonon.supercells_with_displacements displacements += [get_pmg_structure(cell) for cell in supercells] @@ -282,11 +535,7 @@ def generate_frequencies_eigenvectors( supercell_matrix: np.array, displacement: float, cal_anhar_fcs: bool, - fcs_cutoff_radius: list[int], - renorm_phonon: bool, - cal_ther_cond: bool, - ther_cond_mesh: list[int], - ther_cond_temp: list[int], + fcs_cutoff_radius: list[float], sym_reduce: bool, symprec: float, use_symmetrized_structure: str | None, @@ -296,10 +545,20 @@ def generate_frequencies_eigenvectors( total_dft_energy: float, epsilon_static: Matrix3D = None, born: Matrix3D = None, + num_displaced_supercells: int = 0, + anhar_max_order: int = 3, + anhar_fit_methods: Sequence[str] = ("cocktail",), + anhar_alpha_min: int = -12, + random_seed: int | None = 103, **kwargs, ) -> PhononBSDOSDoc: """ - Analyze the phonon runs and summarize the results. + Fit the force constants with pheasy and summarize the phonon results. + + The harmonic force constants give the band structure, density of states + and thermodynamic properties in the output document. With cal_anhar_fcs, + the anharmonic force constants are also fitted and written to files in + the job folder. They are not stored in the output document. Parameters ---------- @@ -309,6 +568,10 @@ def generate_frequencies_eigenvectors( array to describe supercell displacement: float displacement in Angstrom used for supercell computation + cal_anhar_fcs: bool + if True, the anharmonic force constants are fitted as well + fcs_cutoff_radius: list[float] + cutoff radius in Bohr for each force constant order, from second order sym_reduce: bool if True, symmetry will be used in phonopy symprec: float @@ -327,11 +590,30 @@ def generate_frequencies_eigenvectors( The high-frequency dielectric constant born: Matrix3D Born charges - verbose : bool = False - Whether to log error messages. + num_displaced_supercells: int + number of harmonic random displacements requested by the user, + 0 for automatic. Must match the value used to generate them. + anhar_max_order: int + highest anharmonic force constant order, 3 or 4 + anhar_fit_methods: Sequence[str] + "cocktail" (fixed second-order FCs), "one-shot" (all orders fitted + together, written to the one_shot folder), or both + anhar_alpha_min: int + base-10 exponent of the smallest LASSO penalty in the anharmonic fits + random_seed : int | None = 103 + Seed for the harmonic and anharmonic LASSO fits, so that they are + reproducible. kwargs: dict - Additional parameters that are passed to PhononBSDOSDoc.from_forces_born + Further options read by this job, such as npoints_band, filename_bs, + filename_dos, kpoint_density_dos and store_force_constants. """ + _check_anharmonic_settings( + anhar_max_order, + anhar_fit_methods, + cal_anhar_fcs, + fcs_cutoff_radius, + anhar_alpha_min, + ) phonon = _generate_phonon_object( structure, supercell_matrix, @@ -405,34 +687,7 @@ def generate_frequencies_eigenvectors( # separate the dataset into harmonic and anharmonic parts num_har = dataset_disps_array_use.shape[0] if cal_anhar_fcs: - if not ALM: - raise ImportError( - "Error importing ALM. Please ensure the 'alm' library is installed." - ) - - supercell_ph = phonon.supercell - lattice = supercell_ph.cell - positions = supercell_ph.scaled_positions - numbers = supercell_ph.numbers - natom = len(numbers) - - # get the number of free parameters of 2ND FCs from ALM, labeled as n_fp - with ALM(lattice, positions, numbers) as alm: - alm.define(1) - alm.suggest() - n_fp = alm._get_number_of_irred_fc_elements(1) # noqa: SLF001 - - # get the number of displaced supercells based on the - # number of free parameters - num = int(np.ceil(n_fp / (3.0 * natom))) - - # get the number of displaced supercells from phonopy to compared - # with the number of 3, if the number of displaced supercells is - # less than 3, we will use the finite displacement method to generate - # the supercells. Otherwise, we will use the random displacement - # method to generate the supercells. - num_disp_f = len(phonon.displacements) - num_har = int(np.ceil(num * 1.8)) if num_disp_f > 3 else num_disp_f + num_har = _get_num_harmonic_supercells(phonon, num_displaced_supercells) np.save( _DEFAULT_FILE_PATHS["harmonic_displacements"], @@ -475,11 +730,8 @@ def generate_frequencies_eigenvectors( prim = ase_read("POSCAR") supercell = ase_read("SPOSCAR") - # Create the clusters and orbitals for second order force constants - # For the variables: --w, --nbody, they are used to specify the order of the - # force constants. in the near future, we will add the option to specify the - # order of the force constants. And these two variables can be defined by the - # users. + # 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 --dim {int(supercell_matrix[0][0])} " f"{int(supercell_matrix[1][1])} " @@ -523,6 +775,7 @@ def generate_frequencies_eigenvectors( 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: @@ -543,12 +796,14 @@ def generate_frequencies_eigenvectors( subprocess.call(shlex.split(pheasy_cmd_3)) subprocess.call(shlex.split(pheasy_cmd_4)) - # When this code is run on Github tests, it is failing because it is - # not able to find the FORCE_CONSTANTS file. This is because the file is - # somehow getting generated in some temp directory. Can you fix the bug? fc_file = Path(_DEFAULT_FILE_PATHS["force_constants"]) + if cal_anhar_fcs and not fc_file.exists(): + raise RuntimeError( + "The harmonic pheasy fit did not write FORCE_CONSTANTS, so the " + "anharmonic force constants cannot be fitted." + ) - if cal_anhar_fcs and fc_file.exists(): + if cal_anhar_fcs: np.save( _DEFAULT_FILE_PATHS["anharmonic_displacements"], dataset_disps_array_use[num_har:, :, :], @@ -559,85 +814,31 @@ def generate_frequencies_eigenvectors( ) num_anhar = dataset_forces_array_disp.shape[0] - num_har - # We next begin to generate the anharmonic force constants up to fourth - # order using the LASSO method - pheasy_cmd_5 = ( - f"pheasy --dim {int(supercell_matrix[0][0])} " - f"{int(supercell_matrix[1][1])} " - f"{int(supercell_matrix[2][2])} -s -w 4 --symprec " - f"{float(symprec)} " - f"--nbody 2 3 3 --c3 {float(fcs_cutoff_radius[1] / ANGSTROM_TO_BOHR)} " - f"--c4 {float(fcs_cutoff_radius[2] / ANGSTROM_TO_BOHR)}" - ) - - pheasy_cmd_6 = ( - f"pheasy --dim {int(supercell_matrix[0][0])} " - f"{int(supercell_matrix[1][1])} " - f"{int(supercell_matrix[2][2])} -c --symprec " - f"{float(symprec)} -w 4" - ) - pheasy_cmd_7 = ( - f"pheasy --dim {int(supercell_matrix[0][0])} " - f"{int(supercell_matrix[1][1])} " - f"{int(supercell_matrix[2][2])} -w 4 -d --symprec " - f"{float(symprec)} " - f"--ndata {int(num_anhar)} --disp_file " - f"--disp_matrix_file {_DEFAULT_FILE_PATHS['anharmonic_displacements']}" - ) - pheasy_cmd_8 = ( - f"pheasy --dim {int(supercell_matrix[0][0])} " - f"{int(supercell_matrix[1][1])} " - f"{int(supercell_matrix[2][2])} -f -w 4 --fix_fc2 " - f"--symprec {float(symprec)} " - f"--ndata {int(num_anhar)} " - f"--force_matrix_file {_DEFAULT_FILE_PATHS['anharmonic_force_matrix']}" - ) - - subprocess.call(shlex.split(pheasy_cmd_5)) - subprocess.call(shlex.split(pheasy_cmd_6)) - subprocess.call(shlex.split(pheasy_cmd_7)) - subprocess.call(shlex.split(pheasy_cmd_8)) - - # begin to renormzlize the phonon energies - if renorm_phonon: - pheasy_cmd_9 = ( - f"pheasy --dim {int(supercell_matrix[0][0])} " - f"{int(supercell_matrix[1][1])} " - f"{int(supercell_matrix[2][2])} -f -w 4 --fix_fc2 " - f"--hdf5 --symprec {float(symprec)} " - f"--ndata {int(num_anhar)}" - ) - - subprocess.call(shlex.split(pheasy_cmd_9)) - - # write the born charges and dielectric constant to the pheasy format - - # begin to convert the force constants to the phonopy and phono3py format - # for the further lattice thermal conductivity calculations - if cal_ther_cond and fc_file.exists(): - # convert the 2ND order force constants to the phonopy format - fc_phonopy_text = parse_FORCE_CONSTANTS(filename=fc_file) - write_force_constants_to_hdf5(fc_phonopy_text, filename="fc2.hdf5") - - # convert the 3RD order force constants to the phonopy format - - prim_hiphive = ase_read("POSCAR") - supercell_hiphive = ase_read("SPOSCAR") - fcs = HiPhiveForceConstants.read_shengBTE( - supercell_hiphive, "FORCE_CONSTANTS_3RD", prim_hiphive - ) - fcs.write_to_phono3py("fc3.hdf5") - - phono3py_cmd = ( - f"phono3py --dim {int(supercell_matrix[0][0])} " - f"{int(supercell_matrix[1][1])} {int(supercell_matrix[2][2])} " - f"--fc2 --fc3 --br --isotope --wigner " - f"--mesh {ther_cond_mesh[0]} {ther_cond_mesh[1]} {ther_cond_mesh[2]} " - f"--tmin {ther_cond_temp[0]} --tmax {ther_cond_temp[1]} " - f"--tstep {ther_cond_temp[2]}" - ) - - subprocess.call(shlex.split(phono3py_cmd)) + # The cocktail fit runs here, next to the harmonic FORCE_CONSTANTS it + # keeps fixed. The one-shot fit runs in its own folder, because both + # fits write FORCE_CONSTANTS_3RD and fc3.hdf5. + for method in anhar_fit_methods: + work_dir = Path.cwd() + if method == "one-shot": + work_dir = Path(_DEFAULT_FILE_PATHS["one_shot_dir"]).resolve() + work_dir.mkdir(exist_ok=True) + for filename in ( + "POSCAR", + _DEFAULT_FILE_PATHS["anharmonic_displacements"], + _DEFAULT_FILE_PATHS["anharmonic_force_matrix"], + ): + shutil.copy(filename, work_dir / filename) + _run_anharmonic_fit( + method=method, + supercell_matrix=supercell_matrix, + symprec=symprec, + anhar_max_order=anhar_max_order, + fcs_cutoff_radius=fcs_cutoff_radius, + num_anhar=num_anhar, + anhar_alpha_min=anhar_alpha_min, + work_dir=work_dir, + random_seed=random_seed, + ) if fc_file.exists(): # Read the force constants from the output file of pheasy code @@ -678,10 +879,8 @@ def generate_frequencies_eigenvectors( # If imaginary modes are present, we first use the hiphive code to enforce # some symmetry constraints to eliminate the imaginary modes (generally work # for small imaginary modes near Gamma point). If the imaginary modes are - # still present, we will use the pheasy code to generate the force constants - # using a shorter cutoff (10 A) to eliminate the imaginary modes, also we - # just want to remove the imaginary modes near Gamma point. In the future, - # we will only use the pheasy code to do the job. + # still present, a pheasy refit with a shorter cutoff (10 A) follows, but + # its result is currently not used (see the NOTE below). if imaginary_modes: # Define a cluster space using the largest cutoff you can @@ -727,7 +926,11 @@ def generate_frequencies_eigenvectors( ) # Using a shorter cutoff (10 A) to generate the force constants to - # eliminate the imaginary modes near Gamma point in pheasy code + # eliminate the imaginary modes near Gamma point in pheasy code. + # NOTE: the force constants are read back below from + # FORCE_CONSTANTS_short_cutoff, the file the hiPhive step above wrote, so + # the result of this pheasy refit is not used. If the refit succeeds, it + # overwrites FORCE_CONSTANTS in the job folder. if imaginary_modes: pheasy_cmd_11 = ( f"pheasy --dim {int(supercell_matrix[0][0])} " diff --git a/src/atomate2/common/jobs/phonons.py b/src/atomate2/common/jobs/phonons.py index d013269ba7..a148f047d0 100644 --- a/src/atomate2/common/jobs/phonons.py +++ b/src/atomate2/common/jobs/phonons.py @@ -44,7 +44,10 @@ from atomate2.vasp.jobs.base import BaseVaspMaker if TYPE_CHECKING: + from collections.abc import Sequence + from emmet.core.math import Matrix3D + from phonopy.structure.atoms import PhonopyAtoms from atomate2.aims.jobs.base import BaseAimsMaker from atomate2.ase.jobs import AseRelaxMaker @@ -52,6 +55,114 @@ logger = logging.getLogger(__name__) +# CODATA 2018: 1 Angstrom = 1 / 0.529177210903 Bohr +ANGSTROM_TO_BOHR = 1.8897261246257702 + + +def _get_num_irreducible_fcs( + supercell: PhonopyAtoms, + max_order: int, + fcs_cutoff_radius: Sequence[float] | None = None, + nbody: Sequence[int] | None = None, +) -> list[int]: + """ + Count the symmetry-irreducible force constants of each order with ALM. + + Used by the pheasy and hiPhive workflows to size their displacement sets. + + Parameters + ---------- + supercell: PhonopyAtoms + Supercell used for the force constant fit. + max_order: int + Highest force constant order to count, 2, 3 or 4. + fcs_cutoff_radius: Sequence[float] | None + Cutoff radius in Bohr for each order, starting at second order. + Only needed when max_order is larger than 2. + nbody: Sequence[int] | None + Largest number of atoms in a cluster for each order, starting at second + order, as in ALM's define. None keeps all clusters. + + Returns + ------- + list[int] + Number of irreducible force constants for each order, starting at + second order. + """ + try: + from alm import ALM + except ImportError as exc: + raise ImportError( + "Error importing ALM. Please ensure the 'alm' library is installed." + ) from exc + + positions = supercell.scaled_positions + numbers = supercell.numbers + + if max_order == 2: + # TODO: remove ALMODE dependence for 2nd order force constants + # the harmonic count uses no cutoff, so the lattice stays in Angstrom + with ALM(supercell.cell, positions, numbers) as alm: + alm.define(1) + alm.suggest() + return [alm._get_number_of_irred_fc_elements(1)] # noqa: SLF001 + + # ALM expects the lattice in Bohr when cutoff radii in Bohr are given, and + # one cutoff per order and per pair of elements + n_elements = len(set(numbers)) + cutoff_radii = [ + np.full((n_elements, n_elements), radius) + for radius in fcs_cutoff_radius[: max_order - 1] + ] + with ALM(supercell.cell * ANGSTROM_TO_BOHR, positions, numbers) as alm: + alm.define(max_order - 1, cutoff_radii=cutoff_radii, nbody=nbody) + alm.suggest() + return [ + alm._get_number_of_irred_fc_elements(order) # noqa: SLF001 + for order in range(1, max_order) + ] + + +def _get_num_harmonic_supercells( + phonon: Phonopy, num_displaced_supercells: int, multiplier: float = 1.8 +) -> int: + """ + Get the number of displaced supercells used for the harmonic force constants. + + Used by the pheasy and hiPhive workflows. In pheasy, both the displacement + generation and the fit call this function, so the dataset is always split + into its harmonic and anharmonic parts at the same place. + + If phonopy needs three or fewer finite displacements, these are used as + they are. Otherwise random displacements are used. Their number is either + given by the user or set from the number of free second-order force + constants, n_fp. The minimum is ceil(n_fp / (3 * natom)). Because the full + rank of the matrix does not guarantee an accurate fit, multiplier times + this minimum plus one configuration is used. + + Parameters + ---------- + phonon: Phonopy + Phonopy object holding the finite displacements. + num_displaced_supercells: int + Number of random displacements requested by the user, 0 for automatic. + multiplier: float + Number of configurations per minimum configuration set. + + Returns + ------- + int + Number of harmonic displaced supercells. + """ + if len(phonon.displacements) <= 3: + return len(phonon.displacements) + if num_displaced_supercells != 0: + return num_displaced_supercells + (n_fp,) = _get_num_irreducible_fcs(phonon.supercell, 2) + natom = len(phonon.supercell.numbers) + num_disp_sc = int(np.ceil(n_fp / (3.0 * natom))) + return int(np.ceil(num_disp_sc * multiplier)) + 1 + def _get_kpath( structure: Structure, kpath_scheme: str, symprec: float, **kpath_kwargs diff --git a/src/atomate2/vasp/flows/pheasy.py b/src/atomate2/vasp/flows/pheasy.py index fca5039a00..14a4368c4a 100644 --- a/src/atomate2/vasp/flows/pheasy.py +++ b/src/atomate2/vasp/flows/pheasy.py @@ -18,22 +18,24 @@ @dataclass class PhononMaker(BasePhononMaker): - """Maker to calculate harmonic phonons with LASSO-based ML code Pheasy. + """Maker to calculate harmonic phonons and anharmonic FCs with Pheasy. - Calculate the zero-K harmonic phonons of a material and higher-order FCs. - Initially, a tight structural relaxation is performed to obtain a structure - without forces on the atoms. Subsequently, supercells with all atoms displaced - by a small amplitude (generally using 0.01 A) are generated and accurate forces - are computed for these structures for the second order force constants. With the - help of pheasy (LASSO technique), these forces are then converted into a dynamical - matrix. In this Workflow, we separate the harmonic phonon calculations and - anharmonic force constants calculations. To correct for polarization effects, a + Calculate the zero-K harmonic phonons of a material and, optionally, its + anharmonic FCs. Initially, a tight structural relaxation is performed to obtain + a structure without forces on the atoms. Subsequently, displaced supercells + (generally using 0.01 A) are generated and accurate forces are computed for + them. If phonopy needs more than three finite displacements, all atoms are + displaced randomly and pheasy fits the second order force constants with LASSO. + Otherwise, the finite displacements are used with a least-squares fit. These + force constants are then converted into a dynamical matrix. In this Workflow, + we separate the harmonic phonon calculations and anharmonic force constants + calculations. To correct for polarization effects, a correction of the dynamical matrix based on BORN charges can be performed. Finally, phonon densities of states, phonon band structures and thermodynamic properties are computed. For the anharmonic force constants, the supercells with all atoms displaced by a larger amplitude (generally using 0.08 A) are generated and accurate forces are computed for these structures. With the help of pheasy (LASSO technique), - the third- and fourth-order force constants are extracted at once. + the third-order (and optionally fourth-order) force constants are extracted. .. Note:: It is heavily recommended to symmetrize the structure before passing it to @@ -62,35 +64,66 @@ class PhononMaker(BasePhononMaker): number of displacements to be generated using a random-displacement approach for harmonic phonon calculations. The default value is 0 and the number of displacements is automatically determined by the number of atoms in the - supercell and its space group. + supercell and its space group. Not used when phonopy needs at most three + finite displacements. cal_anhar_fcs: bool - if set to True, anharmonic force constants(FCs) up to fourth-order FCs will - be calculated. The default value is False, and only harmonic phonons will - be calculated. + if set to True, anharmonic force constants(FCs) up to order + anhar_max_order will be calculated. The default value is False, and only + harmonic phonons will be calculated. displacement_anhar: float - displacement distance for anharmonic force constants(FCs) up to fourth-order - FCs, for most cases 0.08 A is a good choice, but it can be increased to 0.1 A. + displacement distance for the anharmonic force constants(FCs). The default + is 0.08 A. num_disp_anhar: int number of displacements to be generated using a random-displacement approach - for anharmonic phonon calculations. The default value is 0 and the number of - displacements is automatically determined by the number of atoms in the - supercell, cutoff distance for anharmonic FCs its space group. generally, - 50 large-distance displacements are enough for most cases. + for anharmonic phonon calculations. A non-zero value is used as given. The + default value is 0, and then the number is set so that the fit has 100 force + equations per free force constant, with at least 20 displaced supercells. + Above 600 the job stops. The free force constants are counted with ALM, which + uses its own symmetry search, for this supercell and fcs_cutoff_radius. The + second-order FCs are included when "one-shot" is requested. + anhar_max_order: int + highest order of the anharmonic FCs, 3 or 4. The default of 3 fits the + third-order FCs. 4 also fits fourth-order FCs. + anhar_fit_methods: Sequence[Literal["cocktail", "one-shot"]] + how the anharmonic FCs are fitted to the large-distance dataset. + "cocktail" keeps the second-order FCs fixed to the harmonic fit and fits + only the higher orders. "one-shot" fits the second- and higher-order FCs + together from the same dataset and writes them to the "one_shot" folder. + Both can be requested. The default is ("cocktail",). The cocktail FCs are + written to the job folder, and the one-shot FCs, including fc2.hdf5, to its + one_shot subfolder: FORCE_CONSTANTS_3RD and fc3.hdf5, plus + FORCE_CONSTANTS_4TH and fc4.hdf5 at fourth order. They are not stored in the + output document. The LASSO fits are seeded, so a run is reproducible. With + few displaced supercells, the one-shot result depends on that seed. + anhar_alpha_min: int + base-10 exponent of the smallest LASSO penalty tried by cross-validation + in the anharmonic fits. It must be an integer below -2. The default of -12 + is below pheasy's own default of -6. If the chosen penalty lands on either + end of the search, 10**anhar_alpha_min or pheasy's 1e-2, a warning is + raised. fcs_cutoff_radius: list - cutoff distance for anharmonic force constants(FCs) up to fourth-order FCs. - The default value is [-1, 12, 10], which means that the cutoff distance for - second-order FCs is the Wigner-Seitz cell boundary and the cutoff distance - for third-order FCs is 12 Borh, and the cutoff distance for fourth-order FCs - is 10 Bohr. Generally, the default value is good enough. + 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. min_length: float minimum length of lattice constants will be used to create the supercell, - the default value is 14.0 A. In most cases, the default value is good - enough, but it can be increased for larger supercells. + the default value is 8.0 A. It can be increased for larger supercells. + max_atoms: float | None + maximum number of atoms in the supercell. + force_90_degrees: bool + if set to True, only supercells with three 90 degree angles are allowed. + force_diagonal: bool + if set to True, only diagonal supercell matrices are allowed. The pheasy + commands take the diagonal of the supercell matrix. prefer_90_degrees: bool - if set to True, supercell algorithm will first try to find a supercell - with 3 90 degree angles. + not used by this workflow, which uses force_90_degrees instead. + allow_orthorhombic: bool + not used by this workflow. get_supercell_size_kwargs: dict - kwargs that will be passed to get_supercell_size to determine supercell size + not used by this workflow. use_symmetrized_structure: str allowed strings: "primitive", "conventional", None diff --git a/tests/common/jobs/test_pheasy.py b/tests/common/jobs/test_pheasy.py new file mode 100644 index 0000000000..3bf0bb7067 --- /dev/null +++ b/tests/common/jobs/test_pheasy.py @@ -0,0 +1,304 @@ +"""Tests for the anharmonic force constant fit of the pheasy workflow. + +Forces come from ASE's EMT potential, so the jobs run end to end without any +DFT reference data. +""" + +import warnings +from pathlib import Path + +import h5py +import numpy as np +import pytest +from ase.build import bulk +from ase.calculators.emt import EMT +from jobflow import run_locally +from phonopy.file_IO import parse_FORCE_CONSTANTS +from pymatgen.core import Lattice, Structure +from pymatgen.io.ase import AseAtomsAdaptor + +import atomate2.common.jobs.pheasy as pheasy_jobs +from atomate2.common.jobs.pheasy import ( + _check_lasso_alpha, + _get_num_anharmonic_supercells, + generate_frequencies_eigenvectors, + generate_phonon_displacements, +) +from atomate2.common.jobs.phonons import ( + _generate_phonon_object, + _get_num_irreducible_fcs, +) + +# fcs_cutoff_radius in Bohr. 8 Bohr (4.2 A) covers the first two neighbour +# shells of fcc Cu (2.55 and 3.61 A) and stays inside the 10.8 A supercell. +FCS_CUTOFF_RADIUS = [-1, 8, 8] + +COMMON_KWARGS = { + "supercell_matrix": [[3, 0, 0], [0, 3, 0], [0, 0, 3]], + "displacement": 0.01, + "sym_reduce": True, + "symprec": 1e-3, + "use_symmetrized_structure": None, + "kpath_scheme": "seekpath", + "code": "vasp", +} + +FIT_KWARGS = { + "total_dft_energy": None, + "static_run_job_dir": None, + "static_run_uuid": None, + "born_run_job_dir": None, + "born_run_uuid": None, + "optimization_run_job_dir": None, + "optimization_run_uuid": None, +} + + +def _emt_displacement_data(structures: list[Structure]) -> dict: + forces = [] + for structure in structures: + atoms = AseAtomsAdaptor.get_atoms(structure) + atoms.calc = EMT() + forces.append(atoms.get_forces().tolist()) + return { + "forces": forces, + "displaced_structures": structures, + "dirs": [None] * len(structures), + "uuids": [None] * len(structures), + } + + +def _max_abs_fc(filename: Path, key: str) -> float: + with h5py.File(filename) as f: + return float(np.abs(f[key][:]).max()) + + +def _cu_structure() -> Structure: + return AseAtomsAdaptor.get_structure(bulk("Cu", "fcc", a=3.61, cubic=True)) + + +@pytest.mark.parametrize("num_displaced_supercells", [0, 3]) +def test_harmonic_and_anharmonic_split(tmp_dir, num_displaced_supercells): + """The fit job must split the forces where the displacement job did. + + Before this was fixed, the fit job used one harmonic supercell fewer than + the displacement job made, and it ignored num_displaced_supercells. + """ + # a distorted Cu cell needs more than three finite displacements, so the + # random-displacement path is used for the harmonic force constants + structure = Structure( + Lattice.orthorhombic(3.6, 3.7, 3.8), + ["Cu"] * 4, + [[0, 0, 0], [0.5, 0.5, 0.02], [0.5, 0.03, 0.5], [0.01, 0.5, 0.5]], + ) + kwargs = {**COMMON_KWARGS, "supercell_matrix": [[2, 0, 0], [0, 2, 0], [0, 0, 2]]} + assert len(_generate_phonon_object(structure, **kwargs).displacements) > 3 + anhar_kwargs = { + "cal_anhar_fcs": True, + "fcs_cutoff_radius": [-1, 6, 6], + "num_displaced_supercells": num_displaced_supercells, + } + + job = generate_phonon_displacements( + structure=structure, + displacement_anhar=0.03, + num_disp_anhar=20, + **anhar_kwargs, + **kwargs, + ) + responses = run_locally(job, create_folders=True, ensure_success=True) + displacements = responses[job.uuid][1].output + # the harmonic set, 20 anharmonic supercells, then the undisplaced supercell + num_har = len(displacements) - 20 - 1 + assert num_har == (num_displaced_supercells or 14) + + job = generate_frequencies_eigenvectors( + structure=structure, + displacement_data=_emt_displacement_data(displacements), + **anhar_kwargs, + **FIT_KWARGS, + **kwargs, + ) + run_locally(job, create_folders=True, ensure_success=True) + + (fit_dir,) = (path.parent for path in Path.cwd().glob("job_*/disp_matrix.npy")) + harmonic = np.load(fit_dir / "disp_matrix.npy") + anharmonic = np.load(fit_dir / "disp_matrix_anhar.npy") + assert harmonic.shape[0] == num_har + assert anharmonic.shape[0] == 20 + + # the anharmonic set is drawn with its own seed, so its directions differ + # from those of the random harmonic set + assert not np.allclose(harmonic[0] / 0.01, anharmonic[0] / 0.03) + + +def test_get_num_anharmonic_supercells(monkeypatch): + phonon = _generate_phonon_object(_cu_structure(), **COMMON_KWARGS) + kwargs = { + "supercell": phonon.supercell, + "anhar_max_order": 3, + "fcs_cutoff_radius": [-1, 12, 10], + "anhar_fit_methods": ["one-shot"], + } + + # a value set by the user is used as it is + assert _get_num_anharmonic_supercells(num_disp_anhar=7, **kwargs) == 7 + + # sized from the free force constants, above the floor of 20 + assert _get_num_irreducible_fcs(phonon.supercell, 3, [-1, 12]) == [25, 324] + assert _get_num_anharmonic_supercells(num_disp_anhar=0, **kwargs) == 108 + + # the first cutoff is not used, since pheasy fits fc2 without a cutoff + fc2_cutoff = {**kwargs, "fcs_cutoff_radius": [5, 12, 10]} + assert _get_num_anharmonic_supercells(num_disp_anhar=0, **fc2_cutoff) == 108 + + # fourth order, with the three-body limit this workflow passes to pheasy + order_4 = {**kwargs, "anhar_max_order": 4, "anhar_fit_methods": ["cocktail"]} + assert _get_num_anharmonic_supercells(num_disp_anhar=0, **order_4) == 224 + + # a binary compound, L1_2 Cu3Au, needs one cutoff per pair of elements + cu3au = 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]], + ) + binary = { + **kwargs, + "supercell": _generate_phonon_object(cu3au, **COMMON_KWARGS).supercell, + } + assert _get_num_anharmonic_supercells(num_disp_anhar=0, **binary) == 337 + + # above the ceiling the job stops + monkeypatch.setattr(pheasy_jobs, "_MAX_NUM_DISP_ANHAR", 30) + with pytest.raises(ValueError, match="more than the limit of 30"): + _get_num_anharmonic_supercells(num_disp_anhar=0, **kwargs) + + +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) + + log_file.write_text("- alpha_min: 1e-12\n- alpha_opt: 1.000000e-12\n") + with pytest.warns(UserWarning, match="on the lower bound"): + _check_lasso_alpha(log_file, alpha_min=-12) + + log_file.write_text("- alpha_max: 1e-2\n- alpha_opt: 1.000000e-02\n") + with pytest.warns(UserWarning, match="on the upper bound"): + _check_lasso_alpha(log_file, alpha_min=-12) + + log_file.write_text("Fitting force constants via the ordinary least-square.\n") + with pytest.raises(RuntimeError, match="No LASSO alpha"): + _check_lasso_alpha(log_file, alpha_min=-12) + + log_file.unlink() + with pytest.raises(FileNotFoundError, match="exited without writing"): + _check_lasso_alpha(log_file, alpha_min=-12) + + +def test_anharmonic_fit_cocktail_and_one_shot(tmp_dir): + structure = _cu_structure() + anhar_kwargs = { + "cal_anhar_fcs": True, + "fcs_cutoff_radius": FCS_CUTOFF_RADIUS, + "anhar_fit_methods": ["cocktail", "one-shot"], + } + + job = generate_phonon_displacements( + structure=structure, + num_displaced_supercells=0, + displacement_anhar=0.03, + num_disp_anhar=0, + **anhar_kwargs, + **COMMON_KWARGS, + ) + responses = run_locally(job, create_folders=True, ensure_success=True) + displacements = responses[job.uuid][1].output + # one finite displacement for fcc Cu, the anharmonic set at its floor of + # 20, then the undisplaced supercell + assert len(displacements) == 1 + 20 + 1 + + job = generate_frequencies_eigenvectors( + structure=structure, + displacement_data=_emt_displacement_data(displacements), + **anhar_kwargs, + **FIT_KWARGS, + **COMMON_KWARGS, + ) + # a LASSO penalty on either bound of the search fails the test + with warnings.catch_warnings(): + warnings.filterwarnings("error", message="The LASSO penalty") + run_locally(job, create_folders=True, ensure_success=True) + + (fit_dir,) = (path.parent for path in Path.cwd().glob("job_*/fc3.hdf5")) + one_shot_dir = fit_dir / "one_shot" + + # both fits used LASSO + for folder in (fit_dir, one_shot_dir): + log = (folder / "pheasy_anharmonic_fit.log").read_text() + assert "coordinate descent LASSO" in log + + # the one-shot fc2 agrees with the harmonic fc2, which it does not overwrite + harmonic_fc2 = np.abs(parse_FORCE_CONSTANTS(fit_dir / "FORCE_CONSTANTS")).max() + assert _max_abs_fc(one_shot_dir / "fc2.hdf5", "fc2") == pytest.approx( + harmonic_fc2, rel=0.05 + ) + + # third-order force constants in eV/A^3. With 20 supercells, the one-shot + # value changed by about 10 percent with the LASSO seed in local runs, + # while the cocktail value did not, so only the cocktail value is pinned. + cocktail_fc3 = _max_abs_fc(fit_dir / "fc3.hdf5", "fc3") + assert cocktail_fc3 == pytest.approx(4.551, rel=0.02) + assert _max_abs_fc(one_shot_dir / "fc3.hdf5", "fc3") == pytest.approx( + cocktail_fc3, rel=0.15 + ) + + +def test_anharmonic_fit_fourth_order(tmp_dir): + """The fourth-order fit runs and writes fc4.""" + structure = _cu_structure() + anhar_kwargs = { + "cal_anhar_fcs": True, + "fcs_cutoff_radius": FCS_CUTOFF_RADIUS, + "anhar_max_order": 4, + "anhar_fit_methods": ["cocktail"], + } + + job = generate_phonon_displacements( + structure=structure, + num_displaced_supercells=0, + # the workflow default. At 0.03 A, cross-validation chose the lowest + # penalty, 1e-12. + displacement_anhar=0.08, + num_disp_anhar=20, + **anhar_kwargs, + **COMMON_KWARGS, + ) + responses = run_locally(job, create_folders=True, ensure_success=True) + displacements = responses[job.uuid][1].output + + job = generate_frequencies_eigenvectors( + structure=structure, + displacement_data=_emt_displacement_data(displacements), + **anhar_kwargs, + **FIT_KWARGS, + **COMMON_KWARGS, + ) + # a LASSO penalty on either bound of the search fails the test + with warnings.catch_warnings(): + warnings.filterwarnings("error", message="The LASSO penalty") + run_locally(job, create_folders=True, ensure_success=True) + + (fit_dir,) = (path.parent for path in Path.cwd().glob("job_*/fc4.hdf5")) + assert (fit_dir / "FORCE_CONSTANTS_4TH").exists() + assert ( + "coordinate descent LASSO" + in (fit_dir / "pheasy_anharmonic_fit.log").read_text() + ) + + # force constants in eV/A^3 and eV/A^4 + assert _max_abs_fc(fit_dir / "fc3.hdf5", "fc3") == pytest.approx(4.601, rel=0.02) + assert _max_abs_fc(fit_dir / "fc4.hdf5", "fc4") == pytest.approx(48.92, rel=0.05) diff --git a/tests/vasp/flows/test_pheasy.py b/tests/vasp/flows/test_pheasy.py index 81effe8249..d057def21f 100644 --- a/tests/vasp/flows/test_pheasy.py +++ b/tests/vasp/flows/test_pheasy.py @@ -154,3 +154,68 @@ def test_pheasy_wf_vasp(mock_vasp, clean_dir, si_structure: Structure, test_dir) assert thermo_props["heat_capacity"][1] == pytest.approx(20.0, abs=2.5) # the high-T limit must approach (but not exceed) Dulong-Petit (3R per atom) assert 23.0 < thermo_props["heat_capacity"][2] < 3 * 8.3145 + + +@pytest.mark.parametrize( + ("settings", "match"), + [ + ({"cal_anhar_fcs": True, "anhar_max_order": 5}, "anhar_max_order"), + ({"cal_anhar_fcs": True, "anhar_fit_methods": ["joint"]}, "anhar_fit_methods"), + ({"cal_anhar_fcs": True, "anhar_fit_methods": []}, "anhar_fit_methods"), + ({"cal_anhar_fcs": True, "fcs_cutoff_radius": [-1, -2, 10]}, "positive"), + ({"anhar_alpha_min": -12.0}, "integer below -2"), + ({"cal_anhar_fcs": True, "num_disp_anhar": -1}, "non-negative integer"), + ( + {"cal_anhar_fcs": True, "anhar_fit_methods": ["cocktail", "cocktail"]}, + "without repeats", + ), + ({"anhar_alpha_min": -1}, "integer below -2"), + ( + { + "cal_anhar_fcs": True, + "anhar_max_order": 4, + "fcs_cutoff_radius": [-1, 12], + }, + "positive", + ), + ], +) +def test_pheasy_anharmonic_settings_are_checked(settings, match): + with pytest.raises(ValueError, match=match): + PhononMaker(**settings) + + +def test_pheasy_maker_passes_anharmonic_settings(si_structure: Structure): + """Both pheasy jobs must get the same settings, or the dataset is split wrong.""" + # every value differs from its default, so a dropped keyword is caught + settings = { + "cal_anhar_fcs": True, + "num_displaced_supercells": 5, + "num_disp_anhar": 30, + "displacement_anhar": 0.05, + "anhar_max_order": 4, + "anhar_fit_methods": ["cocktail", "one-shot"], + "anhar_alpha_min": -10, + "fcs_cutoff_radius": [-1, 11, 9], + "displacement": 0.02, + "symprec": 1e-4, + "sym_reduce": False, + } + flow = PhononMaker(**settings).make(structure=si_structure) + jobs = {job.name: job for job in flow.jobs} + displacements = jobs["generate_phonon_displacements"].function_kwargs + fit = jobs["generate_frequencies_eigenvectors"].function_kwargs + for key in ( + "cal_anhar_fcs", + "num_displaced_supercells", + "anhar_max_order", + "anhar_fit_methods", + "fcs_cutoff_radius", + "displacement", + "symprec", + "sym_reduce", + ): + assert displacements[key] == fit[key] == settings[key] + assert displacements["num_disp_anhar"] == 30 + assert displacements["displacement_anhar"] == 0.05 + assert fit["anhar_alpha_min"] == -10 diff --git a/tutorials/pheasy_workflow.ipynb b/tutorials/pheasy_workflow.ipynb index 1a27755a5b..6f6742b043 100644 --- a/tutorials/pheasy_workflow.ipynb +++ b/tutorials/pheasy_workflow.ipynb @@ -46,7 +46,7 @@ "\n", "The workflow has the same basic structure as the phonopy-based phonon workflow and uses [Phonopy](https://doi.org/10.7566/JPSJ.92.012001) to post-process the force constants into band structures, densities of states, and thermodynamic properties.\n", "\n", - "To run this tutorial, the `pheasy` extra must be installed: `pip install 'atomate2[phonons,pheasy]'`. Anharmonic force-constant extraction (`cal_anhar_fcs=True`) additionally requires the `alamode` extra — see the atomate2 VASP documentation for installation hints." + "To run this tutorial, the `pheasy` extra must be installed: `pip install 'atomate2[pheasy]'`. It includes phonopy and ALM. ALM is compiled from source. See the atomate2 VASP documentation if that build fails." ] }, { @@ -83,7 +83,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Then we use the pheasy `PhononMaker` to generate a `Flow`. The maker mirrors the phonopy-based `PhononMaker`; the most important additional switches are `cal_anhar_fcs` (extract anharmonic force constants up to fourth order and renormalize the phonon energies) and `fcs_cutoff_radius` (cutoff radii for the second-, third- and fourth-order force constants). Here we only extract harmonic force constants for silicon. As always, a tight structural relaxation is performed first — make sure it is converged very accurately for production runs." + "Then we use the pheasy `PhononMaker` to generate a `Flow`. The maker mirrors the phonopy-based `PhononMaker`; the most important additional switches are `cal_anhar_fcs` (extract anharmonic force constants, third order by default, set by `anhar_max_order`) and `fcs_cutoff_radius` (cutoff radii in Bohr for the third- and fourth-order force constants, the first entry is not used). Here we only extract harmonic force constants for silicon. As always, a tight structural relaxation is performed first — make sure it is converged very accurately for production runs." ] }, { @@ -117,7 +117,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The flow relaxes the bulk structure, generates a set of randomly displaced supercells, computes their forces with VASP, and then fits the force constants with pheasy. We can visualize the flow first." + "The flow relaxes the bulk structure, generates the displaced supercells (a single finite displacement for silicon), computes their forces with VASP, and then fits the force constants with pheasy. We can visualize the flow first." ] }, { @@ -209,9 +209,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Anharmonic force constants and phonon renormalization\n", + "## Anharmonic force constants\n", "\n", - "To go beyond the harmonic approximation, set `cal_anhar_fcs=True`. Pheasy will then extract third- and fourth-order force constants with LASSO regression (the cutoff radii are controlled by `fcs_cutoff_radius`). With `renorm_phonon=True`, the phonon energies are additionally renormalized at the temperatures given by `renorm_temp`. The lattice thermal conductivity can additionally be computed with `cal_ther_cond=True` (requires phono3py). Both options need more displaced supercells than the harmonic run shown here — pheasy will tell you how many it generated.\n", + "To go beyond the harmonic approximation, set `cal_anhar_fcs=True`. Pheasy will then extract third-order force constants with LASSO regression, and fourth-order ones too with `anhar_max_order=4` (the cutoff radii are controlled by `fcs_cutoff_radius`). The number of anharmonic displaced supercells is set from the number of free force constants, with at least 20. Above 600 the job stops unless `num_disp_anhar` is set.\n", "\n", "The same workflow is also available for forcefields via `from atomate2.forcefields.flows.pheasy import PhononMaker`, which is a cheap way to try it out without DFT." ] From db70c81c014bb7fd9b25808ffad1931c45bf1918 Mon Sep 17 00:00:00 2001 From: Jiongzhi ZHENG Date: Tue, 29 Sep 2026 20:17:01 +0800 Subject: [PATCH 02/28] Update pheasy.py --- src/atomate2/common/jobs/pheasy.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/atomate2/common/jobs/pheasy.py b/src/atomate2/common/jobs/pheasy.py index eb2648b817..1bedf0183a 100644 --- a/src/atomate2/common/jobs/pheasy.py +++ b/src/atomate2/common/jobs/pheasy.py @@ -1,4 +1,5 @@ -"""Jobs for running phonon calculations with phonopy and pheasy.""" +"""Jobs for running phonon，Higher-order FCs calculations with phonopy and pheasy. +and lattice thermal conductivity calculations using SHENGBTE and FOURPHONON.""" from __future__ import annotations From e0e191b9b16ab55a485db7a3ee78faf9a8e6f351 Mon Sep 17 00:00:00 2001 From: Jiongzhi ZHENG Date: Tue, 29 Sep 2026 20:24:40 +0800 Subject: [PATCH 03/28] Update pheasy.py --- src/atomate2/common/jobs/pheasy.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/atomate2/common/jobs/pheasy.py b/src/atomate2/common/jobs/pheasy.py index 1bedf0183a..3bca46db70 100644 --- a/src/atomate2/common/jobs/pheasy.py +++ b/src/atomate2/common/jobs/pheasy.py @@ -840,6 +840,12 @@ def generate_frequencies_eigenvectors( work_dir=work_dir, random_seed=random_seed, ) + + # End of anharmonic fit section + # and begin to calculate the anharmonic phonon priopoerties + # (e.g., phonon lifetimes, linewidths, etc.) using SHENGBTE and FOURPHONON. + # maybe we also need to do the phonon renormalization here. + # TODO: Implement the calculation of anharmonic phonon properties and phonon renormalization here. if fc_file.exists(): # Read the force constants from the output file of pheasy code From 537609c4f15ae035a49cc8807bb65a85e5ce6729 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Tue, 29 Sep 2026 22:27:02 -0700 Subject: [PATCH 04/28] Fix lint in pheasy job comments --- src/atomate2/common/jobs/pheasy.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/atomate2/common/jobs/pheasy.py b/src/atomate2/common/jobs/pheasy.py index 3bca46db70..578ab1e492 100644 --- a/src/atomate2/common/jobs/pheasy.py +++ b/src/atomate2/common/jobs/pheasy.py @@ -1,5 +1,8 @@ -"""Jobs for running phonon，Higher-order FCs calculations with phonopy and pheasy. -and lattice thermal conductivity calculations using SHENGBTE and FOURPHONON.""" +"""Jobs for running phonon and higher-order FC calculations with phonopy and pheasy. + +Lattice thermal conductivity calculations using ShengBTE and FourPhonon are +planned. +""" from __future__ import annotations @@ -840,12 +843,13 @@ def generate_frequencies_eigenvectors( work_dir=work_dir, random_seed=random_seed, ) - + # End of anharmonic fit section - # and begin to calculate the anharmonic phonon priopoerties - # (e.g., phonon lifetimes, linewidths, etc.) using SHENGBTE and FOURPHONON. - # maybe we also need to do the phonon renormalization here. - # TODO: Implement the calculation of anharmonic phonon properties and phonon renormalization here. + # and begin to calculate the anharmonic phonon properties + # (e.g., phonon lifetimes, linewidths, etc.) using ShengBTE and FourPhonon. + # maybe we also need to do the phonon renormalization here. + # TODO: Implement the calculation of anharmonic phonon properties and + # phonon renormalization here. if fc_file.exists(): # Read the force constants from the output file of pheasy code From ac4d511342ba8e2e7cbd10619c50c341f1e2821b Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Wed, 30 Sep 2026 09:09:03 -0700 Subject: [PATCH 05/28] Run the pheasy short-cutoff refit in its own folder with the matrix files --- docs/user/codes/vasp.md | 2 ++ src/atomate2/common/jobs/pheasy.py | 40 +++++++++++++++++++----------- tests/common/jobs/test_pheasy.py | 40 ++++++++++++++++++++++++++++++ 3 files changed, 67 insertions(+), 15 deletions(-) diff --git a/docs/user/codes/vasp.md b/docs/user/codes/vasp.md index ee9e9ba3a7..7edc691cfc 100644 --- a/docs/user/codes/vasp.md +++ b/docs/user/codes/vasp.md @@ -350,6 +350,8 @@ 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`. This workflow was used to build the Materials Project's Harmonic Phonon Database, described in [this preprint](https://chemrxiv.org/doi/full/10.26434/chemrxiv.15004632/v1). +The pheasy fits can fail for a cell that is not in a standard setting, for example the primitive cell of MgO (mp-1265) as the Materials Project serves it. +For such a structure, set `use_symmetrized_structure="primitive"`, as was done for the database. By default, this workflow does not compute anharmonic force constants, but can be extended to using the `cal_anhar_fcs` kwarg. ALM, from the ALAMODE package, counts the free force constants used to size the random displacement sets. diff --git a/src/atomate2/common/jobs/pheasy.py b/src/atomate2/common/jobs/pheasy.py index 578ab1e492..da3aaac57e 100644 --- a/src/atomate2/common/jobs/pheasy.py +++ b/src/atomate2/common/jobs/pheasy.py @@ -72,6 +72,7 @@ "anharmonic_force_matrix": "force_matrix_anhar.npy", "website": "phonon_website.json", "one_shot_dir": "one_shot", + "refit_dir": "short_cutoff_refit", "anharmonic_fit_log": "pheasy_anharmonic_fit.log", } @@ -890,8 +891,7 @@ def generate_frequencies_eigenvectors( # If imaginary modes are present, we first use the hiphive code to enforce # some symmetry constraints to eliminate the imaginary modes (generally work # for small imaginary modes near Gamma point). If the imaginary modes are - # still present, a pheasy refit with a shorter cutoff (10 A) follows, but - # its result is currently not used (see the NOTE below). + # still present, a pheasy refit with a shorter cutoff (10 A) follows. if imaginary_modes: # Define a cluster space using the largest cutoff you can @@ -937,12 +937,19 @@ def generate_frequencies_eigenvectors( ) # Using a shorter cutoff (10 A) to generate the force constants to - # eliminate the imaginary modes near Gamma point in pheasy code. - # NOTE: the force constants are read back below from - # FORCE_CONSTANTS_short_cutoff, the file the hiPhive step above wrote, so - # the result of this pheasy refit is not used. If the refit succeeds, it - # overwrites FORCE_CONSTANTS in the job folder. + # eliminate the imaginary modes near Gamma point in pheasy code. The refit + # runs in its own folder, so FORCE_CONSTANTS of the harmonic fit, which the + # cocktail fit kept fixed, stays unchanged. if imaginary_modes: + refit_dir = Path(_DEFAULT_FILE_PATHS["refit_dir"]).resolve() + refit_dir.mkdir(exist_ok=True) + for filename in ( + "POSCAR", + _DEFAULT_FILE_PATHS["harmonic_displacements"], + _DEFAULT_FILE_PATHS["harmonic_force_matrix"], + ): + shutil.copy(filename, refit_dir / filename) + pheasy_cmd_11 = ( f"pheasy --dim {int(supercell_matrix[0][0])} " f"{int(supercell_matrix[1][1])} " @@ -963,7 +970,8 @@ def generate_frequencies_eigenvectors( f"{int(supercell_matrix[1][1])} " f"{int(supercell_matrix[2][2])} -w 2 -d --symprec " f"{float(symprec)} --c2 10.0 " - f"--ndata {int(num_har)} --disp_file" + f"--ndata {int(num_har)} --disp_file " + f"--disp_matrix_file {_DEFAULT_FILE_PATHS['harmonic_displacements']}" ) phonon.generate_displacements(distance=displacement) @@ -974,7 +982,8 @@ def generate_frequencies_eigenvectors( 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 --rasr BHH --ndata {int(num_har)} " + f"--force_matrix_file {_DEFAULT_FILE_PATHS['harmonic_force_matrix']}" ) else: @@ -983,15 +992,16 @@ def generate_frequencies_eigenvectors( f"{int(supercell_matrix[1][1])} " f"{int(supercell_matrix[2][2])} -f --full_ifc " f"--c2 10.0 -w 2 --symprec {float(symprec)} " - f"--rasr BHH --ndata {int(num_har)}" + f"--rasr BHH --ndata {int(num_har)} " + f"--force_matrix_file {_DEFAULT_FILE_PATHS['harmonic_force_matrix']}" ) - subprocess.call(shlex.split(pheasy_cmd_11)) - subprocess.call(shlex.split(pheasy_cmd_12)) - subprocess.call(shlex.split(pheasy_cmd_13)) - subprocess.call(shlex.split(pheasy_cmd_14)) + for cmd in (pheasy_cmd_11, pheasy_cmd_12, pheasy_cmd_13, pheasy_cmd_14): + subprocess.run(shlex.split(cmd), cwd=refit_dir, check=True) - force_constants = parse_FORCE_CONSTANTS(filename=new_fc_file) + force_constants = parse_FORCE_CONSTANTS( + filename=refit_dir / _DEFAULT_FILE_PATHS["force_constants"] + ) phonon.force_constants = force_constants phonon.symmetrize_force_constants() diff --git a/tests/common/jobs/test_pheasy.py b/tests/common/jobs/test_pheasy.py index 3bf0bb7067..16d749e82b 100644 --- a/tests/common/jobs/test_pheasy.py +++ b/tests/common/jobs/test_pheasy.py @@ -27,6 +27,7 @@ from atomate2.common.jobs.phonons import ( _generate_phonon_object, _get_num_irreducible_fcs, + _run_band_structure_and_plot, ) # fcs_cutoff_radius in Bohr. 8 Bohr (4.2 A) covers the first two neighbour @@ -132,6 +133,45 @@ def test_harmonic_and_anharmonic_split(tmp_dir, num_displaced_supercells): assert not np.allclose(harmonic[0] / 0.01, anharmonic[0] / 0.03) +def test_imaginary_mode_refit(tmp_dir, monkeypatch): + """The short-cutoff refit reads the matrix files and runs in its own folder.""" + + def _report_imaginary_modes(*args, **kwargs): + bs_symm_line, _ = _run_band_structure_and_plot(*args, **kwargs) + return bs_symm_line, True + + monkeypatch.setattr( + pheasy_jobs, "_run_band_structure_and_plot", _report_imaginary_modes + ) + structure = _cu_structure() + kwargs = { + **COMMON_KWARGS, + "supercell_matrix": [[2, 0, 0], [0, 2, 0], [0, 0, 2]], + "cal_anhar_fcs": False, + "fcs_cutoff_radius": FCS_CUTOFF_RADIUS, + } + + job = generate_phonon_displacements( + structure=structure, + num_displaced_supercells=0, + displacement_anhar=0.03, + num_disp_anhar=0, + **kwargs, + ) + responses = run_locally(job, create_folders=True, ensure_success=True) + job = generate_frequencies_eigenvectors( + structure=structure, + displacement_data=_emt_displacement_data(responses[job.uuid][1].output), + **FIT_KWARGS, + **kwargs, + ) + run_locally(job, create_folders=True, ensure_success=True) + + (fit_dir,) = (path.parent for path in Path.cwd().glob("job_*/disp_matrix.npy")) + refit = fit_dir / "short_cutoff_refit" / "FORCE_CONSTANTS" + assert parse_FORCE_CONSTANTS(str(refit)).shape == (32, 32, 3, 3) + + def test_get_num_anharmonic_supercells(monkeypatch): phonon = _generate_phonon_object(_cu_structure(), **COMMON_KWARGS) kwargs = { From c48688747fef88974fe9bc706066e2486e0d2d50 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Wed, 30 Sep 2026 15:46:12 -0700 Subject: [PATCH 06/28] Make pheasy read the supercell of the force data --- pyproject.toml | 3 +- src/atomate2/common/jobs/pheasy.py | 32 +++++++++++++-------- tests/common/jobs/test_pheasy.py | 46 ++++++++++++++++++++++++++++++ 3 files changed, 68 insertions(+), 13 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index affab73559..7d9110d588 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -53,7 +53,8 @@ mp = ["mp-api>=0.37.5"] # ("Force constants shape disagrees with crystal structure setting"). phonons = ["phonopy>=2.43,<5", "seekpath>=2.0.0"] # the pheasy fork adds the --disp_matrix_file and --force_matrix_file options -pheasy = ["atomate2[phonons,alamode]", "hiphive==1.5", "numpy<=2.2", "pheasy @ git+https://gitlab.com/hpsahasrabuddhe/pheasy.git@e67a777ecfde5ab368d71167837b099df75cfd4c"] +# and fixes the --scell option +pheasy = ["atomate2[phonons,alamode]", "hiphive==1.5", "numpy<=2.2", "pheasy @ git+https://gitlab.com/hpsahasrabuddhe/pheasy.git@1728f16642d5f2763f27695cd29ffada0f755d7f"] alamode = ["alm @ git+https://github.com/ttadano/ALM.git@f1d668fdee66e7e7218a04c88daf19d0e14fce0c#subdirectory=python"] hiphive = ["hiphive==1.5", "trainstation>=1.0", "atomate2[phonons,alamode]"] lobster = ["ijson>=3.2.2", "lobsterpy>=0.6.0"] diff --git a/src/atomate2/common/jobs/pheasy.py b/src/atomate2/common/jobs/pheasy.py index da3aaac57e..f882e3580e 100644 --- a/src/atomate2/common/jobs/pheasy.py +++ b/src/atomate2/common/jobs/pheasy.py @@ -314,7 +314,10 @@ def _run_anharmonic_fit( from run to run on the same forces. """ dim = " ".join(str(int(supercell_matrix[i][i])) for i in range(3)) - base = f"pheasy --dim {dim} -w {anhar_max_order} --symprec {float(symprec)}" + base = ( + f"pheasy --scell SPOSCAR --dim {dim} -w {anhar_max_order} " + f"--symprec {float(symprec)}" + ) cutoffs = f"--c3 {float(fcs_cutoff_radius[1] / ANGSTROM_TO_BOHR)}" if anhar_max_order == 4: cutoffs += f" --c4 {float(fcs_cutoff_radius[2] / ANGSTROM_TO_BOHR)}" @@ -631,7 +634,10 @@ def generate_frequencies_eigenvectors( verbose=False, ) - # Write the POSCAR and SPOSCAR files for the input of pheasy code + # Write the POSCAR and SPOSCAR files for the input of pheasy code. pheasy + # reads SPOSCAR (--scell SPOSCAR), so its atoms are in the order of the + # displacement and force matrices. The supercell pheasy builds itself from + # POSCAR can put an atom on a cell face one lattice vector away. supercell = phonon._supercell # noqa: SLF001 write_vasp("POSCAR", get_phonopy_structure(structure)) write_vasp("SPOSCAR", supercell) @@ -738,7 +744,7 @@ def generate_frequencies_eigenvectors( # 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 --dim {int(supercell_matrix[0][0])} " + 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" @@ -747,7 +753,7 @@ def generate_frequencies_eigenvectors( # Create the null space to further reduce the free parameters for # specific force constants and make them physically correct. pheasy_cmd_2 = ( - f"pheasy --dim {int(supercell_matrix[0][0])} " + 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" @@ -756,7 +762,7 @@ def generate_frequencies_eigenvectors( # Generate the Compressive Sensing matrix,i.e., displacement matrix # for the input of machine leaning method.i.e., LASSO, pheasy_cmd_3 = ( - f"pheasy --dim {int(supercell_matrix[0][0])} " + 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)} " @@ -774,7 +780,7 @@ def generate_frequencies_eigenvectors( # constraint, i.e., tag: --rasr BHH, is enforced during the # fitting process. pheasy_cmd_4 = ( - f"pheasy --dim {int(supercell_matrix[0][0])} " + 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)} " @@ -786,7 +792,7 @@ def generate_frequencies_eigenvectors( else: # Calculate the force constants using the least-squred method pheasy_cmd_4 = ( - f"pheasy --dim {int(supercell_matrix[0][0])} " + 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)} " @@ -829,6 +835,7 @@ def generate_frequencies_eigenvectors( work_dir.mkdir(exist_ok=True) for filename in ( "POSCAR", + "SPOSCAR", _DEFAULT_FILE_PATHS["anharmonic_displacements"], _DEFAULT_FILE_PATHS["anharmonic_force_matrix"], ): @@ -945,13 +952,14 @@ def generate_frequencies_eigenvectors( refit_dir.mkdir(exist_ok=True) for filename in ( "POSCAR", + "SPOSCAR", _DEFAULT_FILE_PATHS["harmonic_displacements"], _DEFAULT_FILE_PATHS["harmonic_force_matrix"], ): shutil.copy(filename, refit_dir / filename) pheasy_cmd_11 = ( - f"pheasy --dim {int(supercell_matrix[0][0])} " + f"pheasy --scell SPOSCAR --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)} " @@ -959,14 +967,14 @@ def generate_frequencies_eigenvectors( ) pheasy_cmd_12 = ( - f"pheasy --dim {int(supercell_matrix[0][0])} " + 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)} --c2 10.0 -w 2" ) pheasy_cmd_13 = ( - f"pheasy --dim {int(supercell_matrix[0][0])} " + 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 --symprec " f"{float(symprec)} --c2 10.0 " @@ -978,7 +986,7 @@ def generate_frequencies_eigenvectors( if len(phonon.displacements) > 3: pheasy_cmd_14 = ( - f"pheasy --dim {int(supercell_matrix[0][0])} " + f"pheasy --scell SPOSCAR --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)} " @@ -988,7 +996,7 @@ def generate_frequencies_eigenvectors( else: pheasy_cmd_14 = ( - f"pheasy --dim {int(supercell_matrix[0][0])} " + 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"--c2 10.0 -w 2 --symprec {float(symprec)} " diff --git a/tests/common/jobs/test_pheasy.py b/tests/common/jobs/test_pheasy.py index 16d749e82b..6bb45d0b0f 100644 --- a/tests/common/jobs/test_pheasy.py +++ b/tests/common/jobs/test_pheasy.py @@ -172,6 +172,52 @@ def _report_imaginary_modes(*args, **kwargs): assert parse_FORCE_CONSTANTS(str(refit)).shape == (32, 32, 3, 3) +def test_atom_on_cell_face(tmp_dir): + """An atom on a cell face must not spoil the harmonic fit. + + pheasy used to build its own supercell from POSCAR. Reading POSCAR can move + an atom on a cell face by one lattice vector, and the atoms of that supercell + then no longer matched the order of the force data. + """ + # Cu on the diamond sites of an fcc primitive cell, with the second atom + # just below x = 1 instead of at x = 0, as in a relaxed Si cell. For this + # lattice, reading POSCAR moves that atom by one lattice vector. + structure = Structure( + Lattice([[0, 2.725, 2.725], [2.725, 0, 2.725], [2.725, 2.725, 0]]), + ["Cu", "Cu"], + [[0.25, 0.25, 0.25], [1 - 1e-16, 0, 0]], + ) + kwargs = { + **COMMON_KWARGS, + "supercell_matrix": [[4, 0, 0], [0, 4, 0], [0, 0, 4]], + "cal_anhar_fcs": False, + "fcs_cutoff_radius": FCS_CUTOFF_RADIUS, + } + + job = generate_phonon_displacements( + structure=structure, + num_displaced_supercells=0, + displacement_anhar=0.03, + num_disp_anhar=0, + **kwargs, + ) + responses = run_locally(job, create_folders=True, ensure_success=True) + job = generate_frequencies_eigenvectors( + structure=structure, + displacement_data=_emt_displacement_data(responses[job.uuid][1].output), + **FIT_KWARGS, + **kwargs, + ) + run_locally(job, create_folders=True, ensure_success=True) + + (fit_dir,) = (path.parent for path in Path.cwd().glob("job_*/disp_matrix.npy")) + fc = parse_FORCE_CONSTANTS(str(fit_dir / "FORCE_CONSTANTS")) + disps = np.load(fit_dir / "disp_matrix.npy") + forces = np.load(fit_dir / "force_matrix.npy") + residual = forces + np.einsum("ijab,njb->nia", fc, disps) + assert np.linalg.norm(residual) < 0.05 * np.linalg.norm(forces) + + def test_get_num_anharmonic_supercells(monkeypatch): phonon = _generate_phonon_object(_cu_structure(), **COMMON_KWARGS) kwargs = { From c8734b5b354f887179d4830827b7a5865074933b Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Thu, 24 Sep 2026 08:49:49 -0700 Subject: [PATCH 07/28] Add thermal expansion workflow --- docs/user/codes/vasp.md | 54 ++++ pyproject.toml | 3 +- src/atomate2/common/flows/cte.py | 157 +++++++++++ src/atomate2/common/jobs/cte.py | 102 +++++++ src/atomate2/common/schemas/cte.py | 368 ++++++++++++++++++++++++++ src/atomate2/forcefields/flows/cte.py | 141 ++++++++++ src/atomate2/vasp/flows/cte.py | 85 ++++++ tests/common/jobs/test_cte.py | 285 ++++++++++++++++++++ tests/vasp/flows/test_cte.py | 71 +++++ 9 files changed, 1265 insertions(+), 1 deletion(-) create mode 100644 src/atomate2/common/flows/cte.py create mode 100644 src/atomate2/common/jobs/cte.py create mode 100644 src/atomate2/common/schemas/cte.py create mode 100644 src/atomate2/forcefields/flows/cte.py create mode 100644 src/atomate2/vasp/flows/cte.py create mode 100644 tests/common/jobs/test_cte.py create mode 100644 tests/vasp/flows/test_cte.py diff --git a/docs/user/codes/vasp.md b/docs/user/codes/vasp.md index 7edc691cfc..74f1a90fdd 100644 --- a/docs/user/codes/vasp.md +++ b/docs/user/codes/vasp.md @@ -346,6 +346,7 @@ phonon_flow = PhononMaker(min_length=15.0, store_force_constants=False).make( 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. +It also installs phono3py for the thermal expansion workflow. 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`. @@ -431,6 +432,59 @@ gruneisen_flow = GruneisenMaker( ).make(structure=structure) ``` +### Thermal expansion workflow + +`CTEMaker` calculates the thermal expansion tensor from third-order force constants, with the help of [Pheasy](https://doi.org/10.48550/arXiv.2508.01020) and [phono3py](https://doi.org/10.1088/1361-648X/acd831). +It needs the `pheasy` extra, see the Pheasy section above. + +First, a tight structural relaxation is performed. +The relaxed structure is then passed to the pheasy phonon workflow and to the elastic constant workflow. +The two do not depend on each other, so a workflow manager can run them at the same time. +Neither of them relaxes the structure again, so both use the same structure. +The phonon workflow fits the third-order force constants with LASSO, on randomly displaced supercells with 0.03 Å displacements. +The phonon maker builds its supercells with `min_length=12.0`. +By default, it uses the one-shot fit, which fits the second-order force constants together with them. +The cocktail fit keeps the second-order force constants of the harmonic fit. +The thermal expansion of each fit uses the second-order force constants of that fit. +phono3py then gives the mode Grüneisen tensors on a 12x12x12 q-point mesh. +The thermal expansion tensor follows from the mode heat capacities, the Grüneisen tensors and the elastic compliance. +It is computed for each fit in `anhar_fit_methods` of the phonon maker, from 0 K to 1000 K in steps of 10 K by default. +If a frequency on the mesh is below -0.1 THz, a warning is raised and the thermal expansion of that fit is not computed. + +The mode Grüneisen tensors come from the third-order force constants at the relaxed structure. +The phonon frequencies are not renormalized with temperature. +The workflow uses PBEsol by default. +The stress is more sensitive to ENCUT than the forces are, so check the ENCUT convergence of the elastic tensor for your material. +For metals, set `born_maker=None` in the phonon maker to skip the Born charge calculation. +The `compute_cte` job reads the force constant files from the folder of the pheasy fit, so it must run where that folder can be read. + +A thermal expansion workflow for VASP can be started as follows: +```python +from atomate2.vasp.flows.cte import CTEMaker +from pymatgen.core.structure import Structure + +structure = Structure( + lattice=[[0, 2.13, 2.13], [2.13, 0, 2.13], [2.13, 2.13, 0]], + species=["Mg", "O"], + coords=[[0, 0, 0], [0.5, 0.5, 0.5]], +) + +cte_flow = CTEMaker().make(structure=structure) +``` + +`update_user_incar_settings` changes the INCAR of every VASP job in the flow. +For example, this switches all of them to r2SCAN: +```python +from atomate2.vasp.powerups import update_user_incar_settings + +cte_flow = update_user_incar_settings(cte_flow, {"GGA": None, "METAGGA": "R2SCAN"}) +``` +The Born charge job uses DFPT with `LEPSILON = True`. +Check that your VASP version runs DFPT with r2SCAN before you use this setting. + +The same workflow runs with a force field via `from atomate2.forcefields.flows.cte import CTEMaker`. +`CTEMaker.from_force_field_name` sets one force field for the relaxation, the phonons and the elastic tensor. + ### Quasi-harmonic Workflow Uses the quasi-harmonic approximation with the help of Phonopy to compute thermodynamic properties. diff --git a/pyproject.toml b/pyproject.toml index 7d9110d588..fd6c0aa590 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -54,7 +54,8 @@ mp = ["mp-api>=0.37.5"] phonons = ["phonopy>=2.43,<5", "seekpath>=2.0.0"] # the pheasy fork adds the --disp_matrix_file and --force_matrix_file options # and fixes the --scell option -pheasy = ["atomate2[phonons,alamode]", "hiphive==1.5", "numpy<=2.2", "pheasy @ git+https://gitlab.com/hpsahasrabuddhe/pheasy.git@1728f16642d5f2763f27695cd29ffada0f755d7f"] +# phono3py computes the mode Grueneisen tensors in the thermal expansion workflow +pheasy = ["atomate2[phonons,alamode]", "hiphive==1.5", "numpy<=2.2", "pheasy @ git+https://gitlab.com/hpsahasrabuddhe/pheasy.git@1728f16642d5f2763f27695cd29ffada0f755d7f", "phono3py>=4.5"] alamode = ["alm @ git+https://github.com/ttadano/ALM.git@f1d668fdee66e7e7218a04c88daf19d0e14fce0c#subdirectory=python"] hiphive = ["hiphive==1.5", "trainstation>=1.0", "atomate2[phonons,alamode]"] lobster = ["ijson>=3.2.2", "lobsterpy>=0.6.0"] diff --git a/src/atomate2/common/flows/cte.py b/src/atomate2/common/flows/cte.py new file mode 100644 index 0000000000..d40fd25b26 --- /dev/null +++ b/src/atomate2/common/flows/cte.py @@ -0,0 +1,157 @@ +"""Flow for the thermal expansion from third-order force constants.""" + +from __future__ import annotations + +from abc import ABC, abstractmethod +from dataclasses import dataclass, field +from typing import TYPE_CHECKING + +from jobflow import Flow, Maker +from pymatgen.util.due import Doi, due + +from atomate2.common.jobs.cte import compute_cte + +if TYPE_CHECKING: + from pathlib import Path + + from pymatgen.core.structure import Structure + + from atomate2.common.flows.elastic import BaseElasticMaker + from atomate2.common.flows.pheasy import BasePhononMaker + from atomate2.forcefields.jobs import ForceFieldRelaxMaker + from atomate2.vasp.jobs.base import BaseVaspMaker + + +@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 BaseCTEMaker(Maker, ABC): + """ + Maker to calculate the thermal expansion from third-order force constants. + + A tight structural relaxation is performed first. The relaxed structure is + then passed to the pheasy phonon flow, which fits the second- and + third-order force constants, and to the elastic flow. The two flows do not + depend on each other, so a workflow manager can run them at the same time. + Neither flow may relax the structure again, so that both use the same + structure in the same frame. Finally, phono3py gives the mode + Grueneisen tensors on a q-point mesh, and the thermal expansion tensor + follows from the heat-capacity weighted Grueneisen tensors and the elastic + compliance. The mode Grueneisen tensors come from the third-order force + constants at the relaxed structure. The frequencies are not renormalized + with temperature. + + Parameters + ---------- + name: str + Name of the flows produced by this maker. + bulk_relax_maker: .ForceFieldRelaxMaker, .BaseVaspMaker, or None + A maker to perform a tight relaxation on the bulk. Set to None to skip + the relaxation. + phonon_maker: .BasePhononMaker + The pheasy phonon maker. It must have cal_anhar_fcs=True, + use_symmetrized_structure=None and bulk_relax_maker=None. Its + anhar_fit_methods set which force constants are used for the thermal + expansion. + elastic_maker: .BaseElasticMaker + Maker for the elastic tensor. It must have bulk_relax_maker=None. + temperatures: list[float] + Temperatures in K. + mesh: tuple[int, int, int] | float + q-point mesh for the mode Grueneisen tensors, or a q-point density used + as kppa in pymatgen's Kpoints.automatic_density for the unit cell. + tol_imaginary_modes: float + If a frequency on the mesh is below -tol_imaginary_modes in THz, a + warning is raised and the thermal expansion of that fit is not computed. + min_frequency: float + Modes below this frequency in THz are left out of the thermal expansion. + """ + + name: str = "cte" + bulk_relax_maker: ForceFieldRelaxMaker | BaseVaspMaker | None = None + phonon_maker: BasePhononMaker = None + elastic_maker: BaseElasticMaker = None + temperatures: list[float] = field(default_factory=lambda: list(range(0, 1001, 10))) + mesh: tuple[int, int, int] | float = (12, 12, 12) + tol_imaginary_modes: float = 0.1 + min_frequency: float = 1e-3 + + def __post_init__(self) -> None: + """Check that the phonon and elastic makers fit this workflow.""" + if not self.phonon_maker.cal_anhar_fcs: + raise ValueError( + "The phonon maker needs cal_anhar_fcs=True, since the thermal " + "expansion needs the third-order force constants." + ) + if self.phonon_maker.use_symmetrized_structure is not None: + raise ValueError( + "The phonon maker needs use_symmetrized_structure=None, so that the " + "phonon and elastic calculations use the same frame." + ) + for label, maker in ( + ("phonon", self.phonon_maker), + ("elastic", self.elastic_maker), + ): + if maker.bulk_relax_maker is not None: + raise ValueError( + f"The {label} maker needs bulk_relax_maker=None. Otherwise the " + "structure is relaxed again, and the phonon and elastic " + "calculations do not use the same structure." + ) + + def make(self, structure: Structure, prev_dir: str | Path | None = None) -> Flow: + """ + Make a flow to calculate the thermal expansion. + + Parameters + ---------- + structure: Structure + A pymatgen structure. Start with a structure that is nearly fully + optimized, as the relaxation settings are strict. + prev_dir: str or Path or None + A previous calculation directory to use for copying outputs. + """ + jobs = [] + equilibrium_stress = None + if self.bulk_relax_maker is not None: + bulk_kwargs = {} + if self.prev_calc_dir_argname is not None: + bulk_kwargs[self.prev_calc_dir_argname] = prev_dir + bulk = self.bulk_relax_maker.make(structure, **bulk_kwargs) + jobs.append(bulk) + structure = bulk.output.structure + prev_dir = bulk.output.dir_name + # as in the elastic flow when it runs its own relaxation + equilibrium_stress = bulk.output.output.stress + + phonon_flow = self.phonon_maker.make(structure, prev_dir=prev_dir) + elastic_flow = self.elastic_maker.make( + structure, prev_dir=prev_dir, equilibrium_stress=equilibrium_stress + ) + cte = compute_cte( + phonon_output=phonon_flow.output, + elastic_tensor=elastic_flow.output.elastic_tensor.raw, + elastic_structure=elastic_flow.output.structure, + anhar_fit_methods=self.phonon_maker.anhar_fit_methods, + temperatures=self.temperatures, + mesh=self.mesh, + tol_imaginary_modes=self.tol_imaginary_modes, + min_frequency=self.min_frequency, + symprec=self.phonon_maker.symprec, + ) + jobs += [phonon_flow, elastic_flow, cte] + return Flow(jobs, output=cte.output, name=self.name) + + @property + @abstractmethod + def prev_calc_dir_argname(self) -> str | None: + """Name of the argument that passes the previous calculation directory. + + It differs between codes, so each subclass sets it. + """ diff --git a/src/atomate2/common/jobs/cte.py b/src/atomate2/common/jobs/cte.py new file mode 100644 index 0000000000..982d398513 --- /dev/null +++ b/src/atomate2/common/jobs/cte.py @@ -0,0 +1,102 @@ +"""Jobs for thermal expansion from mode Grueneisen tensors.""" + +from __future__ import annotations + +from pathlib import Path +from typing import TYPE_CHECKING + +from jobflow import job + +from atomate2.common.jobs.gruneisen import PhononDoc, _get_taskdoc_run_dir +from atomate2.common.jobs.pheasy import _ANHARMONIC_FIT_METHODS, _DEFAULT_FILE_PATHS +from atomate2.common.schemas.cte import CTEDocument + +if TYPE_CHECKING: + from collections.abc import Sequence + + from emmet.core.math import MatrixVoigt + from pymatgen.core import Structure + + +@job(output_schema=CTEDocument) +def compute_cte( + phonon_output: PhononDoc, + elastic_tensor: MatrixVoigt, + elastic_structure: Structure, + anhar_fit_methods: Sequence[str] = ("one-shot",), + temperatures: Sequence[float] = tuple(range(0, 1001, 10)), + mesh: tuple[int, int, int] | float = (12, 12, 12), + tol_imaginary_modes: float = 0.1, + min_frequency: float = 1e-3, + symprec: float = 1e-5, +) -> CTEDocument: + """ + Compute the thermal expansion from the pheasy force constants. + + The second- and third-order force constants of each fit method are read + from the folder of the pheasy fit, so this job must run where that folder + can be read. The cocktail fit uses the second-order force constants of the + harmonic fit. The one-shot fit uses its own. + + Parameters + ---------- + phonon_output: PhononDoc + Output document of the pheasy phonon flow, run with cal_anhar_fcs=True. + elastic_tensor: MatrixVoigt + Elastic tensor in GPa and Voigt notation, as ElasticDocument's + elastic_tensor.raw. + elastic_structure: Structure + Structure of the elastic calculation. Its lattice must match the unit + cell of the phonon run, so that both tensors are in the same frame. + anhar_fit_methods: Sequence[str] + Fit methods whose force constants are used, "cocktail" and/or "one-shot". + temperatures: Sequence[float] + Temperatures in K, not negative. + mesh: tuple[int, int, int] | float + q-point mesh, or a q-point density used as kppa in pymatgen's + Kpoints.automatic_density for the unit cell. + tol_imaginary_modes: float + If a frequency on the mesh is below -tol_imaginary_modes in THz, a + warning is raised and the thermal expansion of that fit is not computed. + min_frequency: float + Modes below this frequency in THz are left out. + symprec: float + Symmetry precision passed to phono3py. + + Returns + ------- + CTEDocument + """ + unknown = set(anhar_fit_methods) - set(_ANHARMONIC_FIT_METHODS) + if unknown or not anhar_fit_methods: + raise ValueError( + f"anhar_fit_methods must be a non-empty subset of " + f"{_ANHARMONIC_FIT_METHODS}, not {list(anhar_fit_methods)}." + ) + job_dir_name = _get_taskdoc_run_dir(phonon_output) + if job_dir_name is None: + raise ValueError("The phonon output does not record the pheasy job folder.") + job_dir = Path(job_dir_name) + + # the files written by the pheasy fits + one_shot_dir = job_dir / _DEFAULT_FILE_PATHS["one_shot_dir"] + files = { + "cocktail": ( + job_dir / _DEFAULT_FILE_PATHS["force_constants"], + job_dir / "fc3.hdf5", + ), + "one-shot": (one_shot_dir / "fc2.hdf5", one_shot_dir / "fc3.hdf5"), + } + force_constant_files = {method: files[method] for method in anhar_fit_methods} + + return CTEDocument.from_force_constants( + phonopy_yaml=job_dir / _DEFAULT_FILE_PATHS["phonopy"], + force_constant_files=force_constant_files, + elastic_tensor=elastic_tensor, + structure=elastic_structure, + temperatures=temperatures, + mesh=mesh, + tol_imaginary_modes=tol_imaginary_modes, + min_frequency=min_frequency, + symprec=symprec, + ) diff --git a/src/atomate2/common/schemas/cte.py b/src/atomate2/common/schemas/cte.py new file mode 100644 index 0000000000..a071615926 --- /dev/null +++ b/src/atomate2/common/schemas/cte.py @@ -0,0 +1,368 @@ +"""Schemas for the thermal expansion workflow outputs.""" + +from __future__ import annotations + +import warnings +from pathlib import Path +from typing import TYPE_CHECKING + +import numpy as np +from emmet.core.math import Matrix3D, MatrixVoigt +from emmet.core.structure import StructureMetadata +from pydantic import BaseModel, Field +from pymatgen.core import Structure +from pymatgen.io.phonopy import get_pmg_structure +from pymatgen.io.vasp import Kpoints +from scipy.constants import Boltzmann, Planck + +if TYPE_CHECKING: + from collections.abc import Mapping, Sequence + + from phonopy import Phonopy + from typing_extensions import Self + +# Voigt order of the symmetric 3x3 tensor components, as in pymatgen +_VOIGT_INDICES = ((0, 0), (1, 1), (2, 2), (1, 2), (0, 2), (0, 1)) + + +def get_cte( + frequencies: np.ndarray, + gruneisen_tensors: np.ndarray, + weights: np.ndarray, + elastic_tensor: np.ndarray, + volume: float, + temperatures: Sequence[float], + min_frequency: float = 1e-3, +) -> tuple[np.ndarray, list[np.ndarray | None]]: + """ + Get the thermal expansion tensor from mode Grueneisen tensors. + + The thermal stress of each mode is its heat capacity times its Grueneisen + tensor. The strain that relaxes the summed thermal stress is the thermal + expansion, alpha = S sum(c * gamma) / (N_q * V). Here S is the elastic + compliance, c are the modal heat capacities, N_q is the sum of the q-point + weights and V is the cell volume. The mode Grueneisen tensors are + symmetrized first, since the strain is symmetric. Modes below + min_frequency, including the acoustic modes at Gamma, are left out. + + Parameters + ---------- + frequencies: np.ndarray + Phonon frequencies in THz, with shape (n_qpoints, n_bands). + gruneisen_tensors: np.ndarray + Mode Grueneisen tensors, with shape (n_qpoints, n_bands, 3, 3). + weights: np.ndarray + Weight of each q-point, with shape (n_qpoints,). + elastic_tensor: np.ndarray + Elastic tensor in GPa and Voigt notation, in the same Cartesian frame as + the Grueneisen tensors. + volume: float + Volume of the cell used for the phonons, in Angstrom^3. + temperatures: Sequence[float] + Temperatures in K. + min_frequency: float + Modes below this frequency in THz are left out. + + Returns + ------- + tuple[np.ndarray, list[np.ndarray | None]] + The thermal expansion tensors in 1/K, with shape (n_temperatures, 3, 3), + and the heat-capacity weighted mean Grueneisen tensor at each + temperature. The mean is None where the heat capacity is zero. + """ + frequencies = np.asarray(frequencies, dtype=float) + gruneisen_tensors = np.asarray(gruneisen_tensors, dtype=float) + weights = np.asarray(weights, dtype=float) + + # S in 1/Pa and V in m^3 + compliance = np.linalg.inv(np.asarray(elastic_tensor, dtype=float) * 1e9) + volume_m3 = volume * 1e-30 + + # the Grueneisen tensors of the left-out modes can be large or undefined, + # so zero them + kept = frequencies > min_frequency + gruneisen_tensors = np.where(kept[..., None, None], gruneisen_tensors, 0.0) + gruneisen_tensors = (gruneisen_tensors + np.swapaxes(gruneisen_tensors, -1, -2)) / 2 + gruneisen_voigt = np.stack( + [gruneisen_tensors[..., i, j] for i, j in _VOIGT_INDICES], axis=-1 + ) + energies = Planck * frequencies * 1e12 # J + + alphas, mean_gruneisen = [], [] + for temperature in temperatures: + if temperature <= 0: + heat_capacities = np.zeros_like(frequencies) + else: + x = np.where(kept, energies / (Boltzmann * temperature), 1.0) + # x^2 e^x / (e^x - 1)^2, written with e^-x so that it cannot overflow + heat_capacities = np.where( + kept, Boltzmann * x**2 * np.exp(-x) / np.expm1(-x) ** 2, 0.0 + ) # J/K + weighted = heat_capacities * weights[:, None] + thermal_stress = np.einsum("qb,qbv->v", weighted, gruneisen_voigt) + thermal_stress /= weights.sum() * volume_m3 # Pa/K + + # the compliance gives engineering shear strains, twice the tensor ones + alpha_voigt = compliance @ thermal_stress + alpha = np.empty((3, 3)) + for k, (i, j) in enumerate(_VOIGT_INDICES): + alpha[i, j] = alpha[j, i] = alpha_voigt[k] if k < 3 else alpha_voigt[k] / 2 + alphas.append(alpha) + + total_heat_capacity = weighted.sum() + if total_heat_capacity > 0: + mean = np.einsum("qb,qbij->ij", weighted, gruneisen_tensors) + mean_gruneisen.append(mean / total_heat_capacity) + else: + mean_gruneisen.append(None) + + return np.array(alphas), mean_gruneisen + + +def _expand_born_to_unitcell(phonon: Phonopy) -> np.ndarray: + """Map the Born charges of the phonopy primitive cell onto the unit cell atoms.""" + primitive = phonon.primitive + born = np.asarray(phonon.nac_params["born"]) + if len(primitive) == len(phonon.unitcell): + return born + primitive_indices = [ + primitive.p2p_map[primitive.s2p_map[s]] for s in phonon.supercell.u2s_map + ] + return born[primitive_indices] + + +class CTEResult(BaseModel): + """Thermal expansion from one set of second- and third-order force constants.""" + + fit_method: str = Field( + description='Anharmonic fit that gave the force constants, "cocktail" or ' + '"one-shot".' + ) + lowest_frequency: float = Field( + description="Lowest phonon frequency on the sampling mesh in THz. Imaginary " + "frequencies are negative." + ) + has_imaginary_modes: bool = Field( + description="Whether a frequency on the sampling mesh lies below " + "-tol_imaginary_modes. The thermal expansion is not computed in that case." + ) + thermal_expansion: list[Matrix3D] | None = Field( + None, + description="Thermal expansion tensor in 1/K at each temperature, in the " + "Cartesian frame of the structure.", + ) + volumetric_thermal_expansion: list[float] | None = Field( + None, + description="Volumetric thermal expansion in 1/K at each temperature, the " + "trace of the thermal expansion tensor.", + ) + average_gruneisen: list[Matrix3D | None] | None = Field( + None, + description="Mode Grueneisen tensor averaged with the mode heat capacities " + "at each temperature. None where the heat capacity is zero.", + ) + + +class CTEDocument(StructureMetadata): + """Thermal expansion from mode Grueneisen tensors and the elastic tensor.""" + + structure: Structure | None = Field( + None, description="Structure used for the phonon and elastic calculations." + ) + temperatures: list[float] | None = Field(None, description="Temperatures in K.") + mesh: tuple[int, int, int] | None = Field( + None, description="q-point mesh used for the mode Grueneisen tensors." + ) + elastic_tensor: MatrixVoigt | None = Field( + None, + description="Elastic tensor in GPa and Voigt notation, in the Cartesian " + "frame of the structure.", + ) + min_frequency: float | None = Field( + None, + description="Modes below this frequency in THz are left out of the thermal " + "expansion.", + ) + tol_imaginary_modes: float | None = Field( + None, + description="The thermal expansion of a fit is not computed if a frequency " + "on the mesh is below -tol_imaginary_modes in THz.", + ) + phonon_job_dir: str | None = Field( + None, description="Directory of the pheasy fit that wrote the force constants." + ) + results: list[CTEResult] | None = Field( + None, description="Thermal expansion for each anharmonic fit method." + ) + + @classmethod + def from_force_constants( + cls, + phonopy_yaml: str | Path, + force_constant_files: Mapping[str, tuple[str | Path, str | Path]], + elastic_tensor: MatrixVoigt, + structure: Structure, + temperatures: Sequence[float], + mesh: tuple[int, int, int] | float, + tol_imaginary_modes: float, + min_frequency: float, + symprec: float, + ) -> Self: + """ + Compute the thermal expansion from second- and third-order force constants. + + phono3py gives the frequencies and mode Grueneisen tensors on the q-point + mesh. Only time reversal symmetry is used to reduce the mesh. The + non-analytical term correction is applied to the dynamical matrix when + the phonopy yaml file holds the Born charges and the dielectric tensor. + Only the third-order force constants enter the strain derivative of the + dynamical matrix. The results of each fit are written to + gruneisen_.hdf5 in the current directory. + + Parameters + ---------- + phonopy_yaml: str or Path + phonopy.yaml of the pheasy fit, with the unit cell, the supercell + matrix and the Born charges. + force_constant_files: Mapping + For each fit method, "cocktail" or "one-shot", the files with the + second- and third-order force constants. + elastic_tensor: MatrixVoigt + Elastic tensor in GPa and Voigt notation. + structure: Structure + Structure of the elastic calculation. Its lattice must match the + unit cell in phonopy_yaml, so that both tensors are in the same frame. + temperatures: Sequence[float] + Temperatures in K, not negative. + mesh: tuple[int, int, int] or float + q-point mesh, or a q-point density used as kppa in pymatgen's + Kpoints.automatic_density for the unit cell. + tol_imaginary_modes: float + If a frequency on the mesh is below -tol_imaginary_modes in THz, a + warning is raised and the thermal expansion of that fit is not + computed. + min_frequency: float + Modes below this frequency in THz are left out. + symprec: float + Symmetry precision passed to phono3py. + + Returns + ------- + CTEDocument + """ + import h5py + import phonopy + from phono3py import Phono3py + from phono3py.file_IO import read_fc2_from_hdf5, read_fc3_from_hdf5 + from phono3py.phonon3.gruneisen import Gruneisen + from phonopy.file_IO import parse_FORCE_CONSTANTS + + if not force_constant_files: + raise ValueError("No force constant files were given.") + if min(temperatures) < 0: + raise ValueError("The temperatures must not be negative.") + + phonon = phonopy.load(phonopy_yaml, produce_fc=False, log_level=0) + if not np.allclose(phonon.unitcell.cell, structure.lattice.matrix, atol=1e-5): + raise ValueError( + "The lattice of the elastic calculation differs from the unit cell " + "of the phonon calculation, so the two tensors are not in the same " + "frame." + ) + if list(phonon.unitcell.symbols) != [site.specie.symbol for site in structure]: + raise ValueError( + "The atoms of the elastic calculation differ from the unit cell of " + "the phonon calculation." + ) + + # pheasy writes compact force constants for the unit cell in POSCAR, so + # the unit cell is also the primitive cell here. "P" keeps it as it is. + ph3 = Phono3py( + phonon.unitcell, + supercell_matrix=phonon.supercell_matrix, + primitive_matrix="P", + symprec=symprec, + ) + nac_params = None + if phonon.nac_params is not None: + nac_params = { + **phonon.nac_params, + "born": _expand_born_to_unitcell(phonon), + } + + if isinstance(mesh, int | float | np.number): + kpoints = Kpoints.automatic_density( + structure=get_pmg_structure(ph3.primitive), + kppa=float(mesh), + force_gamma=True, + ) + mesh_numbers = tuple(int(m) for m in kpoints.kpts[0]) + else: + mesh_numbers = tuple(int(m) for m in mesh) + + results = [] + for method, (fc2_file, fc3_file) in force_constant_files.items(): + if Path(fc2_file).suffix == ".hdf5": + fc2 = read_fc2_from_hdf5(fc2_file) + else: + fc2 = parse_FORCE_CONSTANTS(filename=fc2_file) + fc3 = read_fc3_from_hdf5(fc3_file) + + gruneisen = Gruneisen( + fc2, fc3, ph3.supercell, ph3.primitive, nac_params=nac_params + ) + gruneisen.set_sampling_mesh(mesh_numbers, primitive_symmetry=None) + gruneisen.run() + filename = f"gruneisen_{method.replace('-', '_')}" + gruneisen.write(filename=filename) + # the full third-order force constants can take several GB + del gruneisen, fc3 + with h5py.File(f"{filename}.hdf5") as file: + frequencies = file["frequency"][:] + gruneisen_tensors = file["gruneisen_tensor"][:] + weights = file["weight"][:] + + lowest = float(frequencies.min()) + has_imaginary_modes = lowest < -tol_imaginary_modes + result = { + "fit_method": method, + "lowest_frequency": lowest, + "has_imaginary_modes": has_imaginary_modes, + } + if has_imaginary_modes: + warnings.warn( + f"The {method} force constants give a frequency of " + f"{lowest:.3f} THz, below -{tol_imaginary_modes} THz. The " + "thermal expansion is not computed for this fit.", + stacklevel=2, + ) + else: + alphas, mean_gruneisen = get_cte( + frequencies, + gruneisen_tensors, + weights, + np.asarray(elastic_tensor), + ph3.primitive.volume, + temperatures, + min_frequency=min_frequency, + ) + result["thermal_expansion"] = alphas.tolist() + result["volumetric_thermal_expansion"] = np.trace( + alphas, axis1=1, axis2=2 + ).tolist() + result["average_gruneisen"] = [ + None if mean is None else mean.tolist() for mean in mean_gruneisen + ] + results.append(CTEResult(**result)) + + return cls.from_structure( + meta_structure=structure, + structure=structure, + temperatures=list(temperatures), + mesh=mesh_numbers, + elastic_tensor=np.asarray(elastic_tensor).tolist(), + min_frequency=min_frequency, + tol_imaginary_modes=tol_imaginary_modes, + phonon_job_dir=str(Path(phonopy_yaml).parent), + results=results, + ) diff --git a/src/atomate2/forcefields/flows/cte.py b/src/atomate2/forcefields/flows/cte.py new file mode 100644 index 0000000000..ed87a7fd28 --- /dev/null +++ b/src/atomate2/forcefields/flows/cte.py @@ -0,0 +1,141 @@ +"""Define the force field thermal expansion maker.""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import TYPE_CHECKING, Any + +from atomate2.common.flows.cte import BaseCTEMaker +from atomate2.forcefields.flows.elastic import ElasticMaker +from atomate2.forcefields.flows.pheasy import PhononMaker +from atomate2.forcefields.jobs import ForceFieldRelaxMaker, ForceFieldStaticMaker + +if TYPE_CHECKING: + from typing_extensions import Self + + from atomate2.forcefields import MLFF + +_DEFAULT_FORCE_FIELD = "MACE-MP-0" + + +def _get_makers( + force_field_name: str | MLFF | dict, calculator_kwargs: dict | None = None +) -> dict: + """Get the relaxation, phonon and elastic makers for one force field.""" + calculator: dict[str, Any] = { + "force_field_name": force_field_name, + "calculator_kwargs": dict(calculator_kwargs or {}), + } + # the relaxation settings are those of the force field elastic flow + return { + "bulk_relax_maker": ForceFieldRelaxMaker( + relax_cell=True, + relax_kwargs={"fmax": 0.00001}, + fix_symmetry=True, + **calculator, + ), + "phonon_maker": PhononMaker( + min_length=12.0, + bulk_relax_maker=None, + static_energy_maker=None, + phonon_displacement_maker=ForceFieldStaticMaker(**calculator), + cal_anhar_fcs=True, + displacement_anhar=0.03, + anhar_fit_methods=("one-shot",), + ), + "elastic_maker": ElasticMaker( + bulk_relax_maker=None, + elastic_relax_maker=ForceFieldRelaxMaker( + relax_cell=False, + relax_kwargs={"fmax": 0.00001}, + fix_symmetry=True, + **calculator, + ), + ), + } + + +@dataclass +class CTEMaker(BaseCTEMaker): + """ + Maker to calculate the thermal expansion with a force field and pheasy. + + A tight relaxation is performed first. The relaxed structure is then passed + to the pheasy phonon flow, which fits the second- and third-order force + constants, and to the elastic flow. Neither flow relaxes the structure + again. Finally, phono3py gives the mode Grueneisen tensors, and the thermal + expansion tensor follows from them and the elastic tensor. The frequencies + are not renormalized with temperature. + + By default, all steps use MACE-MP-0. Use :obj:`from_force_field_name` to + run every step with another force field. The phonon flow fits the force + constants with the one-shot method from randomly displaced supercells with + 0.03 A displacements, and builds the supercells with min_length=12.0. + + Parameters + ---------- + name: str + Name of the flows produced by this maker. + bulk_relax_maker: .ForceFieldRelaxMaker or None + A maker to perform a tight relaxation on the bulk. Set to None to skip + the relaxation. + phonon_maker: .PhononMaker + The pheasy phonon maker. It must have cal_anhar_fcs=True, + use_symmetrized_structure=None and bulk_relax_maker=None. Its + anhar_fit_methods set which force constants are used for the thermal + expansion. + elastic_maker: .ElasticMaker + Maker for the elastic tensor. It must have bulk_relax_maker=None. + temperatures: list[float] + Temperatures in K. + mesh: tuple[int, int, int] | float + q-point mesh for the mode Grueneisen tensors, or a q-point density used + as kppa in pymatgen's Kpoints.automatic_density for the unit cell. + tol_imaginary_modes: float + If a frequency on the mesh is below -tol_imaginary_modes in THz, a + warning is raised and the thermal expansion of that fit is not computed. + min_frequency: float + Modes below this frequency in THz are left out of the thermal expansion. + """ + + name: str = "cte" + bulk_relax_maker: ForceFieldRelaxMaker | None = field( + default_factory=lambda: _get_makers(_DEFAULT_FORCE_FIELD)["bulk_relax_maker"] + ) + phonon_maker: PhononMaker = field( + default_factory=lambda: _get_makers(_DEFAULT_FORCE_FIELD)["phonon_maker"] + ) + elastic_maker: ElasticMaker = field( + default_factory=lambda: _get_makers(_DEFAULT_FORCE_FIELD)["elastic_maker"] + ) + + @classmethod + def from_force_field_name( + cls, + force_field_name: str | MLFF | dict, + calculator_kwargs: dict | None = None, + **kwargs, + ) -> Self: + """ + Create a thermal expansion maker that uses one force field for all steps. + + Parameters + ---------- + force_field_name: str or .MLFF or dict + The name of the force field. + calculator_kwargs: dict or None + Keyword arguments passed to the force field calculator. + **kwargs + Further keyword arguments passed to CTEMaker. A maker given here + replaces the one built for the force field. + + Returns + ------- + CTEMaker + """ + return cls(**{**_get_makers(force_field_name, calculator_kwargs), **kwargs}) + + @property + def prev_calc_dir_argname(self) -> None: + """Name of the argument that passes the previous calculation directory.""" + return diff --git a/src/atomate2/vasp/flows/cte.py b/src/atomate2/vasp/flows/cte.py new file mode 100644 index 0000000000..9f63c4820c --- /dev/null +++ b/src/atomate2/vasp/flows/cte.py @@ -0,0 +1,85 @@ +"""Define the VASP thermal expansion maker.""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import TYPE_CHECKING + +from atomate2.common.flows.cte import BaseCTEMaker +from atomate2.vasp.flows.core import DoubleRelaxMaker +from atomate2.vasp.flows.elastic import ElasticMaker +from atomate2.vasp.flows.pheasy import PhononMaker +from atomate2.vasp.jobs.core import TightRelaxMaker + +if TYPE_CHECKING: + from atomate2.vasp.jobs.base import BaseVaspMaker + + +@dataclass +class CTEMaker(BaseCTEMaker): + """ + Maker to calculate the thermal expansion with VASP, pheasy and phono3py. + + A tight double relaxation is performed first. The relaxed structure is then + passed to the pheasy phonon flow, which fits the second- and third-order + force constants, and to the elastic flow. Neither flow relaxes the + structure again. Finally, phono3py gives the mode Grueneisen tensors, and + the thermal expansion tensor follows from them and the elastic tensor. The + frequencies are not renormalized with temperature. + + By default, the phonon flow fits the force constants with the one-shot + method from randomly displaced supercells with 0.03 A displacements, and + builds the supercells with min_length=12.0. The phonon flow skips the + static energy calculation, which the thermal expansion does not need. The + elastic flow is the atomate2 elastic flow without its own + relaxation. The stress is more sensitive to ENCUT than the forces are, so + check the ENCUT convergence of the elastic tensor for your material. + + Parameters + ---------- + name: str + Name of the flows produced by this maker. + bulk_relax_maker: .BaseVaspMaker or None + A maker to perform a tight relaxation on the bulk. Set to None to skip + the relaxation. + phonon_maker: .PhononMaker + The pheasy phonon maker. It must have cal_anhar_fcs=True, + use_symmetrized_structure=None and bulk_relax_maker=None. Its + anhar_fit_methods set which force constants are used for the thermal + expansion. + elastic_maker: .ElasticMaker + Maker for the elastic tensor. It must have bulk_relax_maker=None. + temperatures: list[float] + Temperatures in K. + mesh: tuple[int, int, int] | float + q-point mesh for the mode Grueneisen tensors, or a q-point density used + as kppa in pymatgen's Kpoints.automatic_density for the unit cell. + tol_imaginary_modes: float + If a frequency on the mesh is below -tol_imaginary_modes in THz, a + warning is raised and the thermal expansion of that fit is not computed. + min_frequency: float + Modes below this frequency in THz are left out of the thermal expansion. + """ + + name: str = "cte" + bulk_relax_maker: BaseVaspMaker | None = field( + default_factory=lambda: DoubleRelaxMaker.from_relax_maker(TightRelaxMaker()) + ) + phonon_maker: PhononMaker = field( + default_factory=lambda: PhononMaker( + min_length=12.0, + bulk_relax_maker=None, + static_energy_maker=None, + cal_anhar_fcs=True, + displacement_anhar=0.03, + anhar_fit_methods=("one-shot",), + ) + ) + elastic_maker: ElasticMaker = field( + default_factory=lambda: ElasticMaker(bulk_relax_maker=None) + ) + + @property + def prev_calc_dir_argname(self) -> str: + """Name of the argument that passes the previous calculation directory.""" + return "prev_dir" diff --git a/tests/common/jobs/test_cte.py b/tests/common/jobs/test_cte.py new file mode 100644 index 0000000000..8be61f9215 --- /dev/null +++ b/tests/common/jobs/test_cte.py @@ -0,0 +1,285 @@ +"""Tests for the thermal expansion workflow. + +The end-to-end test runs the force field CTEMaker with ASE's EMT potential, so +it needs no DFT reference data. It lives here and not under tests/forcefields, +because pheasy and ALM are only installed in the test-non-ase CI job. +""" + +from pathlib import Path + +import numpy as np +import phono3py.phonon3.gruneisen as gruneisen_module +import phonopy +import pytest +from ase.build import bulk +from jobflow import run_locally +from phonopy import Phonopy +from phonopy.structure.atoms import PhonopyAtoms +from pymatgen.analysis.elasticity import ElasticTensor +from pymatgen.io.ase import AseAtomsAdaptor +from scipy.constants import Boltzmann, Planck +from scipy.spatial.transform import Rotation + +from atomate2.common.jobs.cte import compute_cte +from atomate2.common.schemas.cte import CTEDocument, _expand_born_to_unitcell, get_cte +from atomate2.forcefields.flows.cte import CTEMaker + +EMT = {"@module": "ase.calculators.emt", "@callable": "EMT"} + + +def _heat_capacity(frequency: float, temperature: float) -> float: + x = Planck * frequency * 1e12 / (Boltzmann * temperature) + return Boltzmann * x**2 * np.exp(x) / (np.exp(x) - 1) ** 2 + + +def test_get_cte_cubic(): + """For a cubic crystal with gamma = g * I, alpha = g * Cv / (3 * B * V).""" + c11, c12, c44 = 170.0, 120.0, 75.0 + elastic = np.zeros((6, 6)) + elastic[:3, :3] = c12 + np.fill_diagonal(elastic[:3, :3], c11) + elastic[3:, 3:] = np.eye(3) * c44 + bulk_modulus = (c11 + 2 * c12) / 3 * 1e9 # Pa + + frequencies = np.array([[2.0, 5.0, 7.0], [3.0, 4.0, 6.0]]) + weights = np.array([1, 3]) + gruneisen = 1.7 * np.broadcast_to(np.eye(3), (2, 3, 3, 3)) + volume, temperature = 45.0, 300.0 + + alpha, mean_gruneisen = get_cte( + frequencies, gruneisen, weights, elastic, volume, [temperature] + ) + + heat_capacity = ( + sum( + w * _heat_capacity(f, temperature) + for w, row in zip(weights, frequencies, strict=True) + for f in row + ) + / weights.sum() + ) + expected = 1.7 * heat_capacity / (3 * bulk_modulus * volume * 1e-30) + assert alpha[0] == pytest.approx(expected * np.eye(3), rel=1e-10, abs=1e-20) + assert mean_gruneisen[0] == pytest.approx(1.7 * np.eye(3)) + + +def test_get_cte_rotation(): + """Rotating the inputs must rotate alpha, which checks the shear terms. + + The mode Grueneisen tensors from phono3py are not symmetric, so the input + tensors here are not symmetric either. + """ + # hexagonal elastic tensor, with C66 = (C11 - C12) / 2 + c11, c12, c13, c33, c44 = 350.0, 120.0, 90.0, 400.0, 110.0 + elastic = np.array( + [ + [c11, c12, c13, 0, 0, 0], + [c12, c11, c13, 0, 0, 0], + [c13, c13, c33, 0, 0, 0], + [0, 0, 0, c44, 0, 0], + [0, 0, 0, 0, c44, 0], + [0, 0, 0, 0, 0, (c11 - c12) / 2], + ] + ) + rng = np.random.default_rng(7) + frequencies = rng.uniform(1.0, 10.0, size=(4, 6)) + gruneisen = rng.normal(1.0, 0.5, size=(4, 6, 3, 3)) + weights = np.ones(4) + temperatures = [100.0, 500.0] + + alpha, _ = get_cte(frequencies, gruneisen, weights, elastic, 30.0, temperatures) + rotation = Rotation.from_euler("zxz", [30, 50, 70], degrees=True).as_matrix() + rotated_elastic = ElasticTensor.from_voigt(elastic).rotate(rotation).voigt + rotated_gruneisen = np.einsum("ik,qbkl,jl->qbij", rotation, gruneisen, rotation) + rotated_alpha, _ = get_cte( + frequencies, rotated_gruneisen, weights, rotated_elastic, 30.0, temperatures + ) + + expected = np.einsum("ik,tkl,jl->tij", rotation, alpha, rotation) + assert np.abs(expected[:, 0, 1]).max() > 1e-7 # the shear terms are tested + assert rotated_alpha == pytest.approx(expected, rel=1e-8, abs=1e-15) + + +def test_get_cte_left_out_modes(): + """Modes below min_frequency are left out, and T = 0 gives zero.""" + elastic = np.diag([200.0, 200.0, 200.0, 80.0, 80.0, 80.0]) + frequencies = np.array([[0.0, 0.0, 0.0, 4.0], [-0.5, 3.0, 5.0, 6.0]]) + gruneisen = np.ones((2, 4, 3, 3)) + gruneisen[0, :3] = np.nan # NaN checks that the left-out modes are ignored + + alpha, mean_gruneisen = get_cte( + frequencies, gruneisen, np.ones(2), elastic, 40.0, [0.0, 300.0] + ) + assert np.all(alpha[0] == 0.0) + assert mean_gruneisen[0] is None + assert mean_gruneisen[1] == pytest.approx(np.ones((3, 3))) + + # only the modes at 4, 3, 5 and 6 THz count, over two q-points + heat_capacity = sum(_heat_capacity(f, 300.0) for f in (4.0, 3.0, 5.0, 6.0)) / 2 + thermal_stress = heat_capacity / (40.0 * 1e-30) + expected = np.full((3, 3), thermal_stress / 80e9 / 2) + np.fill_diagonal(expected, thermal_stress / 200e9) + assert alpha[1] == pytest.approx(expected, rel=1e-10) + + +def test_expand_born_to_unitcell(): + """Born charges of the primitive cell are mapped onto the conventional cell.""" + atoms = bulk("MgO", "rocksalt", a=4.2, cubic=True) + unitcell = PhonopyAtoms( + symbols=atoms.get_chemical_symbols(), + cell=atoms.cell[:], + scaled_positions=atoms.get_scaled_positions(), + ) + phonon = Phonopy(unitcell, supercell_matrix=np.eye(3), primitive_matrix="auto") + assert len(phonon.primitive) == 2 + born_primitive = {"Mg": 1.9, "O": -1.9} + phonon.nac_params = { + "born": [np.eye(3) * born_primitive[s] for s in phonon.primitive.symbols], + "dielectric": np.eye(3) * 3.0, + } + + born = _expand_born_to_unitcell(phonon) + assert born.shape == (8, 3, 3) + for symbol, charges in zip(unitcell.symbols, born, strict=True): + assert charges == pytest.approx(np.eye(3) * born_primitive[symbol]) + + +def test_cte_maker_emt(clean_dir, monkeypatch): + """Run the whole force field workflow with EMT forces on fcc Cu.""" + structure = AseAtomsAdaptor.get_structure(bulk("Cu", "fcc", a=3.61, cubic=True)) + maker = CTEMaker.from_force_field_name( + EMT, temperatures=[0, 100, 300], mesh=(8, 8, 8) + ) + # a 2x2x2 supercell of the cubic cell, 32 atoms, and an fc3 cutoff of 6 Bohr + # (3.2 A), which covers the nearest neighbours at 2.55 A + maker.phonon_maker.min_length = 7.0 + maker.phonon_maker.fcs_cutoff_radius = [-1, 6, 6] + maker.phonon_maker.anhar_fit_methods = ["cocktail", "one-shot"] + + flow = maker.make(structure) + # the uuid of the phonon flow output that compute_cte reads + phonon_uuid = flow.jobs[-1].function_kwargs["phonon_output"].uuid + responses = run_locally(flow, create_folders=True, ensure_success=True) + doc = responses[flow.output.uuid][1].output + assert isinstance(doc, CTEDocument) + assert doc.mesh == (8, 8, 8) + assert [result.fit_method for result in doc.results] == ["cocktail", "one-shot"] + + # EMT elastic constants of Cu in GPa + assert doc.elastic_tensor[0][0] == pytest.approx(172.6, rel=0.01) + assert doc.elastic_tensor[0][1] == pytest.approx(115.4, rel=0.01) + assert doc.elastic_tensor[3][3] == pytest.approx(89.9, rel=0.01) + + for result in doc.results: + assert not result.has_imaginary_modes + alpha = np.array(result.thermal_expansion) + assert np.all(alpha[0] == 0.0) + # cubic, so alpha is isotropic in every frame + assert alpha[2] == pytest.approx(alpha[2, 0, 0] * np.eye(3), abs=1e-12) + assert result.volumetric_thermal_expansion[2] == pytest.approx( + 3 * alpha[2, 0, 0] + ) + # linear thermal expansion at 300 K in 1/K, and the mean Grueneisen parameter. + # The one-shot fit has few supercells in this small cell, so its value is only + # compared with the cocktail value. + cocktail, one_shot = doc.results + assert cocktail.thermal_expansion[2][0][0] == pytest.approx(1.859e-5, rel=0.02) + assert np.trace(cocktail.average_gruneisen[2]) / 3 == pytest.approx(2.237, rel=0.02) + assert one_shot.thermal_expansion[2][0][0] == pytest.approx( + cocktail.thermal_expansion[2][0][0], rel=0.2 + ) + assert Path(doc.phonon_job_dir, "one_shot", "fc3.hdf5").exists() + + # the same force constants, with a q-point density instead of a mesh. The + # negative tolerance flags every frequency, to test the imaginary mode check. + phonon_output = responses[phonon_uuid][1].output + job = compute_cte( + phonon_output=phonon_output, + elastic_tensor=doc.elastic_tensor, + elastic_structure=doc.structure, + anhar_fit_methods=["one-shot"], + temperatures=[300], + mesh=100.0, + tol_imaginary_modes=-10.0, + ) + with pytest.warns(UserWarning, match="thermal expansion is not"): + responses = run_locally(job, create_folders=True, ensure_success=True) + flagged_doc = responses[job.uuid][1].output + assert flagged_doc.mesh == (2, 2, 2) + (flagged,) = flagged_doc.results + assert flagged.has_imaginary_modes + assert flagged.thermal_expansion is None + + # the non-analytical term correction with zero Born charges leaves alpha + # unchanged. pheasy only stores Born charges for VASP, so they are added here. + # The phonopy primitive cell has one atom and the unit cell four, so the + # charges must be expanded to four atoms before they reach phono3py. + phonon_job_dir = Path(doc.phonon_job_dir) + phonon = phonopy.load(phonon_job_dir / "phonopy.yaml", produce_fc=False) + assert len(phonon.primitive) == 1 + phonon.nac_params = { + "born": np.zeros((1, 3, 3)), + "dielectric": np.eye(3) * 10.0, + "factor": 14.399652, + } + phonon.save("phonopy_nac.yaml") + + nac_params_used = [] + original_gruneisen = gruneisen_module.Gruneisen + + def _gruneisen(*args, **kwargs): + nac_params_used.append(kwargs["nac_params"]) + return original_gruneisen(*args, **kwargs) + + monkeypatch.setattr(gruneisen_module, "Gruneisen", _gruneisen) + cocktail_files = { + "cocktail": (phonon_job_dir / "FORCE_CONSTANTS", phonon_job_dir / "fc3.hdf5") + } + settings = { + "force_constant_files": cocktail_files, + "elastic_tensor": doc.elastic_tensor, + "temperatures": [300], + "mesh": (8, 8, 8), + "tol_imaginary_modes": 0.1, + "min_frequency": 1e-3, + "symprec": 1e-5, + } + nac_doc = CTEDocument.from_force_constants( + phonopy_yaml="phonopy_nac.yaml", structure=doc.structure, **settings + ) + (nac_params,) = nac_params_used + assert nac_params["born"].shape == (4, 3, 3) + assert nac_params["dielectric"] == pytest.approx(np.eye(3) * 10.0) + assert np.array(nac_doc.results[0].thermal_expansion[0]) == pytest.approx( + np.array(cocktail.thermal_expansion[2]), rel=1e-6, abs=1e-15 + ) + + # a structure with another lattice or other atoms is refused + strained = doc.structure.copy() + strained.apply_strain(0.01) + with pytest.raises(ValueError, match="not in the same frame"): + CTEDocument.from_force_constants( + phonopy_yaml=phonon_job_dir / "phonopy.yaml", structure=strained, **settings + ) + other_atoms = doc.structure.copy() + other_atoms.replace_species({"Cu": "Au"}) + with pytest.raises(ValueError, match="atoms of the elastic calculation"): + CTEDocument.from_force_constants( + phonopy_yaml=phonon_job_dir / "phonopy.yaml", + structure=other_atoms, + **settings, + ) + with pytest.raises(ValueError, match="must not be negative"): + CTEDocument.from_force_constants( + phonopy_yaml=phonon_job_dir / "phonopy.yaml", + structure=doc.structure, + **{**settings, "temperatures": [-10, 300]}, + ) + + +def test_cte_maker_force_field_defaults(): + """The default phonon maker uses min_length=12.0 for the supercells.""" + maker = CTEMaker() + assert maker.phonon_maker.min_length == 12.0 + assert maker.phonon_maker.cal_anhar_fcs + assert maker.phonon_maker.displacement_anhar == 0.03 diff --git a/tests/vasp/flows/test_cte.py b/tests/vasp/flows/test_cte.py new file mode 100644 index 0000000000..95d8ce720c --- /dev/null +++ b/tests/vasp/flows/test_cte.py @@ -0,0 +1,71 @@ +import pytest +from jobflow import Flow +from pymatgen.core.structure import Structure + +from atomate2.vasp.flows.cte import CTEMaker +from atomate2.vasp.flows.elastic import ElasticMaker +from atomate2.vasp.flows.pheasy import PhononMaker + + +def test_cte_maker_vasp_flow(si_structure: Structure): + """The relaxed structure goes to both flows, and their outputs to compute_cte.""" + maker = CTEMaker(temperatures=[100, 300], mesh=(10, 10, 10)) + flow = maker.make(si_structure) + relax, phonon_flow, elastic_flow, cte = flow.jobs + assert isinstance(phonon_flow, Flow) + assert isinstance(elastic_flow, Flow) + + # both flows start from the relaxed structure + relaxed = relax.output.structure + for sub_flow in (phonon_flow, elastic_flow): + structure = sub_flow.jobs[0].function_args[0] + assert structure.uuid == relaxed.uuid + assert structure.attributes == relaxed.attributes + + kwargs = cte.function_kwargs + assert kwargs["phonon_output"].uuid == phonon_flow.output.uuid + assert kwargs["elastic_tensor"].uuid == elastic_flow.output.uuid + assert kwargs["elastic_structure"].uuid == elastic_flow.output.uuid + assert kwargs["anhar_fit_methods"] == ("one-shot",) + assert kwargs["temperatures"] == [100, 300] + assert kwargs["mesh"] == (10, 10, 10) + assert kwargs["symprec"] == maker.phonon_maker.symprec + assert kwargs["min_frequency"] == maker.min_frequency + assert maker.phonon_maker.min_length == 12.0 + + assert kwargs["elastic_tensor"].attributes == ( + ("a", "elastic_tensor"), + ("a", "raw"), + ) + + # the elastic fit gets the stress of the relaxation, as in the elastic flow + fit = next(job for job in elastic_flow.jobs if job.name == "fit_elastic_tensor") + stress = fit.function_kwargs["equilibrium_stress"] + assert stress.uuid == relax.output.uuid + assert stress.attributes == (("a", "output"), ("a", "stress")) + + +@pytest.mark.parametrize( + ("phonon_kwargs", "match"), + [ + ({"bulk_relax_maker": None, "cal_anhar_fcs": False}, "cal_anhar_fcs=True"), + ( + { + "bulk_relax_maker": None, + "cal_anhar_fcs": True, + "use_symmetrized_structure": "primitive", + }, + "use_symmetrized_structure=None", + ), + ({"cal_anhar_fcs": True}, "The phonon maker needs bulk_relax_maker=None"), + ], +) +def test_cte_maker_checks_phonon_maker(phonon_kwargs, match): + phonon_maker = PhononMaker(**phonon_kwargs) + with pytest.raises(ValueError, match=match): + CTEMaker(phonon_maker=phonon_maker) + + +def test_cte_maker_checks_elastic_maker(): + with pytest.raises(ValueError, match="The elastic maker needs bulk_relax_maker"): + CTEMaker(elastic_maker=ElasticMaker()) From 7d25de7884d85addd6c6fa4c56ba4a811bd72180 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Wed, 30 Sep 2026 09:11:40 -0700 Subject: [PATCH 08/28] Standardize the input cell in the thermal expansion workflow --- docs/user/codes/vasp.md | 4 +++- pyproject.toml | 3 ++- src/atomate2/common/flows/cte.py | 22 ++++++++++++++++++++-- src/atomate2/forcefields/flows/cte.py | 7 ++++++- src/atomate2/vasp/flows/cte.py | 7 ++++++- tests/common/jobs/test_cte.py | 7 ++++++- tests/vasp/flows/test_cte.py | 4 +++- 7 files changed, 46 insertions(+), 8 deletions(-) diff --git a/docs/user/codes/vasp.md b/docs/user/codes/vasp.md index 74f1a90fdd..a6891481f9 100644 --- a/docs/user/codes/vasp.md +++ b/docs/user/codes/vasp.md @@ -437,7 +437,9 @@ gruneisen_flow = GruneisenMaker( `CTEMaker` calculates the thermal expansion tensor from third-order force constants, with the help of [Pheasy](https://doi.org/10.48550/arXiv.2508.01020) and [phono3py](https://doi.org/10.1088/1361-648X/acd831). It needs the `pheasy` extra, see the Pheasy section above. -First, a tight structural relaxation is performed. +First, the structure is converted to the standard primitive cell, and a tight structural relaxation is performed. +The pheasy fits can fail for cells that are not in a standard setting. +Set `use_symmetrized_structure="conventional"` to use the standard conventional cell instead. The relaxed structure is then passed to the pheasy phonon workflow and to the elastic constant workflow. The two do not depend on each other, so a workflow manager can run them at the same time. Neither of them relaxes the structure again, so both use the same structure. diff --git a/pyproject.toml b/pyproject.toml index fd6c0aa590..41cf73f8e1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -55,7 +55,8 @@ phonons = ["phonopy>=2.43,<5", "seekpath>=2.0.0"] # the pheasy fork adds the --disp_matrix_file and --force_matrix_file options # and fixes the --scell option # phono3py computes the mode Grueneisen tensors in the thermal expansion workflow -pheasy = ["atomate2[phonons,alamode]", "hiphive==1.5", "numpy<=2.2", "pheasy @ git+https://gitlab.com/hpsahasrabuddhe/pheasy.git@1728f16642d5f2763f27695cd29ffada0f755d7f", "phono3py>=4.5"] +# phonors 0.5 breaks the Grueneisen q-point mesh of phonopy and phono3py 4.5 +pheasy = ["atomate2[phonons,alamode]", "hiphive==1.5", "numpy<=2.2", "pheasy @ git+https://gitlab.com/hpsahasrabuddhe/pheasy.git@1728f16642d5f2763f27695cd29ffada0f755d7f", "phono3py==4.5.0", "phonors<0.5"] alamode = ["alm @ git+https://github.com/ttadano/ALM.git@f1d668fdee66e7e7218a04c88daf19d0e14fce0c#subdirectory=python"] hiphive = ["hiphive==1.5", "trainstation>=1.0", "atomate2[phonons,alamode]"] lobster = ["ijson>=3.2.2", "lobsterpy>=0.6.0"] diff --git a/src/atomate2/common/flows/cte.py b/src/atomate2/common/flows/cte.py index d40fd25b26..95ec25a1e8 100644 --- a/src/atomate2/common/flows/cte.py +++ b/src/atomate2/common/flows/cte.py @@ -4,12 +4,13 @@ from abc import ABC, abstractmethod from dataclasses import dataclass, field -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, Literal from jobflow import Flow, Maker from pymatgen.util.due import Doi, due from atomate2.common.jobs.cte import compute_cte +from atomate2.common.jobs.utils import structure_to_conventional, structure_to_primitive if TYPE_CHECKING: from pathlib import Path @@ -35,7 +36,8 @@ class BaseCTEMaker(Maker, ABC): """ Maker to calculate the thermal expansion from third-order force constants. - A tight structural relaxation is performed first. The relaxed structure is + By default, the structure is first converted to the standard primitive + cell, and a tight structural relaxation follows. The relaxed structure is then passed to the pheasy phonon flow, which fits the second- and third-order force constants, and to the elastic flow. The two flows do not depend on each other, so a workflow manager can run them at the same time. @@ -51,6 +53,12 @@ class BaseCTEMaker(Maker, ABC): ---------- name: str Name of the flows produced by this maker. + use_symmetrized_structure: str or None + Convert the input structure to the standard "primitive" or + "conventional" cell before the relaxation. The phonon and elastic flows + both use the converted structure. The pheasy fits can fail for cells + that are not in a standard setting, so only set this to None for an + input structure that already is. bulk_relax_maker: .ForceFieldRelaxMaker, .BaseVaspMaker, or None A maker to perform a tight relaxation on the bulk. Set to None to skip the relaxation. @@ -74,6 +82,7 @@ class BaseCTEMaker(Maker, ABC): """ name: str = "cte" + use_symmetrized_structure: Literal["primitive", "conventional"] | None = "primitive" bulk_relax_maker: ForceFieldRelaxMaker | BaseVaspMaker | None = None phonon_maker: BasePhononMaker = None elastic_maker: BaseElasticMaker = None @@ -118,6 +127,15 @@ def make(self, structure: Structure, prev_dir: str | Path | None = None) -> Flow A previous calculation directory to use for copying outputs. """ jobs = [] + if self.use_symmetrized_structure == "primitive": + prim_job = structure_to_primitive(structure, self.phonon_maker.symprec) + jobs.append(prim_job) + structure = prim_job.output + elif self.use_symmetrized_structure == "conventional": + conv_job = structure_to_conventional(structure, self.phonon_maker.symprec) + jobs.append(conv_job) + structure = conv_job.output + equilibrium_stress = None if self.bulk_relax_maker is not None: bulk_kwargs = {} diff --git a/src/atomate2/forcefields/flows/cte.py b/src/atomate2/forcefields/flows/cte.py index ed87a7fd28..e2270e8bde 100644 --- a/src/atomate2/forcefields/flows/cte.py +++ b/src/atomate2/forcefields/flows/cte.py @@ -60,7 +60,8 @@ class CTEMaker(BaseCTEMaker): """ Maker to calculate the thermal expansion with a force field and pheasy. - A tight relaxation is performed first. The relaxed structure is then passed + By default, the structure is converted to the standard primitive cell and + relaxed tightly. The relaxed structure is then passed to the pheasy phonon flow, which fits the second- and third-order force constants, and to the elastic flow. Neither flow relaxes the structure again. Finally, phono3py gives the mode Grueneisen tensors, and the thermal @@ -76,6 +77,10 @@ class CTEMaker(BaseCTEMaker): ---------- name: str Name of the flows produced by this maker. + use_symmetrized_structure: str or None + Convert the input structure to the standard "primitive" or + "conventional" cell before the relaxation. The pheasy fits can fail + for cells that are not in a standard setting. bulk_relax_maker: .ForceFieldRelaxMaker or None A maker to perform a tight relaxation on the bulk. Set to None to skip the relaxation. diff --git a/src/atomate2/vasp/flows/cte.py b/src/atomate2/vasp/flows/cte.py index 9f63c4820c..dc7d691b0e 100644 --- a/src/atomate2/vasp/flows/cte.py +++ b/src/atomate2/vasp/flows/cte.py @@ -20,7 +20,8 @@ class CTEMaker(BaseCTEMaker): """ Maker to calculate the thermal expansion with VASP, pheasy and phono3py. - A tight double relaxation is performed first. The relaxed structure is then + By default, the structure is converted to the standard primitive cell, and + a tight double relaxation follows. The relaxed structure is then passed to the pheasy phonon flow, which fits the second- and third-order force constants, and to the elastic flow. Neither flow relaxes the structure again. Finally, phono3py gives the mode Grueneisen tensors, and @@ -39,6 +40,10 @@ class CTEMaker(BaseCTEMaker): ---------- name: str Name of the flows produced by this maker. + use_symmetrized_structure: str or None + Convert the input structure to the standard "primitive" or + "conventional" cell before the relaxation. The pheasy fits can fail + for cells that are not in a standard setting. bulk_relax_maker: .BaseVaspMaker or None A maker to perform a tight relaxation on the bulk. Set to None to skip the relaxation. diff --git a/tests/common/jobs/test_cte.py b/tests/common/jobs/test_cte.py index 8be61f9215..8dfddad6eb 100644 --- a/tests/common/jobs/test_cte.py +++ b/tests/common/jobs/test_cte.py @@ -147,8 +147,12 @@ def test_expand_born_to_unitcell(): def test_cte_maker_emt(clean_dir, monkeypatch): """Run the whole force field workflow with EMT forces on fcc Cu.""" structure = AseAtomsAdaptor.get_structure(bulk("Cu", "fcc", a=3.61, cubic=True)) + # the conventional cell keeps the 4-atom cubic unit cell used below maker = CTEMaker.from_force_field_name( - EMT, temperatures=[0, 100, 300], mesh=(8, 8, 8) + EMT, + use_symmetrized_structure="conventional", + temperatures=[0, 100, 300], + mesh=(8, 8, 8), ) # a 2x2x2 supercell of the cubic cell, 32 atoms, and an fc3 cutoff of 6 Bohr # (3.2 A), which covers the nearest neighbours at 2.55 A @@ -280,6 +284,7 @@ def _gruneisen(*args, **kwargs): def test_cte_maker_force_field_defaults(): """The default phonon maker uses min_length=12.0 for the supercells.""" maker = CTEMaker() + assert maker.use_symmetrized_structure == "primitive" assert maker.phonon_maker.min_length == 12.0 assert maker.phonon_maker.cal_anhar_fcs assert maker.phonon_maker.displacement_anhar == 0.03 diff --git a/tests/vasp/flows/test_cte.py b/tests/vasp/flows/test_cte.py index 95d8ce720c..98c392dde6 100644 --- a/tests/vasp/flows/test_cte.py +++ b/tests/vasp/flows/test_cte.py @@ -11,7 +11,9 @@ def test_cte_maker_vasp_flow(si_structure: Structure): """The relaxed structure goes to both flows, and their outputs to compute_cte.""" maker = CTEMaker(temperatures=[100, 300], mesh=(10, 10, 10)) flow = maker.make(si_structure) - relax, phonon_flow, elastic_flow, cte = flow.jobs + prim, relax, phonon_flow, elastic_flow, cte = flow.jobs + assert prim.name == "structure_to_primitive" + assert relax.jobs[0].function_args[0].uuid == prim.output.uuid assert isinstance(phonon_flow, Flow) assert isinstance(elastic_flow, Flow) From 2347c1844557e4d7c95a86f4de3e0ea17300b5a6 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Wed, 30 Sep 2026 09:42:21 -0700 Subject: [PATCH 09/28] Allow force field jobs without torch installed --- src/atomate2/forcefields/utils.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/atomate2/forcefields/utils.py b/src/atomate2/forcefields/utils.py index d6388d845a..2cd486e43e 100644 --- a/src/atomate2/forcefields/utils.py +++ b/src/atomate2/forcefields/utils.py @@ -496,7 +496,12 @@ def revert_default_dtype() -> Generator[None]: Originally added for use with MACE(Relax|Static)Maker. https://github.com/ACEsuit/mace/issues/328 """ - import torch + try: + import torch + except ImportError: + # force fields that do not use torch, such as EMT, have no dtype to revert + yield + return orig = torch.get_default_dtype() yield From 5bee251f41ad6383e488da18004ad93795d2beaf Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Wed, 30 Sep 2026 10:49:31 -0700 Subject: [PATCH 10/28] Add a thermal expansion tutorial --- .github/workflows/testing.yml | 2 +- docs/user/codes/vasp.md | 1 + tutorials/cte_workflow.ipynb | 355 ++++++++++++++++++++++++++++++++++ tutorials/tutorials.md | 1 + 4 files changed, 358 insertions(+), 1 deletion(-) create mode 100644 tutorials/cte_workflow.ipynb diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index 07e81ea6b8..e34e9e1c01 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -295,7 +295,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 + 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 - name: Test ASE env: diff --git a/docs/user/codes/vasp.md b/docs/user/codes/vasp.md index a6891481f9..31561ca004 100644 --- a/docs/user/codes/vasp.md +++ b/docs/user/codes/vasp.md @@ -486,6 +486,7 @@ Check that your VASP version runs DFPT with r2SCAN before you use this setting. The same workflow runs with a force field via `from atomate2.forcefields.flows.cte import CTEMaker`. `CTEMaker.from_force_field_name` sets one force field for the relaxation, the phonons and the elastic tensor. +The notebook `tutorials/cte_workflow.ipynb` runs it for MgO with MACE-OMAT-0-medium. ### Quasi-harmonic Workflow diff --git a/tutorials/cte_workflow.ipynb b/tutorials/cte_workflow.ipynb new file mode 100644 index 0000000000..d805497acf --- /dev/null +++ b/tutorials/cte_workflow.ipynb @@ -0,0 +1,355 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Thermal expansion workflow with pheasy and a machine-learned potential\n", + "\n", + "This notebook computes the thermal expansion of MgO with the `CTEMaker`\n", + "workflow. The forces come from MACE-OMAT-0-medium, so no DFT is needed." + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "## Background\n", + "\n", + "The workflow computes the thermal expansion tensor from mode Grüneisen tensors.\n", + "The thermal stress of each phonon mode is its heat capacity times its Grüneisen\n", + "tensor. The strain that relaxes the summed thermal stress is the thermal\n", + "expansion:\n", + "\n", + "$$\\alpha = \\frac{S}{N_q V} \\sum_{q\\nu} c_{q\\nu} \\gamma_{q\\nu}$$\n", + "\n", + "Here $S$ is the elastic compliance, $c_{q\\nu}$ and $\\gamma_{q\\nu}$ are the heat\n", + "capacity and Grüneisen tensor of each mode, $N_q$ is the number of q-points and\n", + "$V$ is the cell volume.\n", + "\n", + "The flow has four steps:\n", + "\n", + "1. The structure is converted to the standard primitive cell and relaxed tightly.\n", + "2. The pheasy phonon workflow fits the second- and third-order force constants\n", + " with LASSO. It uses randomly displaced supercells with 0.03 Å displacements.\n", + "3. The elastic workflow fits the elastic tensor on the same relaxed structure.\n", + "4. phono3py computes the mode Grüneisen tensors on a 12x12x12 q-point mesh, and\n", + " the thermal expansion follows from the formula above.\n", + "\n", + "Steps 2 and 3 do not depend on each other, so a workflow manager can run them\n", + "at the same time. The same workflow runs with VASP as\n", + "`atomate2.vasp.flows.cte.CTEMaker`." + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "## Installation\n", + "\n", + "The `pheasy` and `ase` extras are needed, plus MACE and the Materials Project\n", + "client:\n", + "\n", + "```\n", + "pip install 'atomate2[pheasy,ase]'\n", + "pip install 'mace-torch>=0.3.16' mp-api\n", + "```\n", + "\n", + "The `pheasy` extra pulls in pheasy, phono3py and ALM. ALM is compiled from\n", + "source. The hiPhive tutorial and the atomate2 VASP documentation list what the\n", + "build needs." + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "## The potential\n", + "\n", + "MACE-OMAT-0-medium downloads once and caches under `~/.cache/mace`. It is\n", + "released under the Academic Software License.\n", + "\n", + "`float64` is worth the cost here. The third-order force constants come from\n", + "small force differences between displaced supercells, where `float32` noise is\n", + "not negligible." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "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", + "warnings.filterwarnings(\"ignore\")\n", + "\n", + "from mace.calculators import mace_mp # noqa: E402\n", + "\n", + "CALC_KWARGS = {\"model\": \"medium-omat-0\", \"device\": \"cpu\", \"default_dtype\": \"float64\"}\n", + "\n", + "_ = mace_mp(**CALC_KWARGS) # downloads and caches on first use" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "## One calculator for all jobs\n", + "\n", + "Every force field job builds its own calculator. `run_locally` runs all of\n", + "them in this one Python process, and the workflow has a few hundred jobs.\n", + "Loading MACE for each of them grows the memory by about 150 MB per job. This\n", + "cell makes the jobs share one calculator. It is only needed when many jobs run\n", + "in one process, not with a workflow manager." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "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": "7", + "metadata": {}, + "source": [ + "## The structure\n", + "\n", + "The primitive cell of MgO (mp-1265) comes from the Materials Project. Set your\n", + "key first with `export MP_API_KEY=...`.\n", + "\n", + "The Materials Project serves this cell rotated away from the cubic axes, and\n", + "the pheasy fits fail for such a cell. `CTEMaker` therefore converts the input to\n", + "the standard primitive cell before the relaxation. This is the default,\n", + "`use_symmetrized_structure=\"primitive\"`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "from mp_api.client import MPRester\n", + "\n", + "if not os.environ.get(\"MP_API_KEY\"):\n", + " raise OSError(\"export MP_API_KEY before running this cell\")\n", + "\n", + "with MPRester() as mpr:\n", + " structure = mpr.get_structure_by_material_id(\"mp-1265\")\n", + "structure.lattice" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "## Building the workflow\n", + "\n", + "`CTEMaker.from_force_field_name` sets one force field for the relaxation, the\n", + "phonons and the elastic tensor. Any MACE model routes through `mace_mp`, and\n", + "the model name in `calculator_kwargs` selects the weights.\n", + "\n", + "Everything else stays at the workflow defaults. The supercells have\n", + "`min_length=12` Å, which gives 128 atoms for MgO. The third-order force\n", + "constants use the one-shot fit, and the temperatures run from 0 K to 1000 K in\n", + "steps of 10 K." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "from atomate2.forcefields.flows.cte import CTEMaker\n", + "\n", + "maker = CTEMaker.from_force_field_name(\"MACE-MP-0\", calculator_kwargs=CALC_KWARGS)\n", + "flow = maker.make(structure)\n", + "flow.draw_graph().show()" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "## Running the workflow\n", + "\n", + "The flow has a little over 200 jobs, most of them force calculations on the\n", + "displaced supercells. On one Perlmutter CPU node the whole run took about 45\n", + "minutes. Most of that time goes into the pheasy fit of the third-order force\n", + "constants. Expect longer on a laptop.\n", + "\n", + "`create_folders=True` is needed, because the thermal expansion job reads the\n", + "force constants from the folder of the pheasy fit." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12", + "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": "13", + "metadata": {}, + "source": [ + "## Results\n", + "\n", + "The output is a `CTEDocument`. It holds the thermal expansion tensor at each\n", + "temperature for each fit method. The value at 0 K is zero, since the heat\n", + "capacity is zero there. MgO is cubic, so the three diagonal components are\n", + "equal." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "(result,) = doc.results\n", + "temperatures = np.array(doc.temperatures)\n", + "alpha = np.array(result.thermal_expansion)\n", + "i300 = int(np.argmin(np.abs(temperatures - 300)))\n", + "\n", + "fig, ax = plt.subplots()\n", + "ax.plot(temperatures, alpha[:, 0, 0] * 1e6)\n", + "ax.set_xlabel(\"Temperature (K)\")\n", + "ax.set_ylabel(\"Linear thermal expansion (10$^{-6}$ K$^{-1}$)\")\n", + "plt.show()\n", + "\n", + "{\n", + " \"alpha(300 K) in 1/K\": float(alpha[i300, 0, 0]),\n", + " \"mean Grueneisen parameter at 300 K\": float(\n", + " np.trace(result.average_gruneisen[i300]) / 3\n", + " ),\n", + " \"has imaginary modes\": result.has_imaginary_modes,\n", + "}" + ] + }, + { + "cell_type": "markdown", + "id": "15", + "metadata": {}, + "source": [ + "With these settings we get $\\alpha$(300 K) = 9.6e-6 K$^{-1}$ and a mean\n", + "Grüneisen parameter of 1.36. Experiment gives about 1.0e-5 K$^{-1}$ for MgO at\n", + "room temperature." + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "## Checking the fit\n", + "\n", + "If the thermal expansion comes out as zero or very small, check the\n", + "third-order force constants first. All zeros means the LASSO fit dropped them.\n", + "The penalty chosen by cross validation should also lie inside the search\n", + "range, not on one of its bounds." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17", + "metadata": {}, + "outputs": [], + "source": [ + "import re\n", + "from pathlib import Path\n", + "\n", + "import h5py\n", + "\n", + "fit_dir = Path(doc.phonon_job_dir) / \"one_shot\"\n", + "log = (fit_dir / \"pheasy_anharmonic_fit.log\").read_text()\n", + "with h5py.File(fit_dir / \"fc3.hdf5\") as file:\n", + " fc3 = file[\"fc3\"][:]\n", + "\n", + "{\n", + " \"LASSO penalty\": re.findall(r\"alpha_(?:min|max|opt):\\s*\\S+\", log),\n", + " \"max |fc3| in eV/A^3\": float(np.abs(fc3).max()),\n", + "}" + ] + }, + { + "cell_type": "markdown", + "id": "18", + "metadata": {}, + "source": [ + "## Known limitations\n", + "\n", + "The frequencies and Grüneisen tensors are those of the relaxed structure. They\n", + "are not renormalized with temperature, so the result is least reliable at high\n", + "temperature.\n", + "\n", + "A machine-learned potential gives no Born charges, so the non-analytical term\n", + "correction is not applied here. With VASP, the Born charges of the phonon run\n", + "are used.\n", + "\n", + "If a frequency on the q-point mesh is below -0.1 THz, a warning is raised and\n", + "the thermal expansion of that fit is not computed." + ] + } + ], + "metadata": { + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/tutorials/tutorials.md b/tutorials/tutorials.md index e30858a4b9..487c6da5e9 100644 --- a/tutorials/tutorials.md +++ b/tutorials/tutorials.md @@ -18,6 +18,7 @@ pheasy_workflow hiphive_workflow force_fields/phonon_workflow grueneisen_workflow +cte_workflow qha_workflow torchsim_tutorial ``` From f3f113afc26be6803744a3fd79c06cf9b520eb99 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Thu, 1 Oct 2026 06:04:23 -0700 Subject: [PATCH 11/28] Pin pheasy to the commit with the faster sensing matrix --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 41cf73f8e1..979b8bbde6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -56,7 +56,7 @@ phonons = ["phonopy>=2.43,<5", "seekpath>=2.0.0"] # and fixes the --scell option # phono3py computes the mode Grueneisen tensors in the thermal expansion workflow # phonors 0.5 breaks the Grueneisen q-point mesh of phonopy and phono3py 4.5 -pheasy = ["atomate2[phonons,alamode]", "hiphive==1.5", "numpy<=2.2", "pheasy @ git+https://gitlab.com/hpsahasrabuddhe/pheasy.git@1728f16642d5f2763f27695cd29ffada0f755d7f", "phono3py==4.5.0", "phonors<0.5"] +pheasy = ["atomate2[phonons,alamode]", "hiphive==1.5", "numpy<=2.2", "pheasy @ git+https://gitlab.com/hpsahasrabuddhe/pheasy.git@9f24162a4ed0f0ab8911d382fd2617944aade55d", "phono3py==4.5.0", "phonors<0.5"] alamode = ["alm @ git+https://github.com/ttadano/ALM.git@f1d668fdee66e7e7218a04c88daf19d0e14fce0c#subdirectory=python"] hiphive = ["hiphive==1.5", "trainstation>=1.0", "atomate2[phonons,alamode]"] lobster = ["ijson>=3.2.2", "lobsterpy>=0.6.0"] From cecd36e912868ceadd7cf66d4707776c88c5f1cf Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Thu, 1 Oct 2026 07:17:16 -0700 Subject: [PATCH 12/28] Add the five-material and Si convergence sections to the thermal expansion tutorial --- tutorials/cte_workflow.ipynb | 350 ++++++++++++++++++++++++++++++++--- 1 file changed, 327 insertions(+), 23 deletions(-) diff --git a/tutorials/cte_workflow.ipynb b/tutorials/cte_workflow.ipynb index d805497acf..ada8d27c93 100644 --- a/tutorials/cte_workflow.ipynb +++ b/tutorials/cte_workflow.ipynb @@ -5,10 +5,14 @@ "id": "0", "metadata": {}, "source": [ - "# Thermal expansion workflow with pheasy and a machine-learned potential\n", + "# Thermal expansion workflow with pheasy and machine-learned potentials\n", "\n", - "This notebook computes the thermal expansion of MgO with the `CTEMaker`\n", - "workflow. The forces come from MACE-OMAT-0-medium, so no DFT is needed." + "This notebook computes the thermal expansion with the `CTEMaker` workflow. The\n", + "forces come from three MACE potentials, so no DFT is needed.\n", + "\n", + "Section 1 runs five materials with the default settings. MgO is worked through\n", + "step by step. Section 2 converges the third-order cutoff and the supercell size\n", + "for Si, where the default settings are not enough." ] }, { @@ -68,10 +72,23 @@ "id": "3", "metadata": {}, "source": [ - "## The potential\n", + "## The potentials\n", + "\n", + "The notebook uses three MACE foundation models. All three are released under the\n", + "Academic Software License.\n", + "\n", + "- **MACE-OMAT-0-medium.** `mace_mp(model=\"medium-omat-0\")` downloads it once and\n", + " caches it under `~/.cache/mace`.\n", + "- **MACE-MATPES-PBE-0** and **MACE-MATPES-r2SCAN-0.** These are MACE-OMAT-0\n", + " fine-tuned on the MatPES data set at the PBE and r2SCAN levels. They need\n", + " `mace-torch>=0.3.10`. Download the two files from the `mace_matpes_0` release of\n", + " [mace-foundations](https://github.com/ACEsuit/mace-foundations/releases/tag/mace_matpes_0):\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", + " - [MACE-matpes-r2scan-omat-ft.model](https://github.com/ACEsuit/mace-foundations/releases/download/mace_matpes_0/MACE-matpes-r2scan-omat-ft.model)\n", "\n", - "MACE-OMAT-0-medium downloads once and caches under `~/.cache/mace`. It is\n", - "released under the Academic Software License.\n", + "Set the paths to the two downloaded files in the next cell. The MgO example\n", + "below only needs MACE-OMAT-0-medium. Change `device` to `\"cuda\"` to run MACE on a\n", + "GPU.\n", "\n", "`float64` is worth the cost here. The third-order force constants come from\n", "small force differences between displaced supercells, where `float32` noise is\n", @@ -95,9 +112,20 @@ "\n", "from mace.calculators import mace_mp # noqa: E402\n", "\n", - "CALC_KWARGS = {\"model\": \"medium-omat-0\", \"device\": \"cpu\", \"default_dtype\": \"float64\"}\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", + "\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\": \"cpu\", \"default_dtype\": \"float64\"}\n", + "\n", "\n", - "_ = mace_mp(**CALC_KWARGS) # downloads and caches on first use" + "_ = mace_mp(**calc_kwargs(\"OMAT-0-medium\")) # downloads and caches on first use" ] }, { @@ -145,6 +173,17 @@ "cell_type": "markdown", "id": "7", "metadata": {}, + "source": [ + "## 1. Five materials with the default settings\n", + "\n", + "The first part runs MgO with MACE-OMAT-0-medium step by step. The other four\n", + "materials and the other two potentials follow with the same code." + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, "source": [ "## The structure\n", "\n", @@ -160,7 +199,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8", + "id": "9", "metadata": {}, "outputs": [], "source": [ @@ -176,7 +215,7 @@ }, { "cell_type": "markdown", - "id": "9", + "id": "10", "metadata": {}, "source": [ "## Building the workflow\n", @@ -194,28 +233,29 @@ { "cell_type": "code", "execution_count": null, - "id": "10", + "id": "11", "metadata": {}, "outputs": [], "source": [ "from atomate2.forcefields.flows.cte import CTEMaker\n", "\n", - "maker = CTEMaker.from_force_field_name(\"MACE-MP-0\", calculator_kwargs=CALC_KWARGS)\n", + "maker = CTEMaker.from_force_field_name(\n", + " \"MACE-MP-0\", calculator_kwargs=calc_kwargs(\"OMAT-0-medium\")\n", + ")\n", "flow = maker.make(structure)\n", "flow.draw_graph().show()" ] }, { "cell_type": "markdown", - "id": "11", + "id": "12", "metadata": {}, "source": [ "## Running the workflow\n", "\n", "The flow has a little over 200 jobs, most of them force calculations on the\n", - "displaced supercells. On one Perlmutter CPU node the whole run took about 45\n", - "minutes. Most of that time goes into the pheasy fit of the third-order force\n", - "constants. Expect longer on a laptop.\n", + "displaced supercells. On 32 CPU cores of a Perlmutter node the whole run took\n", + "about 11 minutes.\n", "\n", "`create_folders=True` is needed, because the thermal expansion job reads the\n", "force constants from the folder of the pheasy fit." @@ -224,7 +264,7 @@ { "cell_type": "code", "execution_count": null, - "id": "12", + "id": "13", "metadata": {}, "outputs": [], "source": [ @@ -238,7 +278,7 @@ }, { "cell_type": "markdown", - "id": "13", + "id": "14", "metadata": {}, "source": [ "## Results\n", @@ -252,7 +292,7 @@ { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "15", "metadata": {}, "outputs": [], "source": [ @@ -281,7 +321,7 @@ }, { "cell_type": "markdown", - "id": "15", + "id": "16", "metadata": {}, "source": [ "With these settings we get $\\alpha$(300 K) = 9.6e-6 K$^{-1}$ and a mean\n", @@ -291,7 +331,7 @@ }, { "cell_type": "markdown", - "id": "16", + "id": "17", "metadata": {}, "source": [ "## Checking the fit\n", @@ -305,7 +345,7 @@ { "cell_type": "code", "execution_count": null, - "id": "17", + "id": "18", "metadata": {}, "outputs": [], "source": [ @@ -327,7 +367,271 @@ }, { "cell_type": "markdown", - "id": "18", + "id": "19", + "metadata": {}, + "source": [ + "## The other materials and potentials\n", + "\n", + "The same workflow runs for NaCl, KCl, CaO and GaAs, and with all three\n", + "potentials. The function below builds and runs one workflow in its own folder.\n", + "Section 2 uses it as well.\n", + "\n", + "The loop runs 15 workflows, so it is off by default. Set `RUN_ALL = True` to\n", + "run it. On 32 CPU cores of a Perlmutter node each run took between 1 and 26\n", + "minutes." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20", + "metadata": {}, + "outputs": [], + "source": [ + "from pymatgen.core import Structure\n", + "\n", + "from atomate2.common.schemas.cte import CTEDocument\n", + "\n", + "MATERIALS = {\n", + " \"NaCl\": \"mp-22862\",\n", + " \"KCl\": \"mp-23193\",\n", + " \"MgO\": \"mp-1265\",\n", + " \"CaO\": \"mp-2605\",\n", + " \"GaAs\": \"mp-2534\",\n", + "}\n", + "\n", + "\n", + "def run_cte(\n", + " structure: Structure,\n", + " model: str,\n", + " root_dir: str,\n", + " min_length: float | None = None,\n", + " c3: float | None = None,\n", + ") -> CTEDocument:\n", + " \"\"\"Run the thermal expansion workflow and return the CTEDocument.\n", + "\n", + " min_length sets the supercell size in Angstrom. c3 sets the third-order\n", + " cutoff in Bohr. None keeps the workflow default.\n", + " \"\"\"\n", + " maker = CTEMaker.from_force_field_name(\n", + " \"MACE-MP-0\", calculator_kwargs=calc_kwargs(model)\n", + " )\n", + " if min_length is not None:\n", + " maker.phonon_maker.min_length = min_length\n", + " if c3 is not None:\n", + " maker.phonon_maker.fcs_cutoff_radius = [-1, c3, 10]\n", + " flow = maker.make(structure)\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", + "def alpha_300k(doc: CTEDocument) -> float:\n", + " \"\"\"Linear thermal expansion at 300 K in 1/K, from the xx component.\"\"\"\n", + " (result,) = doc.results\n", + " i300 = int(np.argmin(np.abs(np.array(doc.temperatures) - 300)))\n", + " return float(np.array(result.thermal_expansion)[i300, 0, 0])\n", + "\n", + "\n", + "RUN_ALL = False\n", + "if RUN_ALL:\n", + " with MPRester() as mpr:\n", + " structures = {\n", + " name: mpr.get_structure_by_material_id(mp_id)\n", + " for name, mp_id in MATERIALS.items()\n", + " }\n", + " alpha_all = {\n", + " (name, model): alpha_300k(run_cte(structure, model, f\"runs/{name}_{model}\"))\n", + " for name, structure in structures.items()\n", + " for model in MODEL_FILES\n", + " }" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21", + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "\n", + "# alpha(300 K) in 1e-6/K from our runs with the default settings\n", + "ALPHA_DEFAULT = {\n", + " \"OMAT-0-medium\": {\"NaCl\": 35.9, \"KCl\": 39.0, \"MgO\": 9.6, \"CaO\": 12.5, \"GaAs\": 5.9},\n", + " \"MATPES-PBE-0\": {\"NaCl\": 38.1, \"KCl\": 34.8, \"MgO\": 21.5, \"CaO\": 15.2, \"GaAs\": 3.9},\n", + " \"MATPES-r2SCAN-0\": {\n", + " \"NaCl\": 30.9,\n", + " \"KCl\": 37.5,\n", + " \"MgO\": 11.0,\n", + " \"CaO\": 11.7,\n", + " \"GaAs\": 3.7,\n", + " },\n", + "}\n", + "\n", + "pd.DataFrame(ALPHA_DEFAULT)" + ] + }, + { + "cell_type": "markdown", + "id": "22", + "metadata": {}, + "source": [ + "For NaCl, KCl and CaO the three potentials agree within about 30%. For GaAs\n", + "MACE-OMAT-0-medium gives about 50% more than the two MatPES potentials.\n", + "\n", + "MgO stands out. MATPES-PBE-0 gives about twice the value of the other two\n", + "potentials and of experiment. It also gives a much softer MgO. Its bulk modulus\n", + "is 100 GPa, against 154 GPa with MACE-OMAT-0-medium and 163 GPa with\n", + "MATPES-r2SCAN-0. A softer lattice expands more.\n", + "\n", + "We did not converge the settings for these five materials. Section 2 shows how\n", + "to do this for Si, where the default settings are not enough." + ] + }, + { + "cell_type": "markdown", + "id": "23", + "metadata": {}, + "source": [ + "## 2. Converging the cutoff and the supercell for Si\n", + "\n", + "In Si the transverse acoustic modes near the zone boundary have negative\n", + "Grüneisen parameters. Their contribution to the thermal stress partly cancels\n", + "that of the other modes. The thermal expansion is the small difference of two\n", + "larger terms. Small errors in the third-order force constants therefore change\n", + "it a lot. This is why the default settings are not enough for Si.\n", + "\n", + "Two settings control the third-order force constants:\n", + "\n", + "- `phonon_maker.min_length` sets the supercell. For the Si primitive cell, 12, 16\n", + " and 21 Å give the 4x4x4, 5x5x5 and 6x6x6 supercells with 128, 250 and 432\n", + " atoms. The default is 12 Å.\n", + "- `phonon_maker.fcs_cutoff_radius` sets the cutoff radius of each order in Bohr.\n", + " The second entry is the third-order cutoff. The default is\n", + " `[-1, 12, 10]`, so 12 Bohr.\n", + "\n", + "A cutoff should stay below half the shortest distance between periodic images\n", + "of the supercell. Above it an atom starts to interact with its own images. For\n", + "Si this limit is 14.6, 18.3 and 21.9 Bohr for the 4x4x4, 5x5x5 and 6x6x6\n", + "supercells.\n", + "\n", + "A larger cutoff adds third-order force constants. The workflow picks the number\n", + "of displaced supercells so that the fit has about 100 equations per free force\n", + "constant, up to 600 supercells. The cost therefore grows quickly with the\n", + "cutoff. The 6x6x6 run at 18 Bohr used 201 displaced supercells and took about 70\n", + "minutes on 32 CPU cores of a Perlmutter node. At 20 Bohr it used 527 supercells, and the LASSO\n", + "fit needed more than the 220 GB of memory we gave it. We stopped at 18 Bohr.\n", + "\n", + "The scan below runs 11 settings for each potential, so it is off by default." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24", + "metadata": {}, + "outputs": [], + "source": [ + "SI_SCAN = { # min_length in Angstrom: third-order cutoffs in Bohr\n", + " 12: (12, 14, 16),\n", + " 16: (12, 14, 16, 18),\n", + " 21: (12, 14, 16, 18),\n", + "}\n", + "\n", + "RUN_SI_SCAN = False\n", + "if RUN_SI_SCAN:\n", + " with MPRester() as mpr:\n", + " si = mpr.get_structure_by_material_id(\"mp-149\")\n", + " alpha_si = {\n", + " (model, min_length, c3): alpha_300k(\n", + " run_cte(\n", + " si, model, f\"runs/Si_{model}_{min_length}A_{c3}bohr\", min_length, c3\n", + " )\n", + " )\n", + " for model in MODEL_FILES\n", + " for min_length, cutoffs in SI_SCAN.items()\n", + " for c3 in cutoffs\n", + " }" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "25", + "metadata": {}, + "outputs": [], + "source": [ + "# alpha(300 K) in 1e-6/K from our runs, by potential, supercell and cutoff in Bohr\n", + "ALPHA_SI = {\n", + " \"OMAT-0-medium\": {\n", + " \"4x4x4\": {12: 0.12, 14: 0.20, 16: 0.16},\n", + " \"5x5x5\": {12: 0.07, 14: 0.25, 16: 0.26, 18: 0.20},\n", + " \"6x6x6\": {12: 0.15, 14: 0.21, 16: 0.14, 18: 0.20},\n", + " },\n", + " \"MATPES-PBE-0\": {\n", + " \"4x4x4\": {12: 1.63, 14: 1.87, 16: 1.87},\n", + " \"5x5x5\": {12: 1.77, 14: 1.88, 16: 1.95, 18: 1.94},\n", + " \"6x6x6\": {12: 1.92, 14: 2.01, 16: 1.80, 18: 1.78},\n", + " },\n", + " \"MATPES-r2SCAN-0\": {\n", + " \"4x4x4\": {12: 0.68, 14: 1.92, 16: 1.87},\n", + " \"5x5x5\": {12: -0.82, 14: 1.85, 16: 1.94, 18: 1.94},\n", + " \"6x6x6\": {12: 0.99, 14: 2.00, 16: 1.89, 18: 1.96},\n", + " },\n", + "}\n", + "IMAGE_LIMIT = {\"4x4x4\": 14.6, \"5x5x5\": 18.3, \"6x6x6\": 21.9} # Bohr\n", + "\n", + "fig, axes = plt.subplots(1, 3, figsize=(12, 3.6), sharey=True)\n", + "for ax, (model, by_cell) in zip(axes, ALPHA_SI.items(), strict=True):\n", + " for cell, by_cutoff in by_cell.items():\n", + " cutoffs = np.array(list(by_cutoff))\n", + " values = np.array(list(by_cutoff.values()))\n", + " (line,) = ax.plot(cutoffs, values, \"o-\", label=cell)\n", + " beyond = cutoffs > IMAGE_LIMIT[cell]\n", + " ax.plot(\n", + " cutoffs[beyond], values[beyond], \"o\", color=line.get_color(), mfc=\"white\"\n", + " )\n", + " ax.axhline(2.6, color=\"gray\", ls=\"--\", label=\"experiment\")\n", + " ax.set_title(model)\n", + " ax.set_xlabel(\"Third-order cutoff (Bohr)\")\n", + "axes[0].set_ylabel(\"$\\\\alpha$(300 K) (10$^{-6}$ K$^{-1}$)\")\n", + "axes[0].legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "26", + "metadata": {}, + "source": [ + "Open markers are cutoffs above the image limit of the supercell. The dashed\n", + "line is the experimental value of about 2.6e-6 K$^{-1}$ (Okada and Tokumaru,\n", + "J. Appl. Phys. 56, 314 (1984)).\n", + "\n", + "- The default 12 Bohr cutoff is not converged for Si. With MATPES-r2SCAN-0 the\n", + " 5x5x5 supercell even gives a negative value.\n", + "- In the 5x5x5 supercell the value is converged at 16 Bohr. MATPES-PBE-0 and\n", + " MATPES-r2SCAN-0 both give 1.94e-6 K$^{-1}$ at 16 and 18 Bohr.\n", + "- The 6x6x6 supercell at 16 and 18 Bohr gives 1.8e-6 to 2.0e-6 K$^{-1}$. This is\n", + " within about 8% of the 5x5x5 value.\n", + "- MACE-OMAT-0-medium gives about 0.2e-6 K$^{-1}$ for every setting. The settings\n", + " do not change this. It comes from the potential.\n", + "\n", + "For a new material, converge the cutoff in a fixed supercell first. Keep the\n", + "cutoff below the image limit. Then repeat the converged cutoff in a larger\n", + "supercell to check the supercell size." + ] + }, + { + "cell_type": "markdown", + "id": "27", "metadata": {}, "source": [ "## Known limitations\n", From 32d7441ee24c1fab62daaf085c7dd84f65069143 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Thu, 1 Oct 2026 07:48:20 -0700 Subject: [PATCH 13/28] Write the K^-1 units as plain text so the tutorial renders in VS Code and on GitHub --- tutorials/cte_workflow.ipynb | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tutorials/cte_workflow.ipynb b/tutorials/cte_workflow.ipynb index ada8d27c93..850e9aaa48 100644 --- a/tutorials/cte_workflow.ipynb +++ b/tutorials/cte_workflow.ipynb @@ -324,8 +324,8 @@ "id": "16", "metadata": {}, "source": [ - "With these settings we get $\\alpha$(300 K) = 9.6e-6 K$^{-1}$ and a mean\n", - "Grüneisen parameter of 1.36. Experiment gives about 1.0e-5 K$^{-1}$ for MgO at\n", + "With these settings we get $\\alpha$(300 K) = 9.6e-6 K⁻¹ and a mean\n", + "Grüneisen parameter of 1.36. Experiment gives about 1.0e-5 K⁻¹ for MgO at\n", "room temperature." ] }, @@ -612,16 +612,16 @@ "metadata": {}, "source": [ "Open markers are cutoffs above the image limit of the supercell. The dashed\n", - "line is the experimental value of about 2.6e-6 K$^{-1}$ (Okada and Tokumaru,\n", + "line is the experimental value of about 2.6e-6 K⁻¹ (Okada and Tokumaru,\n", "J. Appl. Phys. 56, 314 (1984)).\n", "\n", "- The default 12 Bohr cutoff is not converged for Si. With MATPES-r2SCAN-0 the\n", " 5x5x5 supercell even gives a negative value.\n", "- In the 5x5x5 supercell the value is converged at 16 Bohr. MATPES-PBE-0 and\n", - " MATPES-r2SCAN-0 both give 1.94e-6 K$^{-1}$ at 16 and 18 Bohr.\n", - "- The 6x6x6 supercell at 16 and 18 Bohr gives 1.8e-6 to 2.0e-6 K$^{-1}$. This is\n", + " MATPES-r2SCAN-0 both give 1.94e-6 K⁻¹ at 16 and 18 Bohr.\n", + "- The 6x6x6 supercell at 16 and 18 Bohr gives 1.8e-6 to 2.0e-6 K⁻¹. This is\n", " within about 8% of the 5x5x5 value.\n", - "- MACE-OMAT-0-medium gives about 0.2e-6 K$^{-1}$ for every setting. The settings\n", + "- MACE-OMAT-0-medium gives about 0.2e-6 K⁻¹ for every setting. The settings\n", " do not change this. It comes from the potential.\n", "\n", "For a new material, converge the cutoff in a fixed supercell first. Keep the\n", From c27615d3cbfd42e9bcb6067aafec32a78620576e Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Thu, 1 Oct 2026 15:09:13 -0700 Subject: [PATCH 14/28] Converge the anharmonic LASSO fit with --tol 1e-8 --- src/atomate2/common/jobs/pheasy.py | 6 ++++-- tests/common/jobs/test_pheasy.py | 4 ++++ 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/atomate2/common/jobs/pheasy.py b/src/atomate2/common/jobs/pheasy.py index f882e3580e..18e5d8c2c5 100644 --- a/src/atomate2/common/jobs/pheasy.py +++ b/src/atomate2/common/jobs/pheasy.py @@ -337,10 +337,12 @@ def _run_anharmonic_fit( f"--disp_matrix_file {_DEFAULT_FILE_PATHS['anharmonic_displacements']}" ), # LASSO fit. OLS, pheasy's default, gives dense force constants. - # --std and --rasr are not passed to the anharmonic fit. + # --std and --rasr are not passed to the anharmonic fit. With pheasy's + # default --tol of 1e-4 the fits at small penalties do not converge, and + # the penalty chosen by cross-validation can change between machines. ( f"{base} -f {fix_fc2}-l LASSO --alpha_min {anhar_alpha_min} {seed}" - f"--ndata {int(num_anhar)} --hdf5 " + f"--tol 1e-8 --ndata {int(num_anhar)} --hdf5 " f"--force_matrix_file {_DEFAULT_FILE_PATHS['anharmonic_force_matrix']} " f"-o {log_file}" ), diff --git a/tests/common/jobs/test_pheasy.py b/tests/common/jobs/test_pheasy.py index 6bb45d0b0f..fc1d816d29 100644 --- a/tests/common/jobs/test_pheasy.py +++ b/tests/common/jobs/test_pheasy.py @@ -307,9 +307,13 @@ def test_anharmonic_fit_cocktail_and_one_shot(tmp_dir): # 20, then the undisplaced supercell assert len(displacements) == 1 + 20 + 1 + # the converged cocktail fit picks a penalty of about 1e-12, the default + # lower bound, so the bound is lowered to keep the penalty inside the + # search. In local runs the force constants did not change with the bound. job = generate_frequencies_eigenvectors( structure=structure, displacement_data=_emt_displacement_data(displacements), + anhar_alpha_min=-14, **anhar_kwargs, **FIT_KWARGS, **COMMON_KWARGS, From 42b98c290f3de2660c27a2e821b6833166a52453 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Thu, 1 Oct 2026 17:06:22 -0700 Subject: [PATCH 15/28] Update the thermal expansion tutorial with converged fits and finite-displacement references --- tutorials/cte_workflow.ipynb | 207 +++++++++++++++++++++++++++-------- 1 file changed, 160 insertions(+), 47 deletions(-) diff --git a/tutorials/cte_workflow.ipynb b/tutorials/cte_workflow.ipynb index 850e9aaa48..1bcb8ab8dc 100644 --- a/tutorials/cte_workflow.ipynb +++ b/tutorials/cte_workflow.ipynb @@ -324,9 +324,8 @@ "id": "16", "metadata": {}, "source": [ - "With these settings we get $\\alpha$(300 K) = 9.6e-6 K⁻¹ and a mean\n", - "Grüneisen parameter of 1.36. Experiment gives about 1.0e-5 K⁻¹ for MgO at\n", - "room temperature." + "With these settings we get $\\alpha$(300 K) = 11.2e-6 K⁻¹. Experiment gives\n", + "about 1.0e-5 K⁻¹ for MgO at room temperature." ] }, { @@ -464,18 +463,36 @@ "\n", "# alpha(300 K) in 1e-6/K from our runs with the default settings\n", "ALPHA_DEFAULT = {\n", - " \"OMAT-0-medium\": {\"NaCl\": 35.9, \"KCl\": 39.0, \"MgO\": 9.6, \"CaO\": 12.5, \"GaAs\": 5.9},\n", - " \"MATPES-PBE-0\": {\"NaCl\": 38.1, \"KCl\": 34.8, \"MgO\": 21.5, \"CaO\": 15.2, \"GaAs\": 3.9},\n", + " \"OMAT-0-medium\": {\"NaCl\": 39.4, \"KCl\": 41.6, \"MgO\": 11.2, \"CaO\": 13.9, \"GaAs\": 5.9},\n", + " \"MATPES-PBE-0\": {\"NaCl\": 44.4, \"KCl\": 47.5, \"MgO\": 31.4, \"CaO\": 14.3, \"GaAs\": 4.4},\n", " \"MATPES-r2SCAN-0\": {\n", - " \"NaCl\": 30.9,\n", - " \"KCl\": 37.5,\n", - " \"MgO\": 11.0,\n", - " \"CaO\": 11.7,\n", - " \"GaAs\": 3.7,\n", + " \"NaCl\": 32.8,\n", + " \"KCl\": 17.7,\n", + " \"MgO\": 11.3,\n", + " \"CaO\": 11.1,\n", + " \"GaAs\": 2.5,\n", + " },\n", + "}\n", + "# finite displacements in the largest supercell we ran, 686 to 2662 atoms\n", + "ALPHA_FD_LARGE = {\n", + " \"OMAT-0-medium\": {\"NaCl\": 39.3, \"KCl\": 41.8, \"MgO\": 11.0, \"CaO\": 14.0, \"GaAs\": 6.7},\n", + " \"MATPES-PBE-0\": {\"NaCl\": 40.0, \"KCl\": 44.4, \"MgO\": 39.2, \"CaO\": 15.8, \"GaAs\": 6.7},\n", + " \"MATPES-r2SCAN-0\": {\n", + " \"NaCl\": 33.7,\n", + " \"KCl\": 19.4,\n", + " \"MgO\": 12.3,\n", + " \"CaO\": 11.8,\n", + " \"GaAs\": 4.1,\n", " },\n", "}\n", "\n", - "pd.DataFrame(ALPHA_DEFAULT)" + "pd.concat(\n", + " {\n", + " \"default settings\": pd.DataFrame(ALPHA_DEFAULT),\n", + " \"finite displacements\": pd.DataFrame(ALPHA_FD_LARGE),\n", + " },\n", + " axis=1,\n", + ")" ] }, { @@ -483,16 +500,34 @@ "id": "22", "metadata": {}, "source": [ - "For NaCl, KCl and CaO the three potentials agree within about 30%. For GaAs\n", - "MACE-OMAT-0-medium gives about 50% more than the two MatPES potentials.\n", - "\n", - "MgO stands out. MATPES-PBE-0 gives about twice the value of the other two\n", - "potentials and of experiment. It also gives a much softer MgO. Its bulk modulus\n", - "is 100 GPa, against 154 GPa with MACE-OMAT-0-medium and 163 GPa with\n", - "MATPES-r2SCAN-0. A softer lattice expands more.\n", - "\n", - "We did not converge the settings for these five materials. Section 2 shows how\n", - "to do this for Si, where the default settings are not enough." + "The first set of columns is what the workflow gives with the default settings.\n", + "To check it, we computed the force constants of each run by finite\n", + "displacements with phono3py, with no cutoff and no fit. In the supercell of\n", + "the run the fits agree with finite displacements within 6%. The one exception\n", + "is GaAs with MATPES-r2SCAN-0, at 10%. So the fit is right for the supercell it\n", + "uses.\n", + "\n", + "We then repeated the finite displacements in larger supercells. The second set\n", + "of columns gives the value in the largest one, with 686 to 2662 atoms. Between\n", + "the two largest supercells each value changes by about 3% or less.\n", + "\n", + "- With MACE-OMAT-0-medium the default settings are within 2% of the\n", + " large-supercell value for the four rock-salt materials.\n", + "- With the two MatPES potentials the default settings are off by up to 11% for\n", + " the rock-salt materials, and by 20% for MgO with MATPES-PBE-0.\n", + "- For GaAs the default settings are 11% to 38% too low with all three\n", + " potentials. GaAs has the zinc-blende structure of Si. Section 2 shows the\n", + " same effect for Si in more detail.\n", + "\n", + "The potentials differ much more than these errors. Two results stand out.\n", + "\n", + "- KCl with MATPES-r2SCAN-0 gives 19.4e-6 K⁻¹ in the large supercell. This is\n", + " less than half the value of the other two potentials.\n", + "- MgO with MATPES-PBE-0 gives 39.2e-6 K⁻¹. This is more than three times the\n", + " value of the other two potentials and of experiment. MATPES-PBE-0 also gives\n", + " a much softer MgO. Its bulk modulus is 100 GPa, against 154 GPa with\n", + " MACE-OMAT-0-medium and 163 GPa with MATPES-r2SCAN-0. A softer lattice\n", + " expands more." ] }, { @@ -525,9 +560,9 @@ "A larger cutoff adds third-order force constants. The workflow picks the number\n", "of displaced supercells so that the fit has about 100 equations per free force\n", "constant, up to 600 supercells. The cost therefore grows quickly with the\n", - "cutoff. The 6x6x6 run at 18 Bohr used 201 displaced supercells and took about 70\n", - "minutes on 32 CPU cores of a Perlmutter node. At 20 Bohr it used 527 supercells, and the LASSO\n", - "fit needed more than the 220 GB of memory we gave it. We stopped at 18 Bohr.\n", + "cutoff. The 6x6x6 run at 18 Bohr used 201 displaced supercells. At 20 Bohr it used 527,\n", + "and the LASSO fit needed more than the 220 GB of memory we gave it. We stopped\n", + "at 18 Bohr.\n", "\n", "The scan below runs 11 settings for each potential, so it is off by default." ] @@ -571,22 +606,28 @@ "# alpha(300 K) in 1e-6/K from our runs, by potential, supercell and cutoff in Bohr\n", "ALPHA_SI = {\n", " \"OMAT-0-medium\": {\n", - " \"4x4x4\": {12: 0.12, 14: 0.20, 16: 0.16},\n", - " \"5x5x5\": {12: 0.07, 14: 0.25, 16: 0.26, 18: 0.20},\n", - " \"6x6x6\": {12: 0.15, 14: 0.21, 16: 0.14, 18: 0.20},\n", + " \"4x4x4\": {12: 1.42, 14: 1.90, 16: 1.87},\n", + " \"5x5x5\": {12: 1.53, 14: 2.04, 16: 2.27, 18: 2.06},\n", + " \"6x6x6\": {12: 1.64, 14: 1.99, 16: 2.30, 18: 2.07},\n", " },\n", " \"MATPES-PBE-0\": {\n", - " \"4x4x4\": {12: 1.63, 14: 1.87, 16: 1.87},\n", - " \"5x5x5\": {12: 1.77, 14: 1.88, 16: 1.95, 18: 1.94},\n", - " \"6x6x6\": {12: 1.92, 14: 2.01, 16: 1.80, 18: 1.78},\n", + " \"4x4x4\": {12: 3.25, 14: 4.13, 16: 4.06},\n", + " \"5x5x5\": {12: 3.35, 14: 4.30, 16: 3.91, 18: 3.53},\n", + " \"6x6x6\": {12: 3.43, 14: 4.32, 16: 3.85, 18: 3.48},\n", " },\n", " \"MATPES-r2SCAN-0\": {\n", - " \"4x4x4\": {12: 0.68, 14: 1.92, 16: 1.87},\n", - " \"5x5x5\": {12: -0.82, 14: 1.85, 16: 1.94, 18: 1.94},\n", - " \"6x6x6\": {12: 0.99, 14: 2.00, 16: 1.89, 18: 1.96},\n", + " \"4x4x4\": {12: -0.83, 14: -0.46, 16: -0.49},\n", + " \"5x5x5\": {12: -0.85, 14: -0.32, 16: -0.61, 18: 0.03},\n", + " \"6x6x6\": {12: -0.52, 14: -0.37, 16: -0.64, 18: 0.08},\n", " },\n", "}\n", "IMAGE_LIMIT = {\"4x4x4\": 14.6, \"5x5x5\": 18.3, \"6x6x6\": 21.9} # Bohr\n", + "# finite displacements in the 10x10x10 supercell, from the next section\n", + "ALPHA_FD_CONVERGED = {\n", + " \"OMAT-0-medium\": 2.90,\n", + " \"MATPES-PBE-0\": 3.71,\n", + " \"MATPES-r2SCAN-0\": 0.58,\n", + "}\n", "\n", "fig, axes = plt.subplots(1, 3, figsize=(12, 3.6), sharey=True)\n", "for ax, (model, by_cell) in zip(axes, ALPHA_SI.items(), strict=True):\n", @@ -598,6 +639,9 @@ " ax.plot(\n", " cutoffs[beyond], values[beyond], \"o\", color=line.get_color(), mfc=\"white\"\n", " )\n", + " ax.axhline(\n", + " ALPHA_FD_CONVERGED[model], color=\"black\", ls=\":\", label=\"finite displacements\"\n", + " )\n", " ax.axhline(2.6, color=\"gray\", ls=\"--\", label=\"experiment\")\n", " ax.set_title(model)\n", " ax.set_xlabel(\"Third-order cutoff (Bohr)\")\n", @@ -611,28 +655,97 @@ "id": "26", "metadata": {}, "source": [ - "Open markers are cutoffs above the image limit of the supercell. The dashed\n", - "line is the experimental value of about 2.6e-6 K⁻¹ (Okada and Tokumaru,\n", + "Open markers are cutoffs above the image limit of the supercell. The dotted\n", + "line is the converged finite-displacement value from the next section. The\n", + "dashed line is the experimental value of about 2.6e-6 K⁻¹ (Okada and Tokumaru,\n", "J. Appl. Phys. 56, 314 (1984)).\n", "\n", - "- The default 12 Bohr cutoff is not converged for Si. With MATPES-r2SCAN-0 the\n", - " 5x5x5 supercell even gives a negative value.\n", - "- In the 5x5x5 supercell the value is converged at 16 Bohr. MATPES-PBE-0 and\n", - " MATPES-r2SCAN-0 both give 1.94e-6 K⁻¹ at 16 and 18 Bohr.\n", - "- The 6x6x6 supercell at 16 and 18 Bohr gives 1.8e-6 to 2.0e-6 K⁻¹. This is\n", - " within about 8% of the 5x5x5 value.\n", - "- MACE-OMAT-0-medium gives about 0.2e-6 K⁻¹ for every setting. The settings\n", - " do not change this. It comes from the potential.\n", - "\n", - "For a new material, converge the cutoff in a fixed supercell first. Keep the\n", - "cutoff below the image limit. Then repeat the converged cutoff in a larger\n", - "supercell to check the supercell size." + "- The default 12 Bohr cutoff is far from the converged value. With\n", + " MACE-OMAT-0-medium it gives about half of it. With MATPES-r2SCAN-0 even the\n", + " sign is wrong.\n", + "- Between 12 and 18 Bohr the value goes up and down as each new shell of\n", + " neighbors enters the fit. It does not settle within the cutoffs we could\n", + " afford.\n", + "- From 14 Bohr on, the value at a fixed cutoff changes little between the 5x5x5\n", + " and 6x6x6 supercells. The cutoff limits the result, more than the supercell." ] }, { "cell_type": "markdown", "id": "27", "metadata": {}, + "source": [ + "## Finite displacements as a reference\n", + "\n", + "To find the converged value, we computed the force constants of Si by finite\n", + "displacements with phono3py, in supercells from 4x4x4 to 10x10x10. There is no\n", + "cutoff and no fit. The 10x10x10 supercell has 2000 atoms and needed 6201\n", + "displaced supercells for each potential." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "28", + "metadata": {}, + "outputs": [], + "source": [ + "# alpha(300 K) in 1e-6/K from finite displacements, by potential and supercell\n", + "ALPHA_FD = {\n", + " \"OMAT-0-medium\": {4: 1.74, 5: 1.96, 6: 2.32, 7: 2.64, 8: 2.83, 9: 2.89, 10: 2.90},\n", + " \"MATPES-PBE-0\": {4: 3.75, 5: 3.41, 6: 3.45, 7: 3.58, 8: 3.67, 9: 3.70, 10: 3.71},\n", + " \"MATPES-r2SCAN-0\": {\n", + " 4: -0.53,\n", + " 5: -0.10,\n", + " 6: 0.14,\n", + " 7: 0.44,\n", + " 8: 0.56,\n", + " 9: 0.57,\n", + " 10: 0.58,\n", + " },\n", + "}\n", + "\n", + "fig, ax = plt.subplots(figsize=(5, 3.6))\n", + "for model, by_size in ALPHA_FD.items():\n", + " ax.plot(list(by_size), list(by_size.values()), \"o-\", label=model)\n", + "ax.axhline(2.6, color=\"gray\", ls=\"--\", label=\"experiment\")\n", + "ax.set_xlabel(\"Supercell size n (n x n x n)\")\n", + "ax.set_ylabel(\"$\\\\alpha$(300 K) (10$^{-6}$ K$^{-1}$)\")\n", + "ax.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "29", + "metadata": {}, + "source": [ + "- From 9x9x9 to 10x10x10 the value changes by less than 0.02e-6 K⁻¹. It is\n", + " converged at 2.90e-6 K⁻¹ with MACE-OMAT-0-medium, 3.71e-6 K⁻¹ with\n", + " MATPES-PBE-0 and 0.58e-6 K⁻¹ with MATPES-r2SCAN-0.\n", + "- We also cut the 8x8x8 finite-displacement force constants at the pheasy\n", + " cutoffs and recomputed the thermal expansion. At 16 and 18 Bohr this agrees\n", + " with the 6x6x6 fits within 0.05e-6 K⁻¹. The fits are right for their cutoff.\n", + "- The value only settles in the 9x9x9 and 10x10x10 supercells. They hold\n", + " interactions up to 17 and 19 Å. Each of these MACE potentials sees the\n", + " neighbors within 12 Å of an atom, through two layers with a 6 Å cutoff. So\n", + " the force constants can couple atoms up to 24 Å apart.\n", + "- A pheasy fit with a cutoff of about 19 Å, or 36 Bohr, is far beyond what the\n", + " workflow can do. The fit at 20 Bohr already ran out of memory.\n", + "- MACE-OMAT-0-medium gives 12% more than experiment and MATPES-PBE-0 gives 43%\n", + " more. MATPES-r2SCAN-0 gives about a fifth of the experimental value.\n", + "\n", + "Si is a hard case, because its thermal expansion is a small difference of two\n", + "larger terms. For such a material, compare the result with finite displacements\n", + "in growing supercells. In section 1 the default supercell is within 2% for the\n", + "rock-salt materials with MACE-OMAT-0-medium. With the other potentials, and for\n", + "GaAs, a larger supercell changes the result by up to 38%." + ] + }, + { + "cell_type": "markdown", + "id": "30", + "metadata": {}, "source": [ "## Known limitations\n", "\n", From bc5f382f3c7a9ecfc4af41ddde4cb95fd7bad51c Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Thu, 1 Oct 2026 17:24:51 -0700 Subject: [PATCH 16/28] Update the run times in the thermal expansion tutorial --- tutorials/cte_workflow.ipynb | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/tutorials/cte_workflow.ipynb b/tutorials/cte_workflow.ipynb index 1bcb8ab8dc..ed51c480d0 100644 --- a/tutorials/cte_workflow.ipynb +++ b/tutorials/cte_workflow.ipynb @@ -255,7 +255,7 @@ "\n", "The flow has a little over 200 jobs, most of them force calculations on the\n", "displaced supercells. On 32 CPU cores of a Perlmutter node the whole run took\n", - "about 11 minutes.\n", + "about 8 minutes.\n", "\n", "`create_folders=True` is needed, because the thermal expansion job reads the\n", "force constants from the folder of the pheasy fit." @@ -376,8 +376,7 @@ "Section 2 uses it as well.\n", "\n", "The loop runs 15 workflows, so it is off by default. Set `RUN_ALL = True` to\n", - "run it. On 32 CPU cores of a Perlmutter node each run took between 1 and 26\n", - "minutes." + "run it." ] }, { From 729ef46a23e88d2cdea0edfd3e5bcad5358eec42 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Fri, 2 Oct 2026 07:29:02 -0700 Subject: [PATCH 17/28] Link the phonon database and pass get_supercell_size_kwargs to the pheasy supercell job --- docs/user/codes/vasp.md | 2 +- src/atomate2/common/flows/pheasy.py | 11 ++++++++--- src/atomate2/common/jobs/pheasy.py | 8 ++++++-- src/atomate2/vasp/flows/pheasy.py | 2 +- tests/common/jobs/test_pheasy.py | 23 +++++++++++++++++++++++ tutorials/pheasy_workflow.ipynb | 2 +- 6 files changed, 40 insertions(+), 8 deletions(-) diff --git a/docs/user/codes/vasp.md b/docs/user/codes/vasp.md index 31561ca004..ec9a822d89 100644 --- a/docs/user/codes/vasp.md +++ b/docs/user/codes/vasp.md @@ -350,7 +350,7 @@ It also installs phono3py for the thermal expansion workflow. 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`. -This workflow was used to build the Materials Project's Harmonic Phonon Database, described in [this preprint](https://chemrxiv.org/doi/full/10.26434/chemrxiv.15004632/v1). +This workflow was used to build the [Materials Project's Harmonic Phonon Database](https://next-gen.materialsproject.org/materials?has_props=phonon), described in [this preprint](https://chemrxiv.org/doi/full/10.26434/chemrxiv.15004632/v1). The pheasy fits can fail for a cell that is not in a standard setting, for example the primitive cell of MgO (mp-1265) as the Materials Project serves it. For such a structure, set `use_symmetrized_structure="primitive"`, as was done for the database. diff --git a/src/atomate2/common/flows/pheasy.py b/src/atomate2/common/flows/pheasy.py index 158c90efab..8e3d5f64a9 100644 --- a/src/atomate2/common/flows/pheasy.py +++ b/src/atomate2/common/flows/pheasy.py @@ -6,7 +6,7 @@ from dataclasses import dataclass, field from typing import TYPE_CHECKING, Literal -from pymatgen.util.due import Doi, due +from pymatgen.util.due import Doi, Url, due from atomate2.common.flows.phonons import BasePhononMaker as PurePhonopyMaker from atomate2.common.jobs.pheasy import ( @@ -33,9 +33,13 @@ @due.dcite( - Doi("10.26434/chemrxiv.15004632/v1"), + Url("https://next-gen.materialsproject.org/materials?has_props=phonon"), description="Materials Project's Harmonic Phonon Database.", ) +@due.dcite( + Doi("10.26434/chemrxiv.15004632/v1"), + description="Preprint on the Materials Project's Harmonic Phonon Database.", +) @due.dcite( Doi("10.48550/arXiv.2508.01020"), description="Pheasy code for (an)harmonic force constants.", @@ -148,7 +152,7 @@ class BasePhononMaker(PurePhonopyMaker, ABC): if set to True, only diagonal supercell matrices are allowed. The pheasy commands take the diagonal of the supercell matrix. get_supercell_size_kwargs: dict - not used by this workflow. + kwargs that will be passed to get_supercell_size to determine supercell size use_symmetrized_structure: str allowed strings: "primitive", "conventional", None @@ -400,6 +404,7 @@ def get_supercell_matrix(self, structure: Structure) -> Job | Flow: self.max_atoms, self.force_90_degrees, self.force_diagonal, + **self.get_supercell_size_kwargs, ) @property diff --git a/src/atomate2/common/jobs/pheasy.py b/src/atomate2/common/jobs/pheasy.py index 18e5d8c2c5..6c69ca1176 100644 --- a/src/atomate2/common/jobs/pheasy.py +++ b/src/atomate2/common/jobs/pheasy.py @@ -362,6 +362,7 @@ def get_supercell_size( max_atoms: int, force_90_degrees: bool, force_diagonal: bool, + **kwargs, ) -> list[list[float]]: """ Determine the supercell matrix with pymatgen's CubicSupercellTransformation. @@ -378,14 +379,17 @@ def get_supercell_size( if True, only supercells with three 90 degree angles are allowed force_diagonal: bool if True, only diagonal supercell matrices are allowed + **kwargs: + Additional parameters passed to CubicSupercellTransformation. They + override the defaults angle_tolerance=1e-2 and allow_orthorhombic=False. """ + kwargs = {"angle_tolerance": 1e-2, "allow_orthorhombic": False} | kwargs transformation = CubicSupercellTransformation( min_length=min_length, max_atoms=max_atoms, force_90_degrees=force_90_degrees, force_diagonal=force_diagonal, - angle_tolerance=1e-2, - allow_orthorhombic=False, + **kwargs, ) transformation.apply_transformation(structure=structure) return transformation.transformation_matrix.transpose().tolist() diff --git a/src/atomate2/vasp/flows/pheasy.py b/src/atomate2/vasp/flows/pheasy.py index 14a4368c4a..d934c538dd 100644 --- a/src/atomate2/vasp/flows/pheasy.py +++ b/src/atomate2/vasp/flows/pheasy.py @@ -123,7 +123,7 @@ class PhononMaker(BasePhononMaker): allow_orthorhombic: bool not used by this workflow. get_supercell_size_kwargs: dict - not used by this workflow. + kwargs that will be passed to get_supercell_size to determine supercell size use_symmetrized_structure: str allowed strings: "primitive", "conventional", None diff --git a/tests/common/jobs/test_pheasy.py b/tests/common/jobs/test_pheasy.py index fc1d816d29..58b72dc7a0 100644 --- a/tests/common/jobs/test_pheasy.py +++ b/tests/common/jobs/test_pheasy.py @@ -29,6 +29,7 @@ _get_num_irreducible_fcs, _run_band_structure_and_plot, ) +from atomate2.forcefields.flows.pheasy import PhononMaker # fcs_cutoff_radius in Bohr. 8 Bohr (4.2 A) covers the first two neighbour # shells of fcc Cu (2.55 and 3.61 A) and stays inside the 10.8 A supercell. @@ -260,6 +261,28 @@ def test_get_num_anharmonic_supercells(monkeypatch): _get_num_anharmonic_supercells(num_disp_anhar=0, **kwargs) +def test_get_supercell_size_kwargs(monkeypatch): + received = {} + transformation = pheasy_jobs.CubicSupercellTransformation + + def record_kwargs(**kwargs): + received.update(kwargs) + return transformation(**kwargs) + + monkeypatch.setattr(pheasy_jobs, "CubicSupercellTransformation", record_kwargs) + + # the maker passes get_supercell_size_kwargs on to the job + maker = PhononMaker(get_supercell_size_kwargs={"angle_tolerance": 0.1}) + job = maker.get_supercell_matrix(_cu_structure()) + assert job.function_kwargs == {"angle_tolerance": 0.1} + + # the job passes them to CubicSupercellTransformation. The other default + # is kept. + job.function(*job.function_args, **job.function_kwargs) + assert received["angle_tolerance"] == 0.1 + assert received["allow_orthorhombic"] is False + + def test_check_lasso_alpha(tmp_dir): log_file = Path("pheasy_anharmonic_fit.log") diff --git a/tutorials/pheasy_workflow.ipynb b/tutorials/pheasy_workflow.ipynb index 6f6742b043..91f318cd21 100644 --- a/tutorials/pheasy_workflow.ipynb +++ b/tutorials/pheasy_workflow.ipynb @@ -42,7 +42,7 @@ "\n", "The [Pheasy code](https://doi.org/10.48550/arXiv.2508.01020) extracts (an)harmonic interatomic force constants from a set of displaced supercell calculations using machine-learning (LASSO) regression. Compared to the finite-displacement approach in the standard phonon workflow, far fewer (randomly displaced) supercells are needed, which substantially reduces the number of DFT calculations.\n", "\n", - "This workflow was used to build the Materials Project's Harmonic Phonon Database, described in [this preprint](https://chemrxiv.org/doi/full/10.26434/chemrxiv.15004632/v1).\n", + "This workflow was used to build the [Materials Project's Harmonic Phonon Database](https://next-gen.materialsproject.org/materials?has_props=phonon), described in [this preprint](https://chemrxiv.org/doi/full/10.26434/chemrxiv.15004632/v1).\n", "\n", "The workflow has the same basic structure as the phonopy-based phonon workflow and uses [Phonopy](https://doi.org/10.7566/JPSJ.92.012001) to post-process the force constants into band structures, densities of states, and thermodynamic properties.\n", "\n", From ce11cb8a840dc0ab90c7f8ef85ebe04e30d01149 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Fri, 2 Oct 2026 08:41:34 -0700 Subject: [PATCH 18/28] Note in the docs that the anharmonic pheasy workflow is less tested --- docs/user/codes/vasp.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/user/codes/vasp.md b/docs/user/codes/vasp.md index ec9a822d89..a481fbf326 100644 --- a/docs/user/codes/vasp.md +++ b/docs/user/codes/vasp.md @@ -355,6 +355,12 @@ The pheasy fits can fail for a cell that is not in a standard setting, for examp For such a structure, set `use_symmetrized_structure="primitive"`, as was done for the database. By default, this workflow does not compute anharmonic force constants, but can be extended to using the `cal_anhar_fcs` kwarg. + +```{warning} +The anharmonic part of this workflow has not been tested as widely as the harmonic part. +It might still change in future versions. +``` + ALM, from the ALAMODE package, counts the free force constants used to size the random displacement sets. To install ALAMODE, see their [installation guidelines](https://alamode.readthedocs.io/en/latest/install.html#). From 0ac41b1275b2ac5a44e9b3f7164c97002479a17f Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Fri, 2 Oct 2026 08:56:45 -0700 Subject: [PATCH 19/28] Name the volumetric CTE thermal_expansion as in the QHA document --- src/atomate2/common/schemas/cte.py | 8 ++++---- tests/common/jobs/test_cte.py | 20 ++++++++++---------- tutorials/cte_workflow.ipynb | 4 ++-- 3 files changed, 16 insertions(+), 16 deletions(-) diff --git a/src/atomate2/common/schemas/cte.py b/src/atomate2/common/schemas/cte.py index a071615926..a781dcc408 100644 --- a/src/atomate2/common/schemas/cte.py +++ b/src/atomate2/common/schemas/cte.py @@ -146,12 +146,12 @@ class CTEResult(BaseModel): description="Whether a frequency on the sampling mesh lies below " "-tol_imaginary_modes. The thermal expansion is not computed in that case." ) - thermal_expansion: list[Matrix3D] | None = Field( + thermal_expansion_tensor: list[Matrix3D] | None = Field( None, description="Thermal expansion tensor in 1/K at each temperature, in the " "Cartesian frame of the structure.", ) - volumetric_thermal_expansion: list[float] | None = Field( + thermal_expansion: list[float] | None = Field( None, description="Volumetric thermal expansion in 1/K at each temperature, the " "trace of the thermal expansion tensor.", @@ -346,8 +346,8 @@ def from_force_constants( temperatures, min_frequency=min_frequency, ) - result["thermal_expansion"] = alphas.tolist() - result["volumetric_thermal_expansion"] = np.trace( + result["thermal_expansion_tensor"] = alphas.tolist() + result["thermal_expansion"] = np.trace( alphas, axis1=1, axis2=2 ).tolist() result["average_gruneisen"] = [ diff --git a/tests/common/jobs/test_cte.py b/tests/common/jobs/test_cte.py index 8dfddad6eb..4282a71348 100644 --- a/tests/common/jobs/test_cte.py +++ b/tests/common/jobs/test_cte.py @@ -176,21 +176,21 @@ def test_cte_maker_emt(clean_dir, monkeypatch): for result in doc.results: assert not result.has_imaginary_modes - alpha = np.array(result.thermal_expansion) + alpha = np.array(result.thermal_expansion_tensor) assert np.all(alpha[0] == 0.0) # cubic, so alpha is isotropic in every frame assert alpha[2] == pytest.approx(alpha[2, 0, 0] * np.eye(3), abs=1e-12) - assert result.volumetric_thermal_expansion[2] == pytest.approx( - 3 * alpha[2, 0, 0] - ) + assert result.thermal_expansion[2] == pytest.approx(3 * alpha[2, 0, 0]) # linear thermal expansion at 300 K in 1/K, and the mean Grueneisen parameter. # The one-shot fit has few supercells in this small cell, so its value is only # compared with the cocktail value. cocktail, one_shot = doc.results - assert cocktail.thermal_expansion[2][0][0] == pytest.approx(1.859e-5, rel=0.02) + assert cocktail.thermal_expansion_tensor[2][0][0] == pytest.approx( + 1.859e-5, rel=0.02 + ) assert np.trace(cocktail.average_gruneisen[2]) / 3 == pytest.approx(2.237, rel=0.02) - assert one_shot.thermal_expansion[2][0][0] == pytest.approx( - cocktail.thermal_expansion[2][0][0], rel=0.2 + assert one_shot.thermal_expansion_tensor[2][0][0] == pytest.approx( + cocktail.thermal_expansion_tensor[2][0][0], rel=0.2 ) assert Path(doc.phonon_job_dir, "one_shot", "fc3.hdf5").exists() @@ -212,7 +212,7 @@ def test_cte_maker_emt(clean_dir, monkeypatch): assert flagged_doc.mesh == (2, 2, 2) (flagged,) = flagged_doc.results assert flagged.has_imaginary_modes - assert flagged.thermal_expansion is None + assert flagged.thermal_expansion_tensor is None # the non-analytical term correction with zero Born charges leaves alpha # unchanged. pheasy only stores Born charges for VASP, so they are added here. @@ -254,8 +254,8 @@ def _gruneisen(*args, **kwargs): (nac_params,) = nac_params_used assert nac_params["born"].shape == (4, 3, 3) assert nac_params["dielectric"] == pytest.approx(np.eye(3) * 10.0) - assert np.array(nac_doc.results[0].thermal_expansion[0]) == pytest.approx( - np.array(cocktail.thermal_expansion[2]), rel=1e-6, abs=1e-15 + assert np.array(nac_doc.results[0].thermal_expansion_tensor[0]) == pytest.approx( + np.array(cocktail.thermal_expansion_tensor[2]), rel=1e-6, abs=1e-15 ) # a structure with another lattice or other atoms is refused diff --git a/tutorials/cte_workflow.ipynb b/tutorials/cte_workflow.ipynb index ed51c480d0..d2b6b13f99 100644 --- a/tutorials/cte_workflow.ipynb +++ b/tutorials/cte_workflow.ipynb @@ -301,7 +301,7 @@ "\n", "(result,) = doc.results\n", "temperatures = np.array(doc.temperatures)\n", - "alpha = np.array(result.thermal_expansion)\n", + "alpha = np.array(result.thermal_expansion_tensor)\n", "i300 = int(np.argmin(np.abs(temperatures - 300)))\n", "\n", "fig, ax = plt.subplots()\n", @@ -434,7 +434,7 @@ " \"\"\"Linear thermal expansion at 300 K in 1/K, from the xx component.\"\"\"\n", " (result,) = doc.results\n", " i300 = int(np.argmin(np.abs(np.array(doc.temperatures) - 300)))\n", - " return float(np.array(result.thermal_expansion)[i300, 0, 0])\n", + " return float(np.array(result.thermal_expansion_tensor)[i300, 0, 0])\n", "\n", "\n", "RUN_ALL = False\n", From a754a6aecf8be824a62b23c8592f544a915caf93 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Fri, 2 Oct 2026 09:01:13 -0700 Subject: [PATCH 20/28] Store the supercell matrix in the CTE document --- src/atomate2/common/schemas/cte.py | 4 ++++ tests/common/jobs/test_cte.py | 1 + 2 files changed, 5 insertions(+) diff --git a/src/atomate2/common/schemas/cte.py b/src/atomate2/common/schemas/cte.py index a781dcc408..dcb0a1c080 100644 --- a/src/atomate2/common/schemas/cte.py +++ b/src/atomate2/common/schemas/cte.py @@ -169,6 +169,9 @@ class CTEDocument(StructureMetadata): structure: Structure | None = Field( None, description="Structure used for the phonon and elastic calculations." ) + supercell_matrix: Matrix3D | None = Field( + None, description="Supercell matrix of the phonon calculation." + ) temperatures: list[float] | None = Field(None, description="Temperatures in K.") mesh: tuple[int, int, int] | None = Field( None, description="q-point mesh used for the mode Grueneisen tensors." @@ -358,6 +361,7 @@ def from_force_constants( return cls.from_structure( meta_structure=structure, structure=structure, + supercell_matrix=phonon.supercell_matrix.tolist(), temperatures=list(temperatures), mesh=mesh_numbers, elastic_tensor=np.asarray(elastic_tensor).tolist(), diff --git a/tests/common/jobs/test_cte.py b/tests/common/jobs/test_cte.py index 4282a71348..675a52836d 100644 --- a/tests/common/jobs/test_cte.py +++ b/tests/common/jobs/test_cte.py @@ -167,6 +167,7 @@ def test_cte_maker_emt(clean_dir, monkeypatch): doc = responses[flow.output.uuid][1].output assert isinstance(doc, CTEDocument) assert doc.mesh == (8, 8, 8) + assert np.array(doc.supercell_matrix) == pytest.approx(2 * np.eye(3)) assert [result.fit_method for result in doc.results] == ["cocktail", "one-shot"] # EMT elastic constants of Cu in GPa From fc5329319a42267bd437df8eb8ad7a063fbef21c Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Fri, 2 Oct 2026 09:21:43 -0700 Subject: [PATCH 21/28] Move phono3py into its own extra --- .github/workflows/testing.yml | 2 +- docs/user/codes/vasp.md | 4 ++-- pyproject.toml | 3 ++- tutorials/cte_workflow.ipynb | 12 ++++++------ 4 files changed, 11 insertions(+), 10 deletions(-) diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index e34e9e1c01..7e3db5984b 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -147,7 +147,7 @@ jobs: python -m pip install --upgrade pip mkdir -p ~/.abinit/pseudos cp -r tests/test_data/abinit/pseudos/ONCVPSP-PBE-SR-PDv0.4 ~/.abinit/pseudos - uv pip install .[strict,abinit,approxneb,aims,pheasy,alamode,hiphive] --group tests + uv pip install .[strict,abinit,approxneb,aims,pheasy,alamode,hiphive,phono3py] --group tests uv pip install torch-runstats torch_dftd uv pip install --no-deps nequip==0.5.6 diff --git a/docs/user/codes/vasp.md b/docs/user/codes/vasp.md index a481fbf326..d05c7eb0e0 100644 --- a/docs/user/codes/vasp.md +++ b/docs/user/codes/vasp.md @@ -346,7 +346,6 @@ phonon_flow = PhononMaker(min_length=15.0, store_force_constants=False).make( 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. -It also installs phono3py for the thermal expansion workflow. 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`. @@ -441,7 +440,8 @@ gruneisen_flow = GruneisenMaker( ### Thermal expansion workflow `CTEMaker` calculates the thermal expansion tensor from third-order force constants, with the help of [Pheasy](https://doi.org/10.48550/arXiv.2508.01020) and [phono3py](https://doi.org/10.1088/1361-648X/acd831). -It needs the `pheasy` extra, see the Pheasy section above. +It needs the `pheasy` and `phono3py` extras, `pip install "atomate2[pheasy,phono3py]"`. +For the pheasy extra, see the Pheasy section above. First, the structure is converted to the standard primitive cell, and a tight structural relaxation is performed. The pheasy fits can fail for cells that are not in a standard setting. diff --git a/pyproject.toml b/pyproject.toml index 979b8bbde6..8d9b27680c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -54,9 +54,10 @@ mp = ["mp-api>=0.37.5"] phonons = ["phonopy>=2.43,<5", "seekpath>=2.0.0"] # the pheasy fork adds the --disp_matrix_file and --force_matrix_file options # and fixes the --scell option +pheasy = ["atomate2[phonons,alamode]", "hiphive==1.5", "numpy<=2.2", "pheasy @ git+https://gitlab.com/hpsahasrabuddhe/pheasy.git@9f24162a4ed0f0ab8911d382fd2617944aade55d"] # phono3py computes the mode Grueneisen tensors in the thermal expansion workflow # phonors 0.5 breaks the Grueneisen q-point mesh of phonopy and phono3py 4.5 -pheasy = ["atomate2[phonons,alamode]", "hiphive==1.5", "numpy<=2.2", "pheasy @ git+https://gitlab.com/hpsahasrabuddhe/pheasy.git@9f24162a4ed0f0ab8911d382fd2617944aade55d", "phono3py==4.5.0", "phonors<0.5"] +phono3py = ["atomate2[phonons]", "phono3py==4.5.0", "phonors<0.5"] alamode = ["alm @ git+https://github.com/ttadano/ALM.git@f1d668fdee66e7e7218a04c88daf19d0e14fce0c#subdirectory=python"] hiphive = ["hiphive==1.5", "trainstation>=1.0", "atomate2[phonons,alamode]"] lobster = ["ijson>=3.2.2", "lobsterpy>=0.6.0"] diff --git a/tutorials/cte_workflow.ipynb b/tutorials/cte_workflow.ipynb index d2b6b13f99..92572b3b9f 100644 --- a/tutorials/cte_workflow.ipynb +++ b/tutorials/cte_workflow.ipynb @@ -54,17 +54,17 @@ "source": [ "## Installation\n", "\n", - "The `pheasy` and `ase` extras are needed, plus MACE and the Materials Project\n", - "client:\n", + "The `pheasy`, `phono3py` and `ase` extras are needed, plus MACE and the\n", + "Materials Project client:\n", "\n", "```\n", - "pip install 'atomate2[pheasy,ase]'\n", + "pip install 'atomate2[pheasy,phono3py,ase]'\n", "pip install 'mace-torch>=0.3.16' mp-api\n", "```\n", "\n", - "The `pheasy` extra pulls in pheasy, phono3py and ALM. ALM is compiled from\n", - "source. The hiPhive tutorial and the atomate2 VASP documentation list what the\n", - "build needs." + "The `pheasy` extra pulls in pheasy and ALM, and the `phono3py` extra pulls in\n", + "phono3py. ALM is compiled from source. The hiPhive tutorial and the atomate2\n", + "VASP documentation list what the build needs." ] }, { From cb91f951a51f521ad3823a2839d0fb1730dfbf2a Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Fri, 2 Oct 2026 09:28:44 -0700 Subject: [PATCH 22/28] Note that the thermal expansion workflow needs more testing --- docs/user/codes/vasp.md | 5 +++++ src/atomate2/common/flows/cte.py | 3 +++ src/atomate2/forcefields/flows/cte.py | 3 +++ src/atomate2/vasp/flows/cte.py | 3 +++ 4 files changed, 14 insertions(+) diff --git a/docs/user/codes/vasp.md b/docs/user/codes/vasp.md index d05c7eb0e0..f283cd436d 100644 --- a/docs/user/codes/vasp.md +++ b/docs/user/codes/vasp.md @@ -443,6 +443,11 @@ gruneisen_flow = GruneisenMaker( It needs the `pheasy` and `phono3py` extras, `pip install "atomate2[pheasy,phono3py]"`. For the pheasy extra, see the Pheasy section above. +```{warning} +This workflow is new and has not been tested widely. +It might still change in future versions. +``` + First, the structure is converted to the standard primitive cell, and a tight structural relaxation is performed. The pheasy fits can fail for cells that are not in a standard setting. Set `use_symmetrized_structure="conventional"` to use the standard conventional cell instead. diff --git a/src/atomate2/common/flows/cte.py b/src/atomate2/common/flows/cte.py index 95ec25a1e8..9868a76f33 100644 --- a/src/atomate2/common/flows/cte.py +++ b/src/atomate2/common/flows/cte.py @@ -49,6 +49,9 @@ class BaseCTEMaker(Maker, ABC): constants at the relaxed structure. The frequencies are not renormalized with temperature. + This workflow is new and has not been tested widely. It might still change + in future versions. + Parameters ---------- name: str diff --git a/src/atomate2/forcefields/flows/cte.py b/src/atomate2/forcefields/flows/cte.py index e2270e8bde..d083ddf463 100644 --- a/src/atomate2/forcefields/flows/cte.py +++ b/src/atomate2/forcefields/flows/cte.py @@ -73,6 +73,9 @@ class CTEMaker(BaseCTEMaker): constants with the one-shot method from randomly displaced supercells with 0.03 A displacements, and builds the supercells with min_length=12.0. + This workflow is new and has not been tested widely. It might still change + in future versions. + Parameters ---------- name: str diff --git a/src/atomate2/vasp/flows/cte.py b/src/atomate2/vasp/flows/cte.py index dc7d691b0e..42be710e18 100644 --- a/src/atomate2/vasp/flows/cte.py +++ b/src/atomate2/vasp/flows/cte.py @@ -36,6 +36,9 @@ class CTEMaker(BaseCTEMaker): relaxation. The stress is more sensitive to ENCUT than the forces are, so check the ENCUT convergence of the elastic tensor for your material. + This workflow is new and has not been tested widely. It might still change + in future versions. + Parameters ---------- name: str From ff0582a59fbb933582194d3d8a9fdb474b278933 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Fri, 2 Oct 2026 10:35:13 -0700 Subject: [PATCH 23/28] Run the non-ase tests with -v so each test shows in the CI log --- .github/workflows/testing.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index 7e3db5984b..668d2d0c29 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -168,7 +168,7 @@ jobs: # However this `splitting-algorithm` means that tests cannot depend sensitively on the order they're executed in. run: | micromamba activate a2 - pytest --durations=5 -n auto --splits 3 --group ${{ matrix.split }} --durations-path tests/.pytest-split-durations --splitting-algorithm least_duration --ignore=tests/ase --ignore=tests/openff_md --ignore=tests/openmm_md --ignore=tests/forcefields --ignore=tests/torchsim --cov=atomate2 --cov-report=xml + pytest -v --durations=5 -n auto --splits 3 --group ${{ matrix.split }} --durations-path tests/.pytest-split-durations --splitting-algorithm least_duration --ignore=tests/ase --ignore=tests/openff_md --ignore=tests/openmm_md --ignore=tests/forcefields --ignore=tests/torchsim --cov=atomate2 --cov-report=xml - uses: codecov/codecov-action@v1 From dcde48ddbf17ecf94f77aeb0f40fa7a66bfbbed2 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Fri, 2 Oct 2026 10:51:11 -0700 Subject: [PATCH 24/28] Run the forcefield tests with -v as well --- .github/workflows/testing.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index 668d2d0c29..745a7dd349 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -85,7 +85,7 @@ jobs: # xdist spawns multiple workers each holding a copy of the models. run: | micromamba activate a2 - pytest --durations=5 -n ${{ matrix.dep-group == 'torch-limited' && '1' || 'auto' }} --cov=atomate2 --cov-report=xml tests/forcefields + pytest -v --durations=5 -n ${{ matrix.dep-group == 'torch-limited' && '1' || 'auto' }} --cov=atomate2 --cov-report=xml tests/forcefields - name: Forcefield tutorial if: matrix.dep-group == 'torch-limited' @@ -93,7 +93,7 @@ jobs: MP_API_KEY: ${{ secrets.MP_API_KEY }} run: | micromamba activate a2 - pytest --durations=5 -n auto --nbmake ./tutorials/force_fields + pytest -v --durations=5 -n auto --nbmake ./tutorials/force_fields - uses: codecov/codecov-action@v1 if: matrix.python-version == '3.11' && github.repository == 'materialsproject/atomate2' From 8de4e043c4b34960f300fc71934a4eebf1430586 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Fri, 2 Oct 2026 11:23:15 -0700 Subject: [PATCH 25/28] Move the force field CTE tests to tests/forcefields --- .github/workflows/testing.yml | 8 ++ tests/common/jobs/test_cte.py | 165 +------------------------- tests/forcefields/flows/test_cte.py | 172 ++++++++++++++++++++++++++++ 3 files changed, 183 insertions(+), 162 deletions(-) create mode 100644 tests/forcefields/flows/test_cte.py diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index 745a7dd349..7405467740 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -70,6 +70,14 @@ jobs: cp -r tests/test_data/abinit/pseudos/ONCVPSP-PBE-SR-PDv0.4 ~/.abinit/pseudos uv pip install .[strict,strict-forcefields-${{ matrix.dep-group }},abinit,approxneb,aims] --group tests + # the thermal expansion tests use EMT, so one group is enough to run them + - name: Install pheasy, ALM and phono3py for the thermal expansion tests + if: matrix.dep-group == 'numpy-limited' + run: | + micromamba install -n a2 -c conda-forge compilers boost eigen=3.3 cmake spglib --yes + micromamba activate a2 + uv pip install .[strict,strict-forcefields-${{ matrix.dep-group }},abinit,approxneb,aims,pheasy,phono3py] --group tests + - name: Install pymatgen from master if triggered by pymatgen repo dispatch if: github.event_name == 'repository_dispatch' && github.event.action == 'pymatgen-ci-trigger' run: | diff --git a/tests/common/jobs/test_cte.py b/tests/common/jobs/test_cte.py index 675a52836d..f2173216b5 100644 --- a/tests/common/jobs/test_cte.py +++ b/tests/common/jobs/test_cte.py @@ -1,30 +1,18 @@ -"""Tests for the thermal expansion workflow. +"""Tests for get_cte and the Born charge mapping of the CTE document. -The end-to-end test runs the force field CTEMaker with ASE's EMT potential, so -it needs no DFT reference data. It lives here and not under tests/forcefields, -because pheasy and ALM are only installed in the test-non-ase CI job. +The force field workflow is tested in tests/forcefields/flows/test_cte.py. """ -from pathlib import Path - import numpy as np -import phono3py.phonon3.gruneisen as gruneisen_module -import phonopy import pytest from ase.build import bulk -from jobflow import run_locally from phonopy import Phonopy from phonopy.structure.atoms import PhonopyAtoms from pymatgen.analysis.elasticity import ElasticTensor -from pymatgen.io.ase import AseAtomsAdaptor from scipy.constants import Boltzmann, Planck from scipy.spatial.transform import Rotation -from atomate2.common.jobs.cte import compute_cte -from atomate2.common.schemas.cte import CTEDocument, _expand_born_to_unitcell, get_cte -from atomate2.forcefields.flows.cte import CTEMaker - -EMT = {"@module": "ase.calculators.emt", "@callable": "EMT"} +from atomate2.common.schemas.cte import _expand_born_to_unitcell, get_cte def _heat_capacity(frequency: float, temperature: float) -> float: @@ -142,150 +130,3 @@ def test_expand_born_to_unitcell(): assert born.shape == (8, 3, 3) for symbol, charges in zip(unitcell.symbols, born, strict=True): assert charges == pytest.approx(np.eye(3) * born_primitive[symbol]) - - -def test_cte_maker_emt(clean_dir, monkeypatch): - """Run the whole force field workflow with EMT forces on fcc Cu.""" - structure = AseAtomsAdaptor.get_structure(bulk("Cu", "fcc", a=3.61, cubic=True)) - # the conventional cell keeps the 4-atom cubic unit cell used below - maker = CTEMaker.from_force_field_name( - EMT, - use_symmetrized_structure="conventional", - temperatures=[0, 100, 300], - mesh=(8, 8, 8), - ) - # a 2x2x2 supercell of the cubic cell, 32 atoms, and an fc3 cutoff of 6 Bohr - # (3.2 A), which covers the nearest neighbours at 2.55 A - maker.phonon_maker.min_length = 7.0 - maker.phonon_maker.fcs_cutoff_radius = [-1, 6, 6] - maker.phonon_maker.anhar_fit_methods = ["cocktail", "one-shot"] - - flow = maker.make(structure) - # the uuid of the phonon flow output that compute_cte reads - phonon_uuid = flow.jobs[-1].function_kwargs["phonon_output"].uuid - responses = run_locally(flow, create_folders=True, ensure_success=True) - doc = responses[flow.output.uuid][1].output - assert isinstance(doc, CTEDocument) - assert doc.mesh == (8, 8, 8) - assert np.array(doc.supercell_matrix) == pytest.approx(2 * np.eye(3)) - assert [result.fit_method for result in doc.results] == ["cocktail", "one-shot"] - - # EMT elastic constants of Cu in GPa - assert doc.elastic_tensor[0][0] == pytest.approx(172.6, rel=0.01) - assert doc.elastic_tensor[0][1] == pytest.approx(115.4, rel=0.01) - assert doc.elastic_tensor[3][3] == pytest.approx(89.9, rel=0.01) - - for result in doc.results: - assert not result.has_imaginary_modes - alpha = np.array(result.thermal_expansion_tensor) - assert np.all(alpha[0] == 0.0) - # cubic, so alpha is isotropic in every frame - assert alpha[2] == pytest.approx(alpha[2, 0, 0] * np.eye(3), abs=1e-12) - assert result.thermal_expansion[2] == pytest.approx(3 * alpha[2, 0, 0]) - # linear thermal expansion at 300 K in 1/K, and the mean Grueneisen parameter. - # The one-shot fit has few supercells in this small cell, so its value is only - # compared with the cocktail value. - cocktail, one_shot = doc.results - assert cocktail.thermal_expansion_tensor[2][0][0] == pytest.approx( - 1.859e-5, rel=0.02 - ) - assert np.trace(cocktail.average_gruneisen[2]) / 3 == pytest.approx(2.237, rel=0.02) - assert one_shot.thermal_expansion_tensor[2][0][0] == pytest.approx( - cocktail.thermal_expansion_tensor[2][0][0], rel=0.2 - ) - assert Path(doc.phonon_job_dir, "one_shot", "fc3.hdf5").exists() - - # the same force constants, with a q-point density instead of a mesh. The - # negative tolerance flags every frequency, to test the imaginary mode check. - phonon_output = responses[phonon_uuid][1].output - job = compute_cte( - phonon_output=phonon_output, - elastic_tensor=doc.elastic_tensor, - elastic_structure=doc.structure, - anhar_fit_methods=["one-shot"], - temperatures=[300], - mesh=100.0, - tol_imaginary_modes=-10.0, - ) - with pytest.warns(UserWarning, match="thermal expansion is not"): - responses = run_locally(job, create_folders=True, ensure_success=True) - flagged_doc = responses[job.uuid][1].output - assert flagged_doc.mesh == (2, 2, 2) - (flagged,) = flagged_doc.results - assert flagged.has_imaginary_modes - assert flagged.thermal_expansion_tensor is None - - # the non-analytical term correction with zero Born charges leaves alpha - # unchanged. pheasy only stores Born charges for VASP, so they are added here. - # The phonopy primitive cell has one atom and the unit cell four, so the - # charges must be expanded to four atoms before they reach phono3py. - phonon_job_dir = Path(doc.phonon_job_dir) - phonon = phonopy.load(phonon_job_dir / "phonopy.yaml", produce_fc=False) - assert len(phonon.primitive) == 1 - phonon.nac_params = { - "born": np.zeros((1, 3, 3)), - "dielectric": np.eye(3) * 10.0, - "factor": 14.399652, - } - phonon.save("phonopy_nac.yaml") - - nac_params_used = [] - original_gruneisen = gruneisen_module.Gruneisen - - def _gruneisen(*args, **kwargs): - nac_params_used.append(kwargs["nac_params"]) - return original_gruneisen(*args, **kwargs) - - monkeypatch.setattr(gruneisen_module, "Gruneisen", _gruneisen) - cocktail_files = { - "cocktail": (phonon_job_dir / "FORCE_CONSTANTS", phonon_job_dir / "fc3.hdf5") - } - settings = { - "force_constant_files": cocktail_files, - "elastic_tensor": doc.elastic_tensor, - "temperatures": [300], - "mesh": (8, 8, 8), - "tol_imaginary_modes": 0.1, - "min_frequency": 1e-3, - "symprec": 1e-5, - } - nac_doc = CTEDocument.from_force_constants( - phonopy_yaml="phonopy_nac.yaml", structure=doc.structure, **settings - ) - (nac_params,) = nac_params_used - assert nac_params["born"].shape == (4, 3, 3) - assert nac_params["dielectric"] == pytest.approx(np.eye(3) * 10.0) - assert np.array(nac_doc.results[0].thermal_expansion_tensor[0]) == pytest.approx( - np.array(cocktail.thermal_expansion_tensor[2]), rel=1e-6, abs=1e-15 - ) - - # a structure with another lattice or other atoms is refused - strained = doc.structure.copy() - strained.apply_strain(0.01) - with pytest.raises(ValueError, match="not in the same frame"): - CTEDocument.from_force_constants( - phonopy_yaml=phonon_job_dir / "phonopy.yaml", structure=strained, **settings - ) - other_atoms = doc.structure.copy() - other_atoms.replace_species({"Cu": "Au"}) - with pytest.raises(ValueError, match="atoms of the elastic calculation"): - CTEDocument.from_force_constants( - phonopy_yaml=phonon_job_dir / "phonopy.yaml", - structure=other_atoms, - **settings, - ) - with pytest.raises(ValueError, match="must not be negative"): - CTEDocument.from_force_constants( - phonopy_yaml=phonon_job_dir / "phonopy.yaml", - structure=doc.structure, - **{**settings, "temperatures": [-10, 300]}, - ) - - -def test_cte_maker_force_field_defaults(): - """The default phonon maker uses min_length=12.0 for the supercells.""" - maker = CTEMaker() - assert maker.use_symmetrized_structure == "primitive" - assert maker.phonon_maker.min_length == 12.0 - assert maker.phonon_maker.cal_anhar_fcs - assert maker.phonon_maker.displacement_anhar == 0.03 diff --git a/tests/forcefields/flows/test_cte.py b/tests/forcefields/flows/test_cte.py new file mode 100644 index 0000000000..752dbd7117 --- /dev/null +++ b/tests/forcefields/flows/test_cte.py @@ -0,0 +1,172 @@ +"""Tests for the force field thermal expansion workflow. + +The tests run the force field CTEMaker with ASE's EMT potential, so they need no +DFT reference data and no machine-learned force field. pheasy, ALM and phono3py +are only installed in the numpy-limited forcefield CI job, so the tests are +skipped in the other forcefield jobs. +""" + +from pathlib import Path + +import numpy as np +import phonopy +import pytest +from ase.build import bulk +from jobflow import run_locally +from pymatgen.io.ase import AseAtomsAdaptor + +from atomate2.common.jobs.cte import compute_cte +from atomate2.common.schemas.cte import CTEDocument +from atomate2.forcefields.flows.cte import CTEMaker + +pytest.importorskip("pheasy") +gruneisen_module = pytest.importorskip("phono3py.phonon3.gruneisen") + +EMT = {"@module": "ase.calculators.emt", "@callable": "EMT"} + + +def test_cte_maker_emt(clean_dir, monkeypatch): + """Run the whole force field workflow with EMT forces on fcc Cu.""" + structure = AseAtomsAdaptor.get_structure(bulk("Cu", "fcc", a=3.61, cubic=True)) + # the conventional cell keeps the 4-atom cubic unit cell used below + maker = CTEMaker.from_force_field_name( + EMT, + use_symmetrized_structure="conventional", + temperatures=[0, 100, 300], + mesh=(8, 8, 8), + ) + # a 2x2x2 supercell of the cubic cell, 32 atoms, and an fc3 cutoff of 6 Bohr + # (3.2 A), which covers the nearest neighbours at 2.55 A + maker.phonon_maker.min_length = 7.0 + maker.phonon_maker.fcs_cutoff_radius = [-1, 6, 6] + maker.phonon_maker.anhar_fit_methods = ["cocktail", "one-shot"] + + flow = maker.make(structure) + # the uuid of the phonon flow output that compute_cte reads + phonon_uuid = flow.jobs[-1].function_kwargs["phonon_output"].uuid + responses = run_locally(flow, create_folders=True, ensure_success=True) + doc = responses[flow.output.uuid][1].output + assert isinstance(doc, CTEDocument) + assert doc.mesh == (8, 8, 8) + assert np.array(doc.supercell_matrix) == pytest.approx(2 * np.eye(3)) + assert [result.fit_method for result in doc.results] == ["cocktail", "one-shot"] + + # EMT elastic constants of Cu in GPa + assert doc.elastic_tensor[0][0] == pytest.approx(172.6, rel=0.01) + assert doc.elastic_tensor[0][1] == pytest.approx(115.4, rel=0.01) + assert doc.elastic_tensor[3][3] == pytest.approx(89.9, rel=0.01) + + for result in doc.results: + assert not result.has_imaginary_modes + alpha = np.array(result.thermal_expansion_tensor) + assert np.all(alpha[0] == 0.0) + # cubic, so alpha is isotropic in every frame + assert alpha[2] == pytest.approx(alpha[2, 0, 0] * np.eye(3), abs=1e-12) + assert result.thermal_expansion[2] == pytest.approx(3 * alpha[2, 0, 0]) + # linear thermal expansion at 300 K in 1/K, and the mean Grueneisen parameter. + # The one-shot fit has few supercells in this small cell, so its value is only + # compared with the cocktail value. + cocktail, one_shot = doc.results + assert cocktail.thermal_expansion_tensor[2][0][0] == pytest.approx( + 1.859e-5, rel=0.02 + ) + assert np.trace(cocktail.average_gruneisen[2]) / 3 == pytest.approx(2.237, rel=0.02) + assert one_shot.thermal_expansion_tensor[2][0][0] == pytest.approx( + cocktail.thermal_expansion_tensor[2][0][0], rel=0.2 + ) + assert Path(doc.phonon_job_dir, "one_shot", "fc3.hdf5").exists() + + # the same force constants, with a q-point density instead of a mesh. The + # negative tolerance flags every frequency, to test the imaginary mode check. + phonon_output = responses[phonon_uuid][1].output + job = compute_cte( + phonon_output=phonon_output, + elastic_tensor=doc.elastic_tensor, + elastic_structure=doc.structure, + anhar_fit_methods=["one-shot"], + temperatures=[300], + mesh=100.0, + tol_imaginary_modes=-10.0, + ) + with pytest.warns(UserWarning, match="thermal expansion is not"): + responses = run_locally(job, create_folders=True, ensure_success=True) + flagged_doc = responses[job.uuid][1].output + assert flagged_doc.mesh == (2, 2, 2) + (flagged,) = flagged_doc.results + assert flagged.has_imaginary_modes + assert flagged.thermal_expansion_tensor is None + + # the non-analytical term correction with zero Born charges leaves alpha + # unchanged. pheasy only stores Born charges for VASP, so they are added here. + # The phonopy primitive cell has one atom and the unit cell four, so the + # charges must be expanded to four atoms before they reach phono3py. + phonon_job_dir = Path(doc.phonon_job_dir) + phonon = phonopy.load(phonon_job_dir / "phonopy.yaml", produce_fc=False) + assert len(phonon.primitive) == 1 + phonon.nac_params = { + "born": np.zeros((1, 3, 3)), + "dielectric": np.eye(3) * 10.0, + "factor": 14.399652, + } + phonon.save("phonopy_nac.yaml") + + nac_params_used = [] + original_gruneisen = gruneisen_module.Gruneisen + + def _gruneisen(*args, **kwargs): + nac_params_used.append(kwargs["nac_params"]) + return original_gruneisen(*args, **kwargs) + + monkeypatch.setattr(gruneisen_module, "Gruneisen", _gruneisen) + cocktail_files = { + "cocktail": (phonon_job_dir / "FORCE_CONSTANTS", phonon_job_dir / "fc3.hdf5") + } + settings = { + "force_constant_files": cocktail_files, + "elastic_tensor": doc.elastic_tensor, + "temperatures": [300], + "mesh": (8, 8, 8), + "tol_imaginary_modes": 0.1, + "min_frequency": 1e-3, + "symprec": 1e-5, + } + nac_doc = CTEDocument.from_force_constants( + phonopy_yaml="phonopy_nac.yaml", structure=doc.structure, **settings + ) + (nac_params,) = nac_params_used + assert nac_params["born"].shape == (4, 3, 3) + assert nac_params["dielectric"] == pytest.approx(np.eye(3) * 10.0) + assert np.array(nac_doc.results[0].thermal_expansion_tensor[0]) == pytest.approx( + np.array(cocktail.thermal_expansion_tensor[2]), rel=1e-6, abs=1e-15 + ) + + # a structure with another lattice or other atoms is refused + strained = doc.structure.copy() + strained.apply_strain(0.01) + with pytest.raises(ValueError, match="not in the same frame"): + CTEDocument.from_force_constants( + phonopy_yaml=phonon_job_dir / "phonopy.yaml", structure=strained, **settings + ) + other_atoms = doc.structure.copy() + other_atoms.replace_species({"Cu": "Au"}) + with pytest.raises(ValueError, match="atoms of the elastic calculation"): + CTEDocument.from_force_constants( + phonopy_yaml=phonon_job_dir / "phonopy.yaml", + structure=other_atoms, + **settings, + ) + with pytest.raises(ValueError, match="must not be negative"): + CTEDocument.from_force_constants( + phonopy_yaml=phonon_job_dir / "phonopy.yaml", + structure=doc.structure, + **{**settings, "temperatures": [-10, 300]}, + ) + + +def test_cte_maker_force_field_defaults(): + """The default phonon maker uses min_length=12.0 for the supercells.""" + maker = CTEMaker() + assert maker.use_symmetrized_structure == "primitive" + assert maker.phonon_maker.min_length == 12.0 + assert maker.phonon_maker.cal_anhar_fcs + assert maker.phonon_maker.displacement_anhar == 0.03 From d46d2c9e64dbf88a4ff6cd40ce10bd7847c79465 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Fri, 2 Oct 2026 11:33:29 -0700 Subject: [PATCH 26/28] Skip the force field CTE tests before importing pheasy --- tests/forcefields/flows/test_cte.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/tests/forcefields/flows/test_cte.py b/tests/forcefields/flows/test_cte.py index 752dbd7117..8d4a493df3 100644 --- a/tests/forcefields/flows/test_cte.py +++ b/tests/forcefields/flows/test_cte.py @@ -5,12 +5,17 @@ are only installed in the numpy-limited forcefield CI job, so the tests are skipped in the other forcefield jobs. """ +# ruff: noqa: E402 from pathlib import Path +import pytest + +pytest.importorskip("pheasy") +gruneisen_module = pytest.importorskip("phono3py.phonon3.gruneisen") + import numpy as np import phonopy -import pytest from ase.build import bulk from jobflow import run_locally from pymatgen.io.ase import AseAtomsAdaptor @@ -19,9 +24,6 @@ from atomate2.common.schemas.cte import CTEDocument from atomate2.forcefields.flows.cte import CTEMaker -pytest.importorskip("pheasy") -gruneisen_module = pytest.importorskip("phono3py.phonon3.gruneisen") - EMT = {"@module": "ase.calculators.emt", "@callable": "EMT"} From 0a2256e76a8f9b79c2eacbdc499ebb40035176dd Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Fri, 2 Oct 2026 11:48:43 -0700 Subject: [PATCH 27/28] Move the force field pheasy test to tests/forcefields --- tests/common/jobs/test_pheasy.py | 23 -------------- tests/forcefields/flows/test_pheasy.py | 42 ++++++++++++++++++++++++++ 2 files changed, 42 insertions(+), 23 deletions(-) create mode 100644 tests/forcefields/flows/test_pheasy.py diff --git a/tests/common/jobs/test_pheasy.py b/tests/common/jobs/test_pheasy.py index 58b72dc7a0..fc1d816d29 100644 --- a/tests/common/jobs/test_pheasy.py +++ b/tests/common/jobs/test_pheasy.py @@ -29,7 +29,6 @@ _get_num_irreducible_fcs, _run_band_structure_and_plot, ) -from atomate2.forcefields.flows.pheasy import PhononMaker # fcs_cutoff_radius in Bohr. 8 Bohr (4.2 A) covers the first two neighbour # shells of fcc Cu (2.55 and 3.61 A) and stays inside the 10.8 A supercell. @@ -261,28 +260,6 @@ def test_get_num_anharmonic_supercells(monkeypatch): _get_num_anharmonic_supercells(num_disp_anhar=0, **kwargs) -def test_get_supercell_size_kwargs(monkeypatch): - received = {} - transformation = pheasy_jobs.CubicSupercellTransformation - - def record_kwargs(**kwargs): - received.update(kwargs) - return transformation(**kwargs) - - monkeypatch.setattr(pheasy_jobs, "CubicSupercellTransformation", record_kwargs) - - # the maker passes get_supercell_size_kwargs on to the job - maker = PhononMaker(get_supercell_size_kwargs={"angle_tolerance": 0.1}) - job = maker.get_supercell_matrix(_cu_structure()) - assert job.function_kwargs == {"angle_tolerance": 0.1} - - # the job passes them to CubicSupercellTransformation. The other default - # is kept. - job.function(*job.function_args, **job.function_kwargs) - assert received["angle_tolerance"] == 0.1 - assert received["allow_orthorhombic"] is False - - def test_check_lasso_alpha(tmp_dir): log_file = Path("pheasy_anharmonic_fit.log") diff --git a/tests/forcefields/flows/test_pheasy.py b/tests/forcefields/flows/test_pheasy.py new file mode 100644 index 0000000000..f0114f4e00 --- /dev/null +++ b/tests/forcefields/flows/test_pheasy.py @@ -0,0 +1,42 @@ +"""Tests for the force field pheasy workflow. + +pheasy and ALM are only installed in the numpy-limited forcefield CI job, so the +tests are skipped in the other forcefield jobs. +""" + +import pytest + +pytest.importorskip("pheasy") + +from ase.build import bulk +from pymatgen.core import Structure +from pymatgen.io.ase import AseAtomsAdaptor + +import atomate2.common.jobs.pheasy as pheasy_jobs +from atomate2.forcefields.flows.pheasy import PhononMaker + + +def _cu_structure() -> Structure: + return AseAtomsAdaptor.get_structure(bulk("Cu", "fcc", a=3.61, cubic=True)) + + +def test_get_supercell_size_kwargs(monkeypatch): + received = {} + transformation = pheasy_jobs.CubicSupercellTransformation + + def record_kwargs(**kwargs): + received.update(kwargs) + return transformation(**kwargs) + + monkeypatch.setattr(pheasy_jobs, "CubicSupercellTransformation", record_kwargs) + + # the maker passes get_supercell_size_kwargs on to the job + maker = PhononMaker(get_supercell_size_kwargs={"angle_tolerance": 0.1}) + job = maker.get_supercell_matrix(_cu_structure()) + assert job.function_kwargs == {"angle_tolerance": 0.1} + + # the job passes them to CubicSupercellTransformation. The other default + # is kept. + job.function(*job.function_args, **job.function_kwargs) + assert received["angle_tolerance"] == 0.1 + assert received["allow_orthorhombic"] is False From 3b4d22036507fe406343e23b9738e9ccb4b85687 Mon Sep 17 00:00:00 2001 From: Hrushikesh Sahasrabuddhe <111614145+hrushikesh-s@users.noreply.github.com> Date: Fri, 2 Oct 2026 11:55:22 -0700 Subject: [PATCH 28/28] Run the force field pheasy workflow end to end with EMT --- tests/forcefields/flows/test_pheasy.py | 36 ++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/tests/forcefields/flows/test_pheasy.py b/tests/forcefields/flows/test_pheasy.py index f0114f4e00..3315963a1b 100644 --- a/tests/forcefields/flows/test_pheasy.py +++ b/tests/forcefields/flows/test_pheasy.py @@ -9,11 +9,21 @@ pytest.importorskip("pheasy") from ase.build import bulk +from emmet.core.phonon import ( + PhononBS, + PhononBSDOSDoc, + PhononDOS, + ThermalDisplacementData, +) +from jobflow import run_locally from pymatgen.core import Structure from pymatgen.io.ase import AseAtomsAdaptor import atomate2.common.jobs.pheasy as pheasy_jobs from atomate2.forcefields.flows.pheasy import PhononMaker +from atomate2.forcefields.jobs import ForceFieldRelaxMaker, ForceFieldStaticMaker + +EMT = {"@module": "ase.calculators.emt", "@callable": "EMT"} def _cu_structure() -> Structure: @@ -40,3 +50,29 @@ def record_kwargs(**kwargs): job.function(*job.function_args, **job.function_kwargs) assert received["angle_tolerance"] == 0.1 assert received["allow_orthorhombic"] is False + + +def test_pheasy_wf_force_field(clean_dir): + """Run the harmonic force field pheasy workflow with EMT forces on fcc Cu.""" + # a 2x2x2 supercell of the 4-atom cubic cell + maker = PhononMaker( + min_length=7.0, + use_symmetrized_structure="conventional", + create_thermal_displacements=True, + bulk_relax_maker=ForceFieldRelaxMaker( + force_field_name=EMT, relax_kwargs={"fmax": 0.00001} + ), + static_energy_maker=ForceFieldStaticMaker(force_field_name=EMT), + phonon_displacement_maker=ForceFieldStaticMaker(force_field_name=EMT), + ) + flow = maker.make(_cu_structure()) + responses = run_locally(flow, create_folders=True, ensure_success=True) + ph_doc = responses[flow.jobs[-1].uuid][1].output + + assert isinstance(ph_doc, PhononBSDOSDoc) + assert isinstance(ph_doc.phonon_bandstructure, PhononBS) + assert isinstance(ph_doc.phonon_dos, PhononDOS) + assert isinstance(ph_doc.thermal_displacement_data, ThermalDisplacementData) + assert isinstance(ph_doc.structure, Structure) + assert ph_doc.has_imaginary_modes is False + assert isinstance(ph_doc.force_constants, list)