Source code for pyscf.pbc.cc.kintermediates_rhf

#!/usr/bin/env python
# Copyright 2014-2020 The PySCF Developers. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Authors: James D. McClain
#          Timothy Berkelbach <tim.berkelbach@gmail.com>
#

import numpy as np

from itertools import product
from pyscf import lib
from pyscf.lib import logger
from pyscf.lib.parameters import LOOSE_ZERO_TOL, LARGE_DENOM  # noqa
from pyscf.pbc.lib import kpts_helper
from pyscf.pbc.mp.kmp2 import (get_frozen_mask, get_nocc, get_nmo,
                               padded_mo_coeff, padding_k_idx)  # noqa

#einsum = np.einsum
einsum = lib.einsum

# This is restricted (R)CCSD
# Ref: Hirata et al., J. Chem. Phys. 120, 2581 (2004)

### Eqs. (37)-(39) "kappa"

[docs] def cc_Foo(t1,t2,eris,kconserv): nkpts, nocc, nvir = t1.shape Fki = np.empty((nkpts,nocc,nocc),dtype=t2.dtype) for ki in range(nkpts): kk = ki Fki[ki] = eris.fock[ki,:nocc,:nocc].copy() for kl in range(nkpts): for kc in range(nkpts): kd = kconserv[kk,kc,kl] Soovv = 2*eris.oovv[kk,kl,kc] - eris.oovv[kk,kl,kd].transpose(0,1,3,2) Fki[ki] += einsum('klcd,ilcd->ki',Soovv,t2[ki,kl,kc]) #if ki == kc: kd = kconserv[kk,ki,kl] Soovv = 2*eris.oovv[kk,kl,ki] - eris.oovv[kk,kl,kd].transpose(0,1,3,2) Fki[ki] += einsum('klcd,ic,ld->ki',Soovv,t1[ki],t1[kl]) return Fki
[docs] def cc_Fvv(t1,t2,eris,kconserv): nkpts, nocc, nvir = t1.shape Fac = np.empty((nkpts,nvir,nvir),dtype=t2.dtype) for ka in range(nkpts): kc = ka Fac[ka] = eris.fock[ka,nocc:,nocc:].copy() for kl in range(nkpts): for kk in range(nkpts): kd = kconserv[kk,kc,kl] Soovv = 2*eris.oovv[kk,kl,kc] - eris.oovv[kk,kl,kd].transpose(0,1,3,2) Fac[ka] += -einsum('klcd,klad->ac',Soovv,t2[kk,kl,ka]) #if kk == ka kd = kconserv[ka,kc,kl] Soovv = 2*eris.oovv[ka,kl,kc] - eris.oovv[ka,kl,kd].transpose(0,1,3,2) Fac[ka] += -einsum('klcd,ka,ld->ac',Soovv,t1[ka],t1[kl]) return Fac
[docs] def cc_Fov(t1,t2,eris,kconserv): nkpts, nocc, nvir = t1.shape Fkc = np.empty((nkpts,nocc,nvir),dtype=t2.dtype) Fkc[:] = eris.fock[:,:nocc,nocc:].copy() for kk in range(nkpts): for kl in range(nkpts): Soovv = 2.*eris.oovv[kk,kl,kk] - eris.oovv[kk,kl,kl].transpose(0,1,3,2) Fkc[kk] += einsum('klcd,ld->kc',Soovv,t1[kl]) return Fkc
### Eqs. (40)-(41) "lambda"
[docs] def Loo(t1,t2,eris,kconserv): nkpts, nocc, nvir = t1.shape fov = eris.fock[:,:nocc,nocc:] Lki = cc_Foo(t1,t2,eris,kconserv) for ki in range(nkpts): Lki[ki] += einsum('kc,ic->ki',fov[ki],t1[ki]) for kl in range(nkpts): Lki[ki] += 2*einsum('klic,lc->ki',eris.ooov[ki,kl,ki],t1[kl]) Lki[ki] += -einsum('lkic,lc->ki',eris.ooov[kl,ki,ki],t1[kl]) return Lki
[docs] def Lvv(t1,t2,eris,kconserv): nkpts, nocc, nvir = t1.shape fov = eris.fock[:,:nocc,nocc:] Lac = cc_Fvv(t1,t2,eris,kconserv) for ka in range(nkpts): Lac[ka] += -einsum('kc,ka->ac',fov[ka],t1[ka]) for kk in range(nkpts): Svovv = 2*eris.vovv[ka,kk,ka] - eris.vovv[ka,kk,kk].transpose(0,1,3,2) Lac[ka] += einsum('akcd,kd->ac',Svovv,t1[kk]) return Lac
### Eqs. (42)-(45) "chi"
[docs] def cc_Woooo(t1, t2, eris, kconserv, out=None): nkpts, nocc, nvir = t1.shape Wklij = _new(eris.oooo.shape, t1.dtype, out) for kk in range(nkpts): for kl in range(kk+1): for ki in range(nkpts): kj = kconserv[kk,ki,kl] oooo = einsum('klic,jc->klij',eris.ooov[kk,kl,ki],t1[kj]) oooo += einsum('lkjc,ic->klij',eris.ooov[kl,kk,kj],t1[ki]) oooo += eris.oooo[kk,kl,ki] # ==== Beginning of change ==== # #for kc in range(nkpts): # Wklij[kk,kl,ki] += einsum('klcd,ijcd->klij',eris.oovv[kk,kl,kc],t2[ki,kj,kc]) #Wklij[kk,kl,ki] += einsum('klcd,ic,jd->klij',eris.oovv[kk,kl,ki],t1[ki],t1[kj]) vvoo = eris.oovv[kk,kl].transpose(0,3,4,1,2).reshape(nkpts*nvir,nvir,nocc,nocc) t2t = t2[ki,kj].copy().transpose(0,3,4,1,2) #for kc in range(nkpts): # kd = kconserv[ki,kc,kj] # if kc == ki and kj == kd: # t2t[kc] += einsum('ic,jd->cdij',t1[ki],t1[kj]) t2t[ki] += einsum('ic,jd->cdij',t1[ki],t1[kj]) t2t = t2t.reshape(nkpts*nvir,nvir,nocc,nocc) oooo += einsum('cdkl,cdij->klij',vvoo,t2t) Wklij[kk,kl,ki] = oooo # ===== End of change = ==== # Be careful about making this term only after all the others are created for kl in range(kk+1): for ki in range(nkpts): kj = kconserv[kk,ki,kl] Wklij[kl,kk,kj] = Wklij[kk,kl,ki].transpose(1,0,3,2) return Wklij
[docs] def cc_Wvvvv(t1, t2, eris, kconserv, out=None): Wabcd = _new(eris.vvvv.shape, t1.dtype, out) nkpts, nocc, nvir = t1.shape for ka in range(nkpts): for kb in range(ka+1): for kc in range(nkpts): kd = kconserv[ka,kc,kb] # avoid transpose in loop vvvv = einsum('akcd,kb->abcd', eris.vovv[ka,kb,kc], -t1[kb]) vvvv += einsum('bkdc,ka->abcd', eris.vovv[kb,ka,kd], -t1[ka]) vvvv += eris.vvvv[ka,kb,kc] Wabcd[ka,kb,kc] = vvvv # Be careful: making this term only after all the others are created for kb in range(ka+1): for kc in range(nkpts): kd = kconserv[ka,kc,kb] Wabcd[kb,ka,kd] = Wabcd[ka,kb,kc].transpose(1,0,3,2) return Wabcd
[docs] def cc_Wvoov(t1, t2, eris, kconserv, out=None): Wakic = _new(eris.voov.shape, t1.dtype, out) nkpts, nocc, nvir = t1.shape for ka in range(nkpts): for kk in range(nkpts): voov_i = einsum('xakdc,xid->xakic',eris.vovv[ka,kk,:],t1[:]) voov_i -= einsum('xlkic,la->xakic',eris.ooov[ka,kk,:],t1[ka]) voov_i += eris.voov[ka,kk,:] for ki in range(nkpts): kc = kconserv[ka,ki,kk] #for kl in range(nkpts): # # kl - kd + kk = kc # # => kd = kl - kc + kk # kd = kconserv[kl,kc,kk] # Soovv = 2*eris.oovv[kl,kk,kd] - eris.oovv[kl,kk,kc].transpose(0,1,3,2) # Wakic[ka,kk,ki] += 0.5*einsum('lkdc,ilad->akic',Soovv,t2[ki,kl,ka]) # Wakic[ka,kk,ki] -= 0.5*einsum('lkdc,ilda->akic',eris.oovv[kl,kk,kd],t2[ki,kl,kd]) #Wakic[ka,kk,ki] -= einsum('lkdc,id,la->akic',eris.oovv[ka,kk,ki],t1[ki],t1[ka]) kd = kconserv[ka,kc,kk] tau = t2[:,ki,ka].copy() tau[ka] += 2*einsum('id,la->liad',t1[kd],t1[ka]) oovv_tmp = np.array(eris.oovv[kk,:,kc]) voov_i[ki] -= 0.5*einsum('xklcd,xliad->akic',oovv_tmp,tau) Soovv_tmp = 2*oovv_tmp - eris.oovv[:,kk,kc].transpose(0,2,1,3,4) voov_i[ki] += 0.5*einsum('xklcd,xilad->akic',Soovv_tmp,t2[ki,:,ka]) Wakic[ka,kk,:] = voov_i[:] return Wakic
[docs] def cc_Wvovo(t1, t2, eris, kconserv, out=None): nkpts, nocc, nvir = t1.shape Wakci = _new((nkpts,nkpts,nkpts,nvir,nocc,nvir,nocc), t1.dtype, out) for ka in range(nkpts): for kk in range(nkpts): for kc in range(nkpts): ki = kconserv[ka,kc,kk] vovo = einsum('akcd,id->akci',eris.vovv[ka,kk,kc],t1[ki]) vovo -= einsum('klic,la->akci',eris.ooov[kk,ka,ki],t1[ka]) vovo += np.asarray(eris.ovov[kk,ka,ki]).transpose(1,0,3,2) # ==== Beginning of change ==== # #for kl in range(nkpts): # kd = kconserv[kl,kc,kk] # Wakci[ka,kk,kc] -= 0.5*einsum('lkcd,ilda->akci',eris.oovv[kl,kk,kc],t2[ki,kl,kd]) #Wakci[ka,kk,kc] -= einsum('lkcd,id,la->akci',eris.oovv[ka,kk,kc],t1[ki],t1[ka]) oovvf = eris.oovv[:,kk,kc].reshape(nkpts*nocc,nocc,nvir,nvir) t2f = t2[:,ki,ka].copy() #This is a tau like term #for kl in range(nkpts): # kd = kconserv[kl,kc,kk] # if ki == kd and kl == ka: # t2f[kl] += 2*einsum('id,la->liad',t1[ki],t1[ka]) kd = kconserv[ka,kc,kk] t2f[ka] += 2*einsum('id,la->liad',t1[kd],t1[ka]) t2f = t2f.reshape(nkpts*nocc,nocc,nvir,nvir) vovo -= 0.5*einsum('lkcd,liad->akci',oovvf,t2f) Wakci[ka,kk,kc] = vovo # ===== End of change = ==== return Wakci
[docs] def Wooov(t1, t2, eris, kconserv, out=None): nkpts, nocc, nvir = t1.shape Wklid = _new(eris.ooov.shape, t1.dtype, out) for kk in range(nkpts): for kl in range(nkpts): for ki in range(nkpts): ooov = einsum('ic,klcd->klid',t1[ki],eris.oovv[kk,kl,ki]) ooov += eris.ooov[kk,kl,ki] Wklid[kk,kl,ki] = ooov return Wklid
[docs] def Wvovv(t1, t2, eris, kconserv, out=None): nkpts, nocc, nvir = t1.shape Walcd = _new(eris.vovv.shape, t1.dtype, out) for ka in range(nkpts): for kl in range(nkpts): for kc in range(nkpts): vovv = einsum('ka,klcd->alcd', -t1[ka], eris.oovv[ka,kl,kc]) vovv += eris.vovv[ka,kl,kc] Walcd[ka,kl,kc] = vovv return Walcd
[docs] def W1ovvo(t1, t2, eris, kconserv, out=None): nkpts, nocc, nvir = t1.shape Wkaci = _new((nkpts,nkpts,nkpts,nocc,nvir,nvir,nocc), t1.dtype, out) for kk in range(nkpts): for ka in range(nkpts): for kc in range(nkpts): ki = kconserv[kk,kc,ka] # ovvo[kk,ka,kc,ki] => voov[ka,kk,ki,kc] ovvo = np.asarray(eris.voov[ka,kk,ki]).transpose(1,0,3,2).copy() for kl in range(nkpts): kd = kconserv[ki,ka,kl] St2 = 2.*t2[ki,kl,ka] - t2[kl,ki,ka].transpose(1,0,2,3) ovvo += einsum('klcd,ilad->kaci',eris.oovv[kk,kl,kc],St2) ovvo += -einsum('kldc,ilad->kaci',eris.oovv[kk,kl,kd],t2[ki,kl,ka]) Wkaci[kk,ka,kc] = ovvo return Wkaci
[docs] def W2ovvo(t1, t2, eris, kconserv, out=None): nkpts, nocc, nvir = t1.shape Wkaci = _new((nkpts,nkpts,nkpts,nocc,nvir,nvir,nocc), t1.dtype, out) WWooov = Wooov(t1,t2,eris,kconserv) for kk in range(nkpts): for ka in range(nkpts): for kc in range(nkpts): ki = kconserv[kk,kc,ka] ovvo = einsum('la,lkic->kaci',-t1[ka],WWooov[ka,kk,ki]) ovvo += einsum('akdc,id->kaci',eris.vovv[ka,kk,ki],t1[ki]) Wkaci[kk,ka,kc] = ovvo return Wkaci
[docs] def Wovvo(t1, t2, eris, kconserv, out=None): Wovvo = W1ovvo(t1, t2, eris, kconserv, out) for k, w2 in enumerate(W2ovvo(t1, t2, eris, kconserv)): Wovvo[k] = Wovvo[k] + w2 return Wovvo
[docs] def W1ovov(t1, t2, eris, kconserv, out=None): nkpts, nocc, nvir = t1.shape Wkbid = _new(eris.ovov.shape, t1.dtype, out) for kk in range(nkpts): for kb in range(nkpts): for ki in range(nkpts): kd = kconserv[kk,ki,kb] # kk + kl - kc - kd = 0 # => kc = kk - kd + kl ovov = eris.ovov[kk,kb,ki].copy() for kl in range(nkpts): kc = kconserv[kk,kd,kl] ovov -= einsum('klcd,ilcb->kbid',eris.oovv[kk,kl,kc],t2[ki,kl,kc]) Wkbid[kk,kb,ki] = ovov return Wkbid
[docs] def W2ovov(t1, t2, eris, kconserv, out=None): nkpts, nocc, nvir = t1.shape Wkbid = _new((nkpts,nkpts,nkpts,nocc,nvir,nocc,nvir), t1.dtype, out) WWooov = Wooov(t1,t2,eris,kconserv) for kk in range(nkpts): for kb in range(nkpts): for ki in range(nkpts): kd = kconserv[kk,ki,kb] ovov = einsum('klid,lb->kbid',WWooov[kk,kb,ki],-t1[kb]) ovov += einsum('bkdc,ic->kbid',eris.vovv[kb,kk,kd],t1[ki]) Wkbid[kk,kb,ki] = ovov return Wkbid
[docs] def Wovov(t1, t2, eris, kconserv, out=None): Wovov = W1ovov(t1, t2, eris, kconserv, out) for k, w2 in enumerate(W2ovov(t1, t2, eris, kconserv)): Wovov[k] = Wovov[k] + w2 return Wovov
[docs] def Woooo(t1, t2, eris, kconserv, out=None): nkpts, nocc, nvir = t1.shape Wklij = _new(eris.oooo.shape, t1.dtype, out) for kk in range(nkpts): for kl in range(nkpts): for ki in range(nkpts): kj = kconserv[kk,ki,kl] oooo = einsum('klcd,ic,jd->klij',eris.oovv[kk,kl,ki],t1[ki],t1[kj]) oooo += einsum('klid,jd->klij',eris.ooov[kk,kl,ki],t1[kj]) oooo += einsum('lkjc,ic->klij',eris.ooov[kl,kk,kj],t1[ki]) oooo += eris.oooo[kk,kl,ki] for kc in range(nkpts): #kd = kconserv[kk,kc,kl] oooo += einsum('klcd,ijcd->klij',eris.oovv[kk,kl,kc],t2[ki,kj,kc]) Wklij[kk,kl,ki] = oooo return Wklij
[docs] def Wvvvv(t1, t2, eris, kconserv, out=None): nkpts, nocc, nvir = t1.shape Wabcd = _new((nkpts,nkpts,nkpts,nvir,nvir,nvir,nvir), t2.dtype, out) for ka in range(nkpts): for kb in range(nkpts): for kc in range(nkpts): Wabcd[ka,kb,kc] = get_Wvvvv(t1, t2, eris, kconserv, ka, kb, kc) return Wabcd
[docs] def get_Wvvvv(t1, t2, eris, kconserv, ka, kb, kc): kd = kconserv[ka, kc, kb] nkpts, nocc, nvir = t1.shape if getattr(eris, 'Lpv', None) is not None: # Using GDF to generate Wvvvv on the fly Lpv = eris.Lpv Lac = (Lpv[ka,kc][:,nocc:] - einsum('Lkc,ka->Lac', Lpv[ka,kc][:,:nocc], t1[ka])) Lbd = (Lpv[kb,kd][:,nocc:] - einsum('Lkd,kb->Lbd', Lpv[kb,kd][:,:nocc], t1[kb])) vvvv = einsum('Lac,Lbd->abcd', Lac, Lbd) vvvv *= (1. / nkpts) else: vvvv = einsum('klcd,ka,lb->abcd',eris.oovv[ka,kb,kc],t1[ka],t1[kb]) vvvv += einsum('alcd,lb->abcd',eris.vovv[ka,kb,kc],-t1[kb]) vvvv += einsum('bkdc,ka->abcd',eris.vovv[kb,ka,kd],-t1[ka]) vvvv += eris.vvvv[ka,kb,kc] for kk in range(nkpts): kl = kconserv[kc,kk,kd] vvvv += einsum('klcd,klab->abcd', eris.oovv[kk,kl,kc], t2[kk,kl,ka]) return vvvv
[docs] def Wvvvo(t1, t2, eris, kconserv, _Wvvvv=None, out=None): nkpts, nocc, nvir = t1.shape Wabcj = _new((nkpts,nkpts,nkpts,nvir,nvir,nvir,nocc), t1.dtype, out) WW1ovov = W1ovov(t1,t2,eris,kconserv) WW1ovvo = W1ovvo(t1,t2,eris,kconserv) FFov = cc_Fov(t1,t2,eris,kconserv) for ka in range(nkpts): for kb in range(nkpts): for kc in range(nkpts): kj = kconserv[ka,kc,kb] # Wvovo[ka,kl,kc,kj] <= Wovov[kl,ka,kj,kc].transpose(1,0,3,2) vvvo = einsum('alcj,lb->abcj',WW1ovov[kb,ka,kj].transpose(1,0,3,2),-t1[kb]) vvvo += einsum('kbcj,ka->abcj',WW1ovvo[ka,kb,kc],-t1[ka]) # vvvo[ka,kb,kc,kj] <= vovv[kc,kj,ka,kb].transpose(2,3,0,1).conj() vvvo += np.asarray(eris.vovv[kc,kj,ka]).transpose(2,3,0,1).conj() for kl in range(nkpts): # ka + kl - kc - kd = 0 # => kd = ka - kc + kl kd = kconserv[ka,kc,kl] St2 = 2.*t2[kl,kj,kd] - t2[kl,kj,kb].transpose(0,1,3,2) vvvo += einsum('alcd,ljdb->abcj',eris.vovv[ka,kl,kc], St2) vvvo += einsum('aldc,ljdb->abcj',eris.vovv[ka,kl,kd], -t2[kl,kj,kd]) # kb - kc + kl = kd kd = kconserv[kb,kc,kl] vvvo += einsum('bldc,jlda->abcj',eris.vovv[kb,kl,kd], -t2[kj,kl,kd]) # kl + kk - kb - ka = 0 # => kk = kb + ka - kl kk = kconserv[kb,kl,ka] vvvo += einsum('lkjc,lkba->abcj',eris.ooov[kl,kk,kj],t2[kl,kk,kb]) vvvo += einsum('lkjc,lb,ka->abcj',eris.ooov[kb,ka,kj],t1[kb],t1[ka]) vvvo += einsum('lc,ljab->abcj',-FFov[kc],t2[kc,kj,ka]) Wabcj[ka,kb,kc] = vvvo # Check if t1=0 (HF+MBPT(2)) # einsum will check, but don't make vvvv if you can avoid it! if np.any(t1 != 0): for ka in range(nkpts): for kb in range(nkpts): for kc in range(nkpts): kj = kconserv[ka,kc,kb] if _Wvvvv is None: Wvvvv = get_Wvvvv(t1, t2, eris, kconserv, ka, kb, kc) else: Wvvvv = _Wvvvv[ka, kb, kc] Wabcj[ka,kb,kc] = (Wabcj[ka,kb,kc] + einsum('abcd,jd->abcj', Wvvvv, t1[kj])) return Wabcj
[docs] def Wovoo(t1, t2, eris, kconserv, out=None): nkpts, nocc, nvir = t1.shape WW1ovov = W1ovov(t1,t2,eris,kconserv) WWoooo = Woooo(t1,t2,eris,kconserv) WW1ovvo = W1ovvo(t1,t2,eris,kconserv) FFov = cc_Fov(t1,t2,eris,kconserv) Wkbij = _new((nkpts,nkpts,nkpts,nocc,nvir,nocc,nocc), t1.dtype, out) for kk in range(nkpts): for kb in range(nkpts): for ki in range(nkpts): kj = kconserv[kk,ki,kb] ovoo = einsum('kbid,jd->kbij',WW1ovov[kk,kb,ki], t1[kj]) ovoo += einsum('klij,lb->kbij',WWoooo[kk,kb,ki],-t1[kb]) ovoo += einsum('kbcj,ic->kbij',WW1ovvo[kk,kb,ki],t1[ki]) ovoo += np.array(eris.ooov[ki,kj,kk]).transpose(2,3,0,1).conj() for kd in range(nkpts): # kk + kl - ki - kd = 0 # => kl = ki - kk + kd kl = kconserv[ki,kk,kd] St2 = 2.*t2[kl,kj,kd] - t2[kj,kl,kd].transpose(1,0,2,3) ovoo += einsum('klid,ljdb->kbij', eris.ooov[kk,kl,ki], St2) ovoo += einsum('lkid,ljdb->kbij', -eris.ooov[kl,kk,ki],t2[kl,kj,kd]) kl = kconserv[kb,ki,kd] ovoo += einsum('lkjd,libd->kbij', -eris.ooov[kl,kk,kj],t2[kl,ki,kb]) # kb + kk - kd = kc #kc = kconserv[kb,kd,kk] ovoo += einsum('bkdc,jidc->kbij',eris.vovv[kb,kk,kd],t2[kj,ki,kd]) ovoo += einsum('bkdc,jd,ic->kbij',eris.vovv[kb,kk,kj],t1[kj],t1[ki]) ovoo += einsum('kc,ijcb->kbij',FFov[kk],t2[ki,kj,kk]) Wkbij[kk,kb,ki] = ovoo return Wkbij
def _new(shape, dtype, out): if out is None: # Incore: out = np.empty(shape, dtype=dtype) else: assert (out.shape == shape) assert (out.dtype == dtype) return out
[docs] def get_t3p2_imds_slow(cc, t1, t2, eris=None, t3p2_ip_out=None, t3p2_ea_out=None): """For a description of arguments, see `get_t3p2_imds_slow` in the corresponding `kintermediates.py`. """ from pyscf.pbc.cc.kccsd_t_rhf import _get_epqr if eris is None: eris = cc.ao2mo() fock = eris.fock nkpts, nocc, nvir = t1.shape kconserv = cc.khelper.kconserv dtype = np.result_type(t1, t2) fov = fock[:, :nocc, nocc:] #foo = [fock[ikpt, :nocc, :nocc].diagonal() for ikpt in range(nkpts)] #fvv = [fock[ikpt, nocc:, nocc:].diagonal() for ikpt in range(nkpts)] mo_energy_occ = np.array([eris.mo_energy[ki][:nocc] for ki in range(nkpts)]) mo_energy_vir = np.array([eris.mo_energy[ki][nocc:] for ki in range(nkpts)]) mo_e_o = mo_energy_occ mo_e_v = mo_energy_vir # Get location of padded elements in occupied and virtual space nonzero_opadding, nonzero_vpadding = padding_k_idx(cc, kind="split") ccsd_energy = cc.energy(t1, t2, eris) if t3p2_ip_out is None: t3p2_ip_out = np.zeros((nkpts,nkpts,nkpts,nocc,nvir,nocc,nocc),dtype=dtype) Wmcik = t3p2_ip_out if t3p2_ea_out is None: t3p2_ea_out = np.zeros((nkpts,nkpts,nkpts,nvir,nvir,nvir,nocc),dtype=dtype) Wacek = t3p2_ea_out from itertools import product tmp_t3 = np.empty((nkpts, nkpts, nkpts, nkpts, nkpts, nocc, nocc, nocc, nvir, nvir, nvir), dtype = t2.dtype) def get_w(ki, kj, kk, ka, kb, kc): kf = kconserv[ka,ki,kb] ret = lib.einsum('fiba,kjcf->ijkabc', eris.vovv[kf, ki, kb].conj(), t2[kk, kj, kc]) km = kconserv[kc,kk,kb] ret -= lib.einsum('jima,mkbc->ijkabc', eris.ooov[kj, ki, km].conj(), t2[km, kk, kb]) return ret for ki, kj, kk, ka, kb in product(range(nkpts), repeat=5): kc = kpts_helper.get_kconserv3(cc._scf.cell, cc.kpts, [ki, kj, kk, ka, kb]) tmp_t3[ki, kj, kk, ka, kb] = get_w(ki, kj, kk, ka, kb, kc) tmp_t3[ki, kj, kk, ka, kb] += get_w(ki, kk, kj, ka, kc, kb).transpose(0, 2, 1, 3, 5, 4) tmp_t3[ki, kj, kk, ka, kb] += get_w(kj, ki, kk, kb, ka, kc).transpose(1, 0, 2, 4, 3, 5) tmp_t3[ki, kj, kk, ka, kb] += get_w(kj, kk, ki, kb, kc, ka).transpose(2, 0, 1, 5, 3, 4) tmp_t3[ki, kj, kk, ka, kb] += get_w(kk, ki, kj, kc, ka, kb).transpose(1, 2, 0, 4, 5, 3) tmp_t3[ki, kj, kk, ka, kb] += get_w(kk, kj, ki, kc, kb, ka).transpose(2, 1, 0, 5, 4, 3) eijk = _get_epqr([0,nocc,ki,mo_e_o,nonzero_opadding], [0,nocc,kj,mo_e_o,nonzero_opadding], [0,nocc,kk,mo_e_o,nonzero_opadding]) eabc = _get_epqr([0,nvir,ka,mo_e_v,nonzero_vpadding], [0,nvir,kb,mo_e_v,nonzero_vpadding], [0,nvir,kc,mo_e_v,nonzero_vpadding], fac=[-1.,-1.,-1.]) eijkabc = eijk[:, :, :, None, None, None] + eabc[None, None, None, :, :, :] tmp_t3[ki, kj, kk, ka, kb] /= eijkabc pt1 = np.zeros((nkpts, nocc, nvir), dtype=t2.dtype) for ki in range(nkpts): for km, kn, ke in product(range(nkpts), repeat=3): kf = kconserv[km, ke, kn] Soovv = 2. * eris.oovv[km, kn, ke] - eris.oovv[km, kn, kf].transpose(0, 1, 3, 2) St3 = (tmp_t3[ki, km, kn, ki, ke] - tmp_t3[ki, km, kn, ke, ki].transpose(0, 1, 2, 4, 3, 5)) pt1[ki] += lib.einsum('mnef,imnaef->ia', Soovv, St3) pt2 = np.zeros((nkpts, nkpts, nkpts, nocc, nocc, nvir, nvir), dtype=t2.dtype) for ki, kj, ka in product(range(nkpts), repeat=3): kb = kconserv[ki, ka, kj] for km in range(nkpts): for kn in range(nkpts): # (ia,jb) -> (ia,jb) ke = kconserv[km, kj, kn] pt2[ki, kj, ka] += - 2. * lib.einsum('imnabe,mnje->ijab', tmp_t3[ki, km, kn, ka, kb], eris.ooov[km, kn, kj]) pt2[ki, kj, ka] += lib.einsum('imnabe,nmje->ijab', tmp_t3[ki, km, kn, ka, kb], eris.ooov[kn, km, kj]) pt2[ki, kj, ka] += lib.einsum('inmeab,mnje->ijab', tmp_t3[ki, kn, km, ke, ka], eris.ooov[km, kn, kj]) # (ia,jb) -> (jb,ia) ke = kconserv[km, ki, kn] pt2[ki, kj, ka] += - 2. * lib.einsum('jmnbae,mnie->ijab', tmp_t3[kj, km, kn, kb, ka], eris.ooov[km, kn, ki]) pt2[ki, kj, ka] += lib.einsum('jmnbae,nmie->ijab', tmp_t3[kj, km, kn, kb, ka], eris.ooov[kn, km, ki]) pt2[ki, kj, ka] += lib.einsum('jnmeba,mnie->ijab', tmp_t3[kj, kn, km, ke, kb], eris.ooov[km, kn, ki]) # (ia,jb) -> (ia,jb) pt2[ki, kj, ka] += lib.einsum('ijmabe,me->ijab', tmp_t3[ki, kj, km, ka, kb], fov[km]) pt2[ki, kj, ka] -= lib.einsum('ijmaeb,me->ijab', tmp_t3[ki, kj, km, ka, km], fov[km]) # (ia,jb) -> (jb,ia) pt2[ki, kj, ka] += lib.einsum('jimbae,me->ijab', tmp_t3[kj, ki, km, kb, ka], fov[km]) pt2[ki, kj, ka] -= lib.einsum('jimbea,me->ijab', tmp_t3[kj, ki, km, kb, km], fov[km]) for ke in range(nkpts): # (ia,jb) -> (ia,jb) kf = kconserv[km, ke, kb] pt2[ki, kj, ka] += 2. * lib.einsum('ijmaef,bmef->ijab', tmp_t3[ki, kj, km, ka, ke], eris.vovv[kb, km, ke]) pt2[ki, kj, ka] -= lib.einsum('ijmaef,bmfe->ijab', tmp_t3[ki, kj, km, ka, ke], eris.vovv[kb, km, kf]) pt2[ki, kj, ka] -= lib.einsum('imjfae,bmef->ijab', tmp_t3[ki, km, kj, kf, ka], eris.vovv[kb, km, ke]) # (ia,jb) -> (jb,ia) kf = kconserv[km, ke, ka] pt2[ki, kj, ka] += 2. * lib.einsum('jimbef,amef->ijab', tmp_t3[kj, ki, km, kb, ke], eris.vovv[ka, km, ke]) pt2[ki, kj, ka] -= lib.einsum('jimbef,amfe->ijab', tmp_t3[kj, ki, km, kb, ke], eris.vovv[ka, km, kf]) pt2[ki, kj, ka] -= lib.einsum('jmifbe,amef->ijab', tmp_t3[kj, km, ki, kf, kb], eris.vovv[ka, km, ke]) for ki in range(nkpts): ka = ki eia = LARGE_DENOM * np.ones((nocc, nvir), dtype=eris.mo_energy[0].dtype) n0_ovp_ia = np.ix_(nonzero_opadding[ki], nonzero_vpadding[ka]) eia[n0_ovp_ia] = (mo_e_o[ki][:,None] - mo_e_v[ka])[n0_ovp_ia] pt1[ki] /= eia for ki, ka in product(range(nkpts), repeat=2): eia = LARGE_DENOM * np.ones((nocc, nvir), dtype=eris.mo_energy[0].dtype) n0_ovp_ia = np.ix_(nonzero_opadding[ki], nonzero_vpadding[ka]) eia[n0_ovp_ia] = (mo_e_o[ki][:,None] - mo_e_v[ka])[n0_ovp_ia] for kj in range(nkpts): kb = kconserv[ki, ka, kj] ejb = LARGE_DENOM * np.ones((nocc, nvir), dtype=eris.mo_energy[0].dtype) n0_ovp_jb = np.ix_(nonzero_opadding[kj], nonzero_vpadding[kb]) ejb[n0_ovp_jb] = (mo_e_o[kj][:,None] - mo_e_v[kb])[n0_ovp_jb] eijab = eia[:, None, :, None] + ejb[:, None, :] pt2[ki, kj, ka] /= eijab pt1 += t1 pt2 += t2 for ki, kj, kk, ka, kb in product(range(nkpts), repeat=5): kc = kpts_helper.get_kconserv3(cc._scf.cell, cc.kpts, [ki, kj, kk, ka, kb]) km = kconserv[kc, ki, ka] _oovv = eris.oovv[km, ki, kc] Wmcik[km, kb, kk] += 2. * lib.einsum('ijkabc,mica->mbkj', tmp_t3[ki, kj, kk, ka, kb], _oovv) Wmcik[km, kb, kk] -= lib.einsum('jikabc,mica->mbkj', tmp_t3[kj, ki, kk, ka, kb], _oovv) Wmcik[km, kb, kk] -= lib.einsum('kjiabc,mica->mbkj', tmp_t3[kk, kj, ki, ka, kb], _oovv) for ki, kj, kk, ka, kb in product(range(nkpts), repeat=5): kc = kpts_helper.get_kconserv3(cc._scf.cell, cc.kpts, [ki, kj, kk, ka, kb]) ke = kconserv[ki, ka, kk] _oovv = eris.oovv[ki, kk, ka] Wacek[kc, kb, ke] -= 2. * lib.einsum('ijkabc,ikae->cbej', tmp_t3[ki, kj, kk, ka, kb], _oovv) Wacek[kc, kb, ke] += lib.einsum('jikabc,ikae->cbej', tmp_t3[kj, ki, kk, ka, kb], _oovv) Wacek[kc, kb, ke] += lib.einsum('kjiabc,ikae->cbej', tmp_t3[kk, kj, ki, ka, kb], _oovv) delta_ccsd_energy = cc.energy(pt1, pt2, eris) - ccsd_energy lib.logger.info(cc, 'CCSD energy T3[2] correction : %16.12e', delta_ccsd_energy) return delta_ccsd_energy, pt1, pt2, Wmcik, Wacek
def _add_pt2(pt2, nkpts, kconserv, kpt_indices, orb_indices, val): '''Adds term P(ia|jb)[tmp] to pt2. P(ia|jb)(tmp[i,j,a,b]) = tmp[i,j,a,b] + tmp[j,i,b,a] or equivalently for each i,j,a,b, pt2 is defined as pt2[i,j,a,b] += tmp[i,j,a,b] pt2[j,i,b,a] += tmp[i,j,a,b].transpose(1,0,3,2) If pt2 is lower-triangular, only adds the RHS term that contributes to the lower-triangular pt2. Args: pt2 (ndarray or HDF5 dataset): Full or lower triangular T2 array to which one is adding to. kpt_indices (array-like): K-point indices ki, kj, ka. orb_indices (array-like): Array-like of four tuples describing the range for i,j,a,b. An element of None will convert to slice(None,None). val (ndarray): Values to be added to pt2. ''' assert (len(orb_indices) == 4) ki, kj, ka = kpt_indices kb = kconserv[ki,ka,kj] idxi, idxj, idxa, idxb = [slice(None, None) if x is None else slice(x[0],x[1]) for x in orb_indices] if len(pt2.shape) == 7 and pt2.shape[:2] == (nkpts, nkpts): pt2[ki,kj,ka,idxi,idxj,idxa,idxb] += val pt2[kj,ki,kb,idxj,idxi,idxb,idxa] += val.transpose(1,0,3,2) elif len(pt2.shape) == 6 and pt2.shape[:2] == (nkpts*(nkpts+1)//2, nkpts): if ki <= kj: # Add tmp[i,j,a,b] to pt2[i,j,a,b] idx = (kj*(kj+1))//2 + ki pt2[idx,ka,idxi,idxj,idxa,idxb] += val if ki == kj: pt2[idx,kb,idxj,idxi,idxb,idxa] += val.transpose(1,0,3,2) else: # pt2[i,a,j,b] += tmp[j,i,a,b].transpose(1,0,3,2) idx = (ki*(ki+1))//2 + kj pt2[idx,kb,idxj,idxi,idxb,idxa] += val.transpose(1,0,3,2) else: raise ValueError('No known conversion for t2 shape %s' % pt2.shape)
[docs] def get_t3p2_imds(mycc, t1, t2, eris=None, t3p2_ip_out=None, t3p2_ea_out=None): """For a description of arguments, see `get_t3p2_imds_slow` in the corresponding `kintermediates.py`. """ from pyscf.pbc.cc.kccsd_t_rhf import _get_epqr cpu1 = cpu0 = (logger.process_clock(), logger.perf_counter()) if eris is None: eris = mycc.ao2mo() fock = eris.fock nkpts, nocc, nvir = t1.shape cell = mycc._scf.cell kpts = mycc.kpts kconserv = mycc.khelper.kconserv dtype = np.result_type(t1, t2) fov = fock[:, :nocc, nocc:] #foo = np.asarray([fock[ikpt, :nocc, :nocc].diagonal() for ikpt in range(nkpts)]) #fvv = np.asarray([fock[ikpt, nocc:, nocc:].diagonal() for ikpt in range(nkpts)]) mo_energy_occ = np.array([eris.mo_energy[ki][:nocc] for ki in range(nkpts)]) mo_energy_vir = np.array([eris.mo_energy[ki][nocc:] for ki in range(nkpts)]) mo_e_o = mo_energy_occ mo_e_v = mo_energy_vir ccsd_energy = mycc.energy(t1, t2, eris) if t3p2_ip_out is None: t3p2_ip_out = np.zeros((nkpts,nkpts,nkpts,nocc,nvir,nocc,nocc),dtype=dtype) Wmcik = t3p2_ip_out if t3p2_ea_out is None: t3p2_ea_out = np.zeros((nkpts,nkpts,nkpts,nvir,nvir,nvir,nocc),dtype=dtype) Wacek = t3p2_ea_out # Create necessary temporary eris for fast read from pyscf.pbc.cc.kccsd_t_rhf import create_t3_eris, get_data_slices feri_tmp, t2T, eris_vvop, eris_vooo_C = create_t3_eris(mycc, kconserv, [eris.vovv, eris.oovv, eris.ooov, t2]) #t1T = np.array([x.T for x in t1], dtype=np.complex128, order='C') #fvo = np.array([x.T for x in fov], dtype=np.complex128, order='C') cpu1 = logger.timer_debug1(mycc, 'CCSD(T) tmp eri creation', *cpu1) def get_w(ki, kj, kk, ka, kb, kc, a0, a1, b0, b1, c0, c1): '''Wijkabc intermediate as described in Scuseria paper before Pijkabc acts Function copied for `kccsd_t_rhf.py`''' km = kconserv[kc, kk, kb] kf = kconserv[kk, kc, kj] out = einsum('cfjk,abif->abcijk', t2T[kc,kf,kj,c0:c1,:,:,:], eris_vvop[ka,kb,ki,a0:a1,b0:b1,:,nocc:]) out = out - einsum('cbmk,aijm->abcijk', t2T[kc,kb,km,c0:c1,b0:b1,:,:], eris_vooo_C[ka,ki,kj,a0:a1,:,:,:]) return out def get_permuted_w(ki, kj, kk, ka, kb, kc, orb_indices): '''Pijkabc operating on Wijkabc intermediate as described in Scuseria paper Function copied for `kccsd_t_rhf.py`''' a0, a1, b0, b1, c0, c1 = orb_indices out = get_w(ki, kj, kk, ka, kb, kc, a0, a1, b0, b1, c0, c1) out = out + get_w(kj, kk, ki, kb, kc, ka, b0, b1, c0, c1, a0, a1).transpose(2,0,1,5,3,4) out = out + get_w(kk, ki, kj, kc, ka, kb, c0, c1, a0, a1, b0, b1).transpose(1,2,0,4,5,3) out = out + get_w(ki, kk, kj, ka, kc, kb, a0, a1, c0, c1, b0, b1).transpose(0,2,1,3,5,4) out = out + get_w(kk, kj, ki, kc, kb, ka, c0, c1, b0, b1, a0, a1).transpose(2,1,0,5,4,3) out = out + get_w(kj, ki, kk, kb, ka, kc, b0, b1, a0, a1, c0, c1).transpose(1,0,2,4,3,5) return out def get_data(kpt_indices): idx_args = get_data_slices(kpt_indices, task, kconserv) vvop_indices, vooo_indices, t2T_vvop_indices, t2T_vooo_indices = idx_args vvop_data = [eris_vvop[tuple(x)] for x in vvop_indices] vooo_data = [eris_vooo_C[tuple(x)] for x in vooo_indices] t2T_vvop_data = [t2T[tuple(x)] for x in t2T_vvop_indices] t2T_vooo_data = [t2T[tuple(x)] for x in t2T_vooo_indices] data = [vvop_data, vooo_data, t2T_vvop_data, t2T_vooo_data] return data def add_and_permute(kpt_indices, orb_indices, data): '''Performs permutation and addition of t3 temporary arrays.''' ki, kj, kk, ka, kb, kc = kpt_indices a0, a1, b0, b1, c0, c1 = orb_indices tmp_t3Tv_ijk = np.asarray(data[0], dtype=dtype, order='C') tmp_t3Tv_jik = np.asarray(data[1], dtype=dtype, order='C') tmp_t3Tv_kji = np.asarray(data[2], dtype=dtype, order='C') #out_ijk = np.empty(data[0].shape, dtype=dtype, order='C') #drv = _ccsd.libcc.MPICCadd_and_permute_t3T #drv(ctypes.c_int(nocc), ctypes.c_int(nvir), # ctypes.c_int(0), # out_ijk.ctypes.data_as(ctypes.c_void_p), # tmp_t3Tv_ijk.ctypes.data_as(ctypes.c_void_p), # tmp_t3Tv_jik.ctypes.data_as(ctypes.c_void_p), # tmp_t3Tv_kji.ctypes.data_as(ctypes.c_void_p), # mo_offset.ctypes.data_as(ctypes.c_void_p), # slices.ctypes.data_as(ctypes.c_void_p)) return (2.*tmp_t3Tv_ijk - tmp_t3Tv_jik.transpose(0,1,2,4,3,5) - tmp_t3Tv_kji.transpose(0,1,2,5,4,3)) #return out_ijk # Get location of padded elements in occupied and virtual space nonzero_opadding, nonzero_vpadding = padding_k_idx(mycc, kind="split") mem_now = lib.current_memory()[0] max_memory = max(0, mycc.max_memory - mem_now) blkmin = 4 # temporary t3 array is size: nkpts**3 * blksize**3 * nocc**3 * 16 vir_blksize = min(nvir, max(blkmin, int((max_memory*.9e6/16/nocc**3/nkpts**3)**(1./3)))) tasks = [] logger.debug(mycc, 'max_memory %d MB (%d MB in use)', max_memory, mem_now) logger.debug(mycc, 'virtual blksize = %d (nvir = %d)', vir_blksize, nvir) for a0, a1 in lib.prange(0, nvir, vir_blksize): for b0, b1 in lib.prange(0, nvir, vir_blksize): for c0, c1 in lib.prange(0, nvir, vir_blksize): tasks.append((a0,a1,b0,b1,c0,c1)) eaa = [] for ka in range(nkpts): eaa.append(mo_e_o[ka][:, None] - mo_e_v[ka][None, :]) pt1 = np.zeros((nkpts,nocc,nvir), dtype=dtype) pt2 = np.zeros((nkpts,nkpts,nkpts,nocc,nocc,nvir,nvir), dtype=dtype) for ka, kb in product(range(nkpts), repeat=2): for task_id, task in enumerate(tasks): cput2 = (logger.process_clock(), logger.perf_counter()) a0,a1,b0,b1,c0,c1 = task my_permuted_w = np.zeros((nkpts,)*3 + (a1-a0,b1-b0,c1-c0) + (nocc,)*3, dtype=dtype) for ki, kj, kk in product(range(nkpts), repeat=3): # Find momentum conservation condition for triples # amplitude t3ijkabc kc = kpts_helper.get_kconserv3(cell, kpts, [ki, kj, kk, ka, kb]) kpt_indices = [ki,kj,kk,ka,kb,kc] #data = get_data(kpt_indices) my_permuted_w[ki,kj,kk] = get_permuted_w(ki,kj,kk,ka,kb,kc,task) for ki, kj, kk in product(range(nkpts), repeat=3): # eigenvalue denominator: e(i) + e(j) + e(k) eijk = _get_epqr([0,nocc,ki,mo_e_o,nonzero_opadding], [0,nocc,kj,mo_e_o,nonzero_opadding], [0,nocc,kk,mo_e_o,nonzero_opadding]) # Find momentum conservation condition for triples # amplitude t3ijkabc kc = kpts_helper.get_kconserv3(cell, kpts, [ki, kj, kk, ka, kb]) eabc = _get_epqr([a0,a1,ka,mo_e_v,nonzero_vpadding], [b0,b1,kb,mo_e_v,nonzero_vpadding], [c0,c1,kc,mo_e_v,nonzero_vpadding], fac=[-1.,-1.,-1.]) kpt_indices = [ki,kj,kk,ka,kb,kc] eabcijk = (eijk[None,None,None,:,:,:] + eabc[:,:,:,None,None,None]) tmp_t3Tv_ijk = my_permuted_w[ki,kj,kk] tmp_t3Tv_jik = my_permuted_w[kj,ki,kk] tmp_t3Tv_kji = my_permuted_w[kk,kj,ki] Ptmp_t3Tv = add_and_permute(kpt_indices, task, (tmp_t3Tv_ijk,tmp_t3Tv_jik,tmp_t3Tv_kji)) Ptmp_t3Tv /= eabcijk # Contribution to T1 amplitudes if ki == ka and kc == kconserv[kj, kb, kk]: eris_Soovv = (2.*eris.oovv[kj,kk,kb,:,:,b0:b1,c0:c1] - eris.oovv[kj,kk,kc,:,:,c0:c1,b0:b1].transpose(0,1,3,2)) pt1[ka,:,a0:a1] += 0.5*einsum('abcijk,jkbc->ia', Ptmp_t3Tv, eris_Soovv) # Contribution to T2 amplitudes if ki == ka and kc == kconserv[kj, kb, kk]: tmp = einsum('abcijk,ia->jkbc', Ptmp_t3Tv, 0.5*fov[ki,:,a0:a1]) _add_pt2(pt2, nkpts, kconserv, [kj,kk,kb], [None,None,(b0,b1),(c0,c1)], tmp) kd = kconserv[ka,ki,kb] eris_vovv = eris.vovv[kd,ki,kb,:,:,b0:b1,a0:a1] tmp = einsum('abcijk,diba->jkdc', Ptmp_t3Tv, eris_vovv) _add_pt2(pt2, nkpts, kconserv, [kj,kk,kd], [None,None,None,(c0,c1)], tmp) km = kconserv[kc, kk, kb] eris_ooov = eris.ooov[kj,ki,km,:,:,:,a0:a1] tmp = einsum('abcijk,jima->mkbc', Ptmp_t3Tv, eris_ooov) _add_pt2(pt2, nkpts, kconserv, [km,kk,kb], [None,None,(b0,b1),(c0,c1)], -1.*tmp) # Contribution to Wovoo array km = kconserv[ka,ki,kc] eris_oovv = eris.oovv[km,ki,kc,:,:,c0:c1,a0:a1] tmp = einsum('abcijk,mica->mbkj', Ptmp_t3Tv, eris_oovv) Wmcik[km,kb,kk,:,b0:b1,:,:] += tmp # Contribution to Wvvoo array ke = kconserv[ki,ka,kk] eris_oovv = eris.oovv[ki,kk,ka,:,:,a0:a1,:] tmp = einsum('abcijk,ikae->cbej', Ptmp_t3Tv, eris_oovv) Wacek[kc,kb,ke,c0:c1,b0:b1,:,:] -= tmp logger.timer_debug1(mycc, 'EOM-CCSD T3[2] ka,kb,vir=(%d,%d,%d/%d) [total=%d]'% (ka,kb,task_id,len(tasks),nkpts**5), *cput2) for ki in range(nkpts): ka = ki eia = LARGE_DENOM * np.ones((nocc, nvir), dtype=eris.mo_energy[0].dtype) n0_ovp_ia = np.ix_(nonzero_opadding[ki], nonzero_vpadding[ka]) eia[n0_ovp_ia] = (mo_e_o[ki][:,None] - mo_e_v[ka])[n0_ovp_ia] pt1[ki] /= eia for ki, ka in product(range(nkpts), repeat=2): eia = LARGE_DENOM * np.ones((nocc, nvir), dtype=eris.mo_energy[0].dtype) n0_ovp_ia = np.ix_(nonzero_opadding[ki], nonzero_vpadding[ka]) eia[n0_ovp_ia] = (mo_e_o[ki][:,None] - mo_e_v[ka])[n0_ovp_ia] for kj in range(nkpts): kb = kconserv[ki, ka, kj] ejb = LARGE_DENOM * np.ones((nocc, nvir), dtype=eris.mo_energy[0].dtype) n0_ovp_jb = np.ix_(nonzero_opadding[kj], nonzero_vpadding[kb]) ejb[n0_ovp_jb] = (mo_e_o[kj][:,None] - mo_e_v[kb])[n0_ovp_jb] eijab = eia[:, None, :, None] + ejb[:, None, :] pt2[ki, kj, ka] /= eijab pt1 += t1 pt2 += t2 logger.timer(mycc, 'EOM-CCSD(T) imds', *cpu0) delta_ccsd_energy = mycc.energy(pt1, pt2, eris) - ccsd_energy logger.info(mycc, 'CCSD energy T3[2] correction : %16.12e', delta_ccsd_energy) return delta_ccsd_energy, pt1, pt2, Wmcik, Wacek