#  benz and i, R is distance variable for Rvals
memory 1000 mb

molecule dimer {
0 1
 C1     0.000000    0.703980    1.219329
 H2     0.000000    1.251078    2.166931
 C3     0.000000   -0.703980    1.219329
 H4     0.000000   -1.251078    2.166931
 C5     0.000000    1.407959    0.000000
 H6     0.000000    2.502156   -0.000001
 C7     0.000000   -1.407959    0.000000
 H8     0.000000   -2.502156   -0.000001
 C9     0.000000    0.703979   -1.219329
 H10    0.000000    1.251078   -2.166930
 C11    0.000000   -0.703979   -1.219329
 H12    0.000000   -1.251078   -2.166930
--
-1 1
 i   C11   R   C1   A   C7   90.0
}


# set the scan variable of R above, ex: distance, manually/automatically assign
# ----------------------------------------------------------------
def frange(start, stop, step):    ## do not change the section of def frange
	i = start
	while i < stop:
		yield i
		i += step
# ----------------------------------------------------------------
# Rvals=[2.5, 3.0, 4.0]           ## manually assign the variable, ex: distance at 2.5, 3.0 and 4.0 anstrom
Rvals=[]
for j in frange(2.0, 5.0, 0.1):   ##  automatically assign the variable; the endpoint will not be performed
	Rvals.append(round(j,2))  ##  decimal is 2



 
# basis set
# ----------------------------------------------------------------
##set basis <basisset>
basis {

   assign H   def2-tzvppd
   assign Li  def2-tzvppd
   assign Be  def2-tzvppd
   assign C   def2-tzvppd
   assign N   def2-tzvppd
   assign O   def2-tzvppd
   assign F   def2-tzvppd

   assign Na  def2-tzvppd
   assign Mg  def2-tzvppd
   assign P   def2-tzvppd
   assign S   def2-tzvppd
   assign Cl  def2-tzvppd

   assign K   def2-tzvppd
   assign Ca  def2-tzvppd
   assign Br  def2-tzvppd
   assign I   def2-tzvppd

}


set guess sad
set scf_type df 
set maxiter 500
set freeze_core false   ## false for ions
####cp("df-mp2")   ## moved to below


# get the energy at each variable
# ----------------------------------------------------------------------------
# Initialize a blank dictionary of counterpoise corrected energies
# (Need this for the syntax below to work)

ecp = {}
for x in Rvals:                               ## x is the distance between the ion and the center
   dimer.R = math.sqrt(x**2+ 1.4**2 )         ## 1.4 is the center on the ring to the carbon
   ang= math.atan(x/1.4)                      ## Return the arc tangent in radians
   dimer.A = math.degrees(ang)                ## Convert angle x from radians to degrees

   ecp[x] = energy("mp2", bsse_type = "cp")   ## for new version
   #ecp[x] = cp("df-mp2")                     ## for beta version
   e= ecp[x] * psi_hartree2kcalmol
   psi4.print_out("X, E_int [kcal/mol]: %3.1f   %10.6f\n" % (x, e))


# print out all the energies
# ----------------------------------------------------------------------------
psi4.print_out("\n")
psi4.print_out("CP-corrected interaction energies\n\n")
psi4.print_out("        X [Ang]         E_int [kcal/mol]             \n")
psi4.print_out("-----------------------------------------------------\n")

for x in Rvals:
   e = ecp[x] * psi_hartree2kcalmol
   psi4.print_out("        %3.1f            %10.6f\n" % (x, e))