***, correlated electrostatics (CO)_2
memory,10,m

! Thresholds should be set somewhat tighter than defaults 
! if you calculate intermolecular interaction energies:
gthresh,twoint=1.e-16,energy=1.e-9,coeff=1.e-7

core,0

! lower dimer symmetry to the symmetry of monomer with ghost basis on another monomer

set,zsymel=x
geometry={
Q1
Q2,Q1,R
c2,Q2,xc,Q1,thetab
c1,Q1,xc,Q2,180-thetaa,c2,180
o1,Q1,xo,Q2,thetaa,c2,0
o2,Q2,xo,Q1,180-thetab,o1,180
Be,Q1,R/2,c1,180-thetaa,c2,0 ! dummy atom
}

rco=2.132 bohr

mc=12.00000
mo=15.994915
mco=mc+mo

xo=mc*rco/mco
xc=mo*rco/mco

R=8.20 bohr

thetaa=134.23
thetab=45.77


! Nuclear repulsion for dimer
dummy,be
basis=sto-3g
hf;maxit,0
enucab=enuc

basis=vdz ! small test basis

! Choose the method(s):
$method=[mp2,ccsd,qcisd,qcisd(t),xccsd,x3resp,ccsdnor]

! ccsdnor - CCSD without orbital relaxation
! XCCSD and X3RESP - see J. Chem. Phys. 125, 184109 (2006)


! For the final table:
Term=['E(100)','E(1c0)','E(10c)','E(1cc)','E(1)(total)','E(1)(corr)']

do i=1,#method

jset=1
if($method(i).eq.'xccsd'.or.$method(i).eq.'XCCSD')then
jset=3
endif

! A MONOMER  A MONOMER  A MONOMER  A MONOMER  A MONOMER  A MONOMER  A MONOMER  

! Do SCF for the A monomer (of course with the ghost functions on the monomer B)

scfa=2110.2 ! here SCF results for A will be saved

dummy,c2,o2,be
! Accuracy of SCF should also be set tighter than default:
hf;accu,13;save,scfa
enuca=enuc

! Use MATROP to calculate the effective potential for the monomer A:
{matrop
load,EPOT ! read the bare potential for A (v_A)
save,EPOT,5703.2,triang ! save the bare potential for A in the record 5703.2 
! (it will be overwritten in the calculations for the monomer B)
load,DEN,DEN,scfa  ! read the SCF density for A
coul,J,DEN ! calculate the Coulomb operator with the SCF A density
add,WA,1,EPOT,2,J ! calculate the effective potential for A
save,WA,5701.2,triang ! and save it in the record 5701.2
}

! Calculate the correlated density for the monomer A
! and save it to the record 5051.2
if($method(i).eq.'xccsd'.or.$method(i).eq.'XCCSD')then
{ccsd
eom,1.1,prop_order=4,xden=1
start,5001.2;save,5001.2
expec
dm,5051.2
}
else if($method(i).eq.'x3resp'.or.$method(i).eq.'X3RESP')then
{ccsd
eom,1.1,xden=1,x3resp=1
start,5001.2;save,5001.2
expec
dm,5051.2
}
else if($method(i).eq.'ccsdnor'.or.$method(i).eq.'CCSDNOR')then
{ccsd;cphf,1
eom,1.1,nonatorb=1
start,5001.2;save,5001.2
expec
dm,5051.2
}
else
{$method(i);cphf,1
eom,1.1,nonatorb=1
start,5001.2;save,5001.2
expec,relax
dm,5051.2
}
endif

! B MONOMER  B MONOMER  B MONOMER  B MONOMER  B MONOMER  B MONOMER  B MONOMER  

dummy,c1,o1,be
scfb=2111.2 ! here SCF results for B will be saved
hf;accu,13;save,scfb

enucb=enuc

!  Use MATROP to calculate the effective potential for the monomer B:
matrop
load,EPOT
load,DEN,DEN,scfb
coul,J,DEN
add,WB,1,EPOT,2,J
save,WB,5702.2,triang

