Merge pull request #276 from materialsproject/phonon-tweaks

Phonon tweaks

src/atomate2/settings.py

@@ -48,11 +48,6 @@ class Atomate2Settings(BaseSettings):
     VASP_INCAR_UPDATES: dict = Field(
         default_factory=dict, description="Updates to apply to VASP INCAR files."
     )
-    VASP_RELAX_MAX_FORCE: float = Field(
-        0.25,
-        description="Maximum force allowed on each atom for successful structure "
-        "optimization",
-    )
     VASP_VOLUME_CHANGE_WARNING_TOL: float = Field(
         0.2,
         description="Maximum volume change allowed in VASP relaxations before the "

src/atomate2/vasp/flows/phonons.py

@@ -18,6 +18,7 @@
     get_total_energy_per_cell,
     run_phonon_displacements,
 )
+from atomate2.vasp.sets.core import StaticSetGenerator
 
 if TYPE_CHECKING:
     from pathlib import Path
@@ -137,7 +138,11 @@ class PhononMaker(Maker):
     bulk_relax_maker: BaseVaspMaker | None = field(
         default_factory=lambda: DoubleRelaxMaker.from_relax_maker(TightRelaxMaker())
     )
-    static_energy_maker: BaseVaspMaker | None = field(default_factory=StaticMaker)
+    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
@@ -163,34 +168,29 @@ def make(
         Parameters
         ----------
         structure : .Structure
-            A pymatgen structure. Please start with a structure
-            that is nearly fully optimized as the internal optimizers
-            have very strict settings!
+            A pymatgen structure. Please start with a structure that is nearly fully
+            optimized as the internal optimizers have very strict settings!
         prev_vasp_dir : str or Path or None
             A previous vasp 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, 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
+            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
-            Instead of recomputing born charges and epsilon,
-            these values can also be provided.
-            if born, epsilon_static are provided, the born
-            run will be skipped
+            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.
+            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]
+            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]
         """
         if self.use_symmetrized_structure not in [None, "primitive", "conventional"]:
             raise ValueError(
@@ -216,39 +216,35 @@ def make(
 
         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
+        # 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
+            # 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 = self.bulk_relax_maker.make(structure, prev_vasp_dir=prev_vasp_dir)
             jobs.append(bulk)
             structure = bulk.output.structure
+            prev_vasp_dir = bulk.output.dir_name
             optimization_run_job_dir = bulk.output.dir_name
             optimization_run_uuid = bulk.output.uuid
-        else:
-            optimization_run_job_dir = None
-            optimization_run_uuid = None
 
-        # if supercell_matrix is None, supercell size will be determined
-        # after relax maker to ensure that cell lengths are really larger
-        # than threshold
+        # if supercell_matrix is None, supercell size will be determined after relax
+        # maker to ensure that cell lengths are really larger than threshold
         if supercell_matrix is None:
             supercell_job = get_supercell_size(
                 structure,
@@ -259,6 +255,29 @@ def make(
             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 = self.static_energy_maker.make(
+                structure=structure, prev_vasp_dir=prev_vasp_dir
+            )
+            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_vasp_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 phonopy
         displacements = generate_phonon_displacements(
             structure=structure,
@@ -278,34 +297,15 @@ def make(
             structure=structure,
             supercell_matrix=supercell_matrix,
             phonon_maker=self.phonon_displacement_maker,
+            prev_vasp_dir=prev_vasp_dir,
         )
         jobs.append(vasp_displacement_calcs)
 
-        # Computation of static energy
-        if (self.static_energy_maker is not None) and (
-            total_dft_energy_per_formula_unit is None
-        ):
-            static_job = self.static_energy_maker.make(structure=structure)
-            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
-        else:
-            if 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
-            else:
-                total_dft_energy = None
-            static_run_job_dir = None
-            static_run_uuid = None
-
         # 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_job = self.born_maker.make(structure)
+            born_job = self.born_maker.make(structure, prev_vasp_dir=prev_vasp_dir)
             jobs.append(born_job)
 
             # I am not happy how we currently access "born" charges
@@ -314,9 +314,6 @@ def make(
             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
-        else:
-            born_run_job_dir = None
-            born_run_uuid = None
 
         phonon_collect = generate_frequencies_eigenvectors(
             supercell_matrix=supercell_matrix,
@@ -341,7 +338,7 @@ def make(
             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)

src/atomate2/vasp/jobs/phonons.py

@@ -21,6 +21,8 @@
 from atomate2.vasp.sets.core import StaticSetGenerator
 
 if TYPE_CHECKING:
+    from pathlib import Path
+
     import numpy as np
     from emmet.core.math import Matrix3D
     from pymatgen.core import Structure
@@ -255,6 +257,7 @@ def run_phonon_displacements(
     structure: Structure,
     supercell_matrix,
     phonon_maker: BaseVaspMaker = None,
+    prev_vasp_dir: str | Path = None,
 ):
     """
     Run phonon displacements.
@@ -270,6 +273,8 @@ def run_phonon_displacements(
         supercell matrix for meta data
     phonon_maker : .BaseVaspMaker
         A VaspMaker to use to generate the elastic relaxation jobs.
+    prev_vasp_dir : str or Path or None
+        A previous vasp calculation directory to use for copying outputs.
     """
     if phonon_maker is None:
         phonon_maker = PhononDisplacementMaker()
@@ -283,7 +288,7 @@ def run_phonon_displacements(
     }
 
     for i, displacement in enumerate(displacements):
-        phonon_job = phonon_maker.make(displacement)
+        phonon_job = phonon_maker.make(displacement, prev_vasp_dir=prev_vasp_dir)
         phonon_job.append_name(f" {i + 1}/{len(displacements)}")
 
         # we will add some meta data
@@ -341,10 +346,9 @@ class PhononDisplacementMaker(BaseVaspMaker):
     """
 
     name: str = "phonon static"
-
     input_set_generator: VaspInputGenerator = field(
         default_factory=lambda: StaticSetGenerator(
-            user_kpoints_settings={"grid_density": 7000},
+            user_kpoints_settings={"reciprocal_density": 100},
             user_incar_settings={
                 "IBRION": 2,
                 "ISIF": 3,
@@ -355,7 +359,7 @@ class PhononDisplacementMaker(BaseVaspMaker):
                 "ALGO": "Normal",
                 "NSW": 0,
                 "LCHARG": False,
-                "ISMEAR": 0,
             },
+            auto_ispin=True,
         )
     )

tests/test_data/vasp/NaCl_phonons/phonon_static_1_2/inputs/INCAR

@@ -5,7 +5,7 @@ ENCUT = 700
 GGA = Ps
 IBRION = 2
 ISIF = 3
-ISMEAR = 0
+ISMEAR = -5
 ISPIN = 2
 LAECHG = False
 LASPH = True

tests/test_data/vasp/NaCl_phonons/phonon_static_2_2/inputs/INCAR

@@ -5,7 +5,7 @@ ENCUT = 700
 GGA = Ps
 IBRION = 2
 ISIF = 3
-ISMEAR = 0
+ISMEAR = -5
 ISPIN = 2
 LAECHG = False
 LASPH = True

tests/test_data/vasp/Si_phonons_1/phonon_static_1_1/inputs/INCAR

@@ -5,7 +5,7 @@ ENCUT = 700
 GGA = Ps
 IBRION = 2
 ISIF = 3
-ISMEAR = 0
+ISMEAR = -5
 ISPIN = 2
 LAECHG = False
 LASPH = True

tests/test_data/vasp/Si_phonons_2/phonon_static_1_1/inputs/INCAR

@@ -5,7 +5,7 @@ ENCUT = 700
 GGA = Ps
 IBRION = 2
 ISIF = 3
-ISMEAR = 0
+ISMEAR = -5
 ISPIN = 2
 LAECHG = False
 LASPH = True

tests/test_data/vasp/Si_phonons_3/phonon_static_1_1/inputs/INCAR

@@ -5,7 +5,7 @@ ENCUT = 700
 GGA = Ps
 IBRION = 2
 ISIF = 3
-ISMEAR = 0
+ISMEAR = -5
 ISPIN = 2
 LAECHG = False
 LASPH = True

tests/test_data/vasp/Si_phonons_4/phonon_static_1_1/inputs/INCAR

@@ -5,7 +5,7 @@ ENCUT = 700
 GGA = Ps
 IBRION = 2
 ISIF = 3
-ISMEAR = 0
+ISMEAR = -5
 ISPIN = 2
 LAECHG = False
 LASPH = True