Merge branch 'master' of github.com:materialsproject/pymatgen

Author: shyuep

src/pymatgen/analysis/reaction_calculator.py

@@ -10,7 +10,7 @@
 import numpy as np
 from monty.fractions import gcd_float
 from monty.json import MontyDecoder, MSONable
-from uncertainties import ufloat
+from uncertainties import UFloat, ufloat
 
 from pymatgen.core.composition import Composition
 from pymatgen.entries.computed_entries import ComputedEntry
@@ -100,7 +100,7 @@ def __str__(self):
     __repr__ = __str__
 
     @overload
-    def calculate_energy(self, energies: dict[Composition, ufloat]) -> ufloat:
+    def calculate_energy(self, energies: dict[Composition, ufloat]) -> UFloat:
         pass
 
     @overload
@@ -485,10 +485,15 @@ def calculated_reaction_energy_uncertainty(self) -> float:
 
         for entry in self._reactant_entries + self._product_entries:
             comp, factor = entry.composition.get_reduced_composition_and_factor()
-            energy_ufloat = ufloat(entry.energy, entry.correction_uncertainty)
+            energy_ufloat = (
+                ufloat(entry.energy, entry.correction_uncertainty)
+                if entry.correction_uncertainty and not np.isnan(entry.correction_uncertainty)
+                else entry.energy
+            )
             calc_energies[comp] = min(calc_energies.get(comp, float("inf")), energy_ufloat / factor)
 
-        return self.calculate_energy(calc_energies).std_dev
+        ufloat_reaction_energy = self.calculate_energy(calc_energies)
+        return ufloat_reaction_energy.std_dev if isinstance(ufloat_reaction_energy, UFloat) else np.nan
 
     def as_dict(self) -> dict:
         """

src/pymatgen/entries/compatibility.py

@@ -121,16 +121,14 @@ def correct_entry(self, entry):
         """
         new_corr = self.get_correction(entry)
         old_std_dev = entry.correction_uncertainty
-        if np.isnan(old_std_dev):
-            old_std_dev = 0
-        old_corr = ufloat(entry.correction, old_std_dev)
+        old_corr = ufloat(entry.correction, 0 if np.isnan(old_std_dev) else old_std_dev)
         updated_corr = new_corr + old_corr
 
         # if there are no error values available for the corrections applied,
         # set correction uncertainty to not a number
-        uncertainty = np.nan if updated_corr.nominal_value != 0 and updated_corr.std_dev == 0 else updated_corr.std_dev
-
-        entry.energy_adjustments.append(ConstantEnergyAdjustment(updated_corr.nominal_value, uncertainty))
+        entry.energy_adjustments.append(
+            ConstantEnergyAdjustment(updated_corr.nominal_value, updated_corr.std_dev or np.nan)
+        )
 
         return entry
 
@@ -195,7 +193,7 @@ def get_correction(self, entry: AnyComputedEntry) -> ufloat:
             ufloat: 0.0 +/- 0.0 (from uncertainties package)
         """
         if SETTINGS.get("PMG_POTCAR_CHECKS") is False or not self.check_potcar:
-            return ufloat(0.0, 0.0)
+            return ufloat(0.0, np.nan)
 
         potcar_spec = entry.parameters.get("potcar_spec")
         if self.check_hash:
@@ -211,7 +209,7 @@ def get_correction(self, entry: AnyComputedEntry) -> ufloat:
         expected_psp = {self.valid_potcars.get(el.symbol) for el in entry.elements}
         if expected_psp != psp_settings:
             raise CompatibilityError(f"Incompatible POTCAR {psp_settings}, expected {expected_psp}")
-        return ufloat(0.0, 0.0)
+        return ufloat(0.0, np.nan)
 
 
 @cached_class
@@ -249,6 +247,8 @@ def get_correction(self, entry: ComputedEntry | ComputedStructureEntry) -> ufloa
         if rform in self.cpd_energies:
             correction += self.cpd_energies[rform] * comp.num_atoms - entry.uncorrected_energy
 
+        if correction.std_dev == 0:
+            correction = ufloat(correction.nominal_value, np.nan)  # set std_dev to nan if no uncertainty
         return correction
 
 
@@ -286,9 +286,11 @@ def get_correction(self, entry: ComputedEntry | ComputedStructureEntry) -> ufloa
         """
         comp = entry.composition
         if len(comp) == 1:  # Skip element entry
-            return ufloat(0.0, 0.0)
+            return ufloat(0.0, np.nan)
 
-        correction = ufloat(0.0, 0.0)
+        with warnings.catch_warnings():
+            warnings.filterwarnings("ignore", message="Using UFloat")
+            correction = ufloat(0.0, 0.0)
 
         # Check for sulfide corrections
         if Element("S") in comp:
