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