fix cvs-adc2x gradients

Author: maxscheurer

adcc/gradients/TwoParticleDensityMatrix.py

@@ -58,11 +58,12 @@ def __init__(self, spaces):
         ]
         if self.mospaces.has_core_occupied_space:
             self.blocks += [
-                b.cccc, b.ococ, b.oooo,  b.cvcv,
+                b.cccc, b.ococ, b.cvcv,
                 b.ocov, b.cccv, b.cocv, b.ocoo,
                 b.ccco, b.occv, b.ccvv, b.ocvv,
-                b.vvvv,
             ]
+        # make sure we didn't add any block twice!
+        assert len(list(set(self.blocks))) == len(self.blocks)
         self._tensors = {}
 
     @property
@@ -223,7 +224,7 @@ def to_ao_basis(self, refstate_or_coefficients=None):
             if not hasattr(self, "reference_state"):
                 raise ValueError("Argument reference_state is required if no "
                                  "reference_state is stored in the "
-                                 "OneParticleOperator")
+                                 "TwoParticleDensityMatrix")
             return self.__transform_to_ao(self.reference_state)
         else:
             raise TypeError("Argument type not supported.")

adcc/gradients/__init__.py

@@ -159,7 +159,8 @@ def nuclear_gradient(excitation_or_mp):
         if hf.has_core_occupied_space:
             delta_IJ = hf.density.cc
 
-            g2_hf.oooo = -0.25 * einsum("ik,jl->ijkl", delta_ij, delta_ij)
+            g2_hf.oooo = 0.25 * (- einsum("li,jk->ijkl", delta_ij, delta_ij)
+                                 + einsum("ki,jl->ijkl", delta_ij, delta_ij))
             g2_hf.cccc = -0.5 * einsum("IK,JL->IJKL", delta_IJ, delta_IJ)
             g2_hf.ococ = -1.0 * einsum("ik,JL->iJkL", delta_ij, delta_IJ)
 
@@ -168,9 +169,7 @@ def nuclear_gradient(excitation_or_mp):
                 + einsum("ik,JL->iJkL", delta_ij, g1o.cc + 2.0 * delta_IJ)
                 + einsum("ik,JL->iJkL", g1o.oo, delta_IJ)
             )
-            g2_oresp.oooo = 0.25 * (
-                einsum("ik,jl->ijkl", delta_ij, g1o.oo + delta_ij)
-            )
+            g2_oresp.oooo = einsum("ij,kl->kilj", delta_ij, g1o.oo)
             g2_oresp.ovov = einsum("ij,ab->iajb", delta_ij, g1o.vv)
             g2_oresp.cvcv = einsum("IJ,ab->IaJb", delta_IJ, g1o.vv)
             g2_oresp.ocov = 2 * einsum("ik,Ja->iJka", delta_ij, g1o.cv)

adcc/gradients/amplitude_response.py

@@ -642,11 +642,8 @@ def ampl_relaxed_dms_cvs_adc2x(exci):
     g2a.ccvv = - 1.0 * t2ccvv
     g2a.ocvv = - 1.0 * t2ocvv
     g2a.ococ = 1.0 * einsum("iJab,kLab->iJkL", u.pphh, u.pphh)
-    g2a.vvvv = 1.0 * einsum("iJcd,iJab->abcd", u.pphh, u.pphh)
+    g2a.vvvv = 2.0 * einsum("iJcd,iJab->abcd", u.pphh, u.pphh)
 
-    # TODO: remove
-    # g2a.ococ *= 0.0
-    # g2a.vvvv *= 0.0
 
     g1a.co = (
         - 1.0 * einsum('JbKc,ibKc->Ji', g2a.cvcv, hf.ovcv)

adcc/gradients/orbital_response.py

@@ -33,7 +33,6 @@ def orbital_response_rhs(hf, g1a, g2a):
     response given amplitude-relaxed density matrices (method-specific)
     """
     # TODO: only add non-zero blocks to equations!
-
     if hf.has_core_occupied_space:
         ret_ov = -1.0 * (
             + 2.0 * einsum("JiKa,JK->ia", hf.cocv, g1a.cc)
@@ -78,7 +77,7 @@ def orbital_response_rhs(hf, g1a, g2a):
             + 2.0 * einsum('kJIb,kJab->Ia', g2a.occv, hf.ocvv)
             - 1.0 * einsum('abcd,Ibcd->Ia', g2a.vvvv, hf.cvvv)  # cvs-adc2x
             - 2.0 * einsum('jakb,jIkb->Ia', g2a.ovov, hf.ocov)  # cvs-adc2x
-            - 2.0 * einsum('jIlK,lKja->Ia', g2a.ococ, hf.ocov)  # cvs-adc2x
+            + 2.0 * einsum('jIlK,lKja->Ia', g2a.ococ, hf.ocov)  # cvs-adc2x
         )
         ret = AmplitudeVector(cv=ret_cv, ov=ret_ov)
     else:
@@ -112,7 +111,6 @@ def orbital_response_rhs(hf, g1a, g2a):
 
 def energy_weighted_density_matrix(hf, g1o, g2a):
     if hf.has_core_occupied_space:
-        # CVS-ADC0, CVS-ADC1, CVS-ADC2
         w = OneParticleOperator(hf)
         w.cc = - 0.5 * (
             + einsum("JKab,IKab->IJ", g2a.ccvv, hf.ccvv)
@@ -295,13 +293,6 @@ def __matmul__(self, lam):
             ret = AmplitudeVector(cv=ret_cv, ov=ret_ov)
         else:
             ret = AmplitudeVector(ov=ret_ov)
-        # TODO: generalize once other solvent methods are available
-        if self.hf.environment == "pe":
-            # PE contribution to the orbital Hessian
-            ops = self.hf.operators
-            dm = OneParticleOperator(self.hf, is_symmetric=True)
-            dm.ov = lam.ov
-            ret += ops.pe_induction_elec(dm).ov
         return evaluate(ret)
 
 

adcc/testdata/dump_fdiff_gradient.py

@@ -36,9 +36,6 @@
 
 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.]
 
 
 @dataclass
@@ -137,7 +134,7 @@ def main():
         "cvs-adc0",
         "cvs-adc1",
         "cvs-adc2",
-        # "cvs-adc2x",  # TODO: broken
+        "cvs-adc2x",
     ]
     molnames = [
         "h2o",