! Calculate the correlated density for the monomer B 
! and save it to the record 5052.2
if($method(i).eq.'xccsd'.or.$method(i).eq.'XCCSD')then
{ccsd
eom,1.1,prop_order=4,xden=1
start,5002.2;save,5002.2
expec
dm,5052.2
}
else if($method(i).eq.'x3resp'.or.$method(i).eq.'X3RESP')then
{ccsd
eom,1.1,xden=1,x3resp=1
start,5002.2;save,5002.2
expec
dm,5052.2
}
else if($method(i).eq.'ccsdnor'.or.$method(i).eq.'CCSDNOR')then
{ccsd;cphf,1
eom,1.1,nonatorb=1
start,5002.2;save,5002.2
expec
dm,5052.2
}
else
{$method(i);cphf,1
eom,1.1,nonatorb=1
start,5002.2;save,5002.2
expec,relax
dm,5052.2}
endif


! CALCULATION of the ELECTROSTATIC ENERGY

inr=enucab-enuca-enucb

! this is the begin of MATROP 

{matrop
load,TOTDENA,den,5051.2,set=jset  ! read the correlated density for A
load,TOTDENB,den,5052.2,set=jset  ! read the correlated density for B

load,WA,triang,5701.2       ! read the effective potential for A
load,WB,triang,5702.2       ! read the effective potential for B
load,DENA,DEN,scfA          ! read the SCF density for A
load,DENB,DEN,scfB          ! read the SCF density for B
add,CDENA,1,TOTDENA,-1,DENA ! calculate the correlated part of density for A
add,CDENB,1,TOTDENB,-1,DENB ! calculate the correlated part of density for B

! calculate E(100)elst (just for check with SAPT)
! in two ways (see below)
load,EPOT                   ! load bare potential for monomer B
trace,e100B,DENA,EPOT       
trace,e100A,DENB,WA
load,EPOT,triang,5703.2     ! load bare potential for monomer A
trace,e100A1,DENB,EPOT
trace,e100B1,DENA,WB

! (just for check: the trace of the SCF density should be equal to 
! the number of electrons and the trace of the correlated part should be 0)
load,S                      ! load the overlap matrix
trace,na,DENA,S
trace,nb,DENB,S
trace,nap,CDENA,S
trace,nbp,CDENB,S

! The main point - calculation of the correlated electrostatics:

trace,e1n0,CDENA,WB
trace,e10n,CDENB,WA
coul,J,CDENA
trace,e1nn,J,CDENB
}
! this is the end of MATROP 

! Now extract the results:

! the SCF-correlated part (e10n) and the correlated-SCF part (e1n0) are ready

! the correlated-correlated part:
e1nn=2*e1nn

! the SCF-SCF part (calculated in two ways, both should give the same number!)
e100=INR+e100A+e100B
e100p=INR+e100A1+e100B1
if(abs(e100-e100p).gt.1.e-8) then
text, Warning: SCF check failed,$e100,$e100p
endif

! the total correlation contribution:

e1c=e10n+e1n0+e1nn

! the total electrostatic contribution:
e1t=e100+e1c

! hartree - > Kelvin
e100k=e100*tokelvin
e1n0k=e1n0*tokelvin
e10nk=e10n*tokelvin
e1nnk=e1nn*tokelvin
e1ck=e1c*tokelvin
e1tk=e1t*tokelvin

! hartree - > cm-1
e100c=e100*tocm
e1n0c=e1n0*tocm
e10nc=e10n*tocm
e1nnc=e1nn*tocm
e1cc=e1c*tocm
e1tc=e1t*tocm

! Make the final table:

! final results in microhartrees:
n=6
skala=10**n
eau(1)=e100*skala
eau(2)=e1n0*skala
eau(3)=e10n*skala
eau(4)=e1nn*skala
eau(5)=e1t*skala
eau(6)=e1c*skala

eke(1)=e100k
eke(2)=e1n0k
eke(3)=e10nk
eke(4)=e1nnk
eke(5)=e1tk
eke(6)=e1ck

ecm(1)=e100c
ecm(2)=e1n0c
ecm(3)=e10nc
ecm(4)=e1nnc
ecm(5)=e1tc
ecm(6)=e1cc

text,                   FINAL RESULTS

table,Term,eau
digits,,6
title, Electrostatic interaction energy in the $method(i) approximation 
title, energies in hartree/10^$n

table,Term,eke
digits,,6
title, Electrostatic interaction energy in the $method(i) approximation 
title, energies in Kelvin

table,Term,ecm
digits,,6
title, Electrostatic interaction energy in the $method(i) approximation 
title, energies in cm^-1

data,truncate,5051.2
data,truncate,5052.2
enddo