From 1c51748569e5df5369213883f724d1db6e537d99 Mon Sep 17 00:00:00 2001 From: Janosh Riebesell Date: Wed, 5 Jun 2024 18:45:22 -0400 Subject: [PATCH] split pymatgen/io/vasp/sets.py into sub-modules --- pymatgen/io/vasp/sets/__init__.py | 40 + pymatgen/io/vasp/{sets.py => sets/base.py} | 1656 +------------------- pymatgen/io/vasp/sets/lobster.py | 105 ++ pymatgen/io/vasp/sets/matpes.py | 84 + pymatgen/io/vasp/sets/mit.py | 211 +++ pymatgen/io/vasp/sets/mp.py | 853 ++++++++++ pymatgen/io/vasp/sets/mvl.py | 443 ++++++ 7 files changed, 1750 insertions(+), 1642 deletions(-) create mode 100644 pymatgen/io/vasp/sets/__init__.py rename pymatgen/io/vasp/{sets.py => sets/base.py} (50%) create mode 100644 pymatgen/io/vasp/sets/lobster.py create mode 100644 pymatgen/io/vasp/sets/matpes.py create mode 100644 pymatgen/io/vasp/sets/mit.py create mode 100644 pymatgen/io/vasp/sets/mp.py create mode 100644 pymatgen/io/vasp/sets/mvl.py diff --git a/pymatgen/io/vasp/sets/__init__.py b/pymatgen/io/vasp/sets/__init__.py new file mode 100644 index 00000000000..eb3621d706a --- /dev/null +++ b/pymatgen/io/vasp/sets/__init__.py @@ -0,0 +1,40 @@ +from __future__ import annotations + +from pymatgen.io.vasp.sets.base import ( + MODULE_DIR, + BadInputSetWarning, + DictSet, + UserPotcarFunctional, + VaspInputGenerator, + VaspInputSet, + _load_yaml_config, + batch_write_input, + get_structure_from_prev_run, + get_valid_magmom_struct, +) +from pymatgen.io.vasp.sets.lobster import LobsterSet +from pymatgen.io.vasp.sets.matpes import MatPESStaticSet +from pymatgen.io.vasp.sets.mit import MITMDSet, MITNEBSet, MITRelaxSet +from pymatgen.io.vasp.sets.mp import ( + MPAbsorptionSet, + MPHSEBSSet, + MPHSERelaxSet, + MPMDSet, + MPMetalRelaxSet, + MPNMRSet, + MPNonSCFSet, + MPRelaxSet, + MPScanRelaxSet, + MPScanStaticSet, + MPSOCSet, + MPStaticSet, +) +from pymatgen.io.vasp.sets.mvl import ( + MVLElasticSet, + MVLGBSet, + MVLGWSet, + MVLNPTMDSet, + MVLRelax52Set, + MVLScanRelaxSet, + MVLSlabSet, +) diff --git a/pymatgen/io/vasp/sets.py b/pymatgen/io/vasp/sets/base.py similarity index 50% rename from pymatgen/io/vasp/sets.py rename to pymatgen/io/vasp/sets/base.py index b2e383401fe..52f74b0b05a 100644 --- a/pymatgen/io/vasp/sets.py +++ b/pymatgen/io/vasp/sets/base.py @@ -36,9 +36,8 @@ from copy import deepcopy from dataclasses import dataclass, field from glob import glob -from itertools import chain from pathlib import Path -from typing import TYPE_CHECKING, cast +from typing import TYPE_CHECKING, Literal, Union, cast import numpy as np from monty.dev import deprecated @@ -46,31 +45,24 @@ from monty.serialization import loadfn from pymatgen.analysis.structure_matcher import StructureMatcher -from pymatgen.core import Element, PeriodicSite, SiteCollection, Species, Structure +from pymatgen.core import SiteCollection, Structure from pymatgen.io.core import InputGenerator from pymatgen.io.vasp.inputs import Incar, Kpoints, Poscar, Potcar, VaspInput from pymatgen.io.vasp.outputs import Outcar, Vasprun from pymatgen.symmetry.analyzer import SpacegroupAnalyzer from pymatgen.symmetry.bandstructure import HighSymmKpath -from pymatgen.util.due import Doi, due from pymatgen.util.typing import Kpoint if TYPE_CHECKING: - from typing import Any, Literal, Union + from typing import Any from typing_extensions import Self - from pymatgen.util.typing import Vector3D - UserPotcarFunctional = Union[ - Literal["PBE", "PBE_52", "PBE_54", "LDA", "LDA_52", "LDA_54", "PW91", "LDA_US", "PW91_US"], None - ] - -MODULE_DIR = Path(__file__).resolve().parent -# TODO (janosh): replace with following line once PMG is py3.9+ only -# which guarantees absolute path to __file__ -# reason: faster and shorter than pathlib (slow import and less OOP overhead) -# MODULE_DIR = os.path.dirname(__file__) +UserPotcarFunctional = Union[ + Literal["PBE", "PBE_52", "PBE_54", "LDA", "LDA_52", "LDA_54", "PW91", "LDA_US", "PW91_US"], None +] +MODULE_DIR = Path(__file__).resolve().parent.parent def _load_yaml_config(fname): @@ -1196,1525 +1188,6 @@ def primes_less_than(max_val: int) -> list[int]: return res -@due.dcite( - Doi("10.1016/j.commatsci.2011.02.023"), - description="A high-throughput infrastructure for density functional theory calculations", -) -@dataclass -class MITRelaxSet(VaspInputSet): - """ - Standard implementation of VaspInputSet utilizing parameters in the MIT - High-throughput project. - The parameters are chosen specifically for a high-throughput project, - which means in general pseudopotentials with fewer electrons were chosen. - - Args: - structure (Structure): The Structure to create inputs for. If None, the input - set is initialized without a Structure but one must be set separately before - the inputs are generated. - **kwargs: Keywords supported by VaspInputSet. - - Please refer:: - - A Jain, G. Hautier, C. Moore, S. P. Ong, C. Fischer, T. Mueller, - K. A. Persson, G. Ceder. A high-throughput infrastructure for density - functional theory calculations. Computational Materials Science, - 2011, 50(8), 2295-2310. doi:10.1016/j.commatsci.2011.02.023 - """ - - CONFIG = _load_yaml_config("MITRelaxSet") - - -@dataclass -class MPRelaxSet(VaspInputSet): - """ - Implementation of VaspInputSet utilizing parameters in the public - Materials Project. Typically, the pseudopotentials chosen contain more - electrons than the MIT parameters, and the k-point grid is ~50% more dense. - The LDAUU parameters are also different due to the different PSPs used, - which result in different fitted values. - - Args: - structure (Structure): The Structure to create inputs for. If None, the input - set is initialized without a Structure but one must be set separately before - the inputs are generated. - **kwargs: Keywords supported by VaspInputSet. - """ - - CONFIG = _load_yaml_config("MPRelaxSet") - - -@due.dcite( - Doi("10.1021/acs.jpclett.0c02405"), - description="AccurAccurate and Numerically Efficient r2SCAN Meta-Generalized Gradient Approximation", -) -@due.dcite( - Doi("10.1103/PhysRevLett.115.036402"), - description="Strongly Constrained and Appropriately Normed Semilocal Density Functional", -) -@due.dcite( - Doi("10.1103/PhysRevB.93.155109"), - description="Efficient generation of generalized Monkhorst-Pack grids through the use of informatics", -) -@dataclass -class MPScanRelaxSet(VaspInputSet): - """Write a relaxation input set using the accurate and numerically - efficient r2SCAN variant of the Strongly Constrained and Appropriately Normed - (SCAN) metaGGA density functional. - - Notes: - 1. This functional is officially supported in VASP 6.0.0 and above. On older version, - source code may be obtained by contacting the authors of the referenced manuscript. - The original SCAN functional, available from VASP 5.4.3 onwards, maybe used instead - by passing `user_incar_settings={"METAGGA": "SCAN"}` when instantiating this InputSet. - r2SCAN and SCAN are expected to yield very similar results. - - 2. Meta-GGA calculations require POTCAR files that include - information on the kinetic energy density of the core-electrons, - i.e. "PBE_52" or "PBE_54". Make sure the POTCARs include the - following lines (see VASP wiki for more details): - - $ grep kinetic POTCAR - kinetic energy-density - mkinetic energy-density pseudized - kinetic energy density (partial) - - Args: - bandgap (float): Bandgap of the structure in eV. The bandgap is used to - compute the appropriate k-point density and determine the - smearing settings. - Metallic systems (default, bandgap = 0) use a KSPACING value of 0.22 - and Methfessel-Paxton order 2 smearing (ISMEAR=2, SIGMA=0.2). - Non-metallic systems (bandgap > 0) use the tetrahedron smearing - method (ISMEAR=-5, SIGMA=0.05). The KSPACING value is - calculated from the bandgap via Eqs. 25 and 29 of Wisesa, McGill, - and Mueller [1] (see References). Note that if 'user_incar_settings' - or 'user_kpoints_settings' override KSPACING, the calculation from - bandgap is not performed. - vdw (str): set "rVV10" to enable SCAN+rVV10, which is a versatile - van der Waals density functional by combing the SCAN functional - with the rVV10 non-local correlation functional. rvv10 is the only - dispersion correction available for SCAN at this time. - **kwargs: Keywords supported by VaspInputSet. - - References: - [1] P. Wisesa, K.A. McGill, T. Mueller, Efficient generation of - generalized Monkhorst-Pack grids through the use of informatics, - Phys. Rev. B. 93 (2016) 1-10. doi:10.1103/PhysRevB.93.155109. - - References: - James W. Furness, Aaron D. Kaplan, Jinliang Ning, John P. Perdew, and Jianwei Sun. - Accurate and Numerically Efficient r2SCAN Meta-Generalized Gradient Approximation. - The Journal of Physical Chemistry Letters 0, 11 DOI: 10.1021/acs.jpclett.0c02405 - """ - - bandgap: float | None = None - auto_kspacing: bool = True - user_potcar_functional: UserPotcarFunctional = "PBE_54" - auto_ismear: bool = True - CONFIG = _load_yaml_config("MPSCANRelaxSet") - _valid_potcars: Sequence[str] | None = ("PBE_52", "PBE_54") - - def __post_init__(self): - super().__post_init__() - if self.vdw and self.vdw != "rvv10": - warnings.warn("Use of van der waals functionals other than rVV10 with SCAN is not supported at this time. ") - # delete any vdw parameters that may have been added to the INCAR - vdw_par = loadfn(str(MODULE_DIR / "vdW_parameters.yaml")) - for k in vdw_par[self.vdw]: - self._config_dict["INCAR"].pop(k, None) - - -@dataclass -class MPMetalRelaxSet(VaspInputSet): - """ - Implementation of VaspInputSet utilizing parameters in the public - Materials Project, but with tuning for metals. Key things are a denser - k point density, and a. - """ - - CONFIG = MPRelaxSet.CONFIG - - @property - def incar_updates(self) -> dict: - """Updates to the INCAR config for this calculation type.""" - return {"ISMEAR": 1, "SIGMA": 0.2} - - @property - def kpoints_updates(self) -> dict | Kpoints: - """Updates to the kpoints configuration for this calculation type.""" - return {"reciprocal_density": 200} - - -@dataclass -class MPHSERelaxSet(VaspInputSet): - """Same as the MPRelaxSet, but with HSE parameters.""" - - CONFIG = _load_yaml_config("MPHSERelaxSet") - - -@dataclass -class MPStaticSet(VaspInputSet): - """Create input files for a static calculation. - - Args: - structure (Structure): Structure from previous run. - lepsilon (bool): Whether to add static dielectric calculation - lcalcpol (bool): Whether to turn on evaluation of the Berry phase approximations - for electronic polarization - reciprocal_density (int): For static calculations, we usually set the - reciprocal density by volume. This is a convenience arg to change - that, rather than using user_kpoints_settings. Defaults to 100, - which is ~50% more than that of standard relaxation calculations. - small_gap_multiply ([float, float]): If the gap is less than - 1st index, multiply the default reciprocal_density by the 2nd - index. - **kwargs: Keywords supported by MPRelaxSet. - """ - - lepsilon: bool = False - lcalcpol: bool = False - reciprocal_density: int = 100 - small_gap_multiply: tuple[float, float] | None = None - inherit_incar: bool = True - CONFIG = MPRelaxSet.CONFIG - - @property - def incar_updates(self) -> dict: - """Updates to the INCAR config for this calculation type.""" - updates: dict[str, Any] = {"NSW": 0, "ISMEAR": -5, "LCHARG": True, "LORBIT": 11, "LREAL": False} - if self.lepsilon: - # LPEAD=T: numerical evaluation of overlap integral prevents LRF_COMMUTATOR - # errors and can lead to better expt. agreement but produces slightly - # different results - updates.update({"IBRION": 8, "LEPSILON": True, "LPEAD": True, "NSW": 1, "EDIFF": 1e-5}) - - if self.lcalcpol: - updates["LCALCPOL"] = True - return updates - - @property - def kpoints_updates(self) -> dict | Kpoints: - """Updates to the kpoints configuration for this calculation type.""" - factor = 1.0 - if self.bandgap is not None and self.small_gap_multiply and self.bandgap <= self.small_gap_multiply[0]: - factor = self.small_gap_multiply[1] - - # prefer to use k-point scheme from previous run unless lepsilon = True is specified - if self.prev_kpoints and self.prev_kpoints.style == Kpoints.supported_modes.Monkhorst and not self.lepsilon: # type: ignore - kpoints = Kpoints.automatic_density_by_vol( - self.structure, # type: ignore - int(self.reciprocal_density * factor), - self.force_gamma, - ) - k_div = [kp + 1 if kp % 2 == 1 else kp for kp in kpoints.kpts[0]] # type: ignore - return Kpoints.monkhorst_automatic(k_div) # type: ignore - - return {"reciprocal_density": self.reciprocal_density * factor} - - -@dataclass -class MatPESStaticSet(VaspInputSet): - """Create input files for a MatPES static calculation. - - The goal of MatPES is to generate potential energy surface data. This is a distinctly different - from the objectives of the MP static calculations, which aims to obtain primarily accurate - energies and also electronic structure (DOS). For PES data, force accuracy (and to some extent, - stress accuracy) is of paramount importance. - - The default POTCAR versions have been updated to PBE_54 from the old PBE set used in the - MPStaticSet. However, **U values** are still based on PBE. The implicit assumption here is that - the PBE_54 and PBE POTCARs are sufficiently similar that the U values fitted to the old PBE - functional still applies. - - Args: - structure (Structure): The Structure to create inputs for. If None, the input - set is initialized without a Structure but one must be set separately before - the inputs are generated. - xc_functional ('R2SCAN'|'PBE'): Exchange-correlation functional to use. Defaults to 'PBE'. - **kwargs: Keywords supported by VaspInputSet. - """ - - xc_functional: Literal["R2SCAN", "PBE", "PBE+U"] = "PBE" - prev_incar: dict | str | None = None - # These are parameters that we will inherit from any previous INCAR supplied. They are mostly parameters related - # to symmetry and convergence set by Custodian when errors are encountered in a previous run. Given that our goal - # is to have a strictly homogeneous PES data, all other parameters (e.g., ISMEAR, ALGO, etc.) are not inherited. - inherit_incar: list[str] | bool = ( # type: ignore - "LPEAD", - "NGX", - "NGY", - "NGZ", - "SYMPREC", - "IMIX", - "LMAXMIX", - "KGAMMA", - "ISYM", - "NCORE", - "NPAR", - "NELMIN", - "IOPT", - "NBANDS", - "KPAR", - "AMIN", - "NELMDL", - "BMIX", - "AMIX_MAG", - "BMIX_MAG", - ) - CONFIG = _load_yaml_config("MatPESStaticSet") - - def __post_init__(self): - """Validate inputs""" - super().__post_init__() - valid_xc_functionals = ("R2SCAN", "PBE", "PBE+U") - if self.xc_functional.upper() not in valid_xc_functionals: - raise ValueError( - f"Unrecognized xc_functional='{self.xc_functional}'. " - f"Supported exchange-correlation functionals are {valid_xc_functionals}" - ) - - default_potcars = self.CONFIG["PARENT"].replace("PBE", "PBE_").replace("Base", "") # PBE64Base -> PBE_64 - self.user_potcar_functional = self.user_potcar_functional or default_potcars - if self.user_potcar_functional.upper() != default_potcars: - warnings.warn( - f"{self.user_potcar_functional=} is inconsistent with the recommended {default_potcars}.", UserWarning - ) - - if self.xc_functional.upper() == "R2SCAN": - self._config_dict["INCAR"].update({"METAGGA": "R2SCAN", "ALGO": "ALL", "GGA": None}) - if self.xc_functional.upper().endswith("+U"): - self._config_dict["INCAR"]["LDAU"] = True - - -@dataclass -class MPScanStaticSet(MPScanRelaxSet): - """Create input files for a static calculation using the accurate and numerically - efficient r2SCAN variant of the Strongly Constrained and Appropriately Normed - (SCAN) metaGGA functional. - - Args: - structure (Structure): Structure from previous run. - bandgap (float): Bandgap of the structure in eV. The bandgap is used to - compute the appropriate k-point density and determine the smearing settings. - lepsilon (bool): Whether to add static dielectric calculation - lcalcpol (bool): Whether to turn on evaluation of the Berry phase approximations - for electronic polarization. - **kwargs: Keywords supported by MPScanRelaxSet. - """ - - lepsilon: bool = False - lcalcpol: bool = False - inherit_incar: bool = True - auto_kspacing: bool = True - - @property - def incar_updates(self) -> dict: - """Updates to the INCAR config for this calculation type.""" - updates: dict[str, Any] = { - "LREAL": False, - "NSW": 0, - "LORBIT": 11, - "LVHAR": True, - "ISMEAR": -5, - } - - if self.lepsilon: - # LPEAD=T: numerical evaluation of overlap integral prevents - # LRF_COMMUTATOR errors and can lead to better expt. agreement - # but produces slightly different results - updates.update({"IBRION": 8, "LEPSILON": True, "LPEAD": True, "NSW": 1, "NPAR": None}) - - if self.lcalcpol: - updates["LCALCPOL"] = True - - return updates - - -@dataclass -class MPHSEBSSet(VaspInputSet): - """ - Implementation of a VaspInputSet for HSE band structure computations. - - Remember that HSE band structures must be self-consistent in VASP. A band structure - along symmetry lines for instance needs BOTH a uniform grid with appropriate weights - AND a path along the lines with weight 0. - - Thus, the "uniform" mode is just like regular static SCF but allows adding custom - kpoints (e.g., corresponding to known VBM/CBM) to the uniform grid that have zero - weight (e.g., for better gap estimate). - - The "gap" mode behaves just like the "uniform" mode, however, if starting from a - previous calculation, the VBM and CBM k-points will automatically be added to - ``added_kpoints``. - - The "line" mode is just like "uniform" mode, but additionally adds k-points along - symmetry lines with zero weight. - - The "uniform_dense" mode is like "uniform" mode but additionally adds a denser - uniform mesh with zero weight. This can be useful when calculating Fermi surfaces - or BoltzTraP/AMSET electronic transport using hybrid DFT. - - Args: - structure (Structure): Structure to compute - added_kpoints (list): a list of kpoints (list of 3 number list) added to the - run. The k-points are in fractional coordinates - mode (str): "Line" - generate k-points along symmetry lines for bandstructure. - "Uniform" - generate uniform k-points grid. - reciprocal_density (int): k-point density to use for uniform mesh. - copy_chgcar (bool): Whether to copy the CHGCAR of a previous run. - kpoints_line_density (int): k-point density for high symmetry lines - dedos (float): Energy difference used to set NEDOS, based on the total energy - range. - optics (bool): Whether to add LOPTICS (used for calculating optical response). - nbands_factor (float): Multiplicative factor for NBANDS when starting from a - previous calculation. Choose a higher number if you are doing an LOPTICS - calculation. - **kwargs: Keywords supported by VaspInputSet. - """ - - added_kpoints: list[Vector3D] = field(default_factory=list) - mode: str = "gap" - reciprocal_density: float = 50 - copy_chgcar: bool = True - kpoints_line_density: float = 20 - nbands_factor: float = 1.2 - zero_weighted_reciprocal_density: float = 100 - dedos: float = 0.02 - optics: bool = False - CONFIG = MPHSERelaxSet.CONFIG - - def __post_init__(self) -> None: - """Ensure mode is set correctly.""" - super().__post_init__() - - if "reciprocal_density" in self.user_kpoints_settings: - self.reciprocal_density = self.user_kpoints_settings["reciprocal_density"] - - self.mode = self.mode.lower() - supported_modes = ("line", "uniform", "gap", "uniform_dense") - if self.mode not in supported_modes: - raise ValueError(f"Supported modes are: {', '.join(supported_modes)}") - - @property - def kpoints_updates(self) -> dict: - """Updates to the kpoints configuration for this calculation type.""" - kpoints: dict[str, Any] = {"reciprocal_density": self.reciprocal_density, "explicit": True} - - if self.mode == "line": - # add line_density on top of reciprocal density - kpoints["zero_weighted_line_density"] = self.kpoints_line_density - - elif self.mode == "uniform_dense": - kpoints["zero_weighted_reciprocal_density"] = self.zero_weighted_reciprocal_density - - added_kpoints = deepcopy(self.added_kpoints) - if self.prev_vasprun is not None and self.mode == "gap": - bs = self.prev_vasprun.get_band_structure() - if not bs.is_metal(): - added_kpoints.append(bs.get_vbm()["kpoint"].frac_coords) - added_kpoints.append(bs.get_cbm()["kpoint"].frac_coords) - - kpoints["added_kpoints"] = added_kpoints - - return kpoints - - @property - def incar_updates(self) -> dict: - """Updates to the INCAR config for this calculation type.""" - updates = dict(NSW=0, ISMEAR=0, SIGMA=0.05, ISYM=3, LCHARG=False, NELMIN=5) - - if self.mode == "uniform" and len(self.added_kpoints) == 0: - # automatic setting of nedos using the energy range and the energy step - nedos = _get_nedos(self.prev_vasprun, self.dedos) - - # use tetrahedron method for DOS and optics calculations - updates.update({"ISMEAR": -5, "NEDOS": nedos}) - - else: - # if line mode or explicit k-points (gap) can't use ISMEAR=-5 - # use small sigma to avoid partial occupancies for small band gap materials - updates.update({"ISMEAR": 0, "SIGMA": 0.01}) - - if self.prev_vasprun is not None: - # set nbands - nbands = int(np.ceil(self.prev_vasprun.parameters["NBANDS"] * self.nbands_factor)) - updates["NBANDS"] = nbands - - if self.optics: - # LREAL not supported with LOPTICS - updates.update({"LOPTICS": True, "LREAL": False, "CSHIFT": 1e-5}) - - if self.prev_vasprun is not None and self.prev_outcar is not None: - # turn off spin when magmom for every site is smaller than 0.02. - updates["ISPIN"] = _get_ispin(self.prev_vasprun, self.prev_outcar) - - return updates - - -@dataclass -class MPNonSCFSet(VaspInputSet): - """ - Init a MPNonSCFSet. Typically, you would use the classmethod - from_prev_calc to initialize from a previous SCF run. - - Args: - structure (Structure): Structure to compute - mode (str): Line, Uniform or Boltztrap mode supported. - nedos (int): nedos parameter. Default to 2001. - dedos (float): setting nedos=0 and uniform mode in from_prev_calc, - an automatic nedos will be calculated using the total energy range - divided by the energy step dedos - reciprocal_density (int): density of k-mesh by reciprocal - volume (defaults to 100) - kpoints_line_density (int): Line density for Line mode. - optics (bool): whether to add dielectric function - copy_chgcar: Whether to copy the old CHGCAR when starting from a - previous calculation. - nbands_factor (float): Multiplicative factor for NBANDS when starting - from a previous calculation. Choose a higher number if you are - doing an LOPTICS calculation. - small_gap_multiply ([float, float]): When starting from a previous - calculation, if the gap is less than 1st index, multiply the default - reciprocal_density by the 2nd index. - **kwargs: Keywords supported by MPRelaxSet. - """ - - mode: str = "line" - nedos: int = 2001 - dedos: float = 0.005 - reciprocal_density: float = 100 - kpoints_line_density: float = 20 - optics: bool = False - copy_chgcar: bool = True - nbands_factor: float = 1.2 - small_gap_multiply: tuple[float, float] | None = None - inherit_incar: bool = True - CONFIG = MPRelaxSet.CONFIG - - def __post_init__(self): - """Perform inputset validation.""" - super().__post_init__() - - mode = self.mode = self.mode.lower() - - valid_modes = ("line", "uniform", "boltztrap") - if mode not in valid_modes: - raise ValueError( - f"Invalid {mode=}. Supported modes for NonSCF runs are {', '.join(map(repr, valid_modes))}" - ) - - if (mode != "uniform" or self.nedos < 2000) and self.optics: - warnings.warn("It is recommended to use Uniform mode with a high NEDOS for optics calculations.") - - if self.standardize: - warnings.warn( - "Use of standardize=True with from_prev_run is not " - "recommended as there is no guarantee the copied " - "files will be appropriate for the standardized" - " structure. copy_chgcar is enforced to be false." - ) - self.copy_chgcar = False - - @property - def incar_updates(self) -> dict: - """Updates to the INCAR config for this calculation type.""" - updates: dict[str, Any] = {"LCHARG": False, "LORBIT": 11, "LWAVE": False, "NSW": 0, "ISYM": 0, "ICHARG": 11} - - if self.prev_vasprun is not None: - # set NBANDS - n_bands = int(np.ceil(self.prev_vasprun.parameters["NBANDS"] * self.nbands_factor)) - updates["NBANDS"] = n_bands - - # automatic setting of NEDOS using the energy range and the energy step - nedos = _get_nedos(self.prev_vasprun, self.dedos) if self.nedos == 0 else self.nedos - - if self.mode == "uniform": - # use tetrahedron method for DOS and optics calculations - updates.update({"ISMEAR": -5, "ISYM": 2, "NEDOS": nedos}) - - elif self.mode in ("line", "boltztrap"): - # if line mode or explicit k-points (boltztrap) can't use ISMEAR=-5 - # use small sigma to avoid partial occupancies for small band gap materials - # use a larger sigma if the material is a metal - sigma = 0.2 if self.bandgap == 0 or self.bandgap is None else 0.01 - updates.update({"ISMEAR": 0, "SIGMA": sigma}) - - if self.optics: - # LREAL not supported with LOPTICS = True; automatic NEDOS usually - # underestimates, so set it explicitly - updates.update({"LOPTICS": True, "LREAL": False, "CSHIFT": 1e-5, "NEDOS": nedos}) - - if self.prev_vasprun is not None and self.prev_outcar is not None: - # turn off spin when magmom for every site is smaller than 0.02. - updates["ISPIN"] = _get_ispin(self.prev_vasprun, self.prev_outcar) - - updates["MAGMOM"] = None - return updates - - @property - def kpoints_updates(self) -> dict: - """Updates to the kpoints configuration for this calculation type.""" - factor = 1.0 - if self.bandgap is not None and self.small_gap_multiply and self.bandgap <= self.small_gap_multiply[0]: - factor = self.small_gap_multiply[1] - - if self.mode == "line": - return {"line_density": self.kpoints_line_density * factor} - - if self.mode == "boltztrap": - return {"explicit": True, "reciprocal_density": self.reciprocal_density * factor} - - return {"reciprocal_density": self.reciprocal_density * factor} - - -@dataclass -class MPSOCSet(VaspInputSet): - """An input set for running spin-orbit coupling (SOC) calculations. - - Args: - structure (Structure): the structure must have the 'magmom' site - property and each magnetic moment value must have 3 - components. eg: ``magmom = [[0,0,2], ...]`` - saxis (tuple): magnetic moment orientation - copy_chgcar: Whether to copy the old CHGCAR. Defaults to True. - nbands_factor (float): Multiplicative factor for NBANDS. Choose a - higher number if you are doing an LOPTICS calculation. - reciprocal_density (int): density of k-mesh by reciprocal volume. - small_gap_multiply ([float, float]): If the gap is less than - 1st index, multiply the default reciprocal_density by the 2nd - index. - lepsilon (bool): Whether to add static dielectric calculation - lcalcpol (bool): Whether to turn on evaluation of the Berry phase approximations - for electronic polarization - magmom (list[list[float]]): Override for the structure magmoms. - **kwargs: Keywords supported by VaspInputSet. - """ - - saxis: tuple[int, int, int] = (0, 0, 1) - nbands_factor: float = 1.2 - lepsilon: bool = False - lcalcpol: bool = False - reciprocal_density: float = 100 - small_gap_multiply: tuple[float, float] | None = None - magmom: list[Vector3D] | None = None - inherit_incar: bool = True - copy_chgcar: bool = True - CONFIG = MPRelaxSet.CONFIG - - def __post_init__(self): - super().__post_init__() - if ( - self.structure - and not hasattr(self.structure[0], "magmom") - and not isinstance(self.structure[0].magmom, list) - ): - raise ValueError( - "The structure must have the 'magmom' site property and each magnetic " - "moment value must have 3 components. e.g. magmom = [0,0,2]" - ) - - @property - def incar_updates(self) -> dict: - """Updates to the INCAR config for this calculation type.""" - updates = { - "ISYM": -1, - "LSORBIT": "T", - "ICHARG": 11, - "SAXIS": list(self.saxis), - "NSW": 0, - "ISMEAR": -5, - "LCHARG": True, - "LORBIT": 11, - "LREAL": False, - } - - if self.lepsilon: - # LPEAD=T: numerical evaluation of overlap integral prevents LRF_COMMUTATOR - # errors and can lead to better expt. agreement but produces slightly - # different results - updates.update({"IBRION": 8, "LEPSILON": True, "LPEAD": True, "NSW": 1}) - - if self.lcalcpol: - updates["LCALCPOL"] = True - - if self.prev_vasprun is not None: - # set NBANDS - n_bands = int(np.ceil(self.prev_vasprun.parameters["NBANDS"] * self.nbands_factor)) - updates["NBANDS"] = n_bands - return updates - - @property - def kpoints_updates(self) -> dict: - """Updates to the kpoints configuration for this calculation type.""" - factor = 1.0 - if self.bandgap is not None and self.small_gap_multiply and self.bandgap <= self.small_gap_multiply[0]: - factor = self.small_gap_multiply[1] - return {"reciprocal_density": self.reciprocal_density * factor} - - @VaspInputSet.structure.setter # type: ignore - def structure(self, structure: Structure | None) -> None: - if structure is not None: - if self.magmom: - structure = structure.copy(site_properties={"magmom": self.magmom}) - - # magmom has to be 3D for SOC calculation. - if hasattr(structure[0], "magmom"): - if not isinstance(structure[0].magmom, list): - # project magmom to z-axis - structure = structure.copy(site_properties={"magmom": [[0, 0, site.magmom] for site in structure]}) - else: - raise ValueError("Neither the previous structure has magmom property nor magmom provided") - - VaspInputSet.structure.fset(self, structure) # type: ignore - - -@dataclass -class MPNMRSet(VaspInputSet): - """Init a MPNMRSet. - - Args: - structure (Structure): Structure from previous run. - mode (str): The NMR calculation to run - "cs": for Chemical Shift - "efg" for Electric Field Gradient - isotopes (list): list of Isotopes for quadrupole moments - reciprocal_density (int): density of k-mesh by reciprocal volume. - lepsilon (bool): Whether to add static dielectric calculation - lcalcpol (bool): Whether to turn on evaluation of the Berry phase approximations - for electronic polarization - reciprocal_density (int): For static calculations, we usually set the - reciprocal density by volume. This is a convenience arg to change - that, rather than using user_kpoints_settings. Defaults to 100, - which is ~50% more than that of standard relaxation calculations. - small_gap_multiply ([float, float]): If the gap is less than - 1st index, multiply the default reciprocal_density by the 2nd - index. - **kwargs: Keywords supported by MPRelaxSet. - """ - - mode: Literal["cs", "efg"] = "cs" - isotopes: list = field(default_factory=list) - reciprocal_density: int = 100 - small_gap_multiply: tuple[float, float] | None = None - inherit_incar: bool = True - CONFIG = MPRelaxSet.CONFIG - - @property - def incar_updates(self) -> dict: - """Updates to the INCAR config for this calculation type.""" - updates: dict[str, Any] = {"NSW": 0, "ISMEAR": -5, "LCHARG": True, "LORBIT": 11, "LREAL": False} - if self.mode.lower() == "cs": - updates.update( - LCHIMAG=True, - EDIFF=-1.0e-10, - ISYM=0, - LCHARG=False, - LNMR_SYM_RED=True, - NELMIN=10, - NLSPLINE=True, - PREC="ACCURATE", - SIGMA=0.01, - ) - elif self.mode.lower() == "efg": - isotopes = {ist.split("-")[0]: ist for ist in self.isotopes} - quad_efg = [ - float(Species(s.name).get_nmr_quadrupole_moment(isotopes.get(s.name))) - for s in self.structure.species # type: ignore - ] - updates.update( - ALGO="FAST", - EDIFF=-1.0e-10, - ISYM=0, - LCHARG=False, - LEFG=True, - QUAD_EFG=quad_efg, - NELMIN=10, - PREC="ACCURATE", - SIGMA=0.01, - ) - return updates - - @property - def kpoints_updates(self) -> dict: - """Updates to the kpoints configuration for this calculation type.""" - factor = 1.0 - if self.bandgap is not None and self.small_gap_multiply and self.bandgap <= self.small_gap_multiply[0]: - factor = self.small_gap_multiply[1] - return {"reciprocal_density": self.reciprocal_density * factor} - - -@due.dcite( - Doi("10.1149/2.0061602jes"), - description="Elastic Properties of Alkali Superionic Conductor Electrolytes from First Principles Calculations", -) -class MVLElasticSet(VaspInputSet): - """ - MVL denotes VASP input sets that are implemented by the Materials Virtual - Lab (http://materialsvirtuallab.org) for various research. - - This input set is used to calculate elastic constants in VASP. It is used - in the following work:: - - Z. Deng, Z. Wang, I.-H. Chu, J. Luo, S. P. Ong. - “Elastic Properties of Alkali Superionic Conductor Electrolytes - from First Principles Calculations”, J. Electrochem. Soc. - 2016, 163(2), A67-A74. doi: 10.1149/2.0061602jes - - To read the elastic constants, you may use the Outcar class which parses the - elastic constants. - - Args: - structure (pymatgen.Structure): Input structure. - potim (float): POTIM parameter. The default of 0.015 is usually fine, - but some structures may require a smaller step. - **kwargs: Parameters supported by MPRelaxSet. - """ - - potim: float = 0.015 - CONFIG = MPRelaxSet.CONFIG - - @property - def incar_updates(self) -> dict: - """Updates to the INCAR config for this calculation type.""" - return {"IBRION": 6, "NFREE": 2, "POTIM": self.potim, "NPAR": None} - - -@dataclass -class MVLGWSet(VaspInputSet): - """ - MVL denotes VASP input sets that are implemented by the Materials Virtual - Lab (http://materialsvirtuallab.org) for various research. This is a - flexible input set for GW calculations. - - Note that unlike all other input sets in this module, the PBE_54 series of - functional is set as the default. These have much improved performance for - GW calculations. - - A typical sequence is mode="STATIC" -> mode="DIAG" -> mode="GW" -> - mode="BSE". For all steps other than the first one (static), the - recommendation is to use from_prev_calculation on the preceding run in - the series. - - Args: - structure (Structure): Input structure. - mode (str): Supported modes are "STATIC" (default), "DIAG", "GW", - and "BSE". - 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 GW and BSE. This parameter is used by - from_previous_calculation to set nband. - copy_wavecar: Whether to copy the old WAVECAR, WAVEDER and associated - files when starting from a previous calculation. - nbands_factor (int): Multiplicative factor for NBANDS when starting - from a previous calculation. Only applies if mode=="DIAG". - Need to be tested for convergence. - reciprocal_density (int): Density of k-mesh by reciprocal atom. Only - applies if mode=="STATIC". Defaults to 100. - ncores (int): Numbers of cores used for the calculation. VASP will alter - NBANDS if it was not dividable by ncores. Only applies if - mode=="DIAG". - **kwargs: All kwargs supported by VaspInputSet. Typically, - user_incar_settings is a commonly used option. - """ - - reciprocal_density: float = 100 - mode: str = "STATIC" - copy_wavecar: bool = True - nbands_factor: int = 5 - ncores: int = 16 - nbands: int | None = None - force_gamma: bool = True - inherit_incar: bool = True # inherit incar from previous run if available - SUPPORTED_MODES = ("DIAG", "GW", "STATIC", "BSE") - CONFIG = _load_yaml_config("MVLGWSet") - - def __post_init__(self): - """Validate input settings.""" - super().__post_init__() - self.mode = mode = self.mode.upper() - - if mode not in MVLGWSet.SUPPORTED_MODES: - raise ValueError(f"Invalid {mode=}, supported modes are {', '.join(map(repr, MVLGWSet.SUPPORTED_MODES))}") - - @property - def kpoints_updates(self) -> dict: - """Updates to the kpoints configuration for this calculation type.""" - # Generate gamma center k-points mesh grid for GW calc, which is requested - # by GW calculation. - return {"reciprocal_density": self.reciprocal_density} - - @property - def incar_updates(self) -> dict: - """Updates to the INCAR config for this calculation type.""" - updates = {} - nbands = int(self.prev_vasprun.parameters["NBANDS"]) if self.prev_vasprun is not None else None - - if self.mode == "DIAG": - # Default parameters for diagonalization calculation. - updates.update({"ALGO": "Exact", "NELM": 1, "LOPTICS": True, "LPEAD": True}) - if nbands: - nbands = int(np.ceil(nbands * self.nbands_factor / self.ncores) * self.ncores) - - elif self.mode == "GW": - # Default parameters for GW calculation. - updates.update( - {"ALGO": "GW0", "NELM": 1, "NOMEGA": 80, "ENCUTGW": 250, "EDIFF": None, "LOPTICS": None, "LPEAD": None} - ) - elif self.mode == "BSE": - # Default parameters for BSE calculation. - updates.update({"ALGO": "BSE", "ANTIRES": 0, "NBANDSO": 20, "NBANDSV": 20}) - - if nbands: - updates["NBANDS"] = nbands - - return updates - - @classmethod - def from_prev_calc(cls, prev_calc_dir: str, mode: str = "DIAG", **kwargs) -> Self: - """Generate a set of VASP input files for GW or BSE calculations from a - directory of previous Exact Diag VASP run. - - Args: - prev_calc_dir (str): The directory contains the outputs( - vasprun.xml of previous vasp run. - mode (str): Supported modes are "STATIC", "DIAG" (default), "GW", - and "BSE". - **kwargs: All kwargs supported by MVLGWSet, other than structure, - prev_incar and mode, which are determined from the - prev_calc_dir. - """ - input_set = cls(_dummy_structure, mode=mode, **kwargs) - return input_set.override_from_prev_calc(prev_calc_dir=prev_calc_dir) - - -@dataclass -class MVLSlabSet(VaspInputSet): - """Write a set of slab vasp runs, including both slabs (along the c direction) - and orient unit cells (bulk), to ensure the same KPOINTS, POTCAR and INCAR criterion. - - Args: - structure: Structure - k_product: default to 50, kpoint number * length for a & b - directions, also for c direction in bulk calculations - bulk: - auto_dipole: - set_mix: - sort_structure: - **kwargs: Other kwargs supported by VaspInputSet. - """ - - k_product: int = 50 - bulk: bool = False - auto_dipole: bool = False - set_mix: bool = True - CONFIG = MPRelaxSet.CONFIG - - @property - def incar_updates(self) -> dict: - """Updates to the INCAR config for this calculation type.""" - updates = {"EDIFF": 1e-4, "EDIFFG": -0.02, "ENCUT": 400, "ISMEAR": 0, "SIGMA": 0.05, "ISIF": 3} - if not self.bulk: - updates.update({"ISIF": 2, "LVTOT": True, "NELMIN": 8}) - if self.set_mix: - updates.update({"AMIN": 0.01, "AMIX": 0.2, "BMIX": 0.001}) - if self.auto_dipole: - weights = [s.species.weight for s in self.structure] # type: ignore - center_of_mass = np.average(self.structure.frac_coords, weights=weights, axis=0) # type: ignore - updates.update({"IDIPOL": 3, "LDIPOL": True, "DIPOL": center_of_mass}) - return updates - - @property - def kpoints_updates(self): - """Updates to the kpoints configuration for this calculation type. - - k_product, default to 50, is kpoint number * length for a & b - directions, also for c direction in bulk calculations - Automatic mesh & Gamma is the default setting. - """ - # To get input sets, the input structure has to has the same number - # of required parameters as a Structure object (ie. 4). Slab - # attributes aren't going to affect the VASP inputs anyways so - # converting the slab into a structure should not matter - # use k_product to calculate kpoints, k_product = kpts[0][0] * a - lattice_abc = self.structure.lattice.abc - kpt_calc = [ - int(self.k_product / lattice_abc[0] + 0.5), - int(self.k_product / lattice_abc[1] + 0.5), - 1, - ] - - # calculate kpts (c direction) for bulk. (for slab, set to 1) - if self.bulk: - kpt_calc[2] = int(self.k_product / lattice_abc[2] + 0.5) - - return Kpoints(comment="Generated by pymatgen's MVLGBSet", style=Kpoints.supported_modes.Gamma, kpts=[kpt_calc]) - - def as_dict(self, verbosity=2): - """ - Args: - verbosity (int): Verbosity of dict. e.g. whether to include Structure. - - Returns: - dict: MSONable MVLGBSet representation. - """ - dct = MSONable.as_dict(self) - if verbosity == 1: - dct.pop("structure", None) - return dct - - -@dataclass -class MVLGBSet(VaspInputSet): - """Write a vasp input files for grain boundary calculations, slab or bulk. - - Args: - structure (Structure): provide the structure - k_product: Kpoint number * length for a & b directions, also for c direction in - bulk calculations. Default to 40. - slab_mode (bool): Defaults to False. Use default (False) for a bulk supercell. - Use True if you are performing calculations on a slab-like (i.e., surface) - of the GB, for example, when you are calculating the work of separation. - is_metal (bool): Defaults to True. This determines whether an ISMEAR of 1 is - used (for metals) or not (for insulators and semiconductors) by default. - Note that it does *not* override user_incar_settings, which can be set by - the user to be anything desired. - **kwargs: - Other kwargs supported by MPRelaxSet. - """ - - k_product: int = 40 - slab_mode: bool = False - is_metal: bool = True - CONFIG = MPRelaxSet.CONFIG - - @property - def kpoints_updates(self): - """k_product is kpoint number * length for a & b directions, also for c direction - in bulk calculations Automatic mesh & Gamma is the default setting. - """ - # use k_product to calculate kpoints, k_product = kpts[0][0] * a - lengths = self.structure.lattice.abc - kpt_calc = [ - int(self.k_product / lengths[0] + 0.5), - int(self.k_product / lengths[1] + 0.5), - int(self.k_product / lengths[2] + 0.5), - ] - - if self.slab_mode: - kpt_calc[2] = 1 - - return Kpoints(comment="Generated by pymatgen's MVLGBSet", style=Kpoints.supported_modes.Gamma, kpts=[kpt_calc]) - - @property - def incar_updates(self) -> dict: - """Updates to the INCAR config for this calculation type.""" - # The default incar setting is used for metallic system, for - # insulator or semiconductor, ISMEAR need to be changed. - updates = dict(LCHARG=False, NELM=60, PREC="Normal", EDIFFG=-0.02, ICHARG=0, NSW=200, EDIFF=0.0001) - - if self.is_metal: - updates["ISMEAR"] = 1 - updates["LDAU"] = False - - if self.slab_mode: - # for clean grain boundary and bulk relaxation, full optimization - # relaxation (ISIF=3) is used. For slab relaxation (ISIF=2) is used. - updates["ISIF"] = 2 - updates["NELMIN"] = 8 - - return updates - - -@dataclass -class MVLRelax52Set(VaspInputSet): - """ - Implementation of VaspInputSet utilizing the public Materials Project - parameters for INCAR & KPOINTS and VASP's recommended PAW potentials for - POTCAR. - - Keynotes from VASP manual: - 1. Recommended potentials for calculations using vasp.5.2+ - 2. If dimers with short bonds are present in the compound (O2, CO, - N2, F2, P2, S2, Cl2), it is recommended to use the h potentials. - Specifically, C_h, O_h, N_h, F_h, P_h, S_h, Cl_h - 3. Released on Oct 28, 2018 by VASP. Please refer to VASP - Manual 1.2, 1.3 & 10.2.1 for more details. - - Args: - structure (Structure): input structure. - user_potcar_functional (str): choose from "PBE_52" and "PBE_54". - **kwargs: Other kwargs supported by VaspInputSet. - """ - - user_potcar_functional: UserPotcarFunctional = "PBE_52" - CONFIG = _load_yaml_config("MVLRelax52Set") - _valid_potcars: Sequence[str] | None = ("PBE_52", "PBE_54") - - -class MITNEBSet(VaspInputSet): - """Write NEB inputs. - - Note that EDIFF is not on a per atom basis for this input set. - """ - - def __init__(self, structures, unset_encut=False, **kwargs) -> None: - """ - Args: - structures: List of Structure objects. - unset_encut (bool): Whether to unset ENCUT. - **kwargs: Other kwargs supported by VaspInputSet. - """ - if len(structures) < 3: - raise ValueError(f"You need at least 3 structures for an NEB, got {len(structures)}") - kwargs["sort_structure"] = False - super().__init__(structures[0], MITRelaxSet.CONFIG, **kwargs) - self.structures = self._process_structures(structures) - - self.unset_encut = False - if unset_encut: - self._config_dict["INCAR"].pop("ENCUT", None) - - if "EDIFF" not in self._config_dict["INCAR"]: - self._config_dict["INCAR"]["EDIFF"] = self._config_dict["INCAR"].pop("EDIFF_PER_ATOM") - - # NEB specific defaults - defaults = {"IMAGES": len(structures) - 2, "IBRION": 1, "ISYM": 0, "LCHARG": False, "LDAU": False} - self._config_dict["INCAR"].update(defaults) - - @property - def poscar(self): - """Poscar for structure of first end point.""" - return Poscar(self.structures[0]) - - @property - def poscars(self): - """List of Poscars.""" - return [Poscar(s) for s in self.structures] - - @staticmethod - def _process_structures(structures): - """Remove any atom jumps across the cell.""" - input_structures = structures - structures = [input_structures[0]] - for s in input_structures[1:]: - prev = structures[-1] - for idx, site in enumerate(s): - translate = np.round(prev[idx].frac_coords - site.frac_coords) - if np.any(np.abs(translate) > 0.5): - s.translate_sites([idx], translate, to_unit_cell=False) - structures.append(s) - return structures - - def write_input( - self, - output_dir, - make_dir_if_not_present=True, - write_cif=False, - write_path_cif=False, - write_endpoint_inputs=False, - ): - """ - NEB inputs has a special directory structure where inputs are in 00, - 01, 02, .... - - Args: - output_dir (str): Directory to output the VASP input files - make_dir_if_not_present (bool): Set to True if you want the - directory (and the whole path) to be created if it is not - present. - write_cif (bool): If true, writes a cif along with each POSCAR. - write_path_cif (bool): If true, writes a cif for each image. - write_endpoint_inputs (bool): If true, writes input files for - running endpoint calculations. - """ - output_dir = Path(output_dir) - if make_dir_if_not_present and not output_dir.exists(): - output_dir.mkdir(parents=True) - self.incar.write_file(str(output_dir / "INCAR")) - self.kpoints.write_file(str(output_dir / "KPOINTS")) - self.potcar.write_file(str(output_dir / "POTCAR")) - - for idx, poscar in enumerate(self.poscars): - d = output_dir / str(idx).zfill(2) - if not d.exists(): - d.mkdir(parents=True) - poscar.write_file(str(d / "POSCAR")) - if write_cif: - poscar.structure.to(filename=str(d / f"{idx}.cif")) - if write_endpoint_inputs: - end_point_param = MITRelaxSet(self.structures[0], user_incar_settings=self.user_incar_settings) - - for image in ["00", str(len(self.structures) - 1).zfill(2)]: - end_point_param.incar.write_file(str(output_dir / image / "INCAR")) - end_point_param.kpoints.write_file(str(output_dir / image / "KPOINTS")) - end_point_param.potcar.write_file(str(output_dir / image / "POTCAR")) - if write_path_cif: - sites = { - PeriodicSite(site.species, site.frac_coords, self.structures[0].lattice) - for site in chain(*(struct for struct in self.structures)) - } - neb_path = Structure.from_sites(sorted(sites)) - neb_path.to(filename=f"{output_dir}/path.cif") - - -@dataclass -class MITMDSet(VaspInputSet): - """Write a VASP MD run. This DOES NOT do multiple stage runs. - - Args: - structure (Structure): Input structure. - start_temp (float): Starting temperature. - end_temp (float): Final temperature. - nsteps (int): Number of time steps for simulations. NSW parameter. - time_step (float): The time step for the simulation. The POTIM - parameter. Defaults to 2fs. - spin_polarized (bool): Whether to do spin polarized calculations. - The ISPIN parameter. Defaults to False. - **kwargs: Other kwargs supported by VaspInputSet. - """ - - structure: Structure | None = None - start_temp: float = 0.0 - end_temp: float = 300.0 - nsteps: int = 1000 - time_step: float = 2 - spin_polarized: bool = False - CONFIG = MITRelaxSet.CONFIG - - @property - def incar_updates(self): - """Updates to the INCAR config for this calculation type.""" - # MD default settings - return { - "TEBEG": self.start_temp, - "TEEND": self.end_temp, - "NSW": self.nsteps, - "EDIFF_PER_ATOM": 0.000001, - "LSCALU": False, - "LCHARG": False, - "LPLANE": False, - "LWAVE": True, - "ISMEAR": 0, - "NELMIN": 4, - "LREAL": True, - "BMIX": 1, - "MAXMIX": 20, - "NELM": 500, - "NSIM": 4, - "ISYM": 0, - "ISIF": 0, - "IBRION": 0, - "NBLOCK": 1, - "KBLOCK": 100, - "SMASS": 0, - "POTIM": self.time_step, - "PREC": "Low", - "ISPIN": 2 if self.spin_polarized else 1, - "LDAU": False, - "ENCUT": None, - } - - @property - def kpoints_updates(self) -> Kpoints | dict: - """Updates to the kpoints configuration for this calculation type.""" - return Kpoints.gamma_automatic() - - -@dataclass -class MPMDSet(VaspInputSet): - """ - This a modified version of the old MITMDSet pre 2018/03/12. - - This set serves as the basis for the amorphous skyline paper. - - (1) Aykol, M.; Dwaraknath, S. S.; Sun, W.; Persson, K. A. Thermodynamic - Limit for Synthesis of Metastable Inorganic Materials. Sci. Adv. 2018, - 4 (4). - - Class for writing a VASP MD run. This DOES NOT do multiple stage runs. - Precision remains normal, to increase accuracy of stress tensor. - - Args: - structure (Structure): Input structure. - start_temp (int): Starting temperature. - end_temp (int): Final temperature. - nsteps (int): Number of time steps for simulations. NSW parameter. - time_step (float): The time step for the simulation. The POTIM - parameter. Defaults to None, which will set it automatically - to 2.0 fs for non-hydrogen containing structures and 0.5 fs - for hydrogen containing structures. - spin_polarized (bool): Whether to do spin polarized calculations. - The ISPIN parameter. Defaults to False. - **kwargs: Other kwargs supported by VaspInputSet. - """ - - start_temp: float = 0.0 - end_temp: float = 300.0 - nsteps: int = 1000 - time_step: float | None = None - spin_polarized: bool = False - CONFIG = MPRelaxSet.CONFIG - - @property - def incar_updates(self) -> dict: - """Updates to the INCAR config for this calculation type.""" - updates = { - "TEBEG": self.start_temp, - "TEEND": self.end_temp, - "NSW": self.nsteps, - "EDIFF_PER_ATOM": 0.00001, - "LSCALU": False, - "LCHARG": False, - "LPLANE": False, - "LWAVE": True, - "ISMEAR": 0, - "NELMIN": 4, - "LREAL": True, - "BMIX": 1, - "MAXMIX": 20, - "NELM": 500, - "NSIM": 4, - "ISYM": 0, - "ISIF": 0, - "IBRION": 0, - "NBLOCK": 1, - "KBLOCK": 100, - "SMASS": 0, - "PREC": "Normal", - "ISPIN": 2 if self.spin_polarized else 1, - "LDAU": False, - "ADDGRID": True, - "ENCUT": None, - } - if not self.spin_polarized: - updates["MAGMOM"] = None - - if self.time_step is None: - if Element("H") in self.structure.species: # type: ignore - updates.update({"POTIM": 0.5, "NSW": self.nsteps * 4}) - else: - updates["POTIM"] = 2.0 - else: - updates["POTIM"] = self.time_step - - return updates - - @property - def kpoints_updates(self) -> dict | Kpoints: - """Updates to the kpoints configuration for this calculation type.""" - return Kpoints.gamma_automatic() - - -@dataclass -class MVLNPTMDSet(VaspInputSet): - """Write a VASP MD run in NPT ensemble. - - Notes: - To eliminate Pulay stress, the default ENCUT is set to a rather large - value of ENCUT, which is 1.5 * ENMAX. - """ - - start_temp: float = 0.0 - end_temp: float = 300.0 - nsteps: int = 1000 - time_step: float = 2 - spin_polarized: bool = False - CONFIG = MITRelaxSet.CONFIG - - @property - def incar_updates(self) -> dict: - """Updates to the INCAR config for this calculation type.""" - # NPT-AIMD default settings - updates = { - "ALGO": "Fast", - "ISIF": 3, - "LANGEVIN_GAMMA": [10] * self.structure.ntypesp, # type: ignore - "LANGEVIN_GAMMA_L": 1, - "MDALGO": 3, - "PMASS": 10, - "PSTRESS": 0, - "SMASS": 0, - "TEBEG": self.start_temp, - "TEEND": self.end_temp, - "NSW": self.nsteps, - "EDIFF_PER_ATOM": 0.000001, - "LSCALU": False, - "LCHARG": False, - "LPLANE": False, - "LWAVE": True, - "ISMEAR": 0, - "NELMIN": 4, - "LREAL": True, - "BMIX": 1, - "MAXMIX": 20, - "NELM": 500, - "NSIM": 4, - "ISYM": 0, - "IBRION": 0, - "NBLOCK": 1, - "KBLOCK": 100, - "POTIM": self.time_step, - "PREC": "Low", - "ISPIN": 2 if self.spin_polarized else 1, - "LDAU": False, - } - # Set NPT-AIMD ENCUT = 1.5 * VASP_default - enmax = [self.potcar[i].keywords["ENMAX"] for i in range(self.structure.ntypesp)] # type: ignore - updates["ENCUT"] = max(enmax) * 1.5 - return updates - - @property - def kpoints_updates(self) -> Kpoints | dict: - """Updates to the kpoints configuration for this calculation type.""" - return Kpoints.gamma_automatic() - - -@dataclass -class MVLScanRelaxSet(VaspInputSet): - """Write a relax input set using Strongly Constrained and - Appropriately Normed (SCAN) semilocal density functional. - - Notes: - 1. This functional is only available from VASP.5.4.3 upwards. - - 2. Meta-GGA calculations require POTCAR files that include - information on the kinetic energy density of the core-electrons, - i.e. "PBE_52" or "PBE_54". Make sure the POTCAR including the - following lines (see VASP wiki for more details): - - $ grep kinetic POTCAR - kinetic energy-density - mkinetic energy-density pseudized - kinetic energy density (partial) - - Args: - structure (Structure): input structure. - vdw (str): set "rVV10" to enable SCAN+rVV10, which is a versatile - van der Waals density functional by combing the SCAN functional - with the rVV10 non-local correlation functional. - **kwargs: Other kwargs supported by VaspInputSet. - """ - - user_potcar_functional: UserPotcarFunctional = "PBE_52" - _valid_potcars: Sequence[str] | None = ("PBE_52", "PBE_54") - CONFIG = MPRelaxSet.CONFIG - - def __post_init__(self): - super().__post_init__() - if self.user_potcar_functional not in ("PBE_52", "PBE_54"): - raise ValueError("SCAN calculations require PBE_52 or PBE_54!") - - @property - def incar_updates(self) -> dict: - """Updates to the INCAR config for this calculation type.""" - updates = { - "ADDGRID": True, - "EDIFF": 1e-5, - "EDIFFG": -0.05, - "LASPH": True, - "LDAU": False, - "METAGGA": "SCAN", - "NELM": 200, - } - if self.vdw and self.vdw.lower() == "rvv10": - updates["BPARAM"] = 15.7 # This is the correct BPARAM for SCAN+rVV10 - return updates - - -@dataclass -class LobsterSet(VaspInputSet): - """Input set to prepare VASP runs that can be digested by Lobster (See cohp.de). - - Args: - structure (Structure): input structure. - isym (int): ISYM entry for INCAR, only isym=-1 and isym=0 are allowed - ismear (int): ISMEAR entry for INCAR, only ismear=-5 and ismear=0 are allowed - reciprocal_density (int): density of k-mesh by reciprocal volume - user_supplied_basis (dict): dict including basis functions for all elements in - structure, e.g. {"Fe": "3d 3p 4s", "O": "2s 2p"}; if not supplied, a - standard basis is used - address_basis_file (str): address to a file similar to - "BASIS_PBE_54_standard.yaml" in pymatgen.io.lobster.lobster_basis - user_potcar_settings (dict): dict including potcar settings for all elements in - structure, e.g. {"Fe": "Fe_pv", "O": "O"}; if not supplied, a standard basis - is used. - **kwargs: Other kwargs supported by VaspInputSet. - """ - - isym: int = 0 - ismear: int = -5 - reciprocal_density: int | None = None - address_basis_file: str | None = None - user_supplied_basis: dict | None = None - - # newest potcars are preferred - # Choose PBE_54 unless the user specifies a different potcar_functional - user_potcar_functional: UserPotcarFunctional = "PBE_54" - - CONFIG = MPRelaxSet.CONFIG - _valid_potcars: Sequence[str] | None = ("PBE_52", "PBE_54") - - def __post_init__(self): - super().__post_init__() - warnings.warn("Make sure that all parameters are okay! This is a brand new implementation.") - - if self.isym not in (-1, 0): - raise ValueError("Lobster cannot digest WAVEFUNCTIONS with symmetry. isym must be -1 or 0") - if self.ismear not in (-5, 0): - raise ValueError("Lobster usually works with ismear=-5 or ismear=0") - - self._config_dict["POTCAR"]["W"] = "W_sv" - - @property - def kpoints_updates(self) -> dict | Kpoints: - """Updates to the kpoints configuration for this calculation type.""" - # test, if this is okay - return {"reciprocal_density": self.reciprocal_density or 310} - - @property - def incar_updates(self) -> dict: - """Updates to the INCAR config for this calculation type.""" - from pymatgen.io.lobster import Lobsterin - - potcar_symbols = self.potcar_symbols - - # predefined basis! Check if the basis is okay! (charge spilling and bandoverlaps!) - if self.user_supplied_basis is None and self.address_basis_file is None: - basis = Lobsterin.get_basis(structure=self.structure, potcar_symbols=potcar_symbols) # type: ignore - elif self.address_basis_file is not None: - basis = Lobsterin.get_basis( - structure=self.structure, # type: ignore - potcar_symbols=potcar_symbols, - address_basis_file=self.address_basis_file, - ) - elif self.user_supplied_basis is not None: - # test if all elements from structure are in user_supplied_basis - for atom_type in self.structure.symbol_set: # type: ignore - if atom_type not in self.user_supplied_basis: - raise ValueError(f"There are no basis functions for the atom type {atom_type}") - basis = [f"{key} {value}" for key, value in self.user_supplied_basis.items()] - else: - basis = None - - lobsterin = Lobsterin(settingsdict={"basisfunctions": basis}) - nbands = lobsterin._get_nbands(structure=self.structure) # type: ignore - - return { - "EDIFF": 1e-6, - "NSW": 0, - "LWAVE": True, - "ISYM": self.isym, - "NBANDS": nbands, - "IBRION": -1, - "ISMEAR": self.ismear, - "LORBIT": 11, - "ICHARG": 0, - "ALGO": "Normal", - } - - def get_vasprun_outcar(path: str | Path, parse_dos: bool = True, parse_eigen: bool = True) -> tuple[Vasprun, Outcar]: """Get a Vasprun and Outcar from a directory. @@ -2736,10 +1209,7 @@ def get_vasprun_outcar(path: str | Path, parse_dos: bool = True, parse_eigen: bo outcarfile_fullpath = str(path / "OUTCAR.gz") vsfile = vsfile_fullpath if vsfile_fullpath in vruns else max(vruns) outcarfile = outcarfile_fullpath if outcarfile_fullpath in outcars else max(outcars) - return ( - Vasprun(vsfile, parse_dos=parse_dos, parse_eigen=parse_eigen), - Outcar(outcarfile), - ) + return Vasprun(vsfile, parse_dos=parse_dos, parse_eigen=parse_eigen), Outcar(outcarfile) def get_structure_from_prev_run(vasprun, outcar=None) -> Structure: @@ -2819,7 +1289,7 @@ class BadInputSetWarning(UserWarning): def batch_write_input( structures, - vasp_input_set=MPRelaxSet, + vasp_input_set=None, output_dir=".", make_dir_if_not_present=True, subfolder=None, @@ -2860,6 +1330,11 @@ def batch_write_input( **kwargs: Additional kwargs are passed to the vasp_input_set class in addition to structure. """ + if vasp_input_set is None: + from pymatgen.io.vasp.sets.mp import MPRelaxSet + + vasp_input_set = MPRelaxSet + output_dir = Path(output_dir) for idx, site in enumerate(structures): formula = re.sub(r"\s+", "", site.formula) @@ -2940,109 +1415,6 @@ def get_valid_magmom_struct(structure: Structure, inplace: bool = True, spin_mod return ret_struct -@dataclass -class MPAbsorptionSet(VaspInputSet): - """ - 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. - - Args: - structure (Structure): Input structure. - mode (str): Supported modes are "IPA", "RPA" - copy_wavecar (bool): Whether to copy the WAVECAR from a previous run. - Defaults to True. - 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 VaspInputSet. Typically, user_incar_settings is a - commonly used option. - """ - - # CONFIG = _load_yaml_config("MPAbsorptionSet") - - mode: str = "IPA" - copy_wavecar: bool = True - nbands_factor: float = 2 - reciprocal_density: float = 400 - nkred: tuple[int, int, int] | None = None - nedos: int = 2001 - inherit_incar: bool = True - force_gamma: bool = True - CONFIG = MPRelaxSet.CONFIG - nbands: int | None = None - SUPPORTED_MODES = ("IPA", "RPA") - - def __post_init__(self): - """Validate settings""" - super().__post_init__() - self.mode = self.mode.upper() - if self.mode not in MPAbsorptionSet.SUPPORTED_MODES: - raise ValueError(f"{self.mode} not one of the support modes : {MPAbsorptionSet.SUPPORTED_MODES}") - - @property - def kpoints_updates(self) -> dict | Kpoints: - """Updates to the kpoints configuration for this calculation type. - - 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 {"reciprocal_density": self.reciprocal_density} - - @property - def incar_updates(self) -> dict: - """Updates to the INCAR config for this calculation type.""" - updates = { - "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": self.nedos, - } - - if 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 - updates.update({"ALGO": "CHI", "NELM": 1, "NOMEGA": 1000, "EDIFF": None, "LOPTICS": None, "LWAVE": None}) - - if self.prev_vasprun is not None and self.mode == "IPA": - prev_nbands = int(self.prev_vasprun.parameters["NBANDS"]) if self.nbands is None else self.nbands - updates["NBANDS"] = int(np.ceil(prev_nbands * self.nbands_factor)) - - if self.prev_vasprun is not None and self.mode == "RPA": - # Since in the optical calculation, only the q->0 transition is of interest, - # we can reduce the number of q by the factor of the number of kpoints in - # each corresponding 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 lattices, but using NKREDX, NKREDY, NKREDZ are more sensible for - # other lattices. - self.nkred = self.prev_vasprun.kpoints.kpts[0] if self.nkred is None else self.nkred - updates.update({"NKREDX": self.nkred[0], "NKREDY": self.nkred[1], "NKREDZ": self.nkred[2]}) - - return updates - - def _get_ispin(vasprun: Vasprun | None, outcar: Outcar | None) -> int: """Get value of ISPIN depending on the magnetisation in the OUTCAR and vasprun.""" if outcar is not None and outcar.magnetization is not None: diff --git a/pymatgen/io/vasp/sets/lobster.py b/pymatgen/io/vasp/sets/lobster.py new file mode 100644 index 00000000000..0c2d77b9250 --- /dev/null +++ b/pymatgen/io/vasp/sets/lobster.py @@ -0,0 +1,105 @@ +from __future__ import annotations + +import warnings +from dataclasses import dataclass +from typing import TYPE_CHECKING + +from pymatgen.io.vasp.sets.base import UserPotcarFunctional, VaspInputSet +from pymatgen.io.vasp.sets.mp import MPRelaxSet + +if TYPE_CHECKING: + from collections.abc import Sequence + + from pymatgen.io.vasp.inputs import Kpoints + + +@dataclass +class LobsterSet(VaspInputSet): + """Input set to prepare VASP runs that can be digested by Lobster (See cohp.de). + + Args: + structure (Structure): input structure. + isym (int): ISYM entry for INCAR, only isym=-1 and isym=0 are allowed + ismear (int): ISMEAR entry for INCAR, only ismear=-5 and ismear=0 are allowed + reciprocal_density (int): density of k-mesh by reciprocal volume + user_supplied_basis (dict): dict including basis functions for all elements in + structure, e.g. {"Fe": "3d 3p 4s", "O": "2s 2p"}; if not supplied, a + standard basis is used + address_basis_file (str): address to a file similar to + "BASIS_PBE_54_standard.yaml" in pymatgen.io.lobster.lobster_basis + user_potcar_settings (dict): dict including potcar settings for all elements in + structure, e.g. {"Fe": "Fe_pv", "O": "O"}; if not supplied, a standard basis + is used. + **kwargs: Other kwargs supported by VaspInputSet. + """ + + isym: int = 0 + ismear: int = -5 + reciprocal_density: int | None = None + address_basis_file: str | None = None + user_supplied_basis: dict | None = None + + # newest potcars are preferred + # Choose PBE_54 unless the user specifies a different potcar_functional + user_potcar_functional: UserPotcarFunctional = "PBE_54" + + CONFIG = MPRelaxSet.CONFIG + _valid_potcars: Sequence[str] | None = ("PBE_52", "PBE_54") + + def __post_init__(self): + super().__post_init__() + warnings.warn("Make sure that all parameters are okay! This is a brand new implementation.") + + if self.isym not in (-1, 0): + raise ValueError("Lobster cannot digest WAVEFUNCTIONS with symmetry. isym must be -1 or 0") + if self.ismear not in (-5, 0): + raise ValueError("Lobster usually works with ismear=-5 or ismear=0") + + self._config_dict["POTCAR"]["W"] = "W_sv" + + @property + def kpoints_updates(self) -> dict | Kpoints: + """Updates to the kpoints configuration for this calculation type.""" + # test, if this is okay + return {"reciprocal_density": self.reciprocal_density or 310} + + @property + def incar_updates(self) -> dict: + """Updates to the INCAR config for this calculation type.""" + from pymatgen.io.lobster import Lobsterin + + potcar_symbols = self.potcar_symbols + + # predefined basis! Check if the basis is okay! (charge spilling and bandoverlaps!) + if self.user_supplied_basis is None and self.address_basis_file is None: + basis = Lobsterin.get_basis(structure=self.structure, potcar_symbols=potcar_symbols) # type: ignore + elif self.address_basis_file is not None: + basis = Lobsterin.get_basis( + structure=self.structure, # type: ignore + potcar_symbols=potcar_symbols, + address_basis_file=self.address_basis_file, + ) + elif self.user_supplied_basis is not None: + # test if all elements from structure are in user_supplied_basis + for atom_type in self.structure.symbol_set: # type: ignore + if atom_type not in self.user_supplied_basis: + raise ValueError(f"There are no basis functions for the atom type {atom_type}") + basis = [f"{key} {value}" for key, value in self.user_supplied_basis.items()] + else: + basis = None + + lobsterin = Lobsterin(settingsdict={"basisfunctions": basis}) + nbands = lobsterin._get_nbands(structure=self.structure) # type: ignore + + return { + "EDIFF": 1e-6, + "NSW": 0, + "LWAVE": True, + "ISYM": self.isym, + "NBANDS": nbands, + "IBRION": -1, + "ISMEAR": self.ismear, + "LORBIT": 11, + "ICHARG": 0, + "ALGO": "Normal", + } diff --git a/pymatgen/io/vasp/sets/matpes.py b/pymatgen/io/vasp/sets/matpes.py new file mode 100644 index 00000000000..8f1bccb58df --- /dev/null +++ b/pymatgen/io/vasp/sets/matpes.py @@ -0,0 +1,84 @@ +from __future__ import annotations + +import warnings +from dataclasses import dataclass +from typing import TYPE_CHECKING + +from pymatgen.io.vasp.sets.base import VaspInputSet, _load_yaml_config + +if TYPE_CHECKING: + from typing import Literal + + +@dataclass +class MatPESStaticSet(VaspInputSet): + """Create input files for a MatPES static calculation. + + The goal of MatPES is to generate potential energy surface data. This is a distinctly different + from the objectives of the MP static calculations, which aims to obtain primarily accurate + energies and also electronic structure (DOS). For PES data, force accuracy (and to some extent, + stress accuracy) is of paramount importance. + + The default POTCAR versions have been updated to PBE_54 from the old PBE set used in the + MPStaticSet. However, **U values** are still based on PBE. The implicit assumption here is that + the PBE_54 and PBE POTCARs are sufficiently similar that the U values fitted to the old PBE + functional still applies. + + Args: + structure (Structure): The Structure to create inputs for. If None, the input + set is initialized without a Structure but one must be set separately before + the inputs are generated. + xc_functional ('R2SCAN'|'PBE'): Exchange-correlation functional to use. Defaults to 'PBE'. + **kwargs: Keywords supported by VaspInputSet. + """ + + xc_functional: Literal["R2SCAN", "PBE", "PBE+U"] = "PBE" + prev_incar: dict | str | None = None + # These are parameters that we will inherit from any previous INCAR supplied. They are mostly parameters related + # to symmetry and convergence set by Custodian when errors are encountered in a previous run. Given that our goal + # is to have a strictly homogeneous PES data, all other parameters (e.g., ISMEAR, ALGO, etc.) are not inherited. + inherit_incar: list[str] | bool = ( # type: ignore + "LPEAD", + "NGX", + "NGY", + "NGZ", + "SYMPREC", + "IMIX", + "LMAXMIX", + "KGAMMA", + "ISYM", + "NCORE", + "NPAR", + "NELMIN", + "IOPT", + "NBANDS", + "KPAR", + "AMIN", + "NELMDL", + "BMIX", + "AMIX_MAG", + "BMIX_MAG", + ) + CONFIG = _load_yaml_config("MatPESStaticSet") + + def __post_init__(self): + """Validate inputs""" + super().__post_init__() + valid_xc_functionals = ("R2SCAN", "PBE", "PBE+U") + if self.xc_functional.upper() not in valid_xc_functionals: + raise ValueError( + f"Unrecognized xc_functional='{self.xc_functional}'. " + f"Supported exchange-correlation functionals are {valid_xc_functionals}" + ) + + default_potcars = self.CONFIG["PARENT"].replace("PBE", "PBE_").replace("Base", "") # PBE64Base -> PBE_64 + self.user_potcar_functional = self.user_potcar_functional or default_potcars + if self.user_potcar_functional.upper() != default_potcars: + warnings.warn( + f"{self.user_potcar_functional=} is inconsistent with the recommended {default_potcars}.", UserWarning + ) + + if self.xc_functional.upper() == "R2SCAN": + self._config_dict["INCAR"].update({"METAGGA": "R2SCAN", "ALGO": "ALL", "GGA": None}) + if self.xc_functional.upper().endswith("+U"): + self._config_dict["INCAR"]["LDAU"] = True diff --git a/pymatgen/io/vasp/sets/mit.py b/pymatgen/io/vasp/sets/mit.py new file mode 100644 index 00000000000..bd0426f351b --- /dev/null +++ b/pymatgen/io/vasp/sets/mit.py @@ -0,0 +1,211 @@ +from __future__ import annotations + +from dataclasses import dataclass +from itertools import chain +from pathlib import Path + +import numpy as np + +from pymatgen.core import Structure +from pymatgen.core.sites import PeriodicSite +from pymatgen.io.vasp.inputs import Kpoints, Poscar +from pymatgen.io.vasp.sets.base import VaspInputSet, _load_yaml_config +from pymatgen.util.due import Doi, due + + +@due.dcite( + Doi("10.1016/j.commatsci.2011.02.023"), + description="A high-throughput infrastructure for density functional theory calculations", +) +@dataclass +class MITRelaxSet(VaspInputSet): + """ + Standard implementation of VaspInputSet utilizing parameters in the MIT + High-throughput project. + The parameters are chosen specifically for a high-throughput project, + which means in general pseudopotentials with fewer electrons were chosen. + + Args: + structure (Structure): The Structure to create inputs for. If None, the input + set is initialized without a Structure but one must be set separately before + the inputs are generated. + **kwargs: Keywords supported by VaspInputSet. + + Please refer:: + + A Jain, G. Hautier, C. Moore, S. P. Ong, C. Fischer, T. Mueller, + K. A. Persson, G. Ceder. A high-throughput infrastructure for density + functional theory calculations. Computational Materials Science, + 2011, 50(8), 2295-2310. doi:10.1016/j.commatsci.2011.02.023 + """ + + CONFIG = _load_yaml_config("MITRelaxSet") + + +class MITNEBSet(VaspInputSet): + """Write NEB inputs. + + Note that EDIFF is not on a per atom basis for this input set. + """ + + def __init__(self, structures, unset_encut=False, **kwargs) -> None: + """ + Args: + structures: List of Structure objects. + unset_encut (bool): Whether to unset ENCUT. + **kwargs: Other kwargs supported by VaspInputSet. + """ + if len(structures) < 3: + raise ValueError(f"You need at least 3 structures for an NEB, got {len(structures)}") + kwargs["sort_structure"] = False + super().__init__(structures[0], MITRelaxSet.CONFIG, **kwargs) + self.structures = self._process_structures(structures) + + self.unset_encut = False + if unset_encut: + self._config_dict["INCAR"].pop("ENCUT", None) + + if "EDIFF" not in self._config_dict["INCAR"]: + self._config_dict["INCAR"]["EDIFF"] = self._config_dict["INCAR"].pop("EDIFF_PER_ATOM") + + # NEB specific defaults + defaults = {"IMAGES": len(structures) - 2, "IBRION": 1, "ISYM": 0, "LCHARG": False, "LDAU": False} + self._config_dict["INCAR"].update(defaults) + + @property + def poscar(self): + """Poscar for structure of first end point.""" + return Poscar(self.structures[0]) + + @property + def poscars(self): + """List of Poscars.""" + return [Poscar(s) for s in self.structures] + + @staticmethod + def _process_structures(structures): + """Remove any atom jumps across the cell.""" + input_structures = structures + structures = [input_structures[0]] + for s in input_structures[1:]: + prev = structures[-1] + for idx, site in enumerate(s): + translate = np.round(prev[idx].frac_coords - site.frac_coords) + if np.any(np.abs(translate) > 0.5): + s.translate_sites([idx], translate, to_unit_cell=False) + structures.append(s) + return structures + + def write_input( + self, + output_dir, + make_dir_if_not_present=True, + write_cif=False, + write_path_cif=False, + write_endpoint_inputs=False, + ): + """ + NEB inputs has a special directory structure where inputs are in 00, + 01, 02, .... + + Args: + output_dir (str): Directory to output the VASP input files + make_dir_if_not_present (bool): Set to True if you want the + directory (and the whole path) to be created if it is not + present. + write_cif (bool): If true, writes a cif along with each POSCAR. + write_path_cif (bool): If true, writes a cif for each image. + write_endpoint_inputs (bool): If true, writes input files for + running endpoint calculations. + """ + output_dir = Path(output_dir) + if make_dir_if_not_present and not output_dir.exists(): + output_dir.mkdir(parents=True) + self.incar.write_file(str(output_dir / "INCAR")) + self.kpoints.write_file(str(output_dir / "KPOINTS")) + self.potcar.write_file(str(output_dir / "POTCAR")) + + for idx, poscar in enumerate(self.poscars): + d = output_dir / str(idx).zfill(2) + if not d.exists(): + d.mkdir(parents=True) + poscar.write_file(str(d / "POSCAR")) + if write_cif: + poscar.structure.to(filename=str(d / f"{idx}.cif")) + if write_endpoint_inputs: + end_point_param = MITRelaxSet(self.structures[0], user_incar_settings=self.user_incar_settings) + + for image in ["00", str(len(self.structures) - 1).zfill(2)]: + end_point_param.incar.write_file(str(output_dir / image / "INCAR")) + end_point_param.kpoints.write_file(str(output_dir / image / "KPOINTS")) + end_point_param.potcar.write_file(str(output_dir / image / "POTCAR")) + if write_path_cif: + sites = { + PeriodicSite(site.species, site.frac_coords, self.structures[0].lattice) + for site in chain(*(struct for struct in self.structures)) + } + neb_path = Structure.from_sites(sorted(sites)) + neb_path.to(filename=f"{output_dir}/path.cif") + + +@dataclass +class MITMDSet(VaspInputSet): + """Write a VASP MD run. This DOES NOT do multiple stage runs. + + Args: + structure (Structure): Input structure. + start_temp (float): Starting temperature. + end_temp (float): Final temperature. + nsteps (int): Number of time steps for simulations. NSW parameter. + time_step (float): The time step for the simulation. The POTIM + parameter. Defaults to 2fs. + spin_polarized (bool): Whether to do spin polarized calculations. + The ISPIN parameter. Defaults to False. + **kwargs: Other kwargs supported by VaspInputSet. + """ + + structure: Structure | None = None + start_temp: float = 0.0 + end_temp: float = 300.0 + nsteps: int = 1000 + time_step: float = 2 + spin_polarized: bool = False + CONFIG = MITRelaxSet.CONFIG + + @property + def incar_updates(self): + """Updates to the INCAR config for this calculation type.""" + # MD default settings + return { + "TEBEG": self.start_temp, + "TEEND": self.end_temp, + "NSW": self.nsteps, + "EDIFF_PER_ATOM": 0.000001, + "LSCALU": False, + "LCHARG": False, + "LPLANE": False, + "LWAVE": True, + "ISMEAR": 0, + "NELMIN": 4, + "LREAL": True, + "BMIX": 1, + "MAXMIX": 20, + "NELM": 500, + "NSIM": 4, + "ISYM": 0, + "ISIF": 0, + "IBRION": 0, + "NBLOCK": 1, + "KBLOCK": 100, + "SMASS": 0, + "POTIM": self.time_step, + "PREC": "Low", + "ISPIN": 2 if self.spin_polarized else 1, + "LDAU": False, + "ENCUT": None, + } + + @property + def kpoints_updates(self) -> Kpoints | dict: + """Updates to the kpoints configuration for this calculation type.""" + return Kpoints.gamma_automatic() diff --git a/pymatgen/io/vasp/sets/mp.py b/pymatgen/io/vasp/sets/mp.py new file mode 100644 index 00000000000..331b99a14a2 --- /dev/null +++ b/pymatgen/io/vasp/sets/mp.py @@ -0,0 +1,853 @@ +from __future__ import annotations + +import warnings +from copy import deepcopy +from dataclasses import dataclass, field +from typing import TYPE_CHECKING + +import numpy as np +from monty.serialization import loadfn + +from pymatgen.core import MODULE_DIR, Element, Species, Structure +from pymatgen.io.vasp.inputs import Kpoints +from pymatgen.io.vasp.sets.base import UserPotcarFunctional, VaspInputSet, _get_ispin, _get_nedos, _load_yaml_config +from pymatgen.util.due import Doi, due + +if TYPE_CHECKING: + from collections.abc import Sequence + from typing import Any, Literal + + from pymatgen.util.typing import Vector3D + + +@dataclass +class MPRelaxSet(VaspInputSet): + """ + Implementation of VaspInputSet utilizing parameters in the public + Materials Project. Typically, the pseudopotentials chosen contain more + electrons than the MIT parameters, and the k-point grid is ~50% more dense. + The LDAUU parameters are also different due to the different PSPs used, + which result in different fitted values. + + Args: + structure (Structure): The Structure to create inputs for. If None, the input + set is initialized without a Structure but one must be set separately before + the inputs are generated. + **kwargs: Keywords supported by VaspInputSet. + """ + + CONFIG = _load_yaml_config("MPRelaxSet") + + +@due.dcite( + Doi("10.1021/acs.jpclett.0c02405"), + description="AccurAccurate and Numerically Efficient r2SCAN Meta-Generalized Gradient Approximation", +) +@due.dcite( + Doi("10.1103/PhysRevLett.115.036402"), + description="Strongly Constrained and Appropriately Normed Semilocal Density Functional", +) +@due.dcite( + Doi("10.1103/PhysRevB.93.155109"), + description="Efficient generation of generalized Monkhorst-Pack grids through the use of informatics", +) +@dataclass +class MPScanRelaxSet(VaspInputSet): + """Write a relaxation input set using the accurate and numerically + efficient r2SCAN variant of the Strongly Constrained and Appropriately Normed + (SCAN) metaGGA density functional. + + Notes: + 1. This functional is officially supported in VASP 6.0.0 and above. On older version, + source code may be obtained by contacting the authors of the referenced manuscript. + The original SCAN functional, available from VASP 5.4.3 onwards, maybe used instead + by passing `user_incar_settings={"METAGGA": "SCAN"}` when instantiating this InputSet. + r2SCAN and SCAN are expected to yield very similar results. + + 2. Meta-GGA calculations require POTCAR files that include + information on the kinetic energy density of the core-electrons, + i.e. "PBE_52" or "PBE_54". Make sure the POTCARs include the + following lines (see VASP wiki for more details): + + $ grep kinetic POTCAR + kinetic energy-density + mkinetic energy-density pseudized + kinetic energy density (partial) + + Args: + bandgap (float): Bandgap of the structure in eV. The bandgap is used to + compute the appropriate k-point density and determine the + smearing settings. + Metallic systems (default, bandgap = 0) use a KSPACING value of 0.22 + and Methfessel-Paxton order 2 smearing (ISMEAR=2, SIGMA=0.2). + Non-metallic systems (bandgap > 0) use the tetrahedron smearing + method (ISMEAR=-5, SIGMA=0.05). The KSPACING value is + calculated from the bandgap via Eqs. 25 and 29 of Wisesa, McGill, + and Mueller [1] (see References). Note that if 'user_incar_settings' + or 'user_kpoints_settings' override KSPACING, the calculation from + bandgap is not performed. + vdw (str): set "rVV10" to enable SCAN+rVV10, which is a versatile + van der Waals density functional by combing the SCAN functional + with the rVV10 non-local correlation functional. rvv10 is the only + dispersion correction available for SCAN at this time. + **kwargs: Keywords supported by VaspInputSet. + + References: + [1] P. Wisesa, K.A. McGill, T. Mueller, Efficient generation of + generalized Monkhorst-Pack grids through the use of informatics, + Phys. Rev. B. 93 (2016) 1-10. doi:10.1103/PhysRevB.93.155109. + + References: + James W. Furness, Aaron D. Kaplan, Jinliang Ning, John P. Perdew, and Jianwei Sun. + Accurate and Numerically Efficient r2SCAN Meta-Generalized Gradient Approximation. + The Journal of Physical Chemistry Letters 0, 11 DOI: 10.1021/acs.jpclett.0c02405 + """ + + bandgap: float | None = None + auto_kspacing: bool = True + user_potcar_functional: UserPotcarFunctional = "PBE_54" + auto_ismear: bool = True + CONFIG = _load_yaml_config("MPSCANRelaxSet") + _valid_potcars: Sequence[str] | None = ("PBE_52", "PBE_54") + + def __post_init__(self): + super().__post_init__() + if self.vdw and self.vdw != "rvv10": + warnings.warn("Use of van der waals functionals other than rVV10 with SCAN is not supported at this time. ") + # delete any vdw parameters that may have been added to the INCAR + vdw_par = loadfn(str(MODULE_DIR / "vdW_parameters.yaml")) + for k in vdw_par[self.vdw]: + self._config_dict["INCAR"].pop(k, None) + + +@dataclass +class MPMetalRelaxSet(VaspInputSet): + """ + Implementation of VaspInputSet utilizing parameters in the public + Materials Project, but with tuning for metals. Key things are a denser + k point density, and a. + """ + + CONFIG = MPRelaxSet.CONFIG + + @property + def incar_updates(self) -> dict: + """Updates to the INCAR config for this calculation type.""" + return {"ISMEAR": 1, "SIGMA": 0.2} + + @property + def kpoints_updates(self) -> dict | Kpoints: + """Updates to the kpoints configuration for this calculation type.""" + return {"reciprocal_density": 200} + + +@dataclass +class MPHSERelaxSet(VaspInputSet): + """Same as the MPRelaxSet, but with HSE parameters.""" + + CONFIG = _load_yaml_config("MPHSERelaxSet") + + +@dataclass +class MPStaticSet(VaspInputSet): + """Create input files for a static calculation. + + Args: + structure (Structure): Structure from previous run. + lepsilon (bool): Whether to add static dielectric calculation + lcalcpol (bool): Whether to turn on evaluation of the Berry phase approximations + for electronic polarization + reciprocal_density (int): For static calculations, we usually set the + reciprocal density by volume. This is a convenience arg to change + that, rather than using user_kpoints_settings. Defaults to 100, + which is ~50% more than that of standard relaxation calculations. + small_gap_multiply ([float, float]): If the gap is less than + 1st index, multiply the default reciprocal_density by the 2nd + index. + **kwargs: Keywords supported by MPRelaxSet. + """ + + lepsilon: bool = False + lcalcpol: bool = False + reciprocal_density: int = 100 + small_gap_multiply: tuple[float, float] | None = None + inherit_incar: bool = True + CONFIG = MPRelaxSet.CONFIG + + @property + def incar_updates(self) -> dict: + """Updates to the INCAR config for this calculation type.""" + updates: dict[str, Any] = {"NSW": 0, "ISMEAR": -5, "LCHARG": True, "LORBIT": 11, "LREAL": False} + if self.lepsilon: + # LPEAD=T: numerical evaluation of overlap integral prevents LRF_COMMUTATOR + # errors and can lead to better expt. agreement but produces slightly + # different results + updates.update({"IBRION": 8, "LEPSILON": True, "LPEAD": True, "NSW": 1, "EDIFF": 1e-5}) + + if self.lcalcpol: + updates["LCALCPOL"] = True + return updates + + @property + def kpoints_updates(self) -> dict | Kpoints: + """Updates to the kpoints configuration for this calculation type.""" + factor = 1.0 + if self.bandgap is not None and self.small_gap_multiply and self.bandgap <= self.small_gap_multiply[0]: + factor = self.small_gap_multiply[1] + + # prefer to use k-point scheme from previous run unless lepsilon = True is specified + if self.prev_kpoints and self.prev_kpoints.style == Kpoints.supported_modes.Monkhorst and not self.lepsilon: # type: ignore + kpoints = Kpoints.automatic_density_by_vol( + self.structure, # type: ignore + int(self.reciprocal_density * factor), + self.force_gamma, + ) + k_div = [kp + 1 if kp % 2 == 1 else kp for kp in kpoints.kpts[0]] # type: ignore + return Kpoints.monkhorst_automatic(k_div) # type: ignore + + return {"reciprocal_density": self.reciprocal_density * factor} + + +@dataclass +class MPScanStaticSet(MPScanRelaxSet): + """Create input files for a static calculation using the accurate and numerically + efficient r2SCAN variant of the Strongly Constrained and Appropriately Normed + (SCAN) metaGGA functional. + + Args: + structure (Structure): Structure from previous run. + bandgap (float): Bandgap of the structure in eV. The bandgap is used to + compute the appropriate k-point density and determine the smearing settings. + lepsilon (bool): Whether to add static dielectric calculation + lcalcpol (bool): Whether to turn on evaluation of the Berry phase approximations + for electronic polarization. + **kwargs: Keywords supported by MPScanRelaxSet. + """ + + lepsilon: bool = False + lcalcpol: bool = False + inherit_incar: bool = True + auto_kspacing: bool = True + + @property + def incar_updates(self) -> dict: + """Updates to the INCAR config for this calculation type.""" + updates: dict[str, Any] = { + "LREAL": False, + "NSW": 0, + "LORBIT": 11, + "LVHAR": True, + "ISMEAR": -5, + } + + if self.lepsilon: + # LPEAD=T: numerical evaluation of overlap integral prevents + # LRF_COMMUTATOR errors and can lead to better expt. agreement + # but produces slightly different results + updates.update({"IBRION": 8, "LEPSILON": True, "LPEAD": True, "NSW": 1, "NPAR": None}) + + if self.lcalcpol: + updates["LCALCPOL"] = True + + return updates + + +@dataclass +class MPHSEBSSet(VaspInputSet): + """ + Implementation of a VaspInputSet for HSE band structure computations. + + Remember that HSE band structures must be self-consistent in VASP. A band structure + along symmetry lines for instance needs BOTH a uniform grid with appropriate weights + AND a path along the lines with weight 0. + + Thus, the "uniform" mode is just like regular static SCF but allows adding custom + kpoints (e.g., corresponding to known VBM/CBM) to the uniform grid that have zero + weight (e.g., for better gap estimate). + + The "gap" mode behaves just like the "uniform" mode, however, if starting from a + previous calculation, the VBM and CBM k-points will automatically be added to + ``added_kpoints``. + + The "line" mode is just like "uniform" mode, but additionally adds k-points along + symmetry lines with zero weight. + + The "uniform_dense" mode is like "uniform" mode but additionally adds a denser + uniform mesh with zero weight. This can be useful when calculating Fermi surfaces + or BoltzTraP/AMSET electronic transport using hybrid DFT. + + Args: + structure (Structure): Structure to compute + added_kpoints (list): a list of kpoints (list of 3 number list) added to the + run. The k-points are in fractional coordinates + mode (str): "Line" - generate k-points along symmetry lines for bandstructure. + "Uniform" - generate uniform k-points grid. + reciprocal_density (int): k-point density to use for uniform mesh. + copy_chgcar (bool): Whether to copy the CHGCAR of a previous run. + kpoints_line_density (int): k-point density for high symmetry lines + dedos (float): Energy difference used to set NEDOS, based on the total energy + range. + optics (bool): Whether to add LOPTICS (used for calculating optical response). + nbands_factor (float): Multiplicative factor for NBANDS when starting from a + previous calculation. Choose a higher number if you are doing an LOPTICS + calculation. + **kwargs: Keywords supported by VaspInputSet. + """ + + added_kpoints: list[Vector3D] = field(default_factory=list) + mode: str = "gap" + reciprocal_density: float = 50 + copy_chgcar: bool = True + kpoints_line_density: float = 20 + nbands_factor: float = 1.2 + zero_weighted_reciprocal_density: float = 100 + dedos: float = 0.02 + optics: bool = False + CONFIG = MPHSERelaxSet.CONFIG + + def __post_init__(self) -> None: + """Ensure mode is set correctly.""" + super().__post_init__() + + if "reciprocal_density" in self.user_kpoints_settings: + self.reciprocal_density = self.user_kpoints_settings["reciprocal_density"] + + self.mode = self.mode.lower() + supported_modes = ("line", "uniform", "gap", "uniform_dense") + if self.mode not in supported_modes: + raise ValueError(f"Supported modes are: {', '.join(supported_modes)}") + + @property + def kpoints_updates(self) -> dict: + """Updates to the kpoints configuration for this calculation type.""" + kpoints: dict[str, Any] = {"reciprocal_density": self.reciprocal_density, "explicit": True} + + if self.mode == "line": + # add line_density on top of reciprocal density + kpoints["zero_weighted_line_density"] = self.kpoints_line_density + + elif self.mode == "uniform_dense": + kpoints["zero_weighted_reciprocal_density"] = self.zero_weighted_reciprocal_density + + added_kpoints = deepcopy(self.added_kpoints) + if self.prev_vasprun is not None and self.mode == "gap": + bs = self.prev_vasprun.get_band_structure() + if not bs.is_metal(): + added_kpoints.append(bs.get_vbm()["kpoint"].frac_coords) + added_kpoints.append(bs.get_cbm()["kpoint"].frac_coords) + + kpoints["added_kpoints"] = added_kpoints + + return kpoints + + @property + def incar_updates(self) -> dict: + """Updates to the INCAR config for this calculation type.""" + updates = dict(NSW=0, ISMEAR=0, SIGMA=0.05, ISYM=3, LCHARG=False, NELMIN=5) + + if self.mode == "uniform" and len(self.added_kpoints) == 0: + # automatic setting of nedos using the energy range and the energy step + nedos = _get_nedos(self.prev_vasprun, self.dedos) + + # use tetrahedron method for DOS and optics calculations + updates.update({"ISMEAR": -5, "NEDOS": nedos}) + + else: + # if line mode or explicit k-points (gap) can't use ISMEAR=-5 + # use small sigma to avoid partial occupancies for small band gap materials + updates.update({"ISMEAR": 0, "SIGMA": 0.01}) + + if self.prev_vasprun is not None: + # set nbands + nbands = int(np.ceil(self.prev_vasprun.parameters["NBANDS"] * self.nbands_factor)) + updates["NBANDS"] = nbands + + if self.optics: + # LREAL not supported with LOPTICS + updates.update({"LOPTICS": True, "LREAL": False, "CSHIFT": 1e-5}) + + if self.prev_vasprun is not None and self.prev_outcar is not None: + # turn off spin when magmom for every site is smaller than 0.02. + updates["ISPIN"] = _get_ispin(self.prev_vasprun, self.prev_outcar) + + return updates + + +@dataclass +class MPNonSCFSet(VaspInputSet): + """ + Init a MPNonSCFSet. Typically, you would use the classmethod + from_prev_calc to initialize from a previous SCF run. + + Args: + structure (Structure): Structure to compute + mode (str): Line, Uniform or Boltztrap mode supported. + nedos (int): nedos parameter. Default to 2001. + dedos (float): setting nedos=0 and uniform mode in from_prev_calc, + an automatic nedos will be calculated using the total energy range + divided by the energy step dedos + reciprocal_density (int): density of k-mesh by reciprocal + volume (defaults to 100) + kpoints_line_density (int): Line density for Line mode. + optics (bool): whether to add dielectric function + copy_chgcar: Whether to copy the old CHGCAR when starting from a + previous calculation. + nbands_factor (float): Multiplicative factor for NBANDS when starting + from a previous calculation. Choose a higher number if you are + doing an LOPTICS calculation. + small_gap_multiply ([float, float]): When starting from a previous + calculation, if the gap is less than 1st index, multiply the default + reciprocal_density by the 2nd index. + **kwargs: Keywords supported by MPRelaxSet. + """ + + mode: str = "line" + nedos: int = 2001 + dedos: float = 0.005 + reciprocal_density: float = 100 + kpoints_line_density: float = 20 + optics: bool = False + copy_chgcar: bool = True + nbands_factor: float = 1.2 + small_gap_multiply: tuple[float, float] | None = None + inherit_incar: bool = True + CONFIG = MPRelaxSet.CONFIG + + def __post_init__(self): + """Perform inputset validation.""" + super().__post_init__() + + mode = self.mode = self.mode.lower() + + valid_modes = ("line", "uniform", "boltztrap") + if mode not in valid_modes: + raise ValueError( + f"Invalid {mode=}. Supported modes for NonSCF runs are {', '.join(map(repr, valid_modes))}" + ) + + if (mode != "uniform" or self.nedos < 2000) and self.optics: + warnings.warn("It is recommended to use Uniform mode with a high NEDOS for optics calculations.") + + if self.standardize: + warnings.warn( + "Use of standardize=True with from_prev_run is not " + "recommended as there is no guarantee the copied " + "files will be appropriate for the standardized" + " structure. copy_chgcar is enforced to be false." + ) + self.copy_chgcar = False + + @property + def incar_updates(self) -> dict: + """Updates to the INCAR config for this calculation type.""" + updates: dict[str, Any] = {"LCHARG": False, "LORBIT": 11, "LWAVE": False, "NSW": 0, "ISYM": 0, "ICHARG": 11} + + if self.prev_vasprun is not None: + # set NBANDS + n_bands = int(np.ceil(self.prev_vasprun.parameters["NBANDS"] * self.nbands_factor)) + updates["NBANDS"] = n_bands + + # automatic setting of NEDOS using the energy range and the energy step + nedos = _get_nedos(self.prev_vasprun, self.dedos) if self.nedos == 0 else self.nedos + + if self.mode == "uniform": + # use tetrahedron method for DOS and optics calculations + updates.update({"ISMEAR": -5, "ISYM": 2, "NEDOS": nedos}) + + elif self.mode in ("line", "boltztrap"): + # if line mode or explicit k-points (boltztrap) can't use ISMEAR=-5 + # use small sigma to avoid partial occupancies for small band gap materials + # use a larger sigma if the material is a metal + sigma = 0.2 if self.bandgap == 0 or self.bandgap is None else 0.01 + updates.update({"ISMEAR": 0, "SIGMA": sigma}) + + if self.optics: + # LREAL not supported with LOPTICS = True; automatic NEDOS usually + # underestimates, so set it explicitly + updates.update({"LOPTICS": True, "LREAL": False, "CSHIFT": 1e-5, "NEDOS": nedos}) + + if self.prev_vasprun is not None and self.prev_outcar is not None: + # turn off spin when magmom for every site is smaller than 0.02. + updates["ISPIN"] = _get_ispin(self.prev_vasprun, self.prev_outcar) + + updates["MAGMOM"] = None + return updates + + @property + def kpoints_updates(self) -> dict: + """Updates to the kpoints configuration for this calculation type.""" + factor = 1.0 + if self.bandgap is not None and self.small_gap_multiply and self.bandgap <= self.small_gap_multiply[0]: + factor = self.small_gap_multiply[1] + + if self.mode == "line": + return {"line_density": self.kpoints_line_density * factor} + + if self.mode == "boltztrap": + return {"explicit": True, "reciprocal_density": self.reciprocal_density * factor} + + return {"reciprocal_density": self.reciprocal_density * factor} + + +@dataclass +class MPSOCSet(VaspInputSet): + """An input set for running spin-orbit coupling (SOC) calculations. + + Args: + structure (Structure): the structure must have the 'magmom' site + property and each magnetic moment value must have 3 + components. eg: ``magmom = [[0,0,2], ...]`` + saxis (tuple): magnetic moment orientation + copy_chgcar: Whether to copy the old CHGCAR. Defaults to True. + nbands_factor (float): Multiplicative factor for NBANDS. Choose a + higher number if you are doing an LOPTICS calculation. + reciprocal_density (int): density of k-mesh by reciprocal volume. + small_gap_multiply ([float, float]): If the gap is less than + 1st index, multiply the default reciprocal_density by the 2nd + index. + lepsilon (bool): Whether to add static dielectric calculation + lcalcpol (bool): Whether to turn on evaluation of the Berry phase approximations + for electronic polarization + magmom (list[list[float]]): Override for the structure magmoms. + **kwargs: Keywords supported by VaspInputSet. + """ + + saxis: tuple[int, int, int] = (0, 0, 1) + nbands_factor: float = 1.2 + lepsilon: bool = False + lcalcpol: bool = False + reciprocal_density: float = 100 + small_gap_multiply: tuple[float, float] | None = None + magmom: list[Vector3D] | None = None + inherit_incar: bool = True + copy_chgcar: bool = True + CONFIG = MPRelaxSet.CONFIG + + def __post_init__(self): + super().__post_init__() + if ( + self.structure + and not hasattr(self.structure[0], "magmom") + and not isinstance(self.structure[0].magmom, list) + ): + raise ValueError( + "The structure must have the 'magmom' site property and each magnetic " + "moment value must have 3 components. e.g. magmom = [0,0,2]" + ) + + @property + def incar_updates(self) -> dict: + """Updates to the INCAR config for this calculation type.""" + updates = { + "ISYM": -1, + "LSORBIT": "T", + "ICHARG": 11, + "SAXIS": list(self.saxis), + "NSW": 0, + "ISMEAR": -5, + "LCHARG": True, + "LORBIT": 11, + "LREAL": False, + } + + if self.lepsilon: + # LPEAD=T: numerical evaluation of overlap integral prevents LRF_COMMUTATOR + # errors and can lead to better expt. agreement but produces slightly + # different results + updates.update({"IBRION": 8, "LEPSILON": True, "LPEAD": True, "NSW": 1}) + + if self.lcalcpol: + updates["LCALCPOL"] = True + + if self.prev_vasprun is not None: + # set NBANDS + n_bands = int(np.ceil(self.prev_vasprun.parameters["NBANDS"] * self.nbands_factor)) + updates["NBANDS"] = n_bands + return updates + + @property + def kpoints_updates(self) -> dict: + """Updates to the kpoints configuration for this calculation type.""" + factor = 1.0 + if self.bandgap is not None and self.small_gap_multiply and self.bandgap <= self.small_gap_multiply[0]: + factor = self.small_gap_multiply[1] + return {"reciprocal_density": self.reciprocal_density * factor} + + @VaspInputSet.structure.setter # type: ignore + def structure(self, structure: Structure | None) -> None: + if structure is not None: + if self.magmom: + structure = structure.copy(site_properties={"magmom": self.magmom}) + + # magmom has to be 3D for SOC calculation. + if hasattr(structure[0], "magmom"): + if not isinstance(structure[0].magmom, list): + # project magmom to z-axis + structure = structure.copy(site_properties={"magmom": [[0, 0, site.magmom] for site in structure]}) + else: + raise ValueError("Neither the previous structure has magmom property nor magmom provided") + + VaspInputSet.structure.fset(self, structure) # type: ignore + + +@dataclass +class MPNMRSet(VaspInputSet): + """Init a MPNMRSet. + + Args: + structure (Structure): Structure from previous run. + mode (str): The NMR calculation to run + "cs": for Chemical Shift + "efg" for Electric Field Gradient + isotopes (list): list of Isotopes for quadrupole moments + reciprocal_density (int): density of k-mesh by reciprocal volume. + lepsilon (bool): Whether to add static dielectric calculation + lcalcpol (bool): Whether to turn on evaluation of the Berry phase approximations + for electronic polarization + reciprocal_density (int): For static calculations, we usually set the + reciprocal density by volume. This is a convenience arg to change + that, rather than using user_kpoints_settings. Defaults to 100, + which is ~50% more than that of standard relaxation calculations. + small_gap_multiply ([float, float]): If the gap is less than + 1st index, multiply the default reciprocal_density by the 2nd + index. + **kwargs: Keywords supported by MPRelaxSet. + """ + + mode: Literal["cs", "efg"] = "cs" + isotopes: list = field(default_factory=list) + reciprocal_density: int = 100 + small_gap_multiply: tuple[float, float] | None = None + inherit_incar: bool = True + CONFIG = MPRelaxSet.CONFIG + + @property + def incar_updates(self) -> dict: + """Updates to the INCAR config for this calculation type.""" + updates: dict[str, Any] = {"NSW": 0, "ISMEAR": -5, "LCHARG": True, "LORBIT": 11, "LREAL": False} + if self.mode.lower() == "cs": + updates.update( + LCHIMAG=True, + EDIFF=-1.0e-10, + ISYM=0, + LCHARG=False, + LNMR_SYM_RED=True, + NELMIN=10, + NLSPLINE=True, + PREC="ACCURATE", + SIGMA=0.01, + ) + elif self.mode.lower() == "efg": + isotopes = {ist.split("-")[0]: ist for ist in self.isotopes} + quad_efg = [ + float(Species(s.name).get_nmr_quadrupole_moment(isotopes.get(s.name))) + for s in self.structure.species # type: ignore + ] + updates.update( + ALGO="FAST", + EDIFF=-1.0e-10, + ISYM=0, + LCHARG=False, + LEFG=True, + QUAD_EFG=quad_efg, + NELMIN=10, + PREC="ACCURATE", + SIGMA=0.01, + ) + return updates + + @property + def kpoints_updates(self) -> dict: + """Updates to the kpoints configuration for this calculation type.""" + factor = 1.0 + if self.bandgap is not None and self.small_gap_multiply and self.bandgap <= self.small_gap_multiply[0]: + factor = self.small_gap_multiply[1] + return {"reciprocal_density": self.reciprocal_density * factor} + + +@dataclass +class MPMDSet(VaspInputSet): + """ + This a modified version of the old MITMDSet pre 2018/03/12. + + This set serves as the basis for the amorphous skyline paper. + + (1) Aykol, M.; Dwaraknath, S. S.; Sun, W.; Persson, K. A. Thermodynamic + Limit for Synthesis of Metastable Inorganic Materials. Sci. Adv. 2018, + 4 (4). + + Class for writing a VASP MD run. This DOES NOT do multiple stage runs. + Precision remains normal, to increase accuracy of stress tensor. + + Args: + structure (Structure): Input structure. + start_temp (int): Starting temperature. + end_temp (int): Final temperature. + nsteps (int): Number of time steps for simulations. NSW parameter. + time_step (float): The time step for the simulation. The POTIM + parameter. Defaults to None, which will set it automatically + to 2.0 fs for non-hydrogen containing structures and 0.5 fs + for hydrogen containing structures. + spin_polarized (bool): Whether to do spin polarized calculations. + The ISPIN parameter. Defaults to False. + **kwargs: Other kwargs supported by VaspInputSet. + """ + + start_temp: float = 0.0 + end_temp: float = 300.0 + nsteps: int = 1000 + time_step: float | None = None + spin_polarized: bool = False + CONFIG = MPRelaxSet.CONFIG + + @property + def incar_updates(self) -> dict: + """Updates to the INCAR config for this calculation type.""" + updates = { + "TEBEG": self.start_temp, + "TEEND": self.end_temp, + "NSW": self.nsteps, + "EDIFF_PER_ATOM": 0.00001, + "LSCALU": False, + "LCHARG": False, + "LPLANE": False, + "LWAVE": True, + "ISMEAR": 0, + "NELMIN": 4, + "LREAL": True, + "BMIX": 1, + "MAXMIX": 20, + "NELM": 500, + "NSIM": 4, + "ISYM": 0, + "ISIF": 0, + "IBRION": 0, + "NBLOCK": 1, + "KBLOCK": 100, + "SMASS": 0, + "PREC": "Normal", + "ISPIN": 2 if self.spin_polarized else 1, + "LDAU": False, + "ADDGRID": True, + "ENCUT": None, + } + if not self.spin_polarized: + updates["MAGMOM"] = None + + if self.time_step is None: + if Element("H") in self.structure.species: # type: ignore + updates.update({"POTIM": 0.5, "NSW": self.nsteps * 4}) + else: + updates["POTIM"] = 2.0 + else: + updates["POTIM"] = self.time_step + + return updates + + @property + def kpoints_updates(self) -> dict | Kpoints: + """Updates to the kpoints configuration for this calculation type.""" + return Kpoints.gamma_automatic() + + +@dataclass +class MPAbsorptionSet(VaspInputSet): + """ + 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. + + Args: + structure (Structure): Input structure. + mode (str): Supported modes are "IPA", "RPA" + copy_wavecar (bool): Whether to copy the WAVECAR from a previous run. + Defaults to True. + 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 VaspInputSet. Typically, user_incar_settings is a + commonly used option. + """ + + # CONFIG = _load_yaml_config("MPAbsorptionSet") + + mode: str = "IPA" + copy_wavecar: bool = True + nbands_factor: float = 2 + reciprocal_density: float = 400 + nkred: tuple[int, int, int] | None = None + nedos: int = 2001 + inherit_incar: bool = True + force_gamma: bool = True + CONFIG = MPRelaxSet.CONFIG + nbands: int | None = None + SUPPORTED_MODES = ("IPA", "RPA") + + def __post_init__(self): + """Validate settings""" + super().__post_init__() + self.mode = self.mode.upper() + if self.mode not in MPAbsorptionSet.SUPPORTED_MODES: + raise ValueError(f"{self.mode} not one of the support modes : {MPAbsorptionSet.SUPPORTED_MODES}") + + @property + def kpoints_updates(self) -> dict | Kpoints: + """Updates to the kpoints configuration for this calculation type. + + 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 {"reciprocal_density": self.reciprocal_density} + + @property + def incar_updates(self) -> dict: + """Updates to the INCAR config for this calculation type.""" + updates = { + "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": self.nedos, + } + + if 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 + updates.update({"ALGO": "CHI", "NELM": 1, "NOMEGA": 1000, "EDIFF": None, "LOPTICS": None, "LWAVE": None}) + + if self.prev_vasprun is not None and self.mode == "IPA": + prev_nbands = int(self.prev_vasprun.parameters["NBANDS"]) if self.nbands is None else self.nbands + updates["NBANDS"] = int(np.ceil(prev_nbands * self.nbands_factor)) + + if self.prev_vasprun is not None and self.mode == "RPA": + # Since in the optical calculation, only the q->0 transition is of interest, + # we can reduce the number of q by the factor of the number of kpoints in + # each corresponding 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 lattices, but using NKREDX, NKREDY, NKREDZ are more sensible for + # other lattices. + self.nkred = self.prev_vasprun.kpoints.kpts[0] if self.nkred is None else self.nkred + updates.update({"NKREDX": self.nkred[0], "NKREDY": self.nkred[1], "NKREDZ": self.nkred[2]}) + + return updates diff --git a/pymatgen/io/vasp/sets/mvl.py b/pymatgen/io/vasp/sets/mvl.py new file mode 100644 index 00000000000..860dbcb58ff --- /dev/null +++ b/pymatgen/io/vasp/sets/mvl.py @@ -0,0 +1,443 @@ +from __future__ import annotations + +from dataclasses import dataclass +from typing import TYPE_CHECKING + +import numpy as np +from monty.json import MSONable + +from pymatgen.io.vasp.inputs import Kpoints +from pymatgen.io.vasp.sets.base import UserPotcarFunctional, VaspInputSet, _load_yaml_config +from pymatgen.io.vasp.sets.mit import MITRelaxSet +from pymatgen.io.vasp.sets.mp import MPRelaxSet +from pymatgen.util.due import Doi, due + +if TYPE_CHECKING: + from collections.abc import Sequence + + from typing_extensions import Self + + +@dataclass +class MVLGWSet(VaspInputSet): + """ + MVL denotes VASP input sets that are implemented by the Materials Virtual + Lab (http://materialsvirtuallab.org) for various research. This is a + flexible input set for GW calculations. + + Note that unlike all other input sets in this module, the PBE_54 series of + functional is set as the default. These have much improved performance for + GW calculations. + + A typical sequence is mode="STATIC" -> mode="DIAG" -> mode="GW" -> + mode="BSE". For all steps other than the first one (static), the + recommendation is to use from_prev_calculation on the preceding run in + the series. + + Args: + structure (Structure): Input structure. + mode (str): Supported modes are "STATIC" (default), "DIAG", "GW", + and "BSE". + 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 GW and BSE. This parameter is used by + from_previous_calculation to set nband. + copy_wavecar: Whether to copy the old WAVECAR, WAVEDER and associated + files when starting from a previous calculation. + nbands_factor (int): Multiplicative factor for NBANDS when starting + from a previous calculation. Only applies if mode=="DIAG". + Need to be tested for convergence. + reciprocal_density (int): Density of k-mesh by reciprocal atom. Only + applies if mode=="STATIC". Defaults to 100. + ncores (int): Numbers of cores used for the calculation. VASP will alter + NBANDS if it was not dividable by ncores. Only applies if + mode=="DIAG". + **kwargs: All kwargs supported by VaspInputSet. Typically, + user_incar_settings is a commonly used option. + """ + + reciprocal_density: float = 100 + mode: str = "STATIC" + copy_wavecar: bool = True + nbands_factor: int = 5 + ncores: int = 16 + nbands: int | None = None + force_gamma: bool = True + inherit_incar: bool = True # inherit incar from previous run if available + SUPPORTED_MODES = ("DIAG", "GW", "STATIC", "BSE") + CONFIG = _load_yaml_config("MVLGWSet") + + def __post_init__(self): + """Validate input settings.""" + super().__post_init__() + self.mode = mode = self.mode.upper() + + if mode not in MVLGWSet.SUPPORTED_MODES: + raise ValueError(f"Invalid {mode=}, supported modes are {', '.join(map(repr, MVLGWSet.SUPPORTED_MODES))}") + + @property + def kpoints_updates(self) -> dict: + """Updates to the kpoints configuration for this calculation type.""" + # Generate gamma center k-points mesh grid for GW calc, which is requested + # by GW calculation. + return {"reciprocal_density": self.reciprocal_density} + + @property + def incar_updates(self) -> dict: + """Updates to the INCAR config for this calculation type.""" + updates = {} + nbands = int(self.prev_vasprun.parameters["NBANDS"]) if self.prev_vasprun is not None else None + + if self.mode == "DIAG": + # Default parameters for diagonalization calculation. + updates.update({"ALGO": "Exact", "NELM": 1, "LOPTICS": True, "LPEAD": True}) + if nbands: + nbands = int(np.ceil(nbands * self.nbands_factor / self.ncores) * self.ncores) + + elif self.mode == "GW": + # Default parameters for GW calculation. + updates.update( + {"ALGO": "GW0", "NELM": 1, "NOMEGA": 80, "ENCUTGW": 250, "EDIFF": None, "LOPTICS": None, "LPEAD": None} + ) + elif self.mode == "BSE": + # Default parameters for BSE calculation. + updates.update({"ALGO": "BSE", "ANTIRES": 0, "NBANDSO": 20, "NBANDSV": 20}) + + if nbands: + updates["NBANDS"] = nbands + + return updates + + @classmethod + def from_prev_calc(cls, prev_calc_dir: str, mode: str = "DIAG", **kwargs) -> Self: + """Generate a set of VASP input files for GW or BSE calculations from a + directory of previous Exact Diag VASP run. + + Args: + prev_calc_dir (str): The directory contains the outputs( + vasprun.xml of previous vasp run. + mode (str): Supported modes are "STATIC", "DIAG" (default), "GW", + and "BSE". + **kwargs: All kwargs supported by MVLGWSet, other than structure, + prev_incar and mode, which are determined from the + prev_calc_dir. + """ + input_set = cls(None, mode=mode, **kwargs) + return input_set.override_from_prev_calc(prev_calc_dir=prev_calc_dir) + + +@dataclass +class MVLSlabSet(VaspInputSet): + """Write a set of slab vasp runs, including both slabs (along the c direction) + and orient unit cells (bulk), to ensure the same KPOINTS, POTCAR and INCAR criterion. + + Args: + structure: Structure + k_product: default to 50, kpoint number * length for a & b + directions, also for c direction in bulk calculations + bulk: + auto_dipole: + set_mix: + sort_structure: + **kwargs: Other kwargs supported by VaspInputSet. + """ + + k_product: int = 50 + bulk: bool = False + auto_dipole: bool = False + set_mix: bool = True + CONFIG = MPRelaxSet.CONFIG + + @property + def incar_updates(self) -> dict: + """Updates to the INCAR config for this calculation type.""" + updates = {"EDIFF": 1e-4, "EDIFFG": -0.02, "ENCUT": 400, "ISMEAR": 0, "SIGMA": 0.05, "ISIF": 3} + if not self.bulk: + updates.update({"ISIF": 2, "LVTOT": True, "NELMIN": 8}) + if self.set_mix: + updates.update({"AMIN": 0.01, "AMIX": 0.2, "BMIX": 0.001}) + if self.auto_dipole: + weights = [s.species.weight for s in self.structure] # type: ignore + center_of_mass = np.average(self.structure.frac_coords, weights=weights, axis=0) # type: ignore + updates.update({"IDIPOL": 3, "LDIPOL": True, "DIPOL": center_of_mass}) + return updates + + @property + def kpoints_updates(self): + """Updates to the kpoints configuration for this calculation type. + + k_product, default to 50, is kpoint number * length for a & b + directions, also for c direction in bulk calculations + Automatic mesh & Gamma is the default setting. + """ + # To get input sets, the input structure has to has the same number + # of required parameters as a Structure object (ie. 4). Slab + # attributes aren't going to affect the VASP inputs anyways so + # converting the slab into a structure should not matter + # use k_product to calculate kpoints, k_product = kpts[0][0] * a + lattice_abc = self.structure.lattice.abc + kpt_calc = [ + int(self.k_product / lattice_abc[0] + 0.5), + int(self.k_product / lattice_abc[1] + 0.5), + 1, + ] + + # calculate kpts (c direction) for bulk. (for slab, set to 1) + if self.bulk: + kpt_calc[2] = int(self.k_product / lattice_abc[2] + 0.5) + + return Kpoints(comment="Generated by pymatgen's MVLGBSet", style=Kpoints.supported_modes.Gamma, kpts=[kpt_calc]) + + def as_dict(self, verbosity=2): + """ + Args: + verbosity (int): Verbosity of dict. e.g. whether to include Structure. + + Returns: + dict: MSONable MVLGBSet representation. + """ + dct = MSONable.as_dict(self) + if verbosity == 1: + dct.pop("structure", None) + return dct + + +@dataclass +class MVLGBSet(VaspInputSet): + """Write a vasp input files for grain boundary calculations, slab or bulk. + + Args: + structure (Structure): provide the structure + k_product: Kpoint number * length for a & b directions, also for c direction in + bulk calculations. Default to 40. + slab_mode (bool): Defaults to False. Use default (False) for a bulk supercell. + Use True if you are performing calculations on a slab-like (i.e., surface) + of the GB, for example, when you are calculating the work of separation. + is_metal (bool): Defaults to True. This determines whether an ISMEAR of 1 is + used (for metals) or not (for insulators and semiconductors) by default. + Note that it does *not* override user_incar_settings, which can be set by + the user to be anything desired. + **kwargs: + Other kwargs supported by MPRelaxSet. + """ + + k_product: int = 40 + slab_mode: bool = False + is_metal: bool = True + CONFIG = MPRelaxSet.CONFIG + + @property + def kpoints_updates(self): + """k_product is kpoint number * length for a & b directions, also for c direction + in bulk calculations Automatic mesh & Gamma is the default setting. + """ + # use k_product to calculate kpoints, k_product = kpts[0][0] * a + lengths = self.structure.lattice.abc + kpt_calc = [ + int(self.k_product / lengths[0] + 0.5), + int(self.k_product / lengths[1] + 0.5), + int(self.k_product / lengths[2] + 0.5), + ] + + if self.slab_mode: + kpt_calc[2] = 1 + + return Kpoints(comment="Generated by pymatgen's MVLGBSet", style=Kpoints.supported_modes.Gamma, kpts=[kpt_calc]) + + @property + def incar_updates(self) -> dict: + """Updates to the INCAR config for this calculation type.""" + # The default incar setting is used for metallic system, for + # insulator or semiconductor, ISMEAR need to be changed. + updates = dict(LCHARG=False, NELM=60, PREC="Normal", EDIFFG=-0.02, ICHARG=0, NSW=200, EDIFF=0.0001) + + if self.is_metal: + updates["ISMEAR"] = 1 + updates["LDAU"] = False + + if self.slab_mode: + # for clean grain boundary and bulk relaxation, full optimization + # relaxation (ISIF=3) is used. For slab relaxation (ISIF=2) is used. + updates["ISIF"] = 2 + updates["NELMIN"] = 8 + + return updates + + +@dataclass +class MVLRelax52Set(VaspInputSet): + """ + Implementation of VaspInputSet utilizing the public Materials Project + parameters for INCAR & KPOINTS and VASP's recommended PAW potentials for + POTCAR. + + Keynotes from VASP manual: + 1. Recommended potentials for calculations using vasp.5.2+ + 2. If dimers with short bonds are present in the compound (O2, CO, + N2, F2, P2, S2, Cl2), it is recommended to use the h potentials. + Specifically, C_h, O_h, N_h, F_h, P_h, S_h, Cl_h + 3. Released on Oct 28, 2018 by VASP. Please refer to VASP + Manual 1.2, 1.3 & 10.2.1 for more details. + + Args: + structure (Structure): input structure. + user_potcar_functional (str): choose from "PBE_52" and "PBE_54". + **kwargs: Other kwargs supported by VaspInputSet. + """ + + user_potcar_functional: UserPotcarFunctional = "PBE_52" + CONFIG = _load_yaml_config("MVLRelax52Set") + _valid_potcars: Sequence[str] | None = ("PBE_52", "PBE_54") + + +@due.dcite( + Doi("10.1149/2.0061602jes"), + description="Elastic Properties of Alkali Superionic Conductor Electrolytes from First Principles Calculations", +) +class MVLElasticSet(VaspInputSet): + """ + MVL denotes VASP input sets that are implemented by the Materials Virtual + Lab (http://materialsvirtuallab.org) for various research. + + This input set is used to calculate elastic constants in VASP. It is used + in the following work:: + + Z. Deng, Z. Wang, I.-H. Chu, J. Luo, S. P. Ong. + “Elastic Properties of Alkali Superionic Conductor Electrolytes + from First Principles Calculations”, J. Electrochem. Soc. + 2016, 163(2), A67-A74. doi: 10.1149/2.0061602jes + + To read the elastic constants, you may use the Outcar class which parses the + elastic constants. + + Args: + structure (pymatgen.Structure): Input structure. + potim (float): POTIM parameter. The default of 0.015 is usually fine, + but some structures may require a smaller step. + **kwargs: Parameters supported by MPRelaxSet. + """ + + potim: float = 0.015 + CONFIG = MPRelaxSet.CONFIG + + @property + def incar_updates(self) -> dict: + """Updates to the INCAR config for this calculation type.""" + return {"IBRION": 6, "NFREE": 2, "POTIM": self.potim, "NPAR": None} + + +@dataclass +class MVLNPTMDSet(VaspInputSet): + """Write a VASP MD run in NPT ensemble. + + Notes: + To eliminate Pulay stress, the default ENCUT is set to a rather large + value of ENCUT, which is 1.5 * ENMAX. + """ + + start_temp: float = 0.0 + end_temp: float = 300.0 + nsteps: int = 1000 + time_step: float = 2 + spin_polarized: bool = False + CONFIG = MITRelaxSet.CONFIG + + @property + def incar_updates(self) -> dict: + """Updates to the INCAR config for this calculation type.""" + # NPT-AIMD default settings + updates = { + "ALGO": "Fast", + "ISIF": 3, + "LANGEVIN_GAMMA": [10] * self.structure.ntypesp, # type: ignore + "LANGEVIN_GAMMA_L": 1, + "MDALGO": 3, + "PMASS": 10, + "PSTRESS": 0, + "SMASS": 0, + "TEBEG": self.start_temp, + "TEEND": self.end_temp, + "NSW": self.nsteps, + "EDIFF_PER_ATOM": 0.000001, + "LSCALU": False, + "LCHARG": False, + "LPLANE": False, + "LWAVE": True, + "ISMEAR": 0, + "NELMIN": 4, + "LREAL": True, + "BMIX": 1, + "MAXMIX": 20, + "NELM": 500, + "NSIM": 4, + "ISYM": 0, + "IBRION": 0, + "NBLOCK": 1, + "KBLOCK": 100, + "POTIM": self.time_step, + "PREC": "Low", + "ISPIN": 2 if self.spin_polarized else 1, + "LDAU": False, + } + # Set NPT-AIMD ENCUT = 1.5 * VASP_default + enmax = [self.potcar[i].keywords["ENMAX"] for i in range(self.structure.ntypesp)] # type: ignore + updates["ENCUT"] = max(enmax) * 1.5 + return updates + + @property + def kpoints_updates(self) -> Kpoints | dict: + """Updates to the kpoints configuration for this calculation type.""" + return Kpoints.gamma_automatic() + + +@dataclass +class MVLScanRelaxSet(VaspInputSet): + """Write a relax input set using Strongly Constrained and + Appropriately Normed (SCAN) semilocal density functional. + + Notes: + 1. This functional is only available from VASP.5.4.3 upwards. + + 2. Meta-GGA calculations require POTCAR files that include + information on the kinetic energy density of the core-electrons, + i.e. "PBE_52" or "PBE_54". Make sure the POTCAR including the + following lines (see VASP wiki for more details): + + $ grep kinetic POTCAR + kinetic energy-density + mkinetic energy-density pseudized + kinetic energy density (partial) + + Args: + structure (Structure): input structure. + vdw (str): set "rVV10" to enable SCAN+rVV10, which is a versatile + van der Waals density functional by combing the SCAN functional + with the rVV10 non-local correlation functional. + **kwargs: Other kwargs supported by VaspInputSet. + """ + + user_potcar_functional: UserPotcarFunctional = "PBE_52" + _valid_potcars: Sequence[str] | None = ("PBE_52", "PBE_54") + CONFIG = MPRelaxSet.CONFIG + + def __post_init__(self): + super().__post_init__() + if self.user_potcar_functional not in ("PBE_52", "PBE_54"): + raise ValueError("SCAN calculations require PBE_52 or PBE_54!") + + @property + def incar_updates(self) -> dict: + """Updates to the INCAR config for this calculation type.""" + updates = { + "ADDGRID": True, + "EDIFF": 1e-5, + "EDIFFG": -0.05, + "LASPH": True, + "LDAU": False, + "METAGGA": "SCAN", + "NELM": 200, + } + if self.vdw and self.vdw.lower() == "rvv10": + updates["BPARAM"] = 15.7 # This is the correct BPARAM for SCAN+rVV10 + return updates