Merge pull request #2320 from yang-ruoxi/absorption

Add VASP input set for absorption

pymatgen/apps/borg/tests/test_queen.py

@@ -33,7 +33,7 @@ def test_get_data(self):
         drone = VaspToComputedEntryDrone()
         self.queen = BorgQueen(drone, PymatgenTest.TEST_FILES_DIR, 1)
         data = self.queen.get_data()
-        self.assertEqual(len(data), 12)
+        self.assertEqual(len(data), 14)
 
     def test_load_data(self):
         drone = VaspToComputedEntryDrone()

pymatgen/io/vasp/MPAbsorptionSet.yaml

@@ -0,0 +1,115 @@
+# Default VASP settings for frequency dependent dielectrics calculated by independent particle approximation
+# in the Materials Project.
+# Convergence against NBANDS and KPOINTS should be tested.
+# Optional modes are: independent particle approximation (IPA) or Random-phase approximation (RPA)
+
+PARENT: VASPIncarBase
+INCAR:
+  ALGO: Exact
+  EDIFF: 1.0e-8
+  ENCUT: 520
+  IBRION: -1
+  ICHARG: 1
+  ISMEAR: 0
+  SIGMA: 0.01
+  LWAVE: True
+  LREAL: False
+  NELM: 100
+  ISPIN: 2
+  PREC: Accurate
+  LOPTICS: True
+  CSHIFT: 0.1
+  NEDOS: 2001
+KPOINTS:
+  reciprocal_density: 200
+POTCAR_FUNCTIONAL: PBE
+POTCAR:
+  Ac: Ac
+  Ag: Ag
+  Al: Al
+  Ar: Ar
+  As: As
+  Au: Au
+  B: B
+  Ba: Ba_sv
+  Be: Be_sv
+  Bi: Bi
+  Br: Br
+  C: C
+  Ca: Ca_sv
+  Cd: Cd
+  Ce: Ce
+  Cl: Cl
+  Co: Co
+  Cr: Cr_pv
+  Cs: Cs_sv
+  Cu: Cu_pv
+  Dy: Dy_3
+  Er: Er_3
+  Eu: Eu
+  F: F
+  Fe: Fe_pv
+  Ga: Ga_d
+  Gd: Gd
+  Ge: Ge_d
+  H: H
+  He: He
+  Hf: Hf_pv
+  Hg: Hg
+  Ho: Ho_3
+  I: I
+  In: In_d
+  Ir: Ir
+  K: K_sv
+  Kr: Kr
+  La: La
+  Li: Li_sv
+  Lu: Lu_3
+  Mg: Mg_pv
+  Mn: Mn_pv
+  Mo: Mo_pv
+  N: N
+  Na: Na_pv
+  Nb: Nb_pv
+  Nd: Nd_3
+  Ne: Ne
+  Ni: Ni_pv
+  Np: Np
+  O: O
+  Os: Os_pv
+  P: P
+  Pa: Pa
+  Pb: Pb_d
+  Pd: Pd
+  Pm: Pm_3
+  Pr: Pr_3
+  Pt: Pt
+  Pu: Pu
+  Rb: Rb_sv
+  Re: Re_pv
+  Rh: Rh_pv
+  Ru: Ru_pv
+  S: S
+  Sb: Sb
+  Sc: Sc_sv
+  Se: Se
+  Si: Si
+  Sm: Sm_3
+  Sn: Sn_d
+  Sr: Sr_sv
+  Ta: Ta_pv
+  Tb: Tb_3
+  Tc: Tc_pv
+  Te: Te
+  Th: Th
+  Ti: Ti_pv
+  Tl: Tl_d
+  Tm: Tm_3
+  U: U
+  V: V_pv
+  W: W_pv
+  Xe: Xe
+  Y: Y_sv
+  Yb: Yb_2
+  Zn: Zn
+  Zr: Zr_sv

pymatgen/io/vasp/outputs.py

@@ -603,10 +603,10 @@ def optical_absorption_coeff(self):
             def f(freq, real, imag):
                 """
                 The optical absorption coefficient calculated in terms of
-                equation
+                equation, the unit is in cm-1
                 """
-                hbar = 6.582119514e-16  # eV/K
-                coeff = np.sqrt(np.sqrt(real**2 + imag**2) - real) * np.sqrt(2) / hbar * freq
+                hc = 1.23984 * 1e-4  # plank constant times speed of light, in the unit of eV*cm
+                coeff = np.sqrt(np.sqrt(real**2 + imag**2) - real) * np.sqrt(2) / hc * freq
                 return coeff
 
             absorption_coeff = [
@@ -621,7 +621,11 @@ def converged_electronic(self):
             True if electronic step convergence has been reached in the final
             ionic step
         """
-        final_esteps = self.ionic_steps[-1]["electronic_steps"]
+        if self.incar not in ["CHI"]:
+            final_esteps = self.ionic_steps[-1]["electronic_steps"]
+        else:
+            final_esteps = 0
+            # In a response function run there is no ionic steps, there is no scf step
         if "LEPSILON" in self.incar and self.incar["LEPSILON"]:
             i = 1
             to_check = {"e_wo_entrp", "e_fr_energy", "e_0_energy"}

pymatgen/io/vasp/sets.py

@@ -3092,3 +3092,185 @@ def get_valid_magmom_struct(structure, inplace=True, spin_mode="auto"):
     if not inplace:
         return new_struct
     return None
+
+
+class MPAbsorptionSet(MPRelaxSet):
+    """
+    MP input set for generating frequency dependent dielectrics.
+    Two modes are supported: "IPA" or "RPA".
+    A typical sequence is mode="STATIC" -> mode="IPA" -> mode="RPA"(optional)
+    For all steps other than the first one (static), the
+    recommendation is to use from_prev_calculation on the preceding run in
+    the series. It is important to ensure Gamma centred kpoints for the RPA step.
+
+    """
+
+    # CONFIG = _load_yaml_config("MPAbsorptionSet")
+
+    SUPPORTED_MODES = ("IPA", "RPA")
+
+    def __init__(
+        self,
+        structure,
+        mode="IPA",
+        copy_wavecar=True,
+        nbands=None,
+        nbands_factor=2,
+        reciprocal_density=400,
+        nkred=None,
+        nedos=2001,
+        prev_incar=None,
+        **kwargs,
+    ):
+        """
+        Args:
+            structure (Structure): Input structure.
+            prev_incar (Incar/string): Incar file from previous run.
+            mode (str): Supported modes are "IPA", "RPA"
+            nbands (int): For subsequent calculations, it is generally
+                recommended to perform NBANDS convergence starting from the
+                NBANDS of the previous run for DIAG, and to use the exact same
+                NBANDS for RPA. This parameter is used by
+                from_previous_calculation to set nband.
+            nbands_factor (int): Multiplicative factor for NBANDS when starting
+                from a previous calculation. Only applies if mode=="IPA".
+                Need to be tested for convergence.
+            reciprocal_density: the k-points density
+            nkred: the reduced number of kpoints to calculate, equal to the k-mesh. Only applies in "RPA" mode
+                  because of the q->0 limit.
+            nedos: the density of DOS, default: 2001.
+            **kwargs: All kwargs supported by DictSet. Typically, user_incar_settings is a commonly used option.
+        """
+
+        # Initialize the input set (default: IPA absorption)
+        super().__init__(structure, **kwargs)
+
+        self.prev_incar = prev_incar
+        self.nbands = nbands
+        self.reciprocal_density = reciprocal_density
+        self.mode = mode.upper()
+        if self.mode not in MPAbsorptionSet.SUPPORTED_MODES:
+            raise ValueError(f"{self.mode} not one of the support modes : {MPAbsorptionSet.SUPPORTED_MODES}")
+        self.copy_wavecar = copy_wavecar
+        self.nbands_factor = nbands_factor
+        self.nedos = nedos
+        self.nkred = nkred
+        self.kwargs = kwargs
+
+    @property
+    def kpoints(self):
+        """
+        Generate gamma center k-points mesh grid for optical calculation. It is not mandatory for 'ALGO = Exact',
+        but is requested by 'ALGO = CHI' calculation.
+        """
+        return Kpoints.automatic_density_by_vol(self.structure, self.reciprocal_density, force_gamma=True)
+
+    @property
+    def incar(self):
+        """
+        :return: Incar
+        """
+        parent_incar = super().incar
+        absorption_incar = {
+            "ALGO": "Exact",
+            "EDIFF": 1.0e-8,
+            "IBRION": -1,
+            "ICHARG": 1,
+            "ISMEAR": 0,
+            "SIGMA": 0.01,
+            "LWAVE": True,
+            "LREAL": False,  # for small cell it's more efficient to use reciprocal
+            "NELM": 100,
+            "NSW": 0,
+            "LOPTICS": True,
+            "CSHIFT": 0.1,
+            "NEDOS": 2001,
+        }
+        self._config_dict["INCAR"].update(absorption_incar)
+
+        if self.mode == "IPA":
+            # use the incar from previous static calculation
+            if self.prev_incar is not None:
+                incar = Incar(self.prev_incar)
+                # Default parameters for diagonalization calculation.
+                incar.update({"ALGO": "Exact", "LOPTICS": True, "CSHIFT": 0.1, "LWAVE": "True"})
+                incar.update({"NEDOS": self.nedos}) if self.nedos is not None else incar.update({"NEDOS": 2001})
+            else:
+                incar = Incar(parent_incar)
+
+        elif self.mode == "RPA":
+            # Default parameters for the response function calculation. NELM has to be set to 1.
+            # NOMEGA is set to 1000 in order to get smooth spectrum
+            incar = Incar(self.prev_incar) if self.prev_incar is not None else Incar(parent_incar)
+            incar.update({"ALGO": "CHI", "NELM": 1, "NOMEGA": 1000})
+
+            if self.nkred is not None:
+                incar["NKREDX"] = self.nkred[0]
+                incar["NKREDY"] = self.nkred[1]
+                incar["NKREDZ"] = self.nkred[2]
+
+            incar.pop("EDIFF", None)
+            incar.pop("LOPTICS", None)
+            incar.pop("LWAVE", False)
+
+        else:
+            raise Exception("mode has to be from 'IPA' or 'RPA'")
+
+        if self.nbands:
+            incar["NBANDS"] = self.nbands
+
+        # Respect user set INCAR.
+        incar.update(self.kwargs.get("user_incar_settings", {}))
+
+        return incar
+
+    def override_from_prev_calc(self, prev_calc_dir=".", **kwargs):
+        """
+        Update the input set to include settings from a previous calculation.
+        Args:
+            prev_calc_dir (str): The path to the previous calculation directory.
+        Returns:
+            The input set with the settings (structure, k-points, incar, etc)
+            updated using the previous VASP run.
+        """
+        vasprun, outcar = get_vasprun_outcar(prev_calc_dir)
+        self.prev_incar = vasprun.incar
+        self._structure = vasprun.final_structure
+
+        # The number of bands is multiplied by the factor
+        prev_nbands = int(vasprun.parameters["NBANDS"])
+
+        if self.mode.upper() == "IPA":
+            self.nbands = int(np.ceil(prev_nbands * self.nbands_factor))
+
+        # Since in the optical calculation, only the q->0 transition is of interests, we can reduce the number of q by
+        # the factor of the number of kpoints in each corresonding x, y, z directions. This will reduce the
+        # computational work by factor of 1/nkredx*nkredy*nkredz. An isotropic NKRED can be used for cubic
+        # lattice, but using NKREDX, NKREDY, NKREDZ is more sensible for other lattice.
+        if self.mode.upper() == "RPA":
+            # self.nbands = int(np.ceil(self.nbands * self.nbands_factor / self.ncores) * self.ncores)
+            self.nkred = vasprun.kpoints.kpts[0]
+
+        files_to_transfer = {}
+        if self.copy_wavecar:
+            for fname in ("WAVECAR", "WAVEDER"):
+                w = sorted(glob.glob(str(Path(prev_calc_dir) / (fname + "*"))))
+                if w:
+                    files_to_transfer[fname] = str(w[-1])
+
+        self.files_to_transfer.update(files_to_transfer)
+
+        return self
+
+    @classmethod
+    def from_prev_calc(cls, prev_calc_dir, mode, **kwargs):
+        """
+        Generate a set of Vasp input files for absorption calculation
+        Args:
+            prev_calc_dir (str): The directory contains the outputs(
+                vasprun.xml of previous vasp run.
+            mode (str): Supported modes are "IPA", "RPA" (default)
+            **kwargs: All kwargs supported by MPAbsorptionsSet, other than structure
+        """
+        input_set = cls(_dummy_structure, mode, **kwargs)
+        return input_set.override_from_prev_calc(prev_calc_dir=prev_calc_dir)

pymatgen/io/vasp/tests/test_outputs.py

@@ -90,7 +90,7 @@ def test_charge_charge_dielectric(self):
     def test_optical_absorption_coeff(self):
         v = Vasprun(self.TEST_FILES_DIR / "vasprun.BSE.xml.gz")
         absorption_coeff = v.optical_absorption_coeff
-        self.assertEqual(absorption_coeff[1], 24966408728.917931)
+        self.assertEqual(absorption_coeff[1], 0.13254281688694558)
 
     def test_vasprun_with_more_than_two_unlabelled_dielectric_functions(self):
         with self.assertRaises(NotImplementedError):