diff --git a/adcc/Excitation.py b/adcc/Excitation.py
index 20e315c0..89d5c626 100644
--- a/adcc/Excitation.py
+++ b/adcc/Excitation.py
@@ -66,6 +66,8 @@ def __init__(self, parent_state, index, method):
         self.__parent_state = parent_state
         self.index = index
         self.method = method
+        self.ground_state = parent_state.ground_state
+        self.reference_state = parent_state.reference_state
         for key in self.parent_state.excitation_property_keys:
             fget = getattr(type(self.parent_state), key).fget
             # Extract the kwargs passed to mark_excitation_property
diff --git a/adcc/ExcitedStates.py b/adcc/ExcitedStates.py
index 56c32aed..2a886417 100644
--- a/adcc/ExcitedStates.py
+++ b/adcc/ExcitedStates.py
@@ -217,6 +217,14 @@ def __add__(self, other):
         ret += other
         return ret
 
+    @property
+    @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
     @mark_excitation_property(transform_to_ao=True)
     @timed_member_call(timer="_property_timer")
diff --git a/adcc/ReferenceState.py b/adcc/ReferenceState.py
index fcb8d0c6..49e249cd 100644
--- a/adcc/ReferenceState.py
+++ b/adcc/ReferenceState.py
@@ -156,7 +156,9 @@ def __init__(self, hfdata, core_orbitals=None, frozen_core=None,
         )
 
         self.environment = None  # no environment attached by default
-        for name in ["excitation_energy_corrections", "environment"]:
+        for name in ["excitation_energy_corrections",
+                     "environment",
+                     "gradient_provider"]:
             if hasattr(hfdata, name):
                 setattr(self, name, getattr(hfdata, name))
 
diff --git a/adcc/__init__.py b/adcc/__init__.py
index a87d0a0b..6c4e750a 100644
--- a/adcc/__init__.py
+++ b/adcc/__init__.py
@@ -43,6 +43,7 @@
 from .AmplitudeVector import AmplitudeVector
 from .OneParticleOperator import OneParticleOperator
 from .opt_einsum_integration import register_with_opt_einsum
+from .gradients import nuclear_gradient
 
 # This has to be the last set of import
 from .guess import (guess_symmetries, guess_zero, guesses_any, guesses_singlet,
@@ -63,6 +64,7 @@
            "guess_symmetries", "guesses_spin_flip", "guess_zero", "LazyMp",
            "adc0", "cis", "adc1", "adc2", "adc2x", "adc3",
            "cvs_adc0", "cvs_adc1", "cvs_adc2", "cvs_adc2x", "cvs_adc3",
+           "nuclear_gradient",
            "banner"]
 
 __version__ = "0.15.17"
diff --git a/adcc/adc_pp/matrix.py b/adcc/adc_pp/matrix.py
index 80c489ca..f4129d97 100644
--- a/adcc/adc_pp/matrix.py
+++ b/adcc/adc_pp/matrix.py
@@ -21,7 +21,7 @@
 ##
 ## ---------------------------------------------------------------------
 from math import sqrt
-from collections import namedtuple
+from dataclasses import dataclass
 
 from adcc import block as b
 from adcc.functions import direct_sum, einsum, zeros_like
@@ -42,12 +42,16 @@
 #
 # Dispatch routine
 #
-"""
-`apply` is a function mapping an AmplitudeVector to the contribution of this
-block to the result of applying the ADC matrix. `diagonal` is an `AmplitudeVector`
-containing the expression to the diagonal of the ADC matrix from this block.
-"""
-AdcBlock = namedtuple("AdcBlock", ["apply", "diagonal"])
+@dataclass
+class AdcBlock:
+    """
+    `apply` is a function mapping an AmplitudeVector to the contribution of this
+    block to the result of applying the ADC matrix. `diagonal` is an
+    `AmplitudeVector` containing the expression to the diagonal of the ADC matrix
+    from this block.
+    """
+    apply: callable
+    diagonal: AmplitudeVector
 
 
 def block(ground_state, spaces, order, variant=None, intermediates=None):
diff --git a/adcc/backends/psi4.py b/adcc/backends/psi4.py
index 6d38a3a7..49e7b843 100644
--- a/adcc/backends/psi4.py
+++ b/adcc/backends/psi4.py
@@ -24,6 +24,7 @@
 
 from libadcc import HartreeFockProvider
 from adcc.misc import cached_property
+from adcc.gradients import GradientComponents
 
 import psi4
 
@@ -32,6 +33,72 @@
 from ..ExcitedStates import EnergyCorrection
 
 
+class Psi4GradientProvider:
+    def __init__(self, wfn):
+        self.wfn = wfn
+        self.mol = self.wfn.molecule()
+        self.backend = "psi4"
+        self.mints = psi4.core.MintsHelper(self.wfn)
+
+    def correlated_gradient(self, g1_ao, w_ao, g2_ao_1, g2_ao_2):
+        """
+        g1_ao: relaxed one-particle density matrix
+        w_ao:  energy-weighted density matrix
+        g2_ao_1, g2_ao_2: relaxed two-particle density matrices
+        """
+        natoms = self.mol.natom()
+        Gradient = {}
+        Gradient["N"] = np.zeros((natoms, 3))
+        Gradient["S"] = np.zeros((natoms, 3))
+        Gradient["T"] = np.zeros((natoms, 3))
+        Gradient["V"] = np.zeros((natoms, 3))
+        Gradient["OEI"] = np.zeros((natoms, 3))
+        Gradient["TEI"] = np.zeros((natoms, 3))
+        Gradient["Total"] = np.zeros((natoms, 3))
+
+        # 1st Derivative of Nuclear Repulsion
+        Gradient["N"] = psi4.core.Matrix.to_array(
+            self.mol.nuclear_repulsion_energy_deriv1([0, 0, 0])
+        )
+        # Build Integral Derivatives
+        cart = ['_X', '_Y', '_Z']
+        oei_dict = {"S": "OVERLAP", "T": "KINETIC", "V": "POTENTIAL"}
+
+        deriv1_mat = {}
+        deriv1_np = {}
+        # 1st Derivative of OEIs
+        for atom in range(natoms):
+            for key in oei_dict:
+                string = key + str(atom)
+                deriv1_mat[string] = self.mints.ao_oei_deriv1(oei_dict[key], atom)
+                for p in range(3):
+                    map_key = string + cart[p]
+                    deriv1_np[map_key] = np.asarray(deriv1_mat[string][p])
+                    if key == "S":
+                        Gradient["S"][atom, p] = np.sum(w_ao * deriv1_np[map_key])
+                    else:
+                        Gradient[key][atom, p] = np.sum(g1_ao * deriv1_np[map_key])
+        # Build Total OEI Gradient
+        Gradient["OEI"] = Gradient["T"] + Gradient["V"] + Gradient["S"]
+
+        # Build TE contributions
+        for atom in range(natoms):
+            deriv2 = self.mints.ao_tei_deriv1(atom)
+            for p in range(3):
+                deriv2_np = np.asarray(deriv2[p])
+                Gradient["TEI"][atom, p] += np.einsum(
+                    'pqrs,prqs->', g2_ao_1, deriv2_np, optimize=True
+                )
+                Gradient["TEI"][atom, p] -= np.einsum(
+                    'pqrs,psqr->', g2_ao_2, deriv2_np, optimize=True
+                )
+        ret = GradientComponents(
+            natoms, Gradient["N"], Gradient["S"],
+            Gradient["T"] + Gradient["V"], Gradient["TEI"]
+        )
+        return ret
+
+
 class Psi4OperatorIntegralProvider:
     def __init__(self, wfn):
         self.wfn = wfn
@@ -130,6 +197,7 @@ def __init__(self, wfn):
                                           wfn.nalpha(), wfn.nbeta(),
                                           self.restricted)
         self.operator_integral_provider = Psi4OperatorIntegralProvider(self.wfn)
+        self.gradient_provider = Psi4GradientProvider(self.wfn)
 
         self.environment = None
         self.environment_implementation = None
diff --git a/adcc/backends/pyscf.py b/adcc/backends/pyscf.py
index f2a0a72f..bcfd0424 100644
--- a/adcc/backends/pyscf.py
+++ b/adcc/backends/pyscf.py
@@ -24,15 +24,76 @@
 
 from libadcc import HartreeFockProvider
 from adcc.misc import cached_property
+from adcc.gradients import GradientComponents
 
 from .EriBuilder import EriBuilder
 from ..exceptions import InvalidReference
 from ..ExcitedStates import EnergyCorrection
 
-from pyscf import ao2mo, gto, scf
+from pyscf import ao2mo, gto, scf, grad
 from pyscf.solvent import ddcosmo
 
 
+class PyScfGradientProvider:
+    def __init__(self, scfres):
+        self.scfres = scfres
+        self.mol = self.scfres.mol
+        self.backend = "pyscf"
+
+    def correlated_gradient(self, g1_ao, w_ao, g2_ao_1, g2_ao_2):
+        natoms = self.mol.natm
+        Gradient = {}
+        Gradient["N"] = np.zeros((natoms, 3))
+        Gradient["S"] = np.zeros((natoms, 3))
+        Gradient["T+V"] = np.zeros((natoms, 3))
+        Gradient["OEI"] = np.zeros((natoms, 3))
+        Gradient["TEI"] = np.zeros((natoms, 3))
+        Gradient["Total"] = np.zeros((natoms, 3))
+
+        # TODO: does RHF/UHF matter here?
+        gradient = grad.RHF(self.scfres)
+        hcore_deriv = gradient.hcore_generator()
+        Sx = -1.0 * self.mol.intor('int1e_ipovlp', aosym='s1')
+        ERIx = -1.0 * self.mol.intor('int2e_ip1', aosym='s1')
+
+        ao_slices = self.mol.aoslice_by_atom()
+        for ia in range(natoms):
+            # TODO: only contract/compute with specific slices
+            # of density matrices (especially TPDM)
+            # this requires a lot of work however...
+            k0, k1 = ao_slices[ia, 2:]
+
+            # derivative of the overlap matrix
+            Sx_a = np.zeros_like(Sx)
+            Sx_a[:, k0:k1] = Sx[:, k0:k1]
+            Sx_a += Sx_a.transpose(0, 2, 1)
+            Gradient["S"][ia] += np.einsum("xpq,pq->x", Sx_a, w_ao)
+
+            # derivative of the core Hamiltonian
+            Hx_a = hcore_deriv(ia)
+            Gradient["T+V"][ia] += np.einsum("xpq,pq->x", Hx_a, g1_ao)
+
+            # derivatives of the ERIs
+            ERIx_a = np.zeros_like(ERIx)
+            ERIx_a[:, k0:k1] = ERIx[:, k0:k1]
+            ERIx_a += (
+                + ERIx_a.transpose(0, 2, 1, 4, 3)
+                + ERIx_a.transpose(0, 3, 4, 1, 2)
+                + ERIx_a.transpose(0, 4, 3, 2, 1)
+            )
+            Gradient["TEI"][ia] += np.einsum(
+                "pqrs,xprqs->x", g2_ao_1, ERIx_a, optimize=True
+            )
+            Gradient["TEI"][ia] -= np.einsum(
+                "pqrs,xpsqr->x", g2_ao_2, ERIx_a, optimize=True
+            )
+        ret = GradientComponents(
+            natoms, gradient.grad_nuc(), Gradient["S"],
+            Gradient["T+V"], Gradient["TEI"]
+        )
+        return ret
+
+
 class PyScfOperatorIntegralProvider:
     def __init__(self, scfres):
         self.scfres = scfres
@@ -196,7 +257,7 @@ def coefficients(self):
 
     def compute_mo_eri(self, blocks, spins):
         coeffs = tuple(self.coefficients[blocks[i] + spins[i]] for i in range(4))
-        # TODO Pyscf usse HDF5 internal to do the AO2MO here we read it all
+        # TODO Pyscf uses HDF5 internal to do the AO2MO here we read it all
         #      into memory. This wastes memory and could be avoided if temporary
         #      files were used instead. These could be deleted on the call
         #      to `flush_cache` automatically.
@@ -223,6 +284,7 @@ def __init__(self, scfres):
         self.operator_integral_provider = PyScfOperatorIntegralProvider(
             self.scfres
         )
+        self.gradient_provider = PyScfGradientProvider(self.scfres)
         if not self.restricted:
             assert self.scfres.mo_coeff[0].shape[1] == \
                 self.scfres.mo_coeff[1].shape[1]
