mp2grad_uhf_sv1
SIAL MBPT2_GRAD_AO2
#
# This SIAL program computes the second-order correction to the density
# and the second-order correction to the energy-weighted density matrix.
# They are written to file to be used in the final gradient construction.
# Author: Victor Lotrich
# Copyright: University of Florida
# This program is licensed under the GPLv2
#
# BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
#
# Declare indeces
# ---------------
#
index kiter = 1, 100
#
aoindex mu = 1, norb
aoindex nu = 1, norb
aoindex lambda = 1, norb
aoindex sigma = 1, norb
#
moaindex i = baocc, eaocc
moaindex i1= baocc, eaocc
moaindex i2= baocc, eaocc
moaindex i3= baocc, eaocc
#
moaindex a = bavirt, eavirt
moaindex a1= bavirt, eavirt
moaindex a2= bavirt, eavirt
moaindex a3= bavirt, eavirt
#
mobindex j = bbocc, ebocc
mobindex j1= bbocc, ebocc
mobindex j2= bbocc, ebocc
mobindex j3= bbocc, ebocc
#
mobindex b = bbvirt, ebvirt
mobindex b1= bbvirt, ebvirt
mobindex b2= bbvirt, ebvirt
mobindex b3= bbvirt, ebvirt
#
moaindex p = baocc, eavirt
moaindex p1= baocc, eavirt
moaindex p2= baocc, eavirt
moaindex p3= baocc, eavirt
#
mobindex q = bbocc, ebvirt
mobindex q1= bbocc, ebvirt
mobindex q2= bbocc, ebvirt
mobindex q3= bbocc, ebvirt
#
# Arrays used in transformation for AO2 algorithm
# -----------------------------------------------
#
temp Txixi(mu,i1,lambda,i)
temp T1xixi(mu,i1,lambda,i)
temp Txxii(mu,nu,i1,i)
temp Tixxi(i1,nu,lambda,i)
temp Txipi(mu,i,p,i1)
temp Tpipi(p1,i,p,i1)
temp T1pipi(p1,i,p,i1)
temp T2pipi(p1,i,p,i1)
temp Tixai(i,mu,a,i1)
temp Txaii(mu,a,i,i1)
temp Tiaai(i,a1,a,i1)
temp Taaii(a,a1,i,i1)
temp Txaai(mu,a1,a,i)
temp Taaai(a2,a,a1,i)
temp T1aaai(a2,a,a1,i)
temp Txxai(mu,nu,a,i)
#
served Vxixi(mu,i1,lambda,i)
served Vxxii(mu,nu,i1,i)
served Viixx(i1,i,mu,nu)
served Vixxi(i1,nu,lambda,i)
served Vxipi(mu,i,p,i1)
served Vixai(i,mu,a,i1)
served Vxaii(mu,a,i,i1)
#
# Permanent
# ---------
served Viaai(i,a1,a,i1)
served Vaaii(a,a1,i,i1)
served VSpipi(p1,i,p,i1)
served ASpipi(p1,i,p,i1)
#
temp Txjxj(mu,j1,lambda,j)
temp T1xjxj(mu,j1,lambda,j)
temp Txxjj(mu,nu,j1,j)
temp Tjjxx(j1,nu,lambda,j)
temp Tjxxj(j1,nu,lambda,j)
temp Txjqj(mu,j,q,j1)
temp Tqjqj(q1,j,q,j1)
temp T1qjqj(q1,j,q,j1)
temp T2qjqj(q1,j,q,j1)
temp Tjxbj(j,mu,b,j1)
temp Txbjj(mu,b,j,j1)
temp Tjbbj(j,b1,b,j1)
temp Tbbjj(b,b1,j,j1)
temp Txbbj(mu,b1,b,j)
temp Tbbbj(b2,b,b1,j)
temp T1bbbj(b2,b,b1,j)
temp Txxbj(mu,nu,b,j)
#
served Vxjxj(mu,j1,lambda,j)
served Vxxjj(mu,nu,j1,j)
served Vjjxx(j1,nu,lambda,j)
served Vjxxj(j1,nu,lambda,j)
served Vxjqj(mu,j,q,j1)
served Vjxbj(j,mu,b,j1)
served Vxbjj(mu,b,j,j1)
#
# Permanent
# ---------
served Vjbbj(j,b1,b,j1)
served Vbbjj(b,b1,j,j1)
served VSqjqj(q1,j,q,j1)
served ASqjqj(q1,j,q,j1)
#
temp Txixj(mu,i,nu,j)
temp Txiqj(mu,i,q,j)
temp Tpiqj(p,i,q,j)
temp T1piqj(p,i,q,j)
temp T2piqj(p,i,q,j)
temp Tiixx(i,i1,mu,nu)
temp Tiixb(i,i1,mu,b)
temp Tiibb(i,i1,b1,b)
temp Txajj(mu,a,j,j1)
temp Taajj(a,a1,j,j1)
temp Txabj(mu,a,b,j)
temp Tixxj(i,mu,nu,j)
temp Tixbj(i,mu,b,j)
temp Tiabj(i,a,b,j)
temp Taabj(a,a1,b,j)
#
served Vxixj(mu,i,nu,j)
served Vxiqj(mu,i,q,j)
served Vpiqj(p,i,q,j)
served Apiqj(p,i,q,j)
served Viixb(i,i1,mu,b)
served Viibb(i,i1,b1,b)
served Vxajj(mu,a,j,j1)
served Vixxj(i,mu,nu,j)
served Vixbj(i,mu,b,j)
#
# Permanent
# ---------
served Viabj(i,a,b,j)
served Vaajj(a,a1,j,j1)
#
temp Txbii(mu,b,i,i1)
temp Tbbii(b,b1,i,i1)
temp Tjbii(j,b,i,i1)
temp Txbai(mu,b,a,i)
temp Tbbai(b,b1,a,i)
#
served Vxbii(mu,b,i,i1)
#
# Permanent
# ---------
served Vbbii(b,b1,i,i1)
#
# End Arrays used in transformation for AO2 algorithm
# ---------------------------------------------------
#
# Declare one-particle density arrays
# -----------------------------------
#
distributed Pij_aa(i,i1)
distributed Pij_bb(j,j1)
distributed Pab_aa(a,a1)
distributed Pab_bb(b,b1)
distributed Lai_aa(a,i)
distributed Lai_bb(b,j)
distributed Painew_aa(a,i)
distributed Paiold_aa(a,i)
distributed Painew_bb(b,j)
distributed Paiold_bb(b,j)
distributed Wab_aa(a,a1)
distributed Wab_bb(b,b1)
distributed Wij_aa(i,i1)
distributed Wij_bb(j,j1)
distributed Wai_aa(a,i)
distributed Wai_bb(b,j)
distributed P2_ao(mu,nu)
distributed P2A_ao(mu,nu)
distributed P2B_ao(mu,nu)
distributed W2_ao(mu,nu)
distributed Paa_ao(mu,nu)
distributed Pbb_ao(mu,nu)
distributed Whfa(mu,nu)
distributed Whfb(mu,nu)
distributed Dhfa(mu,nu)
distributed Dhfb(mu,nu)
distributed Dhf(mu,nu)
distributed Lxi(mu,i)
distributed Lxj(nu,j)
distributed Yxi(mu,i)
distributed Yxj(nu,j)
local L0xxxi(mu,nu,lambda,i)
local L0xxxj(mu,nu,lambda,j)
temp V0xxxi(mu,nu,lambda,i)
temp V0xxxj(mu,nu,lambda,j)
local LDhfa(mu,nu)
local LDhfb(mu,nu)
local LP2A_ao(mu,nu)
local LP2B_ao(mu,nu)
local La(p,mu)
local Lb(q,mu)
#
# Declare temporary arrays
# ------------------------
#
temp Txxxx(mu,nu,lambda,sigma)
temp TSxxxx(mu,nu,lambda,sigma)
temp T1xxxx(mu,nu,lambda,sigma)
temp Txxxi(mu,nu,lambda,i)
temp T1xxxi(mu,nu,lambda,i)
temp T2xxxi(mu,nu,lambda,i)
temp T3xxxi(mu,nu,lambda,i)
temp T4xxxi(mu,nu,lambda,i)
temp Txxxj(mu,nu,lambda,j)
temp T1xxxj(mu,nu,lambda,j)
temp T2xxxj(mu,nu,lambda,j)
temp T3xxxj(mu,nu,lambda,j)
temp T4xxxj(mu,nu,lambda,j)
temp Txixx(mu,i,nu,lambda)
temp T1xixx(mu,i,nu,lambda)
temp txxix(nu,mu,i,lambda)
temp txxjx(nu,mu,j,lambda)
temp Txjxx(mu,j,nu,lambda)
temp T1xjxx(mu,j,nu,lambda)
temp gaa(mu,i,nu,lambda)
temp gab(mu,i,nu,lambda)
temp gbb(mu,i,nu,lambda)
temp TAxxxi(mu,lambda,nu,i1)
temp TBxxxi(mu,lambda,sigma,i1)
temp TAxxxj(mu,lambda,nu,j1)
temp TBxxxj(mu,lambda,sigma,j1)
temp Tii(i,i1)
temp T1ii(i,i1)
temp Tjj(j,j1)
temp T1jj(j,j1)
temp Taa(a,a1)
temp Tbb(b,b1)
temp Tai(a,i)
temp T1ai(a,i)
temp Tbj(b,j)
temp T1bj(b,j)
temp Tiiaa(i,i1,a,a1)
temp T1iiaa(i,i1,a,a1)
temp Tjjbb(j,j1,b,b1)
temp T1jjbb(j,j1,b,b1)
temp Txi(mu,i)
temp Ixi(mu,i)
temp I1xi(mu,i)
temp Jxi(mu,i)
temp Kxi(mu,i)
temp Txj(mu,j)
temp Ixj(mu,j)
temp I1xj(mu,j)
temp Jxj(mu,j)
temp Kxj(mu,j)
temp Ixx(mu,nu)
temp I1xx(mu,nu)
temp Jxx(mu,nu)
temp J1xx(mu,nu)
temp Kxx(mu,nu)
temp K1xx(mu,nu)
temp Ixa(mu,a)
temp Jxa(mu,a)
temp Ixb(mu,b)
temp Jxb(mu,b)
temp Kxa(mu,a)
temp Kxb(mu,b)
temp Tpq(mu,nu)
temp Txx(mu,nu)
temp T1xx(mu,nu)
temp T2xx(mu,nu)
temp T3xx(mu,nu)
temp T4xx(mu,nu)
temp T5xx(mu,nu)
temp T6xx(mu,nu)
temp T7xx(mu,nu)
temp T8xx(mu,nu)
temp T9xx(mu,nu)
#
# Declare served arrays
# ---------------------
#
temp AOINT(mu,nu,lambda,sigma)
temp dx1(mu,nu,lambda,sigma)
temp dy1(mu,nu,lambda,sigma)
temp dz1(mu,nu,lambda,sigma)
temp dx2(mu,nu,lambda,sigma)
temp dy2(mu,nu,lambda,sigma)
temp dz2(mu,nu,lambda,sigma)
temp dx3(mu,nu,lambda,sigma)
temp dy3(mu,nu,lambda,sigma)
temp dz3(mu,nu,lambda,sigma)
temp dx4(mu,nu,lambda,sigma)
temp dy4(mu,nu,lambda,sigma)
temp dz4(mu,nu,lambda,sigma)
#
# Declare Local arrays
# --------------------
#
local Xa(lambda,i,sigma,mu)
local Xb(lambda,j,sigma,mu)
local xTa(sigma,mu,lambda,i)
local xTb(sigma,mu,lambda,j)
temp D2(mu,lambda,nu,sigma)
temp D3(mu,lambda,nu,sigma)
#
local Lxixi(mu,i1,lambda,i)
local Lxixj(mu,i,lambda,j)
local Lxjxj(mu,j1,lambda,j)
local Lxxii(mu,nu,i1,i)
local Lxxjj(mu,nu,j1,j)
local Lixxi(i1,nu,lambda,i)
local Lixxj(i,nu,lambda,j)
local Ljxxj(j1,nu,lambda,j)
local Lxipi(mu,i,p,i1)
local Lpipi(p1,i,p,i1)
local Lxaii(mu,a,i,i1)
local Laaii(a1,a,i,i1)
local Lixai(i,mu,a,i1)
local Liaai(i,a1,a,i1)
#
local Lxjqj(mu,j,q,j1)
local Lqjqj(q1,j,q,j1)
local Lxbjj(mu,b,j,j1)
local Lbbjj(b1,b,j,j1)
local Ljxbj(j,mu,b,j1)
local Ljbbj(j,b1,b,j1)
local Lxbii(mu,b,i,i1)
local Lbbii(b1,b,i,i1)
local Lxajj(mu,a,j,j1)
local Laajj(a1,a,j,j1)
local Lixbj(i,mu,b,j)
local Liabj(i,a,b,j)
local Liixb(i,i1,mu,b)
local Liibb(i,i1,b1,b)
local Lxiqj(mu,i,q,j)
local Lpiqj(p,i,q,j)
#
local Lxiai(mu,i,a1,i1)
local Lxjbj(mu,j,b1,j1)
local Lxibj(mu,i,b,j)
local L1xixi(mu,i,nu,i1)
local L1xjxj(mu,j,nu,j1)
local L1xixj(mu,i,nu,j)
#
# Arrays and scalars used in iterative computation of Dai
# --------------------------------------------------------
#
distributed D0ai(a,i)
distributed D1ai(a,i)
distributed D2ai(a,i)
distributed D3ai(a,i)
distributed D4ai(a,i)
#
distributed D0bj(b,j)
distributed D1bj(b,j)
distributed D2bj(b,j)
distributed D3bj(b,j)
distributed D4bj(b,j)
#
distributed e1ai(a,i)
distributed e2ai(a,i)
distributed e3ai(a,i)
distributed e4ai(a,i)
distributed e5ai(a,i)
#
distributed e1bj(b,j)
distributed e2bj(b,j)
distributed e3bj(b,j)
distributed e4bj(b,j)
distributed e5bj(b,j)
#
scalar b11
scalar b12
scalar b13
scalar b14
scalar b15
scalar b16
scalar b17
scalar b18
scalar b19
scalar b110
#
scalar b22
scalar b23
scalar b24
scalar b25
scalar b26
scalar b27
scalar b28
scalar b29
scalar b210
#
scalar b33
scalar b34
scalar b35
scalar b36
scalar b37
scalar b38
scalar b39
scalar b310
scalar b44
scalar b45
scalar b46
scalar b47
scalar b48
scalar b49
scalar b410
#
scalar b55
scalar b56
scalar b57
scalar b58
scalar b59
scalar b510
#
scalar b66
scalar b67
scalar b68
scalar b69
scalar b610
#
scalar b77
scalar b78
scalar b79
scalar b710
#
scalar b88
scalar b89
scalar b810
#
scalar b99
scalar b910
#
scalar b1010
#
scalar Tb11
scalar Tb12
scalar Tb13
scalar Tb14
scalar Tb15
scalar Tb16
scalar Tb17
scalar Tb18
scalar Tb19
scalar Tb110
#
scalar Tb22
scalar Tb23
scalar Tb24
scalar Tb25
scalar Tb26
scalar Tb27
scalar Tb28
scalar Tb29
scalar Tb210
#
scalar Tb33
scalar Tb34
scalar Tb35
scalar Tb36
scalar Tb37
scalar Tb38
scalar Tb39
scalar Tb310
scalar Tb44
scalar Tb45
scalar Tb46
scalar Tb47
scalar Tb48
scalar Tb49
scalar Tb410
#
scalar Tb55
scalar Tb56
scalar Tb57
scalar Tb58
scalar Tb59
scalar Tb510
#
scalar Tb66
scalar Tb67
scalar Tb68
scalar Tb69
scalar Tb610
#
scalar Tb77
scalar Tb78
scalar Tb79
scalar Tb710
#
scalar Tb88
scalar Tb89
scalar Tb810
#
scalar Tb99
scalar Tb910
#
scalar Tb1010
#
scalar c1
scalar c2
scalar c3
scalar c4
scalar c5
scalar c6
scalar c7
scalar c8
scalar c9
scalar c10
#
# Declare scalars
# ---------------
#
scalar etemp
scalar esum
scalar esumaa
scalar esumbb
scalar esumab
scalar ecorraa
scalar ecorrbb
scalar ecorrab
scalar ecorrT
scalar enew
scalar eold
scalar ecrit
scalar ediff
scalar mp2_energy
#
# BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
#
# ---------------------------------------------------------------------------------
#
PROC ENERGY
#
esumaa = 0.0
ecorraa = 0.0
PARDO a, a1, i, i1
#
REQUEST VSpipi(a,i,a1,i1) a
Tpipi(a,i,a1,i1) = VSpipi(a,i,a1,i1)
execute energy_denominator Tpipi
#
etemp = Tpipi(a,i,a1,i1)*VSpipi(a,i,a1,i1)
etemp = 0.25*etemp
esumaa += etemp
ENDPARDO a, a1, i, i1
#
ecorrbb = 0.0
esumbb = 0.0
PARDO b, b1, j, j1
#
REQUEST VSqjqj(b,j,b1,j1) b
Tqjqj(b,j,b1,j1) = VSqjqj(b,j,b1,j1)
execute energy_denominator Tqjqj
#
etemp = Tqjqj(b,j,b1,j1)*VSqjqj(b,j,b1,j1)
etemp = 0.25*etemp
esumbb += etemp
#
ENDPARDO b, b1, j, j1
ecorrab = 0.0
esumab = 0.0
PARDO a, b, i, j
#
REQUEST Vpiqj(a,i,b,j) a
Tpiqj(a,i,b,j) = Vpiqj(a,i,b,j)
execute energy_denominator Tpiqj
#
etemp = Tpiqj(a,i,b,j)*Vpiqj(a,i,b,j)
esumab += etemp
#
ENDPARDO a, b, i, j
execute sip_barrier
collective ecorraa += esumaa
execute print_scalar ecorraa
#
collective ecorrbb += esumbb
execute print_scalar ecorrbb
#
collective ecorrab += esumab
execute print_scalar ecorrab
#
mp2_energy = ecorraa
mp2_energy += ecorrbb
mp2_energy += ecorrab
execute print_scalar mp2_energy
#
totenerg = mp2_energy
totenerg += scfeneg
execute print_scalar totenerg
#
ENDPROC ENERGY
#
# ------------------------------------------------------------------------
#
#
# ------------------------------------------------------------------------
#
PROC TRAN_XXII
#
# ------------------------------------------------------------------------
#
execute sip_barrier
#
PARDO mu, nu, lambda
#
allocate L0xxxi(mu,nu,lambda,*)
allocate L0xxxj(mu,nu,lambda,*)
#
DO sigma
#
compute_integrals AOINT(mu,nu,lambda,sigma)
#
DO i
#
V0xxxi(mu,nu,lambda,i) = AOINT(mu,nu,lambda,sigma)*ca(sigma,i)
L0xxxi(mu,nu,lambda,i) += V0xxxi(mu,nu,lambda,i)
#
ENDDO i
#
DO j
#
V0xxxj(mu,nu,lambda,j) = AOINT(mu,nu,lambda,sigma)*cb(sigma,j)
L0xxxj(mu,nu,lambda,j) += V0xxxj(mu,nu,lambda,j)
#
ENDDO j
#
ENDDO sigma
#
DO i
#
txixx(lambda,i,mu,nu) = L0xxxi(mu,nu,lambda,i)
txxix(nu,mu,i,lambda) = L0xxxi(mu,nu,lambda,i)
#
DO i1
#
Txixi(lambda,i,mu,i1) = txixx(lambda,i,mu,nu)*ca(nu,i1)
prepare Vxixi(lambda,i,mu,i1) += Txixi(lambda,i,mu,i1)
#
ENDDO i1
#
DO j1
#
Txixj(lambda,i,mu,j1) = txixx(lambda,i,mu,nu)*cb(nu,j1)
prepare Vxixj(lambda,i,mu,j1) += Txixj(lambda,i,mu,j1)
#
ENDDO j1
#
DO i1
#
Txxii(nu,mu,i,i1) = txxix(nu,mu,i,lambda)*ca(lambda,i1)
prepare Vxxii(nu,mu,i,i1) += Txxii(nu,mu,i,i1)
#
ENDDO i1
#
DO i1
#
Tixxi(i1,nu,lambda,i) = L0xxxi(mu,nu,lambda,i)*ca(mu,i1)
prepare Vixxi(i1,nu,lambda,i) += Tixxi(i1,nu,lambda,i)
#
ENDDO i1
#
ENDDO i
#
DO j
#
txjxx(lambda,j,mu,nu) = L0xxxj(mu,nu,lambda,j)
txxjx(nu,mu,j,lambda) = L0xxxj(mu,nu,lambda,j)
#
DO j1
#
Txjxj(lambda,j,mu,j1) = txjxx(lambda,j,mu,nu)*cb(nu,j1)
prepare Vxjxj(lambda,j,mu,j1) += Txjxj(lambda,j,mu,j1)
#
ENDDO j1
#
DO j1
#
Txxjj(nu,mu,j,j1) = txxjx(nu,mu,j,lambda)*cb(lambda,j1)
prepare Vxxjj(nu,mu,j,j1) += Txxjj(nu,mu,j,j1)
#
ENDDO j1
#
DO j1
#
Tjxxj(j1,nu,lambda,j) = L0xxxj(mu,nu,lambda,j)*cb(mu,j1)
prepare Vjxxj(j1,nu,lambda,j) += Tjxxj(j1,nu,lambda,j)
#
ENDDO j1
#
ENDDO j
#
DO j
#
DO i1
#
Tixxj(i1,nu,lambda,j) = L0xxxj(mu,nu,lambda,j)*ca(mu,i1)
prepare Vixxj(i1,nu,lambda,j) += Tixxj(i1,nu,lambda,j)
#
ENDDO i1
#
ENDDO j
#
deallocate L0xxxi(mu,nu,lambda,*)
deallocate L0xxxj(mu,nu,lambda,*)
#
ENDPARDO mu, nu, lambda
#
execute sip_barrier
execute server_barrier
#
# ------------------------------------------------------------------------
#
ENDPROC TRAN_XXII
#
# ------------------------------------------------------------------------
#
# ------------------------------------------------------------------------
#
PROC TRAN_PIPI
#
# ------------------------------------------------------------------------
#
execute sip_barrier
#
PARDO mu, i, i1
#
allocate Lxipi(mu,i,*,i1)
#
DO nu
#
REQUEST Vxixi(mu,i,nu,i1) i
#
DO p
Txipi(mu,i,p,i1) = Vxixi(mu,i,nu,i1)*ca(nu,p)
Lxipi(mu,i,p,i1) += Txipi(mu,i,p,i1)
ENDDO p
#
ENDDO nu
#
DO p
#
PREPARE Vxipi(mu,i,p,i1) = Lxipi(mu,i,p,i1)
#
ENDDO p
#
deallocate Lxipi(mu,i,*,i1)
#
ENDPARDO mu, i, i1
#
execute sip_barrier
execute server_barrier
# discard Vxixi
#
PARDO p, i, i1
#
allocate Lpipi(*,i,p,i1)
#
DO mu
#
REQUEST Vxipi(mu,i,p,i1) i
REQUEST Vxipi(mu,i1,p,i) i
Txipi(mu,i,p,i1) = Vxipi(mu,i1,p,i)
Txipi(mu,i,p,i1)-= Vxipi(mu,i,p,i1)
Txipi(mu,i,p,i1)*= -1.0
#
DO p1
#
#Tpipi(p1,i,p,i1) = Vxipi(mu,i,p,i1)*ca(mu,p1)
#T1pipi(p1,i,p,i1) = Vxipi(mu,i1,p,i)*ca(mu,p1)
#Tpipi(p1,i,p,i1) -= T1pipi(p1,i,p,i1)
#
Tpipi(p1,i,p,i1) = Txipi(mu,i,p,i1)*ca(mu,p1)
Lpipi(p1,i,p,i1) += Tpipi(p1,i,p,i1)
#
ENDDO p1
#
ENDDO mu
#
DO p1
#
PREPARE VSpipi(p1,i,p,i1) = Lpipi(p1,i,p,i1)
#
ENDDO p1
#
deallocate Lpipi(*,i,p,i1)
#
ENDPARDO p, i, i1
#
execute sip_barrier
execute server_barrier
# DISCARD Vxipi
#
# ------------------------------------------------------------------------
#
ENDPROC TRAN_PIPI
#
# ------------------------------------------------------------------------
#
# ------------------------------------------------------------------------
#
PROC TRAN_AAII
#
# ------------------------------------------------------------------------
#
execute sip_barrier
#
PARDO mu, i, i1
#
allocate Lxaii(mu,*,i,i1)
#
DO nu
#
REQUEST Vxxii(mu,nu,i,i1) i
#
DO a
Txaii(mu,a,i,i1) = Vxxii(mu,nu,i,i1)*ca(nu,a)
Lxaii(mu,a,i,i1) += Txaii(mu,a,i,i1)
ENDDO a
#
ENDDO nu
#
DO a
#
PREPARE Vxaii(mu,a,i,i1) = Lxaii(mu,a,i,i1)
#
ENDDO a
#
deallocate Lxaii(mu,*,i,i1)
#
ENDPARDO mu, i, i1
#
execute sip_barrier
execute server_barrier
#
PARDO a, i, i1
#
allocate Laaii(*,a,i,i1)
#
DO mu
#
REQUEST Vxaii(mu,a,i,i1) i
#
DO a1
#
Taaii(a1,a,i,i1) = Vxaii(mu,a,i,i1)*ca(mu,a1)
Laaii(a1,a,i,i1) += Taaii(a1,a,i,i1)
#
ENDDO a1
#
ENDDO mu
#
DO a1
#
PREPARE Vaaii(a1,a,i,i1) = Laaii(a1,a,i,i1)
#
ENDDO a1
#
deallocate Laaii(*,a,i,i1)
#
ENDPARDO a, i, i1
#
execute sip_barrier
execute server_barrier
# DISCARD Vxaii
#
# ------------------------------------------------------------------------
#
ENDPROC TRAN_AAII
#
# ------------------------------------------------------------------------
#
# ------------------------------------------------------------------------
#
PROC TRAN_IAAI
#
# ------------------------------------------------------------------------
#
execute sip_barrier
#
PARDO mu, i, i1
#
allocate Lixai(i,mu,*,i1)
#
DO nu
#
REQUEST Vixxi(i,mu,nu,i1) i
#
DO a
#
Tixai(i,mu,a,i1) = Vixxi(i,mu,nu,i1)*ca(nu,a)
Lixai(i,mu,a,i1) += Tixai(i,mu,a,i1)
#
ENDDO a
#
ENDDO nu
#
DO a
#
PREPARE Vixai(i,mu,a,i1) = Lixai(i,mu,a,i1)
#
ENDDO a
#
deallocate Lixai(i,mu,*,i1)
#
ENDPARDO mu, i, i1
#
execute sip_barrier
execute server_barrier
# DISCARD Vixxi
#
PARDO a, i, i1
#
allocate Liaai(i,*,a,i1)
#
DO mu
#
REQUEST Vixai(i,mu,a,i1) i
#
DO a1
#
Tiaai(i,a1,a,i1) = Vixai(i,mu,a,i1)*ca(mu,a1)
Liaai(i,a1,a,i1) += Tiaai(i,a1,a,i1)
#
ENDDO a1
#
ENDDO mu
#
DO a1
#
PREPARE Viaai(i,a1,a,i1) = Liaai(i,a1,a,i1)
#
ENDDO a1
#
deallocate Liaai(i,*,a,i1)
#
ENDPARDO a, i, i1
#
execute sip_barrier
execute server_barrier
# DISCARD Vixai
#
# ------------------------------------------------------------------------
#
ENDPROC TRAN_IAAI
#
# ------------------------------------------------------------------------
#
# ------------------------------------------------------------------------
#
PROC TRAN_QJQJ
#
# ------------------------------------------------------------------------
#
execute sip_barrier
#
PARDO mu, j, j1
#
allocate Lxjqj(mu,j,*,j1)
#
DO nu
#
REQUEST Vxjxj(mu,j,nu,j1) j
#
DO q
#
Txjqj(mu,j,q,j1) = Vxjxj(mu,j,nu,j1)*cb(nu,q)
Lxjqj(mu,j,q,j1) += Txjqj(mu,j,q,j1)
#
ENDDO q
#
ENDDO nu
#
DO q
#
PREPARE Vxjqj(mu,j,q,j1) = Lxjqj(mu,j,q,j1)
#
ENDDO q
#
deallocate Lxjqj(mu,j,*,j1)
#
ENDPARDO mu, j, j1
#
execute sip_barrier
execute server_barrier
# DISCARD Vxjxj
#
PARDO q, j, j1
#
allocate Lqjqj(*,j,q,j1)
#
DO mu
#
REQUEST Vxjqj(mu,j,q,j1) j
REQUEST Vxjqj(mu,j1,q,j) j
Txjqj(mu,j,q,j1) = Vxjqj(mu,j1,q,j)
Txjqj(mu,j,q,j1)-= Vxjqj(mu,j,q,j1)
Txjqj(mu,j,q,j1)*= -1.0
#
DO q1
#
#Tqjqj(q1,j,q,j1) = Vxjqj(mu,j,q,j1)*cb(mu,q1)
#T1qjqj(q1,j,q,j1) = Vxjqj(mu,j1,q,j)*cb(mu,q1)
#Tqjqj(q1,j,q,j1) -= T1qjqj(q1,j,q,j1)
Tqjqj(q1,j,q,j1) = Txjqj(mu,j,q,j1)*cb(mu,q1)
Lqjqj(q1,j,q,j1) += Tqjqj(q1,j,q,j1)
#
ENDDO q1
#
ENDDO mu
#
DO q1
#
PREPARE VSqjqj(q1,j,q,j1) = Lqjqj(q1,j,q,j1)
#
ENDDO q1
#
deallocate Lqjqj(*,j,q,j1)
#
ENDPARDO q, j, j1
#
execute sip_barrier
execute server_barrier
# DISCARD Vxjqj
#
# ------------------------------------------------------------------------
#
ENDPROC TRAN_QJQJ
#
# ------------------------------------------------------------------------
#
# ------------------------------------------------------------------------
#
PROC TRAN_BBJJ
#
# ------------------------------------------------------------------------
#
execute sip_barrier
#
PARDO mu, j, j1
#
allocate Lxbjj(mu,*,j,j1)
#
DO nu
#
REQUEST Vxxjj(mu,nu,j,j1) j
#
DO b
#
Txbjj(mu,b,j,j1) = Vxxjj(mu,nu,j,j1)*cb(nu,b)
Lxbjj(mu,b,j,j1) += Txbjj(mu,b,j,j1)
#
ENDDO b
#
ENDDO nu
#
DO b
#
PREPARE Vxbjj(mu,b,j,j1) = Lxbjj(mu,b,j,j1)
#
ENDDO b
#
deallocate Lxbjj(mu,*,j,j1)
#
ENDPARDO mu, j, j1
#
execute sip_barrier
execute server_barrier
#
PARDO b, j, j1
#
allocate Lbbjj(*,b,j,j1)
#
DO mu
#
REQUEST Vxbjj(mu,b,j,j1) j
#
DO b1
#
Tbbjj(b1,b,j,j1) = Vxbjj(mu,b,j,j1)*cb(mu,b1)
Lbbjj(b1,b,j,j1) += Tbbjj(b1,b,j,j1)
#
ENDDO b1
#
ENDDO mu
#
DO b1
#
PREPARE Vbbjj(b1,b,j,j1) = Lbbjj(b1,b,j,j1)
#
ENDDO b1
#
deallocate Lbbjj(*,b,j,j1)
#
ENDPARDO b, j, j1
#
execute sip_barrier
execute server_barrier
# DISCARD Vxbjj
#
# ------------------------------------------------------------------------
#
ENDPROC TRAN_BBJJ
#
# ------------------------------------------------------------------------
#
# ------------------------------------------------------------------------
#
PROC TRAN_JBBJ
#
# ------------------------------------------------------------------------
#
execute sip_barrier
#
PARDO mu, j, j1
#
allocate Ljxbj(j,nu,*,j1)
#
DO nu
#
REQUEST Vjxxj(j,mu,nu,j1) j
#
DO b
#
Tjxbj(j,mu,b,j1) = Vjxxj(j,mu,nu,j1)*cb(nu,b)
Ljxbj(j,mu,b,j1) += Tjxbj(j,mu,b,j1)
#
ENDDO b
#
ENDDO nu
#
DO b
#
PREPARE Vjxbj(j,mu,b,j1) = Ljxbj(j,mu,b,j1)
#
ENDDO b
#
deallocate Ljxbj(j,nu,*,j1)
#
ENDPARDO mu, j, j1
#
execute sip_barrier
execute server_barrier
# DISCARD Vjxxj
#
PARDO b, j, j1
#
allocate Ljbbj(j,*,b,j1)
#
DO mu
#
REQUEST Vjxbj(j,mu,b,j1) j
#
DO b1
#
Tjbbj(j,b1,b,j1) = Vjxbj(j,mu,b,j1)*cb(mu,b1)
Ljbbj(j,b1,b,j1) += Tjbbj(j,b1,b,j1)
#
ENDDO b1
#
ENDDO mu
#
DO b1
#
PREPARE Vjbbj(j,b1,b,j1) = Ljbbj(j,b1,b,j1)
#
ENDDO b1
#
deallocate Ljbbj(j,*,b,j1)
#
ENDPARDO b, j, j1
#
execute sip_barrier
execute server_barrier
# DISCARD Vjxbj
#
# ------------------------------------------------------------------------
#
ENDPROC TRAN_JBBJ
#
# ------------------------------------------------------------------------
#
# ------------------------------------------------------------------------
#
PROC TRAN_BBII
#
# ------------------------------------------------------------------------
#
execute sip_barrier
#
PARDO mu, i, i1
#
allocate Lxbii(mu,*,i,i1)
#
DO nu
#
REQUEST Vxxii(mu,nu,i,i1) i
#
DO b
#
Txbii(mu,b,i,i1) = Vxxii(mu,nu,i,i1)*cb(nu,b)
Lxbii(mu,b,i,i1) += Txbii(mu,b,i,i1)
#
ENDDO b
#
ENDDO nu
#
DO b
#
PREPARE Vxbii(mu,b,i,i1) = Lxbii(mu,b,i,i1)
#
ENDDO b
#
deallocate Lxbii(mu,*,i,i1)
#
ENDPARDO mu, i, i1
#
execute sip_barrier
execute server_barrier
# DISCARD Vxxii
#
PARDO b, i, i1
#
allocate Lbbii(*,b,i,i1)
#
DO mu
#
REQUEST Vxbii(mu,b,i,i1) i
#
DO b1
#
Tbbii(b1,b,i,i1) = Vxbii(mu,b,i,i1)*cb(mu,b1)
Lbbii(b1,b,i,i1) += Tbbii(b1,b,i,i1)
#
ENDDO b1
#
ENDDO mu
#
DO b1
#
PREPARE Vbbii(b1,b,i,i1) = Lbbii(b1,b,i,i1)
#
ENDDO b1
#
deallocate Lbbii(*,b,i,i1)
#
ENDPARDO b, i, i1
#
execute sip_barrier
execute server_barrier
# DISCARD Vxbii
#
# ------------------------------------------------------------------------
#
ENDPROC TRAN_BBII
#
# ------------------------------------------------------------------------
#
# ------------------------------------------------------------------------
#
PROC TRAN_AAJJ
#
# ------------------------------------------------------------------------
#
execute sip_barrier
#
PARDO mu, j, j1
#
allocate Lxajj(mu,*,j,j1)
#
DO nu
#
REQUEST Vxxjj(mu,nu,j,j1) j
#
DO a
#
Txajj(mu,a,j,j1) = Vxxjj(mu,nu,j,j1)*ca(nu,a)
Lxajj(mu,a,j,j1) += Txajj(mu,a,j,j1)
#
ENDDO a
#
ENDDO nu
#
DO a
#
PREPARE Vxajj(mu,a,j,j1) = Lxajj(mu,a,j,j1)
#
ENDDO a
#
deallocate Lxajj(mu,*,j,j1)
#
ENDPARDO mu, j, j1
#
execute sip_barrier
execute server_barrier
# DISCARD Vxxjj
#
PARDO a, j, j1
#
allocate Laajj(*,a,j,j1)
#
DO mu
#
REQUEST Vxajj(mu,a,j,j1) j
#
DO a1
#
Taajj(a1,a,j,j1) = Vxajj(mu,a,j,j1)*ca(mu,a1)
Laajj(a1,a,j,j1) += Taajj(a1,a,j,j1)
#
ENDDO a1
#
ENDDO mu
#
DO a1
#
PREPARE Vaajj(a1,a,j,j1) = Laajj(a1,a,j,j1)
#
ENDDO a1
#
deallocate Laajj(*,a,j,j1)
#
ENDPARDO a, j, j1
#
execute sip_barrier
execute server_barrier
# DISCARD Vxajj
#
# ------------------------------------------------------------------------
#
ENDPROC TRAN_AAJJ
#
# ------------------------------------------------------------------------
#
# ------------------------------------------------------------------------
#
PROC TRAN_IABJ
#
# ------------------------------------------------------------------------
#
execute sip_barrier
#
PARDO mu, i, j
#
allocate Lixbj(i,mu,*,j)
#
DO nu
#
REQUEST Vixxj(i,mu,nu,j) i
#
DO b
#
Tixbj(i,mu,b,j) = Vixxj(i,mu,nu,j)*cb(nu,b)
Lixbj(i,mu,b,j) += Tixbj(i,mu,b,j)
#
ENDDO b
#
ENDDO nu
#
DO b
#
PREPARE Vixbj(i,mu,b,j) = Lixbj(i,mu,b,j)
#
ENDDO b
#
deallocate Lixbj(i,mu,*,j)
#
ENDPARDO mu, i, j
#
execute sip_barrier
execute server_barrier
# DISCARD Vixxj
#
PARDO b, i, j
#
allocate Liabj(i,*,b,j)
#
DO mu
#
REQUEST Vixbj(i,mu,b,j) i
#
DO a
#
Tiabj(i,a,b,j) = Vixbj(i,mu,b,j)*ca(mu,a)
Liabj(i,a,b,j) += Tiabj(i,a,b,j)
#
ENDDO a
#
ENDDO mu
#
DO a
#
PREPARE Viabj(i,a,b,j) = Liabj(i,a,b,j)
#
ENDDO a
#
deallocate Liabj(i,*,b,j)
#
ENDPARDO b, i, j
#
execute sip_barrier
execute server_barrier
# DISCARD Vixbj
#
# ------------------------------------------------------------------------
#
ENDPROC TRAN_IABJ
#
# ------------------------------------------------------------------------
#
# ------------------------------------------------------------------------
#
PROC TRAN_IIBB
#
# ------------------------------------------------------------------------
#
execute sip_barrier
#
PARDO mu, i, i1
#
allocate Liixb(i,i1,mu,*)
#
DO nu
#
REQUEST Viixx(i,i1,mu,nu) i
#
DO b
#
Tiixb(i,i1,mu,b) = Viixx(i,i1,mu,nu)*cb(nu,b)
Liixb(i,i1,mu,b) += Tiixb(i,i1,mu,b)
#
ENDDO b
#
ENDDO nu
#
DO b
#
PREPARE Viixb(i,i1,mu,b) = Liixb(i,i1,mu,b)
#
ENDDO b
#
deallocate Liixb(i,i1,mu,*)
#
ENDPARDO mu, i, i1
#
execute sip_barrier
execute server_barrier
# DISCARD Viixx
#
PARDO b, i, i1
#
allocate Liibb(i,i1,*,b)
#
DO mu
#
REQUEST Viixb(i,i1,mu,b) i
#
DO b1
#
Tiibb(i,i1,b1,b) = Viixb(i,i1,mu,b)*cb(mu,b1)
Liibb(i,i1,b1,b) += Tiibb(i,i1,b1,b)
#
ENDDO b1
#
ENDDO mu
#
DO b1
#
PREPARE Viibb(i,i1,b1,b) = Liibb(i,i1,b1,b)
#
ENDDO b1
#
deallocate Liibb(i,i1,*,b)
#
ENDPARDO b, i, i1
#
execute sip_barrier
execute server_barrier
# DISCARD Viixb
#
# ------------------------------------------------------------------------
#
ENDPROC TRAN_IIBB
#
# ------------------------------------------------------------------------
#
# ------------------------------------------------------------------------
#
PROC TRAN_PIQJ
#
# ------------------------------------------------------------------------
#
execute sip_barrier
#
PARDO mu, i, j
#
allocate Lxiqj(mu,i,*,j)
#
DO nu
#
REQUEST Vxixj(mu,i,nu,j) i
#
DO q
#
Txiqj(mu,i,q,j) = Vxixj(mu,i,nu,j)*cb(nu,q)
Lxiqj(mu,i,q,j) += Txiqj(mu,i,q,j)
#
ENDDO q
#
ENDDO nu
#
DO q
#
PREPARE Vxiqj(mu,i,q,j) = Lxiqj(mu,i,q,j)
#
ENDDO q
#
deallocate Lxiqj(mu,i,*,j)
#
ENDPARDO mu, i, j
#
execute sip_barrier
execute server_barrier
# DISCARD Vxixj
#
PARDO q, i, j
#
allocate Lpiqj(*,i,q,j)
#
DO mu
#
REQUEST Vxiqj(mu,i,q,j) i
#
DO p
#
Tpiqj(p,i,q,j) = Vxiqj(mu,i,q,j)*ca(mu,p)
Lpiqj(p,i,q,j) += Tpiqj(p,i,q,j)
#
ENDDO p
#
ENDDO mu
#
DO p
#
PREPARE Vpiqj(p,i,q,j) = Lpiqj(p,i,q,j)
#
ENDDO p
#
deallocate Lpiqj(*,i,q,j)
#
ENDPARDO q, i, j
#
execute sip_barrier
execute server_barrier
# DISCARD Vxiqj
#
# ------------------------------------------------------------------------
#
ENDPROC TRAN_PIQJ
#
# ------------------------------------------------------------------------
#
# ------------------------------------------------------------------------
#
PROC TRAN_PPPP
#
# ------------------------------------------------------------------------
#
CALL TRAN_XXII
CALL TRAN_PIPI
CALL TRAN_AAII
CALL TRAN_IAAI
#
CALL TRAN_QJQJ
CALL TRAN_BBJJ
CALL TRAN_JBBJ
#
CALL TRAN_PIQJ
CALL TRAN_AAJJ
CALL TRAN_IABJ
#
CALL TRAN_BBII
#
# ------------------------------------------------------------------------
#
ENDPROC TRAN_PPPP
#
# ------------------------------------------------------------------------
#
# ------------------------------------------------------------------------
#
PROC TRAN_UHF
#
CALL TRAN_PPPP
CALL ENERGY
#
ENDPROC TRAN_UHF
#
# ------------------------------------------------------------------------
#
# -----------------------------------------------------------------------------
#
PROC WHFDENS
# ------------
#
PARDO mu, nu, i
#
Ixi(nu,i) = 0.0
#
DO i1
#
I1xi(nu,i) = ca(nu,i1)*fock_a(i,i1)
Ixi(nu,i) -= I1xi(nu,i)
#
ENDDO i1
#
Tpq(mu,nu) = ca(mu,i)*Ixi(nu,i)
PUT Whfa(mu,nu) += Tpq(mu,nu)
#
ENDPARDO mu, nu, i
#
PARDO mu, nu, j
#
Ixj(nu,j) = 0.0
#
DO j1
#
I1xj(nu,j) = cb(nu,j1)*fock_b(j,j1)
Ixj(nu,j) -= I1xj(nu,j)
#
ENDDO j1
#
Tpq(mu,nu) = cb(mu,j)*Ixj(nu,j)
PUT Whfb(mu,nu) += Tpq(mu,nu)
#
ENDPARDO mu, nu, j
execute sip_barrier
#
ENDPROC WHFDENS
# --------------
#
# -----------------------------------------------------------------------------
#
PROC HFDENS
# -----------
#
PARDO mu, nu, i
#
Txi(nu,i) = ca(nu,i)
Tpq(mu,nu) = ca(mu,i)*Txi(nu,i)
PUT Dhfa(mu,nu) += Tpq(mu,nu)
#
ENDPARDO mu, nu, i
#
PARDO mu, nu, j
#
Txj(nu,j) = cb(nu,j)
Tpq(mu,nu) = cb(mu,j)*Txj(nu,j)
PUT Dhfb(mu,nu) += Tpq(mu,nu)
#
ENDPARDO mu, nu, j
#
ENDPROC HFDENS
# --------------
#
PROC D1TRANS
# ------------
#
# Contract with the derivative integrals
# --------------------------------------
#
PARDO mu, nu
#
GET P2A_ao(mu,nu)
GET P2B_ao(mu,nu)
GET DHFA(mu,nu)
GET DHFB(mu,nu)
#
Tpq(mu,nu) = DHFA(mu,nu)
Tpq(mu,nu) += DHFB(mu,nu)
Tpq(mu,nu) += P2A_ao(mu,nu)
Tpq(mu,nu) += P2B_ao(mu,nu)
#
EXECUTE HCONT1 Tpq(mu,nu)
#
ENDPARDO mu, nu
#
# -----------------------------------------------------------------------------
#
ENDPROC D1TRANS
# ---------------
#
# -----------------------------------------------------------------------------
#
PROC S1TRANS
# ------------
#
# Contract with the derivative integrals
# --------------------------------------
#
PARDO mu, nu
#
GET W2_ao(mu,nu)
GET WHFa(mu,nu)
GET WHFb(mu,nu)
#
Tpq(mu,nu) = W2_ao(mu,nu)
Tpq(mu,nu) += WHFa(mu,nu)
Tpq(mu,nu) += WHFb(mu,nu)
#
EXECUTE SCONT1 Tpq(mu,nu)
#
ENDPARDO mu, nu
#
# -----------------------------------------------------------------------------
#
ENDPROC S1TRANS
# ---------------
#
# ---------------------------------------------------------------------------
#
# The following procedure computes the contribution to Lai which
# depends on the VVVO integrals directly.
#
PROC LAIAO1
# -----------
#
# Zero-out half backtransform vvoo integral arrays.
# -------------------------------------------------
#
execute sip_barrier
#
PARDO mu, nu, i, i1
txixi(mu,i,nu,i1) = 0.0
PREPARE Vxixi(mu,i,nu,i1) = txixi(mu,i,nu,i1)
ENDPARDO mu, nu, i, i1
#
PARDO mu, nu, i, j
txixj(mu,i,nu,j) = 0.0
PREPARE Vxixj(mu,i,nu,j) = txixj(mu,i,nu,j)
ENDPARDO mu, nu, i, j
#
PARDO mu, nu, j, j1
txjxj(mu,j,nu,j1) = 0.0
PREPARE Vxjxj(mu,j,nu,j1) = txjxj(mu,j,nu,j1)
ENDPARDO mu, nu, j, j1
#
execute server_barrier
#
# First half backtransform vvoo integrals
# ---------------------------------------
#
# AAAA spin case
# --------------
#
PARDO a1, i, i1
#
allocate Lxiai(*,i,a1,i1)
#
DO a
#
REQUEST VSpipi(a,i,a1,i1) i
Tpipi(a,i,a1,i1) = VSpipi(a,i,a1,i1)
execute energy_denominator Tpipi(a,i,a1,i1)
#
DO mu
#
Txipi(mu,i,a1,i1) = Tpipi(a,i,a1,i1)*ca(mu,a)
Lxiai(mu,i,a1,i1) += Txipi(mu,i,a1,i1)
#
ENDDO mu
#
ENDDO a
#
DO mu
#
PREPARE Vxipi(mu,i,a1,i1) = Lxiai(mu,i,a1,i1)
#
ENDDO mu
#
deallocate Lxiai(*,i,a1,i1)
#
ENDPARDO a1, i, i1
#
execute sip_barrier
execute server_barrier
#
PARDO mu, i, i1
#
allocate L1xixi(mu,i,*,i1)
#
DO a1
#
REQUEST Vxipi(mu,i,a1,i1) i
#
DO nu
#
Txixi(mu,i,nu,i1) = Vxipi(mu,i,a1,i1)*La(a1,nu) # ca(nu,a1)
L1xixi(mu,i,nu,i1) += Txixi(mu,i,nu,i1)
#
ENDDO nu
#
ENDDO a1
#
DO nu
#
PREPARE Vxixi(mu,i,nu,i1) = L1xixi(mu,i,nu,i1)
#
ENDDO nu
#
deallocate L1xixi(mu,i,*,i1)
#
ENDPARDO mu, i, i1
#
# BBBB spin case
# --------------
#
PARDO b1, j, j1
#
allocate Lxjbj(*,j,b1,j1)
#
DO b
#
REQUEST VSqjqj(b,j,b1,j1) j
Tqjqj(b,j,b1,j1) = VSqjqj(b,j,b1,j1)
execute energy_denominator Tqjqj(b,j,b1,j1)
#
DO mu
#
Txjqj(mu,j,b1,j1) = Tqjqj(b,j,b1,j1)*cb(mu,b)
Lxjbj(mu,j,b1,j1) += Txjqj(mu,j,b1,j1)
#
ENDDO mu
#
ENDDO b
#
DO mu
#
PREPARE Vxjqj(mu,j,b1,j1) = Lxjbj(mu,j,b1,j1)
#
ENDDO mu
#
deallocate Lxjbj(*,j,b1,j1)
#
ENDPARDO b1, j, j1
#
execute sip_barrier
execute server_barrier
#
PARDO mu, j, j1
#
allocate L1xjxj(mu,j,*,j1)
#
DO b1
#
REQUEST Vxjqj(mu,j,b1,j1) j
#
DO nu
#
Txjxj(mu,j,nu,j1) = Vxjqj(mu,j,b1,j1)*Lb(b1,nu) # cb(nu,b1)
L1xjxj(mu,j,nu,j1) += Txjxj(mu,j,nu,j1)
#
ENDDO nu
#
ENDDO b1
#
DO nu
#
PREPARE Vxjxj(mu,j,nu,j1) = L1xjxj(mu,j,nu,j1)
#
ENDDO nu
#
deallocate L1xjxj(mu,j,*,j1)
#
ENDPARDO mu, j, j1
#
# AABB spin case
# --------------
#
PARDO b, i, j
#
allocate Lxibj(*,i,b,j)
#
DO a
#
REQUEST Vpiqj(a,i,b,j) i
Tpiqj(a,i,b,j) = Vpiqj(a,i,b,j)
execute energy_denominator Tpiqj(a,i,b,j)
#
DO mu
#
Txiqj(mu,i,b,j) = Tpiqj(a,i,b,j)*ca(mu,a)
Lxibj(mu,i,b,j) += Txiqj(mu,i,b,j)
#
ENDDO mu
#
ENDDO a
#
DO mu
#
PREPARE Vxiqj(mu,i,b,j) = Lxibj(mu,i,b,j)
#
ENDDO mu
#
deallocate Lxibj(*,i,b,j)
#
ENDPARDO b, i, j
#
execute sip_barrier
execute server_barrier
#
PARDO mu, i, j
#
allocate L1xixj(mu,i,*,j)
#
DO b
#
REQUEST Vxiqj(mu,i,b,j) i
#
DO nu
#
Txixj(mu,i,nu,j) = Vxiqj(mu,i,b,j)*Lb(b,nu) # cb(nu,b)
L1xixj(mu,i,nu,j) += Txixj(mu,i,nu,j)
#
ENDDO nu
#
ENDDO b
#
DO nu
#
PREPARE Vxixj(mu,i,nu,j) = L1xixj(mu,i,nu,j)
#
ENDDO nu
#
deallocate L1xixj(mu,i,*,j)
#
ENDPARDO mu, i, j
#
create Lxi
create Lxj
execute sip_barrier
execute server_barrier
#
# Perform intermediate contraction to get contribution to Lxi and Lxj
# -------------------------------------------------------------------
#
PARDO mu, lambda, sigma
#
allocate xa(lambda,*,sigma,mu)
allocate xb(lambda,*,sigma,mu)
allocate xTa(sigma,mu,lambda,*)
allocate xTb(sigma,mu,lambda,*)
#
# Form Xa
# -------
#
DO i
#
T1xixx(sigma,i,lambda,mu) = 0.0
#
DO i1
#
REQUEST Vxixi(lambda,i,sigma,i1) i
REQUEST Vxixi(sigma,i,lambda,i1) i
#
Txixi(lambda,i,sigma,i1) = Vxixi(lambda,i,sigma,i1)
T1xixi(lambda,i,sigma,i1) = Vxixi(sigma,i,lambda,i1)
Txixi(lambda,i,sigma,i1) -= T1xixi(lambda,i,sigma,i1)
#
Txixx(lambda,i,sigma,mu) = Txixi(lambda,i,sigma,i1)*La(i1,mu) # ca(mu,i1)
Txixx(lambda,i,sigma,mu) *= 0.5
xa(lambda,i,sigma,mu) += Txixx(lambda,i,sigma,mu)
#
ENDDO i1
#
DO j
#
REQUEST Vxixj(sigma,i,lambda,j) i
#Txixx(lambda,i,sigma,mu) = Vxixj(sigma,i,lambda,j)*cb(mu,j)
Txixx(sigma,i,lambda,mu) = Vxixj(sigma,i,lambda,j)*Lb(j,mu)
T1xixx(sigma,i,lambda,mu)-= Txixx(sigma,i,lambda,mu)
#
ENDDO j
#
Txixx(lambda,i,sigma,mu) = T1xixx(sigma,i,lambda,mu)
xa(lambda,i,sigma,mu) += Txixx(lambda,i,sigma,mu)
xTa(sigma,mu,lambda,i) = xa(lambda,i,sigma,mu)
#
ENDDO i
#
# Form Xb
# -------
#
DO j
#
T1xxxj(lambda,mu,sigma,j) = 0.0
#
DO j1
#
REQUEST Vxjxj(lambda,j,sigma,j1) j
REQUEST Vxjxj(sigma,j,lambda,j1) j
#
Txjxj(lambda,j,sigma,j1) = Vxjxj(lambda,j,sigma,j1)
T1xjxj(lambda,j,sigma,j1) = Vxjxj(sigma,j,lambda,j1)
Txjxj(lambda,j,sigma,j1) -= T1xjxj(lambda,j,sigma,j1)
#
Txjxx(lambda,j,sigma,mu) = Txjxj(lambda,j,sigma,j1)*Lb(j1,mu) # cb(mu,j1)
Txjxx(lambda,j,sigma,mu) *= 0.5
xb(lambda,j,sigma,mu) += Txjxx(lambda,j,sigma,mu)
#
ENDDO j1
#
DO i
#
REQUEST Vxixj(lambda,i,sigma,j) i
#Txjxx(lambda,j,sigma,mu) = Vxixj(lambda,i,sigma,j)*ca(mu,i)
Txxxj(lambda,mu,sigma,j) = Vxixj(lambda,i,sigma,j)*La(i,mu)
T1xxxj(lambda,mu,sigma,j)-= Txxxj(lambda,mu,sigma,j)
#
ENDDO i
#
Txjxx(lambda,j,sigma,mu) = T1xxxj(lambda,mu,sigma,j)
xb(lambda,j,sigma,mu) += Txjxx(lambda,j,sigma,mu)
xTb(sigma,mu,lambda,j) = xb(lambda,j,sigma,mu)
#
ENDDO j
#
DO nu
#
compute_integrals aoint(nu,sigma,mu,lambda)
#
# Finish Lxi
# ----------
#
DO i
#
#Ixi(nu,i) = xa(lambda,i,sigma,mu)*aoint(mu,lambda,nu,sigma)
Ixi(nu,i) = aoint(nu,sigma,mu,lambda)*xTa(sigma,mu,lambda,i)
PUT Lxi(nu,i) += Ixi(nu,i)
#
ENDDO i
#
# Finish Lxj
# ----------
#
DO j
#
#Ixj(nu,j) = xb(lambda,j,sigma,mu)*aoint(mu,lambda,nu,sigma)
Ixj(nu,j) = aoint(nu,sigma,mu,lambda)*xTb(sigma,mu,lambda,j)
PUT Lxj(nu,j) += Ixj(nu,j)
#
ENDDO j
#
# Start third-term
# ----------------
#
GET Paa_ao(sigma,nu)
GET Pbb_ao(sigma,nu)
GET Paa_ao(sigma,mu)
GET Pbb_ao(sigma,mu)
#
I1xx(sigma,nu) = Paa_ao(sigma,nu)
I1xx(sigma,nu) += Pbb_ao(sigma,nu)
Ixx(mu,lambda) = aoint(nu,sigma,mu,lambda)*I1xx(sigma,nu)
I1xx(nu,lambda) = aoint(nu,sigma,mu,lambda)*Paa_ao(sigma,mu)
#
DO i
#
Ixi(lambda,i) = Ixx(mu,lambda)*ca(mu,i)
I1xi(lambda,i) = I1xx(nu,lambda)*ca(nu,i)
Ixi(lambda,i) -= I1xi(lambda,i)
#
PUT Yxi(lambda,i) += Ixi(lambda,i)
#
ENDDO i
#
I1xx(sigma,nu) = Pbb_ao(sigma,nu)
I1xx(sigma,nu) += Paa_ao(sigma,nu)
Ixx(mu,lambda) = aoint(nu,sigma,mu,lambda)*I1xx(sigma,nu)
I1xx(nu,lambda) = aoint(nu,sigma,mu,lambda)*Pbb_ao(sigma,mu)
#
DO j
#
Ixj(lambda,j) = Ixx(mu,lambda)*cb(mu,j)
I1xj(lambda,j) = I1xx(nu,lambda)*cb(nu,j)
Ixj(lambda,j) -= I1xj(lambda,j)
#
PUT Yxj(lambda,j) += Ixj(lambda,j)
#
ENDDO j
#
ENDDO nu
#
deallocate xa(lambda,*,sigma,mu)
deallocate xb(lambda,*,sigma,mu)
deallocate xTa(sigma,mu,lambda,*)
deallocate xTb(sigma,mu,lambda,*)
#
ENDPARDO mu, lambda, sigma
execute sip_barrier
#
PARDO lambda, a, i
#
GET Yxi(lambda,i)
Tai(a,i) = ca(lambda,a)*Yxi(lambda,i)
Tai(a,i) *= -1.0
PUT Lai_aa(a,i) += tai(a,i)
#
ENDPARDO lambda, a, i
#
PARDO lambda, b, j
#
GET Yxj(lambda,j)
Tbj(b,j) = cb(lambda,b)*Yxj(lambda,j)
Tbj(b,j) *= -1.0
PUT Lai_bb(b,j) += tbj(b,j)
#
ENDPARDO lambda, b, j
#
# Perform final transformation to get contribution to Lai
# -------------------------------------------------------
#
PARDO a, i, nu
#
GET Lxi(nu,i)
tai(a,i) = Lxi(nu,i)*ca(nu,a)
PUT Lai_aa(a,i) += tai(a,i)
#
ENDPARDO a, i, nu
#
PARDO b, j, nu
#
GET Lxj(nu,j)
tbj(b,j) = Lxj(nu,j)*cb(nu,b)
PUT Lai_bb(b,j) += tbj(b,j)
#
ENDPARDO b, j, nu
#
# Done direct contribution of Lai and Lbj
# ---------------------------------------
#
ENDPROC LAIAO1
# --------------
#
# ---------------------------------------------------------------------------
#
# ---------------------------------------------------------------------------
#
# Procedure DPQRSSEP computes the seperable part of the two-particle
# 'density' matrix.
#
PROC DPQRSSEP
# -------------
#
# Get 1-particle pieces
# ---------------------
#
GET DHFa(mu,lambda)
GET DHFa(mu,sigma)
GET DHFa(mu,nu)
GET DHFa(nu,sigma)
GET DHFa(nu,lambda)
GET DHFa(sigma,lambda)
#
GET DHFb(mu,lambda)
GET DHFb(mu,sigma)
GET DHFb(mu,nu)
GET DHFb(nu,sigma)
GET DHFb(nu,lambda)
GET DHFb(sigma,lambda)
#
GET P2A_ao(mu,lambda)
GET P2A_ao(mu,sigma)
GET P2A_ao(mu,nu)
GET P2A_ao(nu,lambda)
GET P2A_ao(nu,sigma)
GET P2A_ao(sigma,lambda)
#
GET P2B_ao(mu,lambda)
GET P2B_ao(mu,nu)
GET P2B_ao(mu,sigma)
GET P2B_ao(nu,lambda)
GET P2B_ao(nu,sigma)
GET P2B_ao(sigma,lambda)
#
# HF only
# -------
Txx(nu,sigma) = DHFa(nu,sigma)
Txx(nu,sigma) += DHFb(nu,sigma)
Txx(nu,sigma) += P2A_ao(nu,sigma)
Txx(nu,sigma) += P2B_ao(nu,sigma)
T1xx(mu,lambda) = DHFa(mu,lambda)
T1xx(mu,lambda) += DHFb(mu,lambda)
Txxxx(mu,lambda,nu,sigma) = T1xx(mu,lambda)^Txx(nu,sigma)
#
Txx(nu,lambda) = DHFa(nu,lambda)
Txx(nu,lambda) += P2A_ao(nu,lambda)
T1xxxx(mu,lambda,nu,sigma) = DHFa(mu,sigma)^Txx(nu,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
Txx(nu,lambda) = DHFb(nu,lambda)
Txx(nu,lambda) += P2B_ao(nu,lambda)
T1xxxx(mu,lambda,nu,sigma) = DHFb(mu,sigma)^Txx(nu,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
Txx(sigma,lambda) = DHFa(sigma,lambda)
Txx(sigma,lambda) += P2A_ao(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma) = DHFa(mu,nu)^Txx(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
Txx(sigma,lambda) = DHFb(sigma,lambda)
Txx(sigma,lambda) += P2B_ao(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma) = DHFb(mu,nu)^Txx(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
# Correlation Only
# ----------------
#
Txx(nu,sigma) = DHFA(nu,sigma)
Txx(nu,sigma) += DHFB(nu,sigma)
T1xxxx(mu,lambda,nu,sigma) = P2A_ao(mu,lambda)^Txx(nu,sigma)
Txxxx(mu,lambda,nu,sigma) += T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = P2B_ao(mu,lambda)^Txx(nu,sigma)
Txxxx(mu,lambda,nu,sigma) += T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = P2A_ao(mu,sigma)^DHFA(nu,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = P2A_ao(mu,nu)^DHFA(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = P2B_ao(mu,sigma)^DHFB(nu,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = P2B_ao(mu,nu)^DHFB(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
Txxxx(mu,lambda,nu,sigma) *= 0.5
#
ENDPROC DPQRSSEP
# ----------------
#
# ---------------------------------------------------------------------------
#
# -----------------------------------------------------------------------------
#
PROC D2TRANS
# ------------
#
allocate LDhfa(*,*)
allocate LDhfb(*,*)
allocate LP2A_ao(*,*)
allocate LP2B_ao(*,*)
execute sip_barrier
DO mu
DO nu
GET DHFa(mu,nu)
GET P2A_ao(mu,nu)
GET DHFb(mu,nu)
GET P2B_ao(mu,nu)
LDhfa(mu,nu) = DHFa(mu,nu)
LDhfb(mu,nu) = DHFb(mu,nu)
LP2A_ao(mu,nu) = P2A_ao(mu,nu)
LP2B_ao(mu,nu) = P2B_ao(mu,nu)
ENDDO nu
ENDDO mu
execute sip_barrier
#
PARDO mu, nu
#
DO lambda
DO sigma
#
IF mu < lambda
IF nu < sigma
#
D2(mu,lambda,nu,sigma) = 0.0
D3(mu,lambda,sigma,nu) = 0.0
#
DO i1
TAxxxi(mu,lambda,nu,i1) = 0.0
TBxxxi(mu,lambda,sigma,i1) = 0.0
DO i
REQUEST Vxixi(mu,i,nu,i1) i
REQUEST Vxixi(lambda,i,nu,i1) i
REQUEST Vxixi(mu,i,sigma,i1) i
REQUEST Vxixi(lambda,i,sigma,i1) i
#
T1xxxi(mu,lambda,nu,i1) = Vxixi(mu,i,nu,i1)*La(i,lambda) # ca(lambda,i)
T2xxxi(mu,lambda,nu,i1) = Vxixi(lambda,i,nu,i1)*ca(mu,i)
#
TAxxxi(mu,lambda,nu,i1) += T1xxxi(mu,lambda,nu,i1)
TAxxxi(mu,lambda,nu,i1) += T2xxxi(mu,lambda,nu,i1) # -
#
T3xxxi(mu,lambda,sigma,i1) = Vxixi(mu,i,sigma,i1)*La(i,lambda) # ca(lambda,i)
T4xxxi(mu,lambda,sigma,i1) = Vxixi(lambda,i,sigma,i1)*ca(mu,i)
#
TBxxxi(mu,lambda,sigma,i1)+= T4xxxi(mu,lambda,sigma,i1)
TBxxxi(mu,lambda,sigma,i1)+= T3xxxi(mu,lambda,sigma,i1) # -
#
ENDDO i
#
#Txxxx(mu,lambda,nu,sigma) = TBxxxi(mu,lambda,sigma,i1)*ca(nu,i1)
#D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
Txxxx(mu,lambda,sigma,nu) = TBxxxi(mu,lambda,sigma,i1)*La(i1,nu)
D3(mu,lambda,sigma,nu) += Txxxx(mu,lambda,sigma,nu)
#
T1xxxx(mu,lambda,nu,sigma) = TAxxxi(mu,lambda,nu,i1)*La(i1,sigma) # ca(sigma,i1)
D2(mu,lambda,nu,sigma) += T1xxxx(mu,lambda,nu,sigma)
ENDDO i1
#
DO j1
TAxxxj(mu,lambda,nu,j1) = 0.0
TBxxxj(mu,lambda,sigma,j1) = 0.0
DO j
REQUEST Vxjxj(mu,j,nu,j1) j
REQUEST Vxjxj(lambda,j,nu,j1) j
REQUEST Vxjxj(mu,j,sigma,j1) j
REQUEST Vxjxj(lambda,j,sigma,j1) j
#
T1xxxj(mu,lambda,nu,j1) = Vxjxj(mu,j,nu,j1)*Lb(j,lambda) # cb(lambda,j)
T2xxxj(mu,lambda,nu,j1) = Vxjxj(lambda,j,nu,j1)*cb(mu,j)
#
TAxxxj(mu,lambda,nu,j1) += T1xxxj(mu,lambda,nu,j1)
TAxxxj(mu,lambda,nu,j1) += T2xxxj(mu,lambda,nu,j1) # -
#
T3xxxj(mu,lambda,sigma,j1) = Vxjxj(mu,j,sigma,j1)*Lb(j,lambda) # cb(lambda,j)
T4xxxj(mu,lambda,sigma,j1) = Vxjxj(lambda,j,sigma,j1)*cb(mu,j)
#
TBxxxj(mu,lambda,sigma,j1)+= T4xxxj(mu,lambda,sigma,j1)
TBxxxj(mu,lambda,sigma,j1)+= T3xxxj(mu,lambda,sigma,j1) # -
#
ENDDO j
#
#Txxxx(mu,lambda,nu,sigma) = TBxxxj(mu,lambda,sigma,j1)*cb(nu,j1)
#D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
Txxxx(mu,lambda,sigma,nu) = TBxxxj(mu,lambda,sigma,j1)*Lb(j1,nu)
D3(mu,lambda,sigma,nu) += Txxxx(mu,lambda,sigma,nu)
#
T1xxxx(mu,lambda,nu,sigma) = TAxxxj(mu,lambda,nu,j1)*Lb(j1,sigma) # cb(sigma,j1)
D2(mu,lambda,nu,sigma) += T1xxxx(mu,lambda,nu,sigma)
ENDDO j1
#
DO j
TAxxxj(mu,lambda,nu,j) = 0.0
TBxxxj(mu,lambda,sigma,j) = 0.0
DO i
REQUEST Vxixj(mu,i,nu,j) j
REQUEST Vxixj(lambda,i,nu,j) j
REQUEST Vxixj(mu,i,sigma,j) j
REQUEST Vxixj(lambda,i,sigma,j) j
#
T1xxxj(mu,lambda,nu,j) = Vxixj(mu,i,nu,j)*La(i,lambda) # ca(lambda,i)
T2xxxj(mu,lambda,nu,j) = Vxixj(lambda,i,nu,j)*ca(mu,i)
#
TAxxxj(mu,lambda,nu,j) += T1xxxj(mu,lambda,nu,j)
TAxxxj(mu,lambda,nu,j) += T2xxxj(mu,lambda,nu,j) # -
#
T3xxxj(mu,lambda,sigma,j) = Vxixj(mu,i,sigma,j)*La(i,lambda) # ca(lambda,i)
T4xxxj(mu,lambda,sigma,j) = Vxixj(lambda,i,sigma,j)*ca(mu,i)
#
TBxxxj(mu,lambda,sigma,j)+= T4xxxj(mu,lambda,sigma,j)
TBxxxj(mu,lambda,sigma,j)+= T3xxxj(mu,lambda,sigma,j) # -
#
ENDDO i
#
#Txxxx(mu,lambda,nu,sigma) = TBxxxj(mu,lambda,sigma,j)*cb(nu,j)
#D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
#D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
Txxxx(mu,lambda,sigma,nu) = TBxxxj(mu,lambda,sigma,j)*Lb(j,nu)
Txxxx(mu,lambda,sigma,nu) *= 2.0
D3(mu,lambda,sigma,nu) += Txxxx(mu,lambda,sigma,nu)
#
T1xxxx(mu,lambda,nu,sigma) = TAxxxj(mu,lambda,nu,j)*Lb(j,sigma) # cb(sigma,j)
T1xxxx(mu,lambda,nu,sigma)*= 2.0
D2(mu,lambda,nu,sigma) += T1xxxx(mu,lambda,nu,sigma)
#D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
ENDDO j
#
# Get the separable part
# ----------------------
#
# Get 1-particle pieces
# ---------------------
#
# HF only
# -------
Txx(nu,sigma) = LDHFa(nu,sigma)
Txx(nu,sigma) += LDHFb(nu,sigma)
Txx(nu,sigma) += LP2A_ao(nu,sigma)
Txx(nu,sigma) += LP2B_ao(nu,sigma)
T1xx(mu,lambda) = LDHFa(mu,lambda)
T1xx(mu,lambda) += LDHFb(mu,lambda)
Txxxx(mu,lambda,nu,sigma) = T1xx(mu,lambda)^Txx(nu,sigma)
#
T2xx(nu,lambda) = LDHFa(nu,lambda)
T2xx(nu,lambda) += LP2A_ao(nu,lambda)
T1xxxx(mu,lambda,nu,sigma) = LDHFa(mu,sigma)^T2xx(nu,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T3xx(nu,lambda) = LDHFb(nu,lambda)
T3xx(nu,lambda) += LP2B_ao(nu,lambda)
T1xxxx(mu,lambda,nu,sigma) = LDHFb(mu,sigma)^T3xx(nu,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T4xx(sigma,lambda) = LDHFa(sigma,lambda)
T4xx(sigma,lambda) += LP2A_ao(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma) = LDHFa(mu,nu)^T4xx(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T5xx(sigma,lambda) = LDHFb(sigma,lambda)
T5xx(sigma,lambda) += LP2B_ao(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma) = LDHFb(mu,nu)^T5xx(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
# Correlation Only
# ----------------
#
T6xx(nu,sigma) = LDHFA(nu,sigma)
T6xx(nu,sigma) += LDHFB(nu,sigma)
T1xxxx(mu,lambda,nu,sigma) = LP2A_ao(mu,lambda)^T6xx(nu,sigma)
Txxxx(mu,lambda,nu,sigma) += T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = LP2B_ao(mu,lambda)^T6xx(nu,sigma)
Txxxx(mu,lambda,nu,sigma) += T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = LP2A_ao(mu,sigma)^LDHFA(nu,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = LP2A_ao(mu,nu)^LDHFA(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = LP2B_ao(mu,sigma)^LDHFB(nu,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = LP2B_ao(mu,nu)^LDHFB(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
Txxxx(mu,lambda,nu,sigma) *= 0.5
TSxxxx(mu,lambda,nu,sigma) = Txxxx(mu,lambda,nu,sigma)
TSxxxx(mu,lambda,nu,sigma) *= 4.0
#
# Add the the nonseparable part
# -----------------------------
Txxxx(mu,lambda,nu,sigma) = D3(mu,lambda,sigma,nu)
TSxxxx(mu,lambda,nu,sigma) += D2(mu,lambda,nu,sigma)
TSxxxx(mu,lambda,nu,sigma) += txxxx(mu,lambda,nu,sigma)
#
# Set up integrals
# ----------------
execute der_int_setup dx1(mu,lambda,nu,sigma)
execute der_int_setup dy1(mu,lambda,nu,sigma)
execute der_int_setup dz1(mu,lambda,nu,sigma)
execute der_int_setup dx2(mu,lambda,nu,sigma)
execute der_int_setup dy2(mu,lambda,nu,sigma)
execute der_int_setup dz2(mu,lambda,nu,sigma)
execute der_int_setup dx3(mu,lambda,nu,sigma)
execute der_int_setup dy3(mu,lambda,nu,sigma)
execute der_int_setup dz3(mu,lambda,nu,sigma)
execute der_int_setup dx4(mu,lambda,nu,sigma)
execute der_int_setup dy4(mu,lambda,nu,sigma)
execute der_int_setup dz4(mu,lambda,nu,sigma)
#
# Compute integral block
# ----------------------
execute compute_derivative_integrals
#
# Contract density with integral derivatives
# ------------------------------------------
execute DCONT2 TSxxxx(mu,lambda,nu,sigma)
#
ENDIF # nu < sigma
ENDIF # mu < lambda
#
ENDDO sigma
ENDDO lambda
#
ENDPARDO mu, nu
#
PARDO mu, nu
#
DO lambda
DO sigma
#
IF mu == lambda
IF nu < sigma
#
D2(mu,lambda,nu,sigma) = 0.0
D3(mu,lambda,sigma,nu) = 0.0
#
DO i1
TAxxxi(mu,lambda,nu,i1) = 0.0
TBxxxi(mu,lambda,sigma,i1) = 0.0
DO i
REQUEST Vxixi(mu,i,nu,i1) i
REQUEST Vxixi(mu,i,sigma,i1) i
#
T1xxxi(mu,lambda,nu,i1) = Vxixi(mu,i,nu,i1)*La(i,lambda) # ca(lambda,i)
TAxxxi(mu,lambda,nu,i1) += T1xxxi(mu,lambda,nu,i1)
T3xxxi(mu,lambda,sigma,i1) = Vxixi(mu,i,sigma,i1)*La(i,lambda) # ca(lambda,i)
TBxxxi(mu,lambda,sigma,i1)+= T3xxxi(mu,lambda,sigma,i1) # -
#
ENDDO i
#
#Txxxx(mu,lambda,nu,sigma) = TBxxxi(mu,lambda,sigma,i1)*ca(nu,i1)
#D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
Txxxx(mu,lambda,sigma,nu) = TBxxxi(mu,lambda,sigma,i1)*La(i1,nu)
D3(mu,lambda,sigma,nu) += Txxxx(mu,lambda,sigma,nu)
#
T1xxxx(mu,lambda,nu,sigma) = TAxxxi(mu,lambda,nu,i1)*La(i1,sigma) # ca(sigma,i1)
D2(mu,lambda,nu,sigma) += T1xxxx(mu,lambda,nu,sigma)
ENDDO i1
#
DO j1
TAxxxj(mu,lambda,nu,j1) = 0.0
TBxxxj(mu,lambda,sigma,j1) = 0.0
DO j
REQUEST Vxjxj(mu,j,nu,j1) j
REQUEST Vxjxj(mu,j,sigma,j1) j
#
T1xxxj(mu,lambda,nu,j1) = Vxjxj(mu,j,nu,j1)*Lb(j,lambda) # cb(lambda,j)
TAxxxj(mu,lambda,nu,j1) += T1xxxj(mu,lambda,nu,j1)
T3xxxj(mu,lambda,sigma,j1) = Vxjxj(mu,j,sigma,j1)*Lb(j,lambda) # cb(lambda,j)
TBxxxj(mu,lambda,sigma,j1)+= T3xxxj(mu,lambda,sigma,j1) # -
#
ENDDO j
#
#Txxxx(mu,lambda,nu,sigma) = TBxxxj(mu,lambda,sigma,j1)*cb(nu,j1)
#D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
Txxxx(mu,lambda,sigma,nu) = TBxxxj(mu,lambda,sigma,j1)*Lb(j1,nu)
D3(mu,lambda,sigma,nu) += Txxxx(mu,lambda,sigma,nu)
#
T1xxxx(mu,lambda,nu,sigma) = TAxxxj(mu,lambda,nu,j1)*Lb(j1,sigma) # cb(sigma,j1)
D2(mu,lambda,nu,sigma) += T1xxxx(mu,lambda,nu,sigma)
ENDDO j1
#
DO j
TAxxxj(mu,lambda,nu,j) = 0.0
TBxxxj(mu,lambda,sigma,j) = 0.0
DO i
REQUEST Vxixj(mu,i,nu,j) j
REQUEST Vxixj(mu,i,sigma,j) j
#
T1xxxj(mu,lambda,nu,j) = Vxixj(mu,i,nu,j)*La(i,lambda) # ca(lambda,i)
TAxxxj(mu,lambda,nu,j) += T1xxxj(mu,lambda,nu,j)
T3xxxj(mu,lambda,sigma,j) = Vxixj(mu,i,sigma,j)*La(i,lambda) # ca(lambda,i)
TBxxxj(mu,lambda,sigma,j)+= T3xxxj(mu,lambda,sigma,j) # -
#
ENDDO i
#
#Txxxx(mu,lambda,nu,sigma) = TBxxxj(mu,lambda,sigma,j)*cb(nu,j)
#D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
#D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
Txxxx(mu,lambda,sigma,nu) = TBxxxj(mu,lambda,sigma,j)*Lb(j,nu)
Txxxx(mu,lambda,sigma,nu)*= 2.0
D3(mu,lambda,sigma,nu) += Txxxx(mu,lambda,sigma,nu)
#
T1xxxx(mu,lambda,nu,sigma) = TAxxxj(mu,lambda,nu,j)*Lb(j,sigma) # cb(sigma,j)
T1xxxx(mu,lambda,nu,sigma)*= 2.0
D2(mu,lambda,nu,sigma) += T1xxxx(mu,lambda,nu,sigma)
#D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
ENDDO j
#
# Get the separable part
# ----------------------
#
# Get 1-particle pieces
# ---------------------
#
# HF only
# -------
Txx(nu,sigma) = LDHFa(nu,sigma)
Txx(nu,sigma) += LDHFb(nu,sigma)
Txx(nu,sigma) += LP2A_ao(nu,sigma)
Txx(nu,sigma) += LP2B_ao(nu,sigma)
T1xx(mu,lambda) = LDHFa(mu,lambda)
T1xx(mu,lambda) += LDHFb(mu,lambda)
Txxxx(mu,lambda,nu,sigma) = T1xx(mu,lambda)^Txx(nu,sigma)
#
T2xx(nu,lambda) = LDHFa(nu,lambda)
T2xx(nu,lambda) += LP2A_ao(nu,lambda)
T1xxxx(mu,lambda,nu,sigma) = LDHFa(mu,sigma)^T2xx(nu,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T3xx(nu,lambda) = LDHFb(nu,lambda)
T3xx(nu,lambda) += LP2B_ao(nu,lambda)
T1xxxx(mu,lambda,nu,sigma) = LDHFb(mu,sigma)^T3xx(nu,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T4xx(sigma,lambda) = LDHFa(sigma,lambda)
T4xx(sigma,lambda) += LP2A_ao(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma) = LDHFa(mu,nu)^T4xx(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T5xx(sigma,lambda) = LDHFb(sigma,lambda)
T5xx(sigma,lambda) += LP2B_ao(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma) = LDHFb(mu,nu)^T5xx(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
# Correlation Only
# ----------------
#
T6xx(nu,sigma) = LDHFA(nu,sigma)
T6xx(nu,sigma) += LDHFB(nu,sigma)
T1xxxx(mu,lambda,nu,sigma) = LP2A_ao(mu,lambda)^T6xx(nu,sigma)
Txxxx(mu,lambda,nu,sigma) += T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = LP2B_ao(mu,lambda)^T6xx(nu,sigma)
Txxxx(mu,lambda,nu,sigma) += T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = LP2A_ao(mu,sigma)^LDHFA(nu,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = LP2A_ao(mu,nu)^LDHFA(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = LP2B_ao(mu,sigma)^LDHFB(nu,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = LP2B_ao(mu,nu)^LDHFB(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
Txxxx(mu,lambda,nu,sigma) *= 0.5
TSxxxx(mu,lambda,nu,sigma) = Txxxx(mu,lambda,nu,sigma)
TSxxxx(mu,lambda,nu,sigma) *= 2.0
#
# Add the the nonseparable part
# -----------------------------
txxxx(mu,lambda,nu,sigma) = D3(mu,lambda,sigma,nu)
TSxxxx(mu,lambda,nu,sigma) += D2(mu,lambda,nu,sigma)
TSxxxx(mu,lambda,nu,sigma) += txxxx(mu,lambda,nu,sigma)
#
# Set up integrals
# ----------------
execute der_int_setup dx1(mu,lambda,nu,sigma)
execute der_int_setup dy1(mu,lambda,nu,sigma)
execute der_int_setup dz1(mu,lambda,nu,sigma)
execute der_int_setup dx2(mu,lambda,nu,sigma)
execute der_int_setup dy2(mu,lambda,nu,sigma)
execute der_int_setup dz2(mu,lambda,nu,sigma)
execute der_int_setup dx3(mu,lambda,nu,sigma)
execute der_int_setup dy3(mu,lambda,nu,sigma)
execute der_int_setup dz3(mu,lambda,nu,sigma)
execute der_int_setup dx4(mu,lambda,nu,sigma)
execute der_int_setup dy4(mu,lambda,nu,sigma)
execute der_int_setup dz4(mu,lambda,nu,sigma)
#
# Compute integral block
# ----------------------
execute compute_derivative_integrals
#
# Contract density with integral derivatives
# ------------------------------------------
execute DCONT2 TSxxxx(mu,lambda,nu,sigma)
#
ENDIF # nu < sigma
ENDIF # mu == lambda
#
ENDDO sigma
ENDDO lambda
#
ENDPARDO mu, nu
#
PARDO mu, nu
#
DO lambda
DO sigma
#
IF mu < lambda
IF nu == sigma
#
D2(mu,lambda,nu,sigma) = 0.0
D3(mu,lambda,sigma,nu) = 0.0
#
DO i1
TAxxxi(mu,lambda,nu,i1) = 0.0
TBxxxi(mu,lambda,nu,i1) = 0.0
DO i
REQUEST Vxixi(mu,i,nu,i1) i
REQUEST Vxixi(lambda,i,nu,i1) i
#
T1xxxi(mu,lambda,nu,i1) = Vxixi(mu,i,nu,i1)*La(i,lambda) # ca(lambda,i)
TAxxxi(mu,lambda,nu,i1) += T1xxxi(mu,lambda,nu,i1)
T2xxxi(mu,lambda,nu,i1) = Vxixi(lambda,i,nu,i1)*ca(mu,i)
TBxxxi(mu,lambda,nu,i1) += T2xxxi(mu,lambda,nu,i1) # -
#
ENDDO i
Txxxx(mu,lambda,nu,sigma) = TBxxxi(mu,lambda,nu,i1)*La(i1,sigma) # ca(sigma,i1)
D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
Txxxx(mu,lambda,nu,sigma) = TAxxxi(mu,lambda,nu,i1)*La(i1,sigma) # ca(sigma,i1)
D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
ENDDO i1
#
DO j1
TAxxxj(mu,lambda,nu,j1) = 0.0
TBxxxj(mu,lambda,nu,j1) = 0.0
DO j
REQUEST Vxjxj(mu,j,nu,j1) j
REQUEST Vxjxj(lambda,j,nu,j1) j
#
T1xxxj(mu,lambda,nu,j1) = Vxjxj(mu,j,nu,j1)*Lb(j,lambda) # cb(lambda,j)
TAxxxj(mu,lambda,nu,j1) += T1xxxj(mu,lambda,nu,j1)
T2xxxj(mu,lambda,nu,j1) = Vxjxj(lambda,j,nu,j1)*cb(mu,j)
TBxxxj(mu,lambda,nu,j1) += T2xxxj(mu,lambda,nu,j1) # -
#
ENDDO j
Txxxx(mu,lambda,nu,sigma) = TBxxxj(mu,lambda,nu,j1)*Lb(j1,sigma) # cb(sigma,j1)
D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
Txxxx(mu,lambda,nu,sigma) = TAxxxj(mu,lambda,nu,j1)*Lb(j1,sigma) # cb(sigma,j1)
D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
ENDDO j1
#
DO j
TAxxxj(mu,lambda,nu,j) = 0.0
TBxxxj(mu,lambda,nu,j) = 0.0
DO i
REQUEST Vxixj(mu,i,nu,j) j
REQUEST Vxixj(lambda,i,nu,j) j
#
T1xxxj(mu,lambda,nu,j) = Vxixj(mu,i,nu,j)*La(i,lambda) # ca(lambda,i)
TAxxxj(mu,lambda,nu,j) += T1xxxj(mu,lambda,nu,j)
T2xxxj(mu,lambda,nu,j) = Vxixj(lambda,i,nu,j)*ca(mu,i)
TBxxxj(mu,lambda,nu,j) += T2xxxj(mu,lambda,nu,j) # -
#
ENDDO i
#
Txxxx(mu,lambda,nu,sigma) = TBxxxj(mu,lambda,nu,j)*Lb(j,sigma) # cb(sigma,j)
D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
#
Txxxx(mu,lambda,nu,sigma) = TAxxxj(mu,lambda,nu,j)*Lb(j,sigma) # cb(sigma,j)
D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
ENDDO j
#
# Get the separable part
# ----------------------
#
# Get 1-particle pieces
# ---------------------
#
# HF only
# -------
Txx(nu,sigma) = LDHFa(nu,sigma)
Txx(nu,sigma) += LDHFb(nu,sigma)
Txx(nu,sigma) += LP2A_ao(nu,sigma)
Txx(nu,sigma) += LP2B_ao(nu,sigma)
T1xx(mu,lambda) = LDHFa(mu,lambda)
T1xx(mu,lambda) += LDHFb(mu,lambda)
Txxxx(mu,lambda,nu,sigma) = T1xx(mu,lambda)^Txx(nu,sigma)
#
T2xx(nu,lambda) = LDHFa(nu,lambda)
T2xx(nu,lambda) += LP2A_ao(nu,lambda)
T1xxxx(mu,lambda,nu,sigma) = LDHFa(mu,sigma)^T2xx(nu,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T3xx(nu,lambda) = LDHFb(nu,lambda)
T3xx(nu,lambda) += LP2B_ao(nu,lambda)
T1xxxx(mu,lambda,nu,sigma) = LDHFb(mu,sigma)^T3xx(nu,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T4xx(sigma,lambda) = LDHFa(sigma,lambda)
T4xx(sigma,lambda) += LP2A_ao(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma) = LDHFa(mu,nu)^T4xx(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T5xx(sigma,lambda) = LDHFb(sigma,lambda)
T5xx(sigma,lambda) += LP2B_ao(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma) = LDHFb(mu,nu)^T5xx(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
# Correlation Only
# ----------------
#
T6xx(nu,sigma) = LDHFA(nu,sigma)
T6xx(nu,sigma) += LDHFB(nu,sigma)
T1xxxx(mu,lambda,nu,sigma) = LP2A_ao(mu,lambda)^T6xx(nu,sigma)
Txxxx(mu,lambda,nu,sigma) += T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = LP2B_ao(mu,lambda)^T6xx(nu,sigma)
Txxxx(mu,lambda,nu,sigma) += T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = LP2A_ao(mu,sigma)^LDHFA(nu,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = LP2A_ao(mu,nu)^LDHFA(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = LP2B_ao(mu,sigma)^LDHFB(nu,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = LP2B_ao(mu,nu)^LDHFB(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
Txxxx(mu,lambda,nu,sigma) *= 0.5
TSxxxx(mu,lambda,nu,sigma) = Txxxx(mu,lambda,nu,sigma)
TSxxxx(mu,lambda,nu,sigma) *= 2.0
#
# Add the the nonseparable part
# -----------------------------
TSxxxx(mu,lambda,nu,sigma) += D2(mu,lambda,nu,sigma)
#
# Set up integrals
# ----------------
execute der_int_setup dx1(mu,lambda,nu,sigma)
execute der_int_setup dy1(mu,lambda,nu,sigma)
execute der_int_setup dz1(mu,lambda,nu,sigma)
execute der_int_setup dx2(mu,lambda,nu,sigma)
execute der_int_setup dy2(mu,lambda,nu,sigma)
execute der_int_setup dz2(mu,lambda,nu,sigma)
execute der_int_setup dx3(mu,lambda,nu,sigma)
execute der_int_setup dy3(mu,lambda,nu,sigma)
execute der_int_setup dz3(mu,lambda,nu,sigma)
execute der_int_setup dx4(mu,lambda,nu,sigma)
execute der_int_setup dy4(mu,lambda,nu,sigma)
execute der_int_setup dz4(mu,lambda,nu,sigma)
#
# Compute integral block
# ----------------------
execute compute_derivative_integrals
#
# Contract density with integral derivatives
# ------------------------------------------
execute DCONT2 TSxxxx(mu,lambda,nu,sigma)
#
ENDIF # nu == sigma
ENDIF # mu < lambda
#
ENDDO sigma
ENDDO lambda
#
ENDPARDO mu, nu
#
PARDO mu, nu
#
DO lambda
DO sigma
#
IF mu == lambda
IF nu == sigma
#
D2(mu,lambda,nu,sigma) = 0.0
#
DO i1
TAxxxi(mu,lambda,nu,i1) = 0.0
DO i
REQUEST Vxixi(mu,i,nu,i1) i
T1xxxi(mu,lambda,nu,i1) = Vxixi(mu,i,nu,i1)*La(i,lambda) # ca(lambda,i)
TAxxxi(mu,lambda,nu,i1) += T1xxxi(mu,lambda,nu,i1)
ENDDO i
Txxxx(mu,lambda,nu,sigma) = TAxxxi(mu,lambda,nu,i1)*La(i1,sigma) # ca(sigma,i1)
D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
ENDDO i1
#
DO j1
TAxxxj(mu,lambda,nu,j1) = 0.0
DO j
REQUEST Vxjxj(mu,j,nu,j1) j
T1xxxj(mu,lambda,nu,j1) = Vxjxj(mu,j,nu,j1)*Lb(j,lambda) # cb(lambda,j)
TAxxxj(mu,lambda,nu,j1) += T1xxxj(mu,lambda,nu,j1)
ENDDO j
Txxxx(mu,lambda,nu,sigma) = TAxxxj(mu,lambda,nu,j1)*Lb(j1,sigma) # cb(sigma,j1)
D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
ENDDO j1
#
DO j
TAxxxj(mu,lambda,nu,j) = 0.0
DO i
REQUEST Vxixj(mu,i,nu,j) i
T1xxxj(mu,lambda,nu,j) = Vxixj(mu,i,nu,j)*La(i,lambda) # ca(lambda,i)
TAxxxj(mu,lambda,nu,j) += T1xxxj(mu,lambda,nu,j)
ENDDO i
Txxxx(mu,lambda,nu,sigma) = TAxxxj(mu,lambda,nu,j)*Lb(j,sigma) # cb(sigma,j)
D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
D2(mu,lambda,nu,sigma) += Txxxx(mu,lambda,nu,sigma)
ENDDO j
#
# Get the separable part
# ----------------------
#
# Get 1-particle pieces
# ---------------------
#
# HF only
# -------
Txx(nu,sigma) = LDHFa(nu,sigma)
Txx(nu,sigma) += LDHFb(nu,sigma)
Txx(nu,sigma) += LP2A_ao(nu,sigma)
Txx(nu,sigma) += LP2B_ao(nu,sigma)
T1xx(mu,lambda) = LDHFa(mu,lambda)
T1xx(mu,lambda) += LDHFb(mu,lambda)
Txxxx(mu,lambda,nu,sigma) = T1xx(mu,lambda)^Txx(nu,sigma)
#
T2xx(nu,lambda) = LDHFa(nu,lambda)
T2xx(nu,lambda) += LP2A_ao(nu,lambda)
T1xxxx(mu,lambda,nu,sigma) = LDHFa(mu,sigma)^T2xx(nu,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T3xx(nu,lambda) = LDHFb(nu,lambda)
T3xx(nu,lambda) += LP2B_ao(nu,lambda)
T1xxxx(mu,lambda,nu,sigma) = LDHFb(mu,sigma)^T3xx(nu,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T4xx(sigma,lambda) = LDHFa(sigma,lambda)
T4xx(sigma,lambda) += LP2A_ao(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma) = LDHFa(mu,nu)^T4xx(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T5xx(sigma,lambda) = LDHFb(sigma,lambda)
T5xx(sigma,lambda) += LP2B_ao(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma) = LDHFb(mu,nu)^T5xx(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
# Correlation Only
# ----------------
#
T6xx(nu,sigma) = LDHFA(nu,sigma)
T6xx(nu,sigma) += LDHFB(nu,sigma)
T1xxxx(mu,lambda,nu,sigma) = LP2A_ao(mu,lambda)^T6xx(nu,sigma)
Txxxx(mu,lambda,nu,sigma) += T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = LP2B_ao(mu,lambda)^T6xx(nu,sigma)
Txxxx(mu,lambda,nu,sigma) += T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = LP2A_ao(mu,sigma)^LDHFA(nu,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = LP2A_ao(mu,nu)^LDHFA(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = LP2B_ao(mu,sigma)^LDHFB(nu,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
T1xxxx(mu,lambda,nu,sigma) = LP2B_ao(mu,nu)^LDHFB(sigma,lambda)
T1xxxx(mu,lambda,nu,sigma)*= 0.5
Txxxx(mu,lambda,nu,sigma) -= T1xxxx(mu,lambda,nu,sigma)
#
Txxxx(mu,lambda,nu,sigma) *= 0.5
TSxxxx(mu,lambda,nu,sigma) = Txxxx(mu,lambda,nu,sigma)
#
# Add the the nonseparable part
# -----------------------------
TSxxxx(mu,lambda,nu,sigma) += D2(mu,lambda,nu,sigma)
#
# Set up integrals
# ----------------
execute der_int_setup dx1(mu,lambda,nu,sigma)
execute der_int_setup dy1(mu,lambda,nu,sigma)
execute der_int_setup dz1(mu,lambda,nu,sigma)
execute der_int_setup dx2(mu,lambda,nu,sigma)
execute der_int_setup dy2(mu,lambda,nu,sigma)
execute der_int_setup dz2(mu,lambda,nu,sigma)
execute der_int_setup dx3(mu,lambda,nu,sigma)
execute der_int_setup dy3(mu,lambda,nu,sigma)
execute der_int_setup dz3(mu,lambda,nu,sigma)
execute der_int_setup dx4(mu,lambda,nu,sigma)
execute der_int_setup dy4(mu,lambda,nu,sigma)
execute der_int_setup dz4(mu,lambda,nu,sigma)
#
# Compute integral block
# ----------------------
execute compute_derivative_integrals
#
# Contract density with integral derivatives
# ------------------------------------------
execute DCONT2 TSxxxx(mu,lambda,nu,sigma)
#
ENDIF # nu == sigma
ENDIF # mu == lambda
#
ENDDO sigma
ENDDO lambda
#
ENDPARDO mu, nu
#
ENDPROC D2TRANS
# ---------------
#
# -----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------------------
#
PROC UPDATE_PAI
# ---------------
#
PARDO a, i
#
GET Painew_aa(a,i)
GET Paiold_aa(a,i)
Tai(a,i) = Painew_aa(a,i)
execute energy_denominator Tai(a,i)
Tai(a,i) -= Paiold_aa(a,i)
#
if kiter == 1
PUT e1ai(a,i) = tai(a,i)
endif
#
if kiter == 2
PUT e2ai(a,i) = tai(a,i)
endif
#
if kiter == 3
PUT e3ai(a,i) = tai(a,i)
endif
#
if kiter == 4
PUT e4ai(a,i) = tai(a,i)
endif
#
if kiter >= 5
PUT e5ai(a,i) = tai(a,i)
endif
#
ENDPARDO a, i
#
PARDO b, j
#
GET Painew_bb(b,j)
GET Paiold_bb(b,j)
Tbj(b,j) = Painew_bb(b,j)
execute energy_denominator Tbj(b,j)
Tbj(b,j) -= Paiold_bb(b,j)
#
if kiter == 1
PUT e1bj(b,j) = tbj(b,j)
endif
#
if kiter == 2
PUT e2bj(b,j) = tbj(b,j)
endif
#
if kiter == 3
PUT e3bj(b,j) = tbj(b,j)
endif
#
if kiter == 4
PUT e4bj(b,j) = tbj(b,j)
endif
#
if kiter >= 5
PUT e5bj(b,j) = tbj(b,j)
endif
#
ENDPARDO b, j
#
ENDPROC UPDATE_PAI
# ------------------
#
# ----------------------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------------------
#
PROC MOVE_PAI
# -------------
#
PARDO a, i
#
GET Paiold_aa(a,i)
tai(a,i) = Paiold_aa(a,i)
#
if kiter == 1
PUT d1ai(a,i) = tai(a,i)
endif
#
if kiter == 2
PUT d2ai(a,i) = Tai(a,i)
endif
#
if kiter == 3
PUT d3ai(a,i) = Tai(a,i)
endif
#
if kiter == 4
PUT d4ai(a,i) = Tai(a,i)
endif
#
if kiter >= 5
PUT d4ai(a,i) = Tai(a,i)
endif
#
ENDPARDO a, i
#
PARDO b, j
#
GET Paiold_bb(b,j)
tbj(b,j) = Paiold_bb(b,j)
#
if kiter == 1
PUT d1bj(b,j) = tbj(b,j)
endif
#
if kiter == 2
PUT d2bj(b,j) = tbj(b,j)
endif
#
if kiter == 3
PUT d3bj(b,j) = tbj(b,j)
endif
#
if kiter == 4
PUT d4bj(b,j) = tbj(b,j)
endif
#
if kiter >= 5
PUT d4bj(b,j) = tbj(b,j)
endif
#
ENDPARDO b, j
#
ENDPROC MOVE_PAI
# ----------------
#
# ----------------------------------------------------------------------------------------
#
# ---------------------------------------------------------------------------------
#
PROC ZERO_DSCALAR
# -----------------
#
# The scalars used (overlap of error arrays and coeffients) in the DIIS
# expansion are zero'd out.
#
# Zero out scalars.
# -----------------
#
b11 = 0.0
b12 = 0.0
b13 = 0.0
b14 = 0.0
b15 = 0.0
b16 = 0.0
b17 = 0.0
b18 = 0.0
b19 = 0.0
b110 = 0.0
#
b22 = 0.0
b23 = 0.0
b24 = 0.0
b25 = 0.0
b26 = 0.0
b27 = 0.0
b28 = 0.0
b29 = 0.0
b210 = 0.0
#
b33 = 0.0
b34 = 0.0
b35 = 0.0
b36 = 0.0
b37 = 0.0
b38 = 0.0
b39 = 0.0
b310 = 0.0
#
b44 = 0.0
b45 = 0.0
b46 = 0.0
b47 = 0.0
b48 = 0.0
b49 = 0.0
b410 = 0.0
#
b55 = 0.0
b56 = 0.0
b57 = 0.0
b58 = 0.0
b59 = 0.0
b510 = 0.0
#
b66 = 0.0
b67 = 0.0
b68 = 0.0
b69 = 0.0
b610 = 0.0
#
b77 = 0.0
b78 = 0.0
b79 = 0.0
b710 = 0.0
#
b88 = 0.0
b89 = 0.0
b810 = 0.0
#
b99 = 0.0
b910 = 0.0
#
b1010 = 0.0
#
Tb11 = 0.0
Tb12 = 0.0
Tb13 = 0.0
Tb14 = 0.0
Tb15 = 0.0
Tb16 = 0.0
Tb17 = 0.0
Tb18 = 0.0
Tb19 = 0.0
Tb110 = 0.0
#
Tb22 = 0.0
Tb23 = 0.0
Tb24 = 0.0
Tb25 = 0.0
Tb26 = 0.0
Tb27 = 0.0
Tb28 = 0.0
Tb29 = 0.0
Tb210 = 0.0
#
Tb33 = 0.0
Tb34 = 0.0
Tb35 = 0.0
Tb36 = 0.0
Tb37 = 0.0
Tb38 = 0.0
Tb39 = 0.0
Tb310 = 0.0
#
Tb44 = 0.0
Tb45 = 0.0
Tb46 = 0.0
Tb47 = 0.0
Tb48 = 0.0
Tb49 = 0.0
Tb410 = 0.0
#
Tb55 = 0.0
Tb56 = 0.0
Tb57 = 0.0
Tb58 = 0.0
Tb59 = 0.0
Tb510 = 0.0
#
Tb66 = 0.0
Tb67 = 0.0
Tb68 = 0.0
Tb69 = 0.0
Tb610 = 0.0
#
Tb77 = 0.0
Tb78 = 0.0
Tb79 = 0.0
Tb710 = 0.0
#
Tb88 = 0.0
Tb89 = 0.0
Tb810 = 0.0
#
Tb99 = 0.0
Tb910 = 0.0
#
Tb1010 = 0.0
#
c1 = 0.0
c2 = 0.0
c3 = 0.0
c4 = 0.0
c5 = 0.0
c6 = 0.0
c7 = 0.0
c8 = 0.0
c9 = 0.0
c10 = 0.0
#
execute sip_barrier
#
ENDPROC ZERO_DSCALAR
# --------------------
#
# ---------------------------------------------------------------------------------
#
# ---------------------------------------------------------------------------------
#
PROC SETUP_DIIS
# ---------------
#
# Put the matrix elements of B into the sip 'B' matrix.
#
execute diis_setup Tb11
execute diis_setup Tb12
execute diis_setup Tb13
execute diis_setup Tb14
execute diis_setup Tb15
execute diis_setup Tb16
execute diis_setup Tb17
execute diis_setup Tb18
execute diis_setup Tb19
execute diis_setup Tb110
execute diis_setup Tb22
execute diis_setup Tb23
execute diis_setup Tb24
execute diis_setup Tb25
execute diis_setup Tb26
execute diis_setup Tb27
execute diis_setup Tb28
execute diis_setup Tb29
execute diis_setup Tb210
execute diis_setup Tb33
execute diis_setup Tb34
execute diis_setup Tb35
execute diis_setup Tb36
execute diis_setup Tb37
execute diis_setup Tb38
execute diis_setup Tb39
execute diis_setup Tb310
execute diis_setup Tb44
execute diis_setup Tb45
execute diis_setup Tb46
execute diis_setup Tb47
execute diis_setup Tb48
execute diis_setup Tb49
execute diis_setup Tb410
execute diis_setup Tb55
execute diis_setup Tb56
execute diis_setup Tb57
execute diis_setup Tb58
execute diis_setup Tb59
execute diis_setup Tb510
execute diis_setup Tb66
execute diis_setup Tb67
execute diis_setup Tb68
execute diis_setup Tb69
execute diis_setup Tb610
execute diis_setup Tb77
execute diis_setup Tb78
execute diis_setup Tb79
execute diis_setup Tb710
execute diis_setup Tb88
execute diis_setup Tb89
execute diis_setup Tb810
execute diis_setup Tb99
execute diis_setup Tb910
execute diis_setup Tb1010
#
execute sip_barrier
#
ENDPROC SETUP_DIIS
# ------------------
#
# ---------------------------------------------------------------------------------
#
# ---------------------------------------------------------------------------------
#
PROC DIIS1
# ----------
#
# Zero out scalars.
# -----------------
#
CALL ZERO_DSCALAR
#
execute sip_barrier
#
# Determine the 'B-matrix'.
# -------------------------
#
# Compute contributions due to Dai amplitudes.
# --------------------------------------------
#
PARDO a, i
#
# REQUEST/GET amplitude data from all 2 previous iterations.
# ----------------------------------------------------------
#
GET e1ai(a,i) # kiter-1 amplitudes
GET e2ai(a,i) # kiter-0 amplitudes
#
# Compute contributions to the 'B-matrix'.
# ----------------------------------------
#
# b1x x=1,2
# ---------
#
etemp = e1ai(a,i)*e1ai(a,i)
b11 += etemp
#
etemp = e1ai(a,i)*e2ai(a,i)
b12 += etemp
#
# b1x x=2,2
# ---------
#
etemp = e2ai(a,i)*e2ai(a,i)
b22 += etemp
#
ENDPARDO a, i
#
# Compute contributions due to Dbj amplitudes.
# --------------------------------------------
#
PARDO b, j
#
# REQUEST/GET amplitude data from all 2 previous iterations.
# ----------------------------------------------------------
#
GET e1bj(b,j) # kiter-1 amplitudes
GET e2bj(b,j) # kiter-0 amplitudes
#
# Compute contributions to the 'B-matrix'.
# ----------------------------------------
#
# b1x x=1,2
# ---------
#
etemp = e1bj(b,j)*e1bj(b,j)
b11 += etemp
#
etemp = e1bj(b,j)*e2bj(b,j)
b12 += etemp
#
# b1x x=2,2
# ---------
#
etemp = e2bj(b,j)*e2bj(b,j)
b22 += etemp
#
ENDPARDO b, j
#
execute sip_barrier
#
# Collectively sum B-matrix elements.
# -----------------------------------
#
collective Tb11 += b11
collective Tb12 += b12
collective Tb22 += b22
execute server_barrier
#
# Now the unique elements of the 'B-matrix' have been computed, the array is filled out
# in the setup_diis.
# --------------------------------------------------------------------------------------
#
# Put the elements of the 'B-matrix', which have been computed as scalars into the
# R-matrix.
# --------------------------------------------------------------------------------
#
CALL SETUP_DIIS
#
execute compute_diis # --> New instruction
#
c1 = Tb11
c2 = Tb22
execute print_scalar c1
execute print_scalar c2
execute server_barrier
execute sip_barrier
#
# Done computing the c-vector.
# ----------------------------
#
# Form Pai_old.
# -------------
#
PARDO a, i
#
# REQUEST/GET amplitude data from all 2 previous iterations.
# ----------------------------------------------------------
#
GET D0ai(a,i) # kiter-2 amplitudes
GET D1ai(a,i) # kiter-1 amplitudes
#
GET e1ai(a,i) # kiter-1 amplitudes
GET e2ai(a,i) # kiter-0 amplitudes
#
# Compute contributions to updated amplitudes --> tai_old.
# --------------------------------------------------------
#
t1ai(a,i) = d0ai(a,i)
t1ai(a,i) += e1ai(a,i)
t1ai(a,i) *= c1
tai(a,i) = t1ai(a,i)
#
t1ai(a,i) = d1ai(a,i)
t1ai(a,i) += e2ai(a,i)
t1ai(a,i) *= c2
tai(a,i) += t1ai(a,i)
#
PUT Paiold_aa(a,i) = tai(a,i)
#
ENDPARDO a, i
#
# Form Dbj_old.
# -------------
#
PARDO b, j
#
# REQUEST/GET amplitude data from all 2 previous iterations.
# ----------------------------------------------------------
#
GET D0bj(b,j) # kiter-2 amplitudes
GET D1bj(b,j) # kiter-1 amplitudes
#
GET e1bj(b,j) # kiter-1 amplitudes
GET e2bj(b,j) # kiter-0 amplitudes
#
# Compute contributions to updated amplitudes --> tbj_old.
# --------------------------------------------------------
#
t1bj(b,j) = d0bj(b,j)
t1bj(b,j) += e1bj(b,j)
t1bj(b,j) *= c1
tbj(b,j) = t1bj(b,j)
#
t1bj(b,j) = d1bj(b,j)
t1bj(b,j) += e2bj(b,j)
t1bj(b,j) *= c2
tbj(b,j) += t1bj(b,j)
#
PUT Paiold_bb(b,j) = tbj(b,j)
#
ENDPARDO b, j
#
execute sip_barrier
#
ENDPROC DIIS1
# -------------
#
# ---------------------------------------------------------------------------------
#
PROC DIIS2
# ----------
#
# Zero out scalars.
# -----------------
#
CALL ZERO_DSCALAR
#
execute server_barrier
execute sip_barrier
#
# Determine the 'B-matrix'.
# -------------------------
#
# Compute contributions due to dai amplitudes.
# --------------------------------------------
#
PARDO a, i
#
# REQUEST/GET amplitude data from all 4 previous iterations.
# ----------------------------------------------------------
#
GET e1ai(a,i) # kiter-2 amplitudes
GET e2ai(a,i) # kiter-1 amplitudes
GET e3ai(a,i) # kiter-0 amplitudes
#
# Compute contributions to the 'B-matrix'.
# ----------------------------------------
#
# b1x x=1,3
# ---------
#
etemp = e1ai(a,i)*e1ai(a,i)
b11 += etemp
#
etemp = e1ai(a,i)*e2ai(a,i)
b12 += etemp
#
etemp = e1ai(a,i)*e3ai(a,i)
b13 += etemp
#
# b1x x=2,3
# ---------
#
etemp = e2ai(a,i)*e2ai(a,i)
b22 += etemp
#
etemp = e2ai(a,i)*e3ai(a,i)
b23 += etemp
#
# b1x x=3,3
# ---------
#
etemp = e3ai(a,i)*e3ai(a,i)
b33 += etemp
#
ENDPARDO a, i
#
# Compute contributions due to dbj amplitudes.
# --------------------------------------------
#
PARDO b, j
#
# REQUEST/GET amplitude data from all 3 previous iterations.
# ----------------------------------------------------------
#
GET e1bj(b,j) # kiter-2 amplitudes
GET e2bj(b,j) # kiter-1 amplitudes
GET e3bj(b,j) # kiter-0 amplitudes
#
# Compute contributions to the 'B-matrix'.
# ----------------------------------------
#
# b1x x=1,3
# ---------
#
etemp = e1bj(b,j)*e1bj(b,j)
b11 += etemp
#
etemp = e1bj(b,j)*e2bj(b,j)
b12 += etemp
#
etemp = e1bj(b,j)*e3bj(b,j)
b13 += etemp
#
# b1x x=2,3
# ---------
#
etemp = e2bj(b,j)*e2bj(b,j)
b22 += etemp
#
etemp = e2bj(b,j)*e3bj(b,j)
b23 += etemp
#
# b1x x=3,3
# ---------
#
etemp = e3bj(b,j)*e3bj(b,j)
b33 += etemp
#
ENDPARDO b, j
#
execute sip_barrier
#
# Collectively sum B-matrix elements.
# -----------------------------------
#
collective Tb11 += b11
collective Tb12 += b12
collective Tb13 += b13
collective Tb22 += b22
collective Tb23 += b23
collective Tb33 += b33
execute server_barrier
#
# Now the unique elements of the 'B-matrix' have been computed and the array filled out.
# --------------------------------------------------------------------------------------
#
# Put the elements of the 'B-matrix', which have been computed as scalars into the
# R-matrix.
# --------------------------------------------------------------------------------
CALL SETUP_DIIS
#
execute compute_diis # --> New instruction
#
c1 = Tb11
c2 = Tb22
c3 = Tb33
execute print_scalar c1
execute print_scalar c2
execute print_scalar c3
execute server_barrier
execute sip_barrier
#
# Done computing the c-vector.
# ----------------------------
#
# Form the updated amplitudes using the c-vector.
# -----------------------#
# Form Dai_old.
# -------------
#
PARDO a, i
#
# REQUEST/GET amplitude data from all 3 previous iterations.
# ----------------------------------------------------------
#
GET D0ai(a,i) # kiter-3 amplitudes
GET D1ai(a,i) # kiter-2 amplitudes
GET D2ai(a,i) # kiter-1 amplitudes
#
GET e1ai(a,i) # kiter-2 amplitudes
GET e2ai(a,i) # kiter-1 amplitudes
GET e3ai(a,i) # kiter-0 amplitudes
#
# Compute contributions to updated amplitudes --> tai_old.
# --------------------------------------------------------
#
t1ai(a,i) = d0ai(a,i)
t1ai(a,i) += e1ai(a,i)
t1ai(a,i) *= c1
tai(a,i) = t1ai(a,i)
#
t1ai(a,i) = d1ai(a,i)
t1ai(a,i) += e2ai(a,i)
t1ai(a,i) *= c2
tai(a,i) += t1ai(a,i)
#
t1ai(a,i) = d2ai(a,i)
t1ai(a,i) += e3ai(a,i)
t1ai(a,i) *= c3
tai(a,i) += t1ai(a,i)
#
PUT Paiold_aa(a,i) = tai(a,i)
#
ENDPARDO a, i
#
# Form tbj_old.
# -------------
#
PARDO b, j
#
# REQUEST/GET amplitude data from all 3 previous iterations.
# ----------------------------------------------------------
#
GET D0bj(b,j) # kiter-3 amplitudes
GET D1bj(b,j) # kiter-2 amplitudes
GET D2bj(b,j) # kiter-1 amplitudes
#
GET e1bj(b,j) # kiter-2 amplitudes
GET e2bj(b,j) # kiter-1 amplitudes
GET e3bj(b,j) # kiter-0 amplitudes
#
# Compute contributions to updated amplitudes --> tbj_old.
# --------------------------------------------------------
#
t1bj(b,j) = d0bj(b,j)
t1bj(b,j) += e1bj(b,j)
t1bj(b,j) *= c1
tbj(b,j) = t1bj(b,j)
#
t1bj(b,j) = d1bj(b,j)
t1bj(b,j) += e2bj(b,j)
t1bj(b,j) *= c2
tbj(b,j) += t1bj(b,j)
#
t1bj(b,j) = d2bj(b,j)
t1bj(b,j) += e3bj(b,j)
t1bj(b,j) *= c3
tbj(b,j) += t1bj(b,j)
#
PUT Paiold_bb(b,j) = tbj(b,j)
#
ENDPARDO b, j
#
execute sip_barrier
#
ENDPROC DIIS2
# -------------
#
# ---------------------------------------------------------------------------------
#
# ---------------------------------------------------------------------------------
#
PROC DIIS3
# ----------
#
# Zero out scalars.
# -----------------
#
CALL ZERO_DSCALAR
#
execute sip_barrier
#
# Determine the 'B-matrix'.
# -------------------------
#
# Compute contributions due to dai amplitudes.
# --------------------------------------------
#
PARDO a, i
#
# REQUEST/GET amplitude data from all 4 previous iterations.
# ----------------------------------------------------------
#
GET e1ai(a,i) # kiter-3 amplitudes
GET e2ai(a,i) # kiter-2 amplitudes
GET e3ai(a,i) # kiter-1 amplitudes
GET e4ai(a,i) # kiter-0 amplitudes
#
# Compute contributions to the 'B-matrix'.
# ----------------------------------------
#
# b1x x=1,4
# ---------
#
etemp = e1ai(a,i)*e1ai(a,i)
b11 += etemp
#
etemp = e1ai(a,i)*e2ai(a,i)
b12 += etemp
#
etemp = e1ai(a,i)*e3ai(a,i)
b13 += etemp
#
etemp = e1ai(a,i)*e4ai(a,i)
b14 += etemp
#
# b1x x=2,4
# ---------
#
etemp = e2ai(a,i)*e2ai(a,i)
b22 += etemp
#
etemp = e2ai(a,i)*e3ai(a,i)
b23 += etemp
#
etemp = e2ai(a,i)*e4ai(a,i)
b24 += etemp
#
# b1x x=3,4
# ---------
#
etemp = e3ai(a,i)*e3ai(a,i)
b33 += etemp
#
etemp = e3ai(a,i)*e4ai(a,i)
b34 += etemp
#
# b1x x=4,4
# ---------
#
etemp = e4ai(a,i)*e4ai(a,i)
b44 += etemp
#
ENDPARDO a, i
#
# Compute contributions due to dbj amplitudes.
# --------------------------------------------
#
PARDO b, j
#
# REQUEST/GET amplitude data from all 3 previous iterations.
# ----------------------------------------------------------
#
GET e1bj(b,j) # kiter-3 amplitudes
GET e2bj(b,j) # kiter-2 amplitudes
GET e3bj(b,j) # kiter-1 amplitudes
GET e4bj(b,j) # kiter-0 amplitudes
#
# Compute contributions to the 'B-matrix'.
# ----------------------------------------
#
# b1x x=1,4
# ---------
#
etemp = e1bj(b,j)*e1bj(b,j)
b11 += etemp
#
etemp = e1bj(b,j)*e2bj(b,j)
b12 += etemp
#
etemp = e1bj(b,j)*e3bj(b,j)
b13 += etemp
#
etemp = e1bj(b,j)*e4bj(b,j)
b14 += etemp
#
# b1x x=2,4
# ---------
#
etemp = e2bj(b,j)*e2bj(b,j)
b22 += etemp
#
etemp = e2bj(b,j)*e3bj(b,j)
b23 += etemp
#
etemp = e2bj(b,j)*e4bj(b,j)
b24 += etemp
#
# b1x x=3,4
# ---------
#
etemp = e3bj(b,j)*e3bj(b,j)
b33 += etemp
#
etemp = e3bj(b,j)*e4bj(b,j)
b34 += etemp
#
# b1x x=4,4
# ---------
#
etemp = e4bj(b,j)*e4bj(b,j)
b44 += etemp
#
ENDPARDO b, j
#
execute sip_barrier
#
# Collectively sum B-matrix elements.
# -----------------------------------
#
collective Tb11 += b11
collective Tb12 += b12
collective Tb13 += b13
collective Tb14 += b14
collective Tb22 += b22
collective Tb23 += b23
collective Tb24 += b24
collective Tb33 += b33
collective Tb34 += b34
collective Tb44 += b44
execute server_barrier
#
# Now the unique elements of the 'B-matrix' have been computed and the array filled out.
# --------------------------------------------------------------------------------------
#
# Put the elements of the 'B-matrix', which have been computed as scalars into the
# R-matrix.
# --------------------------------------------------------------------------------
CALL SETUP_DIIS
#
execute compute_diis # --> New instruction
#
c1 = Tb11
c2 = Tb22
c3 = Tb33
c4 = Tb44
execute print_scalar c1
execute print_scalar c2
execute print_scalar c3
execute print_scalar c4
execute server_barrier
execute sip_barrier
#
# Done computing the c-vector.
# ----------------------------
#
# Form the updated amplitudes using the c-vector.
# -----------------------#
#
# Form Dai_old.
# -------------
#
PARDO a, i
#
# REQUEST/GET amplitude data from all 3 previous iterations.
# ----------------------------------------------------------
#
GET D0ai(a,i) # kiter-4 amplitudes
GET D1ai(a,i) # kiter-3 amplitudes
GET D2ai(a,i) # kiter-2 amplitudes
GET D3ai(a,i) # kiter-1 amplitudes
#
GET e1ai(a,i) # kiter-3 amplitudes
GET e2ai(a,i) # kiter-2 amplitudes
GET e3ai(a,i) # kiter-1 amplitudes
GET e4ai(a,i) # kiter-0 amplitudes
#
# Compute contributions to updated amplitudes --> tai_old.
# --------------------------------------------------------
#
t1ai(a,i) = d0ai(a,i)
t1ai(a,i) += e1ai(a,i)
t1ai(a,i) *= c1
tai(a,i) = t1ai(a,i)
#
t1ai(a,i) = d1ai(a,i)
t1ai(a,i) += e2ai(a,i)
t1ai(a,i) *= c2
tai(a,i) += t1ai(a,i)
#
t1ai(a,i) = d2ai(a,i)
t1ai(a,i) += e3ai(a,i)
t1ai(a,i) *= c3
tai(a,i) += t1ai(a,i)
#
t1ai(a,i) = d3ai(a,i)
t1ai(a,i) += e4ai(a,i)
t1ai(a,i) *= c4
tai(a,i) += t1ai(a,i)
#
PUT Paiold_aa(a,i) = tai(a,i)
#
ENDPARDO a, i
#
# Form tbj_old.
# -------------
#
PARDO b, j
#
# REQUEST/GET amplitude data from all 3 previous iterations.
# ----------------------------------------------------------
#
GET D0bj(b,j) # kiter-4 amplitudes
GET D1bj(b,j) # kiter-3 amplitudes
GET D2bj(b,j) # kiter-2 amplitudes
GET D3bj(b,j) # kiter-1 amplitudes
#
GET e1bj(b,j) # kiter-3 amplitudes
GET e2bj(b,j) # kiter-2 amplitudes
GET e3bj(b,j) # kiter-1 amplitudes
GET e4bj(b,j) # kiter-0 amplitudes
#
# Compute contributions to updated amplitudes --> tbj_old.
# --------------------------------------------------------
#
t1bj(b,j) = d0bj(b,j)
t1bj(b,j) += e1bj(b,j)
t1bj(b,j) *= c1
tbj(b,j) = t1bj(b,j)
#
t1bj(b,j) = d1bj(b,j)
t1bj(b,j) += e2bj(b,j)
t1bj(b,j) *= c2
tbj(b,j) += t1bj(b,j)
#
t1bj(b,j) = d2bj(b,j)
t1bj(b,j) += e3bj(b,j)
t1bj(b,j) *= c3
tbj(b,j) += t1bj(b,j)
#
t1bj(b,j) = d3bj(b,j)
t1bj(b,j) += e4bj(b,j)
t1bj(b,j) *= c4
tbj(b,j) += t1bj(b,j)
#
PUT Paiold_bb(b,j) = tbj(b,j)
#
ENDPARDO b, j
#
execute sip_barrier
#
ENDPROC DIIS3
# -------------
#
# ---------------------------------------------------------------------------------
#
# ---------------------------------------------------------------------------------
#
PROC DIIS4
# ----------
#
# Zero out scalars.
# -----------------
#
CALL ZERO_DSCALAR
#
execute sip_barrier
#
# Determine the 'B-matrix'.
# -------------------------
#
# Compute contributions due to dai amplitudes.
# --------------------------------------------
#
PARDO a, i
#
# REQUEST/GET amplitude data from all 4 previous iterations.
# ----------------------------------------------------------
#
GET e1ai(a,i) # kiter-4 amplitudes
GET e2ai(a,i) # kiter-3 amplitudes
GET e3ai(a,i) # kiter-2 amplitudes
GET e4ai(a,i) # kiter-1 amplitudes
GET e5ai(a,i) # kiter-0 amplitudes
#
# Compute contributions to the 'B-matrix'.
# ----------------------------------------
#
# b1x x=1,5
# ---------
#
etemp = e1ai(a,i)*e1ai(a,i)
b11 += etemp
#
etemp = e1ai(a,i)*e2ai(a,i)
b12 += etemp
#
etemp = e1ai(a,i)*e3ai(a,i)
b13 += etemp
#
etemp = e1ai(a,i)*e4ai(a,i)
b14 += etemp
#
etemp = e1ai(a,i)*e5ai(a,i)
b15 += etemp
#
# b1x x=2,5
# ---------
#
etemp = e2ai(a,i)*e2ai(a,i)
b22 += etemp
#
etemp = e2ai(a,i)*e3ai(a,i)
b23 += etemp
#
etemp = e2ai(a,i)*e4ai(a,i)
b24 += etemp
#
etemp = e2ai(a,i)*e5ai(a,i)
b25 += etemp
#
# b1x x=3,5
# ---------
#
etemp = e3ai(a,i)*e3ai(a,i)
b33 += etemp
#
etemp = e3ai(a,i)*e4ai(a,i)
b34 += etemp
#
etemp = e3ai(a,i)*e5ai(a,i)
b35 += etemp
#
# b1x x=4,5
# ---------
#
etemp = e4ai(a,i)*e4ai(a,i)
b44 += etemp
#
etemp = e4ai(a,i)*e5ai(a,i)
b45 += etemp
#
# b1x x=5,5
# ---------
#
etemp = e5ai(a,i)*e5ai(a,i)
b55 += etemp
#
ENDPARDO a, i
#
# Compute contributions due to dbj amplitudes.
# --------------------------------------------
#
PARDO b, j
#
# REQUEST/GET amplitude data from all 3 previous iterations.
# ----------------------------------------------------------
#
GET e1bj(b,j) # kiter-4 amplitudes
GET e2bj(b,j) # kiter-3 amplitudes
GET e3bj(b,j) # kiter-2 amplitudes
GET e4bj(b,j) # kiter-1 amplitudes
GET e5bj(b,j) # kiter-0 amplitudes
#
# Compute contributions to the 'B-matrix'.
# ----------------------------------------
#
# b1x x=1,5
# ---------
#
etemp = e1bj(b,j)*e1bj(b,j)
b11 += etemp
#
etemp = e1bj(b,j)*e2bj(b,j)
b12 += etemp
#
etemp = e1bj(b,j)*e3bj(b,j)
b13 += etemp
#
etemp = e1bj(b,j)*e4bj(b,j)
b14 += etemp
#
etemp = e1bj(b,j)*e5bj(b,j)
b15 += etemp
#
# b1x x=2,5
# ---------
#
etemp = e2bj(b,j)*e2bj(b,j)
b22 += etemp
#
etemp = e2bj(b,j)*e3bj(b,j)
b23 += etemp
#
etemp = e2bj(b,j)*e4bj(b,j)
b24 += etemp
#
etemp = e2bj(b,j)*e5bj(b,j)
b25 += etemp
#
# b1x x=3,5
# ---------
#
etemp = e3bj(b,j)*e3bj(b,j)
b33 += etemp
#
etemp = e3bj(b,j)*e4bj(b,j)
b34 += etemp
#
etemp = e3bj(b,j)*e5bj(b,j)
b35 += etemp
#
# b1x x=4,5
# ---------
#
etemp = e4bj(b,j)*e4bj(b,j)
b44 += etemp
#
etemp = e4bj(b,j)*e5bj(b,j)
b45 += etemp
#
# b1x x=5,5
# ---------
#
etemp = e5bj(b,j)*e5bj(b,j)
b55 += etemp
#
ENDPARDO b, j
#
execute sip_barrier
#
# Collectively sum B-matrix elements.
# -----------------------------------
#
collective Tb11 += b11
collective Tb12 += b12
collective Tb13 += b13
collective Tb14 += b14
collective Tb15 += b15
collective Tb22 += b22
collective Tb23 += b23
collective Tb24 += b24
collective Tb25 += b25
collective Tb33 += b33
collective Tb34 += b34
collective Tb35 += b35
collective Tb44 += b44
collective Tb45 += b45
collective Tb55 += b55
execute server_barrier
#
# Now the unique elements of the 'B-matrix' have been computed and the array filled out.
# --------------------------------------------------------------------------------------
#
# Put the elements of the 'B-matrix', which have been computed as scalars into the
# R-matrix.
# --------------------------------------------------------------------------------
CALL SETUP_DIIS
#
execute compute_diis # --> New instruction
#
c1 = Tb11
c2 = Tb22
c3 = Tb33
c4 = Tb44
c5 = Tb55
execute print_scalar c1
execute print_scalar c2
execute print_scalar c3
execute print_scalar c4
execute print_scalar c5
execute sip_barrier
#
# Done computing the c-vector.
# ----------------------------
#
# Form the updated amplitudes using the c-vector.
# -----------------------#
# Form Dai_old.
# -------------
#
PARDO a, i
#
# REQUEST/GET amplitude data from all 3 previous iterations.
# ----------------------------------------------------------
#
GET D0ai(a,i) # kiter-5 amplitudes
GET D1ai(a,i) # kiter-4 amplitudes
GET D2ai(a,i) # kiter-3 amplitudes
GET D3ai(a,i) # kiter-2 amplitudes
GET D4ai(a,i) # kiter-1 amplitudes
#
GET e1ai(a,i) # kiter-4 amplitudes
GET e2ai(a,i) # kiter-3 amplitudes
GET e3ai(a,i) # kiter-2 amplitudes
GET e4ai(a,i) # kiter-1 amplitudes
GET e5ai(a,i) # kiter-0 amplitudes
#
# Compute contributions to updated amplitudes --> tai_old.
# --------------------------------------------------------
#
t1ai(a,i) = d0ai(a,i)
t1ai(a,i) += e1ai(a,i)
t1ai(a,i) *= c1
tai(a,i) = t1ai(a,i)
#
t1ai(a,i) = d1ai(a,i)
t1ai(a,i) += e2ai(a,i)
t1ai(a,i) *= c2
tai(a,i) += t1ai(a,i)
#
t1ai(a,i) = d2ai(a,i)
t1ai(a,i) += e3ai(a,i)
t1ai(a,i) *= c3
tai(a,i) += t1ai(a,i)
#
t1ai(a,i) = d3ai(a,i)
t1ai(a,i) += e4ai(a,i)
t1ai(a,i) *= c4
tai(a,i) += t1ai(a,i)
#
t1ai(a,i) = d4ai(a,i)
t1ai(a,i) += e5ai(a,i)
t1ai(a,i) *= c5
tai(a,i) += t1ai(a,i)
#
PUT Paiold_aa(a,i) = tai(a,i)
#
ENDPARDO a, i
#
# Form tbj_old.
# -------------
#
PARDO b, j
#
# REQUEST/GET amplitude data from all 3 previous iterations.
# ----------------------------------------------------------
#
GET D0bj(b,j) # kiter-5 amplitudes
GET D1bj(b,j) # kiter-4 amplitudes
GET D2bj(b,j) # kiter-3 amplitudes
GET D3bj(b,j) # kiter-2 amplitudes
GET D4bj(b,j) # kiter-1 amplitudes
#
GET e1bj(b,j) # kiter-4 amplitudes
GET e2bj(b,j) # kiter-3 amplitudes
GET e3bj(b,j) # kiter-2 amplitudes
GET e4bj(b,j) # kiter-1 amplitudes
GET e5bj(b,j) # kiter-0 amplitudes
#
# Compute contributions to updated amplitudes --> tbj_old.
# --------------------------------------------------------
#
t1bj(b,j) = d0bj(b,j)
t1bj(b,j) += e1bj(b,j)
t1bj(b,j) *= c1
tbj(b,j) = t1bj(b,j)
#
t1bj(b,j) = d1bj(b,j)
t1bj(b,j) += e2bj(b,j)
t1bj(b,j) *= c2
tbj(b,j) += t1bj(b,j)
#
t1bj(b,j) = d2bj(b,j)
t1bj(b,j) += e3bj(b,j)
t1bj(b,j) *= c3
tbj(b,j) += t1bj(b,j)
#
t1bj(b,j) = d3bj(b,j)
t1bj(b,j) += e4bj(b,j)
t1bj(b,j) *= c4
tbj(b,j) += t1bj(b,j)
#
t1bj(b,j) = d4bj(b,j)
t1bj(b,j) += e5bj(b,j)
t1bj(b,j) *= c5
tbj(b,j) += t1bj(b,j)
#
PUT Paiold_bb(b,j) = tbj(b,j)
#
ENDPARDO b, j
#
execute sip_barrier
#
ENDPROC DIIS4
# -------------
#
# ---------------------------------------------------------------------------------
#
# ---------------------------------------------------------------------------------
#
PROC MOVE4
# ----------
#
# 0 --> 1
# ---------------------------------------------------
#
PARDO a, i
GET d1ai(a,i)
PUT d0ai(a,i) = d1ai(a,i)
ENDPARDO a, i
#
PARDO b, j
GET d1bj(b,j)
PUT d0bj(b,j) = d1bj(b,j)
ENDPARDO b, j
#
execute sip_barrier
#
# 2 --> 1
# ---------------------------------------------------
#
PARDO a, i
GET e2ai(a,i)
PUT e1ai(a,i) = e2ai(a,i)
ENDPARDO a, i
#
PARDO b, j
GET e2bj(b,j)
PUT e1bj(b,j) = e2bj(b,j)
ENDPARDO b, j
#
PARDO a, i
GET d2ai(a,i)
PUT d1ai(a,i) = d2ai(a,i)
ENDPARDO a, i
#
PARDO b, j
GET d2bj(b,j)
PUT d1bj(b,j) = d2bj(b,j)
ENDPARDO b, j
#
execute sip_barrier
#
# 3 --> 2
# ---------------------------------------------------
#
PARDO a, i
GET e3ai(a,i)
PUT e2ai(a,i) = e3ai(a,i)
ENDPARDO a, i
#
PARDO b, j
GET e3bj(b,j)
PUT e2bj(b,j) = e3bj(b,j)
ENDPARDO b, j
#
PARDO a, i
GET d3ai(a,i)
PUT d2ai(a,i) = d3ai(a,i)
ENDPARDO a, i
#
PARDO b, j
GET d3bj(b,j)
PUT d2bj(b,j) = d3bj(b,j)
ENDPARDO b, j
#
execute sip_barrier
#
# 4 --> 3
# ---------------------------------------------------
#
PARDO a, i
GET e4ai(a,i)
PUT e3ai(a,i) = e4ai(a,i)
ENDPARDO a, i
#
PARDO b, j
GET e4bj(b,j)
PUT e3bj(b,j) = e4bj(b,j)
ENDPARDO b, j
#
PARDO a, i
GET d4ai(a,i)
PUT d3ai(a,i) = d4ai(a,i)
ENDPARDO a, i
#
PARDO b, j
GET d4bj(b,j)
PUT d3bj(b,j) = d4bj(b,j)
ENDPARDO b, j
#
execute sip_barrier
#
# 5 --> 4
# ---------------------------------------------------
#
PARDO a, i
GET e5ai(a,i)
PUT e4ai(a,i) = e5ai(a,i)
ENDPARDO a, i
#
PARDO b, j
GET e5bj(b,j)
PUT e4bj(b,j) = e5bj(b,j)
ENDPARDO b, j
#
PARDO a, i
GET Paiold_aa(a,i)
PUT d4ai(a,i) = Paiold_aa(a,i)
ENDPARDO a, i
#
PARDO b, j
GET Paiold_bb(b,j)
PUT d4bj(b,j) = Paiold_bb(b,j)
ENDPARDO b, j
#
execute sip_barrier
#
ENDPROC MOVE4
# -------------
#
##############################################################################
#
# START OF MAIN PROGRAM
#
##############################################################################
#
# Transform integrals
# -------------------
#
execute sip_barrier
CALL TRAN_UHF
allocate La(*,*)
allocate Lb(*,*)
DO mu
DO p
La(p,mu) = ca(mu,p)
ENDDO p
ENDDO mu
DO mu
DO q
Lb(q,mu) = cb(mu,q)
ENDDO q
ENDDO mu
#
# Create one-particle arrays
# --------------------------
#
CREATE Pij_aa
CREATE Pij_bb
CREATE Pab_aa
CREATE Pab_bb
CREATE Lai_aa
CREATE Lai_bb
CREATE Painew_aa
CREATE Paiold_aa
CREATE Painew_bb
CREATE Paiold_bb
CREATE Wab_aa
CREATE Wab_bb
CREATE Wij_aa
CREATE Wij_bb
CREATE Wai_aa
CREATE Wai_bb
CREATE P2_ao
CREATE P2A_ao
CREATE P2B_ao
CREATE W2_ao
CREATE Paa_ao
CREATE Pbb_ao
CREATE Yxi
CREATE Yxj
CREATE WHFa
CREATE WHFb
CREATE DHFa
CREATE DHFb
CREATE DHF
execute sip_barrier
#
# Compute the HF contribution to the weighted density matrix
# ----------------------------------------------------------
#
call WHFDENS
execute sip_barrier
#
# First compute the occupied-occupied block of the density correction
# -------------------------------------------------------------------
#
# AAAA/AAAA piece
# ---------------
PARDO a, i, a1, i2
#
REQUEST VSpipi(a,i,a1,i2) i
Tpipi(a,i,a1,i2) = VSpipi(a,i,a1,i2)
execute energy_denominator Tpipi(a,i,a1,i2)
#
DO i1
#
REQUEST VSpipi(a,i1,a1,i2) i1
T1pipi(a,i1,a1,i2) = VSpipi(a,i1,a1,i2)
execute energy_denominator T1pipi(a,i1,a1,i2)
#
Tii(i,i1) = Tpipi(a,i,a1,i2)*T1pipi(a,i1,a1,i2)
Tii(i,i1) *= -0.5
PUT Pij_aa(i,i1) += Tii(i,i1)
#
ENDDO i1
#
ENDPARDO a, i, a1, i2
#
# AABB/AABB piece
# ---------------
PARDO a, i, b, j
#
REQUEST Vpiqj(a,i,b,j) i
Tpiqj(a,i,b,j) = Vpiqj(a,i,b,j)
execute energy_denominator Tpiqj(a,i,b,j)
#
DO i1
#
REQUEST Vpiqj(a,i1,b,j) j
T1piqj(a,i1,b,j) = Vpiqj(a,i1,b,j)
execute energy_denominator T1piqj(a,i1,b,j)
#
Tii(i,i1) = Tpiqj(a,i,b,j)*T1piqj(a,i1,b,j)
Tii(i,i1) *= -1.0
PUT Pij_aa(i,i1) += Tii(i,i1)
#
ENDDO i1
#
ENDPARDO a, i, b, j
#
# BBBB/BBBB piece
# ---------------
PARDO b, j, b1, j2
#
REQUEST VSqjqj(b,j,b1,j2) j
Tqjqj(b,j,b1,j2) = VSqjqj(b,j,b1,j2)
execute energy_denominator Tqjqj(b,j,b1,j2)
#
DO j1
#
REQUEST VSqjqj(b,j1,b1,j2) j1
T1qjqj(b,j1,b1,j2) = VSqjqj(b,j1,b1,j2)
execute energy_denominator T1qjqj(b,j1,b1,j2)
#
Tjj(j,j1) = Tqjqj(b,j,b1,j2)*T1qjqj(b,j1,b1,j2)
Tjj(j,j1) *= -0.5
PUT Pij_bb(j,j1) += Tjj(j,j1)
#
ENDDO j1
#
ENDPARDO b, j, b1, j2
#
# BBAA/BBAA piece
# ---------------
PARDO b, j, a, i
#
REQUEST Vpiqj(a,i,b,j) i
Tpiqj(a,i,b,j) = Vpiqj(a,i,b,j)
execute energy_denominator Tpiqj(a,i,b,j)
#
DO j1
#
REQUEST Vpiqj(a,i,b,j1) i
T1piqj(a,i,b,j1) = Vpiqj(a,i,b,j1)
execute energy_denominator T1piqj(a,i,b,j1)
#
Tjj(j,j1) = Tpiqj(a,i,b,j)*T1piqj(a,i,b,j1)
Tjj(j,j1) *= -1.0
PUT Pij_bb(j,j1) += Tjj(j,j1)
#
ENDDO j1
#
ENDPARDO b, j, a, i
#
# Done compute the occupied-occupied block of the density correction
# -------------------------------------------------------------------
#
# Compute the virtual-virtual block of the density correction
# -------------------------------------------------------------------
#
# AAAA/AAAA piece
# ---------------
#
PARDO a, a2, i, i1
#
REQUEST VSpipi(a,i,a2,i1) i
Tpipi(a,i,a2,i1) = VSpipi(a,i,a2,i1)
execute energy_denominator Tpipi(a,i,a2,i1)
#
DO a1
#
REQUEST VSpipi(a1,i,a2,i1) i
T1pipi(a1,i,a2,i1) = VSpipi(a1,i,a2,i1)
execute energy_denominator T1pipi(a1,i,a2,i1)
#
Taa(a,a1) = Tpipi(a,i,a2,i1)*T1pipi(a1,i,a2,i1)
Taa(a,a1) *= 0.5
PUT Pab_aa(a,a1) += Taa(a,a1)
#
ENDDO a1
#
ENDPARDO a, a2, i, i1
#
# AABB/AABB piece
# ---------------
#
PARDO a, b, i, j
#
REQUEST Vpiqj(a,i,b,j) i
Tpiqj(a,i,b,j) = Vpiqj(a,i,b,j)
execute energy_denominator Tpiqj(a,i,b,j)
#
DO a1
#
REQUEST Vpiqj(a1,i,b,j) i
T1piqj(a1,i,b,j) = Vpiqj(a1,i,b,j)
execute energy_denominator T1piqj(a1,i,b,j)
#
Taa(a,a1) = Tpiqj(a,i,b,j)*T1piqj(a1,i,b,j)
PUT Pab_aa(a,a1) += Taa(a,a1)
#
ENDDO a1
#
ENDPARDO a, b, i, j
#
# BBBB/BBBB piece
# ---------------
#
PARDO b, b2, j, j1
#
REQUEST VSqjqj(b,j,b2,j1) j
Tqjqj(b,j,b2,j1) = VSqjqj(b,j,b2,j1)
execute energy_denominator Tqjqj(b,j,b2,j1)
#
DO b1
#
REQUEST VSqjqj(b1,j,b2,j1) j
T1qjqj(b1,j,b2,j1) = VSqjqj(b1,j,b2,j1)
execute energy_denominator T1qjqj(b1,j,b2,j1)
#
Tbb(b,b1) = Tqjqj(b,j,b2,j1)*T1qjqj(b1,j,b2,j1)
Tbb(b,b1) *= 0.5
PUT Pab_bb(b,b1) += Tbb(b,b1)
#
ENDDO b1
#
ENDPARDO b, b2, j, j1
#
# BBAA/BBAA piece
# ---------------
#
PARDO b, a, j, i
#
REQUEST Vpiqj(a,i,b,j) j
Tpiqj(a,i,b,j) = Vpiqj(a,i,b,j)
execute energy_denominator Tpiqj(a,i,b,j)
#
DO b1
#
REQUEST Vpiqj(a,i,b1,j) j
T1piqj(a,i,b1,j) = Vpiqj(a,i,b1,j)
execute energy_denominator T1piqj(a,i,b1,j)
#
Tbb(b,b1) = Tpiqj(a,i,b,j)*T1piqj(a,i,b1,j)
PUT Pab_bb(b,b1) += Tbb(b,b1)
#
ENDDO b1
#
ENDPARDO b, a, j, i
#
# End compute the virtual-virtual block of the density correction
# -------------------------------------------------------------------
execute sip_barrier
#
# Backtransform Pab to be used in the 'direct' contribution to Lai
# ----------------------------------------------------------------
#
# Transform Pab_aa
# ----------------
PARDO a, a1
#
GET Pab_aa(a,a1)
#
DO mu
#
Ixa(mu,a1) = Pab_aa(a,a1)*ca(mu,a)
#
DO nu
#
Ixx(mu,nu) = Ixa(mu,a1)*ca(nu,a1)
PUT Paa_ao(mu,nu) += Ixx(mu,nu)
#
ENDDO nu
#
ENDDO mu
#
ENDPARDO a, a1
#
# Transform Pab_bb
# ----------------
PARDO b, b1
#
GET Pab_bb(b,b1)
#
DO mu
#
Ixb(mu,b1) = Pab_bb(b,b1)*cb(mu,b)
#
DO nu
#
Ixx(mu,nu) = Ixb(mu,b1)*cb(nu,b1)
PUT Pbb_ao(mu,nu) += Ixx(mu,nu)
#
ENDDO nu
#
ENDDO mu
#
ENDPARDO b, b1
execute sip_barrier
#
# Compute the right-hand side of Eq. 10 --> Lai_aa
# ------------------------------------------------
#
# Compte the 'direct' contributions
# ---------------------------------
CALL LAIAO1
#
# Second-term
# -----------
#
PARDO a, a1, i1, i2
#
REQUEST VSpipi(a,i1,a1,i2) i1
Tpipi(a,i1,a1,i2) = VSpipi(a,i1,a1,i2)
execute energy_denominator Tpipi(a,i1,a1,i2)
#
DO i
#
REQUEST VSpipi(a1,i2,i,i1) i
Tai(a,i) = Tpipi(a,i1,a1,i2)*VSpipi(a1,i2,i,i1)
Tai(a,i) *= 0.5
PUT Lai_aa(a,i) += Tai(a,i)
#
ENDDO i
#
ENDPARDO a, a1, i1, i2
#
PARDO a, b, i1, j
#
REQUEST Vpiqj(a,i1,b,j) j
Tpiqj(a,i1,b,j) = Vpiqj(a,i1,b,j)
execute energy_denominator Tpiqj(a,i1,b,j)
#
DO i
#
REQUEST Vpiqj(i,i1,b,j) i
Tai(a,i) = Tpiqj(a,i1,b,j)*Vpiqj(i,i1,b,j)
PUT Lai_aa(a,i) += Tai(a,i)
#
ENDDO i
#
ENDPARDO a, b, i1, j
#
# Fourth-term
# -----------
#
PARDO i, i1, i2, a
#
REQUEST VSpipi(a,i,i1,i2) i
GET Pij_aa(i1,i2)
#
Tai(a,i) = VSpipi(a,i,i1,i2)*Pij_aa(i1,i2)
Tai(a,i) *= -1.0
PUT Lai_aa(a,i) += Tai(a,i)
#
ENDPARDO i, i1, i2, a
#
PARDO i, j, j1, a
#
REQUEST Vpiqj(a,i,j,j1) i
GET Pij_bb(j,j1)
#
Tai(a,i) = Vpiqj(a,i,j,j1)*Pij_bb(j,j1)
Tai(a,i) *= -1.0
PUT Lai_aa(a,i) += Tai(a,i)
#
ENDPARDO i, j, j1, a
#
# Compute the right-hand side of Eq. 10 --> Lai_bb
# ------------------------------------------------
#
# Second-term
# -----------
#
PARDO b, b1, j1, j2
#
REQUEST VSqjqj(b,j1,b1,j2) j1
Tqjqj(b,j1,b1,j2) = VSqjqj(b,j1,b1,j2)
execute energy_denominator Tqjqj(b,j1,b1,j2)
#
DO j
#
REQUEST VSqjqj(b1,j2,j,j1) j
Tbj(b,j) = Tqjqj(b,j1,b1,j2)*VSqjqj(b1,j2,j,j1)
Tbj(b,j) *= 0.5
PUT Lai_bb(b,j) += Tbj(b,j)
#
ENDDO j
#
ENDPARDO b, b1, j1, j2
#
PARDO a, b, i, j1
#
REQUEST Vpiqj(a,i,b,j1) i
Tpiqj(a,i,b,j1) = Vpiqj(a,i,b,j1)
execute energy_denominator Tpiqj(a,i,b,j1)
#
DO j
#
REQUEST Vpiqj(a,i,j,j1) i
Tbj(b,j) = Tpiqj(a,i,b,j1)*Vpiqj(a,i,j,j1)
PUT Lai_bb(b,j) += Tbj(b,j)
#
ENDDO j
#
ENDPARDO a, b, i, j1
#
# Fourth-term
# -----------
#
PARDO j, j1, j2, b
#
REQUEST VSqjqj(b,j,j1,j2) j
GET Pij_bb(j1,j2)
#
Tbj(b,j) = VSqjqj(b,j,j1,j2)*Pij_bb(j1,j2)
Tbj(b,j) *= -1.0
PUT Lai_bb(b,j) += Tbj(b,j)
#
ENDPARDO j, j1, j2, b
#
PARDO j, i, i1, b
#
REQUEST Vpiqj(i,i1,b,j) i
GET Pij_aa(i,i1)
#
Tbj(b,j) = Vpiqj(i,i1,b,j)*Pij_aa(i,i1)
Tbj(b,j) *= -1.0
PUT Lai_bb(b,j) += Tbj(b,j)
#
ENDPARDO j, i, i1, b
#
# Done compute the right-hand side of Eq. 10 --> Lai_bb
# -----------------------------------------------------
#
execute sip_barrier
#
# Compute the occupied-virtual block of correlated density
# --------------------------------------------------------
#
# Create error arrays used in the DIIS procedure.
# -----------------------------------------------
#
CREATE D0ai
CREATE D1ai
CREATE D2ai
CREATE D3ai
CREATE D4ai
#
CREATE D0bj
CREATE D1bj
CREATE D2bj
CREATE D3bj
CREATE D4bj
#
CREATE e1ai
CREATE e2ai
CREATE e3ai
CREATE e4ai
CREATE e5ai
#
CREATE e1bj
CREATE e2bj
CREATE e3bj
CREATE e4bj
CREATE e5bj
#
# Get initial guess and form a few intermediates
# ----------------------------------------------
#
eold = 0.0
esum = 0.0
ecrit = 0.0000000001
#
PARDO a, a1, i, i1
#
REQUEST VSpipi(a,i,a1,i1) i
REQUEST Vaaii(a,a1,i1,i) i
REQUEST Viaai(i,a,a1,i1) i
#
Tpipi(a,i,a1,i1) = VSpipi(a,i,a1,i1)
T1pipi(a,i,a1,i1) = Vaaii(a,a1,i1,i)
T2pipi(a,i,a1,i1) = Viaai(i,a,a1,i1)
#
Tpipi(a,i,a1,i1) -= T1pipi(a,i,a1,i1)
Tpipi(a,i,a1,i1) += T2pipi(a,i,a1,i1)
PREPARE ASpipi(a,i,a1,i1) = Tpipi(a,i,a1,i1)
#
ENDPARDO a, a1, i, i1
#
PARDO a, b, i, j
#
REQUEST Vpiqj(a,i,b,j) j
REQUEST Viabj(i,a,b,j) j
#
Tpiqj(a,i,b,j) = Vpiqj(a,i,b,j)
T2piqj(a,i,b,j) = Viabj(i,a,b,j)
#
Tpiqj(a,i,b,j) += T2piqj(a,i,b,j)
#
PREPARE Apiqj(a,i,b,j) = Tpiqj(a,i,b,j)
#
ENDPARDO a, b, i, j
#
PARDO b, b1, j, j1
#
REQUEST VSqjqj(b,j,b1,j1) j
REQUEST Vbbjj(b,b1,j1,j) j
REQUEST Vjbbj(j,b,b1,j1) j
#
Tqjqj(b,j,b1,j1) = VSqjqj(b,j,b1,j1)
T1qjqj(b,j,b1,j1) = Vbbjj(b,b1,j1,j)
T2qjqj(b,j,b1,j1) = Vjbbj(j,b,b1,j1)
#
Tqjqj(b,j,b1,j1) -= T1qjqj(b,j,b1,j1)
Tqjqj(b,j,b1,j1) += T2qjqj(b,j,b1,j1)
PREPARE ASqjqj(b,j,b1,j1) = Tqjqj(b,j,b1,j1)
#
ENDPARDO b, b1, j, j1
#
PARDO a, i
#
GET Lai_aa(a,i)
Tai(a,i) = Lai_aa(a,i)
Tai(a,i) *= -1.0
execute energy_denominator Tai(a,i)
etemp = Tai(a,i)*Tai(a,i)
esum += etemp
PUT Paiold_aa(a,i) = Tai(a,i)
#
ENDPARDO a, i
#
PARDO b, j
#
GET Lai_bb(b,j)
Tbj(b,j) = Lai_bb(b,j)
Tbj(b,j) *= -1.0
execute energy_denominator Tbj(b,j)
etemp = Tbj(b,j)*Tbj(b,j)
esum += etemp
PUT Paiold_bb(b,j) = Tbj(b,j)
#
ENDPARDO b, j
#
execute sip_barrier
execute server_barrier
collective eold += esum
#
# Done initial guess
# ------------------
#
# Start iterations
# ----------------
#
DO kiter
#
PARDO a, i
#
GET Lai_aa(a,i)
Tai(a,i) = Lai_aa(a,i)
Tai(a,i) *= -1.0
PUT Painew_aa(a,i) += Tai(a,i)
#
ENDPARDO a, i
#
PARDO b, j
#
GET Lai_bb(b,j)
Tbj(b,j) = Lai_bb(b,j)
Tbj(b,j) *= -1.0
PUT Painew_bb(b,j) += Tbj(b,j)
#
ENDPARDO b, j
#
PARDO a, a1, i, i1
#
#REQUEST VSpipi(a,i,a1,i1) i
#REQUEST Vaaii(a,a1,i1,i) i
#REQUEST Viaai(i,a,a1,i1) i
REQUEST ASpipi(a,i,a1,i1) i
GET Paiold_aa(a1,i1)
#
#Tpipi(a,i,a1,i1) = VSpipi(a,i,a1,i1)
#T1pipi(a,i,a1,i1) = Vaaii(a,a1,i1,i)
#T2pipi(a,i,a1,i1) = Viaai(i,a,a1,i1)
#
#Tpipi(a,i,a1,i1) -= T1pipi(a,i,a1,i1)
#Tpipi(a,i,a1,i1) += T2pipi(a,i,a1,i1)
#
Tai(a,i) = ASpipi(a,i,a1,i1)*Paiold_aa(a1,i1)
PUT Painew_aa(a,i) += Tai(a,i)
#
ENDPARDO a, a1, i, i1
#
PARDO a, b, i, j
#
#REQUEST Vpiqj(a,i,b,j) j
#REQUEST Viabj(i,a,b,j) j
REQUEST Apiqj(a,i,b,j) j
GET Paiold_bb(b,j)
GET Paiold_aa(a,i)
#
#Tpiqj(a,i,b,j) = Vpiqj(a,i,b,j)
#T2piqj(a,i,b,j) = Viabj(i,a,b,j)
#Tpiqj(a,i,b,j) += T2piqj(a,i,b,j)
#
Tai(a,i) = Apiqj(a,i,b,j)*Paiold_bb(b,j)
PUT Painew_aa(a,i) += Tai(a,i)
#
#Tpiqj(a,i,b,j) = Vpiqj(a,i,b,j)
#T2piqj(a,i,b,j) = Viabj(i,a,b,j)
#Tpiqj(a,i,b,j) += T2piqj(a,i,b,j)
#
Tbj(b,j) = Apiqj(a,i,b,j)*Paiold_aa(a,i)
PUT Painew_bb(b,j) += Tbj(b,j)
#
ENDPARDO a, b, i, j
#
PARDO b, b1, j, j1
#
#REQUEST VSqjqj(b,j,b1,j1) j
#REQUEST Vbbjj(b,b1,j1,j) j
#REQUEST Vjbbj(j,b,b1,j1) j
REQUEST ASqjqj(b,j,b1,j1) j
GET Paiold_bb(b1,j1)
#
#Tqjqj(b,j,b1,j1) = VSqjqj(b,j,b1,j1)
#T1qjqj(b,j,b1,j1) = Vbbjj(b,b1,j1,j)
#T2qjqj(b,j,b1,j1) = Vjbbj(j,b,b1,j1)
#
#Tqjqj(b,j,b1,j1) -= T1qjqj(b,j,b1,j1)
#Tqjqj(b,j,b1,j1) += T2qjqj(b,j,b1,j1)
#
Tbj(b,j) = ASqjqj(b,j,b1,j1)*Paiold_bb(b1,j1)
PUT Painew_bb(b,j) += Tbj(b,j)
#
ENDPARDO b, b1, j, j1
#
execute sip_barrier
#
# Update error vector for diis
# ----------------------------
#
CALL UPDATE_PAI
#
execute sip_barrier
#
esum = 0.0
enew = 0.0
PARDO a, i
#
GET Painew_aa(a,i)
Tai(a,i) = Painew_aa(a,i)
execute energy_denominator Tai(a,i)
PUT Paiold_aa(a,i) = Tai(a,i)
etemp = Painew_aa(a,i)*Painew_aa(a,i)
esum += etemp
Tai(a,i) = 0.0
PUT Painew_aa(a,i) = Tai(a,i)
#
ENDPARDO a, i
#
PARDO b, j
#
GET Painew_bb(b,j)
Tbj(b,j) = Painew_bb(b,j)
execute energy_denominator Tbj(b,j)
PUT Paiold_bb(b,j) = Tbj(b,j)
etemp = Painew_bb(b,j)*Painew_bb(b,j)
esum += etemp
Tbj(b,j) = 0.0
PUT Painew_bb(b,j) = Tbj(b,j)
#
ENDPARDO b, j
#
execute sip_barrier
collective enew += esum
#
# Check on convergence
# --------------------
#
IF enew < eold
ediff = eold - enew
IF ediff < ecrit
exit # kiter
ENDIF
ENDIF
#
IF enew > eold
ediff = enew - eold
IF ediff < ecrit
exit # kiter
ENDIF
ENDIF
#
# Reset eold --> enew
# -------------------
#
eold = enew
#
if kiter == 2
#
# Get uptated amplitudes based on DIIS procedure.
# -----------------------------------------------
#
CALL DIIS1
#
endif # kiter == 2
#
if kiter == 3
#
# Get uptated amplitudes based on DIIS procedure.
# -----------------------------------------------
#
CALL DIIS2
#
endif # kiter == 3
#
if kiter == 4
#
# Get uptated amplitudes based on DIIS procedure.
# -----------------------------------------------
#
CALL DIIS3
#
endif # kiter == 4
#
if kiter >= 5
#
# Get uptated amplitudes based on DIIS procedure.
# -----------------------------------------------
#
CALL DIIS4
CALL MOVE4
#
endif # kiter == 5
#
CALL MOVE_PAI
#
ENDDO kiter
#
execute sip_barrier
#
# Delete error arrays used in the DIIS procedure.
# -----------------------------------------------
#
DELETE D0ai
DELETE D1ai
DELETE D2ai
DELETE D3ai
DELETE D4ai
#
DELETE D0bj
DELETE D1bj
DELETE D2bj
DELETE D3bj
DELETE D4bj
#
DELETE e1ai
DELETE e2ai
DELETE e3ai
DELETE e4ai
DELETE e5ai
#
DELETE e1bj
DELETE e2bj
DELETE e3bj
DELETE e4bj
DELETE e5bj
#
# Done compute the occupied-virtual block of correlated density
# -------------------------------------------------------------
#
# Compute the second-order corrections to the energy weighted
# density matrix.
# -----------------------------------------------------------
#
# Compute Wab_aa
# --------------
#
PARDO a, a1, a2
#
GET Pab_aa(a2,a1)
Taa(a,a1) = Pab_aa(a2,a1)*Fock_a(a2,a)
Taa(a,a1) *= -1.0
PUT Wab_aa(a,a1) += Taa(a,a1)
#
ENDPARDO a, a1, a2
#
PARDO a1, a2, i, i1
#
REQUEST VSpipi(a1,i,a2,i1) i
Tpipi(a1,i,a2,i1) = VSpipi(a1,i,a2,i1)
execute energy_denominator Tpipi(a1,i,a2,i1)
#
DO a
#
REQUEST VSpipi(a,i1,a2,i) i
#
Taa(a,a1) = VSpipi(a,i1,a2,i)*Tpipi(a1,i,a2,i1)
Taa(a,a1) *= 0.5
PUT Wab_aa(a,a1) += Taa(a,a1)
#
ENDDO a
#
ENDPARDO a1, a2, i, i1
#
PARDO a1, b, i, j
#
REQUEST Vpiqj(a1,i,b,j) i
Tpiqj(a1,i,b,j) = Vpiqj(a1,i,b,j)
execute energy_denominator Tpiqj(a1,i,b,j)
#
DO a
#
REQUEST Vpiqj(a,i,b,j) i
#
Taa(a,a1) = Vpiqj(a,i,b,j)*Tpiqj(a1,i,b,j)
Taa(a,a1) *= -1.0
PUT Wab_aa(a,a1) += Taa(a,a1)
#
ENDDO a
#
ENDPARDO a1, b, i, j
#
# Done compute Wab_aa
# -------------------
#
# Compute Wab_bb
# --------------
#
PARDO b, b1, b2
#
GET Pab_bb(b2,b1)
Tbb(b,b1) = Pab_bb(b2,b1)*Fock_b(b2,b)
Tbb(b,b1) *= -1.0
PUT Wab_bb(b,b1) += Tbb(b,b1)
#
ENDPARDO b, b1, b2
#
PARDO b1, b2, j, j1
#
REQUEST VSqjqj(b1,j,b2,j1) j1
Tqjqj(b1,j,b2,j1) = VSqjqj(b1,j,b2,j1)
execute energy_denominator Tqjqj(b1,j,b2,j1)
#
DO b
#
REQUEST VSqjqj(b,j1,b2,j) j
#
Tbb(b,b1) = VSqjqj(b,j1,b2,j)*Tqjqj(b1,j,b2,j1)
Tbb(b,b1) *= 0.5
PUT Wab_bb(b,b1) += Tbb(b,b1)
#
ENDDO b
#
ENDPARDO b1, b2, j, j1
#
PARDO a, b1, i, j
#
REQUEST Vpiqj(a,i,b1,j) i
Tpiqj(a,i,b1,j) = Vpiqj(a,i,b1,j)
execute energy_denominator Tpiqj(a,i,b1,j)
#
DO b
#
REQUEST Vpiqj(a,i,b,j) i
#
Tbb(b,b1) = Vpiqj(a,i,b,j)*Tpiqj(a,i,b1,j)
Tbb(b,b1) *= -1.0
PUT Wab_bb(b,b1) += Tbb(b,b1)
#
ENDDO b
#
ENDPARDO a, b1, i, j
#
# PARDO a, a1
#
# GET Wab_aa(a,a1)
# execute dump_block Wab_aa(a,a1)
#
# ENDPARDO a, a1
#
# PARDO b, b1
#
# GET Wab_bb(b,b1)
# execute dump_block Wab_bb(b,b1)
#
# ENDPARDO b, b1
#
# Done compute Wab_bb
# -------------------
#
# Compute Wij_aa
# --------------
#
# Second-term in Eq. 12
# ---------------------
PARDO i, i1, i2
#
GET Pij_aa(i2,i1)
#
T1ii(i,i1) = Pij_aa(i2,i1)*Fock_a(i2,i)
T1ii(i,i1) *= -1.0
PUT Wij_aa(i,i1) += T1ii(i,i1)
#
ENDPARDO i, i1, i2
#
# Fourth-term in Eq. 12
# ---------------------
PARDO a, a1, i1, i2
#
REQUEST VSpipi(a,i1,a1,i2) i1
Tpipi(a,i1,a1,i2) = VSpipi(a,i1,a1,i2)
execute energy_denominator Tpipi(a,i1,a1,i2)
#
DO i
#
REQUEST VSpipi(a,i2,a1,i) i
#
Tii(i,i1) = VSpipi(a,i2,a1,i)*Tpipi(a,i1,a1,i2)
Tii(i,i1) *= 0.5
PUT Wij_aa(i,i1) += Tii(i,i1)
#
ENDDO i
#
ENDPARDO a, a1, i1, i2
#
PARDO a, b, i1, j
#
REQUEST Vpiqj(a,i1,b,j) i1
Tpiqj(a,i1,b,j) = Vpiqj(a,i1,b,j)
execute energy_denominator Tpiqj(a,i1,b,j)
#
DO i
#
REQUEST Vpiqj(a,i,b,j) i
#
Tii(i,i1) = Vpiqj(a,i,b,j)*Tpiqj(a,i1,b,j)
Tii(i,i1) *= -1.0
PUT Wij_aa(i,i1) += Tii(i,i1)
#
ENDDO i
#
ENDPARDO a, b, i1, j
#
# Third-term in Eq. 12
# --------------------
#
# occupied-occupied contribution
# ------------------------------
#
PARDO i, i1, i2, i3
#
REQUEST VSpipi(i,i1,i2,i3) i
GET Pij_aa(i2,i3)
#
Tii(i,i1) = VSpipi(i,i1,i2,i3)*Pij_aa(i2,i3)
Tii(i,i1) *= -1.0
PUT Wij_aa(i,i1) += Tii(i,i1)
#
ENDPARDO i, i1, i2, i3
#
PARDO i, i1, j, j1
#
REQUEST Vpiqj(i,i1,j,j1) i
GET Pij_bb(j,j1)
#
Tii(i,i1) = Vpiqj(i,i1,j,j1)*Pij_bb(j,j1)
Tii(i,i1) *= -1.0
PUT Wij_aa(i,i1) += Tii(i,i1)
#
ENDPARDO i, i1, j, j1
#
# virtual-virtual contribution
# ----------------------------
PARDO i, i1, a, a1
#
REQUEST Vaaii(a,a1,i,i1) i
REQUEST Viaai(i,a1,a,i1) i
GET Pab_aa(a,a1)
#
Tiiaa(i,i1,a,a1) = Vaaii(a,a1,i,i1)
T1iiaa(i,i1,a,a1) = Viaai(i,a1,a,i1)
Tiiaa(i,i1,a,a1) -= T1iiaa(i,i1,a,a1)
#
Tii(i,i1) = Tiiaa(i,i1,a,a1)*Pab_aa(a,a1)
Tii(i,i1) *= -1.0
PUT Wij_aa(i,i1) += Tii(i,i1)
#
ENDPARDO i, i1, a, a1
#
PARDO i, i1, b, b1
#
REQUEST Vbbii(b,b1,i,i1) i
GET Pab_bb(b,b1)
#
Tii(i,i1) = Vbbii(b,b1,i,i1)*Pab_bb(b,b1)
Tii(i,i1) *= -1.0
PUT Wij_aa(i,i1) += Tii(i,i1)
#
ENDPARDO i, i1, b, b1
#
# virtual-occupied contribution --> Needs checked VFL
# -----------------------------
PARDO i, i1, i2, a
#
REQUEST VSpipi(i,i1,a,i2) i
GET Paiold_aa(a,i2)
#
Tii(i,i1) = VSpipi(i,i1,a,i2)*Paiold_aa(a,i2)
Tii(i,i1) *= -1.0
PUT Wij_aa(i,i1) += Tii(i,i1)
T1ii(i1,i) = Tii(i,i1)
PUT Wij_aa(i1,i) += T1ii(i1,i)
#
ENDPARDO i, i1, i2, a
#
PARDO i, i1, j, b
#
REQUEST Vpiqj(i,i1,b,j) i
GET Paiold_bb(b,j)
#
Tii(i,i1) = Vpiqj(i,i1,b,j)*Paiold_bb(b,j)
Tii(i,i1) *= -1.0
PUT Wij_aa(i,i1) += Tii(i,i1)
T1ii(i1,i) = Tii(i,i1)
PUT Wij_aa(i1,i) += T1ii(i1,i)
#
ENDPARDO i, i1, j, b
#
# Compute Wij_bb
# --------------
#
# Second-term in Eq. 12
# ---------------------
PARDO j, j1, j2
#
GET Pij_bb(j2,j1)
#
T1jj(j2,j1) = Pij_bb(j2,j1)
Tjj(j,j1) = T1jj(j2,j1)*Fock_b(j2,j)
Tjj(j,j1) *= -1.0
PUT Wij_bb(j,j1) += Tjj(j,j1)
#
ENDPARDO j, j1, j2
#
# Fourth-term in Eq. 12
# ---------------------
PARDO b, b1, j1, j2
#
REQUEST VSqjqj(b,j1,b1,j2) j1
Tqjqj(b,j1,b1,j2) = VSqjqj(b,j1,b1,j2)
execute energy_denominator Tqjqj(b,j1,b1,j2)
#
DO j
#
REQUEST VSqjqj(b,j2,b1,j) j
#
Tjj(j,j1) = VSqjqj(b,j2,b1,j)*Tqjqj(b,j1,b1,j2)
Tjj(j,j1) *= 0.5
PUT Wij_bb(j,j1) += Tjj(j,j1)
#
ENDDO j
#
ENDPARDO b, b1, j1, j2
#
PARDO a, b, i, j1
#
REQUEST Vpiqj(a,i,b,j1) i
Tpiqj(a,i,b,j1) = Vpiqj(a,i,b,j1)
execute energy_denominator Tpiqj(a,i,b,j1)
#
DO j
#
REQUEST Vpiqj(a,i,b,j) i
#
Tjj(j,j1) = Vpiqj(a,i,b,j)*Tpiqj(a,i,b,j1)
Tjj(j,j1) *= -1.0
PUT Wij_bb(j,j1) += Tjj(j,j1)
#
ENDDO j
#
ENDPARDO a, b, i, j1
#
# Third-term in Eq. 12
# --------------------
#
# occupied-occupied contribution
# ------------------------------
PARDO j, j1, j2, j3
#
REQUEST VSqjqj(j,j1,j2,j3) j
GET Pij_bb(j2,j3)
#
Tjj(j,j1) = VSqjqj(j,j1,j2,j3)*Pij_bb(j2,j3)
Tjj(j,j1) *= -1.0
PUT Wij_bb(j,j1) += Tjj(j,j1)
#
ENDPARDO j, j1, j2, j3
#
PARDO i, i1, j, j1
#
REQUEST Vpiqj(i,i1,j,j1) j
GET Pij_aa(i,i1)
#
Tjj(j,j1) = Vpiqj(i,i1,j,j1)*Pij_aa(i,i1)
Tjj(j,j1) *= -1.0
PUT Wij_bb(j,j1) += Tjj(j,j1)
#
ENDPARDO i, i1, j, j1
#
# virtual-virtual contribution
# ----------------------------
PARDO j, j1, b, b1
#
REQUEST Vbbjj(b,b1,j,j1) j
REQUEST Vjbbj(j,b1,b,j1) j
GET Pab_bb(b,b1)
#
Tjjbb(j,j1,b,b1) = Vbbjj(b,b1,j,j1)
T1jjbb(j,j1,b,b1) = Vjbbj(j,b1,b,j1)
Tjjbb(j,j1,b,b1) -= T1jjbb(j,j1,b,b1)
#
Tjj(j,j1) = Tjjbb(j,j1,b,b1)*Pab_bb(b,b1)
Tjj(j,j1) *= -1.0
PUT Wij_bb(j,j1) += Tjj(j,j1)
#
ENDPARDO j, j1, b, b1
#
PARDO j, j1, a, a1
#
REQUEST Vaajj(a,a1,j,j1) j
GET Pab_aa(a,a1)
#
Tjj(j,j1) = Vaajj(a,a1,j,j1)*Pab_aa(a,a1)
Tjj(j,j1) *= -1.0
PUT Wij_bb(j,j1) += Tjj(j,j1)
#
ENDPARDO j, j1, a, a1
#
# virtual-occupied contribution --> Needs checked VFL
# -----------------------------
PARDO j, j1, j2, b
#
REQUEST VSqjqj(j,j1,b,j2) j
GET Paiold_bb(b,j2)
#
Tjj(j,j1) = VSqjqj(j,j1,b,j2)*Paiold_bb(b,j2)
Tjj(j,j1) *= -1.0
PUT Wij_bb(j,j1) += Tjj(j,j1)
T1jj(j1,j) = Tjj(j,j1)
PUT Wij_bb(j1,j) += T1jj(j1,j)
#
ENDPARDO j, j1, j2, b
#
PARDO j, j1, i, a
#
REQUEST Vpiqj(a,i,j,j1) i
GET Paiold_aa(a,i)
#
Tjj(j,j1) = Vpiqj(a,i,j,j1)*Paiold_aa(a,i)
Tjj(j,j1) *= -1.0
PUT Wij_bb(j,j1) += Tjj(j,j1)
T1jj(j1,j) = Tjj(j,j1)
PUT Wij_bb(j1,j) += T1jj(j1,j)
#
ENDPARDO j, j1, i, a
#
#
# PARDO i, i1
#
# GET Wij_aa(i,i1)
# execute dump_block Wij_aa(i,i1)
#
# ENDPARDO i, i1
#
# PARDO j, j1
#
# GET Wij_bb(j,j1)
# execute dump_block Wij_bb(j,j1)
#
# ENDPARDO j, j1
#
# Compute Wai_aa
# --------------
PARDO a, i, i1
#
GET Paiold_aa(a,i1)
Tai(a,i) = Paiold_aa(a,i1)*Fock_a(i1,i)
Tai(a,i) *= -1.0
PUT Wai_aa(a,i) += Tai(a,i)
#
ENDPARDO a, i, i1
#
PARDO a, a1, i1, i2
#
REQUEST VSpipi(a1,i1,a,i2) i1
Tpipi(a1,i1,a,i2) = VSpipi(a1,i1,a,i2)
execute energy_denominator Tpipi(a1,i1,a,i2)
#
DO i
#
REQUEST VSpipi(a1,i2,i,i1) i
#
Tai(a,i) = Tpipi(a1,i1,a,i2)*VSpipi(a1,i2,i,i1)
Tai(a,i) *= 0.5
PUT Wai_aa(a,i) += Tai(a,i)
#
ENDDO i
#
ENDPARDO a, a1, i1, i2
#
PARDO a, b, j, i2
#
REQUEST Vpiqj(a,i2,b,j) j
Tpiqj(a,i2,b,j) = Vpiqj(a,i2,b,j)
execute energy_denominator Tpiqj(a,i2,b,j)
#
DO i
#
REQUEST Vpiqj(i,i2,b,j) i
#
Tai(a,i) = Tpiqj(a,i2,b,j)*Vpiqj(i,i2,b,j)
Tai(a,i) *= -1.0
PUT Wai_aa(a,i) += Tai(a,i)
#
ENDDO i
#
ENDPARDO a, b, j, i2
#
# Compute Wai_bb
# --------------
PARDO b, j, j1
#
GET Paiold_bb(b,j1)
Tbj(b,j) = Paiold_bb(b,j1)*Fock_b(j1,j)
Tbj(b,j) *= -1.0
PUT Wai_bb(b,j) += Tbj(b,j)
#
ENDPARDO b, j, j1
#
PARDO b, b1, j1, j2
#
REQUEST VSqjqj(b1,j1,b,j2) j1
Tqjqj(b1,j1,b,j2) = VSqjqj(b1,j1,b,j2)
execute energy_denominator Tqjqj(b1,j1,b,j2)
#
DO j
#
REQUEST VSqjqj(b1,j2,j,j1) j
#
Tbj(b,j) = Tqjqj(b1,j1,b,j2)*VSqjqj(b1,j2,j,j1)
Tbj(b,j) *= 0.5
PUT Wai_bb(b,j) += Tbj(b,j)
#
ENDDO j
#
ENDPARDO b, b1, j1, j2
#
PARDO a, b, i, j2
#
REQUEST Vpiqj(a,i,b,j2) i
Tpiqj(a,i,b,j2) = Vpiqj(a,i,b,j2)
execute energy_denominator Tpiqj(a,i,b,j2)
#
DO j
#
REQUEST Vpiqj(a,i,j,j2) i
#
Tbj(b,j) = Tpiqj(a,i,b,j2)*Vpiqj(a,i,j,j2)
Tbj(b,j) *= -1.0
PUT Wai_bb(b,j) += Tbj(b,j)
#
ENDDO j
#
ENDPARDO a, b, i, j2
#
# PARDO a, i
#
# GET Wai_aa(a,i)
# execute dump_block Wai_aa(a,i)
#
# ENDPARDO a, i
#
# PARDO b, j
#
# GET Wai_bb(b,j)
# execute dump_block Wai_bb(b,j)
#
# ENDPARDO b, j
#
# Done compute the second-order corrections to the energy
# weighted density matrix.
# -----------------------------------------------------------
execute sip_barrier
# discard VSpipi
# discard VSqjqj
# discard Vpiqj
# discard Vaaii
# discard Viaai
# discard Vbbjj
# discard Vjbbj
# discard Viabj
# discard Vaajj
# discard Vbbii
#
# Backtransform Ppq --> P2_ao Wpq --> W2_ao
# -----------------------------------------
#
# Transform Pij_aa
# ----------------
PARDO i, i1
#
GET Pij_aa(i,i1)
GET Wij_aa(i,i1)
#
DO mu
#
Ixi(mu,i1) = Pij_aa(i,i1)*ca(mu,i)
Jxi(mu,i1) = Wij_aa(i,i1)*ca(mu,i)
#
DO nu
#
Ixx(mu,nu) = Ixi(mu,i1)*ca(nu,i1)
Jxx(mu,nu) = Jxi(mu,i1)*ca(nu,i1)
PUT P2A_ao(mu,nu) += Ixx(mu,nu)
PUT W2_ao(mu,nu) += Jxx(mu,nu)
#
ENDDO nu
#
ENDDO mu
#
ENDPARDO i, i1
#
# Transform Pij_bb
# ----------------
PARDO j, j1
#
GET Pij_bb(j,j1)
GET Wij_bb(j,j1)
#
DO mu
#
Ixj(mu,j1) = Pij_bb(j,j1)*cb(mu,j)
Jxj(mu,j1) = Wij_bb(j,j1)*cb(mu,j)
#
DO nu
#
Ixx(mu,nu) = Ixj(mu,j1)*cb(nu,j1)
Jxx(mu,nu) = Jxj(mu,j1)*cb(nu,j1)
PUT P2B_ao(mu,nu) += Ixx(mu,nu)
PUT W2_ao(mu,nu) += Jxx(mu,nu)
#
ENDDO nu
#
ENDDO mu
#
ENDPARDO j, j1
#
# Transform Pab_aa
# ----------------
PARDO a, a1
#
GET Pab_aa(a,a1)
GET Wab_aa(a,a1)
#
DO mu
#
Ixa(mu,a1) = Pab_aa(a,a1)*ca(mu,a)
Jxa(mu,a1) = Wab_aa(a,a1)*ca(mu,a)
#
DO nu
#
Ixx(mu,nu) = Ixa(mu,a1)*ca(nu,a1)
Jxx(mu,nu) = Jxa(mu,a1)*ca(nu,a1)
PUT P2A_ao(mu,nu) += Ixx(mu,nu)
PUT W2_ao(mu,nu) += Jxx(mu,nu)
#
ENDDO nu
#
ENDDO mu
#
ENDPARDO a, a1
#
# Transform Pab_bb
# ----------------
PARDO b, b1
#
GET Pab_bb(b,b1)
GET Wab_bb(b,b1)
#
DO mu
#
Ixb(mu,b1) = Pab_bb(b,b1)*cb(mu,b)
Jxb(mu,b1) = Wab_bb(b,b1)*cb(mu,b)
#
DO nu
#
Ixx(mu,nu) = Ixb(mu,b1)*cb(nu,b1)
Jxx(mu,nu) = Jxb(mu,b1)*cb(nu,b1)
PUT P2B_ao(mu,nu) += Ixx(mu,nu)
PUT W2_ao(mu,nu) += Jxx(mu,nu)
#
ENDDO nu
#
ENDDO mu
#
ENDPARDO b, b1
#
# Transform Pai_aa
# ----------------
PARDO a, i
#
GET Paiold_aa(a,i)
GET Wai_aa(a,i)
#
DO mu
#
Ixi(mu,i) = Paiold_aa(a,i)*ca(mu,a)
Jxi(mu,i) = Wai_aa(a,i)*ca(mu,a)
#
DO nu
#
Ixx(mu,nu) = Ixi(mu,i)*ca(nu,i)
Jxx(mu,nu) = Jxi(mu,i)*ca(nu,i)
I1xx(nu,mu) = Ixx(mu,nu)
J1xx(nu,mu) = Jxx(mu,nu)
#
PUT P2A_ao(mu,nu) += Ixx(mu,nu)
PUT W2_ao(mu,nu) += Jxx(mu,nu)
PUT P2A_ao(nu,mu) += I1xx(nu,mu)
PUT W2_ao(nu,mu) += J1xx(nu,mu)
#
ENDDO nu
#
ENDDO mu
#
ENDPARDO a, i
#
# Transform Pai_bb
# ----------------
PARDO b, j
#
GET Paiold_bb(b,j)
GET Wai_bb(b,j)
GET Lai_bb(b,j)
#
DO mu
#
Ixj(mu,j) = Paiold_bb(b,j)*cb(mu,b)
Jxj(mu,j) = Wai_bb(b,j)*cb(mu,b)
#
DO nu
#
Ixx(mu,nu) = Ixj(mu,j)*cb(nu,j)
Jxx(mu,nu) = Jxj(mu,j)*cb(nu,j)
I1xx(nu,mu) = Ixx(mu,nu)
J1xx(nu,mu) = Jxx(mu,nu)
#
PUT P2B_ao(mu,nu) += Ixx(mu,nu)
PUT W2_ao(mu,nu) += Jxx(mu,nu)
PUT P2B_ao(nu,mu) += I1xx(nu,mu)
PUT W2_ao(nu,mu) += J1xx(nu,mu)
#
ENDDO nu
#
ENDDO mu
#
ENDPARDO b, j
#
# Done backtransform Ppq --> P2_ao Wpq --> W2_ao
# ----------------------------------------------
execute sip_barrier
delete Paiold_bb
delete Wai_bb
delete Lai_bb
delete Paiold_aa
delete Wai_aa
delete Lai_aa
execute sip_barrier
#
# Form the HF density
# -------------------
#
CALL HFDENS
execute sip_barrier
#
# Contract the density with the AO basis core Hamiltonian
# -------------------------------------------------------
#
CALL D1TRANS
CALL S1TRANS
#
# Contract the 'two-particle' contributions
# -----------------------------------------
#
CALL D2TRANS
#
# Calculate the mp2 energy and write to JOBARC
# --------------------------------------------
#
#CALL ENERGY
#execute sip_barrier
#
ENDSIAL MBPT2_GRAD_AO2
#
###############################################################################
#