diff --git a/src/atomate2/common/flows/pheasy.py b/src/atomate2/common/flows/pheasy.py new file mode 100644 index 0000000000..94844ad1cd --- /dev/null +++ b/src/atomate2/common/flows/pheasy.py @@ -0,0 +1,406 @@ +"""Flows for calculating phonons, an-harmonic force constants, +and phonon energy renormalization with the pheasy. +""" + +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 atomate2.common.jobs.pheasy import ( + generate_frequencies_eigenvectors, + generate_phonon_displacements, + run_phonon_displacements, +) +from atomate2.common.jobs.phonons import get_supercell_size, get_total_energy_per_cell +from atomate2.common.jobs.utils import structure_to_conventional, structure_to_primitive + +if TYPE_CHECKING: + from pathlib import Path + + from emmet.core.math import Matrix3D + from pymatgen.core.structure import Structure + + from atomate2.aims.jobs.base import BaseAimsMaker + from atomate2.forcefields.jobs import ForceFieldRelaxMaker, ForceFieldStaticMaker + from atomate2.vasp.jobs.base import BaseVaspMaker + +SUPPORTED_CODES = frozenset(("vasp", "aims", "forcefields")) + + +@dataclass +class BasePhononMaker(Maker, ABC): + """ + + Maker to calculate harmonic phonons with a DFT/force field code and the LASSO-based + machine learning code Pheasy. + + Calculate the zero-K harmonic phonons of a material. 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. + With the help of pheasy (LASSO technique), these forces are then converted into a + dynamical matrix. 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. + + .. Note:: + It is heavily recommended to symmetrize the structure before passing it to + this flow. Otherwise, a different space group might be detected and too + many 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. + + Parameters + ---------- + name : str + Name of the flow produced by this maker. + sym_reduce : bool + Whether to reduce the number of deformations using symmetry. + symprec : float + Symmetry precision to use in the + reduction of symmetry to find the primitive/conventional cell + (use_primitive_standard_structure, use_conventional_standard_structure) + and to handle all symmetry-related tasks in pheasy, we recommend to + use the value of 1e-3. + displacement: float + displacement distance for phonons, for most cases 0.01 A is a good choice, + but it can be increased to 0.02 A for heavier elements. + displacement_anharmonic: float + displacement distance for anharmonic force constants(FCs) up to sixth-order + FCs, for most cases 0.08 A is a good choice, but it can be increased to 0.1 A. + num_displaced_supercells: int + 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. + num_displaced_supercells_anharmonic: 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. + 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. + get_supercell_size_kwargs: dict + kwargs that will be passed to get_supercell_size to determine supercell size + use_symmetrized_structure: str + allowed strings: "primitive", "conventional", None + + - "primitive" will enforce to start the phonon computation + from the primitive standard structure + according to Setyawan, W., & Curtarolo, S. (2010). + High-throughput electronic band structure calculations: + Challenges and tools. Computational Materials Science, + 49(2), 299-312. doi:10.1016/j.commatsci.2010.05.010. + This makes it possible to use certain k-path definitions + with this workflow. Otherwise, we must rely on seekpath + - "conventional" will enforce to start the phonon computation + from the conventional standard structure + according to Setyawan, W., & Curtarolo, S. (2010). + High-throughput electronic band structure calculations: + Challenges and tools. Computational Materials Science, + 49(2), 299-312. doi:10.1016/j.commatsci.2010.05.010. + We will however use seekpath and primitive structures + as determined by from phonopy to compute the phonon band structure + bulk_relax_maker: .ForceFieldRelaxMaker, .BaseAimsMaker, .BaseVaspMaker, or None + A maker to perform a tight relaxation on the bulk. + Set to ``None`` to skip the + bulk relaxation + static_energy_maker: .ForceFieldRelaxMaker, .BaseAimsMaker, .BaseVaspMaker, or None + A maker to perform the computation of the DFT energy on the bulk. + Set to ``None`` to skip the + static energy computation + born_maker: .ForceFieldStaticMaker, .BaseAsimsMaker, .BaseVaspMaker, or None + Maker to compute the BORN charges. + phonon_displacement_maker: .ForceFieldStaticMaker, .BaseAimsMaker, .BaseVaspMaker + Maker used to compute the forces for a supercell. + generate_frequencies_eigenvectors_kwargs : dict + Keyword arguments passed to :obj:`generate_frequencies_eigenvectors`. + create_thermal_displacements: bool + Bool that determines if thermal_displacement_matrices are computed + kpath_scheme: str + scheme to generate kpoints. Please be aware that + you can only use seekpath with any kind of cell + Otherwise, please use the standard primitive structure + Available schemes are: + "seekpath", "hinuma", "setyawan_curtarolo", "latimer_munro". + "seekpath" and "hinuma" are the same definition but + seekpath can be used with any kind of unit cell as + it relies on phonopy to handle the relationship + to the primitive cell and not pymatgen + code: str + determines the dft or force field code. + mp_id: str + The mp_id of the material in the Materials Project database. + store_force_constants: bool + if True, force constants will be stored + socket: bool + If True, use the socket for the calculation + + """ + + name: str = "phonon" + sym_reduce: bool = True + symprec: float = 1e-3 + displacement: float = 0.01 + displacement_anharmonic: float = 0.08 + num_displaced_supercells: int = 0 + num_displaced_supercells_anharmonic: int = 0 + min_length: float | None = 14.0 + prefer_90_degrees: bool = True + get_supercell_size_kwargs: dict = field(default_factory=dict) + use_symmetrized_structure: str | None = None + bulk_relax_maker: ForceFieldRelaxMaker | BaseVaspMaker | BaseAimsMaker | None = None + static_energy_maker: ForceFieldRelaxMaker | BaseVaspMaker | BaseAimsMaker | None = ( + None + ) + born_maker: ForceFieldStaticMaker | BaseVaspMaker | None = None + phonon_displacement_maker: ForceFieldStaticMaker | BaseVaspMaker | BaseAimsMaker = ( + None + ) + create_thermal_displacements: bool = False + generate_frequencies_eigenvectors_kwargs: dict = field(default_factory=dict) + kpath_scheme: str = "seekpath" + code: str = None + mp_id: str = None + store_force_constants: bool = True + socket: bool = False + + def make( + self, + structure: Structure, + prev_dir: str | Path | None = None, + born: list[Matrix3D] | None = None, + epsilon_static: Matrix3D | None = None, + total_dft_energy_per_formula_unit: float | None = None, + supercell_matrix: Matrix3D | None = None, + ) -> Flow: + """Make flow to calculate the phonon properties. + + Parameters + ---------- + structure : Structure + A pymatgen structure object. Please start with a structure + that is nearly fully optimized as the internal optimizers + have very strict settings! + prev_dir : str or Path or None + A previous calculation directory to use for copying outputs. + born: Matrix3D + Instead of recomputing born charges and epsilon, these values can also be + provided manually. If born and epsilon_static are provided, the born run + will be skipped it can be provided in the VASP convention with information + for every atom in unit cell. Please be careful when converting structures + within in this workflow as this could lead to errors + epsilon_static: Matrix3D + The high-frequency dielectric constant to use instead of recomputing born + charges and epsilon. If born, epsilon_static are provided, the born run + will be skipped + total_dft_energy_per_formula_unit: float + It has to be given per formula unit (as a result in corresponding Doc). + Instead of recomputing the energy of the bulk structure every time, this + value can also be provided in eV. If it is provided, the static run will be + skipped. This energy is the typical output dft energy of the dft workflow. + No conversion needed. + supercell_matrix: list + Instead of min_length, also a supercell_matrix can be given, e.g. + [[1.0,0.0,0.0],[0.0,1.0,0.0],[0.0,0.0,1.0] + """ + use_symmetrized_structure = self.use_symmetrized_structure + kpath_scheme = self.kpath_scheme + valid_structs = (None, "primitive", "conventional") + if use_symmetrized_structure not in valid_structs: + raise ValueError( + f"Invalid {use_symmetrized_structure=}, use one of {valid_structs}" + ) + + if use_symmetrized_structure != "primitive" and kpath_scheme != "seekpath": + raise ValueError( + f"You can't use {kpath_scheme=} with the primitive standard " + "structure, please use seekpath" + ) + + valid_schemes = ("seekpath", "hinuma", "setyawan_curtarolo", "latimer_munro") + if kpath_scheme not in valid_schemes: + raise ValueError( + f"{kpath_scheme=} is not implemented, use one of {valid_schemes}" + ) + + if self.code is None or self.code not in SUPPORTED_CODES: + raise ValueError( + "The code variable must be passed and it must be a supported code." + f" Supported codes are: {SUPPORTED_CODES}" + ) + + jobs = [] + + # TODO: should this be after or before structural optimization as the + # optimization could change the symmetry we could add a tutorial and point out + # that the structure should be nearly optimized before the phonon workflow + if self.use_symmetrized_structure == "primitive": + # These structures are compatible with many + # of the kpath algorithms that are used for Materials Project + prim_job = structure_to_primitive(structure, self.symprec) + jobs.append(prim_job) + structure = prim_job.output + elif self.use_symmetrized_structure == "conventional": + # it could be beneficial to use conventional standard structures to arrive + # faster at supercells with right angles + conv_job = structure_to_conventional(structure, self.symprec) + jobs.append(conv_job) + structure = conv_job.output + + optimization_run_job_dir = None + optimization_run_uuid = None + + if self.bulk_relax_maker is not None: + # optionally relax the structure + 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 + optimization_run_job_dir = bulk.output.dir_name + optimization_run_uuid = bulk.output.uuid + + # if supercell_matrix is None, supercell size will be determined after relax + # maker to ensure that cell lengths are really larger than threshold. + # Note that If one wants to calculate the lattice thermal conductivity, + # the supercell dimensions should be forced to be diagonal, e.g., + # supercell_matrix = [[2, 0, 0], [0, 2, 0], [0, 0, 2]] + if supercell_matrix is None: + supercell_job = get_supercell_size( + structure, + self.min_length, + self.prefer_90_degrees, + **self.get_supercell_size_kwargs, + ) + jobs.append(supercell_job) + supercell_matrix = supercell_job.output + + # Computation of static energy + total_dft_energy = None + static_run_job_dir = None + static_run_uuid = None + if (self.static_energy_maker is not None) and ( + total_dft_energy_per_formula_unit is None + ): + static_job_kwargs = {} + if self.prev_calc_dir_argname is not None: + static_job_kwargs[self.prev_calc_dir_argname] = prev_dir + static_job = self.static_energy_maker.make( + structure=structure, **static_job_kwargs + ) + jobs.append(static_job) + total_dft_energy = static_job.output.output.energy + static_run_job_dir = static_job.output.dir_name + static_run_uuid = static_job.output.uuid + prev_dir = static_job.output.dir_name + elif total_dft_energy_per_formula_unit is not None: + # to make sure that one can reuse results from Doc + compute_total_energy_job = get_total_energy_per_cell( + total_dft_energy_per_formula_unit, structure + ) + jobs.append(compute_total_energy_job) + total_dft_energy = compute_total_energy_job.output + + # get a phonon object from pheasy code using the random-displacement approach + displacements = generate_phonon_displacements( + structure=structure, + supercell_matrix=supercell_matrix, + displacement=self.displacement, + num_displaced_supercells=self.num_displaced_supercells, + sym_reduce=self.sym_reduce, + symprec=self.symprec, + use_symmetrized_structure=self.use_symmetrized_structure, + kpath_scheme=self.kpath_scheme, + code=self.code, + ) + jobs.append(displacements) + + # perform the phonon displacement calculations + displacement_calcs = run_phonon_displacements( + displacements=displacements.output, + structure=structure, + supercell_matrix=supercell_matrix, + phonon_maker=self.phonon_displacement_maker, + socket=self.socket, + prev_dir_argname=self.prev_calc_dir_argname, + prev_dir=prev_dir, + ) + jobs.append(displacement_calcs) + + # Computation of BORN charges + born_run_job_dir = None + born_run_uuid = None + if self.born_maker is not None and (born is None or epsilon_static is None): + born_kwargs = {} + if self.prev_calc_dir_argname is not None: + born_kwargs[self.prev_calc_dir_argname] = prev_dir + born_job = self.born_maker.make(structure, **born_kwargs) + jobs.append(born_job) + + # I am not happy how we currently access "born" charges + # This is very vasp specific code aims and forcefields + # do not support this at the moment, if this changes we have + # to update this section + epsilon_static = born_job.output.calcs_reversed[0].output.epsilon_static + born = born_job.output.calcs_reversed[0].output.outcar["born"] + born_run_job_dir = born_job.output.dir_name + born_run_uuid = born_job.output.uuid + + phonon_collect = generate_frequencies_eigenvectors( + supercell_matrix=supercell_matrix, + displacement=self.displacement, + sym_reduce=self.sym_reduce, + symprec=self.symprec, + use_symmetrized_structure=self.use_symmetrized_structure, + kpath_scheme=self.kpath_scheme, + code=self.code, + mp_id=self.mp_id, + structure=structure, + displacement_data=displacement_calcs.output, + epsilon_static=epsilon_static, + born=born, + total_dft_energy=total_dft_energy, + static_run_job_dir=static_run_job_dir, + static_run_uuid=static_run_uuid, + born_run_job_dir=born_run_job_dir, + born_run_uuid=born_run_uuid, + optimization_run_job_dir=optimization_run_job_dir, + optimization_run_uuid=optimization_run_uuid, + create_thermal_displacements=self.create_thermal_displacements, + store_force_constants=self.store_force_constants, + **self.generate_frequencies_eigenvectors_kwargs, + ) + + jobs.append(phonon_collect) + + # create a flow including all jobs for a phonon computation + return Flow(jobs, phonon_collect.output) + + @property + @abstractmethod + def prev_calc_dir_argname(self) -> str | None: + """Name of argument informing static maker of previous calculation directory. + + As this differs between different DFT codes (e.g., VASP, CP2K), it + has been left as a property to be implemented by the inheriting class. + + Note: this is only applicable if a relax_maker is specified; i.e., two + calculations are performed for each ordering (relax -> static) + """ diff --git a/src/atomate2/common/jobs/pheasy.py b/src/atomate2/common/jobs/pheasy.py new file mode 100644 index 0000000000..5c47d21105 --- /dev/null +++ b/src/atomate2/common/jobs/pheasy.py @@ -0,0 +1,371 @@ +"""Jobs for running phonon calculations.""" + +from __future__ import annotations + +import contextlib +import logging +import warnings +from typing import TYPE_CHECKING + +import numpy as np +from jobflow import Flow, Response, job +from phonopy import Phonopy +from pymatgen.core import Structure +from pymatgen.io.phonopy import get_phonopy_structure, get_pmg_structure +from pymatgen.phonon.bandstructure import PhononBandStructureSymmLine +from pymatgen.phonon.dos import PhononDos + +from atomate2.common.schemas.pheasy import Forceconstants, PhononBSDOSDoc +from atomate2.common.schemas.phonons import get_factor + +if TYPE_CHECKING: + from pathlib import Path + + from emmet.core.math import Matrix3D + + from atomate2.aims.jobs.base import BaseAimsMaker + from atomate2.forcefields.jobs import ForceFieldStaticMaker + from atomate2.vasp.jobs.base import BaseVaspMaker + + +logger = logging.getLogger(__name__) + + +@job(data=[Structure]) +def generate_phonon_displacements( + structure: Structure, + supercell_matrix: np.array, + num_displaced_supercells: int, + displacement: float, + sym_reduce: bool, + symprec: float, + use_symmetrized_structure: str | None, + kpath_scheme: str, + code: str, +) -> list[Structure]: + """ + + Generate small-distance perturbed structures with phonopy based on two ways: + (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. + + Parameters + ---------- + structure: Structure object + Fully optimized input structure for phonon run + supercell_matrix: np.array + array to describe supercell matrix + displacement: float + displacement in Angstrom (default: 0.01) + num_displaced_supercells: int + number of displaced supercells defined by users + sym_reduce: bool + if True, symmetry will be used to generate displacements + symprec: float + precision to determine symmetry + use_symmetrized_structure: str or None + primitive, conventional or None + kpath_scheme: str + scheme to generate kpath + code: + code to perform the computations + + """ + warnings.warn( + "Initial magnetic moments will not be considered for the determination " + "of the symmetry of the structure and thus will be removed now.", + stacklevel=1, + ) + cell = get_phonopy_structure( + structure.remove_site_property(property_name="magmom") + if "magmom" in structure.site_properties + else structure + ) + factor = get_factor(code) + + # a bit of code repetition here as I currently + # do not see how to pass the phonopy object? + if use_symmetrized_structure == "primitive" and kpath_scheme != "seekpath": + primitive_matrix: np.ndarray | str = np.eye(3) + else: + primitive_matrix = "auto" + + # TARP: THIS IS BAD! Including for discussions sake + if cell.magnetic_moments is not None and primitive_matrix == "auto": + if np.any(cell.magnetic_moments != 0.0): + raise ValueError( + "For materials with magnetic moments, " + "use_symmetrized_structure must be 'primitive'" + ) + cell.magnetic_moments = None + # create the phonopy object to get some information + # for the displacement generation in ALM code. + phonon = Phonopy( + cell, + supercell_matrix, + primitive_matrix=primitive_matrix, + factor=factor, + symprec=symprec, + is_symmetry=sym_reduce, + ) + + # 1. the ALM module is used to determine how many free parameters + # (irreducible force constants) of second order force constants (FCs) + # within the 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 matrix can not always guarantee the accurate result sometimes, you + # may need to displace more random configurations. At least use one or + # two more configurations based on the suggested number of displacements. + + try: + from alm import ALM + except ImportError as e: + logging.exception( + f"Error importing ALM: {e}. 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) + + # 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. + phonon.generate_displacements(distance=displacement) + num_disp_f = len(phonon.displacements) + if num_disp_f > 3: + num_d = int(np.ceil(num * 1.8)) + else: + pass + + logger.info( + "The number of free parameters of Second Order Force Constants is %s", n_fp + ) + logger.info("") + + logger.info( + "The Number of Equations Used to Obtain the 2ND FCs is %s", 3 * natom * num + ) + logger.info("") + + logger.warning( + "Be Careful!!! Full Rank of Matrix cannot always guarantee the correct result\ + sometimes.\n" + "If the total atoms in the supercell are less than 100 and\n" + "lattice constants are less than 10 angstroms,\n" + "I highly suggest displacing more random configurations.\n" + "At least use one or two more configurations based on the suggested\ + number of displacements." + ) + logger.info("") + + if num_disp_f > 3: + if num_displaced_supercells != 0: + phonon.generate_displacements( + distance=displacement, + number_of_snapshots=num_displaced_supercells, + random_seed=103, + ) + else: + phonon.generate_displacements( + distance=displacement, + number_of_snapshots=num_d, + random_seed=103, + ) + else: + pass + + supercells = phonon.supercells_with_displacements + displacements = [get_pmg_structure(cell) for cell in supercells] + + # add the equilibrium structure to the list for calculating + # the residual forces. + displacements.append(get_pmg_structure(phonon.supercell)) + return displacements + + +@job( + output_schema=PhononBSDOSDoc, + data=[PhononDos, PhononBandStructureSymmLine, Forceconstants], +) +def generate_frequencies_eigenvectors( + structure: Structure, + supercell_matrix: np.array, + displacement: float, + sym_reduce: bool, + symprec: float, + use_symmetrized_structure: str | None, + kpath_scheme: str, + code: str, + mp_id: str, + displacement_data: dict[str, list], + total_dft_energy: float, + epsilon_static: Matrix3D = None, + born: Matrix3D = None, + **kwargs, +) -> PhononBSDOSDoc: + """ + Analyze the phonon runs and summarize the results. + + Parameters + ---------- + structure: Structure object + Fully optimized structure used for phonon runs + supercell_matrix: np.array + array to describe supercell + displacement: float + displacement in Angstrom used for supercell computation + sym_reduce: bool + if True, symmetry will be used in phonopy + symprec: float + precision to determine symmetry + use_symmetrized_structure: str + primitive, conventional, None are allowed + kpath_scheme: str + kpath scheme for phonon band structure computation + code: str + code to run computations + displacement_data: dict + outputs from displacements + total_dft_energy: float + total DFT energy in eV per cell + epsilon_static: Matrix3D + The high-frequency dielectric constant + born: Matrix3D + Born charges + kwargs: dict + Additional parameters that are passed to PhononBSDOSDoc.from_forces_born + """ + return PhononBSDOSDoc.from_forces_born( + structure=structure.remove_site_property(property_name="magmom") + if "magmom" in structure.site_properties + else structure, + supercell_matrix=supercell_matrix, + displacement=displacement, + sym_reduce=sym_reduce, + symprec=symprec, + use_symmetrized_structure=use_symmetrized_structure, + kpath_scheme=kpath_scheme, + code=code, + mp_id=mp_id, + displacement_data=displacement_data, + total_dft_energy=total_dft_energy, + epsilon_static=epsilon_static, + born=born, + **kwargs, + ) + + +# I did not directly import this job from the phonon module +# because I modified the job to pass the displaced structures +# to the output. +@job(data=["forces", "displaced_structures"]) +def run_phonon_displacements( + displacements: list[Structure], + structure: Structure, + supercell_matrix: Matrix3D, + phonon_maker: BaseVaspMaker | ForceFieldStaticMaker | BaseAimsMaker = None, + prev_dir: str | Path = None, + prev_dir_argname: str = None, + socket: bool = False, +) -> Flow: + """ + Run phonon displacements. + + Note, this job will replace itself with N displacement calculations, + or a single socket calculation for all displacements. + + Parameters + ---------- + displacements: Sequence + All displacements to calculate + structure: Structure object + Fully optimized structure used for phonon computations. + supercell_matrix: Matrix3D + supercell matrix for meta data + phonon_maker : .BaseVaspMaker or .ForceFieldStaticMaker or .BaseAimsMaker + A maker to use to generate dispacement calculations + prev_dir: str or Path + The previous working directory + prev_dir_argname: str + argument name for the prev_dir variable + socket: bool + If True use the socket-io interface to increase performance + """ + phonon_jobs = [] + outputs: dict[str, list] = { + "displacement_number": [], + "forces": [], + "uuids": [], + "dirs": [], + "displaced_structures": [], + } + phonon_job_kwargs = {} + if prev_dir is not None and prev_dir_argname is not None: + phonon_job_kwargs[prev_dir_argname] = prev_dir + + if socket: + phonon_job = phonon_maker.make(displacements, **phonon_job_kwargs) + info = { + "original_structure": structure, + "supercell_matrix": supercell_matrix, + "displaced_structures": displacements, + } + phonon_job.update_maker_kwargs( + {"_set": {"write_additional_data->phonon_info:json": info}}, dict_mod=True + ) + phonon_jobs.append(phonon_job) + outputs["displacement_number"] = list(range(len(displacements))) + outputs["uuids"] = [phonon_job.output.uuid] * len(displacements) + outputs["dirs"] = [phonon_job.output.dir_name] * len(displacements) + outputs["forces"] = phonon_job.output.output.all_forces + # add the displaced structures, still need to be careful with the order, + # experimental feature + outputs["displaced_structures"] = displacements + else: + for idx, displacement in enumerate(displacements): + if prev_dir is not None: + phonon_job = phonon_maker.make(displacement, prev_dir=prev_dir) + else: + phonon_job = phonon_maker.make(displacement) + phonon_job.append_name(f" {idx + 1}/{len(displacements)}") + + # we will add some meta data + info = { + "displacement_number": idx, + "original_structure": structure, + "supercell_matrix": supercell_matrix, + "displaced_structure": displacement, + } + with contextlib.suppress(Exception): + phonon_job.update_maker_kwargs( + {"_set": {"write_additional_data->phonon_info:json": info}}, + dict_mod=True, + ) + phonon_jobs.append(phonon_job) + outputs["displacement_number"].append(idx) + outputs["uuids"].append(phonon_job.output.uuid) + outputs["dirs"].append(phonon_job.output.dir_name) + outputs["forces"].append(phonon_job.output.output.forces) + outputs["displaced_structures"].append(displacement) + + displacement_flow = Flow(phonon_jobs, outputs) + return Response(replace=displacement_flow) diff --git a/src/atomate2/common/schemas/pheasy.py b/src/atomate2/common/schemas/pheasy.py new file mode 100644 index 0000000000..ecc6f396f6 --- /dev/null +++ b/src/atomate2/common/schemas/pheasy.py @@ -0,0 +1,757 @@ +"""Schemas for phonon documents.""" + +import copy +import logging +import pickle +import subprocess +from pathlib import Path +from typing import Optional, Union + +import numpy as np + +# import lib by jiongzhi zheng +from ase.io import read +from emmet.core.math import Matrix3D +from emmet.core.structure import StructureMetadata +from hiphive import ( + ClusterSpace, + ForceConstantPotential, + ForceConstants, + enforce_rotational_sum_rules +) + +from hiphive.cutoffs import estimate_maximum_cutoff +from hiphive.utilities import extract_parameters +from monty.json import MSONable +from phonopy import Phonopy +from phonopy.file_IO import parse_FORCE_CONSTANTS +from phonopy.interface.vasp import write_vasp +from phonopy.phonon.band_structure import get_band_qpoints_and_path_connections +from phonopy.structure.symmetry import symmetrize_borns_and_epsilon +from pydantic import Field +from pymatgen.core import Structure +from pymatgen.io.phonopy import ( + get_ph_bs_symm_line, + get_ph_dos, + get_phonopy_structure, + get_pmg_structure, +) +from pymatgen.io.vasp import Kpoints +from pymatgen.phonon.bandstructure import PhononBandStructureSymmLine +from pymatgen.phonon.dos import PhononDos +from pymatgen.phonon.plotter import PhononBSPlotter, PhononDosPlotter +from pymatgen.symmetry.bandstructure import HighSymmKpath +from pymatgen.symmetry.kpath import KPathSeek +from typing_extensions import Self + +# import some classmethod directly from phonons +from atomate2.common.schemas.phonons import ( + PhononComputationalSettings, + PhononJobDirs, + PhononUUIDs, + ThermalDisplacementData, + get_factor, +) + +logger = logging.getLogger(__name__) + + +class Forceconstants(MSONable): + """A force constants class.""" + + def __init__(self, force_constants: list[list[Matrix3D]]) -> None: + self.force_constants = force_constants + + +class PhononBSDOSDoc(StructureMetadata, extra="allow"): # type: ignore[call-arg] + """Collection of all data produced by the phonon workflow.""" + + structure: Optional[Structure] = Field( + None, description="Structure of Materials Project." + ) + + phonon_bandstructure: Optional[PhononBandStructureSymmLine] = Field( + None, + description="Phonon band structure object.", + ) + + phonon_dos: Optional[PhononDos] = Field( + None, + description="Phonon density of states object.", + ) + + free_energies: Optional[list[float]] = Field( + None, + description="vibrational part of the free energies in J/mol per " + "formula unit for temperatures in temperature_list", + ) + + heat_capacities: Optional[list[float]] = Field( + None, + description="heat capacities in J/K/mol per " + "formula unit for temperatures in temperature_list", + ) + + internal_energies: Optional[list[float]] = Field( + None, + description="internal energies in J/mol per " + "formula unit for temperatures in temperature_list", + ) + entropies: Optional[list[float]] = Field( + None, + description="entropies in J/(K*mol) per formula unit" + "for temperatures in temperature_list ", + ) + + temperatures: Optional[list[int]] = Field( + None, + description="temperatures at which the vibrational" + " part of the free energies" + " and other properties have been computed", + ) + + total_dft_energy: Optional[float] = Field("total DFT energy per formula unit in eV") + + has_imaginary_modes: Optional[bool] = Field( + None, description="if true, structure has imaginary modes" + ) + + # needed, e.g. to compute Grueneisen parameter etc + force_constants: Optional[Forceconstants] = Field( + None, description="Force constants between every pair of atoms in the structure" + ) + + born: Optional[list[Matrix3D]] = Field( + None, + description="born charges as computed from phonopy. Only for symmetrically " + "different atoms", + ) + + epsilon_static: Optional[Matrix3D] = Field( + None, description="The high-frequency dielectric constant" + ) + + supercell_matrix: Matrix3D = Field("matrix describing the supercell") + primitive_matrix: Matrix3D = Field( + "matrix describing relationship to primitive cell" + ) + + code: str = Field("String describing the code for the computation") + + phonopy_settings: PhononComputationalSettings = Field( + "Field including settings for Phonopy" + ) + + thermal_displacement_data: Optional[ThermalDisplacementData] = Field( + "Includes all data of the computation of the thermal displacements" + ) + + jobdirs: Optional[PhononJobDirs] = Field( + "Field including all relevant job directories" + ) + + uuids: Optional[PhononUUIDs] = Field("Field including all relevant uuids") + + @classmethod + def from_forces_born( + cls, + structure: Structure, + supercell_matrix: np.array, + displacement: float, + sym_reduce: bool, + symprec: float, + use_symmetrized_structure: Union[str, None], + kpath_scheme: str, + code: str, + mp_id: str, + displacement_data: dict[str, list], + total_dft_energy: float, + epsilon_static: Matrix3D = None, + born: Matrix3D = None, + **kwargs, + ) -> Self: + """Generate collection of phonon data. + + Parameters + ---------- + structure: Structure object + supercell_matrix: numpy array describing the supercell + displacement: float + size of displacement in angstrom + sym_reduce: bool + if True, phonopy will use symmetry + symprec: float + precision to determine kpaths, + primitive cells and symmetry in phonopy and pymatgen + use_symmetrized_structure: str + primitive, conventional or None + kpath_scheme: str + kpath scheme to generate phonon band structure + code: str + which code was used for computation + displacement_data: + output of the displacement data + total_dft_energy: float + total energy in eV per cell + epsilon_static: Matrix3D + The high-frequency dielectric constant + born: Matrix3D + born charges + **kwargs: + additional arguments + """ + factor = get_factor(code) + # This opens the opportunity to add support for other codes + # that are supported by phonopy + + cell = get_phonopy_structure(structure) + + if use_symmetrized_structure == "primitive": + primitive_matrix: Union[np.ndarray, str] = np.eye(3) + else: + primitive_matrix = "auto" + + # TARP: THIS IS BAD! Including for discussions sake + if cell.magnetic_moments is not None and primitive_matrix == "auto": + if np.any(cell.magnetic_moments != 0.0): + raise ValueError( + "For materials with magnetic moments, " + "use_symmetrized_structure must be 'primitive'" + ) + cell.magnetic_moments = None + + # Create the phonon object using the phonopy API to write the POSCAR and + # SPOSCAR files for the input of pheasy code. + phonon = Phonopy( + cell, + supercell_matrix, + primitive_matrix=primitive_matrix, + factor=factor, + symprec=symprec, + is_symmetry=sym_reduce, + ) + + # Write the POSCAR and SPOSCAR files for the input of pheasy code + supercell = phonon.get_supercell() + write_vasp("POSCAR", cell) + write_vasp("SPOSCAR", supercell) + + # get the force-displacement dataset from previous calculations + dataset_forces = [np.array(forces) for forces in displacement_data["forces"]] + dataset_forces_array = np.array(dataset_forces) + + # To deduct the residual forces on an equilibrium structure to eliminate the + # fitting error + dataset_forces_array_rr = dataset_forces_array - dataset_forces_array[-1, :, :] + + # force matrix on the displaced structures + dataset_forces_array_disp = dataset_forces_array_rr[:-1, :, :] + dataset_disps = [ + np.array(disps.cart_coords) + for disps in displacement_data["displaced_structures"] + ] + + # get the displacement dataset + dataset_disps_array_rr = np.round( + (dataset_disps - supercell.get_positions()), decimals=16 + ).astype('double') + dataset_disps_array_use = dataset_disps_array_rr[:-1, :, :] + + # get the number of displacements for harmonic phonon calculation + num_har = dataset_disps_array_use.shape[0] + + # save the displacement and force matrix in the current directory + # for the future use by pheasy code + with open("disp_matrix.pkl", "wb") as file: + pickle.dump(dataset_disps_array_use, file) + with open("force_matrix.pkl", "wb") as file: + pickle.dump(dataset_forces_array_disp, file) + + # get the born charges and dielectric constant + if born is not None and epsilon_static is not None: + if len(structure) == len(born): + borns, epsilon = symmetrize_borns_and_epsilon( + ucell=phonon.unitcell, + borns=np.array(born), + epsilon=np.array(epsilon_static), + symprec=symprec, + primitive_matrix=phonon.primitive_matrix, + supercell_matrix=phonon.supercell_matrix, + is_symmetry=kwargs.get("symmetrize_born", True), + ) + else: + raise ValueError( + "Number of born charges does not agree with number of atoms" + ) + if code == "vasp" and not np.all(np.isclose(borns, 0.0)): + phonon.nac_params = { + "born": borns, + "dielectric": epsilon, + "factor": 14.399652, + } + # Other codes could be added here + else: + borns = None + epsilon = None + + prim = read("POSCAR") + supercell = 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. + pheasy_cmd_1 = ( + f'pheasy --dim "{int(supercell_matrix[0][0])}" "{int(supercell_matrix[1][1])}" ' + f'"{int(supercell_matrix[2][2])}" -s -w 2 --symprec "{float(symprec)}" --nbody 2' + ) + + # Create the null space to further reduce the free parameters for + # specific force constants and make them physically correct. + pheasy_cmd_2 = ( + f'pheasy --dim "{int(supercell_matrix[0][0])}" "{int(supercell_matrix[1][1])}" ' + f'"{int(supercell_matrix[2][2])}" -c --symprec "{float(symprec)}" -w 2' + ) + + # Generate the Compressive Sensing matrix,i.e., displacement matrix + # for the input of machine leaning method.i.e., LASSO, + pheasy_cmd_3 = ( + f'pheasy --dim "{int(supercell_matrix[0][0])}" "{int(supercell_matrix[1][1])}" ' + f'"{int(supercell_matrix[2][2])}" -w 2 -d --symprec "{float(symprec)}" ' + f'--ndata "{int(num_har)}" --disp_file' + ) + + # Here we set a criteria to determine which method to use to generate the + # force constants. If the number of displacements is larger than 3, we + # will use the LASSO method to generate the force constants. Otherwise, + # we will use the least-squred method to generate the force constants. + phonon.generate_displacements(distance=displacement) + disps = phonon.displacements + num_judge = len(disps) + + if num_judge > 3: + # Calculate the force constants using the LASSO method due to the + # random-displacement method Obviously, the rotaional invariance + # constraint, i.e., tag: --rasr BHH, is enforced during the + # fitting process. + pheasy_cmd_4 = ( + f'pheasy --dim "{int(supercell_matrix[0][0])}" "{int(supercell_matrix[1][1])}" ' + f'"{int(supercell_matrix[2][2])}" -f --full_ifc -w 2 --symprec "{float(symprec)}" ' + f'-l LASSO --std --rasr BHH --ndata "{int(num_har)}"' + ) + + else: + # Calculate the force constants using the least-squred method + pheasy_cmd_4 = ( + f'pheasy --dim "{int(supercell_matrix[0][0])}" "{int(supercell_matrix[1][1])}" ' + f'"{int(supercell_matrix[2][2])}" -f --full_ifc -w 2 --symprec "{float(symprec)}" ' + f'--rasr BHH --ndata "{int(num_har)}"' + ) + + logger.info("Start running pheasy in cluster") + + subprocess.call(pheasy_cmd_1, shell=True) + subprocess.call(pheasy_cmd_2, shell=True) + subprocess.call(pheasy_cmd_3, shell=True) + subprocess.call(pheasy_cmd_4, shell=True) + + # Read the force constants from the output file of pheasy code + force_constants = parse_FORCE_CONSTANTS(filename="FORCE_CONSTANTS") + phonon.force_constants = force_constants + # symmetrize the force constants to make them physically correct based on + # the space group symmetry of the crystal structure. + phonon.symmetrize_force_constants() + + # with phonopy.load("phonopy.yaml") the phonopy API can be used + phonon.save("phonopy.yaml") + + # get phonon band structure + kpath_dict, kpath_concrete = PhononBSDOSDoc.get_kpath( + structure=get_pmg_structure(phonon.primitive), + kpath_scheme=kpath_scheme, + symprec=symprec, + ) + + npoints_band = kwargs.get("npoints_band", 101) + qpoints, connections = get_band_qpoints_and_path_connections( + kpath_concrete, npoints=kwargs.get("npoints_band", 101) + ) + + # phonon band structures will always be computed + filename_band_yaml = "phonon_band_structure.yaml" + + # TODO: potentially add kwargs to avoid computation of eigenvectors + phonon.run_band_structure( + qpoints, + path_connections=connections, + with_eigenvectors=kwargs.get("band_structure_eigenvectors", False), + is_band_connection=kwargs.get("band_structure_eigenvectors", False), + ) + phonon.write_yaml_band_structure(filename=filename_band_yaml) + bs_symm_line = get_ph_bs_symm_line( + filename_band_yaml, labels_dict=kpath_dict, has_nac=born is not None + ) + new_plotter = PhononBSPlotter(bs=bs_symm_line) + new_plotter.save_plot( + filename=kwargs.get("filename_bs", "phonon_band_structure.pdf"), + units=kwargs.get("units", "THz"), + ) + + # will determine if imaginary modes are present in the structure + imaginary_modes = bs_symm_line.has_imaginary_freq( + tol=kwargs.get("tol_imaginary_modes", 1e-5) + ) + + # 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. + + if imaginary_modes: + # Define a cluster space using the largest cutoff you can + max_cutoff = estimate_maximum_cutoff(supercell) - 0.01 + cutoffs = [max_cutoff] # only second order needed + cs = ClusterSpace(prim, cutoffs) + + # import the phonopy force constants using the correct supercell also + # provided by phonopy + fcs = ForceConstants.read_phonopy(supercell, "FORCE_CONSTANTS") + + # Find the parameters that best fits the force constants given you + # cluster space + parameters = extract_parameters(fcs, cs) + + # Enforce the rotational sum rules + parameters_rot = enforce_rotational_sum_rules( + cs, parameters, ["Huang","Born-Huang"], alpha=1e-6 + ) + + # use the new parameters to make a fcp and then create the force + # constants and write to a phonopy file + fcp = ForceConstantPotential(cs, parameters_rot) + fcs = fcp.get_force_constants(supercell) + fcs.write_to_phonopy("FORCE_CONSTANTS_new", format="text") + + force_constants = parse_FORCE_CONSTANTS(filename="FORCE_CONSTANTS_new") + phonon.force_constants = force_constants + phonon.symmetrize_force_constants() + + phonon.run_band_structure( + qpoints, path_connections=connections, with_eigenvectors=True + ) + phonon.write_yaml_band_structure(filename=filename_band_yaml) + bs_symm_line = get_ph_bs_symm_line( + filename_band_yaml, labels_dict=kpath_dict, has_nac=born is not None + ) + + new_plotter = PhononBSPlotter(bs=bs_symm_line) + + new_plotter.save_plot( + filename=kwargs.get("filename_bs", "phonon_band_structure.pdf"), + units=kwargs.get("units", "THz"), + ) + + # new_plotter.save_plot("phonon_band_structure.eps", + # img_format=kwargs.get("img_format", "eps"), + # units=kwargs.get("units", "THz"),) + + imaginary_modes_hiphive = bs_symm_line.has_imaginary_freq( + tol=kwargs.get("tol_imaginary_modes", 1e-5) + ) + + else: + imaginary_modes_hiphive = False + imaginary_modes = False + + # Using a shorter cutoff (10 A) to generate the force constants to + # eliminate the imaginary modes near Gamma point in phesay code + if imaginary_modes_hiphive: + pheasy_cmd_11 = ( + f'pheasy --dim "{int(supercell_matrix[0][0])}" ' + f'"{int(supercell_matrix[1][1])}" ' + f'"{int(supercell_matrix[2][2])}" ' + f'-s -w 2 --c2 10.0 ' + f'--symprec "{float(symprec)}" ' + f'--nbody 2' + ) + + pheasy_cmd_12 = ( + f'pheasy --dim "{int(supercell_matrix[0][0])}" ' + f'"{int(supercell_matrix[1][1])}" ' + f'"{int(supercell_matrix[2][2])}" ' + f'-c --symprec "{float(symprec)}" ' + f'--c2 10.0 -w 2' + ) + + pheasy_cmd_13 = ( + f'pheasy --dim "{int(supercell_matrix[0][0])}" ' + f'"{int(supercell_matrix[1][1])}" ' + f'"{int(supercell_matrix[2][2])}" ' + f'-w 2 -d --symprec "{float(symprec)}" ' + f'--c2 10.0 --ndata "{int(num_har)}" ' + f'--disp_file' + ) + + phonon.generate_displacements(distance=displacement) + disps = phonon.displacements + num_judge = len(disps) + + if num_judge > 3: + pheasy_cmd_14 = ( + f'pheasy --dim "{int(supercell_matrix[0][0])}" ' + f'"{int(supercell_matrix[1][1])}" ' + f'"{int(supercell_matrix[2][2])}" ' + f'-f --c2 10.0 --full_ifc -w 2 ' + f'--symprec "{float(symprec)}" ' + f'-l LASSO --std --rasr BHH ' + f'--ndata "{int(num_har)}"' + ) + + else: + pheasy_cmd_14 = ( + f'pheasy --dim "{int(supercell_matrix[0][0])}" ' + f'"{int(supercell_matrix[1][1])}" ' + f'"{int(supercell_matrix[2][2])}" ' + f'-f --full_ifc --c2 10.0 -w 2 ' + f'--symprec "{float(symprec)}" ' + f'--rasr BHH --ndata "{int(num_har)}"' + ) + + logger.info("Start running pheasy in cluster") + + subprocess.call(pheasy_cmd_11, shell=True) + subprocess.call(pheasy_cmd_12, shell=True) + subprocess.call(pheasy_cmd_13, shell=True) + subprocess.call(pheasy_cmd_14, shell=True) + + force_constants = parse_FORCE_CONSTANTS(filename="FORCE_CONSTANTS") + phonon.force_constants = force_constants + phonon.symmetrize_force_constants() + + # with phonon.load("phonopy.yaml") the phonopy API can be used + phonon.save("phonopy.yaml") + + # get phonon band structure + kpath_dict, kpath_concrete = cls.get_kpath( + structure=get_pmg_structure(phonon.primitive), + kpath_scheme=kpath_scheme, + symprec=symprec, + ) + + npoints_band = kwargs.get("npoints_band", 101) + qpoints, connections = get_band_qpoints_and_path_connections( + kpath_concrete, npoints=kwargs.get("npoints_band", 101) + ) + + # phonon band structures will always be cmouted + filename_band_yaml = "phonon_band_structure.yaml" + phonon.run_band_structure( + qpoints, path_connections=connections, with_eigenvectors=True + ) + phonon.write_yaml_band_structure(filename=filename_band_yaml) + bs_symm_line = get_ph_bs_symm_line( + filename_band_yaml, labels_dict=kpath_dict, has_nac=born is not None + ) + new_plotter = PhononBSPlotter(bs=bs_symm_line) + + new_plotter.save_plot( + filename=kwargs.get("filename_bs", "phonon_band_structure.pdf"), + units=kwargs.get("units", "THz") + ) + + imaginary_modes_cutoff = bs_symm_line.has_imaginary_freq( + tol=kwargs.get("tol_imaginary_modes", 1e-5) + ) + imaginary_modes = imaginary_modes_cutoff + # new_plotter.save_plot( + # "phonon_band_structure.eps", + # img_format=kwargs.get("img_format", "eps"), + # units=kwargs.get("units", "THz"), + # ) + else: + pass + + # gets data for visualization on website - yaml is also enough + if kwargs.get("band_structure_eigenvectors"): + bs_symm_line.write_phononwebsite("phonon_website.json") + + # get phonon density of states + filename_dos_yaml = "phonon_dos.yaml" + + kpoint_density_dos = kwargs.get("kpoint_density_dos", 7_000) + kpoint = Kpoints.automatic_density( + structure=get_pmg_structure(phonon.primitive), + kppa=kpoint_density_dos, + force_gamma=True, + ) + phonon.run_mesh(kpoint.kpts[0]) + phonon.run_total_dos() + phonon.write_total_dos(filename=filename_dos_yaml) + dos = get_ph_dos(filename_dos_yaml) + new_plotter_dos = PhononDosPlotter() + new_plotter_dos.add_dos(label="total", dos=dos) + new_plotter_dos.save_plot( + filename=kwargs.get("filename_dos", "phonon_dos.pdf"), + units=kwargs.get("units", "THz"), + ) + + # compute vibrational part of free energies per formula unit + temperature_range = np.arange( + kwargs.get("tmin", 0), kwargs.get("tmax", 500), kwargs.get("tstep", 10) + ) + + free_energies = [ + dos.helmholtz_free_energy( + temp=temp, structure=get_pmg_structure(phonon.primitive) + ) + for temp in temperature_range + ] + + entropies = [ + dos.entropy(temp=temp, structure=get_pmg_structure(phonon.primitive)) + for temp in temperature_range + ] + + internal_energies = [ + dos.internal_energy( + temp=temp, structure=get_pmg_structure(phonon.primitive) + ) + for temp in temperature_range + ] + + heat_capacities = [ + dos.cv(temp=temp, structure=get_pmg_structure(phonon.primitive)) + for temp in temperature_range + ] + + # will compute thermal displacement matrices + # for the primitive cell (phonon.primitive!) + # only this is available in phonopy + if kwargs.get("create_thermal_displacements"): + phonon.run_mesh( + kpoint.kpts[0], with_eigenvectors=True, is_mesh_symmetry=False + ) + freq_min_thermal_displacements = kwargs.get( + "freq_min_thermal_displacements", 0.0 + ) + phonon.run_thermal_displacement_matrices( + t_min=kwargs.get("tmin_thermal_displacements", 0), + t_max=kwargs.get("tmax_thermal_displacements", 500), + t_step=kwargs.get("tstep_thermal_displacements", 100), + freq_min=freq_min_thermal_displacements, + ) + + temperature_range_thermal_displacements = np.arange( + kwargs.get("tmin_thermal_displacements", 0), + kwargs.get("tmax_thermal_displacements", 500), + kwargs.get("tstep_thermal_displacements", 100), + ) + for idx, temp in enumerate(temperature_range_thermal_displacements): + phonon.thermal_displacement_matrices.write_cif( + phonon.primitive, idx, filename=f"tdispmat_{temp}K.cif" + ) + _disp_mat = phonon._thermal_displacement_matrices # noqa: SLF001 + tdisp_mat = _disp_mat.thermal_displacement_matrices.tolist() + + tdisp_mat_cif = _disp_mat.thermal_displacement_matrices_cif.tolist() + + else: + tdisp_mat = None + tdisp_mat_cif = None + + formula_units = ( + structure.composition.num_atoms + / structure.composition.reduced_composition.num_atoms + ) + + total_dft_energy_per_formula_unit = ( + total_dft_energy / formula_units if total_dft_energy is not None else None + ) + + return cls.from_structure( + structure=structure, + meta_structure=structure, + phonon_bandstructure=bs_symm_line, + phonon_dos=dos, + free_energies=free_energies, + internal_energies=internal_energies, + heat_capacities=heat_capacities, + entropies=entropies, + temperatures=temperature_range.tolist(), + total_dft_energy=total_dft_energy_per_formula_unit, + has_imaginary_modes=imaginary_modes, + force_constants={"force_constants": phonon.force_constants.tolist()} + if kwargs["store_force_constants"] + else None, + born=borns.tolist() if borns is not None else None, + epsilon_static=epsilon.tolist() if epsilon is not None else None, + supercell_matrix=phonon.supercell_matrix.tolist(), + primitive_matrix=phonon.primitive_matrix.tolist(), + code=code, + mp_id=mp_id, + thermal_displacement_data={ + "temperatures_thermal_displacements": temperature_range_thermal_displacements.tolist(), # noqa: E501 + "thermal_displacement_matrix_cif": tdisp_mat_cif, + "thermal_displacement_matrix": tdisp_mat, + "freq_min_thermal_displacements": freq_min_thermal_displacements, + } + if kwargs.get("create_thermal_displacements") + else None, + jobdirs={ + "displacements_job_dirs": displacement_data["dirs"], + "static_run_job_dir": kwargs["static_run_job_dir"], + "born_run_job_dir": kwargs["born_run_job_dir"], + "optimization_run_job_dir": kwargs["optimization_run_job_dir"], + "taskdoc_run_job_dir": str(Path.cwd()), + }, + uuids={ + "displacements_uuids": displacement_data["uuids"], + "born_run_uuid": kwargs["born_run_uuid"], + "optimization_run_uuid": kwargs["optimization_run_uuid"], + "static_run_uuid": kwargs["static_run_uuid"], + }, + phonopy_settings={ + "npoints_band": npoints_band, + "kpath_scheme": kpath_scheme, + "kpoint_density_dos": kpoint_density_dos, + }, + ) + + @staticmethod + def get_kpath( + structure: Structure, kpath_scheme: str, symprec: float, **kpath_kwargs + ) -> tuple: + """Get high-symmetry points in k-space in phonopy format. + + Parameters + ---------- + structure: Structure Object + kpath_scheme: str + string describing kpath + symprec: float + precision for symmetry determination + **kpath_kwargs: + additional parameters that can be passed to this method as a dict + """ + if kpath_scheme in ("setyawan_curtarolo", "latimer_munro", "hinuma"): + high_symm_kpath = HighSymmKpath( + structure, path_type=kpath_scheme, symprec=symprec, **kpath_kwargs + ) + kpath = high_symm_kpath.kpath + elif kpath_scheme == "seekpath": + high_symm_kpath = KPathSeek(structure, symprec=symprec, **kpath_kwargs) + kpath = high_symm_kpath._kpath # noqa: SLF001 + else: + raise ValueError(f"Unexpected {kpath_scheme=}") + + path = copy.deepcopy(kpath["path"]) + + for set_idx, label_set in enumerate(kpath["path"]): + for lbl_idx, label in enumerate(label_set): + path[set_idx][lbl_idx] = kpath["kpoints"][label] + return kpath["kpoints"], path + \ No newline at end of file diff --git a/src/atomate2/forcefields/flows/pheasy.py b/src/atomate2/forcefields/flows/pheasy.py new file mode 100644 index 0000000000..d9db491809 --- /dev/null +++ b/src/atomate2/forcefields/flows/pheasy.py @@ -0,0 +1,145 @@ +"""Flows for calculating phonons.""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import Literal + +from atomate2 import SETTINGS +from atomate2.common.flows.pheasy import BasePhononMaker +from atomate2.forcefields.jobs import ForceFieldRelaxMaker, ForceFieldStaticMaker + + +@dataclass +class PhononMaker(BasePhononMaker): + """ + Maker to calculate harmonic phonons with a force field. + + Calculate the harmonic phonons of a material. Initially, a tight structural + relaxation is performed to obtain a structure without forces on the atoms. + Subsequently, supercells with one displaced atom are generated and accurate + forces are computed for these structures. With the help of phonopy, these + forces are then converted into a dynamical matrix. To correct for polarization + effects, a correction of the dynamical matrix based on BORN charges can + be performed. The BORN charges can be supplied manually. + Finally, phonon densities of states, phonon band structures + and thermodynamic properties are computed. + + .. Note:: + It is heavily recommended to symmetrize the structure before passing it to + this flow. Otherwise, a different space group might be detected and too + many displacement calculations will be generated. + 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. + + Parameters + ---------- + name : str + Name of the flows produced by this maker. + sym_reduce : bool + Whether to reduce the number of deformations using symmetry. + symprec : float + Symmetry precision to use in the + reduction of symmetry to find the primitive/conventional cell + (use_primitive_standard_structure, use_conventional_standard_structure) + and to handle all symmetry-related tasks in phonopy + displacement: float + displacement distance for phonons + min_length: float + min length of the supercell that will be built + prefer_90_degrees: bool + if set to True, supercell algorithm will first try to find a supercell + with 3 90 degree angles + get_supercell_size_kwargs: dict + kwargs that will be passed to get_supercell_size to determine supercell size + use_symmetrized_structure: str + allowed strings: "primitive", "conventional", None + + - "primitive" will enforce to start the phonon computation + from the primitive standard structure + according to Setyawan, W., & Curtarolo, S. (2010). + High-throughput electronic band structure calculations: + Challenges and tools. Computational Materials Science, + 49(2), 299-312. doi:10.1016/j.commatsci.2010.05.010. + This makes it possible to use certain k-path definitions + with this workflow. Otherwise, we must rely on seekpath + - "conventional" will enforce to start the phonon computation + from the conventional standard structure + according to Setyawan, W., & Curtarolo, S. (2010). + High-throughput electronic band structure calculations: + Challenges and tools. Computational Materials Science, + 49(2), 299-312. doi:10.1016/j.commatsci.2010.05.010. + We will however use seekpath and primitive structures + as determined by from phonopy to compute the phonon band structure + bulk_relax_maker : .ForceFieldRelaxMaker or None + A maker to perform a tight relaxation on the bulk. + Set to ``None`` to skip the + bulk relaxation + static_energy_maker : .ForceFieldRelaxMaker or None + A maker to perform the computation of the DFT energy on the bulk. + Set to ``None`` to skip the + static energy computation + born_maker: .ForceFieldStaticMaker or None + Maker to compute the BORN charges. + phonon_displacement_maker : .ForceFieldStaticMaker or None + Maker used to compute the forces for a supercell. + generate_frequencies_eigenvectors_kwargs : dict + Keyword arguments passed to :obj:`generate_frequencies_eigenvectors`. + create_thermal_displacements: bool + Bool that determines if thermal_displacement_matrices are computed + kpath_scheme: str + scheme to generate kpoints. Please be aware that + you can only use seekpath with any kind of cell + Otherwise, please use the standard primitive structure + Available schemes are: + "seekpath", "hinuma", "setyawan_curtarolo", "latimer_munro". + "seekpath" and "hinuma" are the same definition but + seekpath can be used with any kind of unit cell as + it relies on phonopy to handle the relationship + to the primitive cell and not pymatgen + code: str + determines the dft or force field code. + store_force_constants: bool + if True, force constants will be stored + socket: bool + If True, use the socket for the calculation + """ + + name: str = "phonon" + sym_reduce: bool = True + symprec: float = SETTINGS.PHONON_SYMPREC + displacement: float = 0.01 + min_length: float | None = 20.0 + prefer_90_degrees: bool = True + get_supercell_size_kwargs: dict = field(default_factory=dict) + use_symmetrized_structure: Literal["primitive", "conventional"] | None = None + bulk_relax_maker: ForceFieldRelaxMaker | None = field( + default_factory=lambda: ForceFieldRelaxMaker( + force_field_name="MACE", relax_kwargs={"fmax": 0.00001} + ) + ) + static_energy_maker: ForceFieldStaticMaker | None = field( + default_factory=lambda: ForceFieldStaticMaker(force_field_name="MACE") + ) + phonon_displacement_maker: ForceFieldStaticMaker = field( + default_factory=lambda: ForceFieldStaticMaker(force_field_name="MACE") + ) + create_thermal_displacements: bool = False + generate_frequencies_eigenvectors_kwargs: dict = field(default_factory=dict) + kpath_scheme: str = "seekpath" + store_force_constants: bool = True + code: str = "forcefields" + born_maker: ForceFieldStaticMaker | None = None + + @property + def prev_calc_dir_argname(self) -> None: + """Name of argument informing static maker of previous calculation directory. + + As this differs between different DFT codes (e.g., VASP, CP2K), it + has been left as a property to be implemented by the inheriting class. + + Note: this is only applicable if a relax_maker is specified; i.e., two + calculations are performed for each ordering (relax -> static) + """ + return diff --git a/src/atomate2/vasp/flows/pheasy.py b/src/atomate2/vasp/flows/pheasy.py new file mode 100644 index 0000000000..df777d50ea --- /dev/null +++ b/src/atomate2/vasp/flows/pheasy.py @@ -0,0 +1,155 @@ +"""Define the VASP PhononMaker.""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import TYPE_CHECKING + +from atomate2 import SETTINGS +from atomate2.common.flows.pheasy import BasePhononMaker +from atomate2.vasp.flows.core import DoubleRelaxMaker +from atomate2.vasp.jobs.core import DielectricMaker, StaticMaker, TightRelaxMaker +from atomate2.vasp.jobs.phonons import PhononDisplacementMaker +from atomate2.vasp.sets.core import StaticSetGenerator + +if TYPE_CHECKING: + from atomate2.vasp.jobs.base import BaseVaspMaker + + +@dataclass +class PhononMaker(BasePhononMaker): + """ + Maker to calculate harmonic phonons with VASP and Phonopy. + + Calculate the harmonic phonons of a material. Initially, a tight structural + relaxation is performed to obtain a structure without forces on the atoms. + Subsequently, supercells with one displaced atom are generated and accurate + forces are computed for these structures. With the help of phonopy, these + forces are then converted into a dynamical matrix. 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. + + .. Note:: + It is heavily recommended to symmetrize the structure before passing it to + this flow. Otherwise, a different space group might be detected and too + many displacement calculations will be generated. + 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. + + Parameters + ---------- + name : str = "phonon" + Name of the flows produced by this maker. + sym_reduce : bool = True + Whether to reduce the number of deformations using symmetry. + symprec : float = 1e-4 + Symmetry precision to use in the + reduction of symmetry to find the primitive/conventional cell + (use_primitive_standard_structure, use_conventional_standard_structure) + and to handle all symmetry-related tasks in phonopy + displacement: float = 0.01 + displacement distance for phonons + min_length: float = 20.0 + min length of the supercell that will be built + prefer_90_degrees: bool = True + if set to True, supercell algorithm will first try to find a supercell + with 3 90 degree angles + get_supercell_size_kwargs: dict = {} + kwargs that will be passed to get_supercell_size to determine supercell size + use_symmetrized_structure: str or None = None + allowed strings: "primitive", "conventional", None + + - "primitive" will enforce to start the phonon computation + from the primitive standard structure + according to Setyawan, W., & Curtarolo, S. (2010). + High-throughput electronic band structure calculations: + Challenges and tools. Computational Materials Science, + 49(2), 299-312. doi:10.1016/j.commatsci.2010.05.010. + This makes it possible to use certain k-path definitions + with this workflow. Otherwise, we must rely on seekpath + - "conventional" will enforce to start the phonon computation + from the conventional standard structure + according to Setyawan, W., & Curtarolo, S. (2010). + High-throughput electronic band structure calculations: + Challenges and tools. Computational Materials Science, + 49(2), 299-312. doi:10.1016/j.commatsci.2010.05.010. + We will however use seekpath and primitive structures + as determined by from phonopy to compute the phonon band structure + bulk_relax_maker : .BaseVaspMaker or None + A maker to perform a tight relaxation on the bulk. + Set to ``None`` to skip the + bulk relaxation + static_energy_maker : .BaseVaspMaker or None + A maker to perform the computation of the DFT energy on the bulk. + Set to ``None`` to skip the + static energy computation + born_maker: .BaseVaspMaker or None + Maker to compute the BORN charges. + phonon_displacement_maker : .BaseVaspMaker or None + Maker used to compute the forces for a supercell. + generate_frequencies_eigenvectors_kwargs : dict + Keyword arguments passed to :obj:`generate_frequencies_eigenvectors`. + create_thermal_displacements: bool + Bool that determines if thermal_displacement_matrices are computed + kpath_scheme: str = "seekpath" + scheme to generate kpoints. Please be aware that + you can only use seekpath with any kind of cell + Otherwise, please use the standard primitive structure + Available schemes are: + "seekpath", "hinuma", "setyawan_curtarolo", "latimer_munro". + "seekpath" and "hinuma" are the same definition but + seekpath can be used with any kind of unit cell as + it relies on phonopy to handle the relationship + to the primitive cell and not pymatgen + code: str = "vasp" + determines the DFT code. currently only vasp is implemented. + This keyword might enable the implementation of other codes + in the future + store_force_constants: bool + if True, force constants will be stored + socket: bool + If True, use the socket for the calculation + """ + + name: str = "phonon" + sym_reduce: bool = True + symprec: float = SETTINGS.PHONON_SYMPREC + displacement: float = 0.01 + num_displaced_supercells: int = 0 + min_length: float | None = 14.0 + prefer_90_degrees: bool = True + get_supercell_size_kwargs: dict = field(default_factory=dict) + use_symmetrized_structure: str | None = None + create_thermal_displacements: bool = True + generate_frequencies_eigenvectors_kwargs: dict = field(default_factory=dict) + kpath_scheme: str = "seekpath" + store_force_constants: bool = True + socket: bool = False + code: str = "vasp" + bulk_relax_maker: BaseVaspMaker | None = field( + default_factory=lambda: DoubleRelaxMaker.from_relax_maker(TightRelaxMaker()) + ) + static_energy_maker: BaseVaspMaker | None = field( + default_factory=lambda: StaticMaker( + input_set_generator=StaticSetGenerator(auto_ispin=True) + ) + ) + born_maker: BaseVaspMaker | None = field(default_factory=DielectricMaker) + phonon_displacement_maker: BaseVaspMaker = field( + default_factory=PhononDisplacementMaker + ) + + @property + def prev_calc_dir_argname(self) -> str: + """Name of argument informing static maker of previous calculation directory. + + As this differs between different DFT codes (e.g., VASP, CP2K), it + has been left as a property to be implemented by the inheriting class. + + Note: this is only applicable if a relax_maker is specified; i.e., two + calculations are performed for each ordering (relax -> static) + """ + return "prev_dir" + \ No newline at end of file