26784388084

Committed 01 Jun 2026 09:54PM UTC coverage: 74.937% (-0.01%) from 74.951%

Build # 26784388084

Build Type

Pull #211

github

Committed by

web-flow

Commit Message

Merge bd97e846f into f7be28bad

Pull Request Pull Request #211: Full 2p-ISR(2) difference density

Coverage Stats

1321 of 2062 branches covered (64.06%)

Branch coverage included in aggregate %.

100 of 118 new or added lines in 8 files covered. (84.75%)

16 existing lines in 4 files now uncovered.

8241 of 10698 relevant lines covered (77.03%)

167772.89 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

86.6

/adcc/adc_pp/state_diffdm_2p.py

#!/usr/bin/env python3
## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
## ---------------------------------------------------------------------
##
## Copyright (C) 2026 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 .. import block as b
from ..LazyMp import LazyMp
from ..AdcMethod import IsrMethod
from ..functions import einsum, zeros_like
from ..Intermediates import Intermediates
from ..AmplitudeVector import AmplitudeVector
from ..TwoParticleDensity import TwoParticleDensity
from ..NParticleOperator import OperatorSymmetry
from .util import (
    check_doubles_amplitudes, check_singles_amplitudes, check_triples_amplitudes
)

from typing import Optional, Union


def diffdm_isr0_2p(mp: LazyMp, amplitude: AmplitudeVector,
                   intermediates: Intermediates) -> TwoParticleDensity:
    check_singles_amplitudes([b.o, b.v], amplitude)
    u1 = amplitude.ph

    hf = mp.reference_state
    d_oo = zeros_like(hf.foo)
    d_oo.set_mask("ii", 1)

    dm = TwoParticleDensity(mp, symmetry=OperatorSymmetry.HERMITIAN)

    # one-particle ISR(0) diffdm
    p1_oo = -einsum("ia,la->il", u1, u1).evaluate()
    p1_vv = einsum("ka,kb->ab", u1, u1).evaluate()

    dm.oooo = (
        # N^4: O^4 / N^4: O^4
        - 4.0 * einsum("il,jk->ijkl", p1_oo, d_oo)
    ).antisymmetrise(0, 1).antisymmetrise(2, 3)
    dm.ovov = (
        # N^4: O^2V^2 / N^4: O^2V^2
        - 1.0 * einsum("ja,ib->iajb", u1, u1)
        # N^4: O^2V^2 / N^4: O^2V^2
        + 1.0 * einsum("ab,ij->iajb", p1_vv, d_oo)
    )
    return dm


def diffdm_isr1s_2p(mp: LazyMp, amplitude: AmplitudeVector,
                    intermediates: Intermediates) -> TwoParticleDensity:
    dm = diffdm_isr0_2p(mp, amplitude, intermediates)  # Get ISR(0) result
    u1 = amplitude.ph

    hf = mp.reference_state
    d_oo = zeros_like(hf.foo)
    d_oo.set_mask("ii", 1)

    t2 = mp.t2(b.oovv)
    # TODO move to intermediates!
    # one-particle ISR(0) diffdm
    p1_oo = -einsum("ia,la->il", u1, u1).evaluate()

    # ISR(2) ISR intermediate (TODO Move to intermediates)
    ru1 = einsum("ijab,jb->ia", t2, u1).evaluate()
    # new ones
    ru1_ooov = einsum("kc,ijac->ijka", u1, t2).evaluate()

    dm.oovv += (
        # N^4: O^2V^2 / N^4: O^2V^2
        + 4.0 * (
            einsum("ib,ja->ijab", ru1, u1)
        ).antisymmetrise(0, 1).antisymmetrise(2, 3)
        # N^5: O^3V^2 / N^4: O^2V^2
        + 2.0 * (
            einsum("ijka,kb->ijab", ru1_ooov, u1)
        ).antisymmetrise(2, 3)
        # N^5: O^3V^2 / N^4: O^2V^2
        - 2.0 * (
            einsum("jk,ikab->ijab", p1_oo, t2)
        ).antisymmetrise(0, 1)
    )

    return dm


def diffdm_isr1_2p(mp: LazyMp, amplitude: AmplitudeVector,
                   intermediates: Intermediates) -> TwoParticleDensity:
    dm = diffdm_isr1s_2p(mp, amplitude, intermediates)  # Get ISR(1)-s result

    try:
        # ISR(1)-d
        check_doubles_amplitudes([b.o, b.o, b.v, b.v], amplitude)
        u1, u2 = amplitude.ph, amplitude.pphh

        hf = mp.reference_state
        d_oo = zeros_like(hf.foo)
        d_oo.set_mask("ii", 1)

        dm.ooov += (
            # N^5: O^3V^2 / N^4: O^2V^2
            - 2.0 * einsum("kb,ijab->ijka", u1, u2)
            # N^4: O^2V^2 / N^4: O^2V^2
            - 4.0 * einsum(
                "ja,ik->ijka", einsum("lb,jlab->ja", u1, u2), d_oo
            ).antisymmetrise(0, 1)
        )

        dm.ovvv += (
            # N^5: O^2V^3 / N^4: O^1V^3
            - 2.0 * einsum("ja,ijbc->iabc", u1, u2)
        )
    except ValueError:
        # no doubles contribution
        pass
    return dm


