tweak test infrastructure

adcc/gradients/TwoParticleDensityMatrix.py

@@ -39,8 +39,6 @@ def __init__(self, spaces):
             self.mospaces = spaces.mospaces
         else:
             self.mospaces = spaces
-#        if self.mospaces.has_core_occupied_space:
-#            raise NotImplementedError("TPDM not implemented for CVS.")
         # Set reference_state if possible
         if isinstance(spaces, libadcc.ReferenceState):
             self.reference_state = spaces

adcc/gradients/amplitude_response.py

@@ -109,7 +109,7 @@ def ampl_relaxed_dms_cvs_adc1(exci):
     g1a.cc = - 1.0 * einsum("Ia,Ja->IJ", u.ph, u.ph)
     g1a.vv = + 1.0 * einsum("Ia,Ib->ab", u.ph, u.ph)
 
-    # Pre-requisites for the OC block of the
+    # Prerequisites for the OC block of the
     # orbital response Lagrange multipliers:
     fc = hf.fock(b.cc).diagonal()
     fo = hf.fock(b.oo).diagonal()

adcc/gradients/test_functionality_gradients.py

@@ -22,6 +22,7 @@
 ## ---------------------------------------------------------------------
 import unittest
 import itertools
+
 import adcc
 import adcc.backends
 
@@ -36,9 +37,11 @@
 
 backends = [b for b in adcc.backends.available()
             if b not in ["molsturm", "veloxchem"]]
-molecules = ["h2o"]
-basissets = ["sto3g", "ccpvdz"]
-methods = ["mp2", "adc1", "adc2"]
+
+molecules = gradient_data["molecules"]
+basissets = gradient_data["basissets"]
+methods = gradient_data["methods"]
+
 combinations = list(itertools.product(molecules, basissets, methods, backends))
 
 
@@ -50,19 +53,19 @@ def template_nuclear_gradient(self, molecule, basis, method, backend):
 
         energy_ref = grad_ref["energy"]
         grad_fdiff = grad_ref["gradient"]
+
+        kwargs = grad_ref["config"]
 
         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 = 5
-            state = adcc.run_adc(scfres, method=method,
-                                 n_singlets=10, conv_tol=1e-11)
-            for ee in state.excitations[:n_limit]:
+            state = adcc.run_adc(scfres, method=method, **kwargs)
+            for i, ee in enumerate(state.excitations):
                 grad = adcc.nuclear_gradient(ee)
                 assert_allclose(energy_ref[ee.index], ee.total_energy, atol=1e-10)
                 assert_allclose(
-                    grad_fdiff[ee.index], grad["Total"], atol=1e-7
+                    grad_fdiff[ee.index], grad["Total"], atol=1e-7, err_msg=f'State {i} wrong.'
                 )
         else:
             # MP2 gradients

adcc/testdata/dump_fdiff_gradient.py

@@ -20,6 +20,7 @@
 ## along with adcc. If not, see <http://www.gnu.org/licenses/>.
 ##
 ## ---------------------------------------------------------------------
+from collections import namedtuple
 import itertools
 import adcc
 import numpy as np
@@ -31,8 +32,10 @@
 from static_data import xyz
 
 
-prefactors_5p = np.array([1.0, -8.0, 8.0, -1.0]) / 12.0
-multipliers_5p = [-2, -1, 1, 2]
+# prefactors_5p = np.array([1.0, -8.0, 8.0, -1.0]) / 12.0
+# multipliers_5p = [-2, -1, 1, 2]
+prefactors_9p = [1./280., -4./105., 1./5., -4./5., 4./5., -1./5., 4./105., -1./280.]
+multipliers_9p = [-4., -3., -2., -1., 1., 2., 3., 4.]
 coords_label = ["x", "y", "z"]
 
 
@@ -47,46 +50,47 @@ def _molstring(elems, coords):
 
 def adc_energy(scfres, method, **kwargs):
     state = adcc.run_adc(method=method, data_or_matrix=scfres,
-                         output=None, **kwargs)
+                         output=None,
+                         **kwargs)
     return state.total_energy
 
 
 def mp_energy(scfres, method, **kwargs):
-    level = {
-        "mp2": 2,
-        "mp3": 3,
-    }
+    level = method[-1]
     refstate = adcc.ReferenceState(scfres)
-    return adcc.LazyMp(refstate).energy(level[method])
+    return adcc.LazyMp(refstate).energy(level)
 
 
-def fdiff_gradient(molstring, method, basis, step=1e-4, **kwargs):
-    m = gto.M(atom=molstring, unit='Bohr', basis=basis)
+def fdiff_gradient(molecule, method, basis, step=1e-4, **kwargs):
+    molstring = xyz[molecule.name]
+    m = gto.M(atom=molstring, unit='Bohr', basis=basis,
+              spin=molecule.multiplicity - 1, charge=molecule.charge)
     coords = m.atom_coords().copy()
     elements = m.elements.copy()
 
-    n_grads = kwargs.get("n_singlets", 1)
-    conv_tol = kwargs.get("conv_tol", 1e-10) / 10
+    conv_tol = kwargs.get("conv_tol", 1e-10) / 100
 
     # run unperturbed system
     scfres = adcc.backends.run_hf(
-        'pyscf', molstring, basis, conv_tol=conv_tol, conv_tol_grad=conv_tol
+        'pyscf', molstring, 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)
     else:
         en = mp_energy(scfres, method, **kwargs)
 
     natoms = len(elements)
-    grad = np.zeros((n_grads, natoms, 3))
+    grad = np.zeros((len(en), natoms, 3))
     at_c = list(itertools.product(range(natoms), range(3)))
     for i, c in tqdm(at_c):
-        for f, p in zip(multipliers_5p, prefactors_5p):
+        for f, p in zip(multipliers_9p, prefactors_9p):
             coords_p = coords.copy()
             coords_p[i, c] += f * step
             geom_p = _molstring(elements, coords_p)
             scfres = adcc.backends.run_hf(
-                'pyscf', geom_p, basis, conv_tol=conv_tol, conv_tol_grad=conv_tol
+                'pyscf', geom_p, basis, conv_tol=conv_tol, conv_tol_grad=conv_tol,
+                charge=molecule.charge, multiplicity=molecule.multiplicity
             )
             if "adc" in method:
                 en_pert = adc_energy(scfres, method, **kwargs)
@@ -97,38 +101,60 @@ def fdiff_gradient(molstring, method, basis, step=1e-4, **kwargs):
 
 
 def main():
-    config_excited = {
-        "n_singlets": 5,
-    }
     basissets = [
         "sto3g",
         "ccpvdz",
     ]
     methods = [
-        "mp2",
+        # "mp2",
+        # "adc0",
         "adc1",
-        "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"
     ]
-    molecules = ["h2o", "hf", "formaldehyde"]
-    ret = {}
+    ret = {
+        "basissets": basissets,
+        "methods": methods,
+        "molecules": molecules,
+    }
+    config_excited = {"n_states": 3}
     for molecule in molecules:
-        ret[molecule] = {}
+        ret[molecule.name] = {}
         for basis in basissets:
-            ret[molecule][basis] = {}
+            ret[molecule.name][basis] = {}
             for method in methods:
                 kwargs = {
-                    "conv_tol": 1e-8,
+                    "conv_tol": 1e-10,
                 }
                 if "adc" in method:
                     kwargs.update(config_excited)
-                basename = f"{molecule}_{basis}_{method}"
+                basename = f"{molecule.name}_{basis}_{method}"
                 print(f"Evaluating finite difference gradient for {basename}.")
-                en, grad = fdiff_gradient(xyz[molecule], method, basis, **kwargs)
+
+                # HACK
+                core_orbitals = None
+                if "cvs" in method:
+                    core_orbitals = 1
+                kwargs["core_orbitals"] = core_orbitals
+                en, grad = fdiff_gradient(molecule, method, basis, **kwargs)
                 if isinstance(en, np.ndarray):
                     en = en.tolist()
                 cont = {
                     "energy": en,
                     "gradient": np.squeeze(grad).tolist(),
+                    "config": kwargs,
                 }
                 ret[molecule][basis][method] = cont
     with open("grad_dump.yml", "w") as yamlout: