Merge remote-tracking branch 'upstream/master' into 4158-fix-eq-check

Author: DanielYang59

README.md

@@ -71,7 +71,10 @@ See [GitHub releases](https://github.com/materialsproject/pymatgen/releases), [`
 
 ## Using pymatgen
 
-Please refer to the official [`pymatgen` docs] for tutorials and examples.
+Please refer to the official [`pymatgen` docs] for tutorials and examples. Dr Anubhav Jain (@computron) has also created
+a series of [tutorials](https://github.com/computron/pymatgen_tutorials)
+and [YouTube videos](https://www.youtube.com/playlist?list=PL7gkuUui8u7_M47KrV4tS4pLwhe7mDAjT), which is a good
+resource, especially for beginners.
 
 ## How to cite pymatgen
 

docs/index.md

@@ -24,7 +24,7 @@
 from pymatgen.util.string import Stringify, formula_double_format
 
 if TYPE_CHECKING:
-    from collections.abc import Generator, Iterator
+    from collections.abc import Generator, ItemsView, Iterator
     from typing import Any, ClassVar, Literal
 
     from typing_extensions import Self
@@ -173,7 +173,7 @@ def __init__(self, *args, strict: bool = False, **kwargs) -> None:
         else:
             elem_map = dict(*args, **kwargs)  # type: ignore[assignment]
         elem_amt = {}
-        self._n_atoms = 0
+        self._n_atoms: int | float = 0
         for key, val in elem_map.items():
             if val < -type(self).amount_tolerance and not self.allow_negative:
                 raise ValueError("Amounts in Composition cannot be negative!")
@@ -184,7 +184,7 @@ def __init__(self, *args, strict: bool = False, **kwargs) -> None:
         if strict and not self.valid:
             raise ValueError(f"Composition is not valid, contains: {', '.join(map(str, self.elements))}")
 
-    def __getitem__(self, key: SpeciesLike) -> float:
+    def __getitem__(self, key: SpeciesLike) -> int | float:
         try:
             sp = get_el_sp(key)
             if isinstance(sp, Species):
@@ -202,6 +202,10 @@ def __len__(self) -> int:
     def __iter__(self) -> Iterator[Species | Element | DummySpecies]:
         return iter(self._data)
 
+    def items(self) -> ItemsView[Species | Element | DummySpecies, int | float]:
+        """Returns Dict.items() for the Composition dict."""
+        return self._data.items()
+
     def __contains__(self, key) -> bool:
         try:
             sp = get_el_sp(key)

src/pymatgen/core/composition.py

@@ -8,23 +8,30 @@
 import warnings
 from fnmatch import fnmatch
 from pathlib import Path
+from tempfile import NamedTemporaryFile
 from typing import TYPE_CHECKING, TypeAlias, cast
 
 import numpy as np
 from monty.io import zopen
 from monty.json import MSONable
 
 from pymatgen.core.structure import Composition, DummySpecies, Element, Lattice, Molecule, Species, Structure
+from pymatgen.io.ase import NO_ASE_ERR, AseAtomsAdaptor
+
+if NO_ASE_ERR is None:
+    from ase.io.trajectory import Trajectory as AseTrajectory
+else:
+    AseTrajectory = None
+
 
 if TYPE_CHECKING:
-    from collections.abc import Iterator
+    from collections.abc import Iterator, Sequence
     from typing import Any
 
     from typing_extensions import Self
 
     from pymatgen.util.typing import Matrix3D, PathLike, SitePropsType, Vector3D
 
-
 __author__ = "Eric Sivonxay, Shyam Dwaraknath, Mingjian Wen, Evan Spotte-Smith"
 __version__ = "0.1"
 __date__ = "Jun 29, 2022"
@@ -563,8 +570,6 @@ def from_file(cls, filename: str | Path, constant_lattice: bool = True, **kwargs
             Trajectory: containing the structures or molecules in the file.
         """
         filename = str(Path(filename).expanduser().resolve())
-        is_mol = False
-        molecules = []
         structures = []
 
         if fnmatch(filename, "*XDATCAR*"):
@@ -578,31 +583,24 @@ def from_file(cls, filename: str | Path, constant_lattice: bool = True, **kwargs
             structures = Vasprun(filename).structures
 
         elif fnmatch(filename, "*.traj"):
-            try:
-                from ase.io.trajectory import Trajectory as AseTrajectory
-
-                from pymatgen.io.ase import AseAtomsAdaptor
-
-                ase_traj = AseTrajectory(filename)
-                # Periodic boundary conditions should be the same for all frames so just check the first
-                pbc = ase_traj[0].pbc
+            if NO_ASE_ERR is None:
+                return cls.from_ase(
+                    filename,
+                    constant_lattice=constant_lattice,
+                    store_frame_properties=True,
+                    additional_fields=None,
+                )
+            raise ImportError("ASE is required to read .traj files. pip install ase")
 
-                if any(pbc):
-                    structures = [AseAtomsAdaptor.get_structure(atoms) for atoms in ase_traj]
-                else:
-                    molecules = [AseAtomsAdaptor.get_molecule(atoms) for atoms in ase_traj]
-                    is_mol = True
+        elif fnmatch(filename, "*.json*"):
+            from monty.serialization import loadfn
 
-            except ImportError as exc:
-                raise ImportError("ASE is required to read .traj files. pip install ase") from exc
+            return loadfn(filename, **kwargs)
 
         else:
-            supported_file_types = ("XDATCAR", "vasprun.xml", "*.traj")
+            supported_file_types = ("XDATCAR", "vasprun.xml", "*.traj", ".json")
             raise ValueError(f"Expect file to be one of {supported_file_types}; got {filename}.")
 
-        if is_mol:
-            return cls.from_molecules(molecules, **kwargs)
-
         return cls.from_structures(structures, constant_lattice=constant_lattice, **kwargs)
 
     @staticmethod
@@ -734,3 +732,148 @@ def _get_site_props(self, frames: ValidIndex) -> SitePropsType | None:
                 return [self.site_properties[idx] for idx in frames]
             raise ValueError("Unexpected frames type.")
         raise ValueError("Unexpected site_properties type.")
+
+    @classmethod
+    def from_ase(
+        cls,
+        trajectory: str | Path | AseTrajectory,
+        constant_lattice: bool | None = None,
+        store_frame_properties: bool = True,
+        property_map: dict[str, str] | None = None,
+        lattice_match_tol: float = 1.0e-6,
+        additional_fields: Sequence[str] | None = ["temperature", "velocities"],
+    ) -> Trajectory:
+        """
+        Convert an ASE trajectory to a pymatgen trajectory.
+
+        Args:
+            trajectory (str, .Path, or ASE .Trajectory) : the ASE trajectory, or a file path to it if a str or .Path
+            constant_lattice (bool or None) : if a bool, whether the lattice is constant in the .Trajectory.
+                If `None`, this is determined on the fly.
+            store_frame_properties (bool) : Whether to store pymatgen .Trajectory `frame_properties` as
+                ASE calculator properties. Defaults to True
+            property_map (dict[str,str]) : A mapping between ASE calculator properties and
+                pymatgen .Trajectory `frame_properties` keys. Ex.:
+                    property_map = {"energy": "e_0_energy"}
+                would map `e_0_energy` in the pymatgen .Trajectory `frame_properties`
+                to ASE's `get_potential_energy` function.
+                See `ase.calculators.calculator.all_properties` for a list of acceptable calculator properties.
+            lattice_match_tol (float = 1.0e-6) : tolerance to which lattices are matched if
+                `constant_lattice = None`.
+            additional_fields (Sequence of str, defaults to ["temperature", "velocities"]) :
+                Optional other fields to save in the pymatgen .Trajectory.
+                Valid options are "temperature" and "velocities".
+
+        Returns:
+            pymatgen .Trajectory
+        """
+        if isinstance(trajectory, str | Path):
+            trajectory = AseTrajectory(trajectory, "r")
+
+        property_map = property_map or {
+            "energy": "energy",
+            "forces": "forces",
+            "stress": "stress",
+        }
+        additional_fields = additional_fields or []
+
+        adaptor = AseAtomsAdaptor()
+
+        structures = []
+        frame_properties = []
+        converter = adaptor.get_structure if (is_pbc := any(trajectory[0].pbc)) else adaptor.get_molecule
+
+        for atoms in trajectory:
+            site_properties = {}
+            if "velocities" in additional_fields:
+                site_properties["velocities"] = atoms.get_velocities()
+
+            structures.append(converter(atoms, site_properties=site_properties))
+
+            if store_frame_properties and atoms.calc:
+                props = {v: atoms.calc.get_property(k) for k, v in property_map.items()}
+                if "temperature" in additional_fields:
+                    props["temperature"] = atoms.get_temperature()
+
+                frame_properties.append(props)
+
+        if constant_lattice is None:
+            constant_lattice = all(
+                np.all(np.abs(ref_struct.lattice.matrix - structures[j].lattice.matrix)) < lattice_match_tol
+                for i, ref_struct in enumerate(structures)
+                for j in range(i + 1, len(structures))
+            )
+
+        if is_pbc:
+            return cls.from_structures(structures, constant_lattice=constant_lattice, frame_properties=frame_properties)
+        return cls.from_molecules(
+            structures,
+            constant_lattice=constant_lattice,
+            frame_properties=frame_properties,
+        )
+
+    def to_ase(
+        self,
+        property_map: dict[str, str] | None = None,
+        ase_traj_file: str | Path | None = None,
+    ) -> AseTrajectory:
+        """
+        Convert a pymatgen .Trajectory to an ASE .Trajectory.
+
+        Args:
+            trajectory (pymatgen .Trajectory) : trajectory to convert
+            property_map (dict[str,str]) : A mapping between ASE calculator properties and
+                pymatgen .Trajectory `frame_properties` keys. Ex.:
+                    property_map = {"energy": "e_0_energy"}
+                would map `e_0_energy` in the pymatgen .Trajectory `frame_properties`
+                to ASE's `get_potential_energy` function.
+                See `ase.calculators.calculator.all_properties` for a list of acceptable calculator properties.
+            ase_traj_file (str, Path, or None (default) ) :  If not None, the name of
+                the file to write the ASE trajectory to.
+
+        Returns:
+            ase .Trajectory
+        """
+        if NO_ASE_ERR is not None:
+            raise ImportError("ASE is required to write .traj files. pip install ase")
+
+        from ase.calculators.calculator import all_properties
+        from ase.calculators.singlepoint import SinglePointCalculator
+
+        property_map = property_map or {
+            "energy": "energy",
+            "forces": "forces",
+            "stress": "stress",
+        }
+
+        if (unrecognized_props := set(property_map).difference(set(all_properties))) != set():
+            raise ValueError(f"Unrecognized ASE calculator properties:\n{', '.join(unrecognized_props)}")
+
+        adaptor = AseAtomsAdaptor()
+
+        temp_file = None
+        if ase_traj_file is None:
+            temp_file = NamedTemporaryFile(delete=False)  # noqa: SIM115
+            ase_traj_file = temp_file.name
+
+        frame_props = self.frame_properties or [{} for _ in range(len(self))]
+        for idx, structure in enumerate(self):
+            atoms = adaptor.get_atoms(structure, msonable=False, velocities=structure.site_properties.get("velocities"))
+
+            props: dict[str, Any] = {k: frame_props[idx][v] for k, v in property_map.items() if v in frame_props[idx]}
+
+            # Ensure that `charges` and `magmoms` are not lost from AseAtomsAdaptor
+            for k in ("charges", "magmoms"):
+                if k in atoms.calc.implemented_properties or k in atoms.calc.results:
+                    props[k] = atoms.calc.get_property(k)
+
+            atoms.calc = SinglePointCalculator(atoms=atoms, **props)
+
+            with AseTrajectory(ase_traj_file, "a" if idx > 0 else "w", atoms=atoms) as _traj_file:
+                _traj_file.write()
+
+        ase_traj = AseTrajectory(ase_traj_file, "r")
+        if temp_file is not None:
+            temp_file.close()
+
+        return ase_traj

src/pymatgen/core/trajectory.py

@@ -156,14 +156,13 @@ def get_cbm_vbm(self, tol: float = 1e-4, abs_tol: bool = False, spin: Spin | Non
 
         # Work backwards until tolerance is reached
         i_gap_start = i_fermi
-        while i_gap_start >= 1 and tdos[i_gap_start - 1] <= tol:
+        while i_gap_start >= 1 and tdos[i_gap_start] <= tol:
             i_gap_start -= 1
 
         # Work forwards until tolerance is reached
-        i_gap_end = i_gap_start
+        i_gap_end = i_gap_start + 1
         while i_gap_end < tdos.shape[0] and tdos[i_gap_end] <= tol:
             i_gap_end += 1
-        i_gap_end -= 1
 
         return self.x[i_gap_end], self.x[i_gap_start]
 
@@ -378,16 +377,15 @@ def get_cbm_vbm(self, tol: float = 1e-4, abs_tol: bool = False, spin: Spin | Non
 
         # Work backwards until tolerance is reached
         i_gap_start = i_fermi
-        while i_gap_start >= 1 and tdos[i_gap_start - 1] <= tol:
+        while i_gap_start >= 1 and tdos[i_gap_start] <= tol:
             i_gap_start -= 1
 
         # Work forwards until tolerance is reached
-        i_gap_end = i_gap_start
+        i_gap_end = i_gap_start + 1
         while i_gap_end < tdos.shape[0] and tdos[i_gap_end] <= tol:
             i_gap_end += 1
-        i_gap_end -= 1
 
-        return self.energies[i_gap_end], self.energies[i_gap_start]
+        return (self.energies[min(i_gap_end, len(self.energies) - 1)], self.energies[max(i_gap_start, 0)])
 
     def get_gap(
         self,
@@ -518,15 +516,15 @@ def get_doping(self, fermi_level: float, temperature: float) -> float:
                 (P-type).
         """
         cb_integral = np.sum(
-            self.tdos[self.idx_mid_gap :]
-            * f0(self.energies[self.idx_mid_gap :], fermi_level, temperature)
-            * self.de[self.idx_mid_gap :],
+            self.tdos[max(self.idx_mid_gap, self.idx_vbm + 1) :]
+            * f0(self.energies[max(self.idx_mid_gap, self.idx_vbm + 1) :], fermi_level, temperature)
+            * self.de[max(self.idx_mid_gap, self.idx_vbm + 1) :],
             axis=0,
         )
         vb_integral = np.sum(
-            self.tdos[: self.idx_mid_gap + 1]
-            * f0(-self.energies[: self.idx_mid_gap + 1], -fermi_level, temperature)
-            * self.de[: self.idx_mid_gap + 1],
+            self.tdos[: min(self.idx_mid_gap, self.idx_cbm - 1) + 1]
+            * f0(-self.energies[: min(self.idx_mid_gap, self.idx_cbm - 1) + 1], -fermi_level, temperature)
+            * self.de[: min(self.idx_mid_gap, self.idx_cbm - 1) + 1],
             axis=0,
         )
         return (vb_integral - cb_integral) / (self.volume * self.A_to_cm**3)

src/pymatgen/electronic_structure/dos.py

@@ -607,6 +607,10 @@ def get_content(
                     content += self.get_aims_control_parameter_str(key, output_type, "%s")
             elif key == "vdw_correction_hirshfeld" and value:
                 content += self.get_aims_control_parameter_str(key, "", "%s")
+            elif key == "xc":
+                if "libxc" in value:
+                    content += self.get_aims_control_parameter_str("override_warning_libxc", ".true.", "%s")
+                content += self.get_aims_control_parameter_str(key, value, "%s")
             elif isinstance(value, bool):
                 content += self.get_aims_control_parameter_str(key, str(value).lower(), ".%s.")
             elif isinstance(value, tuple | list):

src/pymatgen/io/aims/inputs.py

@@ -144,6 +144,10 @@ def test_init(self):
             with pytest.raises(TypeError, match=f"{type(val).__name__!r} object is not iterable"):
                 Composition(val)
 
+    def test_init_mixed_valence(self):
+        assert Composition({"Fe3+": 2, "Fe2+": 2, "Li": 4, "O": 16, "P": 4}).formula == "Li4 Fe4 P4 O16"
+        assert Composition({"Fe3+": 2, "Fe": 2, "Li": 4, "O": 16, "P": 4}).formula == "Li4 Fe4 P4 O16"
+
     def test_str_and_repr(self):
         test_cases = [
             (
@@ -621,8 +625,8 @@ def test_negative_compositions(self):
         # test species
         c1 = Composition({"Mg": 1, "Mg2+": -1}, allow_negative=True)
         assert c1.num_atoms == 2
-        assert c1.element_composition == Composition("Mg-1", allow_negative=True)
-        assert c1.average_electroneg == approx(0.655)
+        assert c1.get_el_amt_dict() == {"Mg": 0}
+        assert c1.average_electroneg == approx(1.31)  # correct Mg electronegativity)
 
     def test_special_formulas(self):
         special_formulas = {