Improve CIF checking, support for isotopes, and correct handling of new VASP 6.4.2 POSCAR format incl. slashes in header (#3542)

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

* Fix LAMMPS guess_element

* Refactor CifAssessor to be class method on CifParser + refactor assoc tests

* Fix failing test by adding DummySpecies.A attr

* Fix oxidation state parsing for isotopes

* Revert element hash to be Z if not an isotope, and add new distinct hashes for isotopes

* fix types in doc strings and tweak var names

* fix DummySpecies.A attr, refactor __eq__

* CifParser document defaults for check_cif, cif_assessor_tol, fix doc str format

* rename method CifParser.assess->check

* check D+ oxi state is correctly parsed

* tighten atomic mass test on H isotopes

* del debug prints

* merge TestCifParserAssess into TestCifIO

* Update docstring for DummySpecies.Z and DummySpecies.A; allow DummySpecies.A to return either int or None

* Revert format of dev_scripts/periodic_table.yaml to minimize diff

---------

Co-authored-by: Janosh Riebesell <[email protected]>

Author: esoteric-ephemera

dev_scripts/periodic_table.yaml

@@ -5435,3 +5435,26 @@ Zr:
   Vickers hardness: 903 MN m<sup>-2</sup>
   X: 1.33
   Youngs modulus: 68 GPa
+D:
+  Atomic no: 1
+  Atomic mass: 2.013553212712
+  Atomic mass no: 2
+  Common oxidation states: [-1, 1]
+  Is named isotope: true
+  Name: Deuterium
+  Oxidation states: [-1, -1]
+  Shannon radii:
+    1:
+      II:
+        ? ''
+        : {crystal_radius: 0.04, ionic_radius: -0.1}
+  Electron affinity: 0.754674
+T:
+  Atomic no: 1
+  Atomic mass: 3.0155007134
+  Atomic mass no: 3
+  Common oxidation states: [-1, -1]
+  Is named isotope: true
+  Name: Tritium
+  Oxidation states: [-1, -1]
+  Electron affinity:

dev_scripts/update_pt_data.py

@@ -19,22 +19,11 @@
 
 from pymatgen.core import Element, get_el_sp
 
-
-def test_yaml():
-    with open("periodic_table.yaml") as file:
-        data = yaml.load(file)
-        print(data)
-
-
-def test_json():
-    with open("periodic_table.json") as file:
-        data = json.load(file)
-        print(data)
+ptable_yaml_path = "periodic_table.yaml"
 
 
 def parse_oxi_state():
-    with open("periodic_table.yaml") as file:
-        data = yaml.load(file)
+    data = loadfn(ptable_yaml_path)
     with open("oxidation_states.txt") as file:
         oxi_data = file.read()
     oxi_data = re.sub("[\n\r]", "", oxi_data)
@@ -72,10 +61,9 @@ def parse_oxi_state():
 
 
 def parse_ionic_radii():
-    with open("periodic_table.yaml") as f:
-        data = yaml.load(f)
-    with open("ionic_radii.csv") as f:
-        radii_data = f.read()
+    data = loadfn(ptable_yaml_path)
+    with open("ionic_radii.csv") as file:
+        radii_data = file.read()
     radii_data = radii_data.split("\r")
     header = radii_data[0].split(",")
     for idx in range(1, len(radii_data)):
@@ -99,15 +87,14 @@ def parse_ionic_radii():
                 data[el]["Ionic_radii"] = ionic_radii
         else:
             print(el)
-    with open("periodic_table2.yaml", "w") as f:
-        yaml.dump(data, f)
+    with open("periodic_table2.yaml", "w") as file:
+        yaml.dump(data, file)
 
 
 def parse_radii():
-    with open("periodic_table.yaml") as f:
-        data = yaml.load(f)
-    with open("radii.csv") as f:
-        radii_data = f.read()
+    data = loadfn(ptable_yaml_path)
+    with open("radii.csv") as file:
+        radii_data = file.read()
     radii_data = radii_data.split("\r")
 
     for line in radii_data:
@@ -134,15 +121,14 @@ def parse_radii():
             data[el]["Van der waals radius"] = vdw_radii
         else:
             print(el)
-    with open("periodic_table2.yaml", "w") as f:
-        yaml.dump(data, f)
-    with open("periodic_table.json", "w") as f:
-        json.dump(data, f)
+    with open("periodic_table2.yaml", "w") as file:
+        yaml.dump(data, file)
+    with open("../pymatgen/core/periodic_table.json", "w") as file:
+        json.dump(data, file)
 
 
 def update_ionic_radii():
-    with open("periodic_table.yaml") as f:
-        data = yaml.load(f)
+    data = loadfn(ptable_yaml_path)
 
     for d in data.values():
         if "Ionic_radii" in d:
@@ -154,15 +140,14 @@ def update_ionic_radii():
         if "Ionic_radii_ls" in d:
             d["Ionic radii ls"] = {k: v / 100 for k, v in d["Ionic_radii_ls"].items()}
             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:
-        json.dump(data, f)
+    with open("periodic_table2.yaml", "w") as file:
+        yaml.dump(data, file)
+    with open("../pymatgen/core/periodic_table.json", "w") as file:
+        json.dump(data, file)
 
 
 def parse_shannon_radii():
-    with open("periodic_table.yaml") as f:
-        data = yaml.load(f)
+    data = loadfn(ptable_yaml_path)
 
     from openpyxl import load_workbook
 
@@ -194,22 +179,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:
-        json.dump(data, f)
+    dumpfn(data, ptable_yaml_path)
+    with open("../pymatgen/core/periodic_table.json", "w") as file:
+        json.dump(data, file)
 
 
 def gen_periodic_table():
-    with open("periodic_table.yaml") as f:
-        data = yaml.load(f)
+    data = loadfn(ptable_yaml_path)
 
-    with open("periodic_table.json", "w") as f:
-        json.dump(data, f)
+    with open("../pymatgen/core/periodic_table.json", "w") as file:
+        json.dump(data, file)
 
 
 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,26 +248,38 @@ 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
+        temp_str = re.sub(r"[\s\(\)]", "", r[3].strip("()[]"))
+        # hyphen-like characters used that can't be parsed by .float
+        bytes_rep = temp_str.encode("unicode_escape").replace(b"\\u2212", b"-")
+        ea[r[2]] = float(bytes_rep.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)
+    for key, val in pt.items():
+        val["Electron affinity"] = ea.get(Element(key).long_name)
     dumpfn(pt, "../pymatgen/core/periodic_table.json")
 
 
 def add_ionization_energies():
     """Update the periodic table data file with ground level and ionization energies from NIST."""
 
-    with open("NIST Atomic Ionization Energies Output.html") as f:
-        soup = BeautifulSoup(f.read(), "html.parser")
-    for t in soup.find_all("table"):
-        if "Hydrogen" in t.text:
+    with open("NIST Atomic Ionization Energies Output.html") as file:
+        soup = BeautifulSoup(file.read(), "html.parser")
+    for table in soup.find_all("table"):
+        if "Hydrogen" in table.text:
             break
     data = collections.defaultdict(list)
-    for tr in t.find_all("tr"):
-        row = [td.get_text().strip() for td in tr.find_all("td")]
+    for row in table.find_all("tr"):
+        row = [td.get_text().strip() for td in row.find_all("td")]
         if row:
             Z = int(row[0])
             val = re.sub(r"\s", "", row[8].strip("()[]"))
@@ -294,9 +289,9 @@ def add_ionization_energies():
     print(data[51])
     assert set(data).issuperset(range(1, 93))  # Ensure that we have data for up to U.
     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, [])