def diffdm_isr2d_2p(mp: LazyMp, amplitude: AmplitudeVector,
                    intermediates: Intermediates) -> TwoParticleDensity:
    dm = diffdm_isr1_2p(mp, amplitude, intermediates)  # Get ISR(1) result
    check_doubles_amplitudes([b.o, b.o, b.v, b.v], amplitude)
    u1, u2 = amplitude.ph, amplitude.pphh
    hf = mp.reference_state
    d_oo = zeros_like(hf.foo)
    d_oo.set_mask("ii", 1)

    t2 = mp.t2(b.oovv)
    td2 = mp.td2(b.oovv)
    p0 = mp.mp2_diffdm

    dm.oooo += (
        # N^6: O^4V^2 / N^4: O^2V^2
        + 2.0 * einsum("ijab,klab->ijkl", u2, u2)
        # N^6: O^5V^1 / N^4: O^2V^2
        - 1.0 * einsum("klmc,ijmc->ijkl", einsum("mb,klbc->klmc", u1, t2),
                       einsum("ma,ijac->ijmc", u1, t2))
        + (
            # N^5: O^3V^2 / N^4: O^2V^2
            - 8.0 * einsum("ik,jl->ijkl", einsum("imab,kmab->ik", u2, u2), d_oo)
            # N^4: O^4 / N^4: O^4
            + 4.0 * einsum("il,jk->ijkl", einsum("ia,la->il", u1, u1), p0.oo)
            # N^6: O^5V^1 / N^4: O^2V^2
            + 4.0 * einsum("jlmc,ikmc->ijkl", einsum("lb,jmbc->jlmc", u1, t2),
                           einsum("ia,kmac->ikmc", u1, t2))
            # N^4: O^2V^2 / N^4: O^2V^2
            + 4.0 * einsum("il,jk->ijkl",
                           einsum("lc,ic->il", einsum("nb,lnbc->lc", u1, t2),
                                  einsum("ma,imac->ic", u1, t2)), d_oo)
            # N^5: O^3V^2 / N^4: O^2V^2
            + 4.0 * einsum("ik,jl->ijkl", einsum("kmnc,imnc->ik",
                                                 einsum("nb,kmbc->kmnc", u1, t2),
                                                 einsum("na,imac->imnc", u1, t2)),
                           d_oo)
            # N^5: O^3V^2 / N^4: O^2V^2
            + 2.0 * einsum("ik,jl->ijkl",
                           einsum("inbc,knbc->ik",
                                  einsum("mn,imbc->inbc",
                                         einsum("ia,la->il", u1, u1), t2), t2),
                           d_oo)
        ).antisymmetrise(0, 1).antisymmetrise(2, 3)
        + (
            # N^5: O^3V^2 / N^4: O^2V^2
            + 4.0 * einsum("iklb,jb->ijkl",
                           einsum("ic,klbc->iklb", einsum("ma,imac->ic", u1, t2),
                                  t2), u1)
            # N^6: O^4V^2 / N^4: O^2V^2
            + 2.0 * einsum("jklm,im->ijkl",
                           einsum("jmbc,klbc->jklm", t2, t2),
                           einsum("ia,la->il", u1, u1))
        ).antisymmetrise(0, 1).symmetrise([(0, 2), (1, 3)])
        + (
            # N^4: O^4 / N^4: O^4
            + 4.0 * einsum("il,jk->ijkl", einsum("im,lm->il",
                                                 einsum("ia,la->il", u1, u1),
                                                 p0.oo), d_oo)
            # N^4: O^2V^2 / N^4: O^2V^2
            + 4.0 * einsum("ik,jl->ijkl",
                           einsum("kb,ib->ik",
                                  einsum("mc,kmbc->kb",
                                         einsum("na,mnac->mc", u1, t2), t2), u1),
                           d_oo)
        ).antisymmetrise(0, 1).antisymmetrise(2, 3).symmetrise([(0, 2), (1, 3)])
    )
    dm.ooov += (
        # N^6: O^4V^2 / N^4: O^2V^2
        + 1.0 * einsum("ijkl,la->ijka", einsum("klbc,ijbc->ijkl", u2, t2), u1)
        # N^5: O^3V^2 / N^4: O^2V^2
        + 2.0 * einsum("kb,ijab->ijka", einsum("lc,klbc->kb", u1, u2), t2)
        + (
            # N^5: O^3V^2 / N^4: O^2V^2
            + 2.0 * einsum("jk,ia->ijka", einsum("klbc,jlbc->jk", u2, t2), u1)
            # N^6: O^4V^2 / N^4: O^2V^2
            - 4.0 * einsum("jklb,ilab->ijka", einsum("jc,klbc->jklb", u1, u2), t2)
            # N^5: O^3V^2 / N^4: O^2V^2
            + 2.0 * einsum("ja,ik->ijka",
                           einsum("jlmb,lmab->ja",
                                  einsum("jc,lmbc->jlmb", u1, u2), t2), d_oo)
            # N^4: O^2V^2 / N^4: O^2V^2
            - 4.0 * einsum("ia,jk->ijka",
                           einsum("lb,ilab->ia",
                                  einsum("mc,lmbc->lb", u1, u2), t2), d_oo)
            # N^5: O^3V^2 / N^4: O^2V^2
            + 2.0 * einsum("ia,jk->ijka",
                           einsum("il,la->ia",
                                  einsum("lmbc,imbc->il", u2, t2), u1), d_oo)
            # N^4: O^3V^1 / N^4: O^3V^1
            + 2.0 * einsum("jk,ia->ijka", einsum("kb,jb->jk", u1, p0.ov), u1)
            # N^4: O^3V^1 / N^4: O^3V^1
            + 2.0 * einsum("jk,ia->ijka", einsum("jb,kb->jk", u1, u1), p0.ov)
            # N^4: O^3V^1 / N^4: O^3V^1
            + 2.0 * einsum("ia,jk->ijka",
                           einsum("il,la->ia",
                                  einsum("lb,ib->il", u1, p0.ov), u1), d_oo)
            # N^4: O^3V^1 / N^4: O^3V^1
            + 2.0 * einsum("ia,jk->ijka",
                           einsum("il,la->ia", einsum("ia,la->il", u1, u1),
                                  p0.ov), d_oo)
        ).antisymmetrise(0, 1)
    )
    dm.oovv += (
        # N^4: O^2V^2 / N^4: O^2V^2
        + 4.0 * einsum("ib,ja->ijab",
                       einsum("kc,ikbc->ib", u1, td2),
                       u1).antisymmetrise(0, 1).antisymmetrise(2, 3)
        # N^5: O^3V^2 / N^4: O^2V^2
        + 2.0 * einsum("ijka,kb->ijab",
                       einsum("kc,ijac->ijka", u1, td2), u1).antisymmetrise(2, 3)
        # N^5: O^3V^2 / N^4: O^2V^2
        + 2.0 * einsum("jk,ikab->ijab", einsum("jc,kc->jk", u1, u1), td2)
    )
    dm.ovov += (
        # N^6: O^3V^3 / N^4: O^2V^2
        - 4.0 * einsum("jkac,ikbc->iajb", u2, u2)
        # N^5: O^2V^3 / N^4: O^2V^2
        + 2.0 * einsum("ab,ij->iajb", einsum("klac,klbc->ab", u2, u2), d_oo)
        # N^6: O^4V^2 / N^4: O^2V^2
        + 1.0 * einsum("iklb,jkla->iajb", einsum("ld,ikbd->iklb", u1, t2),
                       einsum("lc,jkac->jkla", u1, t2))
        # N^4: O^2V^2 / N^4: O^2V^2
        + 1.0 * einsum("ab,ij->iajb", einsum("ka,kb->ab", u1, u1), p0.oo)
        # N^4: O^2V^2 / N^4: O^2V^2
        - 1.0 * einsum("ij,ab->iajb", einsum("ic,jc->ij", u1, u1), p0.vv)
        # N^6: O^3V^3 / N^4: O^2V^2
        + 1.0 * einsum("jlac,ilbc->iajb",
                       einsum("kl,jkac->jlac", einsum("kd,ld->kl", u1, u1), t2), t2)
        # N^6: O^4V^2 / N^4: O^2V^2
        + 0.5 * einsum("ijla,lb->iajb",
                       einsum("ijkl,ka->ijla",
                              einsum("ikcd,jlcd->ijkl", t2, t2), u1), u1)
        # N^6: O^4V^2 / N^4: O^2V^2
        + 0.5 * einsum("jklb,ikla->iajb",
                       einsum("jd,klbd->jklb", u1, t2),
                       einsum("ic,klac->ikla", u1, t2))
        # N^4: O^2V^2 / N^4: O^2V^2
        - 1.0 * einsum("ib,ja->iajb",
                       einsum("kd,ikbd->ib", u1, t2), einsum("lc,jlac->ja", u1, t2))
        # N^5: O^2V^3 / N^4: O^2V^2
        + 1.0 * einsum("ab,ij->iajb",
                       einsum("lmad,lmbd->ab",
                              einsum("km,klad->lmad",
                                     einsum("kc,mc->km", u1, u1), t2), t2), d_oo)
        # N^4: O^2V^2 / N^4: O^2V^2
        - 1.0 * einsum("ab,ij->iajb",
                       einsum("lb,la->ab",
                              einsum("kd,klbd->lb", u1, t2),
                              einsum("mc,lmac->la", u1, t2)), d_oo)
        # N^5: O^3V^2 / N^4: O^2V^2
        - 0.5 * einsum("ab,ij->iajb",
                       einsum("klmb,klma->ab",
                              einsum("md,klbd->klmb", u1, t2),
                              einsum("mc,klac->klma", u1, t2)), d_oo)
        + (
            # N^4: O^2V^2 / N^4: O^2V^2
            + 1.0 * einsum("ib,ja->iajb", einsum("ic,bc->ib", u1, p0.vv), u1)
            # N^4: O^2V^2 / N^4: O^2V^2
            - 1.0 * einsum("ib,ja->iajb", einsum("kb,ik->ib", u1, p0.oo), u1)
            # N^5: O^3V^2 / N^4: O^2V^2
            + 2.0 * einsum("ijlb,la->iajb",
                           einsum("jc,ilbc->ijlb",
                                  einsum("kd,jkcd->jc", u1, t2), t2), u1)
            # N^5: O^3V^2 / N^4: O^2V^2
            + 2.0 * einsum("kb,ijka->iajb",
                           einsum("ld,klbd->kb", u1, t2),
                           einsum("ic,jkac->ijka", u1, t2))
            # N^6: O^4V^2 / N^4: O^2V^2
            - 2.0 * einsum("ijkb,ka->iajb",
                           einsum("jklc,ilbc->ijkb",
                                  einsum("kd,jlcd->jklc", u1, t2), t2), u1)
            # N^6: O^4V^2 / N^4: O^2V^2
            - 2.0 * einsum("ijlb,la->iajb",
                           einsum("ijkc,klbc->ijlb",
                                  einsum("id,jkcd->ijkc", u1, t2), t2), u1)
            # N^6: O^3V^3 / N^4: O^2V^2
            - 2.0 * einsum("ikbc,jkac->iajb",
                           einsum("il,klbc->ikbc",
                                  einsum("id,ld->il", u1, u1), t2), t2)
            # N^4: O^2V^2 / N^4: O^2V^2
            - 1.0 * einsum("ib,ja->iajb",
                           einsum("kc,ikbc->ib",
                                  einsum("ld,klcd->kc", u1, t2), t2), u1)
            # N^4: O^2V^2 / N^4: O^2V^2
            - 1.0 * einsum("ab,ij->iajb",
                           einsum("kb,ka->ab",
                                  einsum("kc,bc->kb", u1, p0.vv), u1), d_oo)
            # N^4: O^2V^2 / N^4: O^2V^2
            - 1.0 * einsum("ab,ij->iajb",
                           einsum("mb,ma->ab",
                                  einsum("ld,lmbd->mb",
                                         einsum("kc,klcd->ld", u1, t2),
                                         t2), u1), d_oo)
        ).symmetrise([(0, 2), (1, 3)])
    )
    dm.ovvv += (
        # N^6: O^3V^3 / N^4: O^1V^3
        + 1.0 * einsum("ijka,jkbc->iabc", einsum("id,jkad->ijka", u1, u2), t2)
        # N^5: O^2V^3 / N^4: O^1V^3
        + 2.0 * einsum("ja,ijbc->iabc", einsum("kd,jkad->ja", u1, u2), t2)
        + (
            # N^5: O^2V^3 / N^4: O^1V^3
            + 2.0 * einsum("ac,ib->iabc", einsum("jkad,jkcd->ac", u2, t2), u1)
            # N^6: O^3V^3 / N^4: O^1V^3
            - 4.0 * einsum("ikab,kc->iabc", einsum("jkad,ijbd->ikab", u2, t2), u1)
            # N^4: O^1V^3 / N^4: O^1V^3
            + 2.0 * einsum("ac,ib->iabc", einsum("ja,jc->ac", u1, p0.ov), u1)
            # N^4: O^1V^3 / N^4: O^1V^3
            + 2.0 * einsum("ac,ib->iabc", einsum("ja,jc->ac", u1, u1), p0.ov)
        ).antisymmetrise(2, 3)
    )
    dm.vvvv += (
        # N^6: O^2V^4 / N^4: V^4
        + 2 * einsum("ijab,ijcd->abcd", u2, u2)
        # N^6: O^2V^4 / N^4: V^4
        - 1 * einsum("jkab,jkcd->abcd",
                     einsum("ij,ikab->jkab", einsum("ie,je->ij", u1, u1), t2), t2)
        + (
            # N^4: V^4 / N^4: V^4
            + 4.0 * einsum("ac,bd->abcd", einsum("ia,ic->ac", u1, u1), p0.vv)
            # N^6: O^3V^3 / N^4: V^4
            + 4.0 * einsum("jabc,jd->abcd",
                           einsum("ijbc,ia->jabc",
                                  einsum("jkbe,ikce->ijbc", t2, t2), u1), u1)
        ).antisymmetrise(0, 1).antisymmetrise(2, 3)
        + (
            # N^5: O^1V^4 / N^4: V^4
            + 4.0 * einsum("ka,kbcd->abcd",
                           einsum("ie,ikae->ka", u1, t2),
                           einsum("jb,jkcd->kbcd", u1, t2))
            # N^6: O^3V^3 / N^4: V^4
            + 2.0 * einsum("ibcd,ia->abcd",
                           einsum("ijkb,jkcd->ibcd",
                                  einsum("ie,jkbe->ijkb", u1, t2), t2), u1)
        ).antisymmetrise(0, 1).symmetrise([(0, 2), (1, 3)])
    )
    return dm


