Corrections to NIST dataset

atomdb/datasets/nist/run.py

@@ -19,8 +19,6 @@
 
 import numpy as np
 
-from scipy import constants
-
 import h5py as h5
 
 import csv
@@ -91,6 +89,57 @@ def load_nist_spectra_data(atnum, nelec, datafile):
     return output
 
 
+def get_energy(atnum, nelec, ip, datafile):
+    """Retrieve the ground state energy from NIST spectra database.
+
+    Parameters
+    ----------
+    atnum : int
+        Atomic number.
+    nelec : int
+        Number of electrons.
+    ip : float
+        Ionization potential of the anion (in Ha).
+    datafile : str
+        HDF5 file containing the NIST atomic spectra data (file database_beta_1.3.0.h5).
+
+    Returns
+    -------
+    float or None
+        Ground state energy of the species (in Ha). Returns None if data is unavailable.
+
+    Notes
+    ----- 
+    - Dianion species (charge = -2) return None.
+    - Anion energy (charge = -1) is computed as:  
+        :math:`E_{anion} = E_{neutral} - IP_{anion}`
+    where :math:`IP_{anion}` is obtained from the conceptual-DFT data in the file data/c6cp04533b1.csv.
+    """
+    default_energy = None
+    charge = atnum - nelec
+    hdf5_path = os.path.join(MODULE_DATAPATH, datafile)
+
+    # Return None for dianions (charge = -2)
+    if charge == -2:
+        return default_energy
+
+    # Load NIST atomic spectra database data (contains neutral and cationic species).
+    # For anions, load data for the neutral species (nelec = atnum)
+    nelec = nelec if charge != -1 else atnum
+    spectra_data = load_nist_spectra_data(atnum, nelec, hdf5_path)
+    energies = spectra_data["energy"]
+
+    # Return energy (in Hartree) for neutral and cationic species (charge ≥ 0)
+    if charge >= 0:
+        return energies[0] / CMINV if len(energies) != 0 else default_energy
+
+    # Compute anion energy (charge = -1), ensuring ip is not zero or None
+    if charge == -1 and ip not in [None, 0]:
+        return (energies[0] / CMINV - ip) if len(energies) != 0 else default_energy
+
+    return default_energy
+
+
 def run(elem, charge, mult, nexc, dataset, datapath):
     r"""Parse NIST related data and compile the AtomDB database entry."""
     # Check arguments
@@ -126,23 +175,15 @@ def run(elem, charge, mult, nexc, dataset, datapath):
         polarizability = atom.pold
         dispersion = {"C6": atom.c6}
 
-    #
-    # Get the ground state energy from database_beta_1.3.0.h5.
-    #
-    # Set an energy default value since there is no data for anions in database_beta_1.3.0.h5.
-    energy = None
-    h5path = os.path.join(MODULE_DATAPATH, "database_beta_1.3.0.h5")
-    if charge >= 0:  # neutral or cationic species
-        spectra_data = load_nist_spectra_data(atnum, nelec, h5path)
-        energies = spectra_data["energy"]
-        # Convert energy to Hartree from cm^{-1} if available
-        energy = energies[0] * CMINV if len(energies) != 0 else energy
-
     # Get conceptual-DFT related properties from c6cp04533b1.csv
     # Locate where each table starts: search for "Element" columns
     csvpath = os.path.join(MODULE_DATAPATH, "c6cp04533b1.csv")
     data = list(csv.reader(open(csvpath, "r")))
     tabid = [i for i, row in enumerate(data) if "Element" in row]
+    # Check that the CSV file has not been modified and it has the expected number of tables (3)
+    if len(tabid) != 3:
+        raise ValueError(f"Unexpected conceptual-DFT CSV file format; expected 3 tables, got {len(tabid)}.")
+
     # Assign each conceptual-DFT data table to a variable.
     # Remove empty and header rows
     table_ips = data[tabid[0] : tabid[1]]
@@ -159,6 +200,9 @@ def run(elem, charge, mult, nexc, dataset, datapath):
     colid = table_etas[0].index(str(charge))
     eta = float(table_etas[atnum][colid]) * EV if len(table_etas[atnum][colid]) > 1 else None
 