@@ -345,6 +347,8 @@ def get_correction(self, entry: ComputedEntry | ComputedStructureEntry) -> ufloa
             else:
                 correction += self.oxide_correction["oxide"] * comp["O"]
 
+        if correction.std_dev == 0:
+            correction = ufloat(correction.nominal_value, np.nan)  # set std_dev to nan if no uncertainty
         return correction
 
 
@@ -432,6 +436,8 @@ def get_correction(self, entry: ComputedEntry | ComputedStructureEntry) -> ufloa
                 correction += ufloat(-1 * MU_H2O * nH2O, 0.0)
                 # correction += 0.5 * 2.46 * nH2O  # this is the old way this correction was calculated
 
+        if correction.std_dev == 0:
+            correction = ufloat(correction.nominal_value, np.nan)  # set std_dev to nan if no uncertainty
         return correction
 
 
@@ -537,6 +543,8 @@ def get_correction(self, entry: ComputedEntry | ComputedStructureEntry) -> ufloa
             if sym in u_corr:
                 correction += ufloat(u_corr[sym], u_errors[sym]) * comp[el]
 
+        if correction.std_dev == 0:
+            correction = ufloat(correction.nominal_value, np.nan)  # set std_dev to nan if no uncertainty
         return correction
 
 
@@ -1076,7 +1084,7 @@ def get_adjustments(self, entry: AnyComputedEntry) -> list[CompositionEnergyAdju
             )
 
         # check the POTCAR symbols