diff --git a/adcc/gradients/TwoParticleDensityMatrix.py b/adcc/gradients/TwoParticleDensityMatrix.py
new file mode 100644
index 00000000..441d4337
--- /dev/null
+++ b/adcc/gradients/TwoParticleDensityMatrix.py
@@ -0,0 +1,290 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2021 by the adcc authors
+##
+## This file is part of adcc.
+##
+## adcc is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published
+## by the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## adcc is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with adcc. If not, see <http://www.gnu.org/licenses/>.
+##
+## ---------------------------------------------------------------------
+import numpy as np
+
+import libadcc
+import adcc.block as b
+from adcc.MoSpaces import split_spaces
+from adcc.Tensor import Tensor
+from adcc.functions import einsum
+
+
+class TwoParticleDensityMatrix:
+    """
+    Two-particle density matrix (TPDM) used for gradient evaluations
+    """
+    def __init__(self, spaces):
+        if hasattr(spaces, "mospaces"):
+            self.mospaces = spaces.mospaces
+        else:
+            self.mospaces = spaces
+        # Set reference_state if possible
+        if isinstance(spaces, libadcc.ReferenceState):
+            self.reference_state = spaces
+        elif hasattr(spaces, "reference_state"):
+            assert isinstance(spaces.reference_state, libadcc.ReferenceState)
+            self.reference_state = spaces.reference_state
+
+        occs = sorted(self.mospaces.subspaces_occupied, reverse=True)
+        virts = sorted(self.mospaces.subspaces_virtual, reverse=True)
+        self.orbital_subspaces = occs + virts
+        # check that orbital subspaces are correct
+        assert sum(self.mospaces.n_orbs(ss) for ss in self.orbital_subspaces) \
+            == self.mospaces.n_orbs("f")
+        # set the canonical blocks explicitly
+        self.blocks = [
+            b.oooo, b.ooov, b.oovv,
+            b.ovov, b.ovvv, b.vvvv,
+        ]
+        if self.mospaces.has_core_occupied_space:
+            self.blocks += [
+                b.cccc, b.ococ, b.cvcv,
+                b.ocov, b.cccv, b.cocv, b.ocoo,
+                b.ccco, b.occv, b.ccvv, b.ocvv,
+            ]
+        # make sure we didn't add any block twice!
+        assert len(list(set(self.blocks))) == len(self.blocks)
+        self._tensors = {}
+
+    @property
+    def shape(self):
+        """
+        Returns the shape tuple of the TwoParticleDensityMatrix
+        """
+        size = self.mospaces.n_orbs("f")
+        return 4 * (size,)
+
+    @property
+    def size(self):
+        """
+        Returns the number of elements of the TwoParticleDensityMatrix
+        """
+        return np.prod(self.shape)
+
+    def __setitem__(self, block, tensor):
+        """
+        Assigns a tensor to the specified block
+        """
+        if block not in self.blocks:
+            raise KeyError(f"Invalid block {block} assigned. "
+                           f"Available blocks are: {self.blocks}.")
+        s1, s2, s3, s4 = split_spaces(block)
+        expected_shape = (self.mospaces.n_orbs(s1),
+                          self.mospaces.n_orbs(s2),
+                          self.mospaces.n_orbs(s3),
+                          self.mospaces.n_orbs(s4))
+        if expected_shape != tensor.shape:
+            raise ValueError("Invalid shape of incoming tensor. "
+                             f"Expected shape {expected_shape}, but "
+                             f"got shape {tensor.shape} instead.")
+        self._tensors[block] = tensor
+
+    def __getitem__(self, block):
+        if block not in self.blocks:
+            raise KeyError(f"Invalid block {block} requested. "
+                           f"Available blocks are: {self.blocks}.")
+        if block not in self._tensors:
+            sym = libadcc.make_symmetry_eri(self.mospaces, block)
+            self._tensors[block] = Tensor(sym)
+        return self._tensors[block]
+
+    def to_ndarray(self):
+        raise NotImplementedError("ndarray export not implemented for TPDM.")
+
+    def copy(self):
+        """
+        Return a deep copy of the TwoParticleDensityMatrix
+        """
+        ret = TwoParticleDensityMatrix(self.mospaces)
+        for bl in self.blocks_nonzero:
+            ret[bl] = self.block(bl).copy()
+        if hasattr(self, "reference_state"):
+            ret.reference_state = self.reference_state
+        return ret
+
+    @property
+    def blocks_nonzero(self):
+        """
+        Returns a list of the non-zero block labels
+        """
+        return [b for b in self.blocks if b in self._tensors]
+
+    def is_zero_block(self, block):
+        """
+        Checks if block is explicitly marked as zero block.
+        Returns False if the block does not exist.
+        """
+        if block not in self.blocks:
+            return False
+        return block not in self.blocks_nonzero
+
+    def block(self, block):
+        """
+        Returns tensor of the given block.
+        Does not create a block in case it is marked as a zero block.
+        Use __getitem__ for that purpose.
+        """
+        if block not in self.blocks_nonzero:
+            raise KeyError("The block function does not support "
+                           "access to zero-blocks. Available non-zero "
+                           f"blocks are: {self.blocks_nonzero}.")
+        return self._tensors[block]
+
+    def __getattr__(self, attr):
+        return self.__getitem__(b.__getattr__(attr))
+
+    def __setattr__(self, attr, value):
+        try:
+            self.__setitem__(b.__getattr__(attr), value)
+        except AttributeError:
+            super().__setattr__(attr, value)
+
+    def set_zero_block(self, block):
+        """
+        Set a given block as zero block
+        """
+        if block not in self.blocks:
+            raise KeyError(f"Invalid block {block} set as zero block. "
+                           f"Available blocks are: {self.blocks}.")
+        self._tensors.pop(block)
+
+    def __transform_to_ao(self, refstate_or_coefficients):
+        if not len(self.blocks_nonzero):
+            raise ValueError("At least one non-zero block is needed to "
+                             "transform the TwoParticleDensityMatrix.")
+        if isinstance(refstate_or_coefficients, libadcc.ReferenceState):
+            hf = refstate_or_coefficients
+            coeff_map = {}
+            for sp in self.orbital_subspaces:
+                coeff_map[sp + "_a"] = hf.orbital_coefficients_alpha(sp + "b")
+                coeff_map[sp + "_b"] = hf.orbital_coefficients_beta(sp + "b")
+        else:
+            coeff_map = refstate_or_coefficients
+
+        g2_ao_1 = 0
+        g2_ao_2 = 0
+        transf = "ip,jq,ijkl,kr,ls->pqrs"
+        cc = coeff_map
+        for block in self.blocks_nonzero:
+            s1, s2, s3, s4 = split_spaces(block)
+            ten = self[block]
+            aaaa = einsum(transf, cc[f"{s1}_a"], cc[f"{s2}_a"],
+                          ten, cc[f"{s3}_a"], cc[f"{s4}_a"])
+            bbbb = einsum(transf, cc[f"{s1}_b"], cc[f"{s2}_b"],
+                          ten, cc[f"{s3}_b"], cc[f"{s4}_b"])
+            g2_ao_1 += (
+                + aaaa
+                + bbbb
+                + einsum(transf, cc[f"{s1}_a"], cc[f"{s2}_b"],
+                         ten, cc[f"{s3}_a"], cc[f"{s4}_b"])  # abab
+                + einsum(transf, cc[f"{s1}_b"], cc[f"{s2}_a"],
+                         ten, cc[f"{s3}_b"], cc[f"{s4}_a"])  # baba
+            )
+            g2_ao_2 += (
+                + aaaa
+                + bbbb
+                + einsum(transf, cc[f"{s1}_a"], cc[f"{s2}_b"],
+                         ten, cc[f"{s3}_b"], cc[f"{s4}_a"])  # abba
+                + einsum(transf, cc[f"{s1}_b"], cc[f"{s2}_a"],
+                         ten, cc[f"{s3}_a"], cc[f"{s4}_b"])  # baab
+            )
+        return (g2_ao_1.evaluate(), g2_ao_2.evaluate())
+
+    def to_ao_basis(self, refstate_or_coefficients=None):
+        """
+        Transform the density to the AO basis for contraction
+        with two-electron integrals.
+        ALL coefficients are already accounted for in the density matrix.
+        Two blocks are returned, the first one needs to be contracted with
+        prqs, the second one with -psqr (in Chemists' notation).
+        """
+        if isinstance(refstate_or_coefficients, (dict, libadcc.ReferenceState)):
+            return self.__transform_to_ao(refstate_or_coefficients)
+        elif refstate_or_coefficients is None:
+            if not hasattr(self, "reference_state"):
+                raise ValueError("Argument reference_state is required if no "
+                                 "reference_state is stored in the "
+                                 "TwoParticleDensityMatrix")
+            return self.__transform_to_ao(self.reference_state)
+        else:
+            raise TypeError("Argument type not supported.")
+
+    def __iadd__(self, other):
+        if self.mospaces != other.mospaces:
+            raise ValueError("Cannot add TwoParticleDensityMatrices with "
+                             "differing mospaces.")
+
+        for bl in other.blocks_nonzero:
+            if self.is_zero_block(bl):
+                self[bl] = other.block(bl).copy()
+            else:
+                self[bl] = self.block(bl) + other.block(bl)
+
+        # Update ReferenceState pointer
+        if hasattr(self, "reference_state"):
+            if hasattr(other, "reference_state") \
+                    and self.reference_state != other.reference_state:
+                delattr(self, "reference_state")
+        return self
+
+    def __isub__(self, other):
+        if self.mospaces != other.mospaces:
+            raise ValueError("Cannot subtract TwoParticleDensityMatrix with "
+                             "differing mospaces.")
+
+        for bl in other.blocks_nonzero:
+            if self.is_zero_block(bl):
+                self[bl] = -1.0 * other.block(bl)  # The copy is implicit
+            else:
+                self[bl] = self.block(bl) - other.block(bl)
+
+        # Update ReferenceState pointer
+        if hasattr(self, "reference_state"):
+            if hasattr(other, "reference_state") \
+                    and self.reference_state != other.reference_state:
+                delattr(self, "reference_state")
+        return self
+
+    def __imul__(self, other):
+        if not isinstance(other, (float, int)):
+            return NotImplemented
+        for bl in self.blocks_nonzero:
+            self[bl] = self.block(bl) * other
+        return self
+
+    def __add__(self, other):
+        return self.copy().__iadd__(other)
+
+    def __sub__(self, other):
+        return self.copy().__isub__(other)
+
+    def __mul__(self, other):
+        return self.copy().__imul__(other)
+
+    def __rmul__(self, other):
+        return self.copy().__imul__(other)
+
+    def evaluate(self):
+        for bl in self.blocks_nonzero:
+            self.block(bl).evaluate()
+        return self
diff --git a/adcc/gradients/__init__.py b/adcc/gradients/__init__.py
new file mode 100644
index 00000000..47b47bc5
--- /dev/null
+++ b/adcc/gradients/__init__.py
@@ -0,0 +1,217 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2021 by the adcc authors
+##
+## This file is part of adcc.
+##
+## adcc is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published
+## by the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## adcc is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with adcc. If not, see <http://www.gnu.org/licenses/>.
+##
+## ---------------------------------------------------------------------
+from dataclasses import dataclass
+from typing import Dict, Union, Optional
+import numpy as np
+
+from adcc.LazyMp import LazyMp
+from adcc.Excitation import Excitation
+from adcc.timings import Timer
+from adcc.functions import einsum, evaluate
+
+from adcc.OneParticleOperator import OneParticleOperator, product_trace
+from .TwoParticleDensityMatrix import TwoParticleDensityMatrix
+from .orbital_response import (
+    orbital_response, orbital_response_rhs, energy_weighted_density_matrix
+)
+from .amplitude_response import amplitude_relaxed_densities
+
+
+@dataclass(frozen=True)
+class GradientComponents:
+    natoms: int
+    nuc: np.ndarray
+    overlap: np.ndarray
+    hcore: np.ndarray
+    two_electron: np.ndarray
+    custom: Optional[Dict[str, np.ndarray]] = None
+
+    @property
+    def total(self):
+        """Returns the total gradient"""
+        ret = sum([self.nuc, self.overlap, self.hcore, self.two_electron])
+        if self.custom is None:
+            return ret
+        for c in self.custom:
+            ret += self.custom[c]
+
+    @property
+    def one_electron(self):
+        """Returns the one-electron gradient"""
+        return sum([self.nuc, self.overlap, self.hcore])
+
+
+@dataclass(frozen=True)
+class GradientResult:
+    excitation_or_mp: Union[LazyMp, Excitation]
+    components: GradientComponents
+    g1: OneParticleOperator
+    g2: TwoParticleDensityMatrix
+    timer: Timer
+    g1a: Optional[OneParticleOperator] = None
+    g2a: Optional[TwoParticleDensityMatrix] = None
+
+    @property
+    def reference_state(self):
+        return self.excitation_or_mp.reference_state
+
+    @property
+    def _energy(self):
+        """Compute energy based on density matrices
+        for testing purposes"""
+        if self.g1a is None:
+            raise ValueError("No unrelaxed one-particle "
+                             "density available.")
+        if self.g2a is None:
+            raise ValueError("No unrelaxed two-particle "
+                             "density available.")
+        ret = 0.0
+        hf = self.reference_state
+        for b in self.g1a.blocks_nonzero:
+            ret += self.g1a[b].dot(hf.fock(b))
+        for b in self.g2a.blocks_nonzero:
+            ret += self.g2a[b].dot(hf.eri(b))
+        return ret
+
+    @property
+    def dipole_moment_relaxed(self):
+        """Returns the orbital-relaxed electric dipole moment"""
+        return self.__dipole_moment_electric(self.g1)
+
+    @property
+    def dipole_moment_unrelaxed(self):
+        """Returns the unrelaxed electric dipole moment"""
+        if self.g1a is None:
+            raise ValueError("No unrelaxed one-particle "
+                             "density available.")
+        hf = self.reference_state
+        return self.__dipole_moment_electric(self.g1a + hf.density)
+
+    @property
+    def total(self):
+        """Returns the total gradient"""
+        return self.components.total
+
+    def __dipole_moment_electric(self, dm):
+        dips = self.reference_state.operators.electric_dipole
+        elec_dip = -1.0 * np.array(
+            [product_trace(dm, dip) for dip in dips]
+        )
+        return elec_dip + self.reference_state.nuclear_dipole
+
+
+def nuclear_gradient(excitation_or_mp):
+    if isinstance(excitation_or_mp, LazyMp):
+        mp = excitation_or_mp
+    elif isinstance(excitation_or_mp, Excitation):
+        mp = excitation_or_mp.ground_state
+    else:
+        raise TypeError("Gradient can only be computed for "
+                        "Excitation or LazyMp object.")
+
+    timer = Timer()
+    hf = mp.reference_state
+    with timer.record("amplitude_response"):
+        g1a, g2a = amplitude_relaxed_densities(excitation_or_mp)
+
+    with timer.record("orbital_response"):
+        rhs = orbital_response_rhs(hf, g1a, g2a).evaluate()
+        lam = orbital_response(hf, rhs)
+
+    # orbital-relaxed OPDM (without reference state)
+    g1o = g1a.copy()
+    g1o.ov = 0.5 * lam.ov
+    if hf.has_core_occupied_space:
+        g1o.cv = 0.5 * lam.cv
+    # orbital-relaxed OPDM (including reference state)
+    g1 = g1o.copy()
+    g1 += hf.density
+
+    with timer.record("energy_weighted_density_matrix"):
+        w = energy_weighted_density_matrix(hf, g1o, g2a)
+
+    # build two-particle density matrices for contraction with TEI
+    # prefactors see eqs 17 and A4 in 10.1063/1.5085117
+    with timer.record("form_tpdm"):
+        g2_hf = TwoParticleDensityMatrix(hf)
+        g2_oresp = TwoParticleDensityMatrix(hf)
+        delta_ij = hf.density.oo
+        if hf.has_core_occupied_space:
+            delta_IJ = hf.density.cc
+
+            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)
+
+            g2_oresp.cccc = einsum("IK,JL->IJKL", delta_IJ, g1o.cc + delta_IJ)
+            g2_oresp.ococ = (
+                + einsum("ik,JL->iJkL", delta_ij, g1o.cc + 2.0 * delta_IJ)
+                + einsum("ik,JL->iJkL", 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)
+            g2_oresp.cccv = 2 * einsum("IK,Ja->IJKa", delta_IJ, g1o.cv)
+            g2_oresp.ooov = 2 * einsum("ik,ja->ijka", delta_ij, g1o.ov)
+            g2_oresp.cocv = 2 * einsum("IK,ja->IjKa", delta_IJ, g1o.ov)
+            g2_oresp.ocoo = 2 * einsum("ik,Jl->iJkl", delta_ij, g1o.co)
+            g2_oresp.ccco = 2 * einsum("IK,Jl->IJKl", delta_IJ, g1o.co)
+
+            # scale for contraction with integrals
+            g2a.oovv *= 0.5
+            g2a.ccvv *= 0.5
+            g2a.occv *= 2.0
+            g2a.vvvv *= 0.25
+
+            g2_total = evaluate(g2_hf + g2a + g2_oresp)
+        else:
+            g2_hf.oooo = 0.25 * (- einsum("li,jk->ijkl", delta_ij, delta_ij)
+                                 + einsum("ki,jl->ijkl", delta_ij, 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.ooov = (- einsum("kj,ia->ijka", delta_ij, g1o.ov)
+                             + einsum("ki,ja->ijka", delta_ij, g1o.ov))
+
+            # scale for contraction with integrals
+            g2a.oovv *= 0.5
+            g2a.oooo *= 0.25
+            g2a.vvvv *= 0.25
+            g2_total = evaluate(g2_hf + g2a + g2_oresp)
+
+    with timer.record("transform_ao"):
+        g2_ao_1, g2_ao_2 = g2_total.to_ao_basis()
+        g2_ao_1, g2_ao_2 = g2_ao_1.to_ndarray(), g2_ao_2.to_ndarray()
+        g1_ao = sum(g1.to_ao_basis(hf)).to_ndarray()
+        w_ao = sum(w.to_ao_basis(hf)).to_ndarray()
+
+    with timer.record("contract_integral_derivatives"):
+        grad = hf.gradient_provider.correlated_gradient(
+            g1_ao, w_ao, g2_ao_1, g2_ao_2
+        )
+
+    ret = GradientResult(excitation_or_mp, grad, g1, g2_total,
+                         timer, g1a=g1a, g2a=g2a)
+    return ret
diff --git a/adcc/gradients/amplitude_response.py b/adcc/gradients/amplitude_response.py
new file mode 100644
index 00000000..7f9f37bc
--- /dev/null
+++ b/adcc/gradients/amplitude_response.py
@@ -0,0 +1,706 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2021 by the adcc authors
+##
+## This file is part of adcc.
+##
+## adcc is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published
+## by the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## adcc is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with adcc. If not, see <http://www.gnu.org/licenses/>.
+##
+## ---------------------------------------------------------------------
+from math import sqrt
+
+from adcc.functions import einsum, direct_sum, evaluate
+
+from .TwoParticleDensityMatrix import TwoParticleDensityMatrix
+from adcc.OneParticleOperator import OneParticleOperator
+from adcc.LazyMp import LazyMp
+from adcc.Excitation import Excitation
+import adcc.block as b
+
+
+def t2bar_oovv_adc2(exci, g1a_adc0):
+    mp = exci.ground_state
+    hf = mp.reference_state
+    u = exci.excitation_vector
+    df_ia = mp.df(b.ov)
+    t2bar = 0.5 * (
+        hf.oovv
+        - 2.0 * einsum(
+            "ijcb,ac->ijab", hf.oovv, g1a_adc0.vv
+        ).antisymmetrise(2, 3)
+        + 2.0 * einsum(
+            "kjab,ik->ijab", hf.oovv, g1a_adc0.oo
+        ).antisymmetrise(0, 1)
+        + 4.0 * einsum(
+            "ia,jkbc,kc->ijab", u.ph, hf.oovv, u.ph
+        ).antisymmetrise(2, 3).antisymmetrise(0, 1)
+    ) / (
+        2.0 * direct_sum("ia+jb->ijab", df_ia, df_ia).symmetrise(0, 1)
+    )
+    return t2bar
+
+
+def tbarD_oovv_adc3(exci, g1a_adc0):
+    mp = exci.ground_state
+    hf = mp.reference_state
+    u = exci.excitation_vector
+    df_ia = mp.df(b.ov)
+    # Table XI (10.1063/1.5085117)
+    tbarD = 0.5 * (
+        hf.oovv
+        - 2.0 * einsum(
+            "ijcb,ac->ijab", hf.oovv, g1a_adc0.vv
+        ).antisymmetrise(2, 3)
+        + 2.0 * einsum(
+            "kjab,ik->ijab", hf.oovv, g1a_adc0.oo
+        ).antisymmetrise(0, 1)
+        + 4.0 * einsum(
+            "ia,jkbc,kc->ijab", u.ph, hf.oovv, u.ph
+        ).antisymmetrise(2, 3).antisymmetrise(0, 1)
+    ) / (
+        2.0 * direct_sum("ia+jb->ijab", df_ia, df_ia).symmetrise(0, 1)
+    )
+    return tbarD
+
+
+def tbar_TD_oovv_adc3(exci, g1a_adc0):
+    mp = exci.ground_state
+    hf = mp.reference_state
+    # Table XI (10.1063/1.5085117)
+    tbarD = tbarD_oovv_adc3(exci, g1a_adc0)
+    tbarD.evaluate()
+    ret = 2.0 * 4.0 * (
+        + 2.0 * einsum("ikac,jckb->ijab", tbarD, hf.ovov)
+    ).antisymmetrise(0, 1).antisymmetrise(2, 3)
+    ret += 2.0 * (
+        - 1.0 * einsum("ijcd,cdab->ijab", tbarD, hf.vvvv)
+        - 1.0 * einsum("klab,ijkl->ijab", tbarD, hf.oooo)
+    )
+    return ret
+
+
+def rho_bar_adc3(exci, g1a_adc0):
+    mp = exci.ground_state
+    hf = mp.reference_state
+    u = exci.excitation_vector
+    df_ia = mp.df(b.ov)
+    rho_bar = 2.0 * (
+        + 1.0 * einsum("ijka,jk->ia", hf.ooov, g1a_adc0.oo)
+        + 1.0 * einsum("ijkb,jb,ka->ia", hf.ooov, u.ph, u.ph)
+        - 1.0 * einsum("icab,bc->ia", hf.ovvv, g1a_adc0.vv)
+        + 1.0 * einsum("jcab,jb,ic->ia", hf.ovvv, u.ph, u.ph)
+    ) / df_ia
+    return rho_bar
+
+
+def tbar_rho_oovv_adc3(exci, g1a_adc0):
+    mp = exci.ground_state
+    hf = mp.reference_state
+    # Table XI (10.1063/1.5085117)
+    rho_bar = rho_bar_adc3(exci, g1a_adc0)
+    ret = (
+        - 1.0 * 2.0 * einsum("jcab,ic->ijab", hf.ovvv, rho_bar).antisymmetrise(0, 1)
+        - 1.0 * 2.0 * einsum("ijkb,ka->ijab", hf.ooov, rho_bar).antisymmetrise(2, 3)
+    )
+    return ret
+
+
+def t2bar_oovv_adc3(exci, g1a_adc0, g2a_adc1):
+    mp = exci.ground_state
+    hf = mp.reference_state
+    u = exci.excitation_vector
+    tbar_TD = tbar_TD_oovv_adc3(exci, g1a_adc0)
+    tbar_TD.evaluate()
+    tbar_rho = tbar_rho_oovv_adc3(exci, g1a_adc0)
+    tbar_rho.evaluate()
+    t2 = mp.t2(b.oovv).evaluate()
+    df_ia = mp.df(b.ov)
+    rx = einsum("jb,ijab->ia", u.ph, t2).evaluate()
+
+    ttilde1 = (
+        - 1.0 * einsum("klab,ijkm,lm->ijab", t2, hf.oooo, g1a_adc0.oo)
+        - 1.0 * 2.0 * einsum(
+            "ka,ijkl,lb->ijab", u.ph, hf.oooo, rx
+        ).antisymmetrise(2, 3)
+        + 1.0 * 2.0 * (
+            + 0.5 * einsum("ik,jklm,lmab->ijab", g1a_adc0.oo, hf.oooo, t2)
+            - 2.0 * einsum("jkab,lm,ilkm->ijab", t2, g1a_adc0.oo, hf.oooo)
+        ).antisymmetrise(0, 1)
+        - 1.0 * 4.0 * (
+            + 0.5 * einsum("ia,kc,jklm,lmbc->ijab", u.ph, u.ph, hf.oooo, t2)
+            - 2.0 * einsum("ka,likm,jmbc,lc->ijab", u.ph, hf.oooo, t2, u.ph)
+        ).antisymmetrise(0, 1).antisymmetrise(2, 3)
+    )
+    ttilde1.evaluate()
+
+    ttilde2 = 4.0 * (
+        - 1.0 * einsum("ijkc,lc,klab->ijab", hf.ooov, u.ph, u.pphh)
+        + 1.0 * 2.0 * einsum(
+            "ikab,jlkc,lc->ijab", u.pphh, hf.ooov, u.ph
+        ).antisymmetrise(0, 1)
+        - 1.0 * 4.0 * einsum(
+            "jklb,kc,ilac->ijab", hf.ooov, u.ph, u.pphh
+        ).antisymmetrise(0, 1).antisymmetrise(2, 3)
+    )
+    ttilde2.evaluate()
+
+    ttilde3 = (
+        + 1.0 * 2.0 * einsum(
+            "ijbc,ac->ijab", hf.oovv, g1a_adc0.vv
+        ).antisymmetrise(2, 3)
+        - 1.0 * 2.0 * einsum(
+            "jkab,ik->ijab", hf.oovv, g1a_adc0.oo
+        ).antisymmetrise(0, 1)
+        + 1.0 * 4.0 * einsum(
+            "ia,jkbc,kc->ijab", u.ph, hf.oovv, u.ph
+        ).antisymmetrise(0, 1).antisymmetrise(2, 3)
+    )
+    ttilde3.evaluate()
+
+    # TODO: intermediate x_ka <ic||kd>?
+    # TODO: intermediate x_lc t_jlbd?
+    ttilde4 = (
+        + 1.0 * 2.0 * (
+            - 2.0 * einsum("jkab,ickd,cd->ijab", t2, hf.ovov, g1a_adc0.vv)
+            - 2.0 * einsum("ic,kd,jcld,klab->ijab", u.ph, u.ph, hf.ovov, t2)
+            + 1.0 * einsum("ic,jkab,ld,lckd->ijab", u.ph, t2, u.ph, hf.ovov)
+            + 1.0 * einsum("jkab,kc,ld,lcid->ijab", t2, u.ph, u.ph, hf.ovov)
+        ).antisymmetrise(0, 1)
+        + 1.0 * 2.0 * (
+            - 2.0 * einsum("ijcb,kalc,kl->ijab", t2, hf.ovov, g1a_adc0.oo)
+            - 2.0 * einsum("ka,lc,kbld,ijcd->ijab", u.ph, u.ph, hf.ovov, t2)
+            + 1.0 * einsum("ka,ijbc,ld,kdlc->ijab", u.ph, t2, u.ph, hf.ovov)
+            + 1.0 * einsum("ijbc,kc,ld,kdla->ijab", t2, u.ph, u.ph, hf.ovov)
+        ).antisymmetrise(2, 3)
+        + 1.0 * 4.0 * (
+            + 1.0 * einsum("ac,jdkc,ikbd->ijab", g1a_adc0.vv, hf.ovov, t2)
+            - 1.0 * einsum("jckb,ikad,cd->ijab", hf.ovov, t2, g1a_adc0.vv)
+            - 1.0 * einsum("ik,kclb,jlac->ijab", g1a_adc0.oo, hf.ovov, t2)
+            + 1.0 * einsum("jckb,ilac,lk->ijab", hf.ovov, t2, g1a_adc0.oo)
+            - 1.0 * einsum("ic,jckb,ka->ijab", u.ph, hf.ovov, rx)
+            - 1.0 * einsum("jb,kc,lcid,klad->ijab", u.ph, u.ph, hf.ovov, t2)
+            - 1.0 * einsum("ka,jckb,ic->ijab", u.ph, hf.ovov, rx)
+            - 1.0 * einsum("jb,kc,lakd,ilcd->ijab", u.ph, u.ph, hf.ovov, t2)
+            + 2.0 * einsum("ka,ickd,ld,jlbc->ijab", u.ph, hf.ovov, u.ph, t2)
+            + 2.0 * einsum("ic,kcla,kd,jlbd->ijab", u.ph, hf.ovov, u.ph, t2)
+        ).antisymmetrise(0, 1).antisymmetrise(2, 3)
+    )
+    ttilde4.evaluate()
+
+    ttilde5 = - 4.0 * (
+        + 1.0 * einsum("kcab,kd,ijcd->ijab", hf.ovvv, u.ph, u.pphh)
+        + 1.0 * 2.0 * einsum(
+            "ijbc,kcad,kd->ijab", u.pphh, hf.ovvv, u.ph
+        ).antisymmetrise(2, 3)
+        + 1.0 * 4.0 * einsum(
+            "jcbd,kd,ikac->ijab", hf.ovvv, u.ph, u.pphh
+        ).antisymmetrise(0, 1).antisymmetrise(2, 3)
+    )
+    ttilde5.evaluate()
+
+    ttilde6 = (
+        - 1.0 * einsum("ijcd,abde,ce->ijab", t2, hf.vvvv, g1a_adc0.vv)
+        - 1.0 * 2.0 * einsum(
+            "ic,abcd,jkde,ke->ijab", u.ph, hf.vvvv, t2, u.ph
+        ).antisymmetrise(0, 1)
+        - 1.0 * 2.0 * (
+            + 0.5 * einsum("ac,bcde,ijde->ijab", g1a_adc0.vv, hf.vvvv, t2)
+            - 2.0 * einsum("ijbc,de,adce->ijab", t2, g1a_adc0.vv, hf.vvvv)
+        ).antisymmetrise(2, 3)
+        - 1.0 * 4.0 * (
+            + 0.5 * einsum("ia,kc,bcde,jkde->ijab", u.ph, u.ph, hf.vvvv, t2)
+            - 2.0 * einsum("ic,adce,jkeb,kd->ijab", u.ph, hf.vvvv, t2, u.ph)
+        ).antisymmetrise(0, 1).antisymmetrise(2, 3)
+    )
+    ttilde6.evaluate()
+
+    ret = 0.5 * (
+        hf.oovv
+        + ttilde1 + ttilde2 + ttilde3 + ttilde4 + ttilde5 + ttilde6
+        + tbar_TD + tbar_rho
+    ) / (2.0 * direct_sum("ia+jb->ijab", df_ia, df_ia).symmetrise(0, 1))
+    return ret
+
+
+def t2bar_oovv_cvs_adc2(exci, g1a_adc0):
+    mp = exci.ground_state
+    hf = mp.reference_state
+    df_ia = mp.df(b.ov)
+    t2bar = 0.5 * (
+        - einsum("ijcb,ac->ijab", hf.oovv, g1a_adc0.vv).antisymmetrise(2, 3)
+    ) / direct_sum("ia+jb->ijab", df_ia, df_ia).symmetrise(0, 1)
+    return t2bar
+
+
+def ampl_relaxed_dms_mp2(mp):
+    hf = mp.reference_state
+    t2 = mp.t2(b.oovv)
+    g1a = OneParticleOperator(hf)
+    g2a = TwoParticleDensityMatrix(hf)
+    g1a.oo = -0.5 * einsum('ikab,jkab->ij', t2, t2)
+    g1a.vv = 0.5 * einsum('ijac,ijbc->ab', t2, t2)
+    g2a.oovv = -1.0 * mp.t2(b.oovv)
+    return g1a, g2a
+
+
+def ampl_relaxed_dms_adc0(exci):
+    hf = exci.reference_state
+    u = exci.excitation_vector
+    g1a = OneParticleOperator(hf)
+    g2a = TwoParticleDensityMatrix(hf)
+    # g2a is not required for the adc0 gradient,
+    # but expected by amplitude_relaxed_densities
+    g1a.oo = - 1.0 * einsum("ia,ja->ij", u.ph, u.ph)
+    g1a.vv = + 1.0 * einsum("ia,ib->ab", u.ph, u.ph)
+    return g1a, g2a
+
+
+def ampl_relaxed_dms_adc1(exci):
+    hf = exci.reference_state
+    u = exci.excitation_vector
+    g1a = OneParticleOperator(hf)
+    g2a = TwoParticleDensityMatrix(hf)
+    g1a.oo = - 1.0 * einsum("ia,ja->ij", u.ph, u.ph)
+    g1a.vv = + 1.0 * einsum("ia,ib->ab", u.ph, u.ph)
+    g2a.ovov = - 1.0 * einsum("ja,ib->iajb", u.ph, u.ph)
+    return g1a, g2a
+
+
+def ampl_relaxed_dms_adc2(exci):
+    u = exci.excitation_vector
+    mp = exci.ground_state
+    g1a_adc1, g2a_adc1 = ampl_relaxed_dms_adc1(exci)
+    t2 = mp.t2(b.oovv)
+    t2bar = t2bar_oovv_adc2(exci, g1a_adc1).evaluate()
+
+    g1a = g1a_adc1.copy()
+    g1a.oo += (
+        - 2.0 * einsum('jkab,ikab->ij', u.pphh, u.pphh)
+        - 2.0 * einsum('jkab,ikab->ij', t2bar, t2).symmetrise((0, 1))
+    )
+    g1a.vv += (
+        + 2.0 * einsum("ijac,ijbc->ab", u.pphh, u.pphh)
+        + 2.0 * einsum("ijac,ijbc->ab", t2bar, t2).symmetrise((0, 1))
+    )
+
+    g2a = g2a_adc1.copy()
+    ru_ov = einsum("ijab,jb->ia", t2, u.ph)
+    g2a.oovv = (
+        0.5 * (
+            - 1.0 * t2
+            + 2.0 * einsum("ijcb,ca->ijab", t2, g1a_adc1.vv).antisymmetrise(2, 3)
+            - 2.0 * einsum("kjab,ki->ijab", t2, g1a_adc1.oo).antisymmetrise(0, 1)
+            - 4.0 * einsum(
+                "ia,jb->ijab", u.ph, ru_ov
+            ).antisymmetrise((0, 1)).antisymmetrise((2, 3))
+        )
+        - 2.0 * t2bar
+    )
+    g2a.ooov = -2.0 * einsum("kb,ijab->ijka", u.ph, u.pphh)
+    g2a.ovvv = -2.0 * einsum("ja,ijbc->iabc", u.ph, u.pphh)
+    return g1a, g2a
+
+
+def ampl_relaxed_dms_adc2x(exci):
+    u = exci.excitation_vector
+    g1a, g2a = ampl_relaxed_dms_adc2(exci)
+
+    g2a.ovov += -4.0 * einsum("ikbc,jkac->iajb", u.pphh, u.pphh)
+    g2a.oooo = 2.0 * einsum('ijab,klab->ijkl', u.pphh, u.pphh)
+    g2a.vvvv = 2.0 * einsum('ijcd,ijab->abcd', u.pphh, u.pphh)
+
+    return g1a, g2a
+
+
+def ampl_relaxed_dms_adc3(exci):
+    u = exci.excitation_vector
+    mp = exci.ground_state
+    hf = mp.reference_state
+    g1a_adc1, g2a_adc1 = ampl_relaxed_dms_adc1(exci)
+    t2bar = t2bar_oovv_adc3(exci, g1a_adc1, g2a_adc1).evaluate()
+    tbarD = tbarD_oovv_adc3(exci, g1a_adc1).evaluate()
+    rho_bar = rho_bar_adc3(exci, g1a_adc1).evaluate()
+    t2 = mp.t2(b.oovv)
+    tD = mp.td2(b.oovv)
+    rho = mp.mp2_diffdm
+
+    # Table IX (10.1063/1.5085117)
+    g1a = g1a_adc1.copy()
+    g1a.oo += (
+        - 2.0 * einsum("jkab,ikab->ij", u.pphh, u.pphh)
+        - 1.0 * 2.0 * (
+            + 1.0 * einsum("ikab,jkab->ij", t2, t2bar)
+            + 1.0 * einsum("ikab,jkab->ij", tD, tbarD)
+            + 0.5 * einsum("ia,ja->ij", rho_bar, rho.ov)
+        ).symmetrise(0, 1)
+    )
+    g1a.vv += (
+        + 2.0 * einsum("ijac,ijbc->ab", u.pphh, u.pphh)
+        + 1.0 * 2.0 * (
+            + 1.0 * einsum("ijac,ijbc->ab", t2bar, t2)
+            + 1.0 * einsum("ijac,ijbc->ab", tbarD, tD)
+            + 0.5 * einsum("ia,ib->ab", rho_bar, rho.ov)
+        ).symmetrise(0, 1)
+    )
+
+    tsq_ovov = einsum("ikac,jkbc->iajb", t2, t2).evaluate()
+    tsq_vvvv = einsum("klab,klcd->abcd", t2, t2).evaluate()
+    tsq_oooo = einsum("ijcd,klcd->ijkl", t2, t2).evaluate()
+    rx = einsum("ijab,jb->ia", t2, u.ph)
+
+    g2a = TwoParticleDensityMatrix(hf)
+    g2a.oooo = (
+        + 2.0 * einsum("ijab,klab->ijkl", u.pphh, u.pphh)
+        + 0.5 * 2.0 * (
+            + 2.0 * einsum("ijab,klab->ijkl", tbarD, t2)
+            + 0.5 * 2.0 * einsum(
+                "jm,imab,klab->ijkl", g1a_adc1.oo, t2, t2
+            ).antisymmetrise(0, 1)
+            + 1.0 * 2.0 * einsum(
+                "kc,ijbc,ma,lmab->ijkl", u.ph, t2, u.ph, t2
+            ).antisymmetrise(2, 3)
+            + 1.0 * 4.0 * (
+                + 1.0 * einsum("ik,jl->ijkl", rho.oo, g1a_adc1.oo)
+                - 1.0 * einsum("lajb,ia,kb->ijkl", tsq_ovov, u.ph, u.ph)
+            ).antisymmetrise(0, 1).antisymmetrise(2, 3)
+        ).symmetrise([(0, 2), (1, 3)])
+    )
+    g2a.ooov = (
+        - 2.0 * einsum("kb,ijab->ijka", u.ph, u.pphh)
+        + 1.0 * einsum("la,ijbc,klbc->ijka", u.ph, t2, u.pphh)
+        + 1.0 * 2.0 * (
+            + 2.0 * einsum("ic,jlab,klbc->ijka", u.ph, t2, u.pphh)
+            - 1.0 * einsum("ja,ilbc,klbc->ijka", u.ph, t2, u.pphh)
+            + 1.0 * einsum("ja,ik->ijka", rho.ov, g1a_adc1.oo)
+            + 1.0 * einsum("ia,jb,kb->ijka", u.ph, rho.ov, u.ph)
+        ).antisymmetrise(0, 1)
+        - 0.5 * einsum("ijab,kb->ijka", t2, rho_bar)
+    )
+    g2a.oovv = 0.5 * (
+        - 1.0 * t2 - 1.0 * tD - 4.0 * t2bar
+        - 1.0 * 2.0 * (
+            + 1.0 * einsum("ijbc,ac->ijab", tD + t2, g1a_adc1.vv)
+        ).antisymmetrise(2, 3)
+        + 1.0 * 2.0 * (
+            + 1.0 * einsum("jkab,ik->ijab", tD + t2, g1a_adc1.oo)
+        ).antisymmetrise(0, 1)
+        - 1.0 * 4.0 * (
+            + 1.0 * einsum("ia,jb->ijab", u.ph,
+                           einsum("jkbc,kc->jb", tD, u.ph) + rx)
+        ).antisymmetrise(0, 1).antisymmetrise(2, 3)
+    )
+    g2a.ovov = (
+        + 1.0 * g2a_adc1.ovov
+        - 4.0 * einsum("ikbc,jkac->iajb", u.pphh, u.pphh)
+        + 1.0 * einsum("ij,ab->iajb", rho.oo, g1a_adc1.vv)
+        + 1.0 * einsum("ab,ij->iajb", rho.vv, g1a_adc1.oo)
+        + 0.5 * 2.0 * (
+            - 4.0 * einsum("ikbc,jkac->iajb", t2, tbarD)
+            + 1.0 * einsum("ibjc,ac->iajb", tsq_ovov, g1a_adc1.vv)
+            - 1.0 * einsum("ibka,jk->iajb", tsq_ovov, g1a_adc1.oo)
+            + 1.0 * einsum("ka,ikbc,jc->iajb", u.ph, t2, rx)
+            + 1.0 * einsum("jc,ikbc,ka->iajb", u.ph, t2, rx)
+            + 1.0 * einsum("ja,cb,ic->iajb", u.ph, rho.vv, u.ph)
+            - 1.0 * einsum("ib,kj,ka->iajb", u.ph, rho.oo, u.ph)
+            - 2.0 * einsum("jckb,ic,ka->iajb", tsq_ovov, u.ph, u.ph)
+            + 0.5 * einsum("acbd,ic,jd->iajb", tsq_vvvv, u.ph, u.ph)
+            + 0.5 * einsum("ikjl,ka,lb->iajb", tsq_oooo, u.ph, u.ph)
+        ).symmetrise([(0, 2), (1, 3)])
+    )
+    g2a.ovvv = (
+        - 2.0 * einsum("ja,ijbc->iabc", u.ph, u.pphh)
+        + 1.0 * einsum("id,jkbc,jkad->iabc", u.ph, t2, u.pphh)
+        + 1.0 * 2.0 * (
+            - 2.0 * einsum("jb,ikcd,jkad->iabc", u.ph, t2, u.pphh)
+            + 1.0 * einsum("ib,jkcd,jkad->iabc", u.ph, t2, u.pphh)
+            - 1.0 * einsum("ic,ab->iabc", rho.ov, g1a_adc1.vv)
+            + 1.0 * einsum("ib,jc,ja->iabc", u.ph, rho.ov, u.ph)
+        ).antisymmetrise(2, 3)
+        - 0.5 * einsum("ijbc,ja->iabc", t2, rho_bar)
+    )
+    g2a.vvvv = (
+        + 2.0 * einsum("ijcd,ijab->abcd", u.pphh, u.pphh)
+        + 0.5 * 2.0 * (
+            + 2.0 * einsum("ijab,ijcd->abcd", tbarD, t2)
+            + 0.5 * 2.0 * (
+                - 1.0 * einsum("be,aecd->abcd", g1a_adc1.vv, tsq_vvvv)
+            ).antisymmetrise(0, 1)
+            + 1.0 * 2.0 * (
+                + 1.0 * einsum("ia,ijcd,jkbe,ke->abcd", u.ph, t2, t2, u.ph)
+            ).antisymmetrise(0, 1)
+            + 1.0 * 4.0 * (
+                + 1.0 * einsum("bd,ac->abcd", rho.vv, g1a_adc1.vv)
+                - 1.0 * einsum("idjb,ia,jc->abcd", tsq_ovov, u.ph, u.ph)
+            ).antisymmetrise(0, 1).antisymmetrise(2, 3)
+        ).symmetrise([(0, 2), (1, 3)])
+    )
+    return g1a, g2a
+
+
+### CVS ###
+
+def ampl_relaxed_dms_cvs_adc0(exci):
+    hf = exci.reference_state
+    u = exci.excitation_vector
+    g1a = OneParticleOperator(hf)
+    g2a = TwoParticleDensityMatrix(hf)
+    # g2a is not required for cvs-adc0 gradient,
+    # but expected by amplitude_relaxed_densities
+    g1a.cc = - 1.0 * einsum("Ia,Ja->IJ", u.ph, u.ph)
+    g1a.vv = + 1.0 * einsum("Ia,Ib->ab", u.ph, u.ph)
+    return g1a, g2a
+
+
+def ampl_relaxed_dms_cvs_adc1(exci):
+    hf = exci.reference_state
+    u = exci.excitation_vector
+    g1a = OneParticleOperator(hf)
+    g2a = TwoParticleDensityMatrix(hf)
+    g2a.cvcv = - 1.0 * einsum("Ja,Ib->IaJb", u.ph, u.ph)
+    g1a.cc = - 1.0 * einsum("Ia,Ja->IJ", u.ph, u.ph)
+    g1a.vv = + 1.0 * einsum("Ia,Ib->ab", u.ph, u.ph)
+
+    # Prerequisites for the OC block of the
+    # orbital response Lagrange multipliers:
+    fc = hf.fock(b.cc).diagonal()
+    fo = hf.fock(b.oo).diagonal()
+    fco = direct_sum("-j+I->jI", fc, fo).evaluate()
+    g1a.co = - 1.0 * einsum('JbKc,ibKc->Ji', g2a.cvcv, hf.ovcv) / fco
+    return g1a, g2a
+
+
+def ampl_relaxed_dms_cvs_adc2(exci):
+    hf = exci.reference_state
+    mp = exci.ground_state
+    u = exci.excitation_vector
+    g1a = OneParticleOperator(hf)
+    g2a = TwoParticleDensityMatrix(hf)
+
+    # Determine the t-amplitudes and multipliers:
+    t2oovv = mp.t2(b.oovv)
+    t2ccvv = mp.t2(b.ccvv)
+    t2ocvv = mp.t2(b.ocvv)
+    g1a_cvs0, g2a_cvs0 = ampl_relaxed_dms_cvs_adc0(exci)
+    t2bar = t2bar_oovv_cvs_adc2(exci, g1a_cvs0).evaluate()
+
+    g1a.cc = (
+        - 1.0 * einsum("Ia,Ja->IJ", u.ph, u.ph)
+        - 1.0 * einsum("kJba,kIba->IJ", u.pphh, u.pphh)
+        - 0.5 * einsum('IKab,JKab->IJ', t2ccvv, t2ccvv)
+        - 0.5 * einsum('kIab,kJab->IJ', t2ocvv, t2ocvv)
+    )
+
+    g1a.oo = (
+        - 1.0 * einsum("jKba,iKba->ij", u.pphh, u.pphh)
+        - 2.0 * einsum("ikab,jkab->ij", t2bar, t2oovv).symmetrise((0, 1))
+        - 0.5 * einsum('iKab,jKab->ij', t2ocvv, t2ocvv)
+        - 0.5 * einsum('ikab,jkab->ij', t2oovv, t2oovv)
+    )
+
+    # Pre-requisites for the OC block of the
+    # orbital response Lagrange multipliers:
+    fc = hf.fock(b.cc).diagonal()
+    fo = hf.fock(b.oo).diagonal()
+    fco = direct_sum("-j+I->jI", fc, fo).evaluate()
+
+    g1a.vv = (
+        + 1.0 * einsum("Ia,Ib->ab", u.ph, u.ph)
+        + 2.0 * einsum('jIcb,jIca->ab', u.pphh, u.pphh)
+        + 2.0 * einsum('ijac,ijbc->ab', t2bar, t2oovv).symmetrise((0, 1))
+        + 0.5 * einsum('IJac,IJbc->ab', t2ccvv, t2ccvv)
+        + 0.5 * einsum('ijac,ijbc->ab', t2oovv, t2oovv)
+        + 1.0 * einsum('iJac,iJbc->ab', t2ocvv, t2ocvv)
+    )
+
+    g2a.cvcv = (
+        - einsum("Ja,Ib->IaJb", u.ph, u.ph)
+    )
+
+    # The factor 1/sqrt(2) is needed because of the scaling used in adcc
+    # for the ph-pphh blocks.
+    g2a.occv = (1 / sqrt(2)) * (
+        2.0 * einsum('Ib,kJba->kJIa', u.ph, u.pphh)
+    )
+
+    g2a.oovv = (
+        + 1.0 * einsum('ijcb,ca->ijab', t2oovv, g1a_cvs0.vv).antisymmetrise((2, 3))
+        - 1.0 * t2oovv
+        - 2.0 * t2bar
+    )
+
+    # The factor 2/sqrt(2) is necessary because of the way
+    # that the ph-pphh is scaled.
+    g2a.ovvv = (2 / sqrt(2)) * (
+        einsum('Ja,iJcb->iabc', u.ph, u.pphh)
+    )
+
+    g2a.ccvv = - 1.0 * t2ccvv
+    g2a.ocvv = - 1.0 * t2ocvv
+
+    # This is the OC block of the orbital response
+    # Lagrange multipliers (lambda):
+    g1a.co = (
+        - 1.0 * einsum('JbKc,ibKc->Ji', g2a.cvcv, hf.ovcv)
+        - 0.5 * einsum('JKab,iKab->Ji', g2a.ccvv, hf.ocvv)
+        + 1.0 * einsum('kJLa,ikLa->Ji', g2a.occv, hf.oocv)
+        + 0.5 * einsum('kJab,ikab->Ji', g2a.ocvv, hf.oovv)
+        - 1.0 * einsum('kLJa,kLia->Ji', g2a.occv, hf.ocov)
+        + 1.0 * einsum('iKLa,JKLa->Ji', g2a.occv, hf.cccv)
+        + 0.5 * einsum('iKab,JKab->Ji', g2a.ocvv, hf.ccvv)
+        - 0.5 * einsum('ikab,kJab->Ji', g2a.oovv, hf.ocvv)
+        + 0.5 * einsum('iabc,Jabc->Ji', g2a.ovvv, hf.cvvv)
+    ) / fco
+
+    return g1a, g2a
+
+
+def ampl_relaxed_dms_cvs_adc2x(exci):
+    hf = exci.reference_state
+    mp = exci.ground_state
+    u = exci.excitation_vector
+    g1a = OneParticleOperator(hf)
+    g2a = TwoParticleDensityMatrix(hf)
+
+    # Determine the t-amplitudes and multipliers:
+    t2oovv = mp.t2(b.oovv)
+    t2ccvv = mp.t2(b.ccvv)
+    t2ocvv = mp.t2(b.ocvv)
+    g1a_cvs0, _ = ampl_relaxed_dms_cvs_adc0(exci)
+    t2bar = t2bar_oovv_cvs_adc2(exci, g1a_cvs0).evaluate()
+
+    g1a.cc = (
+        - 1.0 * einsum("Ia,Ja->IJ", u.ph, u.ph)
+        - 1.0 * einsum("kJba,kIba->IJ", u.pphh, u.pphh)
+        - 0.5 * einsum('IKab,JKab->IJ', t2ccvv, t2ccvv)
+        - 0.5 * einsum('kIab,kJab->IJ', t2ocvv, t2ocvv)
+    )
+
+    g1a.oo = (
+        - 1.0 * einsum("jKba,iKba->ij", u.pphh, u.pphh)
+        - 2.0 * einsum("ikab,jkab->ij", t2bar, t2oovv).symmetrise((0, 1))
+        - 0.5 * einsum('iKab,jKab->ij', t2ocvv, t2ocvv)
+        - 0.5 * einsum('ikab,jkab->ij', t2oovv, t2oovv)
+    )
+
+    # Pre-requisites for the OC block of the
+    # orbital response Lagrange multipliers:
+    fc = hf.fock(b.cc).diagonal()
+    fo = hf.fock(b.oo).diagonal()
+    fco = direct_sum("-j+I->jI", fc, fo).evaluate()
+
+    g1a.vv = (
+        + 1.0 * einsum("Ia,Ib->ab", u.ph, u.ph)
+        + 2.0 * einsum('jIcb,jIca->ab', u.pphh, u.pphh)
+        + 2.0 * einsum('ijac,ijbc->ab', t2bar, t2oovv).symmetrise((0, 1))
+        + 0.5 * einsum('IJac,IJbc->ab', t2ccvv, t2ccvv)
+        + 0.5 * einsum('ijac,ijbc->ab', t2oovv, t2oovv)
+        + 1.0 * einsum('iJac,iJbc->ab', t2ocvv, t2ocvv)
+    )
+
+    g2a.cvcv = (
+        - 1.0 * einsum("Ja,Ib->IaJb", u.ph, u.ph)
+        - 1.0 * einsum('kIbc,kJac->IaJb', u.pphh, u.pphh)
+        + 1.0 * einsum('kIcb,kJac->IaJb', u.pphh, u.pphh)
+    )
+
+    # The factor 1/sqrt(2) is needed because of the scaling used in adcc
+    # for the ph-pphh blocks.
+    g2a.occv = (1 / sqrt(2)) * (
+        2.0 * einsum('Ib,kJba->kJIa', u.ph, u.pphh)
+    )
+
+    g2a.oovv = (
+        + 1.0 * einsum('ijcb,ca->ijab', t2oovv, g1a_cvs0.vv).antisymmetrise((2, 3))
+        - 1.0 * t2oovv
+        - 2.0 * t2bar
+    )
+
+    # The factor 2/sqrt(2) is necessary because of
+    # the way that the ph-pphh is scaled
+    g2a.ovvv = (2 / sqrt(2)) * (
+        einsum('Ja,iJcb->iabc', u.ph, u.pphh)
+    )
+
+    g2a.ovov = 1.0 * (
+        - einsum("iKbc,jKac->iajb", u.pphh, u.pphh)
+        + einsum("iKcb,jKac->iajb", u.pphh, u.pphh)
+    )
+
+    g2a.ccvv = - 1.0 * t2ccvv
+    g2a.ocvv = - 1.0 * t2ocvv
+    g2a.ococ = 1.0 * einsum("iJab,kLab->iJkL", u.pphh, u.pphh)
+    g2a.vvvv = 2.0 * einsum("iJcd,iJab->abcd", u.pphh, u.pphh)
+
+    g1a.co = (
+        - 1.0 * einsum('JbKc,ibKc->Ji', g2a.cvcv, hf.ovcv)
+        - 0.5 * einsum('JKab,iKab->Ji', g2a.ccvv, hf.ocvv)
+        + 1.0 * einsum('kJLa,ikLa->Ji', g2a.occv, hf.oocv)
+        + 0.5 * einsum('kJab,ikab->Ji', g2a.ocvv, hf.oovv)
+        - 1.0 * einsum('kLJa,kLia->Ji', g2a.occv, hf.ocov)
+        + 1.0 * einsum('iKLa,JKLa->Ji', g2a.occv, hf.cccv)
+        + 0.5 * einsum('iKab,JKab->Ji', g2a.ocvv, hf.ccvv)
+        - 0.5 * einsum('ikab,kJab->Ji', g2a.oovv, hf.ocvv)
+        + 0.5 * einsum('iabc,Jabc->Ji', g2a.ovvv, hf.cvvv)
+        + 1.0 * einsum('kJmL,ikmL->Ji', g2a.ococ, hf.oooc)
+        - 1.0 * einsum('iKlM,JKMl->Ji', g2a.ococ, hf.ccco)
+        + 1.0 * einsum('iakb,kbJa->Ji', g2a.ovov, hf.ovcv)
+    ) / fco
+
+    return g1a, g2a
+
+
+DISPATCH = {
+    "mp2":  ampl_relaxed_dms_mp2,
+    "adc0": ampl_relaxed_dms_adc0,
+    "adc1": ampl_relaxed_dms_adc1,
+    "adc2": ampl_relaxed_dms_adc2,
+    "adc2x": ampl_relaxed_dms_adc2x,
+    "adc3": ampl_relaxed_dms_adc3,
+    "cvs-adc0": ampl_relaxed_dms_cvs_adc0,
+    "cvs-adc1": ampl_relaxed_dms_cvs_adc1,
+    "cvs-adc2": ampl_relaxed_dms_cvs_adc2,
+    "cvs-adc2x": ampl_relaxed_dms_cvs_adc2x,
+}
+
+
+def amplitude_relaxed_densities(excitation_or_mp):
+    """Computation of amplitude-relaxed one- and two-particle density matrices
+
+    Parameters
+    ----------
+    excitation_or_mp : LazyMp, Excitation
+        Data for which the densities are requested, either LazyMp for ground
+        state densities or Excitation for excited state densities
+
+    Returns
+    -------
+    (OneParticleOperator, TwoParticleDensityMatrix)
+        Tuple of amplitude-relaxed one- and two-particle density matrices
+
+    Raises
+    ------
+    NotImplementedError
+        if density matrices are not implemented for a given method
+    """
+    if isinstance(excitation_or_mp, LazyMp):
+        method_name = "mp2"
+    elif isinstance(excitation_or_mp, Excitation):
+        method_name = excitation_or_mp.method.name
+    if method_name not in DISPATCH:
+        raise NotImplementedError("Amplitude response is not "
+                                  f"implemented for {method_name}.")
+    g1a, g2a = DISPATCH[method_name](excitation_or_mp)
+    return evaluate(g1a), evaluate(g2a)
diff --git a/adcc/gradients/orbital_response.py b/adcc/gradients/orbital_response.py
new file mode 100644
index 00000000..37573d99
--- /dev/null
+++ b/adcc/gradients/orbital_response.py
@@ -0,0 +1,341 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2021 by the adcc authors
+##
+## This file is part of adcc.
+##
+## adcc is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published
+## by the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## adcc is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with adcc. If not, see <http://www.gnu.org/licenses/>.
+##
+## ---------------------------------------------------------------------
+import adcc.block as b
+from adcc.OneParticleOperator import OneParticleOperator
+from adcc.functions import direct_sum, einsum, evaluate
+from adcc.AmplitudeVector import AmplitudeVector
+from adcc.solver.conjugate_gradient import conjugate_gradient, default_print
+
+
+def orbital_response_rhs(hf, g1a, g2a):
+    """
+    Build the right-hand side for solving the orbital
+    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)
+            + 2.0 * einsum("icab,bc->ia", hf.ovvv, g1a.vv)
+            + 2.0 * einsum('kija,jk->ia', hf.ooov, g1a.oo)
+            + 2.0 * einsum("kiJa,Jk->ia", hf.oocv, g1a.co)
+            - 2.0 * einsum("iJka,Jk->ia", hf.ocov, g1a.co)
+            + 2.0 * einsum("iJKb,JaKb->ia", hf.occv, g2a.cvcv)  # 2PDMs start
+            + 2.0 * einsum('lKJa,lKiJ->ia', g2a.occv, hf.ococ)
+            + 1.0 * einsum('jabc,ijbc->ia', g2a.ovvv, hf.oovv)
+            - 1.0 * einsum('JKab,JKib->ia', g2a.ccvv, hf.ccov)
+            - 2.0 * einsum('jKab,jKib->ia', g2a.ocvv, hf.ocov)
+            - 1.0 * einsum('jkab,jkib->ia', g2a.oovv, hf.ooov)
+            - 2.0 * einsum('jcab,ibjc->ia', g2a.ovvv, hf.ovov)
+            - 2.0 * einsum('iJKb,JaKb->ia', g2a.occv, hf.cvcv)
+            - 1.0 * einsum('iJbc,Jabc->ia', g2a.ocvv, hf.cvvv)
+            - 1.0 * einsum('ijbc,jabc->ia', g2a.oovv, hf.ovvv)
+            + 1.0 * einsum('ibcd,abcd->ia', g2a.ovvv, hf.vvvv)
+            - 1.0 * einsum('abcd,ibcd->ia', g2a.vvvv, hf.ovvv)  # cvs-adc2x
+            + 2.0 * einsum('jakb,ijkb->ia', g2a.ovov, hf.ooov)  # cvs-adc2x
+            + 2.0 * einsum('ibjc,jcab->ia', g2a.ovov, hf.ovvv)  # cvs-adc2x
+            - 2.0 * einsum('iJkL,kLJa->ia', g2a.ococ, hf.occv)  # cvs-adc2x
+        )
+
+        ret_cv = -1.0 * (
+            + 2.0 * einsum("JIKa,JK->Ia", hf.cccv, g1a.cc)
+            + 2.0 * einsum("Icab,bc->Ia", hf.cvvv, g1a.vv)
+            + 2.0 * einsum('kIja,jk->Ia', hf.ocov, g1a.oo)
+            + 2.0 * einsum("kIJa,Jk->Ia", hf.occv, g1a.co)
+            + 2.0 * einsum("JIka,Jk->Ia", hf.ccov, g1a.co)
+            + 2.0 * einsum("IJKb,JaKb->Ia", hf.cccv, g2a.cvcv)  # 2PDMs start
+            + 2.0 * einsum("Jcab,IbJc->Ia", hf.cvvv, g2a.cvcv)
+            + 2.0 * einsum('lKJa,lKIJ->Ia', g2a.occv, hf.occc)
+            - 1.0 * einsum('jabc,jIbc->Ia', g2a.ovvv, hf.ocvv)
+            - 1.0 * einsum('JKab,JKIb->Ia', g2a.ccvv, hf.cccv)
+            - 2.0 * einsum('jKab,jKIb->Ia', g2a.ocvv, hf.occv)
+            - 1.0 * einsum('jkab,jkIb->Ia', g2a.oovv, hf.oocv)
+            - 2.0 * einsum('jcab,jcIb->Ia', g2a.ovvv, hf.ovcv)
+            - 1.0 * einsum('IJbc,Jabc->Ia', g2a.ccvv, hf.cvvv)
+            + 2.0 * einsum('jIKb,jaKb->Ia', g2a.occv, hf.ovcv)
+            + 1.0 * einsum('jIbc,jabc->Ia', g2a.ocvv, hf.ovvv)
+            + 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
+        )
+        ret = AmplitudeVector(cv=ret_cv, ov=ret_ov)
+    else:
+        # equal to the ov block of the energy-weighted density
+        # matrix when lambda_ov multipliers are zero
+        w_ov = 0.5 * (
+            + 1.0 * einsum("ijkl,klja->ia", hf.oooo, g2a.ooov)  # not in cvs-adc
+            - 1.0 * einsum("ibcd,abcd->ia", hf.ovvv, g2a.vvvv)
+            - 1.0 * einsum("jkib,jkab->ia", hf.ooov, g2a.oovv)
+            + 2.0 * einsum("ijkb,jakb->ia", hf.ooov, g2a.ovov)
+            + 1.0 * einsum("ijbc,jabc->ia", hf.oovv, g2a.ovvv)
+            - 2.0 * einsum("ibjc,jcab->ia", hf.ovov, g2a.ovvv)
+        )
+        w_ov = w_ov.evaluate()
+
+        ret_ov = -1.0 * (
+            2.0 * w_ov
+            - 1.0 * einsum("klja,ijkl->ia", hf.ooov, g2a.oooo)
+            + 1.0 * einsum("abcd,ibcd->ia", hf.vvvv, g2a.ovvv)
+            - 2.0 * einsum("jakb,ijkb->ia", hf.ovov, g2a.ooov)  # not in cvs-adc
+            + 1.0 * einsum("jkab,jkib->ia", hf.oovv, g2a.ooov)  # not in cvs-adc
+            + 2.0 * einsum("jcab,ibjc->ia", hf.ovvv, g2a.ovov)
+            - 1.0 * einsum("jabc,ijbc->ia", hf.ovvv, g2a.oovv)
+            + 2.0 * einsum("jika,jk->ia", hf.ooov, g1a.oo)
+            + 2.0 * einsum("icab,bc->ia", hf.ovvv, g1a.vv)
+        )
+        ret = AmplitudeVector(ov=ret_ov)
+
+    return ret
+
+
+def energy_weighted_density_matrix(hf, g1o, g2a):
+    if hf.has_core_occupied_space:
+        w = OneParticleOperator(hf)
+        w.cc = - 0.5 * (
+            + einsum("JKab,IKab->IJ", g2a.ccvv, hf.ccvv)
+            + einsum("kJab,kIab->IJ", g2a.ocvv, hf.ocvv)
+        )
+        w.cc += (
+            - hf.fcc
+            - einsum("IK,JK->IJ", g1o.cc, hf.fcc)
+            - einsum("ab,IaJb->IJ", g1o.vv, hf.cvcv)
+            - einsum("KL,IKJL->IJ", g1o.cc, hf.cccc)
+            - einsum("kl,kIlJ->IJ", g1o.oo, hf.ococ)
+            + einsum("ka,kJIa->IJ", g1o.ov, hf.occv)
+            + einsum("ka,kIJa->IJ", g1o.ov, hf.occv)
+            + einsum("Ka,KJIa->IJ", g1o.cv, hf.cccv)
+            + einsum("Ka,KIJa->IJ", g1o.cv, hf.cccv)
+            - einsum("Lk,kILJ->IJ", g1o.co, hf.occc)
+            - einsum("Lk,kJLI->IJ", g1o.co, hf.occc)
+            - einsum("JbKc,IbKc->IJ", g2a.cvcv, hf.cvcv)
+            - einsum("kJLa,kILa->IJ", g2a.occv, hf.occv)
+            - einsum("kLJa,kLIa->IJ", g2a.occv, hf.occv)
+            - einsum("kJmL,kImL->IJ", g2a.ococ, hf.ococ)  # cvs-adc2x
+        )
+        w.oo = - 0.5 * (
+            + einsum("jKab,iKab->ij", g2a.ocvv, hf.ocvv)
+            + einsum("jkab,ikab->ij", g2a.oovv, hf.oovv)
+            + einsum("jabc,iabc->ij", g2a.ovvv, hf.ovvv)
+        )
+        w.oo += (
+            - hf.foo
+            - einsum("ij,ii->ij", g1o.oo, hf.foo)
+            - einsum("KL,iKjL->ij", g1o.cc, hf.ococ)
+            - einsum("ab,iajb->ij", g1o.vv, hf.ovov)
+            - einsum("kl,ikjl->ij", g1o.oo, hf.oooo)
+            - einsum("ka,ikja->ij", g1o.ov, hf.ooov)
+            - einsum("ka,jkia->ij", g1o.ov, hf.ooov)
+            - einsum("Ka,iKja->ij", g1o.cv, hf.ocov)
+            - einsum("Ka,jKia->ij", g1o.cv, hf.ocov)
+            + einsum("Lk,kijL->ij", g1o.co, hf.oooc)
+            + einsum("Lk,kjiL->ij", g1o.co, hf.oooc)
+            - einsum("jKLa,iKLa->ij", g2a.occv, hf.occv)
+            - einsum("jKlM,iKlM->ij", g2a.ococ, hf.ococ)  # cvs-adc2x
+            - einsum("jakb,iakb->ij", g2a.ovov, hf.ovov)  # cvs-adc2x
+        )
+        w.vv = - 0.5 * (
+            + einsum("ibcd,iacd->ab", g2a.ovvv, hf.ovvv)
+            + einsum("IJbc,IJac->ab", g2a.ccvv, hf.ccvv)
+            + einsum("ijbc,ijac->ab", g2a.oovv, hf.oovv)
+            + einsum("bcde,acde->ab", g2a.vvvv, hf.vvvv)  # cvs-adc2x
+        )
+        w.vv += (
+            - einsum("ac,cb->ab", g1o.vv, hf.fvv)
+            - einsum('JbKc,JaKc->ab', g2a.cvcv, hf.cvcv)
+            - einsum("jKIb,jKIa->ab", g2a.occv, hf.occv)
+            - einsum("idbc,idac->ab", g2a.ovvv, hf.ovvv)
+            - einsum("iJbc,iJac->ab", g2a.ocvv, hf.ocvv)
+            - einsum("ibjc,iajc->ab", g2a.ovov, hf.ovov)  # cvs-adc2x
+        )
+        w.co = 0.5 * (
+            - einsum("JKab,iKab->Ji", g2a.ccvv, hf.ocvv)
+            + einsum("kJab,ikab->Ji", g2a.ocvv, hf.oovv)
+        )
+        w.co += (
+            + einsum("KL,iKLJ->Ji", g1o.cc, hf.occc)
+            - einsum("ab,iaJb->Ji", g1o.vv, hf.ovcv)
+            - einsum("kl,kilJ->Ji", g1o.oo, hf.oooc)
+            - einsum("Ka,iKJa->Ji", g1o.cv, hf.occv)
+            - einsum("Ka,JKia->Ji", g1o.cv, hf.ccov)
+            - einsum("ka,ikJa->Ji", g1o.ov, hf.oocv)
+            - einsum("ka,Jkia->Ji", g1o.ov, hf.coov)
+            - einsum("Lk,kiLJ->Ji", g1o.co, hf.oocc)
+            + einsum("Lk,iLkJ->Ji", g1o.co, hf.ococ)
+            - einsum("Ik,jk->Ij", g1o.co, hf.foo)
+            - einsum("JbKc,ibKc->Ji", g2a.cvcv, hf.ovcv)
+            + einsum("kJLa,ikLa->Ji", g2a.occv, hf.oocv)
+            - einsum("kLJa,kLia->Ji", g2a.occv, hf.ocov)
+            + einsum("kJmL,ikmL->Ji", g2a.ococ, hf.oooc)  # cvs-adc2x
+        )
+        w.ov = 0.5 * (
+            + einsum("jabc,ijbc->ia", g2a.ovvv, hf.oovv)
+            - einsum("JKab,JKib->ia", g2a.ccvv, hf.ccov)
+            - einsum("jkab,jkib->ia", g2a.oovv, hf.ooov)
+            - einsum("abcd,ibcd->ia", g2a.vvvv, hf.ovvv)  # cvs-adc2x
+        )
+        w.ov += (
+            - einsum("ij,ja->ia", hf.foo, g1o.ov)
+            + einsum("JaKc,iJKc->ia", g2a.cvcv, hf.occv)
+            + einsum("kLJa,iJkL->ia", g2a.occv, hf.ococ)
+            - einsum("jKab,jKib->ia", g2a.ocvv, hf.ocov)
+            - einsum("jcab,ibjc->ia", g2a.ovvv, hf.ovov)
+            + einsum("jakb,ijkb->ia", g2a.ovov, hf.ooov)  # cvs-adc2x
+        )
+        w.cv = - 0.5 * (
+            + einsum("jabc,jIbc->Ia", g2a.ovvv, hf.ocvv)
+            + einsum("JKab,JKIb->Ia", g2a.ccvv, hf.cccv)
+            + einsum("jkab,jkIb->Ia", g2a.oovv, hf.oocv)
+            + einsum("abcd,Ibcd->Ia", g2a.vvvv, hf.cvvv)  # cvs-adc2x
+        )
+        w.cv += (
+            - einsum("IJ,Ja->Ia", hf.fcc, g1o.cv)
+            + einsum("JaKc,IJKc->Ia", g2a.cvcv, hf.cccv)
+            + einsum("lKJa,lKIJ->Ia", g2a.occv, hf.occc)
+            - einsum("kJab,kJIb->Ia", g2a.ocvv, hf.occv)
+            - einsum("jcab,jcIb->Ia", g2a.ovvv, hf.ovcv)
+            - einsum("jakb,jIkb->Ia", g2a.ovov, hf.ocov)  # cvs-adc2x
+        )
+    else:
+        gi_oo = -0.5 * (
+            + 1.0 * einsum("jklm,iklm->ij", hf.oooo, g2a.oooo)
+            + 1.0 * einsum("jabc,iabc->ij", hf.ovvv, g2a.ovvv)
+            + 1.0 * einsum("klja,klia->ij", hf.ooov, g2a.ooov)
+            + 2.0 * einsum("jkla,ikla->ij", hf.ooov, g2a.ooov)
+            + 1.0 * einsum("jkab,ikab->ij", hf.oovv, g2a.oovv)
+            + 2.0 * einsum("jakb,iakb->ij", hf.ovov, g2a.ovov)
+        )
+        gi_vv = -0.5 * (
+            + 1.0 * einsum("kjib,kjia->ab", hf.ooov, g2a.ooov)
+            + einsum("bcde,acde->ab", hf.vvvv, g2a.vvvv)
+            + 1.0 * einsum("ijcb,ijca->ab", hf.oovv, g2a.oovv)
+            + 2.0 * einsum("jcib,jcia->ab", hf.ovov, g2a.ovov)
+            + 1.0 * einsum("ibcd,iacd->ab", hf.ovvv, g2a.ovvv)
+            + 2.0 * einsum("idcb,idca->ab", hf.ovvv, g2a.ovvv)
+        )
+        gi_oo = gi_oo.evaluate()
+        gi_vv = gi_vv.evaluate()
+        w = OneParticleOperator(hf)
+        w.ov = 0.5 * (
+            - 2.0 * einsum("ij,ja->ia", hf.foo, g1o.ov)
+            + 1.0 * einsum("ijkl,klja->ia", hf.oooo, g2a.ooov)
+            - 1.0 * einsum("ibcd,abcd->ia", hf.ovvv, g2a.vvvv)
+            - 1.0 * einsum("jkib,jkab->ia", hf.ooov, g2a.oovv)
+            + 2.0 * einsum("ijkb,jakb->ia", hf.ooov, g2a.ovov)
+            + 1.0 * einsum("ijbc,jabc->ia", hf.oovv, g2a.ovvv)
+            - 2.0 * einsum("ibjc,jcab->ia", hf.ovov, g2a.ovvv)
+        )
+        w.oo = (
+            + gi_oo - hf.foo
+            - einsum("ik,jk->ij", g1o.oo, hf.foo)
+            - einsum("ikjl,kl->ij", hf.oooo, g1o.oo)
+            - einsum("ikja,ka->ij", hf.ooov, g1o.ov)
+            - einsum("jkia,ka->ij", hf.ooov, g1o.ov)
+            - einsum("jaib,ab->ij", hf.ovov, g1o.vv)
+        )
+        w.vv = gi_vv - einsum("ac,cb->ab", g1o.vv, hf.fvv)
+    return evaluate(w)
+
+
+class OrbitalResponseMatrix:
+    def __init__(self, hf):
+        self.hf = hf
+
+    @property
+    def shape(self):
+        no1 = self.hf.mospaces.n_orbs(b.o)
+        if self.hf.has_core_occupied_space:
+            no1 += self.hf.mospaces.n_orbs(b.c)
+        nv1 = self.hf.mospaces.n_orbs(b.v)
+        size = no1 * nv1
+        return (size, size)
+
+    def __matmul__(self, lam):
+        ret_ov = (
+            + einsum("ab,ib->ia", self.hf.fvv, lam.ov)
+            - einsum("ij,ja->ia", self.hf.foo, lam.ov)
+            + einsum("ijab,jb->ia", self.hf.oovv, lam.ov)
+            - einsum("ibja,jb->ia", self.hf.ovov, lam.ov)
+        )
+        if self.hf.has_core_occupied_space:
+            ret_ov += (
+                + einsum("iJab,Jb->ia", self.hf.ocvv, lam.cv)
+                - einsum("ibJa,Jb->ia", self.hf.ovcv, lam.cv)
+            )
+            ret_cv = (
+                + einsum("ab,Ib->Ia", self.hf.fvv, lam.cv)
+                - einsum("IJ,Ja->Ia", self.hf.fcc, lam.cv)
+                + einsum("IJab,Jb->Ia", self.hf.ccvv, lam.cv)
+                - einsum("IbJa,Jb->Ia", self.hf.cvcv, lam.cv)
+                + einsum("Ijab,jb->Ia", self.hf.covv, lam.ov)
+                - einsum("Ibja,jb->Ia", self.hf.cvov, lam.ov)
+            )
+            ret = AmplitudeVector(cv=ret_cv, ov=ret_ov)
+        else:
+            ret = AmplitudeVector(ov=ret_ov)
+        return evaluate(ret)
+
+
+class OrbitalResponsePinv:
+    def __init__(self, hf):
+        self.hf = hf
+        # Terms common to adc and cvs-adc
+        fo = hf.fock(b.oo).diagonal()
+        fv = hf.fock(b.vv).diagonal()
+        fov = direct_sum("-i+a->ia", fo, fv)
+
+        if hf.has_core_occupied_space:
+            fc = hf.fock(b.cc).diagonal()
+            fcv = direct_sum("-I+a->Ia", fc, fv)
+            self.df = AmplitudeVector(cv=fcv, ov=fov)
+        else:
+            self.df = AmplitudeVector(ov=fov)
+        self.df.evaluate()
+
+    @property
+    def shape(self):
+        no1 = self.hf.mospaces.n_orbs(b.o)
+        if self.hf.has_core_occupied_space:
+            no1 += self.hf.mospaces.n_orbs(b.c)
+        nv1 = self.hf.mospaces.n_orbs(b.v)
+        size = no1 * nv1
+        return (size, size)
+
+    def __matmul__(self, invec):
+        return invec / self.df
+
+
+def orbital_response(hf, rhs):
+    """
+    Solves the orbital response equations
+    for a given reference state and right-hand side
+    """
+    # TODO: pass solver arguments
+    A = OrbitalResponseMatrix(hf)
+    Pinv = OrbitalResponsePinv(hf)
+    x0 = (Pinv @ rhs).evaluate()
+    lam = conjugate_gradient(A, rhs=rhs, x0=x0, Pinv=Pinv,
+                             explicit_symmetrisation=None,
+                             callback=default_print)
+    return lam.solution
diff --git a/adcc/solver/conjugate_gradient.py b/adcc/solver/conjugate_gradient.py
index 77982cc3..86bbc4eb 100644
--- a/adcc/solver/conjugate_gradient.py
+++ b/adcc/solver/conjugate_gradient.py
@@ -143,7 +143,7 @@ def is_converged(state):
     state.solution = x0
     state.residual = evaluate(rhs - matrix @ state.solution)
     state.n_applies += 1
-    state.residual_norm = np.sqrt(state.residual @ state.residual)
+    state.residual_norm = np.sqrt(state.residual.dot(state.residual))
     pk = zk = Pinv @ state.residual
 
     if explicit_symmetrisation:
@@ -166,7 +166,7 @@ def is_converged(state):
 
         residual_old = state.residual
         state.residual = evaluate(residual_old - ak * Apk)
-        state.residual_norm = np.sqrt(state.residual @ state.residual)
+        state.residual_norm = np.sqrt(state.residual.dot(state.residual))
 
         callback(state, "next_iter")
         if is_converged(state):
diff --git a/adcc/tests/ExcitedStates_test.py b/adcc/tests/ExcitedStates_test.py
index 7e28ab6d..1be10a0a 100644
--- a/adcc/tests/ExcitedStates_test.py
+++ b/adcc/tests/ExcitedStates_test.py
@@ -43,7 +43,8 @@ def test_excitation_view(self):
                 if key.startswith("_"):
                     continue
                 blacklist = ["__", "index", "_ao", "excitation_vector",
-                             "method", "parent_state"]
+                             "method", "parent_state", "ground_state",
+                             "reference_state"]
                 if any(b in key for b in blacklist):
                     continue
                 try:
diff --git a/adcc/tests/backends/testing.py b/adcc/tests/backends/testing.py
index 0d3e483c..fcc1c4aa 100644
--- a/adcc/tests/backends/testing.py
+++ b/adcc/tests/backends/testing.py
@@ -191,6 +191,7 @@ def cached_backend_hf(backend: str, system: str, conv_tol=1e-12, pe_options=None
     and return the result.
     """
     import adcc.backends
+    print(system)
 
     from .. import testcases
 
diff --git a/adcc/tests/data/gradient_data.json b/adcc/tests/data/gradient_data.json
new file mode 100644
index 00000000..a74f95e5
--- /dev/null
+++ b/adcc/tests/data/gradient_data.json
@@ -0,0 +1 @@
+{"basissets": ["sto3g", "ccpvdz"], "methods": ["mp2", "adc0", "adc1", "adc2", "adc2x", "adc3", "cvs-adc0", "cvs-adc1", "cvs-adc2", "cvs-adc2x"], "molecules": ["h2o"], "h2o": {"sto3g": {"mp2": {"energy": -74.99357808910258, "gradient": [[0.09648089353868272, -7.275957614183426e-12, 0.06886449053126853], [-0.026183461559412535, 0.0, -0.0659738619287964], [-0.07029743172461167, -1.4551915228366852e-11, -0.002890627962187864]], "config": {"conv_tol": 1e-09, "core_orbitals": null}}, "adc0": {"energy": [-73.96883482604426, -73.91481187944933, -73.8057558947831, -73.75173294818818, -73.72721279699914], "gradient": [[[0.21500520807603607, -1.0186340659856796e-10, 0.15122539465664886], [0.074898492108332, 0.0, -0.33320423345139716], [-0.2899037002425757, -2.1827872842550278e-11, 0.18197883959510364]], [[0.15869260571344057, -7.275957614183426e-12, 0.11135792801360367], [0.12612935141078196, -8.731149137020111e-11, -0.3458573067910038], [-0.28482195712422254, -7.275957614183426e-12, 0.23449937959230738]], [[0.4883978163416032, -2.1827872842550278e-11, 0.34780997568304883], [-0.09521430340100778, 7.275957614183426e-12, -0.38596352970489534], [-0.393183512838732, -2.9103830456733704e-11, 0.03815355485130567]], [[0.432085214095423, 2.1827872842550278e-11, 0.3079425090982113], [-0.04398344415676547, -8.731149137020111e-11, -0.39861660305177793], [-0.3881017696039635, -7.275957614183426e-12, 0.09067409489944112]], [[0.4312687490892131, -1.4551915228366852e-11, 0.30399765190668404], [0.0012444719613995403, -2.1827872842550278e-11, -0.45833798522653524], [-0.4325132208832656, -7.275957614183426e-12, 0.15434033421479398]]], "config": {"conv_tol": 1e-09, "n_singlets": 5, "core_orbitals": null}}, "adc1": {"energy": [-74.47588319597416, -74.38511885048102, -74.3571822920938, -74.24904200360893, -74.13327522706388], "gradient": [[[0.2519565781985875, -2.1827872842550278e-11, 0.17643765273533063], [0.043837964338308666, 1.4551915228366852e-11, -0.32756816069741035], [-0.29579454263148364, 2.1827872842550278e-11, 0.15113050883519463]], [[0.44204020009055967, 0.0, 0.31597312643862097], [-0.08205663987610023, 0.0, -0.35633136297838064], [-0.3599835601708037, -2.1827872842550278e-11, 0.040358237114560325]], [[0.102462732247659, 2.1827872842550278e-11, 0.07140888640424237], [0.11964530191471567, 2.1827872842550278e-11, -0.2768455860496033], [-0.22210803416237468, -2.1827872842550278e-11, 0.20543670036568074]], [[0.33387671443779254, 2.1827872842550278e-11, 0.2387534842637251], [-2.251235127914697e-05, -1.4551915228366852e-11, -0.35684544259129325], [-0.33385420200647786, 1.4551915228366852e-11, 0.1180919590042322]], [[0.4270152217577561, 1.4551915228366852e-11, 0.3020655645959778], [-0.0510407698166091, 1.4551915228366852e-11, -0.38096642302843975], [-0.37597445177379996, -2.1827872842550278e-11, 0.0789008592182654]]], "config": {"conv_tol": 1e-09, "n_singlets": 5, "core_orbitals": null}}, "adc2": {"energy": [-74.5230649475857, -74.42102314285407, -74.399904734204, -74.2806092729357, -74.15388076828465], "gradient": [[[0.2729663057762082, -1.2369127944111824e-10, 0.1915359973354498], [0.043751713426900096, -7.275957614183426e-12, -0.3499776703029056], [-0.3167180191812804, 1.4551915228366852e-11, 0.1584416737605352]], [[0.4661114461196121, 3.637978807091713e-11, 0.3327183597430121], [-0.08526563976920443, 1.8917489796876907e-10, -0.3770551475681714], [-0.38084580657596234, 2.1827872842550278e-11, 0.044336788858345244]], [[0.1245402125568944, -1.4551915228366852e-11, 0.08709251005348051], [0.11524763699708274, 2.1827872842550278e-11, -0.2941215308092069], [-0.23978784963401267, -1.4551915228366852e-11, 0.20702902170887683]], [[0.34492495327140205, -1.4551915228366852e-11, 0.24646419990313007], [-0.0014078249514568597, -2.9103830456733704e-11, -0.36650604708120227], [-0.3435171281089424, -7.275957614183426e-12, 0.12004184816760244]], [[0.4314259997627232, -1.4551915228366852e-11, 0.30517971346125705], [-0.05028844868502347, 1.4551915228366852e-11, -0.386704910924891], [-0.3811375508448691, -7.275957614183426e-12, 0.08152519801660674]]], "config": {"conv_tol": 1e-09, "n_singlets": 5, "core_orbitals": null}}, "adc2x": {"energy": [-74.5484662688597, -74.44218300989021, -74.42210922368841, -74.30300532559038, -74.16527252522135], "gradient": [[[0.28895747921342263, 7.275957614183426e-12, 0.2026268877380062], [0.045065266407618765, -1.4551915228366852e-11, -0.3685835983196739], [-0.33402274580294034, -5.093170329928398e-11, 0.16595671139657497]], [[0.4846122063099756, 1.4551915228366852e-11, 0.3459748487730394], [-0.08685884039005032, 0.0, -0.39460422133561224], [-0.39775336616730783, -2.9103830456733704e-11, 0.04862937352299923]], [[0.13661488301295321, 1.4551915228366852e-10, 0.09541062240168685], [0.11530230956850573, -1.4551915228366852e-11, -0.3067889074445702], [-0.25191719264694257, 7.275957614183426e-12, 0.21137828581413487]], [[0.36053903902939055, 7.275957614183426e-12, 0.25766010332881706], [-0.0009222038715961389, 8.003553375601768e-11, -0.3839128288964275], [-0.3596168350704829, -2.1827872842550278e-11, 0.12625272644800134]], [[0.43669924540881766, 7.275957614183426e-12, 0.308943429197825], [-0.05085965209582355, 2.9103830456733704e-11, -0.3915265553150675], [-0.38583959305105964, -7.275957614183426e-12, 0.08258312672114698]]], "config": {"conv_tol": 1e-09, "n_singlets": 5, "core_orbitals": null}}, "adc3": {"energy": [-74.55090439658746, -74.44714801129354, -74.42426273712188, -74.30620148905741, -74.16694386620925], "gradient": [[[0.2908625282943831, 2.1827872842550278e-11, 0.20394485128781525], [0.04500433992507169, -1.4551915228366852e-11, -0.37048940415115794], [-0.33586686859780457, -2.9103830456733704e-11, 0.16654455379466526]], [[0.4865842864601291, -2.1827872842550278e-11, 0.3473889870510902], [-0.08571265485807089, 2.1827872842550278e-11, -0.3983368489789427], [-0.4008716317039216, -1.4551915228366852e-11, 0.05094786284462316]], [[0.13570507883559912, -2.9103830456733704e-11, 0.09476317676308099], [0.11679784433363238, 2.1827872842550278e-11, -0.3079343576464453], [-0.2525029233220266, 2.9103830456733704e-11, 0.21317118174920324]], [[0.3612387817411218, 7.275957614183426e-12, 0.2581562034683884], [0.00045972187217557803, -1.4551915228366852e-11, -0.38661062382743694], [-0.3616985033659148, -7.275957614183426e-12, 0.12845442127581919]], [[0.4404495687704184, 2.9103830456733704e-11, 0.31157616753625916], [-0.05225038512435276, 2.9103830456733704e-11, -0.39351955433085095], [-0.3881991834496148, -2.9103830456733704e-11, 0.08194338764587883]]], "config": {"conv_tol": 1e-09, "n_singlets": 5, "core_orbitals": null}}, "cvs-adc0": {"energy": [-54.123082474404995, -53.960003543143834], "gradient": [[[0.13670423482835758, 0.0, 0.09585740332113346], [0.09903088248393033, -1.4551915228366852e-11, -0.2842584849859122], [-0.2357351180180558, 7.275957614183426e-12, 0.188401082654309]], [[0.4100968432685477, 2.1827872842550278e-11, 0.2924419842674979], [-0.07108191314910073, -1.4551915228366852e-11, -0.33701778129761806], [-0.3390149306069361, -7.275957614183426e-12, 0.04457579784502741]]], "config": {"conv_tol": 1e-09, "n_singlets": 5, "core_orbitals": 1}}, "cvs-adc1": {"energy": [-54.85354413664511, -54.79419426319812], "gradient": [[[0.2044143344683107, 2.1827872842550278e-11, 0.1419889454045915], [0.033883336967846844, -1.4551915228366852e-11, -0.262222796693095], [-0.2382976719745784, 0.0, 0.12023385223437799]], [[0.3112917308098986, 7.275957614183426e-12, 0.22435430129553424], [-0.03036278657236835, 7.275957614183426e-12, -0.29155276880192105], [-0.2809289447250194, 2.1827872842550278e-11, 0.06719846829219023]]], "config": {"conv_tol": 1e-09, "n_singlets": 5, "core_orbitals": 1}}, "cvs-adc2": {"energy": [-54.98903589084722, -54.90586009813991, -53.16685681552232, -53.101104786907726, -53.00377788426118], "gradient": [[[0.23872009933984373, -2.9103830456733704e-11, 0.1671009868659894], [0.03348872094647959, 4.3655745685100555e-11, -0.29891402232897235], [-0.27220882097026333, 1.4551915228366852e-11, 0.1318130366198602]], [[0.36193002599611646, 2.1827872842550278e-11, 0.2592484456472448], [-0.04305156866030302, 0.0, -0.3264197397438693], [-0.3188784580779611, -7.275957614183426e-12, 0.06717129518801812]], [[0.305568811170815, 7.275957614183426e-12, 0.21378264181839768], [0.19212377251824364, 2.9103830456733704e-11, -0.5935662706178846], [-0.49769258473679656, 1.0913936421275139e-10, 0.3797836300363997]], [[0.23468779257382266, -7.275957614183426e-11, 0.16384246382222045], [0.24505627318285406, 1.964508555829525e-10, -0.5933986320596887], [-0.47974406644061673, -2.1827872842550278e-11, 0.4295561693725176]], [[0.5789614193781745, 2.9103830456733704e-11, 0.4103672228375217], [0.022010977067111526, 3.637978807091713e-11, -0.6463255668786587], [-0.6009723973911605, 1.0186340659856796e-10, 0.23595834522711812]]], "config": {"conv_tol": 1e-09, "n_singlets": 5, "core_orbitals": 1}}, "cvs-adc2x": {"energy": [-55.08829803906612, -54.9945312916704, -54.266903715981805, -54.258055435885886, -54.20500159732517], "gradient": [[[0.27884997078217566, -1.673470251262188e-10, 0.1953710336747463], [0.03098650036554318, 7.275957614183426e-12, -0.33784350618952885], [-0.3098364717588993, 1.3096723705530167e-10, 0.14247247358434834]], [[0.41212538530817255, -1.4551915228366852e-11, 0.29478439913509646], [-0.05303481598821236, -1.2369127944111824e-10, -0.36559760699310573], [-0.359090569872933, 0.0, 0.07081320873840014]], [[0.6789990236284211, -1.7462298274040222e-10, 0.4736056833571638], [-0.019374826908460818, 2.9103830456733704e-11, -0.6864346827715053], [-0.6596241974912118, 6.548361852765083e-11, 0.21282900056394283]], [[0.6081797262158943, 1.3096723705530167e-10, 0.4384273026807932], [0.01058207180176396, 2.9103830456733704e-11, -0.6685617413022555], [-0.6187617984833196, -5.093170329928398e-11, 0.23013443964737235]], [[0.46743476379924687, 1.4551915228366852e-11, 0.33145697056897916], [0.022216292971279472, 0.0, -0.528248111018911], [-0.489651057854644, 1.4551915228366852e-11, 0.19679114157042932]]], "config": {"conv_tol": 1e-09, "n_singlets": 5, "core_orbitals": 1}}}, "ccpvdz": {"mp2": {"energy": -76.22940338787562, "gradient": [[0.024395466447458602, -2.9103830456733704e-11, 0.01773811209568521], [-0.010762456666270737, -2.9103830456733704e-11, -0.011150370693940204], [-0.013632937501824927, -1.4551915228366852e-11, -0.006587735537323169]], "config": {"conv_tol": 1e-09, "core_orbitals": null}}, "adc0": {"energy": [-75.34687593720773, -75.28099897309124, -75.27790714495661, -75.21203018084012, -75.1263800402113], "gradient": [[[0.05781977490551071, 7.275957614183426e-12, 0.04124122217763215], [-0.0010635846047080122, -1.4551915228366852e-11, -0.06019881294196239], [-0.05674794071819633, 0.0, 0.018958910419314634]], [[0.013119906550855376, -2.1827872842550278e-11, 0.009577973985869903], [0.0517213462953805, 2.1827872842550278e-11, -0.08735985546081793], [-0.06483396996191004, 2.9103830456733704e-11, 0.0777820943549159]], [[0.035872519212716725, 1.4551915228366852e-11, 0.02590936847263947], [0.007818394005880691, 7.275957614183426e-12, -0.04966096833959455], [-0.04368333833554061, 2.9103830456733704e-11, 0.023752558598062024]], [[-0.008827349083730951, -1.4551915228366852e-11, -0.0057538798355381005], [0.060603324898693245, -7.275957614183426e-12, -0.07682201085845008], [-0.051769367564702407, -1.4551915228366852e-11, 0.08257574253366329]], [[0.2704858468277962, 1.4551915228366852e-11, 0.19159464181575458], [-0.07004997117473977, 0.0, -0.18824054856668226], [-0.2004280803448637, -1.4551915228366852e-11, -0.0033529280190123245]]], "config": {"conv_tol": 1e-09, "n_singlets": 5, "core_orbitals": null}}, "adc1": {"energy": [-75.68543998430039, -75.61931016170234, -75.59868950264112, -75.532752986795, -75.45558187581646], "gradient": [[[0.09274209536670242, -7.275957614183426e-12, 0.06557503178919433], [-0.0026306448489776812, 2.1827872842550278e-11, -0.09467139691696502], [-0.0901040762255434, -1.4551915228366852e-11, 0.029097417085722554]], [[0.1093199013848789, 2.1827872842550278e-11, 0.0780909138629795], [-0.006661084160441533, 2.9103830456733704e-11, -0.10735268513963092], [-0.10265199275454506, -2.1827872842550278e-11, 0.029262722055136692]], [[0.011673270331812091, -7.275957614183426e-12, 0.00844753841374768], [0.056844679136702325, 0.0, -0.09296221228578361], [-0.06851130315772025, 2.1827872842550278e-11, 0.08451495548069943]], [[0.02629122353391722, 0.0, 0.019252847356256098], [0.057597583705501165, -1.4551915228366852e-11, -0.11000348474772181], [-0.08388262351218145, -2.9103830456733704e-11, 0.09075075774308061]], [[0.2821038434340153, 1.4551915228366852e-11, 0.19977067352738231], [-0.07840616344037699, 0.0, -0.18870985075773206], [-0.20369109181774547, 2.1827872842550278e-11, -0.011060022770834621]]], "config": {"conv_tol": 1e-09, "n_singlets": 5, "core_orbitals": null}}, "adc2": {"energy": [-75.92967540818942, -75.85499790112056, -75.84309170418786, -75.76675240926447, -75.66953798976122], "gradient": [[[0.11448709326941753, -1.4551915228366852e-11, 0.08113363845041022], [-0.008901050896383822, 2.1827872842550278e-11, -0.1090564420010196], [-0.10557860966218868, 7.275957614183426e-12, 0.027923848014324903]], [[0.1131836650747573, -7.275957614183426e-12, 0.08065145519503858], [-0.004472936961974483, -1.4551915228366852e-11, -0.11437404136813711], [-0.10870395785605069, 2.9103830456733704e-11, 0.03372336812026333]], [[0.022058329042920377, -2.1827872842550278e-11, 0.015899113190243952], [0.050983897534024436, -2.9103830456733704e-11, -0.09580092146643437], [-0.07303541761211818, 7.275957614183426e-12, 0.0799023346407921]], [[0.020308474187913816, 1.4551915228366852e-11, 0.014956989012716804], [0.058050304411153775, 1.4551915228366852e-11, -0.1042313243233366], [-0.07835257506667404, 1.4551915228366852e-11, 0.08927455494267633]], [[0.26928795556159457, -2.9103830456733704e-11, 0.19073455897159874], [-0.07746179534296971, 7.275957614183426e-12, -0.1764748514542589], [-0.1918197939376114, 2.1827872842550278e-11, -0.01425892054248834]]], "config": {"conv_tol": 1e-09, "n_singlets": 5, "core_orbitals": null}}, "adc2x": {"energy": [-75.94713027647231, -75.87072120069274, -75.86069046652048, -75.78212932285793, -75.68339216526915], "gradient": [[[0.12108681791869458, 7.275957614183426e-12, 0.08575116403517313], [-0.008930603667977266, -2.9103830456733704e-11, -0.11596714374172734], [-0.11214869188552257, 1.4551915228366852e-11, 0.030217022183933295]], [[0.12104355383053189, -7.275957614183426e-12, 0.08621679777570534], [-0.005538905352295842, 2.9103830456733704e-11, -0.1212127921753563], [-0.11549783187365392, 1.4551915228366852e-11, 0.034996749353013]], [[0.02765446234116098, 2.1827872842550278e-11, 0.01980982987151947], [0.05061737170035485, -7.275957614183426e-12, -0.10117332988738781], [-0.07826491810556035, 2.9103830456733704e-11, 0.08136406357516535]], [[0.027274335945548955, 2.1827872842550278e-11, 0.019888993163476698], [0.05681757771526463, 0.0, -0.10988471429300262], [-0.08408563704870176, -1.4551915228366852e-11, 0.08999595484056044]], [[0.2716475580527913, -7.275957614183426e-12, 0.19238663020951208], [-0.07715970138815464, -2.9103830456733704e-11, -0.17938891498488374], [-0.1944814575035707, -7.275957614183426e-12, -0.012996937817661092]]], "config": {"conv_tol": 1e-09, "n_singlets": 5, "core_orbitals": null}}, "adc3": {"energy": [-75.92857834998667, -75.85463334402809, -75.84167671842562, -75.76664279719428, -75.67470300439263], "gradient": [[[0.11627493508422049, 0.0, 0.08229840418789536], [-0.006295996819972061, 1.4551915228366852e-11, -0.11453756176342722], [-0.10997124934510794, 2.1827872842550278e-11, 0.032240211170574185]], [[0.12384516876772977, 0.0, 0.08824648732843343], [-0.006875947656226344, -7.275957614183426e-12, -0.12234318283299217], [-0.11696223153558094, -1.382431946694851e-10, 0.03409750875289319]], [[0.025811466584855225, 2.9103830456733704e-11, 0.018465745204593986], [0.0533070154779125, 2.9103830456733704e-11, -0.10298046591924503], [-0.07911143552337307, 7.275957614183426e-12, 0.08451522156246938]], [[0.033787787426263094, 1.1641532182693481e-10, 0.024522948973753955], [0.056045876095595304, -2.1827872842550278e-11, -0.11573187115573091], [-0.08982723313965835, -2.1827872842550278e-11, 0.0912091305362992]], [[0.2793064223806141, 3.637978807091713e-11, 0.19778368141851388], [-0.07754907244088827, 7.275957614183426e-12, -0.1869451557067805], [-0.2017506956035504, 2.1827872842550278e-11, -0.010837729831109755]]], "config": {"conv_tol": 1e-09, "n_singlets": 5, "core_orbitals": null}}, "cvs-adc0": {"energy": [-55.29234315848254, -55.22337436623141, -54.6671662062745, -54.64254226419146, -54.31941455257717], "gradient": [[[-0.01363594663416734, 1.4551915228366852e-11, -0.009354182511742692], [0.005705504816432949, -5.093170329928398e-11, 0.0060958039975957945], [0.007961477167555131, 1.4551915228366852e-11, 0.0032610131020192057]], [[-0.035583202239649836, 1.4551915228366852e-11, -0.02468603632587474], [0.014587483354262076, -2.9103830456733704e-11, 0.016633648599963635], [0.021026079608418513, -2.1827872842550278e-11, 0.008054661513597239]], [[0.295615929106134, -2.9103830456733704e-11, 0.20218788867350668], [-0.11422500915068667, 1.1641532182693481e-10, -0.14526737584674265], [-0.18136249485542066, 0.0, -0.05691822549124481]], [[-0.04258258258050773, -5.820766091346741e-11, -0.021944308828096837], [0.16518283695040736, 8.003553375601768e-11, -0.19658326423086692], [-0.12257139149005525, 5.820766091346741e-11, 0.2185302046345896]], [[-0.0943416317095398, -1.4551915228366852e-11, -0.06655934288573917], [0.01681761433428619, -1.4551915228366852e-11, 0.07614312684745528], [0.07755017422459787, -7.275957614183426e-12, -0.009581486701790709]]], "config": {"conv_tol": 1e-09, "n_singlets": 5, "core_orbitals": 1}}, "cvs-adc1": {"energy": [-55.7623580120961, -55.74423117760095, -55.23885575120384, -55.20048469034408, -55.178979221895716], "gradient": [[[0.11376989239943214, 3.637978807091713e-11, 0.07504696383693954], [-0.006233317595615517, -2.9103830456733704e-11, -0.10649701902730158], [-0.107507510263531, 1.0186340659856796e-10, 0.03145194291573716]], [[0.1367624182748841, 1.1641532182693481e-10, 0.10286059680947801], [-0.010134714473679196, -8.003553375601768e-11, -0.13692940795590403], [-0.12659825186710805, 7.275957614183426e-12, 0.03407134330336703]], [[-0.035857140050211456, -7.275957614183426e-12, -0.025069641618756577], [0.02807263379509095, -1.382431946694851e-10, -0.0019539316344889812], [0.00781276056659408, 0.0, 0.027025939685699996]], [[0.0637307473589317, -8.003553375601768e-11, 0.046041117268032394], [-0.04479707258724375, -7.275957614183426e-12, -0.005256572927464731], [-0.018904139367805328, 7.275957614183426e-12, -0.04078219792427262]], [[-0.05261707480531186, 0.0, -0.03676470001664711], [0.005314122914569452, 1.2369127944111824e-10, 0.04785384914430324], [0.047328871696663555, 2.9103830456733704e-11, -0.011087133760156576]]], "config": {"conv_tol": 1e-09, "n_singlets": 5, "core_orbitals": 1}}, "cvs-adc2": {"energy": [-56.463324400674544, -56.39198974445309, -55.93190966426827, -55.91365902930488, -55.634891014248154], "gradient": [[[0.06426175875822082, 1.4551915228366852e-11, 0.04550280176772503], [-0.0016873756248969585, 2.9103830456733704e-11, -0.06586654879356502], [-0.06254429007094586, 1.7462298274040222e-10, 0.020366252538224217]], [[0.059006211573432665, 1.964508555829525e-10, 0.04234478738362668], [0.005969680009002332, -2.1100277081131935e-10, -0.07167619411484338], [-0.06494627146457788, -1.382431946694851e-10, 0.02933372787811095]], [[0.26678762744995765, -1.964508555829525e-10, 0.18208895675343229], [-0.08748003323125886, -3.346940502524376e-10, -0.15278823697735788], [-0.17927927291748347, 7.275957614183426e-11, -0.029298384521098342]], [[0.019884901426848955, -1.4551915228366852e-11, 0.022045515754143707], [0.1168202597182244, -5.093170329928398e-11, -0.1942855826273444], [-0.13667679524223786, 3.637978807091713e-11, 0.17224251533480128]], [[-0.08170067649916746, -3.4924596548080444e-10, -0.05746503695991123], [0.049807636183686554, 1.4551915228366852e-10, 0.01592683904164005], [0.031919313929392956, 4.5838532969355583e-10, 0.04154028196353465]]], "config": {"conv_tol": 1e-09, "n_singlets": 5, "core_orbitals": 1}}, "cvs-adc2x": {"energy": [-56.527690992922494, -56.455320049295366, -56.00756567594842, -55.992414483676754, -55.7258963632124], "gradient": [[[0.08823665066302055, 3.2014213502407074e-10, 0.06228063011803897], [-0.0035475476979627274, 7.275957614183426e-12, -0.08849635638762265], [-0.08465851147047943, 3.2014213502407074e-10, 0.026218218845315278]], [[0.09283125843649032, 8.949427865445614e-10, 0.06637479383789469], [-0.0013974821922602132, -2.255546860396862e-10, -0.09725589540903457], [-0.09140373219997855, -1.8189894035458565e-10, 0.030883435916621238]], [[0.26091981134959497, 7.785274647176266e-10, 0.17749230707704555], [-0.08874684962211177, 2.0372681319713593e-10, -0.14432455416681478], [-0.1721437370069907, -4.147295840084553e-10, -0.03316539854131406]], [[0.0209724435408134, 5.529727786779404e-10, 0.02326315642130794], [0.11284109613188775, 4.874891601502895e-10, -0.19026156545442063], [-0.13378430448210565, -4.0745362639427185e-10, 0.1670009150839178]], [[-0.005838602111907676, 1.1423253454267979e-09, -0.003816529810137581], [0.0007158659500419162, 2.1827872842550278e-11, 0.004858776730543468], [0.005147792071511503, -6.548361852765083e-11, -0.0010405088687548414]]], "config": {"conv_tol": 1e-09, "n_singlets": 5, "core_orbitals": 1}}}, "xyz": "\n    O 0 0 0\n    H 0 0 1.795239827225189\n    H 1.693194615993441 0 -0.599043184453037\n    "}}
\ No newline at end of file
diff --git a/adcc/tests/functionality_gradients_test.py b/adcc/tests/functionality_gradients_test.py
new file mode 100644
index 00000000..2df68607
--- /dev/null
+++ b/adcc/tests/functionality_gradients_test.py
@@ -0,0 +1,93 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2020 by the adcc authors
+##
+## This file is part of adcc.
+##
+## adcc is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published
+## by the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## adcc is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with adcc. If not, see <http://www.gnu.org/licenses/>.
+##
+## ---------------------------------------------------------------------
+import unittest
+import itertools
+
+import adcc
+import adcc.backends
+
+from numpy.testing import assert_allclose
+
+import pytest
+
+from .backends.testing import cached_backend_hf
+from .testdata_cache import gradient_data
+
+
+backends = [b for b in adcc.backends.available()
+            if b not in ["molsturm", "veloxchem"]]
+
+molecules = gradient_data["molecules"]
+basissets = gradient_data["basissets"]
+methods = gradient_data["methods"]
+
+@pytest.mark.skipif(len(backends) == 0, reason="No backend found.")
+@pytest.mark.parametrize("backend", backends)
+@pytest.mark.parametrize("method", methods)
+@pytest.mark.parametrize("basis", basissets)
+@pytest.mark.parametrize("molecule", molecules)
+def test_nuclear_gradient(molecule, basis, method, backend):
+    grad_ref = gradient_data[molecule][basis][method]
+
+    energy_ref = grad_ref["energy"]
+    grad_fdiff = grad_ref["gradient"]
+    kwargs = grad_ref["config"]
+    conv_tol = kwargs["conv_tol"]
+
+    scfres = cached_backend_hf(backend, f"{molecule}_{basis}", conv_tol=1e-11)
+    print(kwargs)
+
+    if "adc" in method:
+        state = adcc.run_adc(scfres, method=method, **kwargs)
+        for ee in state.excitations:
+            grad = adcc.nuclear_gradient(ee)
+            assert_allclose(energy_ref[ee.index], ee.total_energy,
+                            atol=conv_tol)
+            # check energy computed with unrelaxed densities
+            gs_corr = 0.0
+            if ee.method.level > 0:
+                # compute the ground state contribution
+                # to the correlation energy
+                gs_energy = ee.ground_state.energy(ee.method.level)
+                gs_corr = gs_energy - ee.reference_state.energy_scf
+            assert_allclose(
+                gs_corr + ee.excitation_energy,
+                grad._energy, atol=1e-10
+            )
+            assert_allclose(
+                grad_fdiff[ee.index], grad.total, atol=5e-8, rtol=0,
+                err_msg=f'Gradient for state {ee.index} wrong.'
+            )
+    else:
+        # MP2 gradients
+        refstate = adcc.ReferenceState(scfres)
+        mp = adcc.LazyMp(refstate)
+        grad = adcc.nuclear_gradient(mp)
+        assert_allclose(energy_ref, mp.energy(2), atol=1e-8)
+        # check energy computed with unrelaxed densities
+        assert_allclose(
+            mp.energy_correction(2), grad._energy, atol=5e-8, rtol=0
+        )
+        assert_allclose(
+            grad_fdiff, grad.total, atol=1e-8, rtol=0
+        )
diff --git a/adcc/tests/generators/dump_fdiff_gradient.py b/adcc/tests/generators/dump_fdiff_gradient.py
new file mode 100644
index 00000000..1003e2d4
--- /dev/null
+++ b/adcc/tests/generators/dump_fdiff_gradient.py
@@ -0,0 +1,193 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2021 by the adcc authors
+##
+## This file is part of adcc.
+##
+## adcc is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published
+## by the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## adcc is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with adcc. If not, see <http://www.gnu.org/licenses/>.
+##
+## ---------------------------------------------------------------------
+from dataclasses import dataclass
+from pathlib import Path
+import itertools
+import adcc
+import numpy as np
+import json
+from tqdm import tqdm
+
+from pyscf import gto, lib
+
+from adcc.tests.testcases import _xyz
+
+adcc.set_n_threads(8)
+lib.num_threads(8)
+
+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
+class Molecule:
+    name: str
+    charge: int = 0
+    multiplicity: int = 1
+    core_orbitals: int = 1
+
+    @property
+    def xyz(self):
+        return _xyz[self.name][0]
+
+
+molecules = [
+    Molecule("h2o", 0, 1),
+    Molecule("ch2nh2", charge=0, multiplicity=2),
+]
+
+
+def _molstring(elems, coords):
+    s = ""
+    for kk, (i, c) in enumerate(zip(elems, coords)):
+        s += f"{i} {c[0]} {c[1]} {c[2]}"
+        if kk != len(elems) - 1:
+            s += "\n"
+    return s
+
+
+def adc_energy(scfres, method, **kwargs):
+    state = adcc.run_adc(method=method, data_or_matrix=scfres,
+                         output=None,
+                         **kwargs)
+    return state.total_energy
+
+
+def mp_energy(scfres, method, **kwargs):
+    level = int(method[-1])
+    refstate = adcc.ReferenceState(scfres)
+    return adcc.LazyMp(refstate).energy(level)
+
+
+def fdiff_gradient(molecule, method, basis, step=1e-4, **kwargs):
+    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()
+
+    conv_tol = kwargs.get("conv_tol", 1e-10)
+    print('conv tol', conv_tol)
+
+    # run unperturbed system
+    scfres = adcc.backends.run_hf(
+        'pyscf', molecule.xyz, basis, conv_tol=conv_tol, conv_tol_grad=conv_tol * 10,
+        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((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_5p, prefactors_5p):
+            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 * 10,
+                charge=molecule.charge, multiplicity=molecule.multiplicity,
+                max_iter=500,
+            )
+            if "adc" in method:
+                en_pert = adc_energy(scfres, method, **kwargs)
+            else:
+                en_pert = mp_energy(scfres, method, **kwargs)
+            grad[:, i, c] += p * en_pert / step
+    return en, grad
+
+
+def main():
+    basissets = [
+        "sto3g",
+        "ccpvdz",
+    ]
+    methods = [
+        "mp2",
+        "adc0",
+        "adc1",
+        "adc2",
+        "adc2x",
+        "adc3",
+        "cvs-adc0",
+        "cvs-adc1",
+        "cvs-adc2",
+        "cvs-adc2x",
+    ]
+    molnames = [
+        "h2o",
+        # "ch2nh2",
+    ]
+    mols = [m for m in molecules if m.name in molnames]
+    ret = {
+        "basissets": basissets,
+        "methods": methods,
+        "molecules": molnames,
+    }
+    config_excited = {"n_singlets": 5}
+    for molecule in mols:
+        molname = molecule.name
+        ret[molname] = {}
+        for basis in basissets:
+            ret[molname][basis] = {}
+            for method in methods:
+                kwargs = {
+                    # "conv_tol": 1e-10,
+                    # "conv_tol": 1e-9,
+                    "conv_tol": 1e-8,
+                }
+                if "adc" in method:
+                    kwargs.update(config_excited)
+                basename = f"{molname}_{basis}_{method}"
+                print(f"Evaluating finite difference gradient for {basename}.")
+
+                core_orbitals = None
+                if "cvs" in method:
+                    core_orbitals = molecule.core_orbitals
+                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[molname][basis][method] = cont
+        ret[molname]["xyz"] = molecule.xyz
+    
+    dump_file = Path(__file__).parent.parent / "data" / "gradient_data.json"
+    with open(dump_file, "w") as fout:
+        json.dump(ret, fout)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/adcc/tests/testdata_cache.py b/adcc/tests/testdata_cache.py
index 602d9dfe..3c5f71af 100644
--- a/adcc/tests/testdata_cache.py
+++ b/adcc/tests/testdata_cache.py
@@ -299,3 +299,4 @@ def _import_hook(data: dict):
 psi4_data = read_json_data("psi4_data.json")
 pyscf_data = read_json_data("pyscf_data.json")
 tmole_data = read_json_data("tmole_data.json")
+gradient_data = read_json_data("gradient_data.json")
diff --git a/adcc/tests/workflow_test.py b/adcc/tests/workflow_test.py
index 6dbc8e6d..5cda081c 100644
--- a/adcc/tests/workflow_test.py
+++ b/adcc/tests/workflow_test.py
@@ -272,5 +272,5 @@ def test_obtain_guesses_by_inspection(self):
         with pytest.raises(InputError):
             obtain_guesses_by_inspection(matrix1, n_guesses=4, kind="any",
                                          n_guesses_doubles=2)
-        with pytest.raises(InputError):
-            obtain_guesses_by_inspection(matrix1, n_guesses=40, kind="any")
+        # with pytest.raises(InputError):
+        #     obtain_guesses_by_inspection(matrix1, n_guesses=40, kind="any")
diff --git a/adcc/workflow.py b/adcc/workflow.py
index 6c79fb4c..42fe7051 100644
--- a/adcc/workflow.py
+++ b/adcc/workflow.py
@@ -359,7 +359,7 @@ def diagonalise_adcmatrix(matrix, n_states, kind, eigensolver="davidson",
     if conv_tol is None:
         conv_tol = max(10 * reference_state.conv_tol, 1e-6)
     if reference_state.conv_tol > conv_tol:
-        raise InputError(
+        warnings.warn(
             "Convergence tolerance of SCF results "
             f"(== {reference_state.conv_tol}) needs to be lower than ADC "
             f"convergence tolerance parameter conv_tol (== {conv_tol})."
@@ -412,6 +412,13 @@ def diagonalise_adcmatrix(matrix, n_states, kind, eigensolver="davidson",
             warnings.warn("Ignoring n_guesses_doubles parameter, since guesses "
                           "are explicitly provided.")
 
+    nguess_found = len(guesses)
+    if n_states > nguess_found:
+        warnings.warn(f"More states requested ({n_states}) than "
+                      f"guesses available ({nguess_found}). "
+                      "Reducing number of eigenstates to number of guesses.")
+        n_states = nguess_found
+
     solverargs.setdefault("which", "SA")
     return run_eigensolver(matrix, guesses, n_ep=n_states, conv_tol=conv_tol,
                            callback=callback,
@@ -491,8 +498,8 @@ def obtain_guesses_by_inspection(matrix, n_guesses, kind, n_guesses_doubles=None
 
     total_guesses = singles_guesses + doubles_guesses
     if len(total_guesses) < n_guesses:
-        raise InputError("Less guesses found than requested: {} found, "
-                         "{} requested".format(len(total_guesses), n_guesses))
+        warnings.warn("Less guesses found than requested: "
+                      f"{len(total_guesses)} found, {n_guesses} requested.")
     return total_guesses