diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index b562b202756..450952dd863 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -25,6 +25,7 @@ repos: rev: v1.11.1 hooks: - id: mypy + entry: env MYPYPATH=src mypy - repo: https://github.com/codespell-project/codespell rev: v2.3.0 diff --git a/pyproject.toml b/pyproject.toml index cb62f430144..15a58446165 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -174,6 +174,7 @@ repair-wheel-command = "delocate-wheel --require-archs {delocate_archs} -w {dest [tool.ruff] target-version = "py39" line-length = 120 +output-format = "concise" [tool.ruff.lint] select = ["ALL"] diff --git a/src/pymatgen/io/vasp/sets.py b/src/pymatgen/io/vasp/sets.py deleted file mode 100644 index 08ce0dedb18..00000000000 --- a/src/pymatgen/io/vasp/sets.py +++ /dev/null @@ -1,3205 +0,0 @@ -""" -This module defines the VaspInputSet abstract base class and a concrete implementation for the parameters developed -and tested by the core team of pymatgen, including the Materials Virtual Lab, Materials Project and the MIT high -throughput project. The basic concept behind an input set is to specify a scheme to generate a consistent set of VASP -inputs from a structure without further user intervention. This ensures comparability across runs. - -Read the following carefully before implementing new input sets: - -1. 99% of what needs to be done can be done by specifying user_incar_settings to override some of the defaults of - various input sets. Unless there is an extremely good reason to add a new set, **do not** add one. e.g. if you want - to turn the Hubbard U off, just set "LDAU": False as a user_incar_setting. -2. All derivative input sets should inherit appropriate configurations (e.g., from MPRelaxSet), and more often than - not, VaspInputSet should be the superclass. Superclass delegation should be used where possible. In particular, - you are not supposed to implement your own as_dict or from_dict for derivative sets unless you know what you are - doing. Improper overriding the as_dict and from_dict protocols is the major cause of implementation headaches. If - you need an example, look at how the MPStaticSet is initialized. - -The above are recommendations. The following are **UNBREAKABLE** rules: - -1. All input sets must take in a structure, list of structures or None as the first argument. If None, the input set - should perform a stateless initialization and before any output can be written, a structure must be set. -2. user_incar_settings, user_kpoints_settings and user__settings are ABSOLUTE. Any new sets you implement - must obey this. If a user wants to override your settings, you assume he knows what he is doing. Do not - magically override user supplied settings. You can issue a warning if you think the user is wrong. -3. All input sets must save all supplied args and kwargs as instance variables. e.g. self.arg = arg and - self.kwargs = kwargs in the __init__. This ensures the as_dict and from_dict work correctly. -""" - -from __future__ import annotations - -import abc -import itertools -import os -import re -import warnings -from collections.abc import Sequence -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, Any, cast - -import numpy as np -from monty.dev import deprecated -from monty.json import MSONable -from monty.serialization import loadfn - -from pymatgen.analysis.structure_matcher import StructureMatcher -from pymatgen.core import Element, PeriodicSite, SiteCollection, Species, 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 Callable, Literal, Union - - from typing_extensions import Self - - from pymatgen.util.typing import PathLike, Tuple3Ints, Vector3D - - UserPotcarFunctional = Union[ - Literal["PBE", "PBE_52", "PBE_54", "LDA", "LDA_52", "LDA_54", "PW91", "LDA_US", "PW91_US"], None - ] - -MODULE_DIR = os.path.dirname(__file__) - - -def _load_yaml_config(fname): - config = loadfn(f"{MODULE_DIR}/{fname}.yaml") - if "PARENT" in config: - parent_config = _load_yaml_config(config["PARENT"]) - for k, v in parent_config.items(): - if k not in config: - config[k] = v - elif isinstance(v, dict): - v_new = config.get(k, {}) - v_new.update(v) - config[k] = v_new - return config - - -@dataclass -class VaspInputSet(InputGenerator, abc.ABC): - """ - Base class representing a set of VASP input parameters with a structure - supplied as init parameters and initialized from a dict of settings. - This allows arbitrary settings to be input. In general, - this is rarely used directly unless there is a source of settings in yaml - format (e.g., from a REST interface). It is typically used by other - VaspInputSets for initialization. - - Special consideration should be paid to the way the MAGMOM initialization - for the INCAR is done. The initialization differs depending on the type of - structure and the configuration settings. The order in which the magmom is - determined is as follows: - - 1. If the site is specified in user_incar_settings, use that setting. - 2. If the site itself has a magmom setting (i.e. site.properties["magmom"] = float), - that is used. This can be set with structure.add_site_property(). - 3. If the species of the site has a spin setting, that is used. This can be set - with structure.add_spin_by_element(). - 4. If the species itself has a particular setting in the config file, that - is used, e.g. Mn3+ may have a different magmom than Mn4+. - 5. Lastly, the element symbol itself is checked in the config file. If - there are no settings, a default value of 0.6 is used. - - 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. - config_dict (dict): The config dictionary to use. - files_to_transfer (dict): A dictionary of {filename: filepath}. This allows the - transfer of files from a previous calculation. - user_incar_settings (dict): User INCAR settings. This allows a user to override - INCAR settings, e.g. setting a different MAGMOM for various elements or - species. Note that in the new scheme, ediff_per_atom and hubbard_u are no - longer args. Instead, the CONFIG supports EDIFF_PER_ATOM and EDIFF keys. - The former scales with # of atoms, the latter does not. If both are present, - EDIFF is preferred. To force such settings, just supply - user_incar_settings={"EDIFF": 1e-5, "LDAU": False} for example. The keys - 'LDAUU', 'LDAUJ', 'LDAUL' are special cases since pymatgen defines different - values depending on what anions are present in the structure, so these keys - can be defined in one of two ways, e.g. either {"LDAUU":{"O":{"Fe":5}}} to - set LDAUU for Fe to 5 in an oxide, or {"LDAUU":{"Fe":5}} to set LDAUU to 5 - regardless of the input structure. If a None value is given, that key is - unset. For example, {"ENCUT": None} will remove ENCUT from the - incar settings. Finally, KSPACING is a special setting and can be set to - "auto" in which the KSPACING is set automatically based on the band gap. - user_kpoints_settings (dict or Kpoints): Allow user to override kpoints setting - by supplying a dict. e.g. {"reciprocal_density": 1000}. User can also - supply Kpoints object. - user_potcar_settings (dict): Allow user to override POTCARs. e.g. {"Gd": - "Gd_3"}. This is generally not recommended. - constrain_total_magmom (bool): Whether to constrain the total magmom (NUPDOWN in - INCAR) to be the sum of the expected MAGMOM for all species. - sort_structure (bool): Whether to sort the structure (using the default sort - order of electronegativity) before generating input files. Defaults to True, - the behavior you would want most of the time. This ensures that similar - atomic species are grouped together. - user_potcar_functional (str): Functional to use. Default (None) is to use the - functional in the config dictionary. Valid values: "PBE", "PBE_52", - "PBE_54", "LDA", "LDA_52", "LDA_54", "PW91", "LDA_US", "PW91_US". - force_gamma (bool): Force gamma centered kpoint generation. Default (False) is - to use the Automatic Density kpoint scheme, which will use the Gamma - centered generation scheme for hexagonal cells, and Monkhorst-Pack otherwise. - reduce_structure (None/str): Before generating the input files, generate the - reduced structure. Default (None), does not alter the structure. Valid - values: None, "niggli", "LLL". - vdw: Adds default parameters for van-der-Waals functionals supported by VASP to - INCAR. Supported functionals are: DFT-D2, undamped DFT-D3, DFT-D3 with - Becke-Jonson damping, Tkatchenko-Scheffler, Tkatchenko-Scheffler with - iterative Hirshfeld partitioning, MBD@rSC, dDsC, Dion's vdW-DF, DF2, optPBE, - optB88, optB86b and rVV10. - use_structure_charge (bool): If set to True, then the overall charge of the - structure (structure.charge) is used to set the NELECT variable in the - INCAR. Default is False. - standardize (float): Whether to standardize to a primitive standard cell. - Defaults to False. - sym_prec (float): Tolerance for symmetry finding. - international_monoclinic (bool): Whether to use international convention (vs - Curtarolo) for monoclinic. Defaults True. - validate_magmom (bool): Ensure that the missing magmom values are filled in with - the VASP default value of 1.0. - inherit_incar (bool): Whether to inherit INCAR settings from previous - calculation. This might be useful to port Custodian fixes to child jobs but - can also be dangerous e.g. when switching from GGA to meta-GGA or relax to - static jobs. Defaults to True. - auto_kspacing (bool): If true, determines the value of KSPACING from the bandgap - of a previous calculation. - auto_ismear (bool): If true, the values for ISMEAR and SIGMA will be set - automatically depending on the bandgap of the system. If the bandgap is not - known (e.g., there is no previous VASP directory) then ISMEAR=0 and - SIGMA=0.2; if the bandgap is zero (a metallic system) then ISMEAR=2 and - SIGMA=0.2; if the system is an insulator, then ISMEAR=-5 (tetrahedron - smearing). Note, this only works when generating the input set from a - previous VASP directory. - auto_ispin (bool) = False: - If generating input set from a previous calculation, this controls whether - to disable magnetisation (ISPIN = 1) if the absolute value of all magnetic - moments are less than 0.02. - auto_lreal (bool) = False: - If True, automatically use the VASP recommended LREAL based on cell size. - auto_metal_kpoints - If true and the system is metallic, try and use ``reciprocal_density_metal`` - instead of ``reciprocal_density`` for metallic systems. Note, this only works - if the bandgap is not None. - bandgap_tol (float): Tolerance for determining if a system is metallic when - KSPACING is set to "auto". If the bandgap is less than this value, the - system is considered metallic. Defaults to 1e-4 (eV). - bandgap (float): Used for determining KSPACING if KSPACING == "auto" or - ISMEAR if auto_ismear == True. Set automatically when using from_prev_calc. - prev_incar (str or dict): Previous INCAR used for setting parent INCAR when - inherit_incar == True. Set automatically when using from_prev_calc. - prev_kpoints (str or Kpoints): Previous Kpoints. Set automatically when using - from_prev_calc. - """ - - structure: Structure | None = None - config_dict: dict = field(default_factory=dict) - files_to_transfer: dict = field(default_factory=dict) - user_incar_settings: dict = field(default_factory=dict) - user_kpoints_settings: dict = field(default_factory=dict) - user_potcar_settings: dict = field(default_factory=dict) - constrain_total_magmom: bool = False - sort_structure: bool = True - user_potcar_functional: UserPotcarFunctional = None - force_gamma: bool = False - reduce_structure: Literal["niggli", "LLL"] | None = None - vdw: str | None = None - use_structure_charge: bool = False - standardize: bool = False - sym_prec: float = 0.1 - international_monoclinic: bool = True - validate_magmom: bool = True - inherit_incar: bool | list[str] = False - auto_kspacing: bool = False - auto_ismear: bool = False - auto_ispin: bool = False - auto_lreal: bool = False - auto_metal_kpoints: bool = False - bandgap_tol: float = 1e-4 - bandgap: float | None = None - prev_incar: str | dict | None = None - prev_kpoints: str | Kpoints | None = None - _valid_potcars: Sequence[str] | None = None - - def __post_init__(self) -> None: - """Perform validation.""" - user_potcar_functional = self.user_potcar_functional - if (valid_potcars := self._valid_potcars) and user_potcar_functional not in valid_potcars: - raise ValueError(f"Invalid {user_potcar_functional=}, must be one of {valid_potcars}") - - if hasattr(self, "CONFIG"): - self.config_dict = self.CONFIG - - self._config_dict = deepcopy(self.config_dict) - - # These have been left to stay consistent with previous API - self.user_incar_settings = self.user_incar_settings or {} - self.user_kpoints_settings = self.user_kpoints_settings or {} - - self.vdw = self.vdw.lower() if isinstance(self.vdw, str) else self.vdw - if self.user_incar_settings.get("KSPACING") and self.user_kpoints_settings: - # self.user_kpoints_settings will never be `None` because it is set to - # an empty dict if it is `None`. - warnings.warn( - "You have specified KSPACING and also supplied KPOINTS " - "settings. KSPACING only has effect when there is no " - "KPOINTS file. Since both settings were given, pymatgen" - "will generate a KPOINTS file and ignore KSPACING." - "Remove the `user_kpoints_settings` argument to enable KSPACING.", - BadInputSetWarning, - ) - - if self.vdw: - vdw_par = loadfn(f"{MODULE_DIR}/vdW_parameters.yaml") - if vdw_param := vdw_par.get(self.vdw): - self._config_dict["INCAR"].update(vdw_param) - else: - raise KeyError( - f"Invalid or unsupported van-der-Waals functional. Supported functionals are {', '.join(vdw_par)}." - ) - # 'or' case reads the POTCAR_FUNCTIONAL from the .yaml - self.user_potcar_functional: UserPotcarFunctional = self.user_potcar_functional or self._config_dict.get( - "POTCAR_FUNCTIONAL", "PBE" - ) - - # Warn if a user is overriding POTCAR_FUNCTIONAL - if self.user_potcar_functional != self._config_dict.get("POTCAR_FUNCTIONAL", "PBE"): - warnings.warn( - "Overriding the POTCAR functional is generally not recommended " - " as it significantly affect the results of calculations and " - "compatibility with other calculations done with the same " - "input set. Note that some POTCAR symbols specified in " - "the configuration file may not be available in the selected " - "functional.", - BadInputSetWarning, - ) - - if self.user_potcar_settings: - warnings.warn( - "Overriding POTCARs is generally not recommended as it " - "significantly affect the results of calculations and " - "compatibility with other calculations done with the same " - "input set. In many instances, it is better to write a " - "subclass of a desired input set and override the POTCAR in " - "the subclass to be explicit on the differences.", - BadInputSetWarning, - ) - for key, val in self.user_potcar_settings.items(): - self._config_dict["POTCAR"][key] = val - - if not isinstance(self.structure, Structure): - self._structure: Structure | None = None - else: - self.structure = self.structure - - if isinstance(self.prev_incar, (Path, str)): - self.prev_incar = Incar.from_file(self.prev_incar) - - if isinstance(self.prev_kpoints, (Path, str)): - self.prev_kpoints = Kpoints.from_file(self.prev_kpoints) - - self.prev_vasprun: Vasprun | None = None - self.prev_outcar: Outcar | None = None - self._ispin: Literal[1, 2] | None = None - - def __str__(self) -> str: - return type(self).__name__ - - def __repr__(self) -> str: - return type(self).__name__ - - def write_input( - self, - output_dir: str, - make_dir_if_not_present: bool = True, - include_cif: bool | str = False, - potcar_spec: bool = False, - zip_output: bool | str = False, - ) -> None: - """Write a set of VASP input to a directory. - - 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. - include_cif (bool): Whether to write a CIF file in the output - directory for easier opening by VESTA. - potcar_spec (bool): Instead of writing the POTCAR, write a "POTCAR.spec". - This is intended to help sharing an input set with people who might - not have a license to specific Potcar files. Given a "POTCAR.spec", - the specific POTCAR file can be re-generated using pymatgen with the - "generate_potcar" function in the pymatgen CLI. - zip_output (bool): If True, output will be zipped into a file with the - same name as the InputSet (e.g., MPStaticSet.zip). - """ - vasp_input = self.get_input_set(potcar_spec=potcar_spec) - - cif_name = None - if include_cif: - struct = vasp_input["POSCAR"].structure - cif_name = f"{output_dir}/{struct.formula.replace(' ', '')}.cif" - - vasp_input.write_input( - output_dir=output_dir, - make_dir_if_not_present=make_dir_if_not_present, - cif_name=cif_name, - zip_name=f"{type(self).__name__}.zip" if zip_output else None, - files_to_transfer=self.files_to_transfer, - ) - - def as_dict(self, verbosity: int = 2) -> dict: - """ - Args: - verbosity: Verbosity for generated dict. If 1, structure is - excluded. - - Returns: - dict: MSONable VaspInputSet representation. - """ - dct = MSONable.as_dict(self) - if verbosity == 1: - dct.pop("structure", None) - return dct - - @property # type: ignore[no-redef] - def structure(self) -> Structure | None: # noqa: F811 - """Structure.""" - return self._structure - - @structure.setter - def structure(self, structure: Structure | None) -> None: - if not hasattr(self, "_config_dict"): - self._structure = structure - return - - if isinstance(structure, SiteCollection): # could be Structure or Molecule - if self.user_potcar_functional == "PBE_54" and "W" in structure.symbol_set: - # When using 5.4 POTCARs, default Tungsten POTCAR to W_Sv but still allow user to override - self.user_potcar_settings = {"W": "W_sv", **(self.user_potcar_settings or {})} - if self.reduce_structure: - structure = structure.get_reduced_structure(self.reduce_structure) - if self.sort_structure: - structure = structure.get_sorted_structure() - if self.validate_magmom: - get_valid_magmom_struct(structure, spin_mode="auto") - - struct_has_Yb = any(specie.symbol == "Yb" for site in structure for specie in site.species) - potcar_settings = self._config_dict.get("POTCAR", {}) - if self.user_potcar_settings: - potcar_settings.update(self.user_potcar_settings) - uses_Yb_2_psp = potcar_settings.get("Yb", None) == "Yb_2" - if struct_has_Yb and uses_Yb_2_psp: - warnings.warn( - "The structure contains Ytterbium (Yb) and this InputSet uses the Yb_2 PSP.\n" - "Yb_2 is known to often give bad results since Yb has oxidation state 3+ in most compounds.\n" - "See https://github.com/materialsproject/pymatgen/issues/2968 for details.", - BadInputSetWarning, - ) - if self.standardize and self.sym_prec: - structure = standardize_structure( - structure, - sym_prec=self.sym_prec, - international_monoclinic=self.international_monoclinic, - ) - self._structure = structure - - def get_input_set( - self, - structure: Structure | None = None, - prev_dir: PathLike | None = None, - potcar_spec: bool = False, - ) -> VaspInput: - """Get a VASP input set. - - Note, if both ``structure`` and ``prev_dir`` are set, then the structure - specified will be preferred over the final structure from the last VASP run. - - Args: - structure (Structure): A structure. - prev_dir (PathLike): A previous directory to generate the input set from. - potcar_spec (bool): Instead of generating a Potcar object, use a list of - potcar symbols. This will be written as a "POTCAR.spec" file. This is - intended to help sharing an input set with people who might not have a - license to specific Potcar files. Given a "POTCAR.spec", the specific - POTCAR file can be re-generated using pymatgen with the - "generate_potcar" function in the pymatgen CLI. - - Returns: - VaspInput: A VASP input object. - """ - if structure is None and prev_dir is None and self.structure is None: - raise ValueError("Either structure or prev_dir must be set") - - self._set_previous(prev_dir) - - if structure is not None: - self.structure = structure - - return VaspInput( - incar=self.incar, - kpoints=self.kpoints, - poscar=self.poscar, - potcar="\n".join(self.potcar_symbols) if potcar_spec else self.potcar, - potcar_spec=potcar_spec, - ) - - @deprecated(get_input_set, deadline=(2026, 6, 6)) - def get_vasp_input(self, structure: Structure | None = None) -> Self: - """Get a VaspInput object. - - Returns: - VaspInput. - """ - return self.get_input_set(structure=structure) - - @property - def incar_updates(self) -> dict: - """Updates to the INCAR config for this calculation type.""" - return {} - - @property - def kpoints_updates(self) -> dict: - """Updates to the kpoints configuration for this calculation type. - - Note, these updates will be ignored if the user has set user_kpoint_settings. - - Returns: - dict or Kpoints: A dictionary of updates to apply to the KPOINTS config - or a Kpoints object. - """ - return {} - - def _set_previous(self, prev_dir: PathLike | None = None) -> None: - """Load previous calculation outputs.""" - if prev_dir is None: - return - - vasprun, outcar = get_vasprun_outcar(prev_dir) - self.prev_vasprun = vasprun - self.prev_outcar = outcar - self.prev_incar = vasprun.incar - self.prev_kpoints = Kpoints.from_dict(vasprun.kpoints.as_dict()) - - if vasprun.efermi is None: - # VASP doesn't output efermi in vasprun if IBRION = 1 - vasprun.efermi = outcar.efermi - - bs = vasprun.get_band_structure(efermi="smart") - self.bandgap = 0 if bs.is_metal() else bs.get_band_gap()["energy"] - if self.auto_ispin: - # Turn off spin when magmom for every site is smaller than 0.02. - self._ispin = _get_ispin(vasprun, outcar) - - self.structure = get_structure_from_prev_run(vasprun, outcar) - - @property - def incar(self) -> Incar: - """The INCAR.""" - if self.structure is None: - raise RuntimeError("No structure is associated with the input set!") - - prev_incar: dict[str, Any] = {} - if self.inherit_incar is True and self.prev_incar: - prev_incar = cast(dict[str, Any], self.prev_incar) - elif isinstance(self.inherit_incar, (list, tuple)) and self.prev_incar: - prev_incar = { - k: cast(dict[str, Any], self.prev_incar)[k] for k in self.inherit_incar if k in self.prev_incar - } - - incar_updates = self.incar_updates - settings = dict(self._config_dict["INCAR"]) - auto_updates: dict[str, Any] = {} - if self.auto_ispin and (self._ispin is not None): - auto_updates["ISPIN"] = self._ispin - - # Breaking change - order in which settings applied inconsistent with atomate2 - # apply updates from input set generator to SETTINGS - # _apply_incar_updates(settings, incar_updates) - - # Apply user incar settings to SETTINGS not to INCAR - _apply_incar_updates(settings, self.user_incar_settings) - - # Generate INCAR - structure = self.structure - comp = structure.composition - elements = sorted((el for el in comp.elements if comp[el] > 0), key=lambda e: e.X) - most_electro_neg = elements[-1].symbol - poscar = Poscar(structure) - hubbard_u = settings.get("LDAU", False) - incar = Incar() - - for key, setting in settings.items(): - if key == "MAGMOM": - mag = [] - for site in structure: - if uic_magmom := self.user_incar_settings.get("MAGMOM", {}).get(site.species_string): - mag.append(uic_magmom) - elif hasattr(site, "magmom"): - mag.append(site.magmom) - elif getattr(site.specie, "spin", None) is not None: - mag.append(site.specie.spin) - elif str(site.specie) in setting: - if site.specie.symbol == "Co" and setting[str(site.specie)] <= 1.0: - warnings.warn( - "Co without an oxidation state is initialized as low spin by default in Pymatgen. " - "If this default behavior is not desired, please set the spin on the magmom on the " - "site directly to ensure correct initialization." - ) - mag.append(setting.get(str(site.specie))) - else: - if site.specie.symbol == "Co": - warnings.warn( - "Co without an oxidation state is initialized as low spin by default in Pymatgen. " - "If this default behavior is not desired, please set the spin on the magmom on the " - "site directly to ensure correct initialization." - ) - mag.append(setting.get(site.specie.symbol, 0.6)) - incar[key] = mag - - elif key in {"LDAUU", "LDAUJ", "LDAUL"}: - if hubbard_u: - if hasattr(structure[0], key.lower()): - m = {site.specie.symbol: getattr(site, key.lower()) for site in structure} - incar[key] = [m[sym] for sym in poscar.site_symbols] - # Lookup specific LDAU if specified for most_electroneg atom - elif most_electro_neg in setting and isinstance(setting[most_electro_neg], dict): - incar[key] = [setting[most_electro_neg].get(sym, 0) for sym in poscar.site_symbols] - # Else, use fallback LDAU value if it exists - else: - incar[key] = [ - setting.get(sym, 0) if isinstance(setting.get(sym, 0), (float, int)) else 0 - for sym in poscar.site_symbols - ] - - elif key.startswith("EDIFF") and key != "EDIFFG": - if "EDIFF" not in settings and key == "EDIFF_PER_ATOM": - incar["EDIFF"] = float(setting) * len(structure) - else: - incar["EDIFF"] = float(settings["EDIFF"]) - - elif key == "KSPACING" and self.auto_kspacing: - # Default to metal if no prev calc available - bandgap = 0 if self.bandgap is None else self.bandgap - incar[key] = auto_kspacing(bandgap, self.bandgap_tol) - - else: - incar[key] = setting - - has_u = hubbard_u and sum(incar["LDAUU"]) > 0 - if not has_u: - for key in list(incar): - if key.startswith("LDAU"): - del incar[key] - - # Modify LMAXMIX if you have d or f electrons present. Note that if the user - # explicitly sets LMAXMIX in settings it will override this logic. - # Previously, this was only set if Hubbard U was enabled as per the VASP manual - # but following an investigation it was determined that this would lead to a - # significant difference between SCF -> NonSCF even without Hubbard U enabled. - # Thanks to Andrew Rosen for investigating and reporting. - if "LMAXMIX" not in settings: - # contains f-electrons - if any(el.Z > 56 for el in structure.composition): - incar["LMAXMIX"] = 6 - # contains d-electrons - elif any(el.Z > 20 for el in structure.composition): - incar["LMAXMIX"] = 4 - - # Warn user about LASPH for +U, meta-GGAs, hybrids, and vdW-DF - if not incar.get("LASPH", False) and ( - incar.get("METAGGA") - or incar.get("LHFCALC", False) - or incar.get("LDAU", False) - or incar.get("LUSE_VDW", False) - ): - warnings.warn("LASPH = True should be set for +U, meta-GGAs, hybrids, and vdW-DFT", BadInputSetWarning) - - # Apply previous INCAR settings, be careful not to override user_incar_settings - # or the settings updates from the specific input set implementations - # also skip LDAU/MAGMOM as structure may have changed. - skip = list(self.user_incar_settings) + list(incar_updates) - skip += ["MAGMOM", "NUPDOWN", "LDAUU", "LDAUL", "LDAUJ"] - _apply_incar_updates(incar, prev_incar, skip=skip) - - if self.constrain_total_magmom: - nupdown = sum(mag if abs(mag) > 0.6 else 0 for mag in incar["MAGMOM"]) - if abs(nupdown - round(nupdown)) > 1e-5: - warnings.warn( - "constrain_total_magmom was set to True, but the sum of MAGMOM " - "values is not an integer. NUPDOWN is meant to set the spin " - "multiplet and should typically be an integer. You are likely " - "better off changing the values of MAGMOM or simply setting " - "NUPDOWN directly in your INCAR settings.", - UserWarning, - stacklevel=1, - ) - auto_updates["NUPDOWN"] = nupdown - - if self.use_structure_charge: - auto_updates["NELECT"] = self.nelect - - # Check that ALGO is appropriate - if incar.get("LHFCALC", False) is True and incar.get("ALGO", "Normal") not in ["Normal", "All", "Damped"]: - warnings.warn( - "Hybrid functionals only support Algo = All, Damped, or Normal.", - BadInputSetWarning, - ) - - if self.auto_ismear: - if self.bandgap is None: - # Don't know if we are a metal or insulator so set ISMEAR and SIGMA to - # be safe with the most general settings - auto_updates.update(ISMEAR=0, SIGMA=0.2) - elif self.bandgap <= self.bandgap_tol: - auto_updates.update(ISMEAR=2, SIGMA=0.2) # metal - else: - auto_updates.update(ISMEAR=-5, SIGMA=0.05) # insulator - - if self.auto_lreal: - auto_updates.update(LREAL=_get_recommended_lreal(structure)) - - # Apply updates from auto options, careful not to override user_incar_settings - _apply_incar_updates(incar, auto_updates, skip=list(self.user_incar_settings)) - - # Apply updates from input set generator to INCAR - _apply_incar_updates(incar, incar_updates, skip=list(self.user_incar_settings)) - - # Finally, re-apply `self.user_incar_settings` to make sure any accidentally - # overwritten settings are changed back to the intended values. - # skip dictionary parameters to avoid dictionaries appearing in the INCAR - _apply_incar_updates(incar, self.user_incar_settings, skip=["LDAUU", "LDAUJ", "LDAUL", "MAGMOM"]) - - # Remove unused INCAR parameters - _remove_unused_incar_params(incar, skip=list(self.user_incar_settings)) - - kpoints = self.kpoints - if kpoints is not None: - # Unset KSPACING as we are using a KPOINTS file - incar.pop("KSPACING", None) - - elif "KSPACING" in incar and "KSPACING" not in self.user_incar_settings and "KSPACING" in prev_incar: - # Prefer to inherit KSPACING from previous INCAR if it exists - # TODO: Is that we actually want to do? Copied from current pymatgen inputsets - incar["KSPACING"] = prev_incar["KSPACING"] - - # Ensure adequate number of KPOINTS are present for the tetrahedron method - # (ISMEAR=-5). If KSPACING is in the INCAR file the number of kpoints is not - # known before calling VASP, but a warning is raised when the KSPACING value is - # > 0.5 (2 reciprocal Angstrom). An error handler in Custodian is available to - # correct overly large KSPACING values (small number of kpoints) if necessary. - if kpoints is not None and np.prod(kpoints.kpts) < 4 and incar.get("ISMEAR", 0) == -5: - incar["ISMEAR"] = 0 - - if incar.get("KSPACING", 0) > 0.5 and incar.get("ISMEAR", 0) == -5: - warnings.warn( - "Large KSPACING value detected with ISMEAR = -5. Ensure that VASP " - "generates an adequate number of KPOINTS, lower KSPACING, or " - "set ISMEAR = 0", - BadInputSetWarning, - ) - - ismear = incar.get("ISMEAR", 1) - sigma = incar.get("SIGMA", 0.2) - if ( - all(elem.is_metal for elem in structure.composition) - and incar.get("NSW", 0) > 0 - and (ismear < 0 or (ismear == 0 and sigma > 0.05)) - ): - msg = "" - if ismear < 0: - msg = f"Relaxation of likely metal with ISMEAR < 0 ({ismear})." - elif ismear == 0 and sigma > 0.05: - msg = f"ISMEAR = 0 with a small SIGMA ({sigma}) detected." - warnings.warn( - f"{msg} See VASP recommendations on ISMEAR for metals (https://www.vasp.at/wiki/index.php/ISMEAR).", - BadInputSetWarning, - stacklevel=1, - ) - - return incar - - @property - def poscar(self) -> Poscar: - """Poscar.""" - if self.structure is None: - raise RuntimeError("No structure is associated with the input set!") - - site_properties = self.structure.site_properties - return Poscar( - self.structure, - velocities=site_properties.get("velocities"), - predictor_corrector=site_properties.get("predictor_corrector"), - predictor_corrector_preamble=self.structure.properties.get("predictor_corrector_preamble"), - lattice_velocities=self.structure.properties.get("lattice_velocities"), - ) - - @property - def potcar_functional(self) -> UserPotcarFunctional: - """The functional used for POTCAR generation.""" - return self.user_potcar_functional - - @property - def nelect(self) -> float: - """The default number of electrons for a given structure.""" - if self.structure is None: - raise RuntimeError("No structure is associated with the input set!") - - n_electrons_by_element = {p.element: p.nelectrons for p in self.potcar} - n_elect = sum( - num_atoms * n_electrons_by_element[el.symbol] for el, num_atoms in self.structure.composition.items() - ) - - return n_elect - (self.structure.charge if self.use_structure_charge else 0) - - @property - def kpoints(self) -> Kpoints | None: - """The KPOINTS file.""" - if self.structure is None: - raise RuntimeError("No structure is associated with the input set!") - - if ( - self.user_incar_settings.get("KSPACING") is not None - or self.incar_updates.get("KSPACING") is not None - or self._config_dict["INCAR"].get("KSPACING") is not None - ) and self.user_kpoints_settings == {}: - # If KSPACING specified then always use this over k-points - return None - - # Use user setting if set otherwise default to base config settings - kpoints_updates = self.kpoints_updates - if self.user_kpoints_settings != {}: - kconfig = deepcopy(self.user_kpoints_settings) - elif isinstance(kpoints_updates, Kpoints): - return kpoints_updates - elif kpoints_updates != {}: - kconfig = kpoints_updates - else: - kconfig = deepcopy(self._config_dict.get("KPOINTS", {})) - - if isinstance(kconfig, Kpoints): - return kconfig - - explicit = ( - kconfig.get("explicit") - or len(kconfig.get("added_kpoints", [])) > 0 - or "zero_weighted_reciprocal_density" in kconfig - or "zero_weighted_line_density" in kconfig - ) - # Handle length generation first as this doesn't support any additional options - if kconfig.get("length"): - if explicit: - raise ValueError( - "length option cannot be used with explicit k-point generation, " - "added_kpoints, or zero weighted k-points." - ) - # If length is in kpoints settings use Kpoints.automatic - return Kpoints.automatic(kconfig["length"]) - - base_kpoints = None - if kconfig.get("line_density"): - # Handle line density generation - kpath = HighSymmKpath(self.structure, **kconfig.get("kpath_kwargs", {})) - frac_k_points, k_points_labels = kpath.get_kpoints( - line_density=kconfig["line_density"], coords_are_cartesian=False - ) - base_kpoints = Kpoints( - comment="Non SCF run along symmetry lines", - style=Kpoints.supported_modes.Reciprocal, - num_kpts=len(frac_k_points), - kpts=frac_k_points, - labels=k_points_labels, - kpts_weights=[1] * len(frac_k_points), - ) - - elif kconfig.get("grid_density") or kconfig.get("reciprocal_density"): - # Handle regular weighted k-point grid generation - if kconfig.get("grid_density"): - base_kpoints = Kpoints.automatic_density(self.structure, int(kconfig["grid_density"]), self.force_gamma) - elif kconfig.get("reciprocal_density"): - density = kconfig["reciprocal_density"] - base_kpoints = Kpoints.automatic_density_by_vol(self.structure, density, self.force_gamma) - - if explicit and base_kpoints is not None: - sga = SpacegroupAnalyzer(self.structure, symprec=self.sym_prec) - mesh = sga.get_ir_reciprocal_mesh(base_kpoints.kpts[0]) - base_kpoints = Kpoints( - comment="Uniform grid", - style=Kpoints.supported_modes.Reciprocal, - num_kpts=len(mesh), - kpts=tuple(i[0] for i in mesh), - kpts_weights=[i[1] for i in mesh], - ) - else: - # If not explicit that means no other options have been specified - # so we can return the k-points as is - return base_kpoints - - zero_weighted_kpoints = None - if kconfig.get("zero_weighted_line_density"): - # zero_weighted k-points along line mode path - kpath = HighSymmKpath(self.structure) - frac_k_points, k_points_labels = kpath.get_kpoints( - line_density=kconfig["zero_weighted_line_density"], - coords_are_cartesian=False, - ) - zero_weighted_kpoints = Kpoints( - comment="Hybrid run along symmetry lines", - style=Kpoints.supported_modes.Reciprocal, - num_kpts=len(frac_k_points), - kpts=frac_k_points, - labels=k_points_labels, - kpts_weights=[0] * len(frac_k_points), - ) - elif kconfig.get("zero_weighted_reciprocal_density"): - zero_weighted_kpoints = Kpoints.automatic_density_by_vol( - self.structure, kconfig["zero_weighted_reciprocal_density"], self.force_gamma - ) - sga = SpacegroupAnalyzer(self.structure, symprec=self.sym_prec) - mesh = sga.get_ir_reciprocal_mesh(zero_weighted_kpoints.kpts[0]) - zero_weighted_kpoints = Kpoints( - comment="Uniform grid", - style=Kpoints.supported_modes.Reciprocal, - num_kpts=len(mesh), - kpts=tuple(i[0] for i in mesh), - kpts_weights=[0 for _ in mesh], - ) - - added_kpoints = None - if kconfig.get("added_kpoints"): - points: list = kconfig.get("added_kpoints", []) - added_kpoints = Kpoints( - comment="Specified k-points only", - style=Kpoints.supported_modes.Reciprocal, - num_kpts=len(points), - kpts=points, - labels=["user-defined"] * len(points), - kpts_weights=[0] * len(points), - ) - - if base_kpoints and not (added_kpoints or zero_weighted_kpoints): - return base_kpoints - if added_kpoints and not (base_kpoints or zero_weighted_kpoints): - return added_kpoints - - # Sanity check - if "line_density" in kconfig and zero_weighted_kpoints: - raise ValueError("Cannot combine line_density and zero weighted k-points options") - if zero_weighted_kpoints and not base_kpoints: - raise ValueError("Zero weighted k-points must be used with reciprocal_density or grid_density options") - if not (base_kpoints or zero_weighted_kpoints or added_kpoints): - raise ValueError( - "Invalid k-point generation algo. Supported Keys are 'grid_density' " - "for Kpoints.automatic_density generation, 'reciprocal_density' for " - "KPoints.automatic_density_by_vol generation, 'length' for " - "Kpoints.automatic generation, 'line_density' for line mode generation," - " 'added_kpoints' for specific k-points to include, " - " 'zero_weighted_reciprocal_density' for a zero weighted uniform mesh," - " or 'zero_weighted_line_density' for a zero weighted line mode mesh." - ) - - return _combine_kpoints(base_kpoints, zero_weighted_kpoints, added_kpoints) - - @property - def potcar(self) -> Potcar: - """The input set's POTCAR.""" - if self.structure is None: - raise RuntimeError("No structure is associated with the input set!") - - user_potcar_functional = self.user_potcar_functional - potcar = Potcar(self.potcar_symbols, functional=user_potcar_functional) - - # Warn if the selected POTCARs do not correspond to the chosen user_potcar_functional - for p_single in potcar: - if user_potcar_functional not in p_single.identify_potcar()[0]: - warnings.warn( - f"POTCAR data with symbol {p_single.symbol} is not known by pymatgen to " - f"correspond with the selected {user_potcar_functional=}. This POTCAR " - f"is known to correspond with functionals {p_single.identify_potcar(mode='data')[0]}. " - "Please verify that you are using the right POTCARs!", - BadInputSetWarning, - ) - - return potcar - - @property - def potcar_symbols(self) -> list[str]: - """List of POTCAR symbols.""" - elements = self.poscar.site_symbols - potcar_symbols = [] - settings = self._config_dict["POTCAR"] - - if isinstance(settings[elements[-1]], dict): - for el in elements: - potcar_symbols.append(settings[el]["symbol"] if el in settings else el) - else: - for el in elements: - potcar_symbols.append(settings.get(el, el)) - - return potcar_symbols - - def estimate_nbands(self) -> int: - """Estimate the number of bands that VASP will initialize a - calculation with by default. Note that in practice this - can depend on # of cores (if not set explicitly). - - Note that this formula is slightly different than the formula on the VASP wiki - (as of July 2023). This is because the formula in the source code (`main.F`) is - slightly different than what is on the wiki. - """ - if self.structure is None: - raise RuntimeError("No structure is associated with the input set!") - - n_ions = len(self.structure) - - # As per the VASP source, if non-spin polarized ignore n_mag - if self.incar["ISPIN"] == 1: - n_mag = 0 - - # Otherwise set equal to sum of total magmoms - else: - n_mag = sum(self.incar["MAGMOM"]) - n_mag = np.floor((n_mag + 1) / 2) - - possible_val_1 = np.floor((self.nelect + 2) / 2) + max(np.floor(n_ions / 2), 3) - possible_val_2 = np.floor(self.nelect * 0.6) - - n_bands = max(possible_val_1, possible_val_2) + n_mag - - if self.incar.get("LNONCOLLINEAR") is True: - n_bands = n_bands * 2 - - if n_par := self.incar.get("NPAR"): - n_bands = (np.floor((n_bands + n_par - 1) / n_par)) * n_par - - return int(n_bands) - - def override_from_prev_calc(self, prev_calc_dir: PathLike = ".") -> Self: - """Update the input set to include settings from a previous calculation. - - Args: - prev_calc_dir (PathLike): The path to the previous calculation directory. - - Returns: - VaspInputSet: A new input set with settings (Structure, k-points, incar, etc) - updated using the previous VASP run. - """ - self._set_previous(prev_calc_dir) - - 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." - ) - - files_to_transfer = {} - if getattr(self, "copy_chgcar", False): - chgcars = sorted(glob(str(Path(prev_calc_dir) / "CHGCAR*"))) - if chgcars: - files_to_transfer["CHGCAR"] = str(chgcars[-1]) - - if getattr(self, "copy_wavecar", False): - for fname in ("WAVECAR", "WAVEDER", "WFULL"): - if wavecar_files := sorted(glob(str(Path(prev_calc_dir) / (f"{fname}*")))): - if fname == "WFULL": - for wavecar_file in wavecar_files: - fname = Path(wavecar_file).name - fname = fname.split(".")[0] - files_to_transfer[fname] = wavecar_file - else: - files_to_transfer[fname] = str(wavecar_files[-1]) - - self.files_to_transfer.update(files_to_transfer) - return self - - @classmethod - def from_prev_calc(cls, prev_calc_dir: PathLike, **kwargs) -> Self: - """Generate a set of VASP input files for static calculations from a - directory of previous VASP run. - - Args: - prev_calc_dir (PathLike): Directory containing the outputs( - vasprun.xml and OUTCAR) of previous VASP run. - **kwargs: All kwargs supported by MPStaticSet, other than prev_incar - and prev_structure and prev_kpoints which are determined from - the prev_calc_dir. - """ - input_set = cls(_dummy_structure, **kwargs) - return input_set.override_from_prev_calc(prev_calc_dir=prev_calc_dir) - - def calculate_ng( - self, - max_prime_factor: int = 7, - must_inc_2: bool = True, - custom_encut: float | None = None, - custom_prec: str | None = None, - ) -> tuple: - """ - Calculate the NGX, NGY, and NGZ values using the information available in the INCAR and POTCAR - This is meant to help with making initial guess for the FFT grid so we can interact with the Charge density API. - - Args: - max_prime_factor (int): the valid prime factors of the grid size in each direction - VASP has many different setting for this to handle many compiling options. - For typical MPI options all prime factors up to 7 are allowed - must_inc_2 (bool): Whether 2 must be a prime factor of the result. Defaults to True. - custom_encut (float | None): Calculates the FFT grid parameters using a custom - ENCUT that may be different from what is generated by the input set. Defaults to None. - Do *not* use this unless you know what you are doing. - custom_prec (str | None): Calculates the FFT grid parameters using a custom prec - that may be different from what is generated by the input set. Defaults to None. - Do *not* use this unless you know what you are doing. - """ - # TODO throw error for Ultrasoft potentials - - _RYTOEV = 13.605826 - _AUTOA = 0.529177249 - - # TODO Only do this for VASP 6 for now. Older version require more advanced logic - - if custom_encut is not None: - encut = custom_encut - elif self.incar.get("ENCUT", 0) > 0: - encut = self.incar["ENCUT"] # get the ENCUT val - else: - encut = max(i_species.enmax for i_species in self.get_vasp_input()["POTCAR"]) - - # PREC=Normal is VASP default - PREC = self.incar.get("PREC", "Normal") if custom_prec is None else custom_prec - - # Check for unsupported / invalid PREC tags - if PREC[0].lower() in {"l", "m", "h"}: - raise NotImplementedError( - "PREC = LOW/MEDIUM/HIGH from VASP 4.x and not supported, Please use NORMA/SINGLE/ACCURATE" - ) - if PREC[0].lower() not in {"a", "s", "n", "l", "m", "h"}: - raise ValueError(f"{PREC=} does not exist. If this is no longer correct, please update this code.") - - CUTOFF = [ - np.sqrt(encut / _RYTOEV) / (2 * np.pi / (anorm / _AUTOA)) for anorm in self.poscar.structure.lattice.abc - ] - - # TODO This only works in VASP 6.x - _WFACT = 4 if PREC[0].lower() in {"a", "s"} else 3 - - def next_g_size(cur_g_size): - g_size = int(_WFACT * cur_g_size + 0.5) - return next_num_with_prime_factors(g_size, max_prime_factor, must_inc_2) - - ng_vec = [*map(next_g_size, CUTOFF)] - - # TODO This works for VASP 5.x and 6.x - finer_g_scale = 2 if PREC[0].lower() in {"a", "n"} else 1 - - return ng_vec, [ng_ * finer_g_scale for ng_ in ng_vec] - - @staticmethod - def from_directory(directory: PathLike, optional_files: dict | None = None) -> VaspInput: - """Load a set of VASP inputs from a directory. - - Note that only the standard INCAR, POSCAR, POTCAR and KPOINTS files are read - unless optional_filenames is specified. - - Args: - directory: Directory to read VASP inputs from. - optional_files: Optional files to read in as well as a dict of {filename: Object class}. - Object class must have a static/class method from_file. - """ - directory = Path(directory) - objs = {"INCAR": Incar, "KPOINTS": Kpoints, "POSCAR": Poscar, "POTCAR": Potcar} - - inputs = {} - for name, obj in objs.items(): - if (directory / name).exists(): - inputs[name.upper()] = obj.from_file(directory / name) # type: ignore[attr-defined] - else: - # Handle the case where there is no KPOINTS file - inputs[name.upper()] = None - - optional_inputs = {} - if optional_files is not None: - for name, obj in optional_files.items(): - optional_inputs[str(name)] = obj.from_file(directory / name) # type: ignore[attr-defined] - - return VaspInput( - incar=inputs["INCAR"], - kpoints=inputs["KPOINTS"], - poscar=inputs["POSCAR"], - potcar=inputs["POTCAR"], - optional_files=optional_inputs, # type: ignore[arg-type] - ) - - def _get_nedos(self, dedos: float) -> int: - """Automatic setting of NEDOS using the energy range and the energy step.""" - if self.prev_vasprun is None or self.prev_vasprun.eigenvalues is None: - return 2000 - - emax = max(eigs.max() for eigs in self.prev_vasprun.eigenvalues.values()) - emin = min(eigs.min() for eigs in self.prev_vasprun.eigenvalues.values()) - return int((emax - emin) / dedos) - - -# Create VaspInputGenerator alias to follow atomate2 terminology -VaspInputGenerator = VaspInputSet - - -@deprecated(VaspInputSet, deadline=(2025, 12, 31)) -class DictSet(VaspInputSet): - """Alias for VaspInputSet.""" - - -# Helper functions to determine valid FFT grids for VASP -def next_num_with_prime_factors(n: int, max_prime_factor: int, must_inc_2: bool = True) -> int: - """Get the next number greater than or equal to n that only has the desired prime factors. - - Args: - n (int): Initial guess at the grid density - max_prime_factor (int): the maximum prime factor - must_inc_2 (bool): 2 must be a prime factor of the result - - Returns: - int: first product of the prime_factors that is >= n - """ - if max_prime_factor < 2: - raise ValueError("Must choose a maximum prime factor greater than 2") - - prime_factors = primes_less_than(max_prime_factor) - for new_val in itertools.count(start=n): - if must_inc_2 and new_val % 2 != 0: - continue - - cur_val_ = new_val - for j in prime_factors: - while cur_val_ % j == 0: - cur_val_ //= j - if cur_val_ == 1: - return new_val - - raise ValueError("No factorable number found, not possible.") - - -def primes_less_than(max_val: int) -> list[int]: - """Get the primes less than or equal to the max value.""" - res = [] - for i in range(2, max_val + 1): - for j in range(2, i): - if i % j == 0: - break - else: - res.append(i) - 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: - 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. - - 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) -> None: - 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(f"{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: - """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 |= {"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 isinstance(self.prev_kpoints, Kpoints) - and self.prev_kpoints.style == Kpoints.supported_modes.Monkhorst - and not self.lepsilon - and self.structure is not None - ): - kpoints = Kpoints.automatic_density_by_vol( - self.structure, - int(self.reciprocal_density * factor), - self.force_gamma, - ) - k_div = cast(Kpoint, tuple(kp + 1 if kp % 2 == 1 else kp for kp in kpoints.kpts[0])) - return Kpoints.monkhorst_automatic(k_div) - - 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: tuple[str, ...] | bool = ( # type: ignore[assignment] - "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) -> None: - """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[str, Any]: - """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 |= {"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[str, Any]: - """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[str, Any]: - """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 |= {"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 |= {"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 |= {"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) -> None: - """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[str, Any]: - """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 |= {"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 |= {"ISMEAR": 0, "SIGMA": sigma} - - if self.optics: - # LREAL not supported with LOPTICS = True; automatic NEDOS usually - # underestimates, so set it explicitly - updates |= {"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[str, Any]: - """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: Tuple3Ints = (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) -> None: - 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[str, Any]: - """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 |= {"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[str, Any]: - """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[misc, union-attr] - 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") - - assert VaspInputSet.structure is not None - VaspInputSet.structure.fset(self, structure) - - -@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[str, Any]: - """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" and self.structure is not None: - isotopes = {ist.split("-")[0]: ist for ist in self.isotopes} - quad_efg = [ - float(Species(sp.name).get_nmr_quadrupole_moment(isotopes.get(sp.name))) - for sp in self.structure.species - ] - 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[str, Any]: - """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[str, Any]: - """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) -> None: - """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[str, Any]: - """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[str, Any]: - """Updates to the INCAR config for this calculation type.""" - updates: dict[str, Any] = {} - 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 |= {"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 |= { - "ALGO": "GW0", - "NELM": 1, - "NOMEGA": 80, - "ENCUTGW": 250, - "EDIFF": None, - "LOPTICS": None, - "LPEAD": None, - } - - elif self.mode == "BSE": - # Default parameters for BSE calculation. - updates |= {"ALGO": "BSE", "ANTIRES": 0, "NBANDSO": 20, "NBANDSV": 20} - - if nbands: - updates["NBANDS"] = nbands - - return updates - - @classmethod - def from_prev_calc(cls, prev_calc_dir: PathLike, 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 (PathLike): 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[str, Any]: - """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 |= {"ISIF": 2, "LVTOT": True, "NELMIN": 8} - if self.set_mix: - updates |= {"AMIN": 0.01, "AMIX": 0.2, "BMIX": 0.001} - if self.auto_dipole and self.structure is not None: - weights = [struct.species.weight for struct in self.structure] - center_of_mass = np.average(self.structure.frac_coords, weights=weights, axis=0) - updates |= {"IDIPOL": 3, "LDIPOL": True, "DIPOL": center_of_mass} - return updates - - @property - def kpoints_updates(self) -> Kpoints: - """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 - assert self.structure is not None - 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=[cast(Kpoint, tuple(kpt_calc))], - ) - - def as_dict(self, verbosity: int = 2) -> dict[str, Any]: - """ - 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) -> Kpoints: - """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 # type: ignore[union-attr] - 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=[cast(Kpoint, tuple(kpt_calc))], - ) - - @property - def incar_updates(self) -> dict[str, Any]: - """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: list[Structure], unset_encut: bool = 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: - """Poscar for structure of first end point.""" - return Poscar(self.structures[0]) - - @property - def poscars(self) -> list[Poscar]: - """List of Poscars.""" - return [Poscar(struct) for struct in self.structures] - - @staticmethod - def _process_structures(structures: list[Structure]) -> list[Structure]: - """Remove any atoms jumping 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: PathLike, - make_dir_if_not_present: bool = True, - write_cif: bool = False, # type: ignore[override] - write_path_cif: bool = False, - write_endpoint_inputs: bool = False, # type: ignore[override] - ) -> None: - """ - NEB inputs has a special directory structure where inputs are in 00, - 01, 02, .... - - Args: - output_dir (PathLike): 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")) - assert self.kpoints is not None - 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")) - assert end_point_param.kpoints is not None - 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) -> dict[str, Any]: - """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: - """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[str, Any]: - """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 and self.structure is not None: - if Element("H") in self.structure.species: - updates |= {"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) -> 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[str, Any]: - """Updates to the INCAR config for this calculation type.""" - # NPT-AIMD default settings - assert self.structure is not None - updates = { - "ALGO": "Fast", - "ISIF": 3, - "LANGEVIN_GAMMA": [10] * self.structure.ntypesp, - "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)] - updates["ENCUT"] = max(enmax) * 1.5 - return updates - - @property - def kpoints_updates(self) -> Kpoints: - """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) -> None: - 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[str, Any]: - """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 # 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 - - # Latest 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) -> None: - 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[str, int]: - """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[str, Any]: - """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[arg-type] - elif self.address_basis_file is not None: - basis = Lobsterin.get_basis( - structure=self.structure, # type: ignore[arg-type] - 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[union-attr] - 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[arg-type] - - 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: PathLike, - parse_dos: bool = True, - parse_eigen: bool = True, -) -> tuple[Vasprun, Outcar]: - """Get a Vasprun and Outcar from a directory. - - Args: - path: Path to get the vasprun.xml and OUTCAR. - parse_dos: Whether to parse dos. Defaults to True. - parse_eigen: Whether to parse eigenvalue. Defaults to True. - - Returns: - Vasprun and Outcar files. - """ - path = Path(path) - vruns = list(glob(str(path / "vasprun.xml*"))) - outcars = list(glob(str(path / "OUTCAR*"))) - - if not vruns or not outcars: - raise ValueError(f"Unable to get vasprun.xml/OUTCAR from prev calculation in {path}") - - vsfile_fullpath = str(path / "vasprun.xml") - 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), - ) - - -def get_structure_from_prev_run(vasprun: Vasprun, outcar: Outcar | None = None) -> Structure: - """Process structure from previous run. - - Args: - vasprun (Vasprun): Vasprun that contains the final structure from previous run. - outcar (Outcar): Outcar that contains the magnetization info from previous run. - - Returns: - Structure: The magmom-decorated structure that can be passed to get VASP input files, e.g. - get_kpoints(). - """ - structure = vasprun.final_structure - - site_properties = {} - # MAGMOM - if vasprun.is_spin: - if outcar and outcar.magnetization: - site_properties["magmom"] = [i["tot"] for i in outcar.magnetization] - else: - site_properties["magmom"] = vasprun.parameters["MAGMOM"] - - # LDAU - if vasprun.parameters.get("LDAU", False): - for key in ("LDAUU", "LDAUJ", "LDAUL"): - vals = vasprun.incar[key] - m = {} - l_val = [] - s = 0 - for site in structure: - if site.specie.symbol not in m: - m[site.specie.symbol] = vals[s] - s += 1 - l_val.append(m[site.specie.symbol]) - if len(l_val) == len(structure): - site_properties |= {key.lower(): l_val} - else: - raise ValueError(f"length of list {l_val} not the same as structure") - - return structure.copy(site_properties=site_properties) - - -def standardize_structure( - structure: Structure, - sym_prec: float = 0.1, - international_monoclinic: bool = True, -) -> Structure: - """Get the symmetrically standardized structure. - - Args: - structure (Structure): The structure. - sym_prec (float): Tolerance for symmetry finding for standardization. - international_monoclinic (bool): Whether to use international - convention (vs Curtarolo) for monoclinic. Defaults True. - - Returns: - The symmetrized structure. - """ - sym_finder = SpacegroupAnalyzer(structure, symprec=sym_prec) - new_structure = sym_finder.get_primitive_standard_structure(international_monoclinic=international_monoclinic) - - # The primitive structure finding has had several bugs in the past - # defend through validation - vpa_old = structure.volume / len(structure) - vpa_new = new_structure.volume / len(new_structure) - - if abs(vpa_old - vpa_new) / vpa_old > 0.02: - raise ValueError(f"Standardizing cell failed! VPA old: {vpa_old}, VPA new: {vpa_new}") - - matcher = StructureMatcher() - if not matcher.fit(structure, new_structure): - raise ValueError("Standardizing cell failed! Old structure doesn't match new.") - - return new_structure - - -class BadInputSetWarning(UserWarning): - """Warning class for bad but legal VASP inputs.""" - - -def batch_write_input( - structures: Sequence[Structure], - vasp_input_set=MPRelaxSet, - output_dir: PathLike = ".", - make_dir_if_not_present: bool = True, - subfolder: Callable | None = None, - sanitize: bool = False, - include_cif: bool = False, - potcar_spec: bool = False, - zip_output: bool = False, - **kwargs, -): - """ - Batch write VASP input for a sequence of structures to - output_dir, following the format output_dir/{group}/{formula}_{number}. - - Args: - structures ([Structure]): Sequence of Structures. - vasp_input_set (VaspInputSet): VaspInputSet class that creates - VASP input files from structures. Note that a class should be - supplied. Defaults to MPRelaxSet. - output_dir (str): Directory to output files. Defaults to current - directory ".". - make_dir_if_not_present (bool): Create the directory if not present. - Defaults to True. - subfolder (callable): Function to create subdirectory name from - structure. Defaults to simply "formula_count". - sanitize (bool): Boolean indicating whether to sanitize the - structure before writing the VASP input files. Sanitized output - are generally easier for viewing and certain forms of analysis. - Defaults to False. - include_cif (bool): Whether to output a CIF as well. CIF files are - generally better supported in visualization programs. - potcar_spec (bool): Instead of writing the POTCAR, write a "POTCAR.spec". - This is intended to help sharing an input set with people who might - not have a license to specific Potcar files. Given a "POTCAR.spec", - the specific POTCAR file can be re-generated using pymatgen with the - "generate_potcar" function in the pymatgen CLI. - zip_output (bool): If True, output will be zipped into a file with the - same name as the InputSet (e.g., MPStaticSet.zip) - **kwargs: Additional kwargs are passed to the vasp_input_set class - in addition to structure. - """ - output_dir = Path(output_dir) - for idx, site in enumerate(structures): - formula = re.sub(r"\s+", "", site.formula) - if subfolder is not None: - subdir = subfolder(site) - d = output_dir / subdir - else: - d = output_dir / f"{formula}_{idx}" - if sanitize: - site = site.copy(sanitize=True) - v = vasp_input_set(site, **kwargs) - v.write_input( - str(d), - make_dir_if_not_present=make_dir_if_not_present, - include_cif=include_cif, - potcar_spec=potcar_spec, - zip_output=zip_output, - ) - - -_dummy_structure = Structure( - [1, 0, 0, 0, 1, 0, 0, 0, 1], - ["I"], - [[0, 0, 0]], - site_properties={"magmom": [[0, 0, 1]]}, -) - - -def get_valid_magmom_struct( - structure: Structure, - inplace: bool = True, - spin_mode: str = "auto", -) -> Structure: - """ - Make sure that the structure has valid magmoms based on the kind of calculation. - - Fill in missing Magmom values. - - Args: - structure: The input structure - inplace: True: edit magmoms of the input structure; False: return new structure - spin_mode: "scalar"/"vector"/"none"/"auto" only first letter (s/v/n) is needed. - dictates how the spin configuration will be determined. - - - auto: read the existing magmom values and decide - - scalar: use a single scalar value (for spin up/down) - - vector: use a vector value for spin-orbit systems - - none: Remove all the magmom information - - Returns: - New structure if inplace is False - """ - default_values = {"s": 1.0, "v": [1.0, 1.0, 1.0], "n": None} - if spin_mode[0].lower() == "a": - mode = "n" - for site in structure: - if "magmom" not in site.properties or site.properties["magmom"] is None: - pass - elif isinstance(site.properties["magmom"], (float, int)): - if mode == "v": - raise TypeError("Magmom type conflict") - mode = "s" - if isinstance(site.properties["magmom"], int): - site.properties["magmom"] = float(site.properties["magmom"]) - elif len(site.properties["magmom"]) == 3: - if mode == "s": - raise TypeError("Magmom type conflict") - mode = "v" - else: - raise TypeError("Unrecognized Magmom Value") - else: - mode = spin_mode[0].lower() - - ret_struct = structure if inplace else structure.copy() - for site in ret_struct: - if mode == "n": - if "magmom" in site.properties: - site.properties.pop("magmom") - elif "magmom" not in site.properties or site.properties["magmom"] is None: - site.properties["magmom"] = default_values[mode] - - 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: Tuple3Ints | 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) -> None: - """Validate settings.""" - super().__post_init__() - self.mode = self.mode.upper() - if self.mode not in type(self).SUPPORTED_MODES: - raise ValueError(f"{self.mode} not one of the support modes : {type(self).SUPPORTED_MODES}") - - @property - def kpoints_updates(self) -> dict[str, float]: - """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[str, Any]: - """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 |= {"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 = ( - cast(tuple[int, int, int], self.prev_vasprun.kpoints.kpts[0]) if self.nkred is None else self.nkred - ) - updates |= {"NKREDX": self.nkred[0], "NKREDY": self.nkred[1], "NKREDZ": self.nkred[2]} - - return updates - - -def _get_ispin(vasprun: Vasprun | None, outcar: Outcar | None) -> Literal[1, 2]: - """Get value of ISPIN depending on the magnetisation in the OUTCAR and vasprun.""" - if outcar is not None and outcar.magnetization is not None: - # Turn off spin when magmom for every site is smaller than 0.02. - site_magmom = np.array([i["tot"] for i in outcar.magnetization]) - return 2 if np.any(np.abs(site_magmom) > 0.02) else 1 - if vasprun is not None: - return 2 if vasprun.is_spin else 1 - return 2 - - -def _get_recommended_lreal(structure: Structure) -> Literal["Auto", False]: - """Get recommended LREAL flag based on the structure.""" - return "Auto" if structure.num_sites > 16 else False - - -def _combine_kpoints(*kpoints_objects: Kpoints | None) -> Kpoints: - """Combine multiple Kpoints objects.""" - _labels: list[list[str]] = [] - _kpoints: list[Sequence[Kpoint]] = [] - _weights = [] - - for kpoints_object in filter(None, kpoints_objects): # type: ignore[var-annotated] - if kpoints_object.style != Kpoints.supported_modes.Reciprocal: - raise ValueError("Can only combine kpoints with style=Kpoints.supported_modes.Reciprocal") - if kpoints_object.labels is None: - _labels.append([""] * len(kpoints_object.kpts)) - else: - _labels.append(kpoints_object.labels) - - _kpoints.append(kpoints_object.kpts) - _weights.append(kpoints_object.kpts_weights) - - labels = np.concatenate(_labels).tolist() - kpoints = np.concatenate(_kpoints).tolist() - weights = np.concatenate(_weights).tolist() - - return Kpoints( - comment="Combined k-points", - style=Kpoints.supported_modes.Reciprocal, - num_kpts=len(kpoints), - kpts=cast(Sequence[Kpoint], kpoints), - labels=labels, - kpts_weights=weights, - ) - - -def _apply_incar_updates(incar: dict | Incar, updates: dict[str, Any], skip: Sequence[str] = ()) -> None: - """ - Apply updates to an INCAR file. - - Args: - incar (Incar): An incar. - updates (dict): Updates to apply. - skip (list of str): Keys to skip. - """ - for k, v in updates.items(): - if k in skip: - continue - - if v is None: - incar.pop(k, None) - else: - incar[k] = v - - -def _remove_unused_incar_params(incar: dict | Incar, skip: Sequence[str] = ()) -> None: - """ - Remove INCAR parameters that are not actively used by VASP. - - Args: - incar (Incar): An incar. - skip (list of str): Keys to skip. - """ - # Turn off IBRION/ISIF/POTIM if NSW = 0 - if incar.get("NSW", 0) == 0: - opt_flags = ("EDIFFG", "IBRION", "ISIF", "POTIM") - for opt_flag in opt_flags: - if opt_flag not in skip: - incar.pop(opt_flag, None) - - # Remove MAGMOM if they aren't used - if incar.get("ISPIN", 1) == 1 and "MAGMOM" not in skip: - incar.pop("MAGMOM", None) - - # Turn off +U flags if +U is not even used - if incar.get("LDAU", False) is False: - ldau_flags = ("LDAUU", "LDAUJ", "LDAUL", "LDAUTYPE") - for ldau_flag in ldau_flags: - if ldau_flag not in skip: - incar.pop(ldau_flag, None) - - -def _get_nedos(vasprun: Vasprun | None, dedos: float) -> int: - """Get NEDOS using the energy range and the energy step, - defaults to 2000. - """ - if vasprun is None: - return 2000 - - if vasprun.eigenvalues is None: - raise RuntimeError("eigenvalues cannot be None.") - - emax = max(eigs.max() for eigs in vasprun.eigenvalues.values()) - emin = min(eigs.min() for eigs in vasprun.eigenvalues.values()) - return int((emax - emin) / dedos) - - -def auto_kspacing(bandgap: float | None, bandgap_tol: float) -> float: - """Set kspacing based on the bandgap.""" - if bandgap is None or bandgap <= bandgap_tol: # metallic - return 0.22 - - rmin = max(1.5, 25.22 - 2.87 * bandgap) # Eq. 25 - kspacing = 2 * np.pi * 1.0265 / (rmin - 1.0183) # Eq. 29 - - # Cap kspacing at a max of 0.44, per internal benchmarking - return min(kspacing, 0.44) diff --git a/src/pymatgen/io/vasp/sets/__init__.py b/src/pymatgen/io/vasp/sets/__init__.py new file mode 100644 index 00000000000..3ec9d5596c7 --- /dev/null +++ b/src/pymatgen/io/vasp/sets/__init__.py @@ -0,0 +1,44 @@ +"""Re-export all VASP input sets, for more convenient imports and to maintain backwards compatible imports +following the split up of sets.py into submodules in #3865. +""" + +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/src/pymatgen/io/vasp/sets/base.py b/src/pymatgen/io/vasp/sets/base.py new file mode 100644 index 00000000000..72664a422eb --- /dev/null +++ b/src/pymatgen/io/vasp/sets/base.py @@ -0,0 +1,1538 @@ +# ruff: noqa: PGH003 +""" +This module defines the VaspInputSet abstract base class and a concrete implementation for the parameters developed +and tested by the core team of pymatgen, including the Materials Virtual Lab, Materials Project and the MIT high +throughput project. The basic concept behind an input set is to specify a scheme to generate a consistent set of VASP +inputs from a structure without further user intervention. This ensures comparability across runs. + +Read the following carefully before implementing new input sets: + +1. 99% of what needs to be done can be done by specifying user_incar_settings to override some of the defaults of + various input sets. Unless there is an extremely good reason to add a new set, **do not** add one. e.g. if you want + to turn the Hubbard U off, just set "LDAU": False as a user_incar_setting. +2. All derivative input sets should inherit appropriate configurations (e.g., from MPRelaxSet), and more often than + not, VaspInputSet should be the superclass. Superclass delegation should be used where possible. In particular, + you are not supposed to implement your own as_dict or from_dict for derivative sets unless you know what you are + doing. Improper overriding the as_dict and from_dict protocols is the major cause of implementation headaches. If + you need an example, look at how the MPStaticSet is initialized. + +The above are recommendations. The following are **UNBREAKABLE** rules: + +1. All input sets must take in a structure, list of structures or None as the first argument. If None, the input set + should perform a stateless initialization and before any output can be written, a structure must be set. +2. user_incar_settings, user_kpoints_settings and user__settings are ABSOLUTE. Any new sets you implement + must obey this. If a user wants to override your settings, you assume he knows what he is doing. Do not + magically override user supplied settings. You can issue a warning if you think the user is wrong. +3. All input sets must save all supplied args and kwargs as instance variables. e.g. self.arg = arg and + self.kwargs = kwargs in the __init__. This ensures the as_dict and from_dict work correctly. +""" + +from __future__ import annotations + +import abc +import itertools +import re +import warnings +from collections.abc import Sequence +from copy import deepcopy +from dataclasses import dataclass, field +from glob import glob +from pathlib import Path +from typing import TYPE_CHECKING, Literal, Union, cast + +import numpy as np +from monty.dev import deprecated +from monty.json import MSONable +from monty.serialization import loadfn + +from pymatgen.analysis.structure_matcher import StructureMatcher +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.typing import Kpoint + +if TYPE_CHECKING: + from typing import Any + + from typing_extensions import Self + + +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): + config = loadfn(MODULE_DIR / (f"{fname}.yaml")) + if "PARENT" in config: + parent_config = _load_yaml_config(config["PARENT"]) + for k, v in parent_config.items(): + if k not in config: + config[k] = v + elif isinstance(v, dict): + v_new = config.get(k, {}) + v_new.update(v) + config[k] = v_new + return config + + +@dataclass +class VaspInputSet(InputGenerator, abc.ABC): + """ + Base class representing a set of VASP input parameters with a structure + supplied as init parameters and initialized from a dict of settings. + This allows arbitrary settings to be input. In general, + this is rarely used directly unless there is a source of settings in yaml + format (e.g., from a REST interface). It is typically used by other + VaspInputSets for initialization. + + Special consideration should be paid to the way the MAGMOM initialization + for the INCAR is done. The initialization differs depending on the type of + structure and the configuration settings. The order in which the magmom is + determined is as follows: + + 1. If the site is specified in user_incar_settings, use that setting. + 2. If the site itself has a magmom setting (i.e. site.properties["magmom"] = float), + that is used. This can be set with structure.add_site_property(). + 3. If the species of the site has a spin setting, that is used. This can be set + with structure.add_spin_by_element(). + 4. If the species itself has a particular setting in the config file, that + is used, e.g. Mn3+ may have a different magmom than Mn4+. + 5. Lastly, the element symbol itself is checked in the config file. If + there are no settings, a default value of 0.6 is used. + + 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. + config_dict (dict): The config dictionary to use. + files_to_transfer (dict): A dictionary of {filename: filepath}. This allows the + transfer of files from a previous calculation. + user_incar_settings (dict): User INCAR settings. This allows a user to override + INCAR settings, e.g. setting a different MAGMOM for various elements or + species. Note that in the new scheme, ediff_per_atom and hubbard_u are no + longer args. Instead, the CONFIG supports EDIFF_PER_ATOM and EDIFF keys. + The former scales with # of atoms, the latter does not. If both are present, + EDIFF is preferred. To force such settings, just supply + user_incar_settings={"EDIFF": 1e-5, "LDAU": False} for example. The keys + 'LDAUU', 'LDAUJ', 'LDAUL' are special cases since pymatgen defines different + values depending on what anions are present in the structure, so these keys + can be defined in one of two ways, e.g. either {"LDAUU":{"O":{"Fe":5}}} to + set LDAUU for Fe to 5 in an oxide, or {"LDAUU":{"Fe":5}} to set LDAUU to 5 + regardless of the input structure. If a None value is given, that key is + unset. For example, {"ENCUT": None} will remove ENCUT from the + incar settings. Finally, KSPACING is a special setting and can be set to + "auto" in which the KSPACING is set automatically based on the band gap. + user_kpoints_settings (dict or Kpoints): Allow user to override kpoints setting + by supplying a dict. e.g. {"reciprocal_density": 1000}. User can also + supply Kpoints object. + user_potcar_settings (dict): Allow user to override POTCARs. e.g. {"Gd": + "Gd_3"}. This is generally not recommended. + constrain_total_magmom (bool): Whether to constrain the total magmom (NUPDOWN in + INCAR) to be the sum of the expected MAGMOM for all species. + sort_structure (bool): Whether to sort the structure (using the default sort + order of electronegativity) before generating input files. Defaults to True, + the behavior you would want most of the time. This ensures that similar + atomic species are grouped together. + user_potcar_functional (str): Functional to use. Default (None) is to use the + functional in the config dictionary. Valid values: "PBE", "PBE_52", + "PBE_54", "LDA", "LDA_52", "LDA_54", "PW91", "LDA_US", "PW91_US". + force_gamma (bool): Force gamma centered kpoint generation. Default (False) is + to use the Automatic Density kpoint scheme, which will use the Gamma + centered generation scheme for hexagonal cells, and Monkhorst-Pack otherwise. + reduce_structure (None/str): Before generating the input files, generate the + reduced structure. Default (None), does not alter the structure. Valid + values: None, "niggli", "LLL". + vdw: Adds default parameters for van-der-Waals functionals supported by VASP to + INCAR. Supported functionals are: DFT-D2, undamped DFT-D3, DFT-D3 with + Becke-Jonson damping, Tkatchenko-Scheffler, Tkatchenko-Scheffler with + iterative Hirshfeld partitioning, MBD@rSC, dDsC, Dion's vdW-DF, DF2, optPBE, + optB88, optB86b and rVV10. + use_structure_charge (bool): If set to True, then the overall charge of the + structure (structure.charge) is used to set the NELECT variable in the + INCAR. Default is False. + standardize (float): Whether to standardize to a primitive standard cell. + Defaults to False. + sym_prec (float): Tolerance for symmetry finding. + international_monoclinic (bool): Whether to use international convention (vs + Curtarolo) for monoclinic. Defaults True. + validate_magmom (bool): Ensure that the missing magmom values are filled in with + the VASP default value of 1.0. + inherit_incar (bool): Whether to inherit INCAR settings from previous + calculation. This might be useful to port Custodian fixes to child jobs but + can also be dangerous e.g. when switching from GGA to meta-GGA or relax to + static jobs. Defaults to True. + auto_kspacing (bool): If true, determines the value of KSPACING from the bandgap + of a previous calculation. + auto_ismear (bool): If true, the values for ISMEAR and SIGMA will be set + automatically depending on the bandgap of the system. If the bandgap is not + known (e.g., there is no previous VASP directory) then ISMEAR=0 and + SIGMA=0.2; if the bandgap is zero (a metallic system) then ISMEAR=2 and + SIGMA=0.2; if the system is an insulator, then ISMEAR=-5 (tetrahedron + smearing). Note, this only works when generating the input set from a + previous VASP directory. + auto_ispin (bool) = False: + If generating input set from a previous calculation, this controls whether + to disable magnetisation (ISPIN = 1) if the absolute value of all magnetic + moments are less than 0.02. + auto_lreal (bool) = False: + If True, automatically use the VASP recommended LREAL based on cell size. + auto_metal_kpoints + If true and the system is metallic, try and use ``reciprocal_density_metal`` + instead of ``reciprocal_density`` for metallic systems. Note, this only works + if the bandgap is not None. + bandgap_tol (float): Tolerance for determining if a system is metallic when + KSPACING is set to "auto". If the bandgap is less than this value, the + system is considered metallic. Defaults to 1e-4 (eV). + bandgap (float): Used for determining KSPACING if KSPACING == "auto" or + ISMEAR if auto_ismear == True. Set automatically when using from_prev_calc. + prev_incar (str or dict): Previous INCAR used for setting parent INCAR when + inherit_incar == True. Set automatically when using from_prev_calc. + prev_kpoints (str or Kpoints): Previous Kpoints. Set automatically when using + from_prev_calc. + """ + + structure: Structure | None = None + config_dict: dict = field(default_factory=dict) + files_to_transfer: dict = field(default_factory=dict) + user_incar_settings: dict = field(default_factory=dict) + user_kpoints_settings: dict = field(default_factory=dict) + user_potcar_settings: dict = field(default_factory=dict) + constrain_total_magmom: bool = False + sort_structure: bool = True + user_potcar_functional: UserPotcarFunctional = None + force_gamma: bool = False + reduce_structure: Literal["niggli", "LLL"] | None = None + vdw: str | None = None + use_structure_charge: bool = False + standardize: bool = False + sym_prec: float = 0.1 + international_monoclinic: bool = True + validate_magmom: bool = True + inherit_incar: bool | list[str] = False + auto_kspacing: bool = False + auto_ismear: bool = False + auto_ispin: bool = False + auto_lreal: bool = False + auto_metal_kpoints: bool = False + bandgap_tol: float = 1e-4 + bandgap: float | None = None + prev_incar: str | dict | None = None + prev_kpoints: str | Kpoints | None = None + _valid_potcars: Sequence[str] | None = None + + def __post_init__(self): + """Perform validation.""" + user_potcar_functional = self.user_potcar_functional + if (valid_potcars := self._valid_potcars) and user_potcar_functional not in valid_potcars: + raise ValueError(f"Invalid {user_potcar_functional=}, must be one of {valid_potcars}") + + if hasattr(self, "CONFIG"): + self.config_dict = self.CONFIG + + self._config_dict = deepcopy(self.config_dict) + + # these have been left to stay consistent with previous API + self.user_incar_settings = self.user_incar_settings or {} + self.user_kpoints_settings = self.user_kpoints_settings or {} + + self.vdw = self.vdw.lower() if isinstance(self.vdw, str) else self.vdw + if self.user_incar_settings.get("KSPACING") and self.user_kpoints_settings: + # self.user_kpoints_settings will never be `None` because it is set to + # an empty dict if it is `None`. + warnings.warn( + "You have specified KSPACING and also supplied KPOINTS " + "settings. KSPACING only has effect when there is no " + "KPOINTS file. Since both settings were given, pymatgen" + "will generate a KPOINTS file and ignore KSPACING." + "Remove the `user_kpoints_settings` argument to enable KSPACING.", + BadInputSetWarning, + ) + + if self.vdw: + vdw_par = loadfn(MODULE_DIR / "vdW_parameters.yaml") + if vdw_param := vdw_par.get(self.vdw): + self._config_dict["INCAR"].update(vdw_param) + else: + raise KeyError( + f"Invalid or unsupported van-der-Waals functional. Supported functionals are {', '.join(vdw_par)}." + ) + # 'or' case reads the POTCAR_FUNCTIONAL from the .yaml + self.user_potcar_functional: UserPotcarFunctional = self.user_potcar_functional or self._config_dict.get( + "POTCAR_FUNCTIONAL", "PBE" + ) + + # warn if a user is overriding POTCAR_FUNCTIONAL + if self.user_potcar_functional != self._config_dict.get("POTCAR_FUNCTIONAL", "PBE"): + warnings.warn( + "Overriding the POTCAR functional is generally not recommended " + " as it significantly affect the results of calculations and " + "compatibility with other calculations done with the same " + "input set. Note that some POTCAR symbols specified in " + "the configuration file may not be available in the selected " + "functional.", + BadInputSetWarning, + ) + + if self.user_potcar_settings: + warnings.warn( + "Overriding POTCARs is generally not recommended as it " + "significantly affect the results of calculations and " + "compatibility with other calculations done with the same " + "input set. In many instances, it is better to write a " + "subclass of a desired input set and override the POTCAR in " + "the subclass to be explicit on the differences.", + BadInputSetWarning, + ) + for key, val in self.user_potcar_settings.items(): + self._config_dict["POTCAR"][key] = val + + if not isinstance(self.structure, Structure): + self._structure = None + else: + self.structure = self.structure + + if isinstance(self.prev_incar, (Path, str)): + self.prev_incar = Incar.from_file(self.prev_incar) + + if isinstance(self.prev_kpoints, (Path, str)): + self.prev_kpoints = Kpoints.from_file(self.prev_kpoints) + + self.prev_vasprun = None + self.prev_outcar = None + self._ispin = None + + @deprecated(message="get_vasp_input will be removed in a future version of pymatgen. Use get_input_set instead.") + def get_vasp_input(self, structure=None) -> VaspInput: + """Get a VaspInput object. + + Returns: + VaspInput. + """ + return self.get_input_set(structure=structure) + + def write_input( + self, + output_dir: str, + make_dir_if_not_present: bool = True, + include_cif: bool | str = False, + potcar_spec: bool = False, + zip_output: bool | str = False, + ) -> None: + """Write a set of VASP input to a directory. + + 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. + include_cif (bool): Whether to write a CIF file in the output + directory for easier opening by VESTA. + potcar_spec (bool): Instead of writing the POTCAR, write a "POTCAR.spec". + This is intended to help sharing an input set with people who might + not have a license to specific Potcar files. Given a "POTCAR.spec", + the specific POTCAR file can be re-generated using pymatgen with the + "generate_potcar" function in the pymatgen CLI. + zip_output (bool): If True, output will be zipped into a file with the + same name as the InputSet (e.g., MPStaticSet.zip). + """ + vasp_input = self.get_input_set(potcar_spec=potcar_spec) + + cif_name = None + if include_cif: + struct = vasp_input["POSCAR"].structure + cif_name = f"{output_dir}/{struct.formula.replace(' ', '')}.cif" + + vasp_input.write_input( + output_dir=output_dir, + make_dir_if_not_present=make_dir_if_not_present, + cif_name=cif_name, + zip_name=f"{type(self).__name__}.zip" if zip_output else None, + files_to_transfer=self.files_to_transfer, + ) + + def as_dict(self, verbosity=2): + """ + Args: + verbosity: Verbosity for generated dict. If 1, structure is + excluded. + + Returns: + dict: MSONable VaspInputSet representation. + """ + dct = MSONable.as_dict(self) + if verbosity == 1: + dct.pop("structure", None) + return dct + + @property # type: ignore + def structure(self) -> Structure: # noqa: F811 + """Structure.""" + return self._structure + + @structure.setter + def structure(self, structure: Structure | None) -> None: + if not hasattr(self, "_config_dict"): + self._structure = structure + return + + if isinstance(structure, SiteCollection): # could be Structure or Molecule + if self.user_potcar_functional == "PBE_54" and "W" in structure.symbol_set: + # when using 5.4 POTCARs, default Tungsten POTCAR to W_Sv but still allow user to override + self.user_potcar_settings = {"W": "W_sv", **(self.user_potcar_settings or {})} + if self.reduce_structure: + structure = structure.get_reduced_structure(self.reduce_structure) + if self.sort_structure: + structure = structure.get_sorted_structure() + if self.validate_magmom: + get_valid_magmom_struct(structure, spin_mode="auto") + + struct_has_Yb = any(specie.symbol == "Yb" for site in structure for specie in site.species) + potcar_settings = self._config_dict.get("POTCAR", {}) + if self.user_potcar_settings: + potcar_settings.update(self.user_potcar_settings) + uses_Yb_2_psp = potcar_settings.get("Yb", None) == "Yb_2" + if struct_has_Yb and uses_Yb_2_psp: + warnings.warn( + "The structure contains Ytterbium (Yb) and this InputSet uses the Yb_2 PSP.\n" + "Yb_2 is known to often give bad results since Yb has oxidation state 3+ in most compounds.\n" + "See https://github.com/materialsproject/pymatgen/issues/2968 for details.", + BadInputSetWarning, + ) + if self.standardize and self.sym_prec: + structure = standardize_structure( + structure, + sym_prec=self.sym_prec, + international_monoclinic=self.international_monoclinic, + ) + self._structure = structure + + def get_input_set( + self, + structure: Structure | None = None, + prev_dir: str | Path | None = None, + potcar_spec: bool = False, + ) -> VaspInput: + """Get a VASP input set. + + Note, if both ``structure`` and ``prev_dir`` are set, then the structure + specified will be preferred over the final structure from the last VASP run. + + Args: + structure (Structure): A structure. + prev_dir (str or Path): A previous directory to generate the input set from. + potcar_spec (bool): Instead of generating a Potcar object, use a list of + potcar symbols. This will be written as a "POTCAR.spec" file. This is + intended to help sharing an input set with people who might not have a + license to specific Potcar files. Given a "POTCAR.spec", the specific + POTCAR file can be re-generated using pymatgen with the + "generate_potcar" function in the pymatgen CLI. + + Returns: + VaspInput: A VASP input object. + """ + if structure is None and prev_dir is None and self.structure is None: + raise ValueError("Either structure or prev_dir must be set") + + self._set_previous(prev_dir) + + if structure is not None: + self.structure = structure + + return VaspInput( + incar=self.incar, + kpoints=self.kpoints, + poscar=self.poscar, + potcar="\n".join(self.potcar_symbols) if potcar_spec else self.potcar, + potcar_spec=potcar_spec, + ) + + @property + def incar_updates(self) -> dict: + """Updates to the INCAR config for this calculation type.""" + return {} + + @property + def kpoints_updates(self) -> dict | Kpoints: + """Updates to the kpoints configuration for this calculation type. + + Note, these updates will be ignored if the user has set user_kpoint_settings. + + Returns: + dict or Kpoints: A dictionary of updates to apply to the KPOINTS config + or a Kpoints object. + """ + return {} + + def _set_previous(self, prev_dir: str | Path | None = None): + """Load previous calculation outputs.""" + if prev_dir is not None: + vasprun, outcar = get_vasprun_outcar(prev_dir) + self.prev_vasprun = vasprun + self.prev_outcar = outcar + self.prev_incar = vasprun.incar + self.prev_kpoints = Kpoints.from_dict(vasprun.kpoints.as_dict()) + + if vasprun.efermi is None: + # VASP doesn't output efermi in vasprun if IBRION = 1 + vasprun.efermi = outcar.efermi + + bs = vasprun.get_band_structure(efermi="smart") + self.bandgap = 0 if bs.is_metal() else bs.get_band_gap()["energy"] + if self.auto_ispin: + # turn off spin when magmom for every site is smaller than 0.02. + self._ispin = _get_ispin(vasprun, outcar) + + self.structure = get_structure_from_prev_run(vasprun, outcar) + + @property + def incar(self) -> Incar: + """The INCAR.""" + if self.structure is None: + raise RuntimeError("No structure is associated with the input set!") + + prev_incar: dict[str, Any] = {} + if self.inherit_incar is True and self.prev_incar: + prev_incar = self.prev_incar # type: ignore + elif isinstance(self.inherit_incar, (list, tuple)) and self.prev_incar: + prev_incar = {k: self.prev_incar[k] for k in self.inherit_incar if k in self.prev_incar} # type: ignore + + incar_updates = self.incar_updates + settings = dict(self._config_dict["INCAR"]) + auto_updates = {} + if self.auto_ispin and (self._ispin is not None): + auto_updates["ISPIN"] = self._ispin + + # breaking change - order in which settings applied inconsistent with atomate2 + # apply updates from input set generator to SETTINGS + # _apply_incar_updates(settings, incar_updates) + + # apply user incar settings to SETTINGS not to INCAR + _apply_incar_updates(settings, self.user_incar_settings) + + # generate incar + structure = self.structure + comp = structure.composition + elements = sorted((el for el in comp.elements if comp[el] > 0), key=lambda e: e.X) + most_electro_neg = elements[-1].symbol + poscar = Poscar(structure) + hubbard_u = settings.get("LDAU", False) + incar = Incar() + + for key, setting in settings.items(): + if key == "MAGMOM": + mag = [] + for site in structure: + if uic_magmom := self.user_incar_settings.get("MAGMOM", {}).get(site.species_string): + mag.append(uic_magmom) + elif hasattr(site, "magmom"): + mag.append(site.magmom) + elif getattr(site.specie, "spin", None) is not None: + mag.append(site.specie.spin) + elif str(site.specie) in setting: + if site.specie.symbol == "Co" and setting[str(site.specie)] <= 1.0: + warnings.warn( + "Co without an oxidation state is initialized as low spin by default in Pymatgen. " + "If this default behavior is not desired, please set the spin on the magmom on the " + "site directly to ensure correct initialization." + ) + mag.append(setting.get(str(site.specie))) + else: + if site.specie.symbol == "Co": + warnings.warn( + "Co without an oxidation state is initialized as low spin by default in Pymatgen. " + "If this default behavior is not desired, please set the spin on the magmom on the " + "site directly to ensure correct initialization." + ) + mag.append(setting.get(site.specie.symbol, 0.6)) + incar[key] = mag + elif key in ("LDAUU", "LDAUJ", "LDAUL"): + if hubbard_u: + if hasattr(structure[0], key.lower()): + m = {site.specie.symbol: getattr(site, key.lower()) for site in structure} + incar[key] = [m[sym] for sym in poscar.site_symbols] + # lookup specific LDAU if specified for most_electroneg atom + elif most_electro_neg in setting and isinstance(setting[most_electro_neg], dict): + incar[key] = [setting[most_electro_neg].get(sym, 0) for sym in poscar.site_symbols] + # else, use fallback LDAU value if it exists + else: + incar[key] = [ + setting.get(sym, 0) if isinstance(setting.get(sym, 0), (float, int)) else 0 + for sym in poscar.site_symbols + ] + elif key.startswith("EDIFF") and key != "EDIFFG": + if "EDIFF" not in settings and key == "EDIFF_PER_ATOM": + incar["EDIFF"] = float(setting) * len(structure) + else: + incar["EDIFF"] = float(settings["EDIFF"]) + elif key == "KSPACING" and self.auto_kspacing: + # default to metal if no prev calc available + bandgap = 0 if self.bandgap is None else self.bandgap + incar[key] = auto_kspacing(bandgap, self.bandgap_tol) + else: + incar[key] = setting + has_u = hubbard_u and sum(incar["LDAUU"]) > 0 + if not has_u: + for key in list(incar): + if key.startswith("LDAU"): + del incar[key] + + # Modify LMAXMIX if you have d or f electrons present. Note that if the user + # explicitly sets LMAXMIX in settings it will override this logic. + # Previously, this was only set if Hubbard U was enabled as per the VASP manual + # but following an investigation it was determined that this would lead to a + # significant difference between SCF -> NonSCF even without Hubbard U enabled. + # Thanks to Andrew Rosen for investigating and reporting. + if "LMAXMIX" not in settings: + # contains f-electrons + if any(el.Z > 56 for el in structure.composition): + incar["LMAXMIX"] = 6 + # contains d-electrons + elif any(el.Z > 20 for el in structure.composition): + incar["LMAXMIX"] = 4 + + # Warn user about LASPH for +U, meta-GGAs, hybrids, and vdW-DF + if not incar.get("LASPH", False) and ( + incar.get("METAGGA") + or incar.get("LHFCALC", False) + or incar.get("LDAU", False) + or incar.get("LUSE_VDW", False) + ): + warnings.warn("LASPH = True should be set for +U, meta-GGAs, hybrids, and vdW-DFT", BadInputSetWarning) + + # apply previous incar settings, be careful not to override user_incar_settings + # or the settings updates from the specific input set implementations + # also skip LDAU/MAGMOM as structure may have changed. + skip = list(self.user_incar_settings) + list(incar_updates) + skip += ["MAGMOM", "NUPDOWN", "LDAUU", "LDAUL", "LDAUJ"] + _apply_incar_updates(incar, prev_incar, skip=skip) + + if self.constrain_total_magmom: + nupdown = sum(mag if abs(mag) > 0.6 else 0 for mag in incar["MAGMOM"]) + if abs(nupdown - round(nupdown)) > 1e-5: + warnings.warn( + "constrain_total_magmom was set to True, but the sum of MAGMOM " + "values is not an integer. NUPDOWN is meant to set the spin " + "multiplet and should typically be an integer. You are likely " + "better off changing the values of MAGMOM or simply setting " + "NUPDOWN directly in your INCAR settings.", + UserWarning, + stacklevel=1, + ) + auto_updates["NUPDOWN"] = nupdown + + if self.use_structure_charge: + auto_updates["NELECT"] = self.nelect + + # Check that ALGO is appropriate + if incar.get("LHFCALC", False) is True and incar.get("ALGO", "Normal") not in ["Normal", "All", "Damped"]: + warnings.warn( + "Hybrid functionals only support Algo = All, Damped, or Normal.", + BadInputSetWarning, + ) + + if self.auto_ismear: + if self.bandgap is None: + # don't know if we are a metal or insulator so set ISMEAR and SIGMA to + # be safe with the most general settings + auto_updates.update(ISMEAR=0, SIGMA=0.2) + elif self.bandgap <= self.bandgap_tol: + auto_updates.update(ISMEAR=2, SIGMA=0.2) # metal + else: + auto_updates.update(ISMEAR=-5, SIGMA=0.05) # insulator + + if self.auto_lreal: + auto_updates.update(LREAL=_get_recommended_lreal(structure)) + + # apply updates from auto options, careful not to override user_incar_settings + _apply_incar_updates(incar, auto_updates, skip=list(self.user_incar_settings)) + + # apply updates from input set generator to INCAR + _apply_incar_updates(incar, incar_updates, skip=list(self.user_incar_settings)) + + # Finally, re-apply `self.user_incar_settings` to make sure any accidentally + # overwritten settings are changed back to the intended values. + # skip dictionary parameters to avoid dictionaries appearing in the INCAR + _apply_incar_updates(incar, self.user_incar_settings, skip=["LDAUU", "LDAUJ", "LDAUL", "MAGMOM"]) + + # Remove unused INCAR parameters + _remove_unused_incar_params(incar, skip=list(self.user_incar_settings)) + + kpoints = self.kpoints + if kpoints is not None: + # unset KSPACING as we are using a KPOINTS file + incar.pop("KSPACING", None) + elif "KSPACING" in incar and "KSPACING" not in self.user_incar_settings and "KSPACING" in prev_incar: + # prefer to inherit KSPACING from previous INCAR if it exists + # TODO: Is that we actually want to do? Copied from current pymatgen inputsets + incar["KSPACING"] = prev_incar["KSPACING"] + + # Ensure adequate number of KPOINTS are present for the tetrahedron method + # (ISMEAR=-5). If KSPACING is in the INCAR file the number of kpoints is not + # known before calling VASP, but a warning is raised when the KSPACING value is + # > 0.5 (2 reciprocal Angstrom). An error handler in Custodian is available to + # correct overly large KSPACING values (small number of kpoints) if necessary. + if kpoints is not None and np.prod(kpoints.kpts) < 4 and incar.get("ISMEAR", 0) == -5: + incar["ISMEAR"] = 0 + + if incar.get("KSPACING", 0) > 0.5 and incar.get("ISMEAR", 0) == -5: + warnings.warn( + "Large KSPACING value detected with ISMEAR = -5. Ensure that VASP " + "generates an adequate number of KPOINTS, lower KSPACING, or " + "set ISMEAR = 0", + BadInputSetWarning, + ) + + ismear = incar.get("ISMEAR", 1) + sigma = incar.get("SIGMA", 0.2) + if ( + all(elem.is_metal for elem in structure.composition) + and incar.get("NSW", 0) > 0 + and (ismear < 0 or (ismear == 0 and sigma > 0.05)) + ): + msg = "" + if ismear < 0: + msg = f"Relaxation of likely metal with ISMEAR < 0 ({ismear})." + elif ismear == 0 and sigma > 0.05: + msg = f"ISMEAR = 0 with a small SIGMA ({sigma}) detected." + warnings.warn( + f"{msg} See VASP recommendations on ISMEAR for metals (https://www.vasp.at/wiki/index.php/ISMEAR).", + BadInputSetWarning, + stacklevel=1, + ) + + return incar + + @property + def poscar(self) -> Poscar: + """Poscar.""" + if self.structure is None: + raise RuntimeError("No structure is associated with the input set!") + site_properties = self.structure.site_properties + return Poscar( + self.structure, + velocities=site_properties.get("velocities"), + predictor_corrector=site_properties.get("predictor_corrector"), + predictor_corrector_preamble=self.structure.properties.get("predictor_corrector_preamble"), + lattice_velocities=self.structure.properties.get("lattice_velocities"), + ) + + @property + def potcar_functional(self) -> UserPotcarFunctional: + """The functional used for POTCAR generation.""" + return self.user_potcar_functional + + @property + def nelect(self) -> float: + """The default number of electrons for a given structure.""" + if self.structure is None: + raise RuntimeError("No structure is associated with the input set!") + + n_electrons_by_element = {p.element: p.nelectrons for p in self.potcar} + n_elect = sum( + num_atoms * n_electrons_by_element[el.symbol] for el, num_atoms in self.structure.composition.items() + ) + + return n_elect - (self.structure.charge if self.use_structure_charge else 0) + + @property + def kpoints(self) -> Kpoints | None: + """The kpoints file.""" + if self.structure is None: + raise RuntimeError("No structure is associated with the input set!") + + if ( + self.user_incar_settings.get("KSPACING", None) is not None + or self.incar_updates.get("KSPACING", None) is not None + or self._config_dict["INCAR"].get("KSPACING", None) is not None + ) and self.user_kpoints_settings == {}: + # If KSPACING specified then always use this over k-points + return None + + # use user setting if set otherwise default to base config settings + kpoints_updates = self.kpoints_updates + if self.user_kpoints_settings != {}: + kconfig = deepcopy(self.user_kpoints_settings) + elif isinstance(kpoints_updates, Kpoints): + return kpoints_updates + elif kpoints_updates != {}: + kconfig = kpoints_updates + else: + kconfig = deepcopy(self._config_dict.get("KPOINTS", {})) + + if isinstance(kconfig, Kpoints): + return kconfig + + explicit = ( + kconfig.get("explicit") + or len(kconfig.get("added_kpoints", [])) > 0 + or "zero_weighted_reciprocal_density" in kconfig + or "zero_weighted_line_density" in kconfig + ) + # handle length generation first as this doesn't support any additional options + if kconfig.get("length"): + if explicit: + raise ValueError( + "length option cannot be used with explicit k-point generation, " + "added_kpoints, or zero weighted k-points." + ) + # If length is in kpoints settings use Kpoints.automatic + return Kpoints.automatic(kconfig["length"]) + + base_kpoints = None + if kconfig.get("line_density"): + # handle line density generation + kpath = HighSymmKpath(self.structure, **kconfig.get("kpath_kwargs", {})) + frac_k_points, k_points_labels = kpath.get_kpoints( + line_density=kconfig["line_density"], coords_are_cartesian=False + ) + base_kpoints = Kpoints( + comment="Non SCF run along symmetry lines", + style=Kpoints.supported_modes.Reciprocal, + num_kpts=len(frac_k_points), + kpts=frac_k_points, + labels=k_points_labels, + kpts_weights=[1] * len(frac_k_points), + ) + elif kconfig.get("grid_density") or kconfig.get("reciprocal_density"): + # handle regular weighted k-point grid generation + if kconfig.get("grid_density"): + base_kpoints = Kpoints.automatic_density(self.structure, int(kconfig["grid_density"]), self.force_gamma) + elif kconfig.get("reciprocal_density"): + density = kconfig["reciprocal_density"] + base_kpoints = Kpoints.automatic_density_by_vol(self.structure, density, self.force_gamma) + if explicit: + sga = SpacegroupAnalyzer(self.structure, symprec=self.sym_prec) + mesh = sga.get_ir_reciprocal_mesh(base_kpoints.kpts[0]) # type: ignore + base_kpoints = Kpoints( + comment="Uniform grid", + style=Kpoints.supported_modes.Reciprocal, + num_kpts=len(mesh), + kpts=tuple(i[0] for i in mesh), + kpts_weights=[i[1] for i in mesh], + ) + else: + # if not explicit that means no other options have been specified + # so we can return the k-points as is + return base_kpoints + + zero_weighted_kpoints = None + if kconfig.get("zero_weighted_line_density"): + # zero_weighted k-points along line mode path + kpath = HighSymmKpath(self.structure) + frac_k_points, k_points_labels = kpath.get_kpoints( + line_density=kconfig["zero_weighted_line_density"], + coords_are_cartesian=False, + ) + zero_weighted_kpoints = Kpoints( + comment="Hybrid run along symmetry lines", + style=Kpoints.supported_modes.Reciprocal, + num_kpts=len(frac_k_points), + kpts=frac_k_points, + labels=k_points_labels, + kpts_weights=[0] * len(frac_k_points), + ) + elif kconfig.get("zero_weighted_reciprocal_density"): + zero_weighted_kpoints = Kpoints.automatic_density_by_vol( + self.structure, kconfig["zero_weighted_reciprocal_density"], self.force_gamma + ) + sga = SpacegroupAnalyzer(self.structure, symprec=self.sym_prec) + mesh = sga.get_ir_reciprocal_mesh(zero_weighted_kpoints.kpts[0]) # type: ignore[arg-type] + zero_weighted_kpoints = Kpoints( + comment="Uniform grid", + style=Kpoints.supported_modes.Reciprocal, + num_kpts=len(mesh), + kpts=tuple(i[0] for i in mesh), + kpts_weights=[0 for _ in mesh], + ) + + added_kpoints = None + if kconfig.get("added_kpoints"): + points: list = kconfig.get("added_kpoints") # type: ignore + added_kpoints = Kpoints( + comment="Specified k-points only", + style=Kpoints.supported_modes.Reciprocal, + num_kpts=len(points), + kpts=points, + labels=["user-defined"] * len(points), + kpts_weights=[0] * len(points), + ) + + if base_kpoints and not (added_kpoints or zero_weighted_kpoints): + return base_kpoints + if added_kpoints and not (base_kpoints or zero_weighted_kpoints): + return added_kpoints + + # do some sanity checking + if "line_density" in kconfig and zero_weighted_kpoints: + raise ValueError("Cannot combine line_density and zero weighted k-points options") + if zero_weighted_kpoints and not base_kpoints: + raise ValueError("Zero weighted k-points must be used with reciprocal_density or grid_density options") + if not (base_kpoints or zero_weighted_kpoints or added_kpoints): + raise ValueError( + "Invalid k-point generation algo. Supported Keys are 'grid_density' " + "for Kpoints.automatic_density generation, 'reciprocal_density' for " + "KPoints.automatic_density_by_vol generation, 'length' for " + "Kpoints.automatic generation, 'line_density' for line mode generation," + " 'added_kpoints' for specific k-points to include, " + " 'zero_weighted_reciprocal_density' for a zero weighted uniform mesh," + " or 'zero_weighted_line_density' for a zero weighted line mode mesh." + ) + + return _combine_kpoints(base_kpoints, zero_weighted_kpoints, added_kpoints) # type: ignore + + @property + def potcar(self) -> Potcar: + """The input set's POTCAR.""" + if self.structure is None: + raise RuntimeError("No structure is associated with the input set!") + + user_potcar_functional = self.user_potcar_functional + potcar = Potcar(self.potcar_symbols, functional=user_potcar_functional) + + # warn if the selected POTCARs do not correspond to the chosen user_potcar_functional + for p_single in potcar: + if user_potcar_functional not in p_single.identify_potcar()[0]: + warnings.warn( + f"POTCAR data with symbol {p_single.symbol} is not known by pymatgen to " + f"correspond with the selected {user_potcar_functional=}. This POTCAR " + f"is known to correspond with functionals {p_single.identify_potcar(mode='data')[0]}. " + "Please verify that you are using the right POTCARs!", + BadInputSetWarning, + ) + + return potcar + + @property + def potcar_symbols(self): + """List of POTCAR symbols.""" + elements = self.poscar.site_symbols + potcar_symbols = [] + settings = self._config_dict["POTCAR"] + + if isinstance(settings[elements[-1]], dict): + for el in elements: + potcar_symbols.append(settings[el]["symbol"] if el in settings else el) + else: + for el in elements: + potcar_symbols.append(settings.get(el, el)) + + return potcar_symbols + + def estimate_nbands(self) -> int: + """Estimate the number of bands that VASP will initialize a + calculation with by default. Note that in practice this + can depend on # of cores (if not set explicitly). + Note that this formula is slightly different than the formula on the VASP wiki + (as of July 2023). This is because the formula in the source code (`main.F`) is + slightly different than what is on the wiki. + """ + if self.structure is None: + raise RuntimeError("No structure is associated with the input set!") + + n_ions = len(self.structure) + + if self.incar["ISPIN"] == 1: # per the VASP source, if non-spin polarized ignore n_mag + n_mag = 0 + else: # otherwise set equal to sum of total magmoms + n_mag = sum(self.incar["MAGMOM"]) + n_mag = np.floor((n_mag + 1) / 2) + + possible_val_1 = np.floor((self.nelect + 2) / 2) + max(np.floor(n_ions / 2), 3) + possible_val_2 = np.floor(self.nelect * 0.6) + + n_bands = max(possible_val_1, possible_val_2) + n_mag + + if self.incar.get("LNONCOLLINEAR") is True: + n_bands = n_bands * 2 + + if n_par := self.incar.get("NPAR"): + n_bands = (np.floor((n_bands + n_par - 1) / n_par)) * n_par + + return int(n_bands) + + def override_from_prev_calc(self, prev_calc_dir="."): + """ + Update the input set to include settings from a previous calculation. + + Args: + prev_calc_dir (str): The path to the previous calculation directory. + + Returns: + VaspInputSet: A new input set with settings (Structure, k-points, incar, etc) + updated using the previous VASP run. + """ + self._set_previous(prev_calc_dir) + + 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." + ) + + files_to_transfer = {} + if getattr(self, "copy_chgcar", False): + chgcars = sorted(glob(str(Path(prev_calc_dir) / "CHGCAR*"))) + if chgcars: + files_to_transfer["CHGCAR"] = str(chgcars[-1]) + + if getattr(self, "copy_wavecar", False): + for fname in ("WAVECAR", "WAVEDER", "WFULL"): + wavecar_files = sorted(glob(str(Path(prev_calc_dir) / (f"{fname}*")))) + if wavecar_files: + if fname == "WFULL": + for wavecar_file in wavecar_files: + fname = Path(wavecar_file).name + fname = fname.split(".")[0] + files_to_transfer[fname] = wavecar_file + else: + files_to_transfer[fname] = str(wavecar_files[-1]) + + self.files_to_transfer.update(files_to_transfer) + return self + + @classmethod + def from_prev_calc(cls, prev_calc_dir: str, **kwargs) -> Self: + """Generate a set of VASP input files for static calculations from a + directory of previous VASP run. + + Args: + prev_calc_dir (str): Directory containing the outputs( + vasprun.xml and OUTCAR) of previous vasp run. + **kwargs: All kwargs supported by MPStaticSet, other than prev_incar + and prev_structure and prev_kpoints which are determined from + the prev_calc_dir. + """ + input_set = cls(_dummy_structure, **kwargs) + return input_set.override_from_prev_calc(prev_calc_dir=prev_calc_dir) + + def __str__(self) -> str: + return type(self).__name__ + + def __repr__(self) -> str: + return type(self).__name__ + + def calculate_ng( + self, + max_prime_factor: int = 7, + must_inc_2: bool = True, + custom_encut: float | None = None, + custom_prec: str | None = None, + ) -> tuple: + """ + Calculates the NGX, NGY, and NGZ values using the information available in the INCAR and POTCAR + This is meant to help with making initial guess for the FFT grid so we can interact with the Charge density API. + + Args: + max_prime_factor (int): the valid prime factors of the grid size in each direction + VASP has many different setting for this to handle many compiling options. + For typical MPI options all prime factors up to 7 are allowed + must_inc_2 (bool): Whether 2 must be a prime factor of the result. Defaults to True. + custom_encut (float | None): Calculates the FFT grid parameters using a custom + ENCUT that may be different from what is generated by the input set. Defaults to None. + Do *not* use this unless you know what you are doing. + custom_prec (str | None): Calculates the FFT grid parameters using a custom prec + that may be different from what is generated by the input set. Defaults to None. + Do *not* use this unless you know what you are doing. + """ + # TODO throw error for Ultrasoft potentials + + _RYTOEV = 13.605826 + _AUTOA = 0.529177249 + + # TODO Only do this for VASP 6 for now. Older version require more advanced logic + + if custom_encut is not None: + encut = custom_encut + elif self.incar.get("ENCUT", 0) > 0: + encut = self.incar["ENCUT"] # get the ENCUT val + else: + encut = max(i_species.enmax for i_species in self.get_vasp_input()["POTCAR"]) + + # PREC=Normal is VASP default + PREC = self.incar.get("PREC", "Normal") if custom_prec is None else custom_prec + + # Check for unsupported / invalid PREC tags + if PREC[0].lower() in {"l", "m", "h"}: + raise NotImplementedError( + "PREC = LOW/MEDIUM/HIGH from VASP 4.x and not supported, Please use NORMA/SINGLE/ACCURATE" + ) + if PREC[0].lower() not in {"a", "s", "n", "l", "m", "h"}: + raise ValueError(f"{PREC=} does not exist. If this is no longer correct, please update this code.") + + CUTOFF = [ + np.sqrt(encut / _RYTOEV) / (2 * np.pi / (anorm / _AUTOA)) for anorm in self.poscar.structure.lattice.abc + ] + + # TODO This only works in VASP 6.x + _WFACT = 4 if PREC[0].lower() in {"a", "s"} else 3 + + def next_g_size(cur_g_size): + g_size = int(_WFACT * cur_g_size + 0.5) + return next_num_with_prime_factors(g_size, max_prime_factor, must_inc_2) + + ng_vec = [*map(next_g_size, CUTOFF)] + + # TODO This works for VASP 5.x and 6.x + finer_g_scale = 2 if PREC[0].lower() in {"a", "n"} else 1 + + return ng_vec, [ng_ * finer_g_scale for ng_ in ng_vec] + + @staticmethod + def from_directory(directory: str | Path, optional_files: dict | None = None) -> VaspInput: + """Load a set of VASP inputs from a directory. + + Note that only the standard INCAR, POSCAR, POTCAR and KPOINTS files are read + unless optional_filenames is specified. + + Parameters + ---------- + directory + Directory to read VASP inputs from. + optional_files + Optional files to read in as well as a dict of {filename: Object class}. + Object class must have a static/class method from_file. + """ + directory = Path(directory) + objs = {"INCAR": Incar, "KPOINTS": Kpoints, "POSCAR": Poscar, "POTCAR": Potcar} + + inputs = {} + for name, obj in objs.items(): + if (directory / name).exists(): + inputs[name.upper()] = obj.from_file(directory / name) # type: ignore[attr-defined] + else: + # handle the case where there is no KPOINTS file + inputs[name.upper()] = None + + optional_inputs = {} + if optional_files is not None: + for name, obj in optional_files.items(): + optional_inputs[str(name)] = obj.from_file(directory / name) # type: ignore[attr-defined] + + return VaspInput( + incar=inputs["INCAR"], + kpoints=inputs["KPOINTS"], + poscar=inputs["POSCAR"], + potcar=inputs["POTCAR"], + optional_files=optional_inputs, # type: ignore[arg-type] + ) + + def _get_nedos(self, dedos: float) -> int: + """Automatic setting of nedos using the energy range and the energy step.""" + if self.prev_vasprun is None: + return 2000 + + emax = max(eigs.max() for eigs in self.prev_vasprun.eigenvalues.values()) + emin = min(eigs.min() for eigs in self.prev_vasprun.eigenvalues.values()) + return int((emax - emin) / dedos) + + +# create VaspInputGenerator alias to follow atomate2 terminology +VaspInputGenerator = VaspInputSet + + +class DictSet(VaspInputSet): + """Alias for VaspInputSet.""" + + def __post_init__(self): + super().__post_init__() + warnings.warn( + "DictSet is deprecated, and will be removed on 2025-12-31. Use VaspInputSet", + category=FutureWarning, + stacklevel=2, + ) + + +# Helper functions to determine valid FFT grids for VASP +def next_num_with_prime_factors(n: int, max_prime_factor: int, must_inc_2: bool = True) -> int: + """Get the next number greater than or equal to n that only has the desired prime factors. + + Args: + n (int): Initial guess at the grid density + max_prime_factor (int): the maximum prime factor + must_inc_2 (bool): 2 must be a prime factor of the result + + Returns: + int: first product of the prime_factors that is >= n + """ + if max_prime_factor < 2: + raise ValueError("Must choose a maximum prime factor greater than 2") + prime_factors = primes_less_than(max_prime_factor) + for new_val in itertools.count(start=n): + if must_inc_2 and new_val % 2 != 0: + continue + cur_val_ = new_val + for j in prime_factors: + while cur_val_ % j == 0: + cur_val_ //= j + if cur_val_ == 1: + return new_val + raise ValueError("No factorable number found, not possible.") + + +def primes_less_than(max_val: int) -> list[int]: + """Get the primes less than or equal to the max value.""" + res = [] + for i in range(2, max_val + 1): + for j in range(2, i): + if i % j == 0: + break + else: + res.append(i) + return res + + +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. + + Args: + path: Path to get the vasprun.xml and OUTCAR. + parse_dos: Whether to parse dos. Defaults to True. + parse_eigen: Whether to parse eigenvalue. Defaults to True. + + Returns: + Vasprun and Outcar files. + """ + path = Path(path) + vruns = list(glob(str(path / "vasprun.xml*"))) + outcars = list(glob(str(path / "OUTCAR*"))) + + if len(vruns) == 0 or len(outcars) == 0: + raise ValueError(f"Unable to get vasprun.xml/OUTCAR from prev calculation in {path}") + vsfile_fullpath = str(path / "vasprun.xml") + 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) + + +def get_structure_from_prev_run(vasprun, outcar=None) -> Structure: + """ + Process structure from previous run. + + Args: + vasprun (Vasprun): Vasprun that contains the final structure from previous run. + outcar (Outcar): Outcar that contains the magnetization info from previous run. + + Returns: + Structure: The magmom-decorated structure that can be passed to get VASP input files, e.g. + get_kpoints(). + """ + structure = vasprun.final_structure + + site_properties = {} + # magmom + if vasprun.is_spin: + if outcar and outcar.magnetization: + site_properties["magmom"] = [i["tot"] for i in outcar.magnetization] + else: + site_properties["magmom"] = vasprun.parameters["MAGMOM"] + # LDAU + if vasprun.parameters.get("LDAU", False): + for key in ("LDAUU", "LDAUJ", "LDAUL"): + vals = vasprun.incar[key] + m = {} + l_val = [] + s = 0 + for site in structure: + if site.specie.symbol not in m: + m[site.specie.symbol] = vals[s] + s += 1 + l_val.append(m[site.specie.symbol]) + if len(l_val) == len(structure): + site_properties.update({key.lower(): l_val}) + else: + raise ValueError(f"length of list {l_val} not the same as structure") + + return structure.copy(site_properties=site_properties) + + +def standardize_structure(structure, sym_prec=0.1, international_monoclinic=True): + """Get the symmetrically standardized structure. + + Args: + structure (Structure): The structure. + sym_prec (float): Tolerance for symmetry finding for standardization. + international_monoclinic (bool): Whether to use international + convention (vs Curtarolo) for monoclinic. Defaults True. + + Returns: + The symmetrized structure. + """ + sym_finder = SpacegroupAnalyzer(structure, symprec=sym_prec) + new_structure = sym_finder.get_primitive_standard_structure(international_monoclinic=international_monoclinic) + + # the primitive structure finding has had several bugs in the past + # defend through validation + vpa_old = structure.volume / len(structure) + vpa_new = new_structure.volume / len(new_structure) + + if abs(vpa_old - vpa_new) / vpa_old > 0.02: + raise ValueError(f"Standardizing cell failed! VPA old: {vpa_old}, VPA new: {vpa_new}") + + matcher = StructureMatcher() + if not matcher.fit(structure, new_structure): + raise ValueError("Standardizing cell failed! Old structure doesn't match new.") + + return new_structure + + +class BadInputSetWarning(UserWarning): + """Warning class for bad but legal VASP inputs.""" + + +def batch_write_input( + structures, + vasp_input_set=None, + output_dir=".", + make_dir_if_not_present=True, + subfolder=None, + sanitize=False, + include_cif=False, + potcar_spec=False, + zip_output=False, + **kwargs, +): + """ + Batch write vasp input for a sequence of structures to + output_dir, following the format output_dir/{group}/{formula}_{number}. + + Args: + structures ([Structure]): Sequence of Structures. + vasp_input_set (VaspInputSet): VaspInputSet class that creates + vasp input files from structures. Note that a class should be + supplied. Defaults to MPRelaxSet. + output_dir (str): Directory to output files. Defaults to current + directory ".". + make_dir_if_not_present (bool): Create the directory if not present. + Defaults to True. + subfolder (callable): Function to create subdirectory name from + structure. Defaults to simply "formula_count". + sanitize (bool): Boolean indicating whether to sanitize the + structure before writing the VASP input files. Sanitized output + are generally easier for viewing and certain forms of analysis. + Defaults to False. + include_cif (bool): Whether to output a CIF as well. CIF files are + generally better supported in visualization programs. + potcar_spec (bool): Instead of writing the POTCAR, write a "POTCAR.spec". + This is intended to help sharing an input set with people who might + not have a license to specific Potcar files. Given a "POTCAR.spec", + the specific POTCAR file can be re-generated using pymatgen with the + "generate_potcar" function in the pymatgen CLI. + zip_output (bool): If True, output will be zipped into a file with the + same name as the InputSet (e.g., MPStaticSet.zip) + **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) + if subfolder is not None: + subdir = subfolder(site) + d = output_dir / subdir + else: + d = output_dir / f"{formula}_{idx}" + if sanitize: + site = site.copy(sanitize=True) + v = vasp_input_set(site, **kwargs) + v.write_input( + str(d), + make_dir_if_not_present=make_dir_if_not_present, + include_cif=include_cif, + potcar_spec=potcar_spec, + zip_output=zip_output, + ) + + +_dummy_structure = Structure( + [1, 0, 0, 0, 1, 0, 0, 0, 1], + ["I"], + [[0, 0, 0]], + site_properties={"magmom": [[0, 0, 1]]}, +) + + +def get_valid_magmom_struct(structure: Structure, inplace: bool = True, spin_mode: str = "auto") -> Structure: + """ + Make sure that the structure has valid magmoms based on the kind of calculation. + + Fill in missing Magmom values. + + Args: + structure: The input structure + inplace: True: edit magmoms of the input structure; False: return new structure + spin_mode: "scalar"/"vector"/"none"/"auto" only first letter (s/v/n) is needed. + dictates how the spin configuration will be determined. + + - auto: read the existing magmom values and decide + - scalar: use a single scalar value (for spin up/down) + - vector: use a vector value for spin-orbit systems + - none: Remove all the magmom information + + Returns: + New structure if inplace is False + """ + default_values = {"s": 1.0, "v": [1.0, 1.0, 1.0], "n": None} + if spin_mode[0].lower() == "a": + mode = "n" + for site in structure: + if "magmom" not in site.properties or site.properties["magmom"] is None: + pass + elif isinstance(site.properties["magmom"], (float, int)): + if mode == "v": + raise TypeError("Magmom type conflict") + mode = "s" + if isinstance(site.properties["magmom"], int): + site.properties["magmom"] = float(site.properties["magmom"]) + elif len(site.properties["magmom"]) == 3: + if mode == "s": + raise TypeError("Magmom type conflict") + mode = "v" + else: + raise TypeError("Unrecognized Magmom Value") + else: + mode = spin_mode[0].lower() + + ret_struct = structure if inplace else structure.copy() + for site in ret_struct: + if mode == "n": + if "magmom" in site.properties: + site.properties.pop("magmom") + elif "magmom" not in site.properties or site.properties["magmom"] is None: + site.properties["magmom"] = default_values[mode] + + return ret_struct + + +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: + # Turn off spin when magmom for every site is smaller than 0.02. + site_magmom = np.array([i["tot"] for i in outcar.magnetization]) + return 2 if np.any(np.abs(site_magmom) > 0.02) else 1 + if vasprun is not None: + return 2 if vasprun.is_spin else 1 + return 2 + + +def _get_recommended_lreal(structure: Structure) -> str | bool: + """Get recommended LREAL flag based on the structure.""" + return "Auto" if structure.num_sites > 16 else False + + +def _combine_kpoints(*kpoints_objects: Sequence[Kpoints]) -> Kpoints: + """Combine multiple Kpoints objects.""" + _labels: list[list[str]] = [] + _kpoints: list[Sequence[Kpoint]] = [] + _weights = [] + + kpoints_obj: Kpoints + for kpoints_obj in kpoints_objects: # type: ignore[assignment] + if kpoints_obj is None: + continue + if kpoints_obj.style != Kpoints.supported_modes.Reciprocal: + raise ValueError("Can only combine kpoints with style=Kpoints.supported_modes.Reciprocal") + if kpoints_obj.labels is None: + _labels.append([""] * len(kpoints_obj.kpts)) + else: + _labels.append(kpoints_obj.labels) + + _kpoints.append(kpoints_obj.kpts) + _weights.append(kpoints_obj.kpts_weights) + + labels = np.concatenate(_labels).tolist() + kpoints = np.concatenate(_kpoints).tolist() + weights = np.concatenate(_weights).tolist() + + return Kpoints( + comment="Combined k-points", + style=Kpoints.supported_modes.Reciprocal, + num_kpts=len(kpoints), + kpts=cast(Sequence[Kpoint], kpoints), + labels=labels, + kpts_weights=weights, + ) + + +def _apply_incar_updates(incar, updates, skip: Sequence[str] = ()) -> None: + """ + Apply updates to an INCAR file. + + Args: + incar (dict): An incar. + updates (dict): Updates to apply. + skip (list of str): Keys to skip. + """ + for k, v in updates.items(): + if k in skip: + continue + + if v is None: + incar.pop(k, None) + else: + incar[k] = v + + +def _remove_unused_incar_params(incar, skip: Sequence[str] = ()) -> None: + """ + Remove INCAR parameters that are not actively used by VASP. + + Args: + incar (dict): An incar. + skip (list of str): Keys to skip. + """ + # Turn off IBRION/ISIF/POTIM if NSW = 0 + opt_flags = ["EDIFFG", "IBRION", "ISIF", "POTIM"] + if incar.get("NSW", 0) == 0: + for opt_flag in opt_flags: + if opt_flag not in skip: + incar.pop(opt_flag, None) + + # Remove MAGMOMs if they aren't used + if incar.get("ISPIN", 1) == 1 and "MAGMOM" not in skip: + incar.pop("MAGMOM", None) + + # Turn off +U flags if +U is not even used + ldau_flags = ["LDAUU", "LDAUJ", "LDAUL", "LDAUTYPE"] + if incar.get("LDAU", False) is False: + for ldau_flag in ldau_flags: + if ldau_flag not in skip: + incar.pop(ldau_flag, None) + + +def _get_nedos(vasprun: Vasprun | None, dedos: float) -> int: + """Get NEDOS using the energy range and the energy step, + defaults to 2000. + """ + if vasprun is None: + return 2000 + + if vasprun.eigenvalues is None: + raise RuntimeError("eigenvalues cannot be None.") + + emax = max(eigs.max() for eigs in vasprun.eigenvalues.values()) + emin = min(eigs.min() for eigs in vasprun.eigenvalues.values()) + return int((emax - emin) / dedos) + + +def auto_kspacing(bandgap, bandgap_tol): + """Set kspacing based on the bandgap.""" + if bandgap is None or bandgap <= bandgap_tol: # metallic + return 0.22 + + rmin = max(1.5, 25.22 - 2.87 * bandgap) # Eq. 25 + kspacing = 2 * np.pi * 1.0265 / (rmin - 1.0183) # Eq. 29 + + # cap kspacing at a max of 0.44, per internal benchmarking + return min(kspacing, 0.44) diff --git a/src/pymatgen/io/vasp/sets/lobster.py b/src/pymatgen/io/vasp/sets/lobster.py new file mode 100644 index 00000000000..07c9d63ebe5 --- /dev/null +++ b/src/pymatgen/io/vasp/sets/lobster.py @@ -0,0 +1,107 @@ +# ruff: noqa: PGH003 + +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/src/pymatgen/io/vasp/sets/matpes.py b/src/pymatgen/io/vasp/sets/matpes.py new file mode 100644 index 00000000000..d78fe22888e --- /dev/null +++ b/src/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 # noqa: PGH003 + "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/src/pymatgen/io/vasp/sets/mit.py b/src/pymatgen/io/vasp/sets/mit.py new file mode 100644 index 00000000000..bd0426f351b --- /dev/null +++ b/src/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/src/pymatgen/io/vasp/sets/mp.py b/src/pymatgen/io/vasp/sets/mp.py new file mode 100644 index 00000000000..12cf6da757d --- /dev/null +++ b/src/pymatgen/io/vasp/sets/mp.py @@ -0,0 +1,862 @@ +# ruff: noqa: PGH003 + +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 Element, Species, Structure +from pymatgen.io.vasp.inputs import Kpoints +from pymatgen.io.vasp.sets.base import ( + MODULE_DIR, + 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 key in vdw_par[self.vdw]: + self._config_dict["INCAR"].pop(key, 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/src/pymatgen/io/vasp/sets/mvl.py b/src/pymatgen/io/vasp/sets/mvl.py new file mode 100644 index 00000000000..49e44dbcf35 --- /dev/null +++ b/src/pymatgen/io/vasp/sets/mvl.py @@ -0,0 +1,445 @@ +# ruff: noqa: PGH003 + +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[union-attr] + 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 diff --git a/tests/io/vasp/test_sets.py b/tests/io/vasp/test_sets.py index 738f994b0bd..4556d74209e 100644 --- a/tests/io/vasp/test_sets.py +++ b/tests/io/vasp/test_sets.py @@ -2121,7 +2121,7 @@ def test_vasp_input_set_alias(): def test_dict_set_alias(): with pytest.warns( - FutureWarning, match="DictSet is deprecated, and will be removed on 2025-12-31\nUse VaspInputSet" + FutureWarning, match="DictSet is deprecated, and will be removed on 2025-12-31. Use VaspInputSet" ): DictSet() assert isinstance(DictSet(), VaspInputSet)