1#!/usr/bin/env python 2# Copyright 2014-2021 The PySCF Developers. All Rights Reserved. 3# 4# Licensed under the Apache License, Version 2.0 (the "License"); 5# you may not use this file except in compliance with the License. 6# You may obtain a copy of the License at 7# 8# http://www.apache.org/licenses/LICENSE-2.0 9# 10# Unless required by applicable law or agreed to in writing, software 11# distributed under the License is distributed on an "AS IS" BASIS, 12# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13# See the License for the specific language governing permissions and 14# limitations under the License. 15 16import tempfile 17import h5py 18import numpy as np 19from pyscf import lib 20 21#einsum = np.einsum 22einsum = lib.einsum 23 24# Ref: Gauss and Stanton, J. Chem. Phys. 103, 3561 (1995) Table III 25 26# Section (a) 27 28def make_tau(t2, t1a, t1b, fac=1, out=None): 29 tmp = einsum('ia,jb->ijab',t1a,t1b) 30 t1t1 = tmp - tmp.transpose(1,0,2,3) - tmp.transpose(0,1,3,2) + tmp.transpose(1,0,3,2) 31 tau1 = t2 + fac*0.50*t1t1 32 return tau1 33 34def cc_Fvv(t1,t2,eris): 35 nocc, nvir = t1.shape 36 fov = eris.fock[:nocc,nocc:] 37 fvv = eris.fock[nocc:,nocc:] 38 eris_vovv = np.array(eris.ovvv).transpose(1,0,3,2) 39 tau_tilde = make_tau(t2,t1,t1,fac=0.5) 40 Fae = ( fvv - 0.5*einsum('me,ma->ae',fov,t1) 41 + einsum('mf,amef->ae',t1,eris_vovv) 42 - 0.5*einsum('mnaf,mnef->ae',tau_tilde,eris.oovv) ) 43 return Fae 44 45def cc_Foo(t1,t2,eris): 46 nocc, nvir = t1.shape 47 fov = eris.fock[:nocc,nocc:] 48 foo = eris.fock[:nocc,:nocc] 49 tau_tilde = make_tau(t2,t1,t1,fac=0.5) 50 Fmi = (foo + 0.5*einsum('me,ie->mi',fov,t1) 51 + einsum('ne,mnie->mi',t1,eris.ooov) 52 + 0.5*einsum('inef,mnef->mi',tau_tilde,eris.oovv) ) 53 return Fmi 54 55def cc_Fov(t1,t2,eris): 56 nocc, nvir = t1.shape 57 fov = eris.fock[:nocc,nocc:] 58 Fme = fov + einsum('nf,mnef->me',t1,eris.oovv) 59 return Fme 60 61def cc_Woooo(t1,t2,eris): 62 tau = make_tau(t2,t1,t1) 63 tmp = einsum('je,mnie->mnij',t1,eris.ooov) 64 Wmnij = eris.oooo + tmp - tmp.transpose(0,1,3,2) 65 Wmnij += 0.25*einsum('ijef,mnef->mnij',tau,eris.oovv) 66 return Wmnij 67 68def cc_Wvvvv(t1,t2,eris): 69 tau = make_tau(t2,t1,t1) 70 #eris_vovv = np.array(eris.ovvv).transpose(1,0,3,2) 71 #tmp = einsum('mb,amef->abef',t1,eris_vovv) 72 #Wabef = eris.vvvv - tmp + tmp.transpose(1,0,2,3) 73 #Wabef += 0.25*einsum('mnab,mnef->abef',tau,eris.oovv) 74 if t1.dtype == np.complex128: ds_type = 'c16' 75 else: ds_type = 'f8' 76 _tmpfile1 = tempfile.NamedTemporaryFile(dir=lib.param.TMPDIR) 77 fimd = h5py.File(_tmpfile1.name) 78 nocc, nvir = t1.shape 79 Wabef = fimd.create_dataset('vvvv', (nvir,nvir,nvir,nvir), ds_type) 80 for a in range(nvir): 81 Wabef[a] = eris.vvvv[a] 82 Wabef[a] -= einsum('mb,mfe->bef',t1,eris.ovvv[:,a,:,:]) 83 Wabef[a] += einsum('m,mbfe->bef',t1[:,a],eris.ovvv) 84 Wabef[a] += 0.25*einsum('mnb,mnef->bef',tau[:,:,a,:],eris.oovv) 85 return Wabef 86 87def cc_Wovvo(t1,t2,eris): 88 eris_ovvo = -np.array(eris.ovov).transpose(0,1,3,2) 89 eris_oovo = -np.array(eris.ooov).transpose(0,1,3,2) 90 Wmbej = eris_ovvo.copy() 91 Wmbej += einsum('jf,mbef->mbej',t1,eris.ovvv) 92 Wmbej += -einsum('nb,mnej->mbej',t1,eris_oovo) 93 Wmbej += -0.5*einsum('jnfb,mnef->mbej',t2,eris.oovv) 94 Wmbej += -einsum('jf,nb,mnef->mbej',t1,t1,eris.oovv) 95 return Wmbej 96 97### Section (b) 98 99def Fvv(t1,t2,eris): 100 ccFov = cc_Fov(t1,t2,eris) 101 Fae = cc_Fvv(t1,t2,eris) - 0.5*einsum('ma,me->ae',t1,ccFov) 102 return Fae 103 104def Foo(t1,t2,eris): 105 ccFov = cc_Fov(t1,t2,eris) 106 Fmi = cc_Foo(t1,t2,eris) + 0.5*einsum('ie,me->mi',t1,ccFov) 107 return Fmi 108 109def Fov(t1,t2,eris): 110 Fme = cc_Fov(t1,t2,eris) 111 return Fme 112 113def Woooo(t1,t2,eris): 114 tau = make_tau(t2,t1,t1) 115 Wmnij = cc_Woooo(t1,t2,eris) + 0.25*einsum('ijef,mnef->mnij',tau,eris.oovv) 116 return Wmnij 117 118def Wvvvv(t1,t2,eris): 119 tau = make_tau(t2,t1,t1) 120 #Wabef = cc_Wvvvv(t1,t2,eris) + 0.25*einsum('mnab,mnef->abef',tau,eris.oovv) 121 if t1.dtype == np.complex128: ds_type = 'c16' 122 else: ds_type = 'f8' 123 _tmpfile1 = tempfile.NamedTemporaryFile(dir=lib.param.TMPDIR) 124 fimd = h5py.File(_tmpfile1.name) 125 nocc, nvir = t1.shape 126 Wabef = fimd.create_dataset('vvvv', (nvir,nvir,nvir,nvir), ds_type) 127 #_cc_Wvvvv = cc_Wvvvv(t1,t2,eris) 128 for a in range(nvir): 129 #Wabef[a] = _cc_Wvvvv[a] 130 Wabef[a] = eris.vvvv[a] 131 Wabef[a] -= einsum('mb,mfe->bef',t1,eris.ovvv[:,a,:,:]) 132 Wabef[a] += einsum('m,mbfe->bef',t1[:,a],eris.ovvv) 133 #Wabef[a] += 0.25*einsum('mnb,mnef->bef',tau[:,:,a,:],eris.oovv) 134 135 #Wabef[a] += 0.25*einsum('mnb,mnef->bef',tau[:,:,a,:],eris.oovv) 136 Wabef[a] += 0.5*einsum('mnb,mnef->bef',tau[:,:,a,:],eris.oovv) 137 return Wabef 138 139def Wovvo(t1,t2,eris): 140 Wmbej = cc_Wovvo(t1,t2,eris) - 0.5*einsum('jnfb,mnef->mbej',t2,eris.oovv) 141 return Wmbej 142 143def Wooov(t1,t2,eris): 144 Wmnie = eris.ooov + einsum('if,mnfe->mnie',t1,eris.oovv) 145 return Wmnie 146 147def Wvovv(t1,t2,eris): 148 eris_vovv = -np.array(eris.ovvv).transpose(1,0,2,3) 149 Wamef = eris_vovv - einsum('na,nmef->amef',t1,eris.oovv) 150 return Wamef 151 152def Wovoo(t1,t2,eris): 153 eris_ovvo = -np.array(eris.ovov).transpose(0,1,3,2) 154 tmp1 = einsum('mnie,jnbe->mbij',eris.ooov,t2) 155 tmp2 = (einsum('ie,mbej->mbij',t1,eris_ovvo) 156 - einsum('ie,njbf,mnef->mbij',t1,t2,eris.oovv)) 157 FFov = Fov(t1,t2,eris) 158 WWoooo = Woooo(t1,t2,eris) 159 tau = make_tau(t2,t1,t1) 160 Wmbij = (eris.ovoo - einsum('me,ijbe->mbij',FFov,t2) 161 - einsum('nb,mnij->mbij',t1,WWoooo) 162 + 0.5 * einsum('mbef,ijef->mbij',eris.ovvv,tau) 163 + tmp1 - tmp1.transpose(0,1,3,2) 164 + tmp2 - tmp2.transpose(0,1,3,2) ) 165 return Wmbij 166 167def Wvvvo(t1,t2,eris,_Wvvvv=None): 168 eris_ovvo = -np.array(eris.ovov).transpose(0,1,3,2) 169 eris_vvvo = -np.array(eris.ovvv).transpose(2,3,1,0).conj() 170 eris_oovo = -np.array(eris.ooov).transpose(0,1,3,2) 171 tmp1 = einsum('mbef,miaf->abei',eris.ovvv,t2) 172 tmp2 = (einsum('ma,mbei->abei',t1,eris_ovvo) 173 - einsum('ma,nibf,mnef->abei',t1,t2,eris.oovv)) 174 FFov = Fov(t1,t2,eris) 175 tau = make_tau(t2,t1,t1) 176 Wabei = eris_vvvo 177 Wabei += -einsum('me,miab->abei',FFov,t2) 178 Wabei += 0.5 * einsum('mnei,mnab->abei',eris_oovo,tau) 179 Wabei += -tmp1 + tmp1.transpose(1,0,2,3) 180 Wabei += -tmp2 + tmp2.transpose(1,0,2,3) 181 nocc,nvir = t1.shape 182 if _Wvvvv is None: 183 _Wvvvv = Wvvvv(t1,t2,eris) 184 for a in range(nvir): 185 Wabei[a] += einsum('if,bef->bei',t1,_Wvvvv[a]) 186 return Wabei 187 188def Wvvvo_incore(t1,t2,eris): 189 eris_ovvo = -np.array(eris.ovov).transpose(0,1,3,2) 190 eris_vvvo = -np.array(eris.ovvv).transpose(2,3,1,0).conj() 191 eris_oovo = -np.array(eris.ooov).transpose(0,1,3,2) 192 tmp1 = einsum('mbef,miaf->abei',eris.ovvv,t2) 193 tmp2 = (einsum('ma,mbei->abei',t1,eris_ovvo) 194 - einsum('ma,nibf,mnef->abei',t1,t2,eris.oovv)) 195 FFov = Fov(t1,t2,eris) 196 WWvvvv = Wvvvv(t1,t2,eris) 197 tau = make_tau(t2,t1,t1) 198 Wabei = (eris_vvvo - einsum('me,miab->abei',FFov,t2) 199 + einsum('if,abef->abei',t1,WWvvvv) 200 + 0.5 * einsum('mnei,mnab->abei',eris_oovo,tau) 201 - tmp1 + tmp1.transpose(1,0,2,3) 202 - tmp2 + tmp2.transpose(1,0,2,3) ) 203 return Wabei 204