add SCI data to molecular data (#293)

python/ffsim/molecular_data.py

@@ -26,6 +26,7 @@
 import pyscf
 import pyscf.cc
 import pyscf.ci
+import pyscf.fci
 import pyscf.mcscf
 import pyscf.mp
 import pyscf.symm
@@ -68,6 +69,9 @@ class MolecularData:
             CCSD t2 amplitudes.
         cisd_energy (float | None): The CISD energy.
         cisd_vec (np.ndarray | None): The CISD state vector.
+        sci_energy (float | None): The SCI energy.
+        sci_vec (tuple[np.ndarray, np.ndarray, np.ndarray] | None): The SCI state
+            vector coefficients, spin alpha strings, and spin beta strings.
         fci_energy (float | None): The FCI energy.
         fci_vec (np.ndarray | None): The FCI state vector.
         dipole_integrals (np.ndarray | None): The dipole integrals.
@@ -104,6 +108,9 @@ class MolecularData:
     # CISD data
     cisd_energy: float | None = None
     cisd_vec: np.ndarray | None = None
+    # SCI data
+    sci_energy: float | None = None
+    sci_vec: tuple[np.ndarray, np.ndarray, np.ndarray] | None = None
     # FCI data
     fci_energy: float | None = None
     fci_vec: np.ndarray | None = None
@@ -220,6 +227,19 @@ def run_cisd(self, *, store_cisd_vec: bool = False) -> None:
         if store_cisd_vec:
             self.cisd_vec = cisd_vec
 
+    def run_sci(self, *, store_sci_vec: bool = False) -> None:
+        """Run SCI and store results."""
+        sci = pyscf.fci.SCI(self.scf.run())
+        sci_energy, sci_vec = sci.kernel(
+            self.one_body_integrals,
+            self.two_body_integrals,
+            norb=self.norb,
+            nelec=self.nelec,
+        )
+        self.sci_energy = sci_energy + self.core_energy
+        if store_sci_vec:
+            self.sci_vec = (sci_vec, *sci_vec._strs)
+
     def run_fci(self, *, store_fci_vec: bool = False) -> None:
         """Run FCI and store results."""
         cas = pyscf.mcscf.CASCI(self.scf.run(), ncas=self.norb, nelecas=self.nelec)
@@ -343,6 +363,8 @@ def as_array_tuple_or_none(val):
             ccsd_t2=arrays_func(data.get("ccsd_t2")),
             cisd_energy=data.get("cisd_energy"),
             cisd_vec=as_array_or_none(data.get("cisd_vec")),
+            sci_energy=data.get("sci_energy"),
+            sci_vec=as_array_tuple_or_none(data.get("sci_vec")),
             fci_energy=data.get("fci_energy"),
             fci_vec=as_array_or_none(data.get("fci_vec")),
             dipole_integrals=as_array_or_none(data.get("dipole_integrals")),

tests/python/molecular_data_test.py

@@ -31,6 +31,7 @@ def _assert_mol_data_equal(
             "np.ndarray | None",
             "np.ndarray | tuple[np.ndarray, np.ndarray] | None",
             "np.ndarray | tuple[np.ndarray, np.ndarray, np.ndarray] | None",
+            "tuple[np.ndarray, np.ndarray, np.ndarray] | None",
         ]:
             if actual is not None:
                 if isinstance(actual, tuple):
@@ -85,11 +86,13 @@ def test_molecular_data_run_methods():
     mol_data.run_mp2()
     mol_data.run_ccsd()
     mol_data.run_cisd()
+    mol_data.run_sci()
     mol_data.run_fci()
 
     np.testing.assert_allclose(mol_data.mp2_energy, -108.58852784026)
     np.testing.assert_allclose(mol_data.ccsd_energy, -108.5933309085008)
     np.testing.assert_allclose(mol_data.cisd_energy, -108.5878344909782)
+    np.testing.assert_allclose(mol_data.sci_energy, -108.59598682615388)
     np.testing.assert_allclose(mol_data.fci_energy, -108.595987350986)
 
 
@@ -108,6 +111,7 @@ def test_json_closed_shell(tmp_path: pathlib.Path):
     mol_data.run_mp2(store_t2=True)
     mol_data.run_ccsd(store_t1=True, store_t2=True)
     mol_data.run_cisd(store_cisd_vec=True)
+    mol_data.run_sci(store_sci_vec=True)
     mol_data.run_fci(store_fci_vec=True)
 
     for compression in [None, "gzip", "bz2", "lzma"]:
@@ -132,6 +136,7 @@ def test_json_open_shell(tmp_path: pathlib.Path):
     mol_data.run_mp2(store_t2=True)
     mol_data.run_ccsd(store_t1=True, store_t2=True)
     mol_data.run_cisd(store_cisd_vec=True)
+    mol_data.run_sci(store_sci_vec=True)
     mol_data.run_fci(store_fci_vec=True)
 
     for compression in [None, "gzip", "bz2", "lzma"]: