1.3. Implementation Examples
This section provides concrete examples of implementing and registering new Hamiltonian blocks using the protocols outlined in the previous sections. The code snippets below are pulled directly from the PyBEST source files.
1.3.1. Understanding the Directory Routing
The target subfolder within eff_ham/blocks/ is determined dynamically by sorting the block’s rank and index groups into occupied (o) or virtual (v) footprints. For the detailed index classification tables and step-by-step path generation rules, see Determine Directory and Create File.
1.3.2. Example 1: Dense Block Implementation (Protocol A)
This is a practical example of how to implement a standard, cacheable block. We are implementing the X_ad block here.
Indices:
a(virtual),d(virtual) translates to thevvdirectory.File Name:
x_ad.py
#
# X_ad = f_ad + L<ak|dc> * t_k^c - (f_md + L<mk|dc> * t_k^c) * t_m^a - L<mk|dc> * t_mk^ac
#
def get_X_ad_block(self, force_update: bool = False) -> DenseTwoIndex:
"""Get or compute the Hamiltonian X_ad block.
Args:
force_update (bool): If True, forces recomputation of X_ad even if cached.
"""
if "X_ad" in self._cache and not force_update:
return self._cache.load("X_ad")
# Load integrals and amplitudes
# NOTE: we assume that all objects are present as
# we ensure their existence when `self` is instantiated
f_vv = self._ham.get("fock_vv")
f_ov = self._ham.get("fock_ov")
e_vovv = self._ham.get("eri_vovv")
e_oovv = self._ham.get("eri_oovv")
t_1 = self.rcc_iodata.t_1
t_2 = self.rcc_iodata.t_2
# f_ad -> f_vv
X_ad = self.init_block("X_ad", alloc=f_vv.copy)
# L<ak|dc> * t_k^c -> ( 2<ab|cd> - <ab|dc> ) * t_b^d
e_vovv.contract("abcd,bd->ac", t_1, out=X_ad, factor=2.0)
e_vovv.contract("abcd,bc->ad", t_1, out=X_ad, factor=-1.0)
# - (f_md + L<mk|dc> * t_k^c) * t_m^a -> - (e_ov + [ 2<ab|cd> - <ab|dc> ] * t_b^d) * t_a^e
#
# - f_md * t_m^a
f_ov.contract("ab,ac->cb", t_1, out=X_ad, factor=-1.0)
# - L<mk|dc> * t_mk^ac = ad -> - ( 2<ab|cd> - <ab|dc> ) * t_ab^ed
e_oovv.contract("abcd,aebd->ec", t_2, out=X_ad, factor=-2.0)
e_oovv.contract("abcd,aebc->ed", t_2, out=X_ad, factor=1.0)
if self.disconnected:
# - L<mk|dc> * t_k^c = md
tmp = e_oovv.contract("abcd,bd->ac", t_1, factor=2.0)
e_oovv.contract("abcd,bc->ad", t_1, out=tmp, factor=-1.0)
# t_m^a -> md * ma = ad
tmp.contract("ab,ac->cb", t_1, out=X_ad, factor=-1.0)
return X_ad
1.3.3. Example 2: Special Block Implementation (Protocol B)
This is an example of implementing Special blocks. Because Special blocks bypass the cache entirely regardless of the linear algebra factory, they only require the Protocol B implementation. We are implementing the X_iAjb_sigma and X_iAJB_sigma blocks.
Indices:
i, a, j, b(occupied, virtual, occupied, virtual) translates to theovovdirectory.File Name:
x_iajb.py
Notice how both distinct spin blocks are implemented in the same lowercase file, but their exact casing is preserved in their respective method signatures. Also note how they compute directly into the developer-provided sigma output tensor.
#
# - <jk|ba> * t_kj^ab - 1/4<ij||db> * (t_ij^db - t_ij^bd)
#
def get_X_iAjb_sigma_block_cholesky(self, sigma: DenseFourIndex) -> None:
"""Compute the diagonal term involving the oovv block for sigma_iAjb.
This corresponds to the 'X_iAjb_sigma' special block contraction.
Args:
sigma (DenseFourIndex): The current value of the iAjb sigma diagonal.
"""
# Load integrals and amplitudes
# NOTE: we assume that all objects are present as
# we ensure their existence when `self` is instantiated
e_oovv = self._ham.get("eri_oovv")
t_2 = self.rcc_iodata.t_2
# - <jk|ba> * t_kj^ab
# jkba,kajb->ajb
# abcd,bdac->dac
tmp_ajb = e_oovv.contract("abcd,bdac->dac", t_2, factor=-1.0)
tmp_ajb.expand("bcd->abcd", out=sigma)
# -1/4<ij||db> * (t_ij^db - t_ij^bd) -> =-1/4(<ij|db> - <ij|bd>) * (t_ij^db - t_ij^bd)
# (1) -1/4<ij|db> * t_ij^db
# ijdb,idjb->ijb
# abcd,acbd->abd
tmp_ijb = e_oovv.contract("abcd,acbd->abd", t_2, factor=-0.25)
# (2) -1/4<ij|db> * (-t_ij^bd)
# ijdb,ibjd->ijb
# abcd,adbc->abc
e_oovv.contract("abcd,adbc->abc", t_2, factor=0.25, out=tmp_ijb)
# (3) -1/4(-<ij|bd>) * t_ij^db
# ijbd,idjb->ijb
# abcd,adbc->abc
e_oovv.contract("abcd,adbc->abc", t_2, factor=0.25, out=tmp_ijb)
# (4) -1/4(-<ij|bd>) * (-t_ij^bd)
# ijbd,ibjd->ijb
# abcd,acbd->abd
e_oovv.contract("abcd,acbd->abd", t_2, factor=-0.25, out=tmp_ijb)
tmp_ijb.expand("acd->abcd", out=sigma)
#
# -1/4 <jk||ba> * (t_kj^ab - t_kj^ba) - <ij||db> * t_ij^db
#
def get_X_iAJB_sigma_block_cholesky(self, sigma: DenseFourIndex) -> None:
"""Compute the diagonal term involving the oovv block for sigma_iAJB.
This corresponds to the 'X_iAJB_sigma' special block contraction.
Args:
sigma (DenseFourIndex): The current value of the iAJB sigma diagonal.
"""
# Load integrals and amplitudes
# NOTE: we assume that all objects are present as
# we ensure their existence when `self` is instantiated
e_oovv = self._ham.get("eri_oovv")
t_2 = self.rcc_iodata.t_2
# -1/4 <jk||ba> * (t_kj^ab - t_kj^ba) -> -1/4(<jk|ba> - <jk|ab>) * (t_kj^ab - t_kj^ba)
# (1) -1/4<jk|ba> * t_kj^ab
# jkba,kajb->ajb
# abcd,bdac->dac
tmp_ajb = e_oovv.contract("abcd,bdac->dac", t_2, factor=-0.25)
# (2) -1/4<jk|ba> * (-t_kj^ba)
# jkba,kbja->ajb
# abcd,bcad->dac
e_oovv.contract("abcd,bcad->dac", t_2, factor=0.25, out=tmp_ajb)
# (3) -1/4(-<jk|ab>) * t_kj^ab
# jkab,kajb->ajb
# abcd,bcad->dac
e_oovv.contract("abcd,bcad->dac", t_2, factor=0.25, out=tmp_ajb)
# (4) -1/4(-<jk|ab>) * (-t_kj^ba)
# jkab,kbja->ajb
# abcd,bdac->dac
e_oovv.contract("abcd,bdac->dac", t_2, factor=-0.25, out=tmp_ajb)
# expand ajb->iajb
tmp_ajb.expand("bcd->abcd", out=sigma)
# - <ij||db> * t_ij^db -> -(<ij|db> - <ij|bd>) * t_ij^db
# (1) - <ij|db> * t_ij^db
# ijdb,idjb->ijb
# abcd,acbd->abd
tmp_ijb = e_oovv.contract("abcd,acbd->abd", t_2, factor=-1.0)
# (2) -(-<ij|bd>) * t_ij^db
# ijbd,idjb->ijb
# abcd,adbc->abc
e_oovv.contract("abcd,adbc->abc", t_2, out=tmp_ijb)
# expand ijb->iajb
tmp_ijb.expand("acd->abcd", out=sigma)
1.3.4. Example 3: Cholesky Block Implementation with Helper Functions (Protocols A & B)
This is an example of a Cholesky block (belonging to GeneralEffHamCholeskyBlocks). Because its execution depends on the active linear algebra factory, we must provide both the Dense (Protocol A) and the On-the-Fly Cholesky (Protocol B) implementations within the same file.
This example also demonstrates how to use helper functions (e.g., get_t1_t2_X_ajbc). The dynamic loader binds them automatically because their names contain the exact block label.
Indices:
a, j, b, c(virtual, occupied, virtual, virtual) translates to thevovvdirectory.File Name:
x_ajbc.py
# X_ajbc = <ab|cj> + <ab|cd> * t_j^d - <am|cj> * t_m^b - <bm|jc> * t_m^a
# - <am|cd> * t_j^d * t_m^b - <mb|cd> * t_j^d * t_m^a + <mk|jc> * t_k^a * t_m^b
# - f_mc * t_jm^ba + <am|cd> * (t_jm^bd - t_jm^db) + <am||cd> * t_jm^bd - <bm|dc> * t_jm^da
# + <mk|cj> * t_km^ba - <km|cd> * t_kj^ab * t_m^d - <km||cd> * t_kj^ab * t_m^d
# - <km|cd> * t_k^a * (t_jm^bd - t_jm ^db) - <km||cd> * t_k^a * t_jm^bd
# + <km|cd> * t_m^b * t_jk^da + <km|cd> * t_j^d * t_mk^ba + <km|cd> * t_k^a * t_j^d * t_m^b
def get_X_ajbc_block(self, force_update: bool = False) -> DenseFourIndex:
"""Get or compute the Hamiltonian X_ajbc block.
Args:
force_update (bool): If True, forces recomputation of X_ajbc even if cached.
Returns:
DenseFourIndex: The X_ajbc block.
"""
if "X_ajbc" in self._cache and not force_update:
return self._cache.load("X_ajbc")
# Load integrals and amplitudes
# NOTE: we assume that all objects are present as
# we ensure their existence when `self` is instantiated
fock_ov = self._ham.get("fock_ov")
e_ooov = self._ham.get("eri_ooov")
e_oovo = self._ham.get("eri_oovo")
e_oovv = self._ham.get("eri_oovv")
e_vovo = self._ham.get("eri_vovo")
e_voov = self._ham.get("eri_voov")
e_ovvv = self._ham.get("eri_ovvv")
e_vvvv = self._ham.get("eri_vvvv")
t_1 = self.rcc_iodata.t_1
t_2 = self.rcc_iodata.t_2
e_vovv = self._ham.get("eri_vovv")
# <ab|cj> -> <aj|cb>
# abcd->abdc
X_ajbc = self.init_block("X_ajbc", alloc=self.alloc(e_vovv, "abcd->abdc"))
# <ab|cd> * t_j^d
# abcd,jd -> ajbc
# abcd,ed -> aebc
e_vvvv.contract("abcd,ed->aebc", t_1, factor=1.0, out=X_ajbc)
# - <am|cj> * t_m^b
# amcj,mb -> ajbc
# abcd,be -> adec
e_vovo.contract("abcd,be->adec", t_1, factor=-1.0, out=X_ajbc)
# - <bm|jc> * t_m^a
# bmjc,ma -> ajbc
# abcd,be -> ecad
e_voov.contract("abcd,be->ecad", t_1, factor=-1.0, out=X_ajbc)
# - f_mc * t_jm^ba
# jbma,mc -> ajbc
# abcd,ce -> dabe
t_2.contract("abcd,ce->dabe", fock_ov, factor=-1.0, out=X_ajbc)
# <am|cd> * (t_jm^bd - t_jm^db)
# amcd,jbmd -> ajbc
# abcd,efbd -> aefc
e_vovv.contract("abcd,efbd->aefc", t_2, factor=2.0, out=X_ajbc)
# amcd,jdmb -> ajbc
# abcd,edbf -> aefc
e_vovv.contract("abcd,edbf->aefc", t_2, factor=-1.0, out=X_ajbc)
# <am||cd> * t_jm^bd -> (<am|cd> - <am|dc>) * t_jm^bd
# First part is done by factor=2.0 in line 776 as these contractions are identical
# amdc,jbmd -> ajbc
# abcd,efbc -> aefd
e_vovv.contract("abcd,efbc->aefd", t_2, factor=-1.0, out=X_ajbc)
# - <bm|dc> * t_jm^da
# bmdc,jdma -> ajbc
# abcd,ecbf -> fead
e_vovv.contract("abcd,ecbf->fead", t_2, factor=-1.0, out=X_ajbc)
# <mk|cj> * t_km^ba
# mkcj,kbma -> ajbc
# abcd,beaf -> fdec
e_oovo.contract("abcd,beaf->fdec", t_2, factor=1.0, out=X_ajbc)
if self.disconnected:
#
# a) get T1T2 part of xajbc intermediate
#
self.get_t1_t2_X_ajbc(X_ajbc, t_1, t_2)
# - <am|cd> * t_j^d * t_m^b
# amcd,jd -> ajmc > ajmc,mb -> ajbc
# abcd,ed -> aebc > abcd,ce -> abed
tmp = e_vovv.contract("abcd,ed->aebc", t_1, factor=-1.0)
tmp.contract("abcd,ce->abed", t_1, out=X_ajbc)
del tmp
# - <mb|cd> * t_j^d * t_m^a
# mbcd,jd -> mjbc > mjbc,ma -> ajbc
# abcd,ed -> aebc > abcd,ae -> ebcd
tmp = e_ovvv.contract("abcd,ed->aebc", t_1, factor=-1.0)
tmp.contract("abcd,ae->ebcd", t_1, out=X_ajbc)
del tmp
# <mk|jc> * t_k^a * t_m^b
# mkjc,ka -> ajmc > ajmc,mb -> ajbc
# abcd,be -> ecad > abcd,ce -> abed
tmp = e_ooov.contract("abcd,be->ecad", t_1, factor=1.0)
tmp.contract("abcd,ce->abed", t_1, out=X_ajbc)
del tmp
# <km|cd> * t_k^a * t_j^d * t_m^b
# kmcd,ka -> admc > admc,jd -> ajmc > ajmc,mb -> ajbc
tmp = e_oovv.contract("abcd,ae->edbc", t_1)
tmp = tmp.contract("abcd,eb->aecd", t_1)
tmp.contract("abcd,ce->abed", t_1, out=X_ajbc)
return X_ajbc
# - <km|cd> * t_kj^ab * t_m^d - <km||cd> * t_kj^ab * t_m^d
# - <km|cd> * t_k^a * (t_jm^bd - t_jm ^db) - <km||cd> * t_k^a * t_jm^bd
# + <km|cd> * t_m^b * t_jk^db + <km|cd> * t_j^d * t_mk^ba
def get_t1_t2_X_ajbc(
self,
X_ajbc: DenseFourIndex,
t_1: DenseTwoIndex,
t_2: DenseFourIndex,
) -> None:
"""Calculate the T_1 * T_2 contribution to the effective Hamiltonian term for the X_ajbc intermediate.
Args:
X_ajbc (DenseFourIndex): Output array representing the X_ajbc intermediate
t_1 (DenseTwoIndex): T_1 amplidutes
t_2 (DenseFourIndex): T_2 amplidutes
"""
# Load integrals and amplitudes
# NOTE: we assume that all objects are present as
# we ensure their existence when `self` is instantiated
e_oovv = self._ham.get("eri_oovv")
# - <km|cd> * t_kj^ab * t_m^d
# - <km||cd> * t_kj^ab * t_m^d
# > - ( 2 <km|cd> - <km|dc> ) * t_kj^ab * t_m^d
# (1) kmcd,md -> kc
# abcd,bd -> ac
# (2) kc,kajb -> ajbc > kajb,kc -> ajbc
# abcd,ae -> bcde
tmp = e_oovv.contract("abcd,bd->ac", t_1, factor=-2.0)
# <km|dc> * t_kj^ab * t_m^d
# (1) kmdc,md -> kc
# abcd,bc -> ad
# (2) kajb,kc -> ajbc
# abcd,ae -> bcde
e_oovv.contract("abcd,bc->ad", t_1, out=tmp)
t_2.contract("abcd,ae->bcde", tmp, out=X_ajbc)
del tmp
# - <km|cd> * t_k^a * (t_jm^bd - t_jm ^db)
# kmcd,ka -> amdc > amdc,jbmd -> ajbc
# amdc,jdmb -> ajbc
tmp = e_oovv.contract("abcd,ae->ebdc", t_1, factor=-1.0)
tmp.contract("abcd,efbc->aefd", t_2, out=X_ajbc, factor=2.0)
tmp.contract("abcd,ecbf->aefd", t_2, out=X_ajbc, factor=-1.0)
# - <km||cd> * t_k^a * t_jm^bd -> - (<km|cd> - <km|dc>) * t_k^a * t_jm^bd
# kmcd,ka -> amdc > amdc,jbmd -> ajbc
# abcd,ae->ebdc done 5 lines above, stored as tmp
# abcd,efbc->aefd done 5 lines above, thus factor=2.0
del tmp
# kmdc,ka -> amdc > amdc,jbmd -> ajbc
tmp = e_oovv.contract("abcd,ae->ebcd", t_1, factor=1.0)
tmp.contract("abcd,efbc->aefd", t_2, out=X_ajbc)
del tmp
# <km|cd> * t_m^b * t_jk^da
# (1) kmcd,mb -> bkdc
# abcd,be -> eadc
# (2) bkdc,jdka -> ajbc
# abcd,ecbf -> fead
tmp = e_oovv.contract("abcd,be->eadc", t_1)
tmp.contract("abcd,ecbf->fead", t_2, out=X_ajbc)
# <km|cd> * t_j^d * t_mk^ba
# kmcd,jd -> kjmc > kjmc,mbka -> ajbc
tmp = e_oovv.contract("abcd,ed->aebc", t_1)
tmp.contract("abcd,ceaf->fbed", t_2, out=X_ajbc)
del tmp
# ajbc Cholesky case
# (<ab|cj> + <ab|cd> * t_j^d - <am|cj> * t_m^b - <bm|jc> * t_m^a
# - <am|cd> * t_j^d * t_m^b - <mb|cd> * t_j^d * t_m^a + <mk|jc> * t_k^a * t_m^b
# - f_mc * t_jm^ba + <am|cd> * (t_jm^bd - t_jm^db) + <am||cd> * t_jm^bd - <bm|dc> * t_jm^da
# + <mk|cj> * t_km^ba - <km|cd> * t_kj^ab * t_m^d - <km||cd> * t_kj^ab * t_m^d
# - <km|cd> * t_k^a * (t_jm^bd - t_jm ^db) - <km||cd> * t_k^a * t_jm^bd
# + <km|cd> * t_m^b * t_jk^da + <km|cd> * t_j^d * t_mk^ba + <km|cd> * t_k^a * t_j^d * t_m^b)
# * r_i^c
def get_X_ajbc_block_cholesky(
self, bvec: DenseTwoIndex, sigma: DenseFourIndex
) -> None:
"""Compute the Hamiltonian X_ajbc block using Cholesky representation.
Args:
bvec (DenseTwoIndex): Current approximation to the CI doubles coefficients
sigma (DenseFourIndex): Output sigma vector
"""
# Load integrals and amplitudes
# NOTE: we assume that all objects are present as
# we ensure their existence when `self` is instantiated
e_nnvo = self._ham.get("eri_nnvo")
e_nnvv = self._ham.get("eri_nnvv")
e_vovv = self._ham.get("eri_vovv")
oooo = self.get_range("oooo")
ovoo = self.get_range("ovoo")
ooov = self.get_range("ooov")
ovov = self.get_range("ovov")
f_ov = self._ham.get("fock_ov")
t_1 = self.rcc_iodata.t_1
t_2 = self.rcc_iodata.t_2
# <pq|cj> * r_i^c
# pqcj,ic->ipjq
# abcd,ec->eadb
tmp_ipjq_cj = e_nnvo.contract("abcd,ec->eadb", bvec)
# <ab|cj> * r_i^c
# tmp_ipjq -> tmp_iajb
tmp_ipjq_cj.contract("abcd->abcd", out=sigma, **ovov)
# (- <am|cj> * r_i^c) * t_m^b
# tmp_ipjq -> tmp_iajm
# iajm,mb -> iajb
# abcd,de -> abce
tmp_ipjq_cj.contract("abcd,de->abce", t_1, out=sigma, factor=-1.0, **ovoo)
# (- <bm|jc> * r_i^c) * t_m^a => (- <mb|cj> * r_i^c) * t_m^a
# tmp_ipjq -> tmp_imjb
# imjb,ma -> iajb
# abcd,be -> aebc
tmp_ipjq_cj.contract("abcd,be->aebc", t_1, out=sigma, factor=-1.0, **ooov)
# (<mk|cj> * r_i^c) * t_km^ba
# tmp_ipjq -> tmp_imjk
# imjk,kbma -> iajb
# abcd,debf -> afce
tmp_ipjq_cj.contract("abcd,debf->afce", t_2, out=sigma, **oooo)
# (- f_mc * t_jm^ba) * r_i^c
# mc,ic -> im
# ab,cb -> ca
tmp = f_ov.contract("ab,cb->ca", bvec)
# jbma,im -> iajb
# abcd,ec -> edab
t_2.contract("abcd,ec->edab", tmp, out=sigma, factor=-1.0)
del tmp
# <pq|cd> * r_i^c * t_j^d
# pqcd,ic,jd->ipjq
# abcd,ec,fd->eafb
tmp_ipjq_cd = e_nnvv.contract("abcd,ec,fd->eafb", t_1, bvec)
# (<ab|cd> * r_i^c) * t_j^d
# tmp_ipjq -> tmp_iajb
tmp_ipjq_cd.contract("abcd->abcd", out=sigma, **ovov)
# <am|cd> * (t_jm^bd - t_jm^db) * r_i^c
# amcd,ic -> iadm
# abcd,ec -> eadb
tmp = e_vovv.contract("abcd,ec->eadb", bvec)
# iadm,jdmb -> iajb
# abcd,ecdf -> abef
tmp.contract("abcd,ecdf->abef", t_2, out=sigma, factor=-1.0)
# iadm,jbmd -> iajb
# abcd,efdc -> abef
tmp.contract("abcd,efdc->abef", t_2, out=sigma, factor=2.0)
del tmp
# <am||cd> * t_jm^bd * r_i^c-> (<am|cd> - <am|dc>) * t_jm^bd * r_i^c
# iadm,jbmd -> iajb
# First part is done by factor=2.0 in the contraction above as these contractions are identical
# Exchange part
# amdc,ic -> iadm
# abcd,ed -> eacb
tmp = e_vovv.contract("abcd,ed->eacb", bvec)
# iadm,jbmd -> iajb
# abcd,efdc -> abef
tmp.contract("abcd,efdc->abef", t_2, out=sigma, factor=-1.0)
# - <bm|dc> * t_jm^da * r_i^c
# bmdc,ic -> ibdm
# abcd,ed -> eacb
# Same as exchange part from above
# ibdm,jdma -> iajb
# abcd,ebdf -> afeb
tmp.contract("abcd,ebdf->afeb", t_2, out=sigma, factor=-1.0)
del tmp
if self.disconnected:
# (<mk|jc> * r_i^c) * t_k^a * t_m^b => (<km|cj> * r_i^c) * t_k^a * t_m^b
# tmp_ipjq -> tmp_ikjm
# ikjm,ka -> iajm
# abcd,be -> aecd
tmp = tmp_ipjq_cj.contract("abcd,be->aecd", t_1, **oooo)
# iajm,mb -> iajb
# abcd,ce -> abce
tmp.contract("abcd,ce->abce", t_1, out=sigma)
del tmp
del tmp_ipjq_cj
# (- <am|cd> * t_j^d * t_m^b) * r_i^c
# tmp_ipjq -> tmp_iajm
# iajm,mb -> iajb
# abcd,de -> abce
tmp_ipjq_cd.contract(
"abcd,de->abce", t_1, out=sigma, factor=-1.0, **ovoo
)
# (- <mb|cd> * t_j^d * t_m^a) * r_i^c
# tmp_ipjq -> tmp_imjb
# imjb,ma -> iajb
# abcd,be -> aecd
tmp_ipjq_cd.contract(
"abcd,be->aecd", t_1, out=sigma, factor=-1.0, **ooov
)
# (<km|cd> * t_k^a * t_j^d * t_m^b) * r_i^c
# tmp_ipjq -> tmp_ikjm
# ikjm,ka -> iajm
# abcd,be -> aecd
tmp = tmp_ipjq_cd.contract("abcd,be->aecd", t_1, **oooo)
# iajm,mb -> iajb
# abcd,de -> abce
tmp.contract("abcd,de->abce", t_1, out=sigma)
del tmp
del tmp_ipjq_cd
# a) get T1T2 part of xajbc intermediate
self.get_t1_t2_X_ajbc_cholesky(sigma, bvec, t_1, t_2)
else:
# If disconnected is False, we delete the temporary objects
del tmp_ipjq_cj
del tmp_ipjq_cd
# (- <km|cd> * t_kj^ab * t_m^d - <km||cd> * t_kj^ab * t_m^d
# - <km|cd> * t_k^a * (t_jm^bd - t_jm ^db) - <km||cd> * t_k^a * t_jm^bd
# + <km|cd> * t_m^b * t_jk^db + <km|cd> * t_j^d * t_mk^ba)
# * r_i^c
def get_t1_t2_X_ajbc_cholesky(
self,
sigma: DenseFourIndex,
bvec: DenseTwoIndex,
t_1: DenseTwoIndex,
t_2: DenseFourIndex,
) -> None:
"""Calculate the T_1 * T_2 contribution to the effective Hamiltonian term for the X_ajbc intermediate.
**Arguments:**
:sigma: (DenseFourIndex)
The output sigma vector
:bvec: (DenseTwoIndex)
The current approximation to the CI doubles coefficient
:t_1: (DenseTwoIndex)
T_1 amplidutes
:t_2: (DenseFourIndex)
T_2 amplidutes
"""
# Load integrals and amplitudes
# NOTE: we assume that all objects are present as
# we ensure their existence when `self` is instantiated
e_oovv = self._ham.get("eri_oovv")
# <km|cd> * r_i^c
# kmcd,ic->ikdm
# abcd,ec->eadb
tmp_ikdm = e_oovv.contract("abcd,ec->eadb", bvec)
# (- <km|cd> * t_kj^ab * t_m^d) * r_i^c
# (- <km||cd> * t_kj^ab * t_m^d) * r_i^c
# > r_i^c * [- ( 2 <km|cd> - <km|dc> ) * t_kj^ab * t_m^d]
# kmdc,ic -> ikdm
# abcd,ed -> eacb
tmp_exchange = e_oovv.contract("abcd,ed->eacb", bvec)
# ikdm,md -> ik
# abcd,dc -> ab
tmp = tmp_exchange.contract("abcd,dc->ab", t_1)
# ikdm,md -> ik
# abcd,dc -> ab
tmp_ikdm.contract("abcd,dc->ab", t_1, out=tmp, factor=-2.0)
# kajb,ik -> iajb
# abcd,ea -> eabc
t_2.contract("abcd,ea->eabc", tmp, out=sigma)
del tmp
# [- <km|cd> * t_k^a * (t_jm^bd - t_jm ^db)] * r_i^c
# (- <km||cd> * t_k^a * t_jm^bd) * r_i^c -> - [(<km|cd> - <km|dc>) * t_k^a * t_jm^bd] * r_i^c
# > r_i^c * [- ( 2<km|cd> - <km|dc> ) * t_k^a * t_jm^bd]
# kmdc,ic -> ikdm
# abcd,ed -> eacb
# Same as tmp_exchange above
# ikdm,ka -> iadm
# abcd,be -> aebc
tmp_exchange = tmp_exchange.contract("abcd,be->aebc", t_1)
# ikdm,ka -> iadm
# abcd,be -> aebc
tmp_ikdm.contract("abcd,be->aebc", t_1, out=tmp_exchange, factor=-2.0)
# iadm,jbmd -> iajb
# abcd,efdc -> abef
tmp_exchange.contract("abcd,efdc->abef", t_2, out=sigma)
del tmp_exchange
# r_i^c * [- <km|cd> * t_k^a * (-t_jm^db)]
# ikdm,ka -> iadm
# abcd,be -> aebc
tmp = tmp_ikdm.contract("abcd,be->aebc", t_1)
# iadm,jdmb -> iajb
# abcd,ecdf -> abef
tmp.contract("abcd,ecdf->abef", t_2, out=sigma)
del tmp
# (<km|cd> * t_m^b * t_jk^da) * r_i^c
# ikdm,mb -> ikdb
# abcd,de -> abce
tmp = tmp_ikdm.contract("abcd,de->abce", t_1)
# ikdb,jdka -> iajb
# abcd,ecbf -> afed
tmp.contract("abcd,ecbf->afed", t_2, out=sigma)
del tmp
# (<km|cd> * t_j^d * t_mk^ba) * r_i^c
# ikdm,jd -> ikjm
# abcd,ec -> abed
tmp = tmp_ikdm.contract("abcd,ec->abed", t_1)
# ikjm,mbka -> iajb
# abcd,debf -> afce
tmp.contract("abcd,debf->afce", t_2, out=sigma)
del tmp
del tmp_ikdm