def diffdm_isr2_2p(mp: LazyMp, amplitude: AmplitudeVector,
                   intermediates: Intermediates) -> TwoParticleDensity:
    dm = diffdm_isr2d_2p(mp, amplitude, intermediates)  # Get ISR(1) result
    # evaluate additional contributions from the S-T block
    # if the vector has a triples component
    try:
        check_triples_amplitudes([b.o, b.o, b.o, b.v, b.v, b.v], amplitude)
    except ValueError:
        return dm
    u1, u3 = amplitude.ph, amplitude.ppphhh
    # N^6: O^3V^3 / N^6: O^3V^3
    dm.oovv += - 6 * einsum("kc,ijkabc->ijab", u1, u3)
    return dm


# dict controlling the dispatch of the state_diffdm function
DISPATCH = {
    "isr0": diffdm_isr0_2p,
    "isr1s": diffdm_isr1s_2p,
    "isr1": diffdm_isr1_2p,
    "isr2d": diffdm_isr2d_2p,
    "isr2": diffdm_isr2_2p,
}


def state_diffdm_2p(method: Union[str, IsrMethod], ground_state: LazyMp,
                    amplitude: AmplitudeVector,
                    intermediates: Optional[Intermediates] = None):
    """
    Compute the two-particle difference density matrix of an excited state
    in the MO basis.

    Parameters
    ----------
    method : str, IsrMethod
        The method to use for the computation (e.g. "isr2")
    ground_state : LazyMp
        The ground state upon which the excitation was based
    amplitude : AmplitudeVector
        The amplitude vector
    intermediates : adcc.Intermediates
        Intermediates from the ADC calculation to reuse
    """
    if not isinstance(method, IsrMethod):
        method = IsrMethod(method)
    if not isinstance(ground_state, LazyMp):
        raise TypeError("ground_state should be a LazyMp object.")
    if not isinstance(amplitude, AmplitudeVector):
        raise TypeError("amplitude should be an AmplitudeVector object.")
    if intermediates is None:
        intermediates = Intermediates(ground_state)

    if method.name not in DISPATCH:
        raise NotImplementedError("state_diffdm_2p is not implemented "
                                  f"for {method.name}.")
    else:
        ret = DISPATCH[method.name](ground_state, amplitude, intermediates)
        return ret.evaluate()

