Source code for pyscf.cc.uccsd_lambda

#!/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.
#
# Author: Qiming Sun <osirpt.sun@gmail.com>
#


import numpy
from pyscf import lib
from pyscf.lib import logger
from pyscf.cc import uccsd
from pyscf.cc import ccsd_lambda

einsum = lib.einsum

[docs] def kernel(mycc, eris=None, t1=None, t2=None, l1=None, l2=None, max_cycle=50, tol=1e-8, verbose=logger.INFO): if eris is None: eris = mycc.ao2mo() return ccsd_lambda.kernel(mycc, eris, t1, t2, l1, l2, max_cycle, tol, verbose, make_intermediates, update_lambda)
# l2, t2 as ijab
[docs] def make_intermediates(mycc, t1, t2, eris): t1a, t1b = t1 t2aa, t2ab, t2bb = t2 nocca, nvira = t1a.shape noccb, nvirb = t1b.shape fooa = eris.focka[:nocca,:nocca] fova = eris.focka[:nocca,nocca:] fvoa = eris.focka[nocca:,:nocca] fvva = eris.focka[nocca:,nocca:] foob = eris.fockb[:noccb,:noccb] fovb = eris.fockb[:noccb,noccb:] fvob = eris.fockb[noccb:,:noccb] fvvb = eris.fockb[noccb:,noccb:] tauaa, tauab, taubb = uccsd.make_tau(t2, t1, t1) ovov = numpy.asarray(eris.ovov) ovov = ovov - ovov.transpose(0,3,2,1) OVOV = numpy.asarray(eris.OVOV) OVOV = OVOV - OVOV.transpose(0,3,2,1) ovOV = numpy.asarray(eris.ovOV) v1a = fvva - einsum('ja,jb->ba', fova, t1a) v1b = fvvb - einsum('ja,jb->ba', fovb, t1b) v1a += einsum('jcka,jkbc->ba', ovov, tauaa) * .5 v1a -= einsum('jaKC,jKbC->ba', ovOV, tauab) * .5 v1a -= einsum('kaJC,kJbC->ba', ovOV, tauab) * .5 v1b += einsum('jcka,jkbc->ba', OVOV, taubb) * .5 v1b -= einsum('kcJA,kJcB->BA', ovOV, tauab) * .5 v1b -= einsum('jcKA,jKcB->BA', ovOV, tauab) * .5 v2a = fooa + einsum('ib,jb->ij', fova, t1a) v2b = foob + einsum('ib,jb->ij', fovb, t1b) v2a += einsum('ibkc,jkbc->ij', ovov, tauaa) * .5 v2a += einsum('ibKC,jKbC->ij', ovOV, tauab) v2b += einsum('ibkc,jkbc->ij', OVOV, taubb) * .5 v2b += einsum('kcIB,kJcB->IJ', ovOV, tauab) ovoo = numpy.asarray(eris.ovoo) ovoo = ovoo - ovoo.transpose(2,1,0,3) OVOO = numpy.asarray(eris.OVOO) OVOO = OVOO - OVOO.transpose(2,1,0,3) OVoo = numpy.asarray(eris.OVoo) ovOO = numpy.asarray(eris.ovOO) v2a -= numpy.einsum('ibkj,kb->ij', ovoo, t1a) v2a += numpy.einsum('KBij,KB->ij', OVoo, t1b) v2b -= numpy.einsum('ibkj,kb->ij', OVOO, t1b) v2b += numpy.einsum('kbIJ,kb->IJ', ovOO, t1a) v5a = fvoa + numpy.einsum('kc,jkbc->bj', fova, t2aa) v5a += numpy.einsum('KC,jKbC->bj', fovb, t2ab) v5b = fvob + numpy.einsum('kc,jkbc->bj', fovb, t2bb) v5b += numpy.einsum('kc,kJcB->BJ', fova, t2ab) tmp = fova - numpy.einsum('kdlc,ld->kc', ovov, t1a) tmp += numpy.einsum('kcLD,LD->kc', ovOV, t1b) v5a += einsum('kc,kb,jc->bj', tmp, t1a, t1a) tmp = fovb - numpy.einsum('kdlc,ld->kc', OVOV, t1b) tmp += numpy.einsum('ldKC,ld->KC', ovOV, t1a) v5b += einsum('kc,kb,jc->bj', tmp, t1b, t1b) v5a -= einsum('lckj,klbc->bj', ovoo, t2aa) * .5 v5a -= einsum('LCkj,kLbC->bj', OVoo, t2ab) v5b -= einsum('LCKJ,KLBC->BJ', OVOO, t2bb) * .5 v5b -= einsum('lcKJ,lKcB->BJ', ovOO, t2ab) oooo = numpy.asarray(eris.oooo) OOOO = numpy.asarray(eris.OOOO) ooOO = numpy.asarray(eris.ooOO) woooo = einsum('icjl,kc->ikjl', ovoo, t1a) wOOOO = einsum('icjl,kc->ikjl', OVOO, t1b) wooOO = einsum('icJL,kc->ikJL', ovOO, t1a) wooOO += einsum('JCil,KC->ilJK', OVoo, t1b) woooo += (oooo - oooo.transpose(0,3,2,1)) * .5 wOOOO += (OOOO - OOOO.transpose(0,3,2,1)) * .5 wooOO += ooOO.copy() woooo += einsum('icjd,klcd->ikjl', ovov, tauaa) * .25 wOOOO += einsum('icjd,klcd->ikjl', OVOV, taubb) * .25 wooOO += einsum('icJD,kLcD->ikJL', ovOV, tauab) v4ovvo = einsum('jbld,klcd->jbck', ovov, t2aa) v4ovvo += einsum('jbLD,kLcD->jbck', ovOV, t2ab) v4ovvo += numpy.asarray(eris.ovvo) v4ovvo -= numpy.asarray(eris.oovv).transpose(0,3,2,1) v4OVVO = einsum('jbld,klcd->jbck', OVOV, t2bb) v4OVVO += einsum('ldJB,lKdC->JBCK', ovOV, t2ab) v4OVVO += numpy.asarray(eris.OVVO) v4OVVO -= numpy.asarray(eris.OOVV).transpose(0,3,2,1) v4OVvo = einsum('ldJB,klcd->JBck', ovOV, t2aa) v4OVvo += einsum('JBLD,kLcD->JBck', OVOV, t2ab) v4OVvo += numpy.asarray(eris.OVvo) v4ovVO = einsum('jbLD,KLCD->jbCK', ovOV, t2bb) v4ovVO += einsum('jbld,lKdC->jbCK', ovov, t2ab) v4ovVO += numpy.asarray(eris.ovVO) v4oVVo = einsum('jdLB,kLdC->jBCk', ovOV, t2ab) v4oVVo -= numpy.asarray(eris.ooVV).transpose(0,3,2,1) v4OvvO = einsum('lbJD,lKcD->JbcK', ovOV, t2ab) v4OvvO -= numpy.asarray(eris.OOvv).transpose(0,3,2,1) woovo = einsum('ibck,jb->ijck', v4ovvo, t1a) wOOVO = einsum('ibck,jb->ijck', v4OVVO, t1b) wOOvo = einsum('IBck,JB->IJck', v4OVvo, t1b) wOOvo -= einsum('IbcK,jb->IKcj', v4OvvO, t1a) wooVO = einsum('ibCK,jb->ijCK', v4ovVO, t1a) wooVO -= einsum('iBCk,JB->ikCJ', v4oVVo, t1b) woovo += ovoo.conj().transpose(3,2,1,0) * .5 wOOVO += OVOO.conj().transpose(3,2,1,0) * .5 wooVO += OVoo.conj().transpose(3,2,1,0) wOOvo += ovOO.conj().transpose(3,2,1,0) woovo -= einsum('iclk,jlbc->ikbj', ovoo, t2aa) woovo += einsum('LCik,jLbC->ikbj', OVoo, t2ab) wOOVO -= einsum('iclk,jlbc->ikbj', OVOO, t2bb) wOOVO += einsum('lcIK,lJcB->IKBJ', ovOO, t2ab) wooVO -= einsum('iclk,lJcB->ikBJ', ovoo, t2ab) wooVO += einsum('LCik,JLBC->ikBJ', OVoo, t2bb) wooVO -= einsum('icLK,jLcB->ijBK', ovOO, t2ab) wOOvo -= einsum('ICLK,jLbC->IKbj', OVOO, t2ab) wOOvo += einsum('lcIK,jlbc->IKbj', ovOO, t2aa) wOOvo -= einsum('IClk,lJbC->IJbk', OVoo, t2ab) wvvvo = einsum('jack,jb->back', v4ovvo, t1a) wVVVO = einsum('jack,jb->back', v4OVVO, t1b) wVVvo = einsum('JAck,JB->BAck', v4OVvo, t1b) wVVvo -= einsum('jACk,jb->CAbk', v4oVVo, t1a) wvvVO = einsum('jaCK,jb->baCK', v4ovVO, t1a) wvvVO -= einsum('JacK,JB->caBK', v4OvvO, t1b) wvvvo += einsum('lajk,jlbc->back', .25*ovoo, tauaa) wVVVO += einsum('lajk,jlbc->back', .25*OVOO, taubb) wVVvo -= einsum('LAjk,jLcB->BAck', OVoo, tauab) wvvVO -= einsum('laJK,lJbC->baCK', ovOO, tauab) w3a = numpy.einsum('jbck,jb->ck', v4ovvo, t1a) w3a += numpy.einsum('JBck,JB->ck', v4OVvo, t1b) w3b = numpy.einsum('jbck,jb->ck', v4OVVO, t1b) w3b += numpy.einsum('jbCK,jb->CK', v4ovVO, t1a) wovvo = v4ovvo wOVVO = v4OVVO wovVO = v4ovVO wOVvo = v4OVvo woVVo = v4oVVo wOvvO = v4OvvO wovvo += lib.einsum('jbld,kd,lc->jbck', ovov, t1a, -t1a) wOVVO += lib.einsum('jbld,kd,lc->jbck', OVOV, t1b, -t1b) wovVO += lib.einsum('jbLD,KD,LC->jbCK', ovOV, t1b, -t1b) wOVvo += lib.einsum('ldJB,kd,lc->JBck', ovOV, t1a, -t1a) woVVo += lib.einsum('jdLB,kd,LC->jBCk', ovOV, t1a, t1b) wOvvO += lib.einsum('lbJD,KD,lc->JbcK', ovOV, t1b, t1a) wovvo -= einsum('jblk,lc->jbck', ovoo, t1a) wOVVO -= einsum('jblk,lc->jbck', OVOO, t1b) wovVO -= einsum('jbLK,LC->jbCK', ovOO, t1b) wOVvo -= einsum('JBlk,lc->JBck', OVoo, t1a) woVVo += einsum('LBjk,LC->jBCk', OVoo, t1b) wOvvO += einsum('lbJK,lc->JbcK', ovOO, t1a) if nvira > 0 and nocca > 0: ovvv = numpy.asarray(eris.get_ovvv()) ovvv = ovvv - ovvv.transpose(0,3,2,1) v1a -= numpy.einsum('jabc,jc->ba', ovvv, t1a) v5a += einsum('kdbc,jkcd->bj', ovvv, t2aa) * .5 woovo += einsum('idcb,kjbd->ijck', ovvv, tauaa) * .25 wovvo += einsum('jbcd,kd->jbck', ovvv, t1a) wvvvo -= ovvv.conj().transpose(3,2,1,0) * .5 wvvvo += einsum('jacd,kjbd->cabk', ovvv, t2aa) wvvVO += einsum('jacd,jKdB->caBK', ovvv, t2ab) ovvv = tmp = None if nvirb > 0 and noccb > 0: OVVV = numpy.asarray(eris.get_OVVV()) OVVV = OVVV - OVVV.transpose(0,3,2,1) v1b -= numpy.einsum('jabc,jc->ba', OVVV, t1b) v5b += einsum('KDBC,JKCD->BJ', OVVV, t2bb) * .5 wOOVO += einsum('idcb,kjbd->ijck', OVVV, taubb) * .25 wOVVO += einsum('jbcd,kd->jbck', OVVV, t1b) wVVVO -= OVVV.conj().transpose(3,2,1,0) * .5 wVVVO += einsum('jacd,kjbd->cabk', OVVV, t2bb) wVVvo += einsum('JACD,kJbD->CAbk', OVVV, t2ab) OVVV = tmp = None if nvirb > 0 and nocca > 0: OVvv = numpy.asarray(eris.get_OVvv()) v1a += numpy.einsum('JCba,JC->ba', OVvv, t1b) v5a += einsum('KDbc,jKcD->bj', OVvv, t2ab) wOOvo += einsum('IDcb,kJbD->IJck', OVvv, tauab) wOVvo += einsum('JBcd,kd->JBck', OVvv, t1a) wOvvO -= einsum('JDcb,KD->JbcK', OVvv, t1b) wvvVO -= OVvv.conj().transpose(3,2,1,0) wvvvo -= einsum('KDca,jKbD->cabj', OVvv, t2ab) wvvVO -= einsum('KDca,JKBD->caBJ', OVvv, t2bb) wVVvo += einsum('KAcd,jKdB->BAcj', OVvv, t2ab) OVvv = tmp = None if nvira > 0 and noccb > 0: ovVV = numpy.asarray(eris.get_ovVV()) v1b += numpy.einsum('jcBA,jc->BA', ovVV, t1a) v5b += einsum('kdBC,kJdC->BJ', ovVV, t2ab) wooVO += einsum('idCB,jKdB->ijCK', ovVV, tauab) wovVO += einsum('jbCD,KD->jbCK', ovVV, t1b) woVVo -= einsum('jdCB,kd->jBCk', ovVV, t1a) wVVvo -= ovVV.conj().transpose(3,2,1,0) wVVVO -= einsum('kdCA,kJdB->CABJ', ovVV, t2ab) wVVvo -= einsum('kdCA,jkbd->CAbj', ovVV, t2aa) wvvVO += einsum('kaCD,kJbD->baCJ', ovVV, t2ab) ovVV = tmp = None w3a += v5a w3b += v5b w3a += lib.einsum('cb,jb->cj', v1a, t1a) w3b += lib.einsum('cb,jb->cj', v1b, t1b) w3a -= lib.einsum('jk,jb->bk', v2a, t1a) w3b -= lib.einsum('jk,jb->bk', v2b, t1b) class _IMDS: pass imds = _IMDS() imds.ftmp = lib.H5TmpFile() dtype = numpy.result_type(t2ab, eris.vvvv).char imds.woooo = imds.ftmp.create_dataset('woooo', (nocca,nocca,nocca,nocca), dtype) imds.wooOO = imds.ftmp.create_dataset('wooOO', (nocca,nocca,noccb,noccb), dtype) imds.wOOOO = imds.ftmp.create_dataset('wOOOO', (noccb,noccb,noccb,noccb), dtype) imds.wovvo = imds.ftmp.create_dataset('wovvo', (nocca,nvira,nvira,nocca), dtype) imds.wOVVO = imds.ftmp.create_dataset('wOVVO', (noccb,nvirb,nvirb,noccb), dtype) imds.wovVO = imds.ftmp.create_dataset('wovVO', (nocca,nvira,nvirb,noccb), dtype) imds.wOVvo = imds.ftmp.create_dataset('wOVvo', (noccb,nvirb,nvira,nocca), dtype) imds.woVVo = imds.ftmp.create_dataset('woVVo', (nocca,nvirb,nvirb,nocca), dtype) imds.wOvvO = imds.ftmp.create_dataset('wOvvO', (noccb,nvira,nvira,noccb), dtype) imds.woovo = imds.ftmp.create_dataset('woovo', (nocca,nocca,nvira,nocca), dtype) imds.wOOVO = imds.ftmp.create_dataset('wOOVO', (noccb,noccb,nvirb,noccb), dtype) imds.wOOvo = imds.ftmp.create_dataset('wOOvo', (noccb,noccb,nvira,nocca), dtype) imds.wooVO = imds.ftmp.create_dataset('wooVO', (nocca,nocca,nvirb,noccb), dtype) imds.wvvvo = imds.ftmp.create_dataset('wvvvo', (nvira,nvira,nvira,nocca), dtype) imds.wVVVO = imds.ftmp.create_dataset('wVVVO', (nvirb,nvirb,nvirb,noccb), dtype) imds.wVVvo = imds.ftmp.create_dataset('wVVvo', (nvirb,nvirb,nvira,nocca), dtype) imds.wvvVO = imds.ftmp.create_dataset('wvvVO', (nvira,nvira,nvirb,noccb), dtype) imds.woooo[:] = woooo imds.wOOOO[:] = wOOOO imds.wooOO[:] = wooOO imds.wovvo[:] = wovvo imds.wOVVO[:] = wOVVO imds.wovVO[:] = wovVO imds.wOVvo[:] = wOVvo imds.woVVo[:] = woVVo imds.wOvvO[:] = wOvvO imds.woovo[:] = woovo imds.wOOVO[:] = wOOVO imds.wOOvo[:] = wOOvo imds.wooVO[:] = wooVO imds.wvvvo[:] = wvvvo imds.wVVVO[:] = wVVVO imds.wVVvo[:] = wVVvo imds.wvvVO[:] = wvvVO imds.v1a = v1a imds.v1b = v1b imds.v2a = v2a imds.v2b = v2b imds.w3a = w3a imds.w3b = w3b imds.ftmp.flush() return imds
# update L1, L2
[docs] def update_lambda(mycc, t1, t2, l1, l2, eris, imds): time0 = logger.process_clock(), logger.perf_counter() log = logger.Logger(mycc.stdout, mycc.verbose) t1a, t1b = t1 t2aa, t2ab, t2bb = t2 l1a, l1b = l1 l2aa, l2ab, l2bb = l2 nocca, nvira = t1a.shape noccb, nvirb = t1b.shape u1a = numpy.zeros_like(l1a) u1b = numpy.zeros_like(l1b) u2aa = numpy.zeros_like(l2aa) u2ab = numpy.zeros_like(l2ab) u2bb = numpy.zeros_like(l2bb) mo_ea_o = eris.mo_energy[0][:nocca] mo_ea_v = eris.mo_energy[0][nocca:] + mycc.level_shift mo_eb_o = eris.mo_energy[1][:noccb] mo_eb_v = eris.mo_energy[1][noccb:] + mycc.level_shift fova = eris.focka[:nocca,nocca:] fovb = eris.fockb[:noccb,noccb:] v1a = imds.v1a - numpy.diag(mo_ea_v) v1b = imds.v1b - numpy.diag(mo_eb_v) v2a = imds.v2a - numpy.diag(mo_ea_o) v2b = imds.v2b - numpy.diag(mo_eb_o) mvv = einsum('klca,klcb->ba', l2aa, t2aa) * .5 mvv+= einsum('lKaC,lKbC->ba', l2ab, t2ab) mVV = einsum('klca,klcb->ba', l2bb, t2bb) * .5 mVV+= einsum('kLcA,kLcB->BA', l2ab, t2ab) moo = einsum('kicd,kjcd->ij', l2aa, t2aa) * .5 moo+= einsum('iKdC,jKdC->ij', l2ab, t2ab) mOO = einsum('kicd,kjcd->ij', l2bb, t2bb) * .5 mOO+= einsum('kIcD,kJcD->IJ', l2ab, t2ab) #m3 = lib.einsum('ijcd,cdab->ijab', l2, eris.vvvv) * .5 m3aa, m3ab, m3bb = mycc._add_vvvv(None, (l2aa.conj(),l2ab.conj(),l2bb.conj()), eris) m3aa = m3aa.conj() m3ab = m3ab.conj() m3bb = m3bb.conj() m3aa += lib.einsum('klab,ikjl->ijab', l2aa, numpy.asarray(imds.woooo)) m3bb += lib.einsum('klab,ikjl->ijab', l2bb, numpy.asarray(imds.wOOOO)) m3ab += lib.einsum('kLaB,ikJL->iJaB', l2ab, numpy.asarray(imds.wooOO)) ovov = numpy.asarray(eris.ovov) ovov = ovov - ovov.transpose(0,3,2,1) OVOV = numpy.asarray(eris.OVOV) OVOV = OVOV - OVOV.transpose(0,3,2,1) ovOV = numpy.asarray(eris.ovOV) mvv1 = einsum('jc,jb->bc', l1a, t1a) + mvv mVV1 = einsum('jc,jb->bc', l1b, t1b) + mVV moo1 = einsum('ic,kc->ik', l1a, t1a) + moo mOO1 = einsum('ic,kc->ik', l1b, t1b) + mOO if nvira > 0 and nocca > 0: ovvv = numpy.asarray(eris.get_ovvv()) ovvv = ovvv - ovvv.transpose(0,3,2,1) tmp = lib.einsum('ijcd,kd->ijck', l2aa, t1a) m3aa -= lib.einsum('kbca,ijck->ijab', ovvv, tmp) tmp = einsum('ic,jbca->jiba', l1a, ovvv) tmp+= einsum('kiab,jk->ijab', l2aa, v2a) tmp-= einsum('ik,kajb->ijab', moo1, ovov) u2aa += tmp - tmp.transpose(1,0,2,3) u1a += numpy.einsum('iacb,bc->ia', ovvv, mvv1) ovvv = tmp = None if nvirb > 0 and noccb > 0: OVVV = numpy.asarray(eris.get_OVVV()) OVVV = OVVV - OVVV.transpose(0,3,2,1) tmp = lib.einsum('ijcd,kd->ijck', l2bb, t1b) m3bb -= lib.einsum('kbca,ijck->ijab', OVVV, tmp) tmp = einsum('ic,jbca->jiba', l1b, OVVV) tmp+= einsum('kiab,jk->ijab', l2bb, v2b) tmp-= einsum('ik,kajb->ijab', mOO1, OVOV) u2bb += tmp - tmp.transpose(1,0,2,3) u1b += numpy.einsum('iaCB,BC->ia', OVVV, mVV1) OVVV = tmp = None if nvirb > 0 and nocca > 0: OVvv = numpy.asarray(eris.get_OVvv()) tmp = lib.einsum('iJcD,KD->iJcK', l2ab, t1b) m3ab -= lib.einsum('KBca,iJcK->iJaB', OVvv, tmp) tmp = einsum('ic,JAcb->JibA', l1a, OVvv) tmp-= einsum('kIaB,jk->IjaB', l2ab, v2a) tmp-= einsum('IK,jaKB->IjaB', mOO1, ovOV) u2ab += tmp.transpose(1,0,2,3) u1b += numpy.einsum('iacb,bc->ia', OVvv, mvv1) OVvv = tmp = None if nvira > 0 and noccb > 0: ovVV = numpy.asarray(eris.get_ovVV()) tmp = lib.einsum('iJdC,kd->iJCk', l2ab, t1a) m3ab -= lib.einsum('kaCB,iJCk->iJaB', ovVV, tmp) tmp = einsum('IC,jbCA->jIbA', l1b, ovVV) tmp-= einsum('iKaB,JK->iJaB', l2ab, v2b) tmp-= einsum('ik,kaJB->iJaB', moo1, ovOV) u2ab += tmp u1a += numpy.einsum('iaCB,BC->ia', ovVV, mVV1) ovVV = tmp = None tauaa, tauab, taubb = uccsd.make_tau(t2, t1, t1) tmp = lib.einsum('ijcd,klcd->ijkl', l2aa, tauaa) ovov = numpy.asarray(eris.ovov) ovov = ovov - ovov.transpose(0,3,2,1) m3aa += lib.einsum('kalb,ijkl->ijab', ovov, tmp) * .25 tmp = lib.einsum('ijcd,klcd->ijkl', l2bb, taubb) OVOV = numpy.asarray(eris.OVOV) OVOV = OVOV - OVOV.transpose(0,3,2,1) m3bb += lib.einsum('kalb,ijkl->ijab', OVOV, tmp) * .25 tmp = lib.einsum('iJcD,kLcD->iJkL', l2ab, tauab) ovOV = numpy.asarray(eris.ovOV) m3ab += lib.einsum('kaLB,iJkL->iJaB', ovOV, tmp) * .5 tmp = lib.einsum('iJdC,lKdC->iJKl', l2ab, tauab) m3ab += lib.einsum('laKB,iJKl->iJaB', ovOV, tmp) * .5 u1a += numpy.einsum('ijab,jb->ia', m3aa, t1a) u1a += numpy.einsum('iJaB,JB->ia', m3ab, t1b) u1b += numpy.einsum('IJAB,JB->IA', m3bb, t1b) u1b += numpy.einsum('jIbA,jb->IA', m3ab, t1a) u2aa += m3aa u2bb += m3bb u2ab += m3ab u2aa += ovov.transpose(0,2,1,3) u2bb += OVOV.transpose(0,2,1,3) u2ab += ovOV.transpose(0,2,1,3) fov1 = fova + numpy.einsum('kcjb,kc->jb', ovov, t1a) fov1+= numpy.einsum('jbKC,KC->jb', ovOV, t1b) tmp = numpy.einsum('ia,jb->ijab', l1a, fov1) tmp+= einsum('kica,jbck->ijab', l2aa, imds.wovvo) tmp+= einsum('iKaC,jbCK->ijab', l2ab, imds.wovVO) tmp = tmp - tmp.transpose(1,0,2,3) u2aa += tmp - tmp.transpose(0,1,3,2) fov1 = fovb + numpy.einsum('kcjb,kc->jb', OVOV, t1b) fov1+= numpy.einsum('kcJB,kc->JB', ovOV, t1a) tmp = numpy.einsum('ia,jb->ijab', l1b, fov1) tmp+= einsum('kica,jbck->ijab', l2bb, imds.wOVVO) tmp+= einsum('kIcA,JBck->IJAB', l2ab, imds.wOVvo) tmp = tmp - tmp.transpose(1,0,2,3) u2bb += tmp - tmp.transpose(0,1,3,2) fov1 = fovb + numpy.einsum('kcjb,kc->jb', OVOV, t1b) fov1+= numpy.einsum('kcJB,kc->JB', ovOV, t1a) u2ab += numpy.einsum('ia,JB->iJaB', l1a, fov1) u2ab += einsum('iKaC,JBCK->iJaB', l2ab, imds.wOVVO) u2ab += einsum('kica,JBck->iJaB', l2aa, imds.wOVvo) u2ab += einsum('kIaC,jBCk->jIaB', l2ab, imds.woVVo) u2ab += einsum('iKcA,JbcK->iJbA', l2ab, imds.wOvvO) fov1 = fova + numpy.einsum('kcjb,kc->jb', ovov, t1a) fov1+= numpy.einsum('jbKC,KC->jb', ovOV, t1b) u2ab += numpy.einsum('ia,jb->jiba', l1b, fov1) u2ab += einsum('kIcA,jbck->jIbA', l2ab, imds.wovvo) u2ab += einsum('KICA,jbCK->jIbA', l2bb, imds.wovVO) ovoo = numpy.asarray(eris.ovoo) ovoo = ovoo - ovoo.transpose(2,1,0,3) OVOO = numpy.asarray(eris.OVOO) OVOO = OVOO - OVOO.transpose(2,1,0,3) OVoo = numpy.asarray(eris.OVoo) ovOO = numpy.asarray(eris.ovOO) tmp = einsum('ka,jbik->ijab', l1a, ovoo) tmp+= einsum('ijca,cb->ijab', l2aa, v1a) tmp+= einsum('ca,icjb->ijab', mvv1, ovov) u2aa -= tmp - tmp.transpose(0,1,3,2) tmp = einsum('ka,jbik->ijab', l1b, OVOO) tmp+= einsum('ijca,cb->ijab', l2bb, v1b) tmp+= einsum('ca,icjb->ijab', mVV1, OVOV) u2bb -= tmp - tmp.transpose(0,1,3,2) u2ab -= einsum('ka,JBik->iJaB', l1a, OVoo) u2ab += einsum('iJaC,CB->iJaB', l2ab, v1b) u2ab -= einsum('ca,icJB->iJaB', mvv1, ovOV) u2ab -= einsum('KA,ibJK->iJbA', l1b, ovOO) u2ab += einsum('iJcA,cb->iJbA', l2ab, v1a) u2ab -= einsum('CA,ibJC->iJbA', mVV1, ovOV) u1a += fova u1b += fovb u1a += einsum('ib,ba->ia', l1a, v1a) u1a -= einsum('ja,ij->ia', l1a, v2a) u1b += einsum('ib,ba->ia', l1b, v1b) u1b -= einsum('ja,ij->ia', l1b, v2b) u1a += numpy.einsum('jb,iabj->ia', l1a, eris.ovvo) u1a -= numpy.einsum('jb,ijba->ia', l1a, eris.oovv) u1a += numpy.einsum('JB,iaBJ->ia', l1b, eris.ovVO) u1b += numpy.einsum('jb,iabj->ia', l1b, eris.OVVO) u1b -= numpy.einsum('jb,ijba->ia', l1b, eris.OOVV) u1b += numpy.einsum('jb,iabj->ia', l1a, eris.OVvo) u1a -= einsum('kjca,ijck->ia', l2aa, imds.woovo) u1a -= einsum('jKaC,ijCK->ia', l2ab, imds.wooVO) u1b -= einsum('kjca,ijck->ia', l2bb, imds.wOOVO) u1b -= einsum('kJcA,IJck->IA', l2ab, imds.wOOvo) u1a -= einsum('ikbc,back->ia', l2aa, imds.wvvvo) u1a -= einsum('iKbC,baCK->ia', l2ab, imds.wvvVO) u1b -= einsum('IKBC,BACK->IA', l2bb, imds.wVVVO) u1b -= einsum('kIcB,BAck->IA', l2ab, imds.wVVvo) u1a += numpy.einsum('jiba,bj->ia', l2aa, imds.w3a) u1a += numpy.einsum('iJaB,BJ->ia', l2ab, imds.w3b) u1b += numpy.einsum('JIBA,BJ->IA', l2bb, imds.w3b) u1b += numpy.einsum('jIbA,bj->IA', l2ab, imds.w3a) tmpa = t1a + numpy.einsum('kc,kjcb->jb', l1a, t2aa) tmpa += numpy.einsum('KC,jKbC->jb', l1b, t2ab) tmpa -= einsum('bd,jd->jb', mvv1, t1a) tmpa -= einsum('lj,lb->jb', moo, t1a) tmpb = t1b + numpy.einsum('kc,kjcb->jb', l1b, t2bb) tmpb += numpy.einsum('kc,kJcB->JB', l1a, t2ab) tmpb -= einsum('bd,jd->jb', mVV1, t1b) tmpb -= einsum('lj,lb->jb', mOO, t1b) u1a += numpy.einsum('jbia,jb->ia', ovov, tmpa) u1a += numpy.einsum('iaJB,JB->ia', ovOV, tmpb) u1b += numpy.einsum('jbia,jb->ia', OVOV, tmpb) u1b += numpy.einsum('jbIA,jb->IA', ovOV, tmpa) u1a -= numpy.einsum('iajk,kj->ia', ovoo, moo1) u1a -= numpy.einsum('iaJK,KJ->ia', ovOO, mOO1) u1b -= numpy.einsum('iajk,kj->ia', OVOO, mOO1) u1b -= numpy.einsum('IAjk,kj->IA', OVoo, moo1) tmp = fova - numpy.einsum('kbja,jb->ka', ovov, t1a) tmp += numpy.einsum('kaJB,JB->ka', ovOV, t1b) u1a -= lib.einsum('ik,ka->ia', moo, tmp) u1a -= lib.einsum('ca,ic->ia', mvv, tmp) tmp = fovb - numpy.einsum('kbja,jb->ka', OVOV, t1b) tmp += numpy.einsum('jbKA,jb->KA', ovOV, t1a) u1b -= lib.einsum('ik,ka->ia', mOO, tmp) u1b -= lib.einsum('ca,ic->ia', mVV, tmp) eia = lib.direct_sum('i-j->ij', mo_ea_o, mo_ea_v) eIA = lib.direct_sum('i-j->ij', mo_eb_o, mo_eb_v) u1a /= eia u1b /= eIA u2aa /= lib.direct_sum('ia+jb->ijab', eia, eia) u2ab /= lib.direct_sum('ia+jb->ijab', eia, eIA) u2bb /= lib.direct_sum('ia+jb->ijab', eIA, eIA) time0 = log.timer_debug1('update l1 l2', *time0) return (u1a,u1b), (u2aa,u2ab,u2bb)
if __name__ == '__main__': from pyscf import gto from pyscf import scf from pyscf.cc import gccsd mol = gto.Mole() mol.atom = [ [8 , (0. , 0. , 0.)], [1 , (0. , -0.757 , 0.587)], [1 , (0. , 0.757 , 0.587)]] mol.basis = '631g' mol.spin = 2 mol.build() mf = scf.UHF(mol).run() mycc = gccsd.GCCSD(scf.addons.convert_to_ghf(mf)) eris = mycc.ao2mo() mycc.kernel() l1, l2 = mycc.solve_lambda(mycc.t1, mycc.t2, eris=eris) l1ref = mycc.spin2spatial(l1, mycc.mo_coeff.orbspin) l2ref = mycc.spin2spatial(l2, mycc.mo_coeff.orbspin) mycc = uccsd.UCCSD(mf) eris = mycc.ao2mo() mycc.kernel() conv, l1, l2 = kernel(mycc, eris, mycc.t1, mycc.t2, tol=1e-8) print(abs(l1[0]-l1ref[0]).max()) print(abs(l1[1]-l1ref[1]).max()) print(abs(l2[0]-l2ref[0]).max()) print(abs(l2[1]-l2ref[1]).max()) print(abs(l2[2]-l2ref[2]).max())