-        # this should return ufloat(0, 0) or raise a CompatibilityError or ValueError
+        # this should return ufloat(0, np.nan) or raise a CompatibilityError or ValueError
         if entry.parameters.get("software", "vasp") == "vasp":
             pc = PotcarCorrection(
                 MPRelaxSet,

src/pymatgen/entries/computed_entries.py

@@ -363,8 +363,10 @@ def correction(self) -> float:
         Returns:
             float: the total energy correction / adjustment applied to the entry in eV.
         """
-        # adds to ufloat(0.0, 0.0) to ensure that no corrections still result in ufloat object
-        corr = ufloat(0.0, 0.0) + sum(ufloat(ea.value, ea.uncertainty) for ea in self.energy_adjustments)
+        # either sum of adjustments or ufloat with nan std_dev, so that no corrections still result in ufloat object:
+        corr = sum(ufloat(ea.value, ea.uncertainty) for ea in self.energy_adjustments if ea.value) or ufloat(
+            0.0, np.nan
+        )
         return corr.nominal_value
 
     @correction.setter
@@ -386,11 +388,11 @@ def correction_uncertainty(self) -> float:
         Returns:
             float: the uncertainty of the energy adjustments applied to the entry in eV.
         """
-        # adds to ufloat(0.0, 0.0) to ensure that no corrections still result in ufloat object
-        unc = ufloat(0.0, 0.0) + sum(
-            (ufloat(ea.value, ea.uncertainty) if not np.isnan(ea.uncertainty) else ufloat(ea.value, 0))
+        # either sum of adjustments or ufloat with nan std_dev, so that no corrections still result in ufloat object:
+        unc = sum(
+            (ufloat(ea.value, ea.uncertainty) if not np.isnan(ea.uncertainty) and ea.value else ufloat(ea.value, 0))
             for ea in self.energy_adjustments
-        )
+        ) or ufloat(0.0, np.nan)
 
         if unc.nominal_value != 0 and unc.std_dev == 0:
             return np.nan

src/pymatgen/io/jdftx/inputs.py

@@ -20,7 +20,7 @@ class is written.
 import scipy.constants as const
 from monty.json import MSONable
 
-from pymatgen.core import Structure
+from pymatgen.core import Lattice, Structure
 from pymatgen.core.periodic_table import Element
 from pymatgen.core.units import bohr_to_ang
 from pymatgen.io.jdftx.generic_tags import AbstractTag, BoolTagContainer, DumpTagContainer, MultiformatTag, TagContainer
@@ -788,17 +788,8 @@ def from_jdftxinfile(cls, jdftxinfile: JDFTXInfile, sort_structure: bool = False
         Returns:
             JDFTXStructure: The created JDFTXStructure object.
         """
-        jl = jdftxinfile["lattice"]
-        lattice = np.zeros([3, 3])
-        for i in range(3):
-            for j in range(3):
-                lattice[i][j] += float(jl[f"R{i}{j}"])
-        if "latt-scale" in jdftxinfile:
-            latt_scale = np.array([jdftxinfile["latt-scale"][x] for x in ["s0", "s1", "s2"]])
-            lattice *= latt_scale
-        lattice = lattice.T  # convert to row vector format
-        lattice *= const.value("Bohr radius") * 10**10  # Bohr radius in Ang; convert to Ang
 
+        lattice = _infile_to_pmg_lattice(jdftxinfile)
         atomic_symbols = [x["species-id"] for x in jdftxinfile["ion"]]
         coords = np.array([[x["x0"], x["x1"], x["x2"]] for x in jdftxinfile["ion"]])
         coords *= const.value("Bohr radius") * 10**10  # Bohr radius in Ang; convert to Ang
@@ -1016,3 +1007,58 @@ def _multi_getattr(varbase: Any, varname: str):
     for var in varlist:
         varbase = getattr(varbase, var)
     return varbase
+
+
+def _infile_to_pmg_lattice(jdftxinfile: JDFTXInfile) -> Lattice:
+    jl = jdftxinfile["lattice"]
+    if "R00" not in jl:
+        return _infile_special_format_to_pmg_lattice(jdftxinfile)
+    return _infile_matrix_format_to_pmg_lattice(jdftxinfile)
+
+
+def _infile_special_format_to_pmg_lattice(jdftxinfile: JDFTXInfile) -> Lattice:
+    jl = jdftxinfile["lattice"]
+    if "modification" in jl:
+        raise NotImplementedError("Special case lattices with modification not implemented yet")
+    # Second boolean to check if latt-scale was passed but does nothing
+    if ("latt-scale" in jl) and (not all(np.isclose(float(x), 1) in jl for x in ["s0", "s1", "s2"])):
+        raise NotImplementedError("latt-scale for special case lattices not implemented yet")
+    if "Monoclinic" in jl:
+        lattice = Lattice.monoclinic(jl["a"], jl["b"], jl["c"], jl["beta"])
+    elif "Orthorhombic" in jl:
+        lattice = Lattice.orthorhombic(jl["a"], jl["b"], jl["c"])
+    elif "Cubic" in jl:
+        lattice = Lattice.cubic(jl["a"])
+    elif "Hexagonal" in jl:
+        lattice = Lattice.hexagonal(jl["a"], jl["c"])
+    elif "Rhombohedral" in jl:
+        lattice = Lattice.rhombohedral(jl["a"], jl["alpha"])
+    elif "Tetragonal" in jl:
+        lattice = Lattice.tetragonal(jl["a"], jl["c"])
+    elif "Triclinic" in jl:
+        lattice = Lattice.from_parameters(jl["a"], jl["b"], jl["c"], jl["alpha"], jl["beta"], jl["gamma"])
+    else:
+        raise ValueError(f"Unable to convert JDFTX lattice tag {jl} to pymatgen Lattice object")
+    return lattice
+
+
+def _infile_matrix_format_to_pmg_lattice(jdftxinfile: JDFTXInfile) -> Lattice:
+    """Convert JDFTX lattice tag to pymatgen Lattice object.
+
+    Args:
+        jl (dict[str, Any]): JDFTX lattice tag.
+
+    Returns:
+        Lattice: Pymatgen Lattice object.
+    """
+    _lattice = np.zeros([3, 3])
+    jl = jdftxinfile["lattice"]
+    for i in range(3):
+        for j in range(3):
+            _lattice[i][j] += float(jl[f"R{i}{j}"])
+    if "latt-scale" in jdftxinfile:
+        latt_scale = np.array([jdftxinfile["latt-scale"][x] for x in ["s0", "s1", "s2"]])
+        _lattice *= latt_scale
+    _lattice = _lattice.T  # convert to row vector format
+    _lattice *= const.value("Bohr radius") * 10**10  # Bohr radius in Ang; convert to Ang
+    return Lattice(_lattice)

src/pymatgen/symmetry/analyzer.py

@@ -771,95 +771,79 @@ def get_conventional_standard_structure(
                     latt2 = Lattice([m[tp2[0]], m[tp2[1]], m[2]])
                     lengths = latt2.lengths
                     angles = latt2.angles
-                    if angles[0] > 90:
-                        # if the angle is > 90 we invert a and b to get
-                        # an angle < 90
-                        a, b, c, alpha, beta, gamma = Lattice([-m[tp2[0]], -m[tp2[1]], m[2]]).parameters
-                        transf = np.zeros(shape=(3, 3))
+                    alpha_degrees = angles[0]
+                    if alpha_degrees == 90:
+                        continue
+
+                    transf = np.zeros(shape=(3, 3))
+                    transf[2][2] = 1
+
+                    if alpha_degrees > 90:
+                        # if the angle is > 90 we invert a and b to get an angle < 90
+                        a, b, c, alpha_degrees, beta, gamma = Lattice([-m[tp2[0]], -m[tp2[1]], m[2]]).parameters
                         transf[0][tp2[0]] = -1
                         transf[1][tp2[1]] = -1
-                        transf[2][2] = 1
-                        alpha = math.pi * alpha / 180
-                        new_matrix = [
-                            [a, 0, 0],
-                            [0, b, 0],
-                            [0, c * cos(alpha), c * sin(alpha)],
-                        ]
-                        continue
 
-                    if angles[0] < 90:
-                        transf = np.zeros(shape=(3, 3))
+                    elif alpha_degrees < 90:  # otherwise no inversion
                         transf[0][tp2[0]] = 1
                         transf[1][tp2[1]] = 1
-                        transf[2][2] = 1
                         a, b, c = lengths
-                        alpha = math.pi * angles[0] / 180
-                        new_matrix = [
-                            [a, 0, 0],
-                            [0, b, 0],
-                            [0, c * cos(alpha), c * sin(alpha)],
-                        ]
+                        alpha_degrees = angles[0]
+
+                    alpha_radians = math.pi * alpha_degrees / 180
+                    new_matrix = [
+                        [a, 0, 0],
+                        [0, b, 0],
+                        [0, c * cos(alpha_radians), c * sin(alpha_radians)],
+                    ]
 
                 if new_matrix is None:
-                    # this if is to treat the case
-                    # where alpha==90 (but we still have a monoclinic sg
+                    # this if is to treat the case where alpha==90 (but we still have a monoclinic sg), with C unchanged
+                    # transf is already defined as [[0,0,0],[0,0,0],[0,0,1]], and sorted_dic only contains a & b:
                     new_matrix = [[a, 0, 0], [0, b, 0], [0, 0, c]]
-                    transf = np.zeros(shape=(3, 3))
-                    transf[2] = [0, 0, 1]  # see issue #1929
                     for idx, dct in enumerate(sorted_dic):
                         transf[idx][dct["orig_index"]] = 1
+
             # if not C-setting
             else:
                 # try all permutations of the axis
-                # keep the ones with the non-90 angle=alpha
-                # and b<c
+                # keep the ones with the non-90 angle=alpha and b<c
                 new_matrix = None
 
                 for tp3 in itertools.permutations(list(range(3)), 3):
                     m = lattice.matrix
                     a, b, c, alpha, beta, gamma = Lattice([m[tp3[0]], m[tp3[1]], m[tp3[2]]]).parameters
+                    if alpha == 90 or b >= c:
+                        continue
+                    transf = np.zeros(shape=(3, 3))
+                    transf[2][tp3[2]] = 1
+
                     if alpha > 90 and b < c:
                         a, b, c, alpha, beta, gamma = Lattice([-m[tp3[0]], -m[tp3[1]], m[tp3[2]]]).parameters
-                        transf = np.zeros(shape=(3, 3))
                         transf[0][tp3[0]] = -1
                         transf[1][tp3[1]] = -1
-                        transf[2][tp3[2]] = 1
-                        alpha = math.pi * alpha / 180
-                        new_matrix = [
-                            [a, 0, 0],
-                            [0, b, 0],
-                            [0, c * cos(alpha), c * sin(alpha)],
-                        ]
-                        continue
 
-                    if alpha < 90 and b < c:
-                        transf = np.zeros(shape=(3, 3))
+                    elif alpha < 90 and b < c:
                         transf[0][tp3[0]] = 1
                         transf[1][tp3[1]] = 1
-                        transf[2][tp3[2]] = 1
-                        alpha = math.pi * alpha / 180
-                        new_matrix = [
-                            [a, 0, 0],
-                            [0, b, 0],
-                            [0, c * cos(alpha), c * sin(alpha)],
-                        ]
 
-                if new_matrix is None:
-                    # this if is to treat the case
-                    # where alpha==90 (but we still have a monoclinic sg
+                    alpha = math.pi * alpha / 180
                     new_matrix = [
-                        [sorted_lengths[0], 0, 0],
-                        [0, sorted_lengths[1], 0],
-                        [0, 0, sorted_lengths[2]],
+                        [a, 0, 0],
+                        [0, b, 0],
+                        [0, c * cos(alpha), c * sin(alpha)],
                     ]
-                    transf = np.zeros(shape=(3, 3))
-                    for idx, dct in enumerate(sorted_dic):
+
+                if new_matrix is None:
+                    # this if is to treat the case where alpha==90 (but we still have a monoclinic sg)
+                    new_matrix = list(np.zeros(shape=(3, 3)))
+                    for idx, dct in enumerate(sorted_dic):  # vectors sorted by length (a<=b<=c)
+                        new_matrix[idx] = dct["vec"]
                         transf[idx][dct["orig_index"]] = 1
 
             if international_monoclinic:
                 # The above code makes alpha the non-right angle.
-                # The following will convert to proper international convention
-                # that beta is the non-right angle.
+                # The following will convert to proper international convention that beta is the non-right angle
                 op = [[0, 1, 0], [1, 0, 0], [0, 0, -1]]
                 transf = np.dot(op, transf)
                 new_matrix = np.dot(op, new_matrix)