1	#!/usr/bin/env python3
2	## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
3	## ---------------------------------------------------------------------
4	##
5	## Copyright (C) 2026 by the adcc authors
6	##
7	## This file is part of adcc.
8	##
9	## adcc is free software: you can redistribute it and/or modify
10	## it under the terms of the GNU General Public License as published
11	## by the Free Software Foundation, either version 3 of the License, or
12	## (at your option) any later version.
13	##
14	## adcc is distributed in the hope that it will be useful,
15	## but WITHOUT ANY WARRANTY; without even the implied warranty of
16	## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17	## GNU General Public License for more details.
18	##
19	## You should have received a copy of the GNU General Public License
20	## along with adcc. If not, see <http://www.gnu.org/licenses/>.
21	##
22	## ---------------------------------------------------------------------
23	from .. import block as b	2✔
24	from ..LazyMp import LazyMp	2✔
25	from ..AdcMethod import IsrMethod	2✔
26	from ..functions import einsum, zeros_like	2✔
27	from ..Intermediates import Intermediates	2✔
28	from ..AmplitudeVector import AmplitudeVector	2✔
29	from ..TwoParticleDensity import TwoParticleDensity	2✔
30	from ..NParticleOperator import OperatorSymmetry	2✔
31	from .util import (	2✔
32	check_doubles_amplitudes, check_singles_amplitudes, check_triples_amplitudes
33	)
34
35	from typing import Optional, Union	2✔
36
37
38	def diffdm_isr0_2p(mp: LazyMp, amplitude: AmplitudeVector,	2✔
39	intermediates: Intermediates) -> TwoParticleDensity:
40	check_singles_amplitudes([b.o, b.v], amplitude)	2✔
41	u1 = amplitude.ph	2✔
42
43	hf = mp.reference_state	2✔
44	d_oo = zeros_like(hf.foo)	2✔
45	d_oo.set_mask("ii", 1)	2✔
46
47	dm = TwoParticleDensity(mp, symmetry=OperatorSymmetry.HERMITIAN)	2✔
48
49	# one-particle ISR(0) diffdm
50	p1_oo = -einsum("ia,la->il", u1, u1).evaluate()	2✔
51	p1_vv = einsum("ka,kb->ab", u1, u1).evaluate()	2✔
52
53	dm.oooo = (	2✔
54	# N^4: O^4 / N^4: O^4
55	- 4.0 * einsum("il,jk->ijkl", p1_oo, d_oo)
56	).antisymmetrise(0, 1).antisymmetrise(2, 3)
57	dm.ovov = (	2✔
58	# N^4: O^2V^2 / N^4: O^2V^2
59	- 1.0 * einsum("ja,ib->iajb", u1, u1)
60	# N^4: O^2V^2 / N^4: O^2V^2
61	+ 1.0 * einsum("ab,ij->iajb", p1_vv, d_oo)
62	)
63	return dm	2✔
64
65
66	def diffdm_isr1s_2p(mp: LazyMp, amplitude: AmplitudeVector,	2✔
67	intermediates: Intermediates) -> TwoParticleDensity:
68	dm = diffdm_isr0_2p(mp, amplitude, intermediates) # Get ISR(0) result	2✔
69	u1 = amplitude.ph	2✔
70
71	hf = mp.reference_state	2✔
72	d_oo = zeros_like(hf.foo)	2✔
73	d_oo.set_mask("ii", 1)	2✔
74
75	t2 = mp.t2(b.oovv)	2✔
76	# TODO move to intermediates!
77	# one-particle ISR(0) diffdm
78	p1_oo = -einsum("ia,la->il", u1, u1).evaluate()	2✔
79
80	# ISR(2) ISR intermediate (TODO Move to intermediates)
81	ru1 = einsum("ijab,jb->ia", t2, u1).evaluate()	2✔
82	# new ones
83	ru1_ooov = einsum("kc,ijac->ijka", u1, t2).evaluate()	2✔
84
85	dm.oovv += (	2✔
86	# N^4: O^2V^2 / N^4: O^2V^2
87	+ 4.0 * (
88	einsum("ib,ja->ijab", ru1, u1)
89	).antisymmetrise(0, 1).antisymmetrise(2, 3)
90	# N^5: O^3V^2 / N^4: O^2V^2
91	+ 2.0 * (
92	einsum("ijka,kb->ijab", ru1_ooov, u1)
93	).antisymmetrise(2, 3)
94	# N^5: O^3V^2 / N^4: O^2V^2
95	- 2.0 * (
96	einsum("jk,ikab->ijab", p1_oo, t2)
97	).antisymmetrise(0, 1)
98	)
99
100	return dm	2✔
101
102
103	def diffdm_isr1_2p(mp: LazyMp, amplitude: AmplitudeVector,	2✔
104	intermediates: Intermediates) -> TwoParticleDensity:
105	dm = diffdm_isr1s_2p(mp, amplitude, intermediates) # Get ISR(1)-s result	2✔
106
107	try:	2✔
108	# ISR(1)-d
109	check_doubles_amplitudes([b.o, b.o, b.v, b.v], amplitude)	2✔
110	u1, u2 = amplitude.ph, amplitude.pphh	2✔
111
112	hf = mp.reference_state	2✔
113	d_oo = zeros_like(hf.foo)	2✔
114	d_oo.set_mask("ii", 1)	2✔
115
116	dm.ooov += (	2✔
117	# N^5: O^3V^2 / N^4: O^2V^2
118	- 2.0 * einsum("kb,ijab->ijka", u1, u2)
119	# N^4: O^2V^2 / N^4: O^2V^2
120	- 4.0 * einsum(
121	"ja,ik->ijka", einsum("lb,jlab->ja", u1, u2), d_oo
122	).antisymmetrise(0, 1)
123	)
124
125	dm.ovvv += (	2✔
126	# N^5: O^2V^3 / N^4: O^1V^3
127	- 2.0 * einsum("ja,ijbc->iabc", u1, u2)
128	)
129	except ValueError:	×
130	# no doubles contribution
131	pass	×
132	return dm	2✔
133
134
135	def diffdm_isr2d_2p(mp: LazyMp, amplitude: AmplitudeVector,	2✔
136	intermediates: Intermediates) -> TwoParticleDensity:
137	dm = diffdm_isr1_2p(mp, amplitude, intermediates) # Get ISR(1) result	2✔
138	check_doubles_amplitudes([b.o, b.o, b.v, b.v], amplitude)	2✔
139	u1, u2 = amplitude.ph, amplitude.pphh	2✔
140	hf = mp.reference_state	2✔
141	d_oo = zeros_like(hf.foo)	2✔
142	d_oo.set_mask("ii", 1)	2✔
143
144	t2 = mp.t2(b.oovv)	2✔
145	td2 = mp.td2(b.oovv)	2✔
146	p0 = mp.mp2_diffdm	2✔
147
148	dm.oooo += (	2✔
149	# N^6: O^4V^2 / N^4: O^2V^2
150	+ 2.0 * einsum("ijab,klab->ijkl", u2, u2)
151	# N^6: O^5V^1 / N^4: O^2V^2
152	- 1.0 * einsum("klmc,ijmc->ijkl", einsum("mb,klbc->klmc", u1, t2),
153	einsum("ma,ijac->ijmc", u1, t2))
154	+ (
155	# N^5: O^3V^2 / N^4: O^2V^2
156	- 8.0 * einsum("ik,jl->ijkl", einsum("imab,kmab->ik", u2, u2), d_oo)
157	# N^4: O^4 / N^4: O^4
158	+ 4.0 * einsum("il,jk->ijkl", einsum("ia,la->il", u1, u1), p0.oo)
159	# N^6: O^5V^1 / N^4: O^2V^2
160	+ 4.0 * einsum("jlmc,ikmc->ijkl", einsum("lb,jmbc->jlmc", u1, t2),
161	einsum("ia,kmac->ikmc", u1, t2))
162	# N^4: O^2V^2 / N^4: O^2V^2
163	+ 4.0 * einsum("il,jk->ijkl",
164	einsum("lc,ic->il", einsum("nb,lnbc->lc", u1, t2),
165	einsum("ma,imac->ic", u1, t2)), d_oo)
166	# N^5: O^3V^2 / N^4: O^2V^2
167	+ 4.0 * einsum("ik,jl->ijkl", einsum("kmnc,imnc->ik",
168	einsum("nb,kmbc->kmnc", u1, t2),
169	einsum("na,imac->imnc", u1, t2)),
170	d_oo)
171	# N^5: O^3V^2 / N^4: O^2V^2
172	+ 2.0 * einsum("ik,jl->ijkl",
173	einsum("inbc,knbc->ik",
174	einsum("mn,imbc->inbc",
175	einsum("ia,la->il", u1, u1), t2), t2),
176	d_oo)
177	).antisymmetrise(0, 1).antisymmetrise(2, 3)
178	+ (
179	# N^5: O^3V^2 / N^4: O^2V^2
180	+ 4.0 * einsum("iklb,jb->ijkl",
181	einsum("ic,klbc->iklb", einsum("ma,imac->ic", u1, t2),
182	t2), u1)
183	# N^6: O^4V^2 / N^4: O^2V^2
184	+ 2.0 * einsum("jklm,im->ijkl",
185	einsum("jmbc,klbc->jklm", t2, t2),
186	einsum("ia,la->il", u1, u1))
187	).antisymmetrise(0, 1).symmetrise([(0, 2), (1, 3)])
188	+ (
189	# N^4: O^4 / N^4: O^4
190	+ 4.0 * einsum("il,jk->ijkl", einsum("im,lm->il",
191	einsum("ia,la->il", u1, u1),
192	p0.oo), d_oo)
193	# N^4: O^2V^2 / N^4: O^2V^2
194	+ 4.0 * einsum("ik,jl->ijkl",
195	einsum("kb,ib->ik",
196	einsum("mc,kmbc->kb",
197	einsum("na,mnac->mc", u1, t2), t2), u1),
198	d_oo)
199	).antisymmetrise(0, 1).antisymmetrise(2, 3).symmetrise([(0, 2), (1, 3)])
200	)
201	dm.ooov += (	2✔
202	# N^6: O^4V^2 / N^4: O^2V^2
203	+ 1.0 * einsum("ijkl,la->ijka", einsum("klbc,ijbc->ijkl", u2, t2), u1)
204	# N^5: O^3V^2 / N^4: O^2V^2
205	+ 2.0 * einsum("kb,ijab->ijka", einsum("lc,klbc->kb", u1, u2), t2)
206	+ (
207	# N^5: O^3V^2 / N^4: O^2V^2
208	+ 2.0 * einsum("jk,ia->ijka", einsum("klbc,jlbc->jk", u2, t2), u1)
209	# N^6: O^4V^2 / N^4: O^2V^2
210	- 4.0 * einsum("jklb,ilab->ijka", einsum("jc,klbc->jklb", u1, u2), t2)
211	# N^5: O^3V^2 / N^4: O^2V^2
212	+ 2.0 * einsum("ja,ik->ijka",
213	einsum("jlmb,lmab->ja",
214	einsum("jc,lmbc->jlmb", u1, u2), t2), d_oo)
215	# N^4: O^2V^2 / N^4: O^2V^2
216	- 4.0 * einsum("ia,jk->ijka",
217	einsum("lb,ilab->ia",
218	einsum("mc,lmbc->lb", u1, u2), t2), d_oo)
219	# N^5: O^3V^2 / N^4: O^2V^2
220	+ 2.0 * einsum("ia,jk->ijka",
221	einsum("il,la->ia",
222	einsum("lmbc,imbc->il", u2, t2), u1), d_oo)
223	# N^4: O^3V^1 / N^4: O^3V^1
224	+ 2.0 * einsum("jk,ia->ijka", einsum("kb,jb->jk", u1, p0.ov), u1)
225	# N^4: O^3V^1 / N^4: O^3V^1
226	+ 2.0 * einsum("jk,ia->ijka", einsum("jb,kb->jk", u1, u1), p0.ov)
227	# N^4: O^3V^1 / N^4: O^3V^1
228	+ 2.0 * einsum("ia,jk->ijka",
229	einsum("il,la->ia",
230	einsum("lb,ib->il", u1, p0.ov), u1), d_oo)
231	# N^4: O^3V^1 / N^4: O^3V^1
232	+ 2.0 * einsum("ia,jk->ijka",
233	einsum("il,la->ia", einsum("ia,la->il", u1, u1),
234	p0.ov), d_oo)
235	).antisymmetrise(0, 1)
236	)
237	dm.oovv += (	2✔
238	# N^4: O^2V^2 / N^4: O^2V^2
239	+ 4.0 * einsum("ib,ja->ijab",
240	einsum("kc,ikbc->ib", u1, td2),
241	u1).antisymmetrise(0, 1).antisymmetrise(2, 3)
242	# N^5: O^3V^2 / N^4: O^2V^2
243	+ 2.0 * einsum("ijka,kb->ijab",
244	einsum("kc,ijac->ijka", u1, td2), u1).antisymmetrise(2, 3)
245	# N^5: O^3V^2 / N^4: O^2V^2
246	+ 2.0 * einsum("jk,ikab->ijab", einsum("jc,kc->jk", u1, u1), td2)
247	)
248	dm.ovov += (	2✔
249	# N^6: O^3V^3 / N^4: O^2V^2
250	- 4.0 * einsum("jkac,ikbc->iajb", u2, u2)
251	# N^5: O^2V^3 / N^4: O^2V^2
252	+ 2.0 * einsum("ab,ij->iajb", einsum("klac,klbc->ab", u2, u2), d_oo)
253	# N^6: O^4V^2 / N^4: O^2V^2
254	+ 1.0 * einsum("iklb,jkla->iajb", einsum("ld,ikbd->iklb", u1, t2),
255	einsum("lc,jkac->jkla", u1, t2))
256	# N^4: O^2V^2 / N^4: O^2V^2
257	+ 1.0 * einsum("ab,ij->iajb", einsum("ka,kb->ab", u1, u1), p0.oo)
258	# N^4: O^2V^2 / N^4: O^2V^2
259	- 1.0 * einsum("ij,ab->iajb", einsum("ic,jc->ij", u1, u1), p0.vv)
260	# N^6: O^3V^3 / N^4: O^2V^2
261	+ 1.0 * einsum("jlac,ilbc->iajb",
262	einsum("kl,jkac->jlac", einsum("kd,ld->kl", u1, u1), t2), t2)
263	# N^6: O^4V^2 / N^4: O^2V^2
264	+ 0.5 * einsum("ijla,lb->iajb",
265	einsum("ijkl,ka->ijla",
266	einsum("ikcd,jlcd->ijkl", t2, t2), u1), u1)
267	# N^6: O^4V^2 / N^4: O^2V^2
268	+ 0.5 * einsum("jklb,ikla->iajb",
269	einsum("jd,klbd->jklb", u1, t2),
270	einsum("ic,klac->ikla", u1, t2))
271	# N^4: O^2V^2 / N^4: O^2V^2
272	- 1.0 * einsum("ib,ja->iajb",
273	einsum("kd,ikbd->ib", u1, t2), einsum("lc,jlac->ja", u1, t2))
274	# N^5: O^2V^3 / N^4: O^2V^2
275	+ 1.0 * einsum("ab,ij->iajb",
276	einsum("lmad,lmbd->ab",
277	einsum("km,klad->lmad",
278	einsum("kc,mc->km", u1, u1), t2), t2), d_oo)
279	# N^4: O^2V^2 / N^4: O^2V^2
280	- 1.0 * einsum("ab,ij->iajb",
281	einsum("lb,la->ab",
282	einsum("kd,klbd->lb", u1, t2),
283	einsum("mc,lmac->la", u1, t2)), d_oo)
284	# N^5: O^3V^2 / N^4: O^2V^2
285	- 0.5 * einsum("ab,ij->iajb",
286	einsum("klmb,klma->ab",
287	einsum("md,klbd->klmb", u1, t2),
288	einsum("mc,klac->klma", u1, t2)), d_oo)
289	+ (
290	# N^4: O^2V^2 / N^4: O^2V^2
291	+ 1.0 * einsum("ib,ja->iajb", einsum("ic,bc->ib", u1, p0.vv), u1)
292	# N^4: O^2V^2 / N^4: O^2V^2
293	- 1.0 * einsum("ib,ja->iajb", einsum("kb,ik->ib", u1, p0.oo), u1)
294	# N^5: O^3V^2 / N^4: O^2V^2
295	+ 2.0 * einsum("ijlb,la->iajb",
296	einsum("jc,ilbc->ijlb",
297	einsum("kd,jkcd->jc", u1, t2), t2), u1)
298	# N^5: O^3V^2 / N^4: O^2V^2
299	+ 2.0 * einsum("kb,ijka->iajb",
300	einsum("ld,klbd->kb", u1, t2),
301	einsum("ic,jkac->ijka", u1, t2))
302	# N^6: O^4V^2 / N^4: O^2V^2
303	- 2.0 * einsum("ijkb,ka->iajb",
304	einsum("jklc,ilbc->ijkb",
305	einsum("kd,jlcd->jklc", u1, t2), t2), u1)
306	# N^6: O^4V^2 / N^4: O^2V^2
307	- 2.0 * einsum("ijlb,la->iajb",
308	einsum("ijkc,klbc->ijlb",
309	einsum("id,jkcd->ijkc", u1, t2), t2), u1)
310	# N^6: O^3V^3 / N^4: O^2V^2
311	- 2.0 * einsum("ikbc,jkac->iajb",
312	einsum("il,klbc->ikbc",
313	einsum("id,ld->il", u1, u1), t2), t2)
314	# N^4: O^2V^2 / N^4: O^2V^2
315	- 1.0 * einsum("ib,ja->iajb",
316	einsum("kc,ikbc->ib",
317	einsum("ld,klcd->kc", u1, t2), t2), u1)
318	# N^4: O^2V^2 / N^4: O^2V^2
319	- 1.0 * einsum("ab,ij->iajb",
320	einsum("kb,ka->ab",
321	einsum("kc,bc->kb", u1, p0.vv), u1), d_oo)
322	# N^4: O^2V^2 / N^4: O^2V^2
323	- 1.0 * einsum("ab,ij->iajb",
324	einsum("mb,ma->ab",
325	einsum("ld,lmbd->mb",
326	einsum("kc,klcd->ld", u1, t2),
327	t2), u1), d_oo)
328	).symmetrise([(0, 2), (1, 3)])
329	)
330	dm.ovvv += (	2✔
331	# N^6: O^3V^3 / N^4: O^1V^3
332	+ 1.0 * einsum("ijka,jkbc->iabc", einsum("id,jkad->ijka", u1, u2), t2)
333	# N^5: O^2V^3 / N^4: O^1V^3
334	+ 2.0 * einsum("ja,ijbc->iabc", einsum("kd,jkad->ja", u1, u2), t2)
335	+ (
336	# N^5: O^2V^3 / N^4: O^1V^3
337	+ 2.0 * einsum("ac,ib->iabc", einsum("jkad,jkcd->ac", u2, t2), u1)
338	# N^6: O^3V^3 / N^4: O^1V^3
339	- 4.0 * einsum("ikab,kc->iabc", einsum("jkad,ijbd->ikab", u2, t2), u1)
340	# N^4: O^1V^3 / N^4: O^1V^3
341	+ 2.0 * einsum("ac,ib->iabc", einsum("ja,jc->ac", u1, p0.ov), u1)
342	# N^4: O^1V^3 / N^4: O^1V^3
343	+ 2.0 * einsum("ac,ib->iabc", einsum("ja,jc->ac", u1, u1), p0.ov)
344	).antisymmetrise(2, 3)
345	)
346	dm.vvvv += (	2✔
347	# N^6: O^2V^4 / N^4: V^4
348	+ 2 * einsum("ijab,ijcd->abcd", u2, u2)
349	# N^6: O^2V^4 / N^4: V^4
350	- 1 * einsum("jkab,jkcd->abcd",
351	einsum("ij,ikab->jkab", einsum("ie,je->ij", u1, u1), t2), t2)
352	+ (
353	# N^4: V^4 / N^4: V^4
354	+ 4.0 * einsum("ac,bd->abcd", einsum("ia,ic->ac", u1, u1), p0.vv)
355	# N^6: O^3V^3 / N^4: V^4
356	+ 4.0 * einsum("jabc,jd->abcd",
357	einsum("ijbc,ia->jabc",
358	einsum("jkbe,ikce->ijbc", t2, t2), u1), u1)
359	).antisymmetrise(0, 1).antisymmetrise(2, 3)
360	+ (
361	# N^5: O^1V^4 / N^4: V^4
362	+ 4.0 * einsum("ka,kbcd->abcd",
363	einsum("ie,ikae->ka", u1, t2),
364	einsum("jb,jkcd->kbcd", u1, t2))
365	# N^6: O^3V^3 / N^4: V^4
366	+ 2.0 * einsum("ibcd,ia->abcd",
367	einsum("ijkb,jkcd->ibcd",
368	einsum("ie,jkbe->ijkb", u1, t2), t2), u1)
369	).antisymmetrise(0, 1).symmetrise([(0, 2), (1, 3)])
370	)
371	return dm	2✔
372
373
374	def diffdm_isr2_2p(mp: LazyMp, amplitude: AmplitudeVector,	2✔
375	intermediates: Intermediates) -> TwoParticleDensity:
376	dm = diffdm_isr2d_2p(mp, amplitude, intermediates) # Get ISR(1) result	2✔
377	# evaluate additional contributions from the S-T block
378	# if the vector has a triples component
379	try:	2✔
380	check_triples_amplitudes([b.o, b.o, b.o, b.v, b.v, b.v], amplitude)	2✔
381	except ValueError:	2✔
382	return dm	2✔
NEW 383	u1, u3 = amplitude.ph, amplitude.ppphhh	×
384	# N^6: O^3V^3 / N^6: O^3V^3
NEW 385	dm.oovv += - 6 * einsum("kc,ijkabc->ijab", u1, u3)	×
NEW 386	return dm	×
387
388
389	# dict controlling the dispatch of the state_diffdm function
390	DISPATCH = {	2✔
391	"isr0": diffdm_isr0_2p,
392	"isr1s": diffdm_isr1s_2p,
393	"isr1": diffdm_isr1_2p,
394	"isr2d": diffdm_isr2d_2p,
395	"isr2": diffdm_isr2_2p,
396	}
397
398
399	def state_diffdm_2p(method: Union[str, IsrMethod], ground_state: LazyMp,	2✔
400	amplitude: AmplitudeVector,
401	intermediates: Optional[Intermediates] = None):
402	"""
403	Compute the two-particle difference density matrix of an excited state
404	in the MO basis.
405
406	Parameters
407	----------
408	method : str, IsrMethod
409	The method to use for the computation (e.g. "isr2")
410	ground_state : LazyMp
411	The ground state upon which the excitation was based
412	amplitude : AmplitudeVector
413	The amplitude vector
414	intermediates : adcc.Intermediates
415	Intermediates from the ADC calculation to reuse
416	"""
417	if not isinstance(method, IsrMethod):	2!
418	method = IsrMethod(method)	×
419	if not isinstance(ground_state, LazyMp):	2!
420	raise TypeError("ground_state should be a LazyMp object.")	×
421	if not isinstance(amplitude, AmplitudeVector):	2!
422	raise TypeError("amplitude should be an AmplitudeVector object.")	×
423	if intermediates is None:	2✔
424	intermediates = Intermediates(ground_state)	2✔
425
426	if method.name not in DISPATCH:	2!
427	raise NotImplementedError("state_diffdm_2p is not implemented "	×
428	f"for {method.name}.")
429	else:
430	ret = DISPATCH[method.name](ground_state, amplitude, intermediates)	2✔
431	return ret.evaluate()	2✔

adc-connect / adcc / 26784388084

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous