PhononMaker add options to calculate_pdos and save force constant…

…s to file (#1008)

* add calculate_pdos and make force constance can be save

* add calculate_pdos and make force constance can be save

* Update phonons.py

* Apply pre-commit hooks for code formatting and linting

* Fix issues found by pre-commit hooks

* add force_constants2file parameter

* - Renamed "force_constants2file" to "force_constants_filename"
- Add "create_force_constants_file" in "generate_frequencies_eigenvectors_kwargs"
- Move "force_constants_filename" and "calculate_pdos" to "generate_frequencies_eigenvectors_kwargs"

* fix doc string placement

* refactor from_forces_born handling create_force_constants_file

---------

Co-authored-by: Yaotang Zhang <[email protected]>
Co-authored-by: Janosh Riebesell <[email protected]>

src/atomate2/common/flows/phonons.py

@@ -107,6 +107,12 @@ class BasePhononMaker(Maker, ABC):
         Maker used to compute the forces for a supercell.
     generate_frequencies_eigenvectors_kwargs : dict
         Keyword arguments passed to :obj:`generate_frequencies_eigenvectors`.
+        - create_force_constants_file: bool
+            If True, a force constants file will be created
+        - force_constants_filename: str
+            If store_force_constants is True, the file name to store the force constants
+        - calculate_pdos: bool
+            If True, the projected phonon density of states will be calculated
     create_thermal_displacements: bool
         Bool that determines if thermal_displacement_matrices are computed
     kpath_scheme: str
@@ -146,7 +152,14 @@ class BasePhononMaker(Maker, ABC):
         None
     )
     create_thermal_displacements: bool = True
-    generate_frequencies_eigenvectors_kwargs: dict = field(default_factory=dict)
+    generate_frequencies_eigenvectors_kwargs: dict = field(
+        default_factory=lambda: {
+            "create_force_constants_file": False,
+            "force_constants_filename": "FORCE_CONSTANTS",
+            "calculate_pdos": False,
+        }
+    )
+
     kpath_scheme: str = "seekpath"
     code: str = None
     store_force_constants: bool = True

src/atomate2/common/schemas/phonons.py

@@ -346,11 +346,28 @@ def from_forces_born(
         phonon.produce_force_constants(forces=set_of_forces)
 
         filename_phonopy_yaml = kwargs.get("filename_phonopy_yaml", "phonopy.yaml")
+        create_force_constants_file = kwargs.get("create_force_constants_file", False)
+        force_constants_filename = kwargs.get(
+            "force_constants_filename", "FORCE_CONSTANTS"
+        )
         # if kwargs.get("filename_phonopy_yaml") is None:
         #    kwargs["filename_phonopy_yaml"] = "phonopy.yaml"
 
         # with phonopy.load("phonopy.yaml") the phonopy API can be used
-        phonon.save(filename_phonopy_yaml)
+        phonon.save(
+            filename_phonopy_yaml,
+            settings={
+                "force_constants": kwargs.get(
+                    "store_force_constants", not create_force_constants_file
+                )
+            },
+        )
+        if create_force_constants_file:
+            from phonopy.file_IO import write_FORCE_CONSTANTS
+
+            write_FORCE_CONSTANTS(  # save force_constants to text file
+                phonon.force_constants, filename=force_constants_filename
+            )
 
         # get phonon band structure
         kpath_dict, kpath_concrete = PhononBSDOSDoc.get_kpath(
@@ -361,7 +378,7 @@ def from_forces_born(
 
         npoints_band = kwargs.get("npoints_band", 101)
         qpoints, connections = get_band_qpoints_and_path_connections(
-            kpath_concrete, npoints=kwargs.get("npoints_band", 101)
+            kpath_concrete, npoints=npoints_band
         )
 
         # phonon band structures will always be computed
@@ -406,6 +423,20 @@ def from_forces_born(
             kppa=kpoint_density_dos,
             force_gamma=True,
         )
+
+        # projected dos
+        if kwargs.get("calculate_pdos", False):
+            phonon.run_mesh(
+                kpoint.kpts[0], with_eigenvectors=True, is_mesh_symmetry=False
+            )
+            phonon_dos_sigma = kwargs.get("phonon_dos_sigma")
+            dos_use_tetrahedron_method = kwargs.get("dos_use_tetrahedron_method", True)
+            phonon.run_projected_dos(
+                sigma=phonon_dos_sigma,
+                use_tetrahedron_method=dos_use_tetrahedron_method,
+            )
+            phonon.write_projected_dos()
+
         phonon.run_mesh(kpoint.kpts[0])
         phonon_dos_sigma = kwargs.get("phonon_dos_sigma")
         dos_use_tetrahedron_method = kwargs.get("dos_use_tetrahedron_method", True)