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

molecule dimer {
0 1
 C1     1.098175     1.565263     0.000000  
 H2     0.631555     2.557430     0.000000  
 H3     1.735324     1.460846     0.892506  
 H4     1.735324     1.460846    -0.892506  
 C5     0.000000     0.514178     0.000000  
 O6    -1.206163     0.792232     0.000000  
 N7     0.430929    -0.785971     0.000000  
 H8     1.423106    -0.979872     0.000000  
 C9    -0.527406    -1.884499     0.000000  
 H10   -1.171471    -1.836020     0.890708  
 H11   -1.171471    -1.836020    -0.890708  
 H12    0.025820    -2.832924     0.000000  
--
-1 1
 i   H8   R   N7   180.0   C5   0.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, 4.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 basis_guess true 
set maxiter 1000 
set mom_start 31 
set diis false 
set soscf true 
set soscf_max_iter 15 
set soscf_conv 1.e-4 
set damping_percentage 20.0
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 R in Rvals:
   dimer.R = R
   ecp[R] = energy("mp2", bsse_type = "cp")   ## for new version
   #ecp[R] = cp("df-mp2")                     ## for beta version
   e= ecp[R] * psi_hartree2kcalmol
   psi4.print_out("R, E_int [kcal/mol]: %3.1f   %10.6f\n" % (R, e))


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

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