diff --git a/adcc/Excitation.py b/adcc/Excitation.py
index 1419e4f6..0ec515db 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 2dd656fe..c4914a91 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 a3274c61..c34f71cd 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 87b136ae..0c00e722 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 96393b1c..a49fb43e 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
@@ -116,6 +183,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 067e5a82..e2546997 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
@@ -111,7 +172,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.
@@ -138,6 +199,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..bffdaa28
--- /dev/null
+++ b/adcc/gradients/TwoParticleDensityMatrix.py
@@ -0,0 +1,289 @@
+#!/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 .
+##
+## ---------------------------------------------------------------------
+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.oooo, b.cvcv,
+ b.ocov, b.cccv, b.cocv, b.ocoo,
+ b.ccco, b.occv, b.ccvv, b.ocvv,
+ b.vvvv,
+ ]
+ 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 "
+ "OneParticleOperator")
+ 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..875684a8
--- /dev/null
+++ b/adcc/gradients/__init__.py
@@ -0,0 +1,218 @@
+#!/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 .
+##
+## ---------------------------------------------------------------------
+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("ik,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 = 0.25 * (
+ einsum("ik,jl->ijkl", delta_ij, g1o.oo + delta_ij)
+ )
+ 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..9e0447b5
--- /dev/null
+++ b/adcc/gradients/amplitude_response.py
@@ -0,0 +1,710 @@
+#!/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 .
+##
+## ---------------------------------------------------------------------
+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 ?
+ # 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 = 1.0 * einsum("iJcd,iJab->abcd", u.pphh, u.pphh)
+
+ # TODO: remove
+ # g2a.ococ *= 0.0
+ # g2a.vvvv *= 0.0
+
+ g1a.co = (
+ - 1.0 * einsum('JbKc,ibKc->Ji', g2a.cvcv, hf.ovcv)
+ - 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..0b1b2ab0
--- /dev/null
+++ b/adcc/gradients/orbital_response.py
@@ -0,0 +1,349 @@
+#!/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 .
+##
+## ---------------------------------------------------------------------
+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:
+ # CVS-ADC0, CVS-ADC1, CVS-ADC2
+ 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)
+ # TODO: generalize once other solvent methods are available
+ if self.hf.environment == "pe":
+ # PE contribution to the orbital Hessian
+ ops = self.hf.operators
+ dm = OneParticleOperator(self.hf, is_symmetric=True)
+ dm.ov = lam.ov
+ ret += ops.pe_induction_elec(dm).ov
+ return evaluate(ret)
+
+
+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/gradients/test_functionality_gradients.py b/adcc/gradients/test_functionality_gradients.py
new file mode 100644
index 00000000..827ee57c
--- /dev/null
+++ b/adcc/gradients/test_functionality_gradients.py
@@ -0,0 +1,98 @@
+#!/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 .
+##
+## ---------------------------------------------------------------------
+import unittest
+import itertools
+
+import adcc
+import adcc.backends
+
+from numpy.testing import assert_allclose
+
+import pytest
+
+from ..misc import expand_test_templates
+from adcc.backends.testing import cached_backend_hf
+from adcc.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"]
+
+combinations = list(itertools.product(molecules, basissets, methods, backends))
+
+
+@pytest.mark.skipif(len(backends) == 0, reason="No backend found.")
+@expand_test_templates(combinations)
+class TestNuclearGradients(unittest.TestCase):
+ def template_nuclear_gradient(self, molecule, basis, method, backend):
+ grad_ref = gradient_data[molecule][basis][method]
+
+ energy_ref = grad_ref["energy"]
+ grad_fdiff = grad_ref["gradient"]
+ kwargs = grad_ref["config"]
+ conv_tol = kwargs["conv_tol"]
+
+ scfres = cached_backend_hf(backend, molecule, basis, conv_tol=1e-13)
+
+ if "adc" in method:
+ # TODO: convergence needs to be very very tight...
+ # so we want to make sure all vectors are tightly converged
+ n_limit = 2
+ kwargs["n_singlets"] = kwargs["n_singlets"] + 6
+ state = adcc.run_adc(scfres, method=method, **kwargs)
+ for ee in state.excitations[:n_limit]:
+ 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=1e-7,
+ 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=1e-8
+ )
+ assert_allclose(
+ grad_fdiff, grad.total, atol=1e-8
+ )
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/test_ExcitedStates.py b/adcc/test_ExcitedStates.py
index b4d1e0b3..cbb1aebf 100644
--- a/adcc/test_ExcitedStates.py
+++ b/adcc/test_ExcitedStates.py
@@ -42,7 +42,8 @@ def base_test(self, system, method, kind):
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/test_workflow.py b/adcc/test_workflow.py
index 8cb0f306..34be4482 100644
--- a/adcc/test_workflow.py
+++ b/adcc/test_workflow.py
@@ -215,7 +215,10 @@ def test_diagonalise_adcmatrix(self):
assert res.converged
assert res.eigenvalues == approx(ref_singlets[:3])
- with pytest.raises(InputError): # Too low tolerance
+ # with pytest.raises(InputError): # Too low tolerance
+ with pytest.warns(
+ UserWarning, match="needs to be lower than ADC convergence tolerance"
+ ):
res = diagonalise_adcmatrix(matrix, n_states=9, kind="singlet",
eigensolver="davidson",
conv_tol=1e-14)
@@ -266,5 +269,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/testdata/cache.py b/adcc/testdata/cache.py
index 6810dca5..3429b820 100644
--- a/adcc/testdata/cache.py
+++ b/adcc/testdata/cache.py
@@ -273,3 +273,4 @@ def read_yaml_data(fname):
tmole_data = read_yaml_data("tmole_dump.yml")
psi4_data = read_yaml_data("psi4_dump.yml")
pyscf_data = read_yaml_data("pyscf_dump.yml")
+gradient_data = read_yaml_data("grad_dump.yml")
diff --git a/adcc/testdata/dump_fdiff_gradient.py b/adcc/testdata/dump_fdiff_gradient.py
new file mode 100644
index 00000000..9c548dac
--- /dev/null
+++ b/adcc/testdata/dump_fdiff_gradient.py
@@ -0,0 +1,185 @@
+#!/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 .
+##
+## ---------------------------------------------------------------------
+from dataclasses import dataclass
+import itertools
+import adcc
+import numpy as np
+import yaml
+from tqdm import tqdm
+
+from pyscf import gto, lib
+
+from static_data 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]
+
+
+molecules = [
+ Molecule("h2o", 0, 1),
+]
+
+
+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) / 100
+
+ # run unperturbed system
+ scfres = adcc.backends.run_hf(
+ 'pyscf', molecule.xyz, basis, conv_tol=conv_tol, conv_tol_grad=conv_tol,
+ charge=molecule.charge, multiplicity=molecule.multiplicity
+ )
+ if "adc" in method:
+ en = adc_energy(scfres, method, **kwargs)
+ ngrad = len(en)
+ else:
+ en = mp_energy(scfres, method, **kwargs)
+ ngrad = 1
+
+ natoms = len(elements)
+ grad = np.zeros((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", # TODO: broken
+ ]
+ molnames = [
+ "h2o",
+ ]
+ mols = [m for m in molecules if m.name in molnames]
+ ret = {
+ "basissets": basissets,
+ "methods": methods,
+ "molecules": molnames,
+ }
+ config_excited = {"n_singlets": 3}
+ 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,
+ }
+ 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
+ with open("grad_dump.yml", "w") as yamlout:
+ yaml.safe_dump(ret, yamlout)
+
+
+if __name__ == "__main__":
+ main()
diff --git a/adcc/testdata/grad_dump.yml b/adcc/testdata/grad_dump.yml
new file mode 100644
index 00000000..d2d20eb5
--- /dev/null
+++ b/adcc/testdata/grad_dump.yml
@@ -0,0 +1,622 @@
+basissets:
+- sto3g
+- ccpvdz
+h2o:
+ ccpvdz:
+ adc0:
+ config:
+ conv_tol: 1.0e-10
+ core_orbitals: null
+ n_singlets: 3
+ energy:
+ - -75.34687593732811
+ - -75.28099897302941
+ - -75.27790714502405
+ gradient:
+ - - - 0.05781772690534126
+ - 3.2741809263825417e-10
+ - 0.041246956439863425
+ - - -0.0010630962424329482
+ - 9.022187441587448e-10
+ - -0.060199278326763306
+ - - -0.05675463087391108
+ - 1.9063008949160576e-09
+ - 0.018952317521325313
+ - - - 0.013116452595568262
+ - 4.0745362639427185e-10
+ - 0.009581870828696992
+ - - 0.05172148926794762
+ - -6.475602276623249e-10
+ - -0.08735981993959285
+ - - -0.0648379422709695
+ - 1.7462298274040222e-10
+ - 0.07777794662979431
+ - - - 0.035870219457137864
+ - 4.2928149923682213e-10
+ - 0.02591435170324985
+ - - 0.007818772384780459
+ - 1.2078089639544487e-09
+ - -0.04966128067462705
+ - - -0.04368899230030365
+ - 2.2264430299401283e-09
+ - 0.0237469249550486
+ adc1:
+ config:
+ conv_tol: 1.0e-10
+ core_orbitals: null
+ n_singlets: 3
+ energy:
+ - -75.68543998437532
+ - -75.61931016176936
+ - -75.59868950258485
+ gradient:
+ - - - 0.09274009053478949
+ - -9.240466170012951e-10
+ - 0.06558011189918034
+ - - -0.0026302200931240804
+ - 1.0186340659856796e-09
+ - -0.09467176448379178
+ - - -0.09010987027431838
+ - -6.83940015733242e-10
+ - 0.02909164885932114
+ - - - 0.10931798897217959
+ - -8.221832104027271e-10
+ - 0.078095546028635
+ - - -0.006660707054834347
+ - 1.0550138540565968e-09
+ - -0.10735300996020669
+ - - -0.10265728181548184
+ - -7.566995918750763e-10
+ - 0.029257460359076504
+ - - - 0.01167026567418361
+ - 1.1568772606551647e-09
+ - 0.008451193309156224
+ - - 0.056844835431547835
+ - 7.421476766467094e-10
+ - -0.09296221967088059
+ - - -0.06851509937405353
+ - -1.1496013030409813e-09
+ - 0.08451102519757114
+ adc2:
+ config:
+ conv_tol: 1.0e-10
+ core_orbitals: null
+ n_singlets: 3
+ energy:
+ - -75.9296754082655
+ - -75.85499790115733
+ - -75.84309170418106
+ gradient:
+ - - - 0.11448501809354639
+ - -1.1496013030409813e-09
+ - 0.08113867529755225
+ - - -0.008900642889784649
+ - -3.637978807091713e-11
+ - -0.10905679840652738
+ - - -0.10558437401778065
+ - -2.837623469531536e-10
+ - 0.02791811802308075
+ - - - 0.1131814870168455
+ - -1.3169483281672e-09
+ - 0.08065583834832069
+ - - -0.0044726152045768686
+ - 4.3655745685100555e-11
+ - -0.11437428731733235
+ - - -0.10870887064811541
+ - -3.128661774098873e-10
+ - 0.03371844370121835
+ - - - 0.022055651032133028
+ - -1.4260876923799515e-09
+ - 0.01590312836196972
+ - - 0.05098412500228733
+ - 7.275957614183426e-11
+ - -0.09580103961343411
+ - - -0.07303977584524546
+ - -5.165929906070232e-10
+ - 0.07989790782448836
+ adc2x:
+ config:
+ conv_tol: 1.0e-10
+ core_orbitals: null
+ n_singlets: 3
+ energy:
+ - -75.94713027654686
+ - -75.87072120070411
+ - -75.8606904664802
+ gradient:
+ - - - 0.12108468737278599
+ - 3.5652192309498787e-10
+ - 0.08575623746583005
+ - - -0.00893019756040303
+ - -9.89530235528946e-10
+ - -0.11596749575255672
+ - - -0.11215448959410423
+ - -1.2369127944111824e-10
+ - 0.030211256875190884
+ - - - 0.12104129993531387
+ - 6.039044819772243e-10
+ - 0.08622117596678436
+ - - -0.005538591824006289
+ - -1.1350493878126144e-09
+ - -0.12121302656305488
+ - - -0.11550270780571736
+ - 5.820766091346741e-11
+ - 0.03499184886459261
+ - - - 0.027651791664538905
+ - 8.512870408594608e-10
+ - 0.019813940722087864
+ - - 0.05061761011893395
+ - 4.43833414465189e-10
+ - -0.10117346070910571
+ - - -0.07826940130325966
+ - 2.3283064365386963e-10
+ - 0.08135951995063806
+ adc3:
+ config:
+ conv_tol: 1.0e-10
+ core_orbitals: null
+ n_singlets: 3
+ energy:
+ - -75.9285783500615
+ - -75.8546333440395
+ - -75.84167671840466
+ gradient:
+ - - - 0.11627274371130625
+ - 1.673470251262188e-09
+ - 0.08230358791479375
+ - - -0.006295582810707856
+ - 3.92901711165905e-10
+ - -0.11453791718668072
+ - - -0.10997716135170776
+ - -9.604264050722122e-10
+ - 0.03223432630329626
+ - - - 0.12384291545458836
+ - 1.5425030142068863e-09
+ - 0.0882510380979511
+ - - -0.006875615137687419
+ - 3.4924596548080444e-10
+ - -0.1223434514249675
+ - - -0.1169673006006633
+ - -9.822542779147625e-10
+ - 0.03409241016197484
+ - - - 0.025808624421188142
+ - -1.6007106751203537e-10
+ - 0.01846984199801227
+ - - 0.05330723668157589
+ - 2.473825588822365e-10
+ - -0.1029805622601998
+ - - -0.0791158616921166
+ - -1.229636836796999e-09
+ - 0.08451071724266512
+ cvs-adc0:
+ config:
+ conv_tol: 1.0e-10
+ core_orbitals: 1
+ n_singlets: 3
+ energy:
+ - -55.29234315857762
+ - -55.223374366273546
+ - -54.667166206322406
+ gradient:
+ - - - -0.01364794644905487
+ - -1.4551915228366852e-11
+ - -0.00933599537529517
+ - - 0.005706506060960237
+ - -1.7462298274040222e-10
+ - 0.006095205535530113
+ - - 0.007941438860143535
+ - 2.1100277081131935e-10
+ - 0.003240785052184947
+ - - - -0.03559545424650423
+ - -1.6007106751203537e-10
+ - -0.02466860027925577
+ - - 0.014588374950108118
+ - -7.275957614183426e-11
+ - 0.016633203136734664
+ - - 0.02100707723002415
+ - 2.9103830456733704e-10
+ - 0.008035392405872699
+ - - - 0.2956047719635535
+ - -8.731149137020111e-11
+ - 0.20220450891793007
+ - - -0.11422410366503755
+ - -1.3096723705530167e-10
+ - -0.14526787615614012
+ - - -0.18138067025574856
+ - 2.473825588822365e-10
+ - -0.05693663733109133
+ cvs-adc1:
+ config:
+ conv_tol: 1.0e-10
+ core_orbitals: 1
+ n_singlets: 3
+ energy:
+ - -55.76235801207001
+ - -55.744231177685215
+ - -55.23885575129546
+ gradient:
+ - - - 0.11375777773355367
+ - -8.149072527885437e-10
+ - 0.07506357444799505
+ - - -0.00623249256022973
+ - -1.7462298274040222e-10
+ - -0.10649733465106692
+ - - -0.10752528410375817
+ - 1.0331859812140465e-09
+ - 0.03143375746731181
+ - - - 0.1367510933050653
+ - -5.966285243630409e-10
+ - 0.10287793751922436
+ - - -0.010133748408406973
+ - 1.964508555829525e-10
+ - -0.13692998980695847
+ - - -0.12661734437278938
+ - 6.912159733474255e-10
+ - 0.03405204997397959
+ - - - -0.03586818584881257
+ - -3.7834979593753815e-10
+ - -0.025053270866919775
+ - - 0.028073501482140273
+ - 3.055902197957039e-10
+ - -0.001954444458533544
+ - - 0.007794684926921036
+ - 8.658389560878277e-10
+ - 0.027007713106286246
+ cvs-adc2:
+ config:
+ conv_tol: 1.0e-10
+ core_orbitals: 1
+ n_singlets: 3
+ energy:
+ - -56.46332440074828
+ - -56.39198974450135
+ - -55.93190966430339
+ gradient:
+ - - - 0.06425007266079774
+ - 3.8562575355172157e-10
+ - 0.045520459178078454
+ - - -0.001686401323240716
+ - -6.548361852765083e-11
+ - -0.06586711463023676
+ - - -0.06256367078458425
+ - -3.055902197957039e-10
+ - 0.020346652701846324
+ - - - 0.05899445136310533
+ - -1.7462298274040222e-10
+ - 0.04236197158752475
+ - - 0.005970590587821789
+ - -2.9103830456733704e-11
+ - -0.07167668161127949
+ - - -0.06496504214737797
+ - -5.529727786779404e-10
+ - 0.029314707753655966
+ - - - 0.26677660722634755
+ - -1.3096723705530167e-10
+ - 0.18210556224221364
+ - - -0.08747911818500143
+ - 6.766640581190586e-10
+ - -0.15278875653166324
+ - - -0.17929748843016569
+ - -1.6007106751203537e-10
+ - -0.02931680792971747
+ mp2:
+ config:
+ conv_tol: 1.0e-10
+ core_orbitals: null
+ energy: -76.22940338787355
+ gradient:
+ - - 0.024395443680987228
+ - 5.020410753786564e-10
+ - 0.0177381719040568
+ - - -0.010762449142930564
+ - -3.2741809263825417e-10
+ - -0.01115037479030434
+ - - -0.01363299589138478
+ - -8.294591680169106e-10
+ - -0.006587800489796791
+ sto3g:
+ adc0:
+ config:
+ conv_tol: 1.0e-10
+ core_orbitals: null
+ n_singlets: 3
+ energy:
+ - -73.96883482604383
+ - -73.91481187944892
+ - -73.80575589478212
+ gradient:
+ - - - 0.21500520483095897
+ - 1.7462298274040222e-10
+ - 0.15122539390722523
+ - - 0.07489849283592775
+ - 1.673470251262188e-10
+ - -0.333204234681034
+ - - -0.2899036974267801
+ - 5.093170329928398e-11
+ - 0.18197884086112026
+ - - - 0.1586926045420114
+ - 3.2014213502407074e-10
+ - 0.11135792706045322
+ - - 0.1261293514544377
+ - 2.1100277081131935e-10
+ - -0.3458573076277389
+ - - -0.28482195568358293
+ - 1.7462298274040222e-10
+ - 0.23449938038538676
+ - - - 0.4883978156649391
+ - 3.128661774098873e-10
+ - 0.3478099736967124
+ - - -0.09521430350287119
+ - 1.2369127944111824e-10
+ - -0.38596352994500194
+ - - -0.39318351200199686
+ - 1.382431946694851e-10
+ - 0.038153556030010805
+ adc1:
+ config:
+ conv_tol: 1.0e-10
+ core_orbitals: null
+ n_singlets: 3
+ energy:
+ - -74.47588319597367
+ - -74.38511885048034
+ - -74.35718229209289
+ gradient:
+ - - - 0.25195657482981915
+ - -1.673470251262188e-10
+ - 0.17643765209504636
+ - - 0.04383796465117484
+ - 2.6921043172478676e-10
+ - -0.32756816265464295
+ - - -0.29579454003396677
+ - -5.456968210637569e-10
+ - 0.151130509802897
+ - - - 0.44204019877361134
+ - -1.6007106751203537e-10
+ - 0.31597312413941836
+ - - -0.08205663986882428
+ - 2.6193447411060333e-10
+ - -0.35633136333490256
+ - - -0.35998355938500026
+ - -5.311449058353901e-10
+ - 0.04035823843150865
+ - - - 0.10246272853692062
+ - -9.458744898438454e-11
+ - 0.07140888737194473
+ - - 0.11964530187105993
+ - 2.837623469531536e-10
+ - -0.2768455888144672
+ - - -0.22210803085181396
+ - -4.802132025361061e-10
+ - 0.20543670067127096
+ adc2:
+ config:
+ conv_tol: 1.0e-10
+ core_orbitals: null
+ n_singlets: 3
+ energy:
+ - -74.52306494758524
+ - -74.4210231428533
+ - -74.39990473420053
+ gradient:
+ - - - 0.27296630301862024
+ - 1.964508555829525e-10
+ - 0.19153599717537872
+ - - 0.04375171427091118
+ - 3.8562575355172157e-10
+ - -0.3499776718817884
+ - - -0.31671801683842205
+ - 1.8917489796876907e-10
+ - 0.15844167445175117
+ - - - 0.4661114452901529
+ - 2.1827872842550278e-10
+ - 0.3327183580331621
+ - - -0.08526563921623165
+ - 4.656612873077393e-10
+ - -0.37705514822300756
+ - - -0.3808458057610551
+ - 1.382431946694851e-10
+ - 0.044336789804219734
+ - - - 0.12454020935547305
+ - 1.5279510989785194e-10
+ - 0.08709251107939053
+ - - 0.11524763766647084
+ - 4.2928149923682213e-10
+ - -0.29412153316661716
+ - - -0.2397878467090777
+ - -2.9103830456733704e-11
+ - 0.20702902177436044
+ adc2x:
+ config:
+ conv_tol: 1.0e-10
+ core_orbitals: null
+ n_singlets: 3
+ energy:
+ - -74.54846626885896
+ - -74.4421830098894
+ - -74.4221092236318
+ gradient:
+ - - - 0.2889574766741134
+ - 1.4551915228366852e-10
+ - 0.20262688743241597
+ - - 0.045065267200698145
+ - 4.3655745685100555e-11
+ - -0.3685835993819637
+ - - -0.3340227430453524
+ - -1.4551915228366852e-11
+ - 0.16595671254617628
+ - - - 0.48461220582248643
+ - 2.1100277081131935e-10
+ - 0.3459748471796047
+ - - -0.08685884004808031
+ - 6.548361852765083e-11
+ - -0.39460422148840735
+ - - -0.39775336504681036
+ - -8.003553375601768e-11
+ - 0.048629374672600534
+ - - - 0.1366148896777304
+ - 2.1827872842550278e-10
+ - 0.09541060922492761
+ - - 0.11530230973585276
+ - 5.093170329928398e-11
+ - -0.30678889842238277
+ - - -0.25191719870781526
+ - -4.3655745685100555e-11
+ - 0.2113782895321492
+ adc3:
+ config:
+ conv_tol: 1.0e-10
+ core_orbitals: null
+ n_singlets: 3
+ energy:
+ - -74.55090439658672
+ - -74.44714801129275
+ - -74.4242627370736
+ gradient:
+ - - - 0.2908625259224209
+ - -4.001776687800884e-10
+ - 0.2039448505747714
+ - - 0.04500434071087511
+ - -1.1641532182693481e-10
+ - -0.3704894050679286
+ - - -0.3358668661385309
+ - 3.346940502524376e-10
+ - 0.16654455448588124
+ - - - 0.48658428593626013
+ - -4.656612873077393e-10
+ - 0.34738898534851614
+ - - -0.08571265437058173
+ - -8.003553375601768e-11
+ - -0.3983368489643908
+ - - -0.40087163097632583
+ - 3.637978807091713e-10
+ - 0.05094786360859871
+ - - - 0.13570507940312382
+ - -3.92901711165905e-10
+ - 0.09476317036023829
+ - - 0.11679784574516816
+ - -1.0186340659856796e-10
+ - -0.3079343543649884
+ - - -0.2525029244279722
+ - 3.41970007866621e-10
+ - 0.2131711840411299
+ cvs-adc0:
+ config:
+ conv_tol: 1.0e-10
+ core_orbitals: 1
+ n_singlets: 3
+ energy:
+ - -54.12308247440466
+ - -53.96000354314296
+ gradient:
+ - - - 0.13670422898576362
+ - 4.220055416226387e-10
+ - 0.09585740162583534
+ - - 0.09903088509599911
+ - -2.546585164964199e-10
+ - -0.284258487423358
+ - - -0.23573511330323527
+ - 3.5652192309498787e-10
+ - 0.18840108531730948
+ - - - 0.4100968398124678
+ - 3.5652192309498787e-10
+ - 0.2924419815244619
+ - - -0.07108191118459217
+ - -1.7462298274040222e-10
+ - -0.3370177824544953
+ - - -0.3390149278129684
+ - 3.128661774098873e-10
+ - 0.04457580037706066
+ cvs-adc1:
+ config:
+ conv_tol: 1.0e-10
+ core_orbitals: 1
+ n_singlets: 3
+ energy:
+ - -54.85354413664419
+ - -54.7941942631979
+ gradient:
+ - - - 0.20441432864026865
+ - -2.0372681319713593e-10
+ - 0.14198894474975532
+ - - 0.03388333945622435
+ - -2.9103830456733704e-10
+ - -0.2622227993415436
+ - - -0.23829766743438086
+ - -2.9831426218152046e-10
+ - 0.12023385422071442
+ - - - 0.31129172734654276
+ - -8.731149137020111e-11
+ - 0.22435429733013734
+ - - -0.03036278401123127
+ - -3.637978807091713e-10
+ - -0.2915527691147872
+ - - -0.28092894267319934
+ - -3.128661774098873e-10
+ - 0.06719847115891753
+ cvs-adc2:
+ config:
+ conv_tol: 1.0e-10
+ core_orbitals: 1
+ n_singlets: 3
+ energy:
+ - -54.989035890846495
+ - -54.905860098139534
+ - -53.1668568155216
+ gradient:
+ - - - 0.23872009358456125
+ - -3.8562575355172157e-10
+ - 0.1671009848068934
+ - - 0.03348872340575326
+ - -2.1100277081131935e-10
+ - -0.29891402418434154
+ - - -0.27220881724497303
+ - 5.093170329928398e-11
+ - 0.13181303915916942
+ - - - 0.361930021950684
+ - -4.220055416226387e-10
+ - 0.25924844172550365
+ - - -0.04305156636110041
+ - -3.2014213502407074e-10
+ - -0.3264197400130797
+ - - -0.3188784557278268
+ - 7.275957614183426e-11
+ - 0.06717129818571266
+ - - - 0.30556880148651544
+ - -4.220055416226387e-10
+ - 0.21378263917722506
+ - - 0.19212377665098757
+ - -3.2014213502407074e-10
+ - -0.5935662729316391
+ - - -0.4976925785158528
+ - 2.1827872842550278e-11
+ - 0.3797836338781053
+ mp2:
+ config:
+ conv_tol: 1.0e-10
+ core_orbitals: null
+ energy: -74.99357808910254
+ gradient:
+ - - 0.096480893320404
+ - -2.764863893389702e-10
+ - 0.06886449074954726
+ - - -0.02618346233794
+ - 1.964508555829525e-10
+ - -0.06597386256908067
+ - - -0.07029743114981102
+ - 8.731149137020111e-11
+ - -0.0028906281804665923
+ xyz: "\n O 0 0 0\n H 0 0 1.795239827225189\n H 1.693194615993441 0 -0.599043184453037\n\
+ \ "
+methods:
+- mp2
+- adc0
+- adc1
+- adc2
+- adc2x
+- adc3
+- cvs-adc0
+- cvs-adc1
+- cvs-adc2
+molecules:
+- h2o
diff --git a/adcc/testdata/static_data.py b/adcc/testdata/static_data.py
index 80720d20..53150dd5 100644
--- a/adcc/testdata/static_data.py
+++ b/adcc/testdata/static_data.py
@@ -22,6 +22,7 @@
## ---------------------------------------------------------------------
import os
+
# all coordinates in Bohr
xyz = {
"h2o": """
diff --git a/adcc/workflow.py b/adcc/workflow.py
index a290032e..3b45bd90 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})."
@@ -406,6 +406,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,
@@ -485,8 +492,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