***,He-H2 correlated electrostatics from Molpro
memory,4,m

! Thresholds should be set somewhat tighter than defaults 
! if you calculate intermolecular interaction energies:
gthresh,twoint=1.e-15,energy=1.e-9,coeff=1.e-6 ! minimal thresholds good for 2002.6 (thanks to Song Huajie <hjsongmoru@tom.com>) - changed on web 24.09.2003
!gthresh,twoint=1.e-15,energy=1.e-7,coeff=1.e-6 ! this works OK with 2002.7

! just a small test basis
basis=avdz

! Set geometry
r=4.0 bohr
rhh=1.4 bohr

geometry={
he,,  0.,   0. ,   r he
h1,,  -rhh/2,   0. ,   0. h 
h2,,  rhh/2,   0. ,   0. h}

! Nuclear repulsion term  (only necessary for E(100))
inr=2*2/sqrt((rhh/2)**2+r**2)

! Choose the method(s):
method=[ccsd] ! available methods: ccsd,qcisd,mp2

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

do i=1,#method

! 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,h1,h2
! Accuracy of SCF should also be set tighter than default:
hf;accu,13;save,scfa

! 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
$method(i);core,0;cphf,1
expec
save,density=5051.2

if($method(i).eq.CCSD.and.$version.lt.2002.5) then
data,copy,5050.2,5051.2 ! with relaxation contribution (better)
data,copy,5040.2,5041.2 ! nonrelaxed (worse, but can be 
! calculated with frozen core)
endif

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

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

!  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
$method(i);core,0;cphf,1
expec
save,density=5052.2

if($method(i).eq.CCSD.and.$version.lt.2002.5) then
data,copy,5050.2,5052.2
data,copy,5040.2,5042.2
endif

! CALCULATION of the ELECTROSTATIC ENERGY

! this is the begin of MATROP 

if($method(i).eq.CCSD.and.$version.lt.2002.5) then

matrop
load,TOTDENA,triang,5051.2  ! read the correlated density for A
load,TOTDENB,triang,5052.2  ! 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

else ! dependent on molpro version

matrop
load,TOTDENA,den,5051.2  ! read the correlated density for A
load,TOTDENB,den,5052.2  ! 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 

endif ! dependent on molpro version

! 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 approximation 
title, energies in hartree/10^$n

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

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

enddo