#! Example potential energy surface scan and CP-correction for Ar2
#Based on a Psi4 example from the wiki
#http://psicode.org/psi4manual/master/tutorial.html#analysis-of-intermolecular-interactions

memory 40 gb

molecule dimer {
  Ar 
  Ar 1 R
}

#Define a monomer with a ghost atom. We need this for the CP correction.
#Only one monomer need be defined as, for this system, the other is 
#identical by symmetry.

molecule mono {
   Ar 
  @Ar 1 R
}

Rvals=[3.0, 3.5, 3.75, 4.0, 4.25, 4.5, 4.75, 5.0, 5.5, 6.0]

# For testing only:
#Rvals=[3.0, 3.5]
#molecules=[dimer,mono]


set basis aug-cc-pVDZ
set freeze_core True

# Initialize a blank dictionary of counterpoise corrected energies
# (Need this for the syntax below to work)
e_d_hf  = {}   # for dimer HF energies
e_d_mp2 = {}   # dimer MP2
e_d_cc  = {}   # dimer CCSD(T)
e_m_hf  = {}   # And same for the monomer
e_m_mp2 = {}
e_m_cc  = {}

# The CCSD(T) lines have been commented out below as I could not get that to work.

for R in Rvals:
   dimer.R = R
   energy('ccsd(t)',molecule=dimer)
   e_d_hf[R]  = get_variable('SCF TOTAL ENERGY')
   e_d_mp2[R] = get_variable('MP2 TOTAL ENERGY')
   e_d_cc[R]  = get_variable('CCSD(T) TOTAL ENERGY')
   
   #Remove tmp files etc or else the mono jobs will fail.
   clean()
   
   mono.R = R
   energy('ccsd(t)',molecule=mono)
   e_m_hf[R]  = get_variable('SCF TOTAL ENERGY')
   e_m_mp2[R] = get_variable('MP2 TOTAL ENERGY')
   e_m_cc[R]  = get_variable('CCSD(T) TOTAL ENERGY')
   
   #And another clean or else the next job will fail.
   clean()


psi4.print_out("\n")
psi4.print_out("CP-corrected HF/MP2/CCSD(T)/aug-cc-pVDZ interaction energies\n\n")
psi4.print_out("        R [Ang]     HF         MP2          CCSD(T)  [kJ/mol] \n")
psi4.print_out("--------------------------------------------------------------\n")
for R in Rvals:
   e_hf  = (e_d_hf[R] -2.0*e_m_hf[R])   * psi_hartree2kJmol
   e_mp2 = (e_d_mp2[R] -2.0*e_m_mp2[R]) * psi_hartree2kJmol
   e_cc  = (e_d_cc[R] -2.0*e_m_cc[R])   * psi_hartree2kJmol
   psi4.print_out("        %3.1f     %10.6f     %10.6f     %10.6f  \n" % (R, e_hf,e_mp2,e_cc))