Housekeeping (#1270)

* bump ase to 3.26 + add MTKNPT as default for NPT

* improve VASP documentation around run_vasp_kwargs, custodian settings, and magmom initialization

* remove deprecated forcefield makers that were slated for deletion on 1 January, 2025

* Suppress warnings in file client by default

* Fallback methods in MPMorph's automatic volume determination

Author: esoteric-ephemera

docs/user/codes/vasp.md

@@ -42,6 +42,66 @@ The most important settings to consider are:
   NCORE, KPAR etc.
 - `VASP_VDW_KERNEL_DIR`: The path to the VASP Van der Waals kernel.
 
+## FAQs
+
+<b>How can I update the Custodian handlers used in a VASP job?</b>
+- Every `Maker` which derives from `BaseVaspMaker` (see below) has a `run_vasp_kwargs` kwarg.
+So, for example, to run a `StaticMaker` with only the `VaspErrorHandler` as the custodian handler, you would do this:
+
+```py
+from atomate2.vasp.jobs.core import StaticMaker
+from custodian.vasp.handlers import VaspErrorHandler
+
+maker = StaticMaker(run_vasp_kwargs={"handlers": [VaspErrorHandler]})
+```
+
+<b>How can I change the other Custodian settings used to run VASP?</b>
+- These can be set through `run_vasp_kwargs.custodian_kwargs`:
+
+```py
+maker = StaticMaker(
+    run_vasp_kwargs={
+        "custodian_kwargs": {
+            "max_errors_per_job": 5,
+            "monitor_freq": 100,
+        }
+    }
+)
+```
+For all possible `custodian_kwargs`, see the [`Custodian` class](https://github.com/materialsproject/custodian/blob/aa02baf5bc2a1883c5f8a8b6808340eeae324a99/src/custodian/custodian.py#L68).
+NB: You cannot set the following four `Custodian` fields using `custodian_kwargs`: `handlers`, `jobs`, `validators`, `max_errors`, and `scratch_dir`.
+
+<b>Can I change how VASP is run for each per-`Maker`/`Job`?</b>
+Yes! Still using the `run_vasp_kwargs` and either the `vasp_cmd`, which represents the path to (and including) `vasp_std`, and `vasp_gamma_cmd`, which is the path to `vasp_gam`:
+
+```py
+maker = StaticMaker(run_vasp_kwargs={"vasp_cmd": "/path/to/vasp_std"})
+```
+
+<b>How can I use non-colinear VASP?</b>
+The same method as before applies, you can simply set:
+
+```py
+vasp_ncl_path = "/path/to/vasp_ncl"
+maker = StaticMaker(
+    run_vasp_kwargs={"vasp_cmd": vasp_ncl_path, "vasp_gamma_cmd": vasp_ncl_path}
+)
+```
+
+<b>How can I update the magnetic moments (MAGMOM) tag used to start a calculation?</b>
+You can specify MAGMOM using a `dict` of defined values, such as:
+
+```py
+from pymatgen.io.vasp.sets import MPRelaxSet
+
+vis = MPRelaxSet(user_incar_settings={"MAGMOM": {"In": 0.5, "Ga": 0.5, "As": -0.5}})
+```
+You can also specify different magnetic moments for different oxidation states, such as `{"Ga3+": 0.25}`.
+However, note that `"Ga0+"`, which has been assigned zero-valued oxidation state, is distinct from `"Ga"`, which has not been assigned an oxidation state.
+
+Alternatively, MAGMOM can be set by giving a structure assigned magnetic moments: `structure.add_site_property("magmom", list[float])`.
+This will override the default MAGMOM settings of a VASP input set.
+
 (vasp_workflows)=
 
 ## List of VASP workflows

pyproject.toml

@@ -51,7 +51,7 @@ defects = [
     "python-ulid>=2.7",
 ]
 forcefields = [
-    "ase>=3.25.0",
+    "ase>=3.26.0",
     "calorine>=3.0",
     "chgnet>=0.2.2",
     "mace-torch>=0.3.3",
@@ -63,7 +63,7 @@ forcefields = [
     "torchdata<=0.7.1",                            # TODO: remove when issue fixed
 ]
 approxneb = ["pymatgen-analysis-diffusion>=2024.7.15"]
-ase = ["ase>=3.25.0"]
+ase = ["ase>=3.26.0"]
 ase-ext = ["tblite>=0.3.0; platform_system=='Linux'"]
 openmm = [
     "mdanalysis>=2.8.0",
@@ -95,7 +95,7 @@ tests = [
 ]
 strict = [
     "PyYAML==6.0.2",
-    "ase==3.25.0",
+    "ase==3.26.0",
     "cclib==1.8.1",
     "click==8.2.1",
     "custodian==2025.8.13",

src/atomate2/ase/md.py

@@ -4,6 +4,7 @@
 
 import contextlib
 import io
+import logging
 import os
 import sys
 import time
@@ -41,6 +42,8 @@
 
     from atomate2.ase.schemas import AseMoleculeTaskDoc, AseStructureTaskDoc
 
+logger = logging.getLogger(__name__)
+
 
 class MDEnsemble(Enum):
     """Define known MD ensembles."""
@@ -60,12 +63,13 @@ class DynamicsPresets(Enum):
     nvt_nose_hoover = "ase.md.npt.NPT"
     npt_berendsen = "ase.md.nptberendsen.NPTBerendsen"
     npt_nose_hoover = "ase.md.npt.NPT"  # noqa: PIE796
+    npt_nose_hoover_chain = "ase.md.nose_hoover_chain.MTKNPT"
 
 
 default_dynamics = {
     MDEnsemble.nve: "velocityverlet",
     MDEnsemble.nvt: "langevin",
-    MDEnsemble.npt: "nose-hoover",
+    MDEnsemble.npt: "nose-hoover-chain",
 }
 
 _valid_dynamics: dict[MDEnsemble, set[str]] = {}
@@ -251,18 +255,23 @@ def _get_ensemble_defaults(self) -> None:
             if (
                 (
                     isinstance(self.dynamics, DynamicsPresets)
-                    and DynamicsPresets(self.dynamics) == DynamicsPresets.npt_berendsen
+                    and DynamicsPresets(self.dynamics)
+                    == DynamicsPresets.npt_nose_hoover
                 )
                 or (
                     isinstance(self.dynamics, type)
                     and issubclass(self.dynamics, MolecularDynamics)
-                    and self.dynamics.__name__ == "NPTBerendsen"
+                    and self.dynamics.__name__ == "NPT"
                 )
-                or (isinstance(self.dynamics, str) and self.dynamics == "berendsen")
+                or (isinstance(self.dynamics, str) and self.dynamics == "nose-hoover")
             ):
-                stress_kwarg = "pressure_au"
-            else:
+                logger.warning(
+                    "The `NPT` module in ASE is no longer recommended."
+                    "Users are advised to switch to Nose-Hoover chain / MTKNPT."
+                )
                 stress_kwarg = "externalstress"
+            else:
+                stress_kwarg = "pressure_au"
 
             self.ase_md_kwargs[stress_kwarg] = self.ase_md_kwargs.get(
                 stress_kwarg, self.p_schedule[0] * 1e3 * units.bar
@@ -340,8 +349,8 @@ def run_ase(
             if self.dynamics not in _valid_dynamics[self.ensemble]:
                 raise ValueError(
                     f"{self.dynamics} thermostat not available for "
-                    f"{self.ensemble.value}."
-                    f"Available {self.ensemble.value} thermostats are:"
+                    f"{self.ensemble.value}. "
+                    f"Available {self.ensemble.value} thermostats are: "
                     " ".join(_valid_dynamics[self.ensemble])
                 )
 
@@ -391,7 +400,10 @@ def run_ase(
         def _callback(dyn: MolecularDynamics = md_runner) -> None:
             if self.ensemble == MDEnsemble.nve:
                 return
-            dyn.set_temperature(temperature_K=self.t_schedule[dyn.nsteps])
+            if hasattr(dyn, "_temperature_K"):
+                dyn._temperature_K = self.t_schedule[dyn.nsteps]  # noqa: SLF001
+            else:
+                dyn.set_temperature(temperature_K=self.t_schedule[dyn.nsteps])
             if self.ensemble == MDEnsemble.nvt:
                 return
 

src/atomate2/ase/utils.py

@@ -465,6 +465,22 @@ def relax(
                 write_atoms.calc = self.calculator
             else:
                 write_atoms = atoms
+
+            # ase==3.26.0 change: writing FixAtoms and FixCartesian
+            # constraints to extxyz supported.
+            # Write only these constraints to extxyz
+            if len(write_atoms.constraints) > 0:
+                from ase.constraints import FixAtoms, FixCartesian
+
+                write_atoms_constraints = [
+                    cons
+                    for cons in write_atoms.constraints
+                    if isinstance(cons, FixAtoms | FixCartesian)
+                ]
+                del write_atoms.constraints
+                for cons in write_atoms_constraints:
+                    write_atoms.set_constraint(cons)
+
             ase_write(
                 final_atoms_object_file, write_atoms, format="extxyz", append=True
             )

src/atomate2/common/jobs/mpmorph.py

@@ -16,13 +16,17 @@
 import os
 from pathlib import Path
 from tempfile import TemporaryDirectory
+from typing import TYPE_CHECKING
 
 import numpy as np
 import pandas as pd
 from jobflow import Job
 from pymatgen.core import Composition, Molecule, Structure
 from pymatgen.io.packmol import PackmolBoxGen
 
+if TYPE_CHECKING:
+    from pymatgen.core import Element, Species
+
 _DEFAULT_AVG_VOL_FILE = Path("~/.cache/atomate2").expanduser() / "db_avg_vols.json.gz"
 if not _DEFAULT_AVG_VOL_FILE.parents[0].exists():
     os.makedirs(_DEFAULT_AVG_VOL_FILE.parents[0], exist_ok=True)
@@ -80,14 +84,11 @@ def get_average_volume_from_mp_api(
     with MPRester(api_key=mp_api_key) as mpr:
         comp_entries = mpr.get_entries(composition.reduced_formula, inc_structure=True)
 
-    vols = None
-
-    if len(comp_entries) > 0:
-        vols = [
-            entry.structure.volume / entry.structure.num_sites for entry in comp_entries
-        ]
+    vols = [
+        entry.structure.volume / entry.structure.num_sites for entry in comp_entries
+    ]
 
-    else:
+    if not vols:
         # Find all Materials project entries containing the elements in the
         # desired composition to estimate starting volume.
         with MPRester() as mpr:
@@ -104,6 +105,27 @@ def get_average_volume_from_mp_api(
 
         vols = [entry.structure.volume / entry.structure.num_sites for entry in entries]
 
+    # Fallback: mix atomic volume by relative weight in composition
+    if not vols:
+        by_comp: dict[Element | Species, list[float]] = {
+            ele: [] for ele in composition.elements
+        }
+        for entry in _entries:
+            if len(entry.composition.elements) == 1:
+                by_comp[entry.composition.elements[0]].append(
+                    entry.structure.volume / entry.structure.num_sites
+                )
+        vols = [
+            coeff * np.mean(by_comp[ele]) / composition.num_atoms
+            for ele, coeff in composition.items()
+        ]
+
+        if any(not v for v in by_comp.values()):
+            raise ValueError(
+                "No unary data for "
+                f"{', '.join(str(k) for k, v in by_comp.items() if not v)}."
+            )
+
     return np.mean(vols)
 
 
@@ -250,10 +272,10 @@ def get_entry_from_dict(chem_env: str) -> dict | None:
             return {k: data[k].squeeze() for k in ("avg_vol", "count")}
         return None
 
-    chem_env_key = _get_chem_env_key_from_composition(
+    full_chem_env_key = _get_chem_env_key_from_composition(
         composition, ignore_oxi_states=ignore_oxi_states
     )
-    if (avg_vol := get_entry_from_dict(chem_env_key)) is not None:
+    if (avg_vol := get_entry_from_dict(full_chem_env_key)) is not None:
         return avg_vol["avg_vol"]
 
     vols = []
@@ -269,6 +291,19 @@ def get_entry_from_dict(chem_env: str) -> dict | None:
                 vols.append(avg_vol["avg_vol"] * avg_vol["count"])
                 counts += avg_vol["count"]
 
+    # Fallback, relative weight of monatomic volumes
+    if counts == 0:
+        by_comp = {ele: get_entry_from_dict(ele.name) for ele in composition.elements}
+        if any(v is None for v in by_comp.values()):
+            raise ValueError(
+                "No unary data for "
+                f"{', '.join(str(k) for k, v in by_comp.items() if v is None)}"
+            )
+        return (
+            sum(coeff * by_comp[ele]["avg_vol"] for ele, coeff in composition.items())
+            / composition.num_atoms
+        )
+
     return sum(vols) / counts