implement isrmatrix_vector_product

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 <http://www.gnu.org/licenses/>.
+##
+## ---------------------------------------------------------------------
+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

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 <http://www.gnu.org/licenses/>.
+##
+## ---------------------------------------------------------------------
+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)