+    # Get the ground state energy from NIST database_beta_1.3.0.h5.
+    energy = get_energy(atnum, nelec, ip, "database_beta_1.3.0.h5")
+
     # Return Species instance
     fields = dict(
         elem=elem,

atomdb/species.py

@@ -77,9 +77,13 @@ def wrapper(self):
             "at_radius": "at_radius",
             "polarizability": "pold",
             "dispersion_c6": "c6",
+            "atmass": "mass",
         }
 
+        if name == "atmass":
+            return getattr(Element(self._data.elem), namemap[name])
         if name in namemap:
+            # Only return Element property if neutral, otherwise None
             charge = self._data.atnum - self._data.nelec
             return getattr(Element(self._data.elem), namemap[name]) if charge == 0 else None
 
@@ -476,7 +480,7 @@ def atmass(self):
             https://github.com/theochem/AtomDB/blob/master/atomdb/data/data_info.csv
 
         """
-        return Element(self._data.elem).mass
+        pass
 
     @scalar
     def cov_radius(self):

atomdb/test/test_nist.py

@@ -36,7 +36,7 @@
             "at_radius": {"wc": 0.529192875 * ANGSTROM, "cr": 0.53 * ANGSTROM},
             "polarizability": {"crc": 0.666793 * ANGSTROM**3, "chu": 4.5, "sn": 4.50711},
             "dispersion_c6": {"chu": 6.499026705},
-            "energy": -109678.77174307 * CMINV,
+            "energy": -109678.77174307 / CMINV,
             "ip": 13.598443 * EV,
             "mu": -7.18 * EV,
             "eta": 12.84 * EV,
@@ -70,7 +70,7 @@
             "at_radius": {"wc": 0.62 * ANGSTROM, "cr": 0.67 * ANGSTROM},
             "polarizability": {"crc": 1.76 * ANGSTROM**3, "chu": 12.0, "sn": 11.3},
             "dispersion_c6": {"chu": 46.6},
-            "energy": -8308396.1899999995 * CMINV,
+            "energy": -8308396.1899999995 / CMINV,
             "ip": 11.2603 * EV,
             "mu": -6.26 * EV,
             "eta": 10 * EV,
@@ -92,7 +92,7 @@
             "at_radius": None,
             "polarizability": None,
             "dispersion_c6": None,
-            "energy": -8217575.77 * CMINV,
+            "energy": -8217575.77 / CMINV,
             "ip": 24.3833 * EV,
             "mu": -17.82 * EV,
             "eta": 13.12 * EV,
@@ -114,13 +114,35 @@
             "at_radius": None,
             "polarizability": None,
             "dispersion_c6": None,
-            "energy": None,
+            "energy": (-8308396.1899999995 / CMINV) - (1.262118 * EV),
             "ip": 1.262118 * EV,
             "mu": 2.08 * EV,
             "eta": 8.02 * EV,
         },
         id="C anion",
     ),
+    pytest.param(
+        {
+            "dataset": "nist",
+            "elem": "C",
+            "atnum": 6,
+            "obasis_name": None,
+            "nelec": 8,
+            "nspin": 2,
+            "nexc": 0,
+            "atmass": {"stb": 12.0106 * AMU, "nist": 12 * AMU},
+            "cov_radius": None,
+            "vdw_radius": None,
+            "at_radius": None,
+            "polarizability": None,
+            "dispersion_c6": None,
+            "energy": None,
+            "ip": None,
+            "mu": 8.41 * EV,
+            "eta": 0 * EV,
+        },
+        id="C dianion",
+    ),
 ]
 
 
@@ -136,8 +158,4 @@ def test_nist_data(case):
 
     for attr, value in case.items():
         data_value = getattr(sp, attr)
-        # try:
-        #     data_value = getattr(sp, attr)
-        # except AttributeError:
-        #     data_value = getattr(sp._data, attr)
-        assert data_value == value, f"{elem} {attr} is not as expected."
+        assert data_value == pytest.approx(value), f"{elem} {attr} is not as expected."