From 0cf13fc9fe891b4df110a5e06fa994ae6c9124e6 Mon Sep 17 00:00:00 2001 From: apapapostolou Date: Thu, 11 Aug 2022 17:31:25 +0200 Subject: [PATCH 1/5] implement isrmatrix_vector_product --- adcc/adc_pp/isrmatrix_vector_product.py | 152 +++++++++++++++++++ adcc/adc_pp/test_isrmatrix_vector_product.py | 75 +++++++++ 2 files changed, 227 insertions(+) create mode 100644 adcc/adc_pp/isrmatrix_vector_product.py create mode 100644 adcc/adc_pp/test_isrmatrix_vector_product.py diff --git a/adcc/adc_pp/isrmatrix_vector_product.py b/adcc/adc_pp/isrmatrix_vector_product.py new file mode 100644 index 00000000..4afc37d1 --- /dev/null +++ b/adcc/adc_pp/isrmatrix_vector_product.py @@ -0,0 +1,152 @@ +from adcc import AdcMethod +from adcc import LazyMp +from adcc import block as b +from adcc.functions import einsum +from adcc.AmplitudeVector import AmplitudeVector + + +def isrmvp_adc0(ground_state, dip, vec): + assert isinstance(vec, AmplitudeVector) + ph = ( + + 1.0 * einsum('ic,ac->ia', vec.ph, dip.vv) + - 1.0 * einsum('ka,ki->ia', vec.ph, dip.oo) + ) + return AmplitudeVector(ph=ph) + + +def isrmvp_adc2(ground_state, dip, vec): + assert isinstance(vec, AmplitudeVector) + if dip.is_symmetric: + dip_vo = dip.ov.transpose((1, 0)) + else: + dip_vo = dip.vo.copy() + p0 = ground_state.mp2_diffdm + t2 = ground_state.t2(b.oovv) + + ph = ( + # product of the ph diagonal block with the singles block of the vector + # zeroth order + + 1.0 * einsum('ic,ac->ia', vec.ph, dip.vv) + - 1.0 * einsum('ka,ki->ia', vec.ph, dip.oo) + # second order + # (2,1) + - 1.0 * einsum('ic,jc,aj->ia', vec.ph, p0.ov, dip_vo) + - 1.0 * einsum('ka,kb,bi->ia', vec.ph, p0.ov, dip_vo) + - 1.0 * einsum('ic,ja,jc->ia', vec.ph, p0.ov, dip.ov) # h.c. + - 1.0 * einsum('ka,ib,kb->ia', vec.ph, p0.ov, dip.ov) # h.c. + # (2,2) + - 0.25 * einsum('ic,mnef,mnaf,ec->ia', vec.ph, t2, t2, dip.vv) + - 0.25 * einsum('ic,mnef,mncf,ae->ia', vec.ph, t2, t2, dip.vv) # h.c. + # (2,3) + - 0.5 * einsum('ic,mnce,mnaf,ef->ia', vec.ph, t2, t2, dip.vv) + + 1.0 * einsum('ic,mncf,jnaf,jm->ia', vec.ph, t2, t2, dip.oo) + # (2,4) + + 0.25 * einsum('ka,mnef,inef,km->ia', vec.ph, t2, t2, dip.oo) + + 0.25 * einsum('ka,mnef,knef,mi->ia', vec.ph, t2, t2, dip.oo) # h.c. + # (2,5) + - 1.0 * einsum('ka,knef,indf,ed->ia', vec.ph, t2, t2, dip.vv) + + 0.5 * einsum('ka,knef,imef,mn->ia', vec.ph, t2, t2, dip.oo) + # (2,6) + + 0.5 * einsum('kc,knef,inaf,ec->ia', vec.ph, t2, t2, dip.vv) + - 0.5 * einsum('kc,mncf,inaf,km->ia', vec.ph, t2, t2, dip.oo) + + 0.5 * einsum('kc,inef,kncf,ae->ia', vec.ph, t2, t2, dip.vv) # h.c. + - 0.5 * einsum('kc,mnaf,kncf,mi->ia', vec.ph, t2, t2, dip.oo) # h.c. + # (2,7) + - 1.0 * einsum('kc,kncf,imaf,mn->ia', vec.ph, t2, t2, dip.oo) + + 1.0 * einsum('kc,knce,inaf,ef->ia', vec.ph, t2, t2, dip.vv) + + # product of the ph-2p2h coupling block with the doubles block of the vector + + 0.5 * ( + - 2.0 * einsum('ilad,ld->ia', vec.pphh, dip.ov) + + 2.0 * einsum('ilad,lndf,fn->ia', vec.pphh, t2, dip_vo) + + 2.0 * einsum('ilca,lc->ia', vec.pphh, dip.ov) + - 2.0 * einsum('ilca,lncf,fn->ia', vec.pphh, t2, dip_vo) + - 2.0 * einsum('klad,kled,ei->ia', vec.pphh, t2, dip_vo) + - 2.0 * einsum('ilcd,nlcd,an->ia', vec.pphh, t2, dip_vo) + ) + ) + + pphh = ( + # product of the 2p2h-ph coupling block with the singles block of the vector + + 0.5 * ( + ( + - 1.0 * einsum('ia,bj->ijab', vec.ph, dip_vo) + + 1.0 * einsum('ia,jnbf,nf->ijab', vec.ph, t2, dip.ov) + + 1.0 * einsum('ja,bi->ijab', vec.ph, dip_vo) + - 1.0 * einsum('ja,inbf,nf->ijab', vec.ph, t2, dip.ov) + + 1.0 * einsum('ib,aj->ijab', vec.ph, dip_vo) + - 1.0 * einsum('ib,jnaf,nf->ijab', vec.ph, t2, dip.ov) + - 1.0 * einsum('jb,ai->ijab', vec.ph, dip_vo) + + 1.0 * einsum('jb,inaf,nf->ijab', vec.ph, t2, dip.ov) + ).antisymmetrise(0,1).antisymmetrise(2,3) + +( + - 1.0 * einsum('ka,ijeb,ke->ijab', vec.ph, t2, dip.ov) + + 1.0 * einsum('kb,ijea,ke->ijab', vec.ph, t2, dip.ov) + ).antisymmetrise(2,3) + +( + - 1.0 * einsum('ic,njab,nc->ijab', vec.ph, t2, dip.ov) + + 1.0 * einsum('jc,niab,nc->ijab', vec.ph, t2, dip.ov) + ).antisymmetrise(0,1) + ) + + # product of the 2p2h diagonal block with the doubles block of the vector + + 0.5 * ( + ( + + 2.0 * einsum('ijcb,ac->ijab', vec.pphh, dip.vv) + - 2.0 * einsum('ijca,bc->ijab', vec.pphh, dip.vv) + ).antisymmetrise(2,3) + +( + - 2.0 * einsum('kjab,ki->ijab', vec.pphh, dip.oo) + + 2.0 * einsum('kiab,kj->ijab', vec.pphh, dip.oo) + ).antisymmetrise(0,1) + ) + ) + return AmplitudeVector(ph=ph, pphh=pphh) + + +DISPATCH = { + "adc0": isrmvp_adc0, + "adc1": isrmvp_adc0, # identical to ADC(0) + "adc2": isrmvp_adc2, +} + + +def isrmatrix_vector_product(method, ground_state, dips, vec): + """Compute the matrix-vector product of an ISR one-particle operator + for the provided ADC method. + The product was derived using the original equations from the work of + Schirmer and Trofimov (J. Schirmer and A. B. Trofimov, “Intermediate state + representation approach to physical properties of electronically excited + molecules,” J. Chem. Phys. 120, 11449–11464 (2004).). + + Parameters + ---------- + method: str, AdcMethod + The method to use for the computation of the matrix-vector product + ground_state : adcc.LazyMp + The MP ground state + dips : OneParticleOperator or list of OneParticleOperator + One-particle matrix elements associated with the dipole operator + vec: AmplitudeVector + A vector with singles and doubles block + + Returns + ------- + AmplitudeVector or list of AmplitudeVector + """ + if not isinstance(method, AdcMethod): + method = AdcMethod(method) + if method.name not in DISPATCH: + raise NotImplementedError(f"isrmatrix_vector_product is not implemented for {method.name}.") + if not isinstance(ground_state, LazyMp): + raise TypeError("ground_state should be a LazyMp object.") + unpack = False + if not isinstance(dips, list): + unpack = True + dips = [dips] + + ret = [DISPATCH[method.name](ground_state, dip, vec) for dip in dips] + if unpack: + assert len(ret) == 1 + ret = ret[0] + return ret diff --git a/adcc/adc_pp/test_isrmatrix_vector_product.py b/adcc/adc_pp/test_isrmatrix_vector_product.py new file mode 100644 index 00000000..3015308f --- /dev/null +++ b/adcc/adc_pp/test_isrmatrix_vector_product.py @@ -0,0 +1,75 @@ +import unittest +import numpy as np +from adcc.misc import expand_test_templates +from adcc.testdata.cache import cache +from adcc.adc_pp.isrmatrix_vector_product import isrmatrix_vector_product +from adcc.OneParticleOperator import product_trace +from adcc.adc_pp.state2state_transition_dm import state2state_transition_dm + +methods = ["adc0", "adc1", "adc2"] + +@expand_test_templates(methods) +class TestIsrmatrixVectorProduct(unittest.TestCase): + def base_test(self, system, method, kind, op_kind): + 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 isrmatrix_vector_product or + # via the state-to-state transition density matrices + # (the second method serves as a reference here) + + for excitation1 in state.excitations: + product_vecs = isrmatrix_vector_product( + method, mp, dips, excitation1.excitation_vector + ) + for excitation2 in state.excitations: + s2s_tdm = [excitation2.excitation_vector @ pv for pv in product_vecs] + tdm = state2state_transition_dm( + state.property_method, + state.ground_state, + 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) + + # + # General + # + def template_h2o_sto3g_singlets_elec(self, method): + self.base_test("h2o_sto3g", method, "singlet", "electric") + + def template_h2o_def2tzvp_singlets_elec(self, method): + self.base_test("h2o_def2tzvp", method, "singlet", "electric") + + def template_h2o_sto3g_triplets_elec(self, method): + self.base_test("h2o_sto3g", method, "triplet", "electric") + + def template_h2o_def2tzvp_triplets_elec(self, method): + self.base_test("h2o_def2tzvp", method, "triplet", "electric") + + def template_cn_sto3g_elec(self, method): + self.base_test("cn_sto3g", method, "state", "electric") + + def template_cn_ccpvdz_elec(self, method): + self.base_test("cn_ccpvdz", method, "state", "electric") + + def template_h2o_sto3g_singlets_mag(self, method): + self.base_test("h2o_sto3g", method, "singlet", "magnetic") + + def template_h2o_sto3g_triplets_mag(self, method): + self.base_test("h2o_sto3g", method, "triplet", "magnetic") + + def template_cn_sto3g_mag(self, method): + self.base_test("cn_sto3g", method, "state", "magnetic") From 996f244f17649ac3cdb99c42b0af09141bc7c76c Mon Sep 17 00:00:00 2001 From: apapapostolou Date: Fri, 30 Sep 2022 17:45:05 +0200 Subject: [PATCH 2/5] restructure isrmatrix vector product --- adcc/IsrMatrix.py | 158 ++++++++++++++ adcc/adc_pp/bmatrix.py | 216 +++++++++++++++++++ adcc/adc_pp/isrmatrix_vector_product.py | 152 ------------- adcc/adc_pp/test_isrmatrix_vector_product.py | 75 ------- adcc/test_IsrMatrix.py | 118 ++++++++++ 5 files changed, 492 insertions(+), 227 deletions(-) create mode 100644 adcc/IsrMatrix.py create mode 100644 adcc/adc_pp/bmatrix.py delete mode 100644 adcc/adc_pp/isrmatrix_vector_product.py delete mode 100644 adcc/adc_pp/test_isrmatrix_vector_product.py create mode 100644 adcc/test_IsrMatrix.py diff --git a/adcc/IsrMatrix.py b/adcc/IsrMatrix.py new file mode 100644 index 00000000..29d79aa2 --- /dev/null +++ b/adcc/IsrMatrix.py @@ -0,0 +1,158 @@ +#!/usr/bin/env python3 +## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab +## --------------------------------------------------------------------- +## +## Copyright (C) 2018 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 .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, dips, 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. + dips : adcc.OneParticleOperator or list of adcc.OneParticleOperator objects + One-particle matrix elements associated with the dipole 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(dips, list): + self.dips = [dips] + else: + self.dips = dips + + 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, dip, block.split("_"), + order=order, variant=variant) + for block, order in self.block_orders.items() if order is not None + } for dip in self.dips] + # 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. + """ + 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(dip.is_symmetric for dip in self.dips): + return self.matvec(v) + else: + diffv = [dip.ov + dip.vo.transpose((1, 0)) for dip in self.dips] + # 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 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..b2ff42a2 --- /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) 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 . +## +## --------------------------------------------------------------------- +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, dip, spaces, order, variant=None): + """ + Gets ground state, one-particle matrix elements associated + with the dipole 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, dip) + + +# +# 0th order main +# +def block_ph_ph_0(ground_state, dip): + def apply(ampl): + return AmplitudeVector(ph=( + + 1.0 * einsum('ic,ac->ia', ampl.ph, dip.vv) + - 1.0 * einsum('ka,ki->ia', ampl.ph, dip.oo) + )) + return AdcBlock(apply) + + +def block_pphh_pphh_0(ground_state, dip): + def apply(ampl): + return AmplitudeVector(pphh=0.5 * ( + ( + + 2.0 * einsum('ijcb,ac->ijab', ampl.pphh, dip.vv) + - 2.0 * einsum('ijca,bc->ijab', ampl.pphh, dip.vv) + ).antisymmetrise(2, 3) + + ( + - 2.0 * einsum('kjab,ki->ijab', ampl.pphh, dip.oo) + + 2.0 * einsum('kiab,kj->ijab', ampl.pphh, dip.oo) + ).antisymmetrise(0, 1) + )) + return AdcBlock(apply) + + +# +# 0th order coupling +# +def block_ph_pphh_0(ground_state, dip): + return AdcBlock(lambda ampl: 0, 0) + + +def block_pphh_ph_0(ground_state, dip): + 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, dip): + if dip.is_symmetric: + dip_vo = dip.ov.transpose((1, 0)) + else: + dip_vo = dip.vo.copy() + t2 = ground_state.t2(b.oovv) + + def apply(ampl): + return AmplitudeVector(ph=0.5 * ( + - 2.0 * einsum('ilad,ld->ia', ampl.pphh, dip.ov) + + 2.0 * einsum('ilad,lndf,fn->ia', ampl.pphh, t2, dip_vo) + + 2.0 * einsum('ilca,lc->ia', ampl.pphh, dip.ov) + - 2.0 * einsum('ilca,lncf,fn->ia', ampl.pphh, t2, dip_vo) + - 2.0 * einsum('klad,kled,ei->ia', ampl.pphh, t2, dip_vo) + - 2.0 * einsum('ilcd,nlcd,an->ia', ampl.pphh, t2, dip_vo) + )) + return AdcBlock(apply) + + +def block_pphh_ph_1(ground_state, dip): + if dip.is_symmetric: + dip_vo = dip.ov.transpose((1, 0)) + else: + dip_vo = dip.vo.copy() + t2 = ground_state.t2(b.oovv) + + def apply(ampl): + return AmplitudeVector(pphh=0.5 * ( + ( + - 1.0 * einsum('ia,bj->ijab', ampl.ph, dip_vo) + + 1.0 * einsum('ia,jnbf,nf->ijab', ampl.ph, t2, dip.ov) + + 1.0 * einsum('ja,bi->ijab', ampl.ph, dip_vo) + - 1.0 * einsum('ja,inbf,nf->ijab', ampl.ph, t2, dip.ov) + + 1.0 * einsum('ib,aj->ijab', ampl.ph, dip_vo) + - 1.0 * einsum('ib,jnaf,nf->ijab', ampl.ph, t2, dip.ov) + - 1.0 * einsum('jb,ai->ijab', ampl.ph, dip_vo) + + 1.0 * einsum('jb,inaf,nf->ijab', ampl.ph, t2, dip.ov) + ).antisymmetrise(0, 1).antisymmetrise(2, 3) + + ( + - 1.0 * einsum('ka,ijeb,ke->ijab', ampl.ph, t2, dip.ov) + + 1.0 * einsum('kb,ijea,ke->ijab', ampl.ph, t2, dip.ov) + ).antisymmetrise(2, 3) + + ( + - 1.0 * einsum('ic,njab,nc->ijab', ampl.ph, t2, dip.ov) + + 1.0 * einsum('jc,niab,nc->ijab', ampl.ph, t2, dip.ov) + ).antisymmetrise(0, 1) + )) + return AdcBlock(apply) + + +# +# 2nd order main +# +def block_ph_ph_2(ground_state, dip): + if dip.is_symmetric: + dip_vo = dip.ov.transpose((1, 0)) + else: + dip_vo = dip.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, dip.vv) + - 1.0 * einsum('ka,ki->ia', ampl.ph, dip.oo) + # 2nd order + # (2,1) + - 1.0 * einsum('ic,jc,aj->ia', ampl.ph, p0.ov, dip_vo) + - 1.0 * einsum('ka,kb,bi->ia', ampl.ph, p0.ov, dip_vo) + - 1.0 * einsum('ic,ja,jc->ia', ampl.ph, p0.ov, dip.ov) # h.c. + - 1.0 * einsum('ka,ib,kb->ia', ampl.ph, p0.ov, dip.ov) # h.c. + # (2,2) + - 0.25 * einsum('ic,mnef,mnaf,ec->ia', ampl.ph, t2, t2, dip.vv) + - 0.25 * einsum('ic,mnef,mncf,ae->ia', ampl.ph, t2, t2, dip.vv) # h.c. + # (2,3) + - 0.5 * einsum('ic,mnce,mnaf,ef->ia', ampl.ph, t2, t2, dip.vv) + + 1.0 * einsum('ic,mncf,jnaf,jm->ia', ampl.ph, t2, t2, dip.oo) + # (2,4) + + 0.25 * einsum('ka,mnef,inef,km->ia', ampl.ph, t2, t2, dip.oo) + + 0.25 * einsum('ka,mnef,knef,mi->ia', ampl.ph, t2, t2, dip.oo) # h.c. + # (2,5) + - 1.0 * einsum('ka,knef,indf,ed->ia', ampl.ph, t2, t2, dip.vv) + + 0.5 * einsum('ka,knef,imef,mn->ia', ampl.ph, t2, t2, dip.oo) + # (2,6) + + 0.5 * einsum('kc,knef,inaf,ec->ia', ampl.ph, t2, t2, dip.vv) + - 0.5 * einsum('kc,mncf,inaf,km->ia', ampl.ph, t2, t2, dip.oo) + + 0.5 * einsum('kc,inef,kncf,ae->ia', ampl.ph, t2, t2, dip.vv) # h.c. + - 0.5 * einsum('kc,mnaf,kncf,mi->ia', ampl.ph, t2, t2, dip.oo) # h.c. + # (2,7) + - 1.0 * einsum('kc,kncf,imaf,mn->ia', ampl.ph, t2, t2, dip.oo) + + 1.0 * einsum('kc,knce,inaf,ef->ia', ampl.ph, t2, t2, dip.vv) + )) + return AdcBlock(apply) diff --git a/adcc/adc_pp/isrmatrix_vector_product.py b/adcc/adc_pp/isrmatrix_vector_product.py deleted file mode 100644 index 4afc37d1..00000000 --- a/adcc/adc_pp/isrmatrix_vector_product.py +++ /dev/null @@ -1,152 +0,0 @@ -from adcc import AdcMethod -from adcc import LazyMp -from adcc import block as b -from adcc.functions import einsum -from adcc.AmplitudeVector import AmplitudeVector - - -def isrmvp_adc0(ground_state, dip, vec): - assert isinstance(vec, AmplitudeVector) - ph = ( - + 1.0 * einsum('ic,ac->ia', vec.ph, dip.vv) - - 1.0 * einsum('ka,ki->ia', vec.ph, dip.oo) - ) - return AmplitudeVector(ph=ph) - - -def isrmvp_adc2(ground_state, dip, vec): - assert isinstance(vec, AmplitudeVector) - if dip.is_symmetric: - dip_vo = dip.ov.transpose((1, 0)) - else: - dip_vo = dip.vo.copy() - p0 = ground_state.mp2_diffdm - t2 = ground_state.t2(b.oovv) - - ph = ( - # product of the ph diagonal block with the singles block of the vector - # zeroth order - + 1.0 * einsum('ic,ac->ia', vec.ph, dip.vv) - - 1.0 * einsum('ka,ki->ia', vec.ph, dip.oo) - # second order - # (2,1) - - 1.0 * einsum('ic,jc,aj->ia', vec.ph, p0.ov, dip_vo) - - 1.0 * einsum('ka,kb,bi->ia', vec.ph, p0.ov, dip_vo) - - 1.0 * einsum('ic,ja,jc->ia', vec.ph, p0.ov, dip.ov) # h.c. - - 1.0 * einsum('ka,ib,kb->ia', vec.ph, p0.ov, dip.ov) # h.c. - # (2,2) - - 0.25 * einsum('ic,mnef,mnaf,ec->ia', vec.ph, t2, t2, dip.vv) - - 0.25 * einsum('ic,mnef,mncf,ae->ia', vec.ph, t2, t2, dip.vv) # h.c. - # (2,3) - - 0.5 * einsum('ic,mnce,mnaf,ef->ia', vec.ph, t2, t2, dip.vv) - + 1.0 * einsum('ic,mncf,jnaf,jm->ia', vec.ph, t2, t2, dip.oo) - # (2,4) - + 0.25 * einsum('ka,mnef,inef,km->ia', vec.ph, t2, t2, dip.oo) - + 0.25 * einsum('ka,mnef,knef,mi->ia', vec.ph, t2, t2, dip.oo) # h.c. - # (2,5) - - 1.0 * einsum('ka,knef,indf,ed->ia', vec.ph, t2, t2, dip.vv) - + 0.5 * einsum('ka,knef,imef,mn->ia', vec.ph, t2, t2, dip.oo) - # (2,6) - + 0.5 * einsum('kc,knef,inaf,ec->ia', vec.ph, t2, t2, dip.vv) - - 0.5 * einsum('kc,mncf,inaf,km->ia', vec.ph, t2, t2, dip.oo) - + 0.5 * einsum('kc,inef,kncf,ae->ia', vec.ph, t2, t2, dip.vv) # h.c. - - 0.5 * einsum('kc,mnaf,kncf,mi->ia', vec.ph, t2, t2, dip.oo) # h.c. - # (2,7) - - 1.0 * einsum('kc,kncf,imaf,mn->ia', vec.ph, t2, t2, dip.oo) - + 1.0 * einsum('kc,knce,inaf,ef->ia', vec.ph, t2, t2, dip.vv) - - # product of the ph-2p2h coupling block with the doubles block of the vector - + 0.5 * ( - - 2.0 * einsum('ilad,ld->ia', vec.pphh, dip.ov) - + 2.0 * einsum('ilad,lndf,fn->ia', vec.pphh, t2, dip_vo) - + 2.0 * einsum('ilca,lc->ia', vec.pphh, dip.ov) - - 2.0 * einsum('ilca,lncf,fn->ia', vec.pphh, t2, dip_vo) - - 2.0 * einsum('klad,kled,ei->ia', vec.pphh, t2, dip_vo) - - 2.0 * einsum('ilcd,nlcd,an->ia', vec.pphh, t2, dip_vo) - ) - ) - - pphh = ( - # product of the 2p2h-ph coupling block with the singles block of the vector - + 0.5 * ( - ( - - 1.0 * einsum('ia,bj->ijab', vec.ph, dip_vo) - + 1.0 * einsum('ia,jnbf,nf->ijab', vec.ph, t2, dip.ov) - + 1.0 * einsum('ja,bi->ijab', vec.ph, dip_vo) - - 1.0 * einsum('ja,inbf,nf->ijab', vec.ph, t2, dip.ov) - + 1.0 * einsum('ib,aj->ijab', vec.ph, dip_vo) - - 1.0 * einsum('ib,jnaf,nf->ijab', vec.ph, t2, dip.ov) - - 1.0 * einsum('jb,ai->ijab', vec.ph, dip_vo) - + 1.0 * einsum('jb,inaf,nf->ijab', vec.ph, t2, dip.ov) - ).antisymmetrise(0,1).antisymmetrise(2,3) - +( - - 1.0 * einsum('ka,ijeb,ke->ijab', vec.ph, t2, dip.ov) - + 1.0 * einsum('kb,ijea,ke->ijab', vec.ph, t2, dip.ov) - ).antisymmetrise(2,3) - +( - - 1.0 * einsum('ic,njab,nc->ijab', vec.ph, t2, dip.ov) - + 1.0 * einsum('jc,niab,nc->ijab', vec.ph, t2, dip.ov) - ).antisymmetrise(0,1) - ) - - # product of the 2p2h diagonal block with the doubles block of the vector - + 0.5 * ( - ( - + 2.0 * einsum('ijcb,ac->ijab', vec.pphh, dip.vv) - - 2.0 * einsum('ijca,bc->ijab', vec.pphh, dip.vv) - ).antisymmetrise(2,3) - +( - - 2.0 * einsum('kjab,ki->ijab', vec.pphh, dip.oo) - + 2.0 * einsum('kiab,kj->ijab', vec.pphh, dip.oo) - ).antisymmetrise(0,1) - ) - ) - return AmplitudeVector(ph=ph, pphh=pphh) - - -DISPATCH = { - "adc0": isrmvp_adc0, - "adc1": isrmvp_adc0, # identical to ADC(0) - "adc2": isrmvp_adc2, -} - - -def isrmatrix_vector_product(method, ground_state, dips, vec): - """Compute the matrix-vector product of an ISR one-particle operator - for the provided ADC method. - The product was derived using the original equations from the work of - Schirmer and Trofimov (J. Schirmer and A. B. Trofimov, “Intermediate state - representation approach to physical properties of electronically excited - molecules,” J. Chem. Phys. 120, 11449–11464 (2004).). - - Parameters - ---------- - method: str, AdcMethod - The method to use for the computation of the matrix-vector product - ground_state : adcc.LazyMp - The MP ground state - dips : OneParticleOperator or list of OneParticleOperator - One-particle matrix elements associated with the dipole operator - vec: AmplitudeVector - A vector with singles and doubles block - - Returns - ------- - AmplitudeVector or list of AmplitudeVector - """ - if not isinstance(method, AdcMethod): - method = AdcMethod(method) - if method.name not in DISPATCH: - raise NotImplementedError(f"isrmatrix_vector_product is not implemented for {method.name}.") - if not isinstance(ground_state, LazyMp): - raise TypeError("ground_state should be a LazyMp object.") - unpack = False - if not isinstance(dips, list): - unpack = True - dips = [dips] - - ret = [DISPATCH[method.name](ground_state, dip, vec) for dip in dips] - if unpack: - assert len(ret) == 1 - ret = ret[0] - return ret diff --git a/adcc/adc_pp/test_isrmatrix_vector_product.py b/adcc/adc_pp/test_isrmatrix_vector_product.py deleted file mode 100644 index 3015308f..00000000 --- a/adcc/adc_pp/test_isrmatrix_vector_product.py +++ /dev/null @@ -1,75 +0,0 @@ -import unittest -import numpy as np -from adcc.misc import expand_test_templates -from adcc.testdata.cache import cache -from adcc.adc_pp.isrmatrix_vector_product import isrmatrix_vector_product -from adcc.OneParticleOperator import product_trace -from adcc.adc_pp.state2state_transition_dm import state2state_transition_dm - -methods = ["adc0", "adc1", "adc2"] - -@expand_test_templates(methods) -class TestIsrmatrixVectorProduct(unittest.TestCase): - def base_test(self, system, method, kind, op_kind): - 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 isrmatrix_vector_product or - # via the state-to-state transition density matrices - # (the second method serves as a reference here) - - for excitation1 in state.excitations: - product_vecs = isrmatrix_vector_product( - method, mp, dips, excitation1.excitation_vector - ) - for excitation2 in state.excitations: - s2s_tdm = [excitation2.excitation_vector @ pv for pv in product_vecs] - tdm = state2state_transition_dm( - state.property_method, - state.ground_state, - 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) - - # - # General - # - def template_h2o_sto3g_singlets_elec(self, method): - self.base_test("h2o_sto3g", method, "singlet", "electric") - - def template_h2o_def2tzvp_singlets_elec(self, method): - self.base_test("h2o_def2tzvp", method, "singlet", "electric") - - def template_h2o_sto3g_triplets_elec(self, method): - self.base_test("h2o_sto3g", method, "triplet", "electric") - - def template_h2o_def2tzvp_triplets_elec(self, method): - self.base_test("h2o_def2tzvp", method, "triplet", "electric") - - def template_cn_sto3g_elec(self, method): - self.base_test("cn_sto3g", method, "state", "electric") - - def template_cn_ccpvdz_elec(self, method): - self.base_test("cn_ccpvdz", method, "state", "electric") - - def template_h2o_sto3g_singlets_mag(self, method): - self.base_test("h2o_sto3g", method, "singlet", "magnetic") - - def template_h2o_sto3g_triplets_mag(self, method): - self.base_test("h2o_sto3g", method, "triplet", "magnetic") - - def template_cn_sto3g_mag(self, method): - self.base_test("cn_sto3g", method, "state", "magnetic") diff --git a/adcc/test_IsrMatrix.py b/adcc/test_IsrMatrix.py new file mode 100644 index 00000000..9037a6bf --- /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) 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 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 + ) From 64fd2943f040de6182ff3004c94d05f26b631f2e Mon Sep 17 00:00:00 2001 From: apapapostolou Date: Tue, 4 Oct 2022 13:05:24 +0200 Subject: [PATCH 3/5] make suggested changes --- adcc/IsrMatrix.py | 42 ++++++++++---- adcc/adc_pp/bmatrix.py | 128 ++++++++++++++++++++--------------------- adcc/test_IsrMatrix.py | 2 +- 3 files changed, 95 insertions(+), 77 deletions(-) diff --git a/adcc/IsrMatrix.py b/adcc/IsrMatrix.py index 29d79aa2..95f4e70d 100644 --- a/adcc/IsrMatrix.py +++ b/adcc/IsrMatrix.py @@ -2,7 +2,7 @@ ## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab ## --------------------------------------------------------------------- ## -## Copyright (C) 2018 by the adcc authors +## Copyright (C) 2022 by the adcc authors ## ## This file is part of adcc. ## @@ -27,6 +27,7 @@ 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 @@ -41,7 +42,7 @@ class IsrMatrix(AdcMatrixlike): "adc3": dict(ph_ph=3, ph_pphh=2, pphh_ph=2, pphh_pphh=1), # noqa: E501 } - def __init__(self, method, hf_or_mp, dips, block_orders=None): + 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. @@ -52,8 +53,9 @@ def __init__(self, method, hf_or_mp, dips, block_orders=None): Method to use. hf_or_mp : adcc.ReferenceState or adcc.LazyMp HF reference or MP ground state. - dips : adcc.OneParticleOperator or list of adcc.OneParticleOperator objects - One-particle matrix elements associated with the dipole operator. + 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. @@ -69,10 +71,14 @@ def __init__(self, method, hf_or_mp, dips, block_orders=None): if not isinstance(method, AdcMethod): method = AdcMethod(method) - if not isinstance(dips, list): - self.dips = [dips] + if not isinstance(operators, list): + self.operators = [operators] else: - self.dips = dips + self.operators = operators + 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 @@ -109,10 +115,11 @@ def __init__(self, method, hf_or_mp, dips, block_orders=None): if self.is_core_valence_separated: variant = "cvs" blocks = [{ - block: ppbmatrix.block(self.ground_state, dip, block.split("_"), - order=order, variant=variant) + 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 dip in self.dips] + } 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 @@ -123,6 +130,9 @@ 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()) @@ -135,10 +145,10 @@ def matvec(self, v): def rmatvec(self, v): # Hermitian operators - if all(dip.is_symmetric for dip in self.dips): + if all(op.is_symmetric for op in self.operators): return self.matvec(v) else: - diffv = [dip.ov + dip.vo.transpose((1, 0)) for dip in self.dips] + 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 [ @@ -150,6 +160,14 @@ def rmatvec(self, v): 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): diff --git a/adcc/adc_pp/bmatrix.py b/adcc/adc_pp/bmatrix.py index b2ff42a2..a795f9f8 100644 --- a/adcc/adc_pp/bmatrix.py +++ b/adcc/adc_pp/bmatrix.py @@ -2,7 +2,7 @@ ## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab ## --------------------------------------------------------------------- ## -## Copyright (C) 2020 by the adcc authors +## Copyright (C) 2022 by the adcc authors ## ## This file is part of adcc. ## @@ -40,10 +40,10 @@ AdcBlock = namedtuple("AdcBlock", ["apply"]) -def block(ground_state, dip, spaces, order, variant=None): +def block(ground_state, operator, spaces, order, variant=None): """ Gets ground state, one-particle matrix elements associated - with the dipole operator, spaces (ph, pphh and so on) + 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. @@ -70,31 +70,31 @@ def block(ground_state, dip, spaces, order, variant=None): if fn not in globals(): raise ValueError("Could not dispatch: " f"spaces={spaces} order={order} variant=variant") - return globals()[fn](ground_state, dip) + return globals()[fn](ground_state, operator) # # 0th order main # -def block_ph_ph_0(ground_state, dip): +def block_ph_ph_0(ground_state, op): def apply(ampl): return AmplitudeVector(ph=( - + 1.0 * einsum('ic,ac->ia', ampl.ph, dip.vv) - - 1.0 * einsum('ka,ki->ia', ampl.ph, dip.oo) + + 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, dip): +def block_pphh_pphh_0(ground_state, op): def apply(ampl): return AmplitudeVector(pphh=0.5 * ( ( - + 2.0 * einsum('ijcb,ac->ijab', ampl.pphh, dip.vv) - - 2.0 * einsum('ijca,bc->ijab', ampl.pphh, dip.vv) + + 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, dip.oo) - + 2.0 * einsum('kiab,kj->ijab', ampl.pphh, dip.oo) + - 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) @@ -103,11 +103,11 @@ def apply(ampl): # # 0th order coupling # -def block_ph_pphh_0(ground_state, dip): +def block_ph_pphh_0(ground_state, op): return AdcBlock(lambda ampl: 0, 0) -def block_pphh_ph_0(ground_state, dip): +def block_pphh_ph_0(ground_state, op): return AdcBlock(lambda ampl: 0, 0) @@ -120,51 +120,51 @@ def block_pphh_ph_0(ground_state, dip): # # 1st order coupling # -def block_ph_pphh_1(ground_state, dip): - if dip.is_symmetric: - dip_vo = dip.ov.transpose((1, 0)) +def block_ph_pphh_1(ground_state, op): + if op.is_symmetric: + op_vo = op.ov.transpose((1, 0)) else: - dip_vo = dip.vo.copy() + 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, dip.ov) - + 2.0 * einsum('ilad,lndf,fn->ia', ampl.pphh, t2, dip_vo) - + 2.0 * einsum('ilca,lc->ia', ampl.pphh, dip.ov) - - 2.0 * einsum('ilca,lncf,fn->ia', ampl.pphh, t2, dip_vo) - - 2.0 * einsum('klad,kled,ei->ia', ampl.pphh, t2, dip_vo) - - 2.0 * einsum('ilcd,nlcd,an->ia', ampl.pphh, t2, dip_vo) + - 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, dip): - if dip.is_symmetric: - dip_vo = dip.ov.transpose((1, 0)) +def block_pphh_ph_1(ground_state, op): + if op.is_symmetric: + op_vo = op.ov.transpose((1, 0)) else: - dip_vo = dip.vo.copy() + 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, dip_vo) - + 1.0 * einsum('ia,jnbf,nf->ijab', ampl.ph, t2, dip.ov) - + 1.0 * einsum('ja,bi->ijab', ampl.ph, dip_vo) - - 1.0 * einsum('ja,inbf,nf->ijab', ampl.ph, t2, dip.ov) - + 1.0 * einsum('ib,aj->ijab', ampl.ph, dip_vo) - - 1.0 * einsum('ib,jnaf,nf->ijab', ampl.ph, t2, dip.ov) - - 1.0 * einsum('jb,ai->ijab', ampl.ph, dip_vo) - + 1.0 * einsum('jb,inaf,nf->ijab', ampl.ph, t2, dip.ov) + - 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, dip.ov) - + 1.0 * einsum('kb,ijea,ke->ijab', ampl.ph, t2, dip.ov) + - 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, dip.ov) - + 1.0 * einsum('jc,niab,nc->ijab', ampl.ph, t2, dip.ov) + - 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) @@ -173,44 +173,44 @@ def apply(ampl): # # 2nd order main # -def block_ph_ph_2(ground_state, dip): - if dip.is_symmetric: - dip_vo = dip.ov.transpose((1, 0)) +def block_ph_ph_2(ground_state, op): + if op.is_symmetric: + op_vo = op.ov.transpose((1, 0)) else: - dip_vo = dip.vo.copy() + 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, dip.vv) - - 1.0 * einsum('ka,ki->ia', ampl.ph, dip.oo) + + 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, dip_vo) - - 1.0 * einsum('ka,kb,bi->ia', ampl.ph, p0.ov, dip_vo) - - 1.0 * einsum('ic,ja,jc->ia', ampl.ph, p0.ov, dip.ov) # h.c. - - 1.0 * einsum('ka,ib,kb->ia', ampl.ph, p0.ov, dip.ov) # h.c. + - 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, dip.vv) - - 0.25 * einsum('ic,mnef,mncf,ae->ia', ampl.ph, t2, t2, dip.vv) # h.c. + - 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, dip.vv) - + 1.0 * einsum('ic,mncf,jnaf,jm->ia', ampl.ph, t2, t2, dip.oo) + - 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, dip.oo) - + 0.25 * einsum('ka,mnef,knef,mi->ia', ampl.ph, t2, t2, dip.oo) # h.c. + + 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, dip.vv) - + 0.5 * einsum('ka,knef,imef,mn->ia', ampl.ph, t2, t2, dip.oo) + - 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, dip.vv) - - 0.5 * einsum('kc,mncf,inaf,km->ia', ampl.ph, t2, t2, dip.oo) - + 0.5 * einsum('kc,inef,kncf,ae->ia', ampl.ph, t2, t2, dip.vv) # h.c. - - 0.5 * einsum('kc,mnaf,kncf,mi->ia', ampl.ph, t2, t2, dip.oo) # h.c. + + 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, dip.oo) - + 1.0 * einsum('kc,knce,inaf,ef->ia', ampl.ph, t2, t2, dip.vv) + - 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 index 9037a6bf..eb6c99e9 100644 --- a/adcc/test_IsrMatrix.py +++ b/adcc/test_IsrMatrix.py @@ -2,7 +2,7 @@ ## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab ## --------------------------------------------------------------------- ## -## Copyright (C) 2020 by the adcc authors +## Copyright (C) 2022 by the adcc authors ## ## This file is part of adcc. ## From c275e5a9e3cf1487541162b75e5d320860f4ccd2 Mon Sep 17 00:00:00 2001 From: apapapostolou Date: Thu, 6 Oct 2022 11:12:58 +0200 Subject: [PATCH 4/5] update tpa example --- examples/water/tpa.py | 17 +++++------------ 1 file changed, 5 insertions(+), 12 deletions(-) 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) From 306788d5a2638d85e70f6369b8580c9fc4db4f38 Mon Sep 17 00:00:00 2001 From: apapapostolou Date: Wed, 12 Oct 2022 13:10:06 +0200 Subject: [PATCH 5/5] fix list reference --- adcc/IsrMatrix.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/adcc/IsrMatrix.py b/adcc/IsrMatrix.py index 95f4e70d..2ba35ad0 100644 --- a/adcc/IsrMatrix.py +++ b/adcc/IsrMatrix.py @@ -74,7 +74,7 @@ def __init__(self, method, hf_or_mp, operators, block_orders=None): if not isinstance(operators, list): self.operators = [operators] else: - self.operators = operators + 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 "