+    for key, val in pt.items():
+        del val["Ionization energy"]
+        val["Ionization energies"] = data.get(Element(key).long_name, [])
     dumpfn(pt, "../pymatgen/core/periodic_table.json")
 
 

pymatgen/analysis/diffraction/core.py

@@ -88,7 +88,7 @@ def get_plot(
             two_theta_range (tuple[float, float]): Range of two_thetas to calculate in degrees.
                 Defaults to (0, 90). Set to None if you want all diffracted beams within the limiting
                 sphere of radius 2 / wavelength.
-            annotate_peaks (str or None): Whether and how to annotate the peaks
+            annotate_peaks (str | None): Whether and how to annotate the peaks
                 with hkl indices. Default is 'compact', i.e. show short
                 version (oriented vertically), e.g. 100. If 'full', show
                 long version, e.g. (1, 0, 0). If None, do not show anything.
@@ -162,7 +162,7 @@ def show_plot(self, structure: Structure, **kwargs):
                 two_thetas to calculate in degrees. Defaults to (0, 90). Set to
                 None if you want all diffracted beams within the limiting
                 sphere of radius 2 / wavelength.
-            annotate_peaks (str or None): Whether and how to annotate the peaks
+            annotate_peaks (str | None): Whether and how to annotate the peaks
                 with hkl indices. Default is 'compact', i.e. show short
                 version (oriented vertically), e.g. 100. If 'full', show
                 long version, e.g. (1, 0, 0). If None, do not show anything.
@@ -180,7 +180,7 @@ def plot_structures(self, structures, fontsize=6, **kwargs):
                 two_thetas to calculate in degrees. Defaults to (0, 90). Set to
                 None if you want all diffracted beams within the limiting
                 sphere of radius 2 / wavelength.
-            annotate_peaks (str or None): Whether and how to annotate the peaks
+            annotate_peaks (str | None): Whether and how to annotate the peaks
                 with hkl indices. Default is 'compact', i.e. show short
                 version (oriented vertically), e.g. 100. If 'full', show
                 long version, e.g. (1, 0, 0). If None, do not show anything.

pymatgen/analysis/structure_matcher.py

@@ -1092,9 +1092,9 @@ def get_transformation(self, struct1, struct2):
             struct2 (Structure): Structure to transform.
 
         Returns:
-            supercell (numpy.ndarray(3, 3)): supercell matrix
-            vector (numpy.ndarray(3)): fractional translation vector
-            mapping (list(int or None)):
+            supercell (np.array(3, 3)): supercell matrix
+            vector (np.array(3)): fractional translation vector
+            mapping (list[int | None]):
                 The first len(struct1) items of the mapping vector are the
                 indices of struct1's corresponding sites in struct2 (or None
                 if there is no corresponding site), and the other items are