Add zero-electron species to gaussian dataset (#138)

* Add zero-electron species to gaussian dataset

* Add proton msg file under gaussian dataset test files.

* Add special case for intensive properties with zero-electron components

---------

Co-authored-by: Gabriela Sanchez-Diaz <[email protected]>

.gitignore

@@ -146,8 +146,10 @@ cython_debug/
 .spyproject
 .vscode/
 
-# Raw data files
-atomdb/datasets/*/db/
+# Data files
+#*.msg
+repo_data.txt
+atomdb/datasets/*/raw/
 
 # Generated documentation
 docs/source/api/

atomdb/__main__.py

@@ -19,7 +19,6 @@
 
 from atomdb import compile_species, load
 
-
 # Initialize command line argument parser
 #
 parser = ArgumentParser(prog="atomdb", description="Compile or query an AtomDB entry.")
@@ -58,7 +57,7 @@
 
     args = parser.parse_args()
 
-    if args.compile:
+    if args.compile_species:
 
         compile_species(args.elem, args.charge, args.mult, args.exc, args.dataset)
 

atomdb/datasets/gaussian/run.py

@@ -18,35 +18,25 @@
 import os
 
 import numpy as np
-
-from gbasis.wrappers import from_iodata
-
 from gbasis.evals.density import evaluate_density as eval_dens
-from gbasis.evals.density import evaluate_density_gradient
-from gbasis.evals.density import evaluate_density_hessian
 from gbasis.evals.density import evaluate_posdef_kinetic_energy_density as eval_pd_ked
-from gbasis.evals.eval_deriv import evaluate_deriv_basis
 from gbasis.evals.eval import evaluate_basis
-
+from gbasis.wrappers import from_iodata
+from grid.atomgrid import AtomGrid
 from grid.onedgrid import UniformInteger
 from grid.rtransform import ExpRTransform
-from grid.atomgrid import AtomGrid
-
 from iodata import load_one
 
 import atomdb
-
-from atomdb.periodic import Element
 from atomdb.datasets.tools import (
+    eval_orb_ked,
     eval_orbs_density,
     eval_orbs_radial_d_density,
     eval_orbs_radial_dd_density,
-    eval_orb_ked,
     eval_radial_d_density,
     eval_radial_dd_density,
-    eval_orb_ked,
 )
-
+from atomdb.periodic import Element
 
 __all__ = [
     "run",
@@ -123,7 +113,11 @@ def run(elem, charge, mult, nexc, dataset, datapath):
 
     # Load data from fchk
     obasis_name = "def2-svpd"
-    data = _load_fchk(atnum, elem, nelec, mult, obasis_name, datapath)
+    if nelec == 0:
+        # For zero-electron case. use 1-electron case as a base
+        data = _load_fchk(atnum, elem, 1, 2, obasis_name, datapath)
+    else:
+        data = _load_fchk(atnum, elem, nelec, mult, obasis_name, datapath)
 
     # Unrestricted Hartree-Fock SCF results
     energy = data.energy
@@ -137,7 +131,7 @@ def run(elem, charge, mult, nexc, dataset, datapath):
     coeffs_b = mo_coeffs[:, norba:]
 
     # check for inconsistencies in filenames
-    if not np.allclose(np.array([n_up, n_dn]), np.array([sum(occs_up), sum(occs_dn)])):
+    if nelec != 0 and not np.allclose([n_up, n_dn], [sum(occs_up), sum(occs_dn)]):
         raise ValueError(f"Inconsistent data in fchk file for N: {atnum}, M: {mult} CH: {charge}")
 
     # Prepare data for computing Species properties
@@ -234,13 +228,38 @@ def run(elem, charge, mult, nexc, dataset, datapath):
 
     # Conceptual-DFT properties (WIP)
     # NOTE: Only the alpha component of the MOs is used bellow
+    # NOTE: Handle zero-electron case here
     mo_energy_occ_up = mo_e_up[:n_up]
     mo_energy_virt_up = mo_e_up[n_up:]
-    ip = -mo_energy_occ_up[-1]  # - energy_HOMO_alpha
-    ea = -mo_energy_virt_up[0]  # - energy_LUMO_alpha
+    ip = -mo_energy_occ_up[-1] if nelec != 0 else 0 # - energy_HOMO_alpha
+    ea = -mo_energy_virt_up[0] if nelec != 0 else 0 # - energy_LUMO_alpha
     mu = None
     eta = None
 
+    # Set appropriate values to zero for zero-electron case
+    if nelec == 0:
+        energy=0.0,
+        mo_e_up[...] = 0
+        mo_e_dn[...] = 0
+        occs_up[...] = 0
+        occs_dn[...] = 0
+        # Density
+        orb_dens_avg_up[...] = 0
+        orb_dens_avg_dn[...] = 0
+        dens_avg_tot[...] = 0
+        # Density gradient
+        orb_d_dens_avg_up[...] = 0
+        orb_d_dens_avg_dn[...] = 0
+        d_dens_avg_tot[...] = 0
+        # Density laplacian
+        orb_dd_dens_avg_up[...] = 0
+        orb_dd_dens_avg_dn[...] = 0
+        dd_dens_avg_tot[...] = 0
+        # KED
+        orb_ked_avg_up[...] = 0
+        orb_ked_avg_dn[...] = 0
+        ked_avg_tot[...] = 0
+
     # Return Species instance
     fields = dict(
         elem=elem,
@@ -261,6 +280,7 @@ def run(elem, charge, mult, nexc, dataset, datapath):
         mo_occs_a=occs_up,
         mo_occs_b=occs_dn,
         ip=ip,
+        # ea=ea,
         mu=mu,
         eta=eta,
         rs=rs,

atomdb/promolecule.py

@@ -248,7 +248,7 @@ def nspin(self, p=1):
         def f(atom):
             return atom.nspin
 
-        return _intensive_property(self.atoms, self.coeffs, f, p=1)
+        return _intensive_property_exclude_zero(self.atoms, self.coeffs, f, p=1)
 
     def mult(self, p=1):
         r"""Compute the multiplicity of the promolecule.
@@ -701,7 +701,7 @@ def make_promolecule(
             promol._extend(*(min(good_combs, key=itemgetter(0))[1:]))
         else:
             raise ValueError(
-                "Unable to construct species with non-integer" "charge/spin from database entries"
+                "Unable to construct species with non-integer charge/spin from database entries"
             )
 
     # Return Promolecule instance
@@ -787,6 +787,32 @@ def _intensive_property(atoms, coeffs, f, p=1):
     )
 
 
+def _intensive_property_exclude_zero(atoms, coeffs, f, p=1):
+    r"""Helper function for computing intensive properties (excluding zero-electron proatoms).
+
+    Parameters
+    ----------
+    atoms: list of Species
+        Species instances.
+    coeffs: np.ndarray((N,), dtype=float)
+        Coefficients of each species.
+    f: callable
+        Property function.
+    p: int, default=1 (linear mean)
+        Type of mean used in the computation.
+
+    Returns
+    -------
+    prop: float
+        Intensive property.
+    """
+    # P-mean of each atom's property value
+    return (
+            sum(coeff * f(atom) ** p for atom, coeff in zip(atoms, coeffs) if atom.nelec != 0)
+                / sum(coeff for atom, coeff in zip(atoms, coeffs) if atom.nelec != 0) ** (1 / p)
+    )
+
+
 def _radial_vector_outer_triu(radii):
     r"""Evaluate the outer products of a set of radial unit vectors.
 

atomdb/test/test_promolecule.py

@@ -149,6 +149,22 @@
         },
         id="Be floating point charge/integer mult (neg)",
     ),
+    pytest.param(
+        {
+            "atnums": [3, 1],
+            "charges": [0, 1],
+            "mults": [2, 1],
+            "coords": np.asarray(
+                [
+                    [0.0, 0.0, 0.0],
+                    [0.0, 0.0, 1.0],
+                ],
+                dtype=float,
+            ),
+            "dataset": "gaussian",
+        },
+        id="LiH with zero electrons on the hydrogen center",
+    ),
 ]