try to finally make damn tests work :(

adcc/ExcitedStates.py

@@ -222,6 +222,8 @@ def __add__(self, other):
     @mark_excitation_property()
     def total_energy(self):
         # TODO: excitation_energy_uncorrected for PE-ADC with postSCF
+        if self.method.level == 0:
+            return self.excitation_energy + self.reference_state.energy_scf
         return self.excitation_energy + self.ground_state.energy(self.method.level)
 
     @cached_property

adcc/gradients/test_functionality_gradients.py

@@ -22,6 +22,7 @@
 ## ---------------------------------------------------------------------
 import unittest
 import itertools
+import numpy as np
 
 import adcc
 import adcc.backends
@@ -52,20 +53,23 @@ def template_nuclear_gradient(self, molecule, basis, method, backend):
         grad_ref = gradient_data[molecule][basis][method]
 
         energy_ref = grad_ref["energy"]
-        grad_fdiff = grad_ref["gradient"]
-
+        grad_fdiff = grad_ref["gradient"] 
         kwargs = grad_ref["config"]
+        conv_tol = kwargs["conv_tol"]
 
         scfres = cached_backend_hf(backend, molecule, basis, conv_tol=1e-13)
         if "adc" in method:
             # TODO: convergence needs to be very very tight...
             # so we want to make sure all vectors are tightly converged
+            n_limit = 2 # kwargs["n_singlets"]
+            kwargs["n_singlets"] = kwargs["n_singlets"] + 2
             state = adcc.run_adc(scfres, method=method, **kwargs)
-            for i, ee in enumerate(state.excitations):
+            for ee in state.excitations[:n_limit]:
                 grad = adcc.nuclear_gradient(ee)
-                assert_allclose(energy_ref[ee.index], ee.total_energy, atol=1e-10)
+                assert_allclose(energy_ref[ee.index], ee.total_energy, atol=conv_tol)
                 assert_allclose(
-                    grad_fdiff[ee.index], grad["Total"], atol=1e-7, err_msg=f'State {i} wrong.'
+                    grad_fdiff[ee.index], grad["Total"], atol=1e-7,
+                    err_msg=f'Gradient for state {ee.index} wrong.'
                 )
         else:
             # MP2 gradients

adcc/testdata/dump_fdiff_gradient.py

@@ -27,9 +27,12 @@
 import yaml
 from tqdm import tqdm
 
-from pyscf import gto
+from pyscf import gto, lib
 
-from static_data import xyz
+from static_data import molecules
+
+adcc.set_n_threads(8)
+lib.num_threads(8)
 
 
 # prefactors_5p = np.array([1.0, -8.0, 8.0, -1.0]) / 12.0
@@ -56,14 +59,13 @@ def adc_energy(scfres, method, **kwargs):
 
 
 def mp_energy(scfres, method, **kwargs):
-    level = method[-1]
+    level = int(method[-1])
     refstate = adcc.ReferenceState(scfres)
     return adcc.LazyMp(refstate).energy(level)
 
 
 def fdiff_gradient(molecule, method, basis, step=1e-4, **kwargs):
-    molstring = xyz[molecule.name]
-    m = gto.M(atom=molstring, unit='Bohr', basis=basis,
+    m = gto.M(atom=molecule.xyz, unit='Bohr', basis=basis,
               spin=molecule.multiplicity - 1, charge=molecule.charge)
     coords = m.atom_coords().copy()
     elements = m.elements.copy()
@@ -72,16 +74,18 @@ def fdiff_gradient(molecule, method, basis, step=1e-4, **kwargs):
 
     # run unperturbed system
     scfres = adcc.backends.run_hf(
-        'pyscf', molstring, basis, conv_tol=conv_tol, conv_tol_grad=conv_tol,
+        'pyscf', molecule.xyz, basis, conv_tol=conv_tol, conv_tol_grad=conv_tol,
         charge=molecule.charge, multiplicity=molecule.multiplicity
     )
     if "adc" in method:
         en = adc_energy(scfres, method, **kwargs)
+        ngrad = len(en)
     else:
         en = mp_energy(scfres, method, **kwargs)
+        ngrad = 1
 
     natoms = len(elements)
-    grad = np.zeros((len(en), natoms, 3))
+    grad = np.zeros((ngrad, natoms, 3))
     at_c = list(itertools.product(range(natoms), range(3)))
     for i, c in tqdm(at_c):
         for f, p in zip(multipliers_9p, prefactors_9p):
@@ -106,47 +110,43 @@ def main():
         "ccpvdz",
     ]
     methods = [
-        # "mp2",
-        # "adc0",
+        "mp2",
+        "adc0",
         "adc1",
-        # "adc2",
-        # "cvs-adc0",
-        # "cvs-adc1",
-        # "cvs-adc2",
+        "adc2",
+        "cvs-adc0",
+        "cvs-adc1",
+        "cvs-adc2",
         # "cvs-adc2x",
     ]
-    Molecule = namedtuple("Molecule", ["name", "charge", "multiplicity"])
-    molecules = [
-        # Molecule("h2o", 0, 1),
-        # "h2s",
-        Molecule("cn", 0, 2),
-        # "ch2nh2",
-        # "hf",
-        # "formaldehyde"
+    molnames = [
+        "h2o",
+        "h2s",
     ]
+    mols = [m for m in molecules if m.name in molnames]
     ret = {
         "basissets": basissets,
         "methods": methods,
-        "molecules": molecules,
+        "molecules": molnames,
     }
-    config_excited = {"n_states": 3}
-    for molecule in molecules:
-        ret[molecule.name] = {}
+    config_excited = {"n_singlets": 3}
+    for molecule in mols:
+        molname = molecule.name
+        ret[molname] = {}
         for basis in basissets:
-            ret[molecule.name][basis] = {}
+            ret[molname][basis] = {}
             for method in methods:
                 kwargs = {
                     "conv_tol": 1e-10,
                 }
                 if "adc" in method:
                     kwargs.update(config_excited)
-                basename = f"{molecule.name}_{basis}_{method}"
+                basename = f"{molname}_{basis}_{method}"
                 print(f"Evaluating finite difference gradient for {basename}.")
 
-                # HACK
                 core_orbitals = None
                 if "cvs" in method:
-                    core_orbitals = 1
+                    core_orbitals = molecule.core_orbitals
                 kwargs["core_orbitals"] = core_orbitals
                 en, grad = fdiff_gradient(molecule, method, basis, **kwargs)
                 if isinstance(en, np.ndarray):
@@ -156,7 +156,7 @@ def main():
                     "gradient": np.squeeze(grad).tolist(),
                     "config": kwargs,
                 }
-                ret[molecule][basis][method] = cont
+                ret[molname][basis][method] = cont
     with open("grad_dump.yml", "w") as yamlout:
         yaml.safe_dump(ret, yamlout)