Add CIF checking for incorrect/missing elements, support for isotopes (primarily Deuterium and Tritium), and correct handling of new VASP POSCAR formats

Author: esoteric-ephemera

dev_scripts/periodic_table.yaml

@@ -21,20 +21,17 @@
 
 
 def test_yaml():
-    with open("periodic_table.yaml") as file:
-        data = yaml.load(file)
-        print(data)
+    loadfn("periodic_table.yaml")
 
 
 def test_json():
-    with open("periodic_table.json") as file:
+    with open("../pymatgen/core/periodic_table.json") as file:
         data = json.load(file)
         print(data)
 
 
 def parse_oxi_state():
-    with open("periodic_table.yaml") as file:
-        data = yaml.load(file)
+    data = loadfn("periodic_table.yaml")
     with open("oxidation_states.txt") as file:
         oxi_data = file.read()
     oxi_data = re.sub("[\n\r]", "", oxi_data)
@@ -72,8 +69,7 @@ def parse_oxi_state():
 
 
 def parse_ionic_radii():
-    with open("periodic_table.yaml") as f:
-        data = yaml.load(f)
+    data = loadfn("periodic_table.yaml")
     with open("ionic_radii.csv") as f:
         radii_data = f.read()
     radii_data = radii_data.split("\r")
@@ -104,8 +100,7 @@ def parse_ionic_radii():
 
 
 def parse_radii():
-    with open("periodic_table.yaml") as f:
-        data = yaml.load(f)
+    data = loadfn("periodic_table.yaml")
     with open("radii.csv") as f:
         radii_data = f.read()
     radii_data = radii_data.split("\r")
@@ -136,13 +131,12 @@ def parse_radii():
             print(el)
     with open("periodic_table2.yaml", "w") as f:
         yaml.dump(data, f)
-    with open("periodic_table.json", "w") as f:
+    with open("../pymatgen/core/periodic_table.json", "w") as f:
         json.dump(data, f)
 
 
 def update_ionic_radii():
-    with open("periodic_table.yaml") as f:
-        data = yaml.load(f)
+    data = loadfn("periodic_table.yaml")
 
     for d in data.values():
         if "Ionic_radii" in d:
@@ -156,13 +150,12 @@ def update_ionic_radii():
             del d["Ionic_radii_ls"]
     with open("periodic_table2.yaml", "w") as f:
         yaml.dump(data, f)
-    with open("periodic_table.json", "w") as f:
+    with open("../pymatgen/core/periodic_table.json", "w") as f:
         json.dump(data, f)
 
 
 def parse_shannon_radii():
-    with open("periodic_table.yaml") as f:
-        data = yaml.load(f)
+    data = loadfn("periodic_table.yaml")
 
     from openpyxl import load_workbook
 
@@ -194,22 +187,20 @@ def parse_shannon_radii():
         if el in data:
             data[el]["Shannon radii"] = dict(radii[el])
 
-    with open("periodic_table.yaml", "w") as f:
-        yaml.safe_dump(data, f)
-    with open("periodic_table.json", "w") as f:
+    dumpfn(data, "periodic_table.yaml")
+    with open("../pymatgen/core/periodic_table.json", "w") as f:
         json.dump(data, f)
 
 
 def gen_periodic_table():
-    with open("periodic_table.yaml") as f:
-        data = yaml.load(f)
+    data = loadfn("periodic_table.yaml")
 
-    with open("periodic_table.json", "w") as f:
+    with open("../pymatgen/core/periodic_table.json", "w") as f:
         json.dump(data, f)
 
 
 def gen_iupac_ordering():
-    periodic_table = loadfn("periodic_table.json")
+    periodic_table = loadfn("../pymatgen/core/periodic_table.json")
     order = [
         ([18], range(6, 0, -1)),  # noble gasses
         ([1], range(7, 1, -1)),  # alkali metals
@@ -265,12 +256,24 @@ def add_electron_affinities():
             row.append(td.get_text().strip())
         data.append(row)
     data.pop(0)
-    ea = {int(r[0]): float(re.sub(r"[\s\(\)]", "", r[3].strip("()[]"))) for r in data}
-    assert set(ea).issuperset(range(1, 93))  # Ensure that we have data for up to U.
+
+    ea = {}
+    max_Z = max(Element(element).Z for element in Element.__members__)
+    for r in data:
+        # don't want superheavy elements or less common isotopes
+        if int(r[0]) > max_Z or r[2] in ea:
+            continue
+        tempstr = re.sub(r"[\s\(\)]", "", r[3].strip("()[]"))
+        # hyphen-like characters used that can't be parsed by .float
+        bytesrep = tempstr.encode("unicode_escape").replace(b"\\u2212", b"-")
+        ea[r[2]] = float(bytesrep.decode("unicode_escape"))
+
+    Z_set = {Element.from_name(element).Z for element in ea}
+    assert Z_set.issuperset(range(1, 93))  # Ensure that we have data for up to U.
     print(ea)
     pt = loadfn("../pymatgen/core/periodic_table.json")
     for k, v in pt.items():
-        v["Electron affinity"] = ea.get(Element(k).Z)
+        v["Electron affinity"] = ea.get(Element(k).long_name)
     dumpfn(pt, "../pymatgen/core/periodic_table.json")
 
 
@@ -296,12 +299,12 @@ def add_ionization_energies():
     pt = loadfn("../pymatgen/core/periodic_table.json")
     for k, v in pt.items():
         del v["Ionization energy"]
-        v["Ionization energies"] = data.get(Element(k).Z, [])
+        v["Ionization energies"] = data.get(Element(k).long_name, [])
     dumpfn(pt, "../pymatgen/core/periodic_table.json")
 
 
 if __name__ == "__main__":
-    # parse_shannon_radii()
+    parse_shannon_radii()
     # add_ionization_energies()
     add_electron_affinities()
     # gen_periodic_table()

dev_scripts/update_pt_data.py

@@ -128,6 +128,7 @@ def __init__(self, *args, strict: bool = False, **kwargs) -> None:
             raise ValueError("float('NaN') is not a valid Composition, did you mean str('NaN')?")
         else:
             elem_map = dict(*args, **kwargs)  # type: ignore
+
         elem_amt = {}
         self._n_atoms = 0
         for key, val in elem_map.items():

pymatgen/core/composition.py

@@ -54,6 +54,7 @@ def __init__(self, symbol: SpeciesLike) -> None:
             Z (int): Atomic number.
             symbol (str): Element symbol.
             long_name (str): Long name for element. E.g., "Hydrogen".
+            A (int) : Atomic mass number (number of protons plus number of neutrons).
             atomic_radius_calculated (float): Calculated atomic radius for the element. This is the empirical value.
                 Data is obtained from http://wikipedia.org/wiki/Atomic_radii_of_the_elements_(data_page).
             van_der_waals_radius (float): Van der Waals radius for the element. This is the empirical value determined
@@ -101,12 +102,26 @@ def __init__(self, symbol: SpeciesLike) -> None:
         # Store key variables for quick access
         self.Z = data["Atomic no"]
 
+        self._is_named_isotope = data.get("Is named isotope", False)
+        if self._is_named_isotope:
+            for sym in _pt_data:
+                if _pt_data[sym]["Atomic no"] == self.Z and not _pt_data[sym].get("Is named isotope", False):
+                    self.symbol = sym
+                    break
+
         at_r = data.get("Atomic radius", "no data")
         if str(at_r).startswith("no data"):
             self._atomic_radius = None
         else:
             self._atomic_radius = Length(at_r, "ang")
         self._atomic_mass = Mass(data["Atomic mass"], "amu")
+
+        self._atomic_mass_number = None
+        self.A = None
+        if data.get("Atomic mass no"):
+            self.A = data.get("Atomic mass no")
+            self._atomic_mass_number = Mass(data["Atomic mass no"], "amu")
+
         self.long_name = data["Name"]
         self._data = data
 
@@ -140,6 +155,14 @@ def atomic_mass(self) -> FloatWithUnit:
         """
         return self._atomic_mass
 
+    @property
+    def atomic_mass_number(self) -> FloatWithUnit | None:
+        """
+        Returns:
+            float: The atomic mass of the element in amu.
+        """
+        return self._atomic_mass_number
+
     def __getattr__(self, item: str) -> Any:
         """Key access to available element data.
 
@@ -454,7 +477,7 @@ def ground_state_term_symbol(self):
         return J_sorted_terms[-1][0]
 
     def __eq__(self, other: object) -> bool:
-        return isinstance(self, Element) and isinstance(other, Element) and self.Z == other.Z
+        return isinstance(self, Element) and isinstance(other, Element) and self.Z == other.Z and self.A == other.A
 
     def __hash__(self) -> int:
         return self.Z
@@ -482,17 +505,19 @@ def __lt__(self, other):
         return self.symbol < other.symbol
 
     @staticmethod
-    def from_Z(Z: int) -> Element:
+    def from_Z(Z: int, A: int | None = None) -> Element:
         """Get an element from an atomic number.
 
         Args:
-            Z (int): Atomic number
+            Z (int): Atomic number (number of protons)
+            A (int or None) : Atomic mass number (number of protons + neutrons)
 
         Returns:
             Element with atomic number Z.
         """
         for sym, data in _pt_data.items():
-            if data["Atomic no"] == Z:
+            amn = data.get("Atomic mass no") if A else None
+            if data["Atomic no"] == Z and amn == A:
                 return Element(sym)
         raise ValueError(f"Unexpected atomic number {Z=}")
 
@@ -506,6 +531,9 @@ def from_name(name: str) -> Element:
         Returns:
             Element with the name 'name'
         """
+        # to accommodate the British Enlgish speaking world
+        GBE_to_AmE = {"aluminium": "aluminum", "caesium": "cesium"}
+        name = GBE_to_AmE.get(name.lower(), name)
         for sym, data in _pt_data.items():
             if data["Name"] == name.capitalize():
                 return Element(sym)
@@ -765,6 +793,8 @@ class Element(ElementBase):
     # necessary to preserve backwards compatibility with a time when Element is
     # a regular object that is constructed with Element(symbol).
     H = "H"
+    D = "D"
+    T = "T"
     He = "He"
     Li = "Li"
     Be = "Be"
@@ -1407,6 +1437,8 @@ def get_el_sp(obj: int | SpeciesLike) -> Element | Species | DummySpecies:
             that can be determined.
     """
     if isinstance(obj, (Element, Species, DummySpecies)):
+        if hasattr(obj, "_is_named_isotope") and obj._is_named_isotope:
+            return Element(obj.name) if isinstance(obj, Element) else Species(obj.name)
         return obj
 
     try:

pymatgen/core/periodic_table.json

@@ -334,7 +334,6 @@ def __init__(
                     species = Composition({get_el_sp(species): 1})  # type: ignore
                 except TypeError:
                     species = Composition(species)
-
             total_occu = species.num_atoms
             if total_occu > 1 + Composition.amount_tolerance:
                 raise ValueError("Species occupancies sum to more than 1!")