diff --git a/adcc/IsrMatrix.py b/adcc/IsrMatrix.py
new file mode 100644
index 00000000..2ba35ad0
--- /dev/null
+++ b/adcc/IsrMatrix.py
@@ -0,0 +1,176 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2022 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 libadcc
+
+from .AdcMatrix import AdcMatrixlike
+from .LazyMp import LazyMp
+from .adc_pp import bmatrix as ppbmatrix
+from .timings import Timer, timed_member_call
+from .AdcMethod import AdcMethod
+from .OneParticleOperator import OneParticleOperator
+from .AmplitudeVector import AmplitudeVector
+
+
+class IsrMatrix(AdcMatrixlike):
+ # Default perturbation-theory orders for the matrix blocks (== standard ADC-PP).
+ default_block_orders = {
+ # ph_ph=0, ph_pphh=None, pphh_ph=None, pphh_pphh=None),
+ "adc0": dict(ph_ph=0, ph_pphh=None, pphh_ph=None, pphh_pphh=None), # noqa: E501
+ "adc1": dict(ph_ph=1, ph_pphh=None, pphh_ph=None, pphh_pphh=None), # noqa: E501
+ "adc2": dict(ph_ph=2, ph_pphh=1, pphh_ph=1, pphh_pphh=0), # noqa: E501
+ "adc2x": dict(ph_ph=2, ph_pphh=1, pphh_ph=1, pphh_pphh=1), # noqa: E501
+ "adc3": dict(ph_ph=3, ph_pphh=2, pphh_ph=2, pphh_pphh=1), # noqa: E501
+ }
+
+ def __init__(self, method, hf_or_mp, operators, block_orders=None):
+ """
+ Initialise an ISR matrix of a given one-particle operator
+ for the provided ADC method.
+
+ Parameters
+ ----------
+ method : str or adcc.AdcMethod
+ Method to use.
+ hf_or_mp : adcc.ReferenceState or adcc.LazyMp
+ HF reference or MP ground state.
+ operators : adcc.OneParticleOperator or list of adcc.OneParticleOperator
+ objects
+ One-particle matrix elements associated with a one-particle operator.
+ block_orders : optional
+ The order of perturbation theory to employ for each matrix block.
+ If not set, defaults according to the selected ADC method are chosen.
+ """
+ if isinstance(hf_or_mp, (libadcc.ReferenceState,
+ libadcc.HartreeFockSolution_i)):
+ hf_or_mp = LazyMp(hf_or_mp)
+ if not isinstance(hf_or_mp, LazyMp):
+ raise TypeError("hf_or_mp is not a valid object. It needs to be "
+ "either a LazyMp, a ReferenceState or a "
+ "HartreeFockSolution_i.")
+
+ if not isinstance(method, AdcMethod):
+ method = AdcMethod(method)
+
+ if not isinstance(operators, list):
+ self.operators = [operators]
+ else:
+ self.operators = operators.copy()
+ if not all(isinstance(op, OneParticleOperator) for op in self.operators):
+ raise TypeError("operators is not a valid object. It needs to be "
+ "either an OneParticleOperator or a list of "
+ "OneParticleOperator objects.")
+
+ self.timer = Timer()
+ self.method = method
+ self.ground_state = hf_or_mp
+ self.reference_state = hf_or_mp.reference_state
+ self.mospaces = hf_or_mp.reference_state.mospaces
+ self.is_core_valence_separated = method.is_core_valence_separated
+ self.ndim = 2
+ self.extra_terms = []
+
+ # Determine orders of PT in the blocks
+ if block_orders is None:
+ block_orders = self.default_block_orders[method.base_method.name]
+ else:
+ tmp_orders = self.default_block_orders[method.base_method.name].copy()
+ tmp_orders.update(block_orders)
+ block_orders = tmp_orders
+
+ # Sanity checks on block_orders
+ for block in block_orders.keys():
+ if block not in ("ph_ph", "ph_pphh", "pphh_ph", "pphh_pphh"):
+ raise ValueError(f"Invalid block order key: {block}")
+ if block_orders["ph_pphh"] != block_orders["pphh_ph"]:
+ raise ValueError("ph_pphh and pphh_ph should always have "
+ "the same order")
+ if block_orders["ph_pphh"] is not None \
+ and block_orders["pphh_pphh"] is None:
+ raise ValueError("pphh_pphh cannot be None if ph_pphh isn't.")
+ self.block_orders = block_orders
+
+ # Build the blocks
+ with self.timer.record("build"):
+ variant = None
+ if self.is_core_valence_separated:
+ variant = "cvs"
+ blocks = [{
+ block: ppbmatrix.block(self.ground_state, operator,
+ block.split("_"), order=order,
+ variant=variant)
+ for block, order in self.block_orders.items() if order is not None
+ } for operator in self.operators]
+ # TODO Rename to self.block in 0.16.0
+ self.blocks_ph = [{
+ b: bl[b].apply for b in bl
+ } for bl in blocks]
+
+ @timed_member_call()
+ def matvec(self, v):
+ """
+ Compute the matrix-vector product of the ISR one-particle
+ operator and return the result.
+
+ If a list of OneParticleOperator objects was passed to the class
+ instantiation operator, a list of AmplitudeVector objects is returned.
+ """
+ ret = [
+ sum(block(v) for block in bl_ph.values())
+ for bl_ph in self.blocks_ph
+ ]
+ if len(ret) == 1:
+ return ret[0]
+ else:
+ return ret
+
+ def rmatvec(self, v):
+ # Hermitian operators
+ if all(op.is_symmetric for op in self.operators):
+ return self.matvec(v)
+ else:
+ diffv = [op.ov + op.vo.transpose((1, 0)) for op in self.operators]
+ # anti-Hermitian operators
+ if all(dv.dot(dv) < 1e-12 for dv in diffv):
+ return [
+ AmplitudeVector(ph=-1.0 * mv.ph, pphh=-1.0 * mv.pphh)
+ for mv in self.matvec(v)
+ ]
+ # operators without any symmetry
+ else:
+ return NotImplemented
+
+ def __matmul__(self, other):
+ """
+ If the matrix-vector product is to be calculated with multiple vectors,
+ a list of AmplitudeVector objects is returned.
+
+ Consequently, applying the ISR matrix with multiple operators to multiple
+ vectors results in an N_vectors x N_operators 2D list of AmplitudeVector
+ objects.
+ """
+ if isinstance(other, AmplitudeVector):
+ return self.matvec(other)
+ if isinstance(other, list):
+ if all(isinstance(elem, AmplitudeVector) for elem in other):
+ return [self.matvec(ov) for ov in other]
+ return NotImplemented
diff --git a/adcc/adc_pp/bmatrix.py b/adcc/adc_pp/bmatrix.py
new file mode 100644
index 00000000..a795f9f8
--- /dev/null
+++ b/adcc/adc_pp/bmatrix.py
@@ -0,0 +1,216 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2022 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 collections import namedtuple
+from adcc import block as b
+from adcc.functions import einsum
+from adcc.AmplitudeVector import AmplitudeVector
+
+
+__all__ = ["block"]
+
+
+#
+# Dispatch routine
+#
+
+"""
+`apply` is a function mapping an AmplitudeVector to the contribution of this
+block to the result of applying the ADC matrix.
+"""
+AdcBlock = namedtuple("AdcBlock", ["apply"])
+
+
+def block(ground_state, operator, spaces, order, variant=None):
+ """
+ Gets ground state, one-particle matrix elements associated
+ with a one-particle operator, spaces (ph, pphh and so on)
+ and the perturbation theory order for the block,
+ variant is "cvs" or sth like that.
+
+ The matrix-vector product was derived up to second order
+ using the original equations from
+ J. Schirmer and A. B. Trofimov, J. Chem. Phys. 120, 11449–11464 (2004).
+ """
+ if isinstance(variant, str):
+ variant = [variant]
+ elif variant is None:
+ variant = []
+
+ if ground_state.has_core_occupied_space and "cvs" not in variant:
+ raise ValueError("Cannot run a general (non-core-valence approximated) "
+ "ADC method on top of a ground state with a "
+ "core-valence separation.")
+ if not ground_state.has_core_occupied_space and "cvs" in variant:
+ raise ValueError("Cannot run a core-valence approximated ADC method on "
+ "top of a ground state without a "
+ "core-valence separation.")
+
+ fn = "_".join(["block"] + variant + spaces + [str(order)])
+
+ if fn not in globals():
+ raise ValueError("Could not dispatch: "
+ f"spaces={spaces} order={order} variant=variant")
+ return globals()[fn](ground_state, operator)
+
+
+#
+# 0th order main
+#
+def block_ph_ph_0(ground_state, op):
+ def apply(ampl):
+ return AmplitudeVector(ph=(
+ + 1.0 * einsum('ic,ac->ia', ampl.ph, op.vv)
+ - 1.0 * einsum('ka,ki->ia', ampl.ph, op.oo)
+ ))
+ return AdcBlock(apply)
+
+
+def block_pphh_pphh_0(ground_state, op):
+ def apply(ampl):
+ return AmplitudeVector(pphh=0.5 * (
+ (
+ + 2.0 * einsum('ijcb,ac->ijab', ampl.pphh, op.vv)
+ - 2.0 * einsum('ijca,bc->ijab', ampl.pphh, op.vv)
+ ).antisymmetrise(2, 3)
+ + (
+ - 2.0 * einsum('kjab,ki->ijab', ampl.pphh, op.oo)
+ + 2.0 * einsum('kiab,kj->ijab', ampl.pphh, op.oo)
+ ).antisymmetrise(0, 1)
+ ))
+ return AdcBlock(apply)
+
+
+#
+# 0th order coupling
+#
+def block_ph_pphh_0(ground_state, op):
+ return AdcBlock(lambda ampl: 0, 0)
+
+
+def block_pphh_ph_0(ground_state, op):
+ return AdcBlock(lambda ampl: 0, 0)
+
+
+#
+# 1st order main
+#
+block_ph_ph_1 = block_ph_ph_0
+
+
+#
+# 1st order coupling
+#
+def block_ph_pphh_1(ground_state, op):
+ if op.is_symmetric:
+ op_vo = op.ov.transpose((1, 0))
+ else:
+ op_vo = op.vo.copy()
+ t2 = ground_state.t2(b.oovv)
+
+ def apply(ampl):
+ return AmplitudeVector(ph=0.5 * (
+ - 2.0 * einsum('ilad,ld->ia', ampl.pphh, op.ov)
+ + 2.0 * einsum('ilad,lndf,fn->ia', ampl.pphh, t2, op_vo)
+ + 2.0 * einsum('ilca,lc->ia', ampl.pphh, op.ov)
+ - 2.0 * einsum('ilca,lncf,fn->ia', ampl.pphh, t2, op_vo)
+ - 2.0 * einsum('klad,kled,ei->ia', ampl.pphh, t2, op_vo)
+ - 2.0 * einsum('ilcd,nlcd,an->ia', ampl.pphh, t2, op_vo)
+ ))
+ return AdcBlock(apply)
+
+
+def block_pphh_ph_1(ground_state, op):
+ if op.is_symmetric:
+ op_vo = op.ov.transpose((1, 0))
+ else:
+ op_vo = op.vo.copy()
+ t2 = ground_state.t2(b.oovv)
+
+ def apply(ampl):
+ return AmplitudeVector(pphh=0.5 * (
+ (
+ - 1.0 * einsum('ia,bj->ijab', ampl.ph, op_vo)
+ + 1.0 * einsum('ia,jnbf,nf->ijab', ampl.ph, t2, op.ov)
+ + 1.0 * einsum('ja,bi->ijab', ampl.ph, op_vo)
+ - 1.0 * einsum('ja,inbf,nf->ijab', ampl.ph, t2, op.ov)
+ + 1.0 * einsum('ib,aj->ijab', ampl.ph, op_vo)
+ - 1.0 * einsum('ib,jnaf,nf->ijab', ampl.ph, t2, op.ov)
+ - 1.0 * einsum('jb,ai->ijab', ampl.ph, op_vo)
+ + 1.0 * einsum('jb,inaf,nf->ijab', ampl.ph, t2, op.ov)
+ ).antisymmetrise(0, 1).antisymmetrise(2, 3)
+ + (
+ - 1.0 * einsum('ka,ijeb,ke->ijab', ampl.ph, t2, op.ov)
+ + 1.0 * einsum('kb,ijea,ke->ijab', ampl.ph, t2, op.ov)
+ ).antisymmetrise(2, 3)
+ + (
+ - 1.0 * einsum('ic,njab,nc->ijab', ampl.ph, t2, op.ov)
+ + 1.0 * einsum('jc,niab,nc->ijab', ampl.ph, t2, op.ov)
+ ).antisymmetrise(0, 1)
+ ))
+ return AdcBlock(apply)
+
+
+#
+# 2nd order main
+#
+def block_ph_ph_2(ground_state, op):
+ if op.is_symmetric:
+ op_vo = op.ov.transpose((1, 0))
+ else:
+ op_vo = op.vo.copy()
+ p0 = ground_state.mp2_diffdm
+ t2 = ground_state.t2(b.oovv)
+
+ def apply(ampl):
+ return AmplitudeVector(ph=(
+ # 0th order
+ + 1.0 * einsum('ic,ac->ia', ampl.ph, op.vv)
+ - 1.0 * einsum('ka,ki->ia', ampl.ph, op.oo)
+ # 2nd order
+ # (2,1)
+ - 1.0 * einsum('ic,jc,aj->ia', ampl.ph, p0.ov, op_vo)
+ - 1.0 * einsum('ka,kb,bi->ia', ampl.ph, p0.ov, op_vo)
+ - 1.0 * einsum('ic,ja,jc->ia', ampl.ph, p0.ov, op.ov) # h.c.
+ - 1.0 * einsum('ka,ib,kb->ia', ampl.ph, p0.ov, op.ov) # h.c.
+ # (2,2)
+ - 0.25 * einsum('ic,mnef,mnaf,ec->ia', ampl.ph, t2, t2, op.vv)
+ - 0.25 * einsum('ic,mnef,mncf,ae->ia', ampl.ph, t2, t2, op.vv) # h.c.
+ # (2,3)
+ - 0.5 * einsum('ic,mnce,mnaf,ef->ia', ampl.ph, t2, t2, op.vv)
+ + 1.0 * einsum('ic,mncf,jnaf,jm->ia', ampl.ph, t2, t2, op.oo)
+ # (2,4)
+ + 0.25 * einsum('ka,mnef,inef,km->ia', ampl.ph, t2, t2, op.oo)
+ + 0.25 * einsum('ka,mnef,knef,mi->ia', ampl.ph, t2, t2, op.oo) # h.c.
+ # (2,5)
+ - 1.0 * einsum('ka,knef,indf,ed->ia', ampl.ph, t2, t2, op.vv)
+ + 0.5 * einsum('ka,knef,imef,mn->ia', ampl.ph, t2, t2, op.oo)
+ # (2,6)
+ + 0.5 * einsum('kc,knef,inaf,ec->ia', ampl.ph, t2, t2, op.vv)
+ - 0.5 * einsum('kc,mncf,inaf,km->ia', ampl.ph, t2, t2, op.oo)
+ + 0.5 * einsum('kc,inef,kncf,ae->ia', ampl.ph, t2, t2, op.vv) # h.c.
+ - 0.5 * einsum('kc,mnaf,kncf,mi->ia', ampl.ph, t2, t2, op.oo) # h.c.
+ # (2,7)
+ - 1.0 * einsum('kc,kncf,imaf,mn->ia', ampl.ph, t2, t2, op.oo)
+ + 1.0 * einsum('kc,knce,inaf,ef->ia', ampl.ph, t2, t2, op.vv)
+ ))
+ return AdcBlock(apply)
diff --git a/adcc/test_IsrMatrix.py b/adcc/test_IsrMatrix.py
new file mode 100644
index 00000000..eb6c99e9
--- /dev/null
+++ b/adcc/test_IsrMatrix.py
@@ -0,0 +1,118 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2022 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 numpy as np
+
+from adcc.IsrMatrix import IsrMatrix
+from adcc.testdata.cache import cache
+from adcc.OneParticleOperator import product_trace
+from adcc.adc_pp.state2state_transition_dm import state2state_transition_dm
+from adcc.misc import expand_test_templates
+
+
+testcases = [("h2o_sto3g", "singlet"), ("cn_sto3g", "state")]
+methods = ["adc0", "adc1", "adc2"]
+operator_kinds = ["electric", "magnetic"]
+
+
+@expand_test_templates(list(itertools.product(testcases, methods, operator_kinds)))
+class TestIsrMatrix(unittest.TestCase):
+ def template_matrix_vector_product(self, case, method, op_kind):
+ (system, kind) = case
+ state = cache.adc_states[system][method][kind]
+ mp = state.ground_state
+ if op_kind == "electric": # example of a symmetric operator
+ dips = state.reference_state.operators.electric_dipole
+ elif op_kind == "magnetic": # example of an asymmetric operator
+ dips = state.reference_state.operators.magnetic_dipole
+ else:
+ raise NotImplementedError(
+ "Tests are only implemented for"
+ "electric and magnetic dipole operators."
+ )
+
+ # computing Y_m @ B @ Y_n yields the state-to-state
+ # transition dipole moments (n->m) (for n not equal to m)
+ # they can either be obtained using the matvec method of the IsrMatrix
+ # class or via the state-to-state transition density matrices
+ # (the second method serves as a reference here)
+
+ matrix = IsrMatrix(method, mp, dips)
+ for excitation1 in state.excitations:
+ resv = matrix @ excitation1.excitation_vector
+ for excitation2 in state.excitations:
+ s2s_tdm = [excitation2.excitation_vector @ v for v in resv]
+ tdm = state2state_transition_dm(
+ state.property_method, mp,
+ excitation1.excitation_vector,
+ excitation2.excitation_vector,
+ state.matrix.intermediates
+ )
+ s2s_tdm_ref = np.array([product_trace(tdm, dip) for dip in dips])
+ np.testing.assert_allclose(s2s_tdm, s2s_tdm_ref, atol=1e-12)
+
+
+class TestIsrMatrixInterface(unittest.TestCase):
+ def test_matvec(self):
+ system = "h2o_sto3g"
+ method = "adc2"
+ kind = "singlet"
+
+ state = cache.adc_states[system][method][kind]
+ mp = state.ground_state
+ dips = state.reference_state.operators.electric_dipole
+ magdips = state.reference_state.operators.magnetic_dipole
+ vecs = [exc.excitation_vector for exc in state.excitations[:2]]
+
+ matrix_ref = IsrMatrix(method, mp, dips)
+ refv1, refv2 = matrix_ref @ vecs
+
+ # multiple vectors
+ for i, dip in enumerate(dips):
+ matrix = IsrMatrix(method, mp, dip)
+ resv1, resv2 = matrix @ vecs
+ diffs = [refv1[i] - resv1, refv2[i] - resv2]
+ for j in range(2):
+ assert diffs[j].ph.dot(diffs[j].ph) < 1e-12
+ assert diffs[j].pphh.dot(diffs[j].pphh) < 1e-12
+
+ # compute Y_n @ B @ Y_n with matvec and rmatvec
+ # symmetric operators
+ for vec1 in vecs:
+ for vec2 in vecs:
+ resl = [vec1 @ mvprod for mvprod in matrix_ref.matvec(vec2)]
+ resr = [mvprod @ vec2 for mvprod in matrix_ref.rmatvec(vec1)]
+ np.testing.assert_allclose(
+ np.array(resl), np.array(resr), atol=1e-12
+ )
+
+ # anti-symmetric operators
+ magmatrix = IsrMatrix(method, mp, magdips)
+ for vec1 in vecs:
+ for vec2 in vecs:
+ resl = [vec1 @ mvprod for mvprod in magmatrix.matvec(vec2)]
+ resr = [mvprod @ vec2 for mvprod in magmatrix.rmatvec(vec1)]
+ np.testing.assert_allclose(
+ np.array(resl), np.array(resr), atol=1e-12
+ )
diff --git a/examples/water/tpa.py b/examples/water/tpa.py
index 842c2489..be46ae84 100644
--- a/examples/water/tpa.py
+++ b/examples/water/tpa.py
@@ -11,8 +11,7 @@
from adcc.solver import IndexSymmetrisation
from adcc.solver.conjugate_gradient import conjugate_gradient, default_print
from adcc.adc_pp.modified_transition_moments import modified_transition_moments
-from adcc.adc_pp.state2state_transition_dm import state2state_transition_dm
-from adcc.OneParticleOperator import product_trace
+from adcc.IsrMatrix import IsrMatrix
class ShiftedMat(adcc.AdcMatrix):
@@ -44,6 +43,7 @@ def __matmul__(self, other):
# setup modified transition moments
dips = state.reference_state.operators.electric_dipole
rhss = modified_transition_moments("adc2", state.ground_state, dips)
+isrmatrix = IsrMatrix("adc2", state.ground_state, dips)
S = np.zeros((len(state.excitation_energy), 3, 3))
for f, ee in enumerate(state.excitation_energy):
@@ -63,20 +63,13 @@ def __matmul__(self, other):
explicit_symmetrisation=explicit_symmetrisation
)
response.append(res)
+ right_vec = isrmatrix @ state.excitation_vector[f]
for mu in range(3):
for nu in range(mu, 3):
- tdm_mu_f = state2state_transition_dm(
- "adc2", matrix.ground_state, response[mu].solution,
- state.excitation_vector[f]
- )
- tdm_nu_f = state2state_transition_dm(
- "adc2", matrix.ground_state, response[nu].solution,
- state.excitation_vector[f]
- )
# compute the matrix element
S[f, mu, nu] = (
- product_trace(tdm_mu_f, dips[nu])
- + product_trace(tdm_nu_f, dips[mu])
+ response[mu].solution @ right_vec[nu]
+ + response[nu].solution @ right_vec[mu]
)
S[f, nu, mu] = S[f, mu, nu]
print("Two-Photon Matrix for state", f)