Hi Psi4 guys,
I am trying to calculate the interaction energies for a water dimer for various different separations rx_value between the two molecules by creating a list of values inside rx_values to loop over.
When I set rx_values = [2.10000] my output is as follows (seems fine):
-------------------------------------------------------------------------
SAPT2 energies using aug-cc-pVTZ basis
-------------------------------------------------------------------------
a.u. | kJ/mol
-------------------------------------------------------------------------
Rx | elst exch ind disp total
-------------------------------------------------------------------------
3.96842 -90.72709 207.06986 -10.52348 -26.90564 78.91366
When I set rx_values = [2.20000] my output is follows (seems fine):
-------------------------------------------------------------------------
SAPT2 energies using aug-cc-pVTZ basis
-------------------------------------------------------------------------
a.u. | kJ/mol
-------------------------------------------------------------------------
Rx | elst exch ind disp total
-------------------------------------------------------------------------
4.15740 -66.05530 150.62750 -8.10631 -21.68163 54.78426
However, when I try to calculate the energies for both of these values in a loop the second value and values after that are incorrect, as shown when I set rx_values = [2.10000, 2.20000] as the second value is not the same as previously calculated, this appears to be because when it is ran as a loop the geometry is incorrectly calculated for all separations after the first one but I do not understand why this is happening.
-------------------------------------------------------------------------
SAPT2 energies using aug-cc-pVTZ basis
-------------------------------------------------------------------------
a.u. | kJ/mol
-------------------------------------------------------------------------
Rx | elst exch ind disp total
-------------------------------------------------------------------------
3.96842 -90.72709 207.06986 -10.52348 -26.90564 78.91366
4.15740 -3.61529 5.01672 -0.68081 -2.66735 -1.94674
This is really frustrating and I cannot find out what is wrong, my input file is shown below.
# Input file for calculating various intermolecular energies in Psi4
# Configuration is for a water dimer using SAPT2 as of writing
# Define 4 threads and 8 GB to be used
set_num_threads(4)
memory 8 GB
# Define the geometries of the atoms within the water dimer
# The x-coordinate of the oxygen atom of water2 is at distance rx from oxygen atom of water1
molecule water_dimer {
# water1
O 0.00000 0.00000 -0.28053
H 0.00005 0.76774 0.28052
H 0.00000 -0.76784 0.28041
--
# water2
# rx is the x separation between the oxygen of water1 and the oxygen of water2
# h1rx is the separation between the oxygen of water1 and the first hydrrogen of water2
# h1rx is defined as h1rx = rx + 0.00005
# h2rx is the separation between the oxygen of water1 and the second hydrogen of water2
# h2rx is defined as h2rx = rx
O rx 0.00000 0.28053
H h1rx -0.76774 -0.28052
H h2rx 0.76784 -0.28041
# Define units for molecule geometries to be in Ångströms
units angstrom
}
# Define values of rx to perform calculations for, given in terms of the Bohr radius (atomic units)
rx_values = [2.10000, 2.20000]
# Applying settings
set {
# Choosing a basis set
# Chosen basis family is aug-cc-pVTZ, the MP2 auxilary basis can be specified with aug-cc-pvtz-ri
basis aug-cc-pVTZ
}
# Initialise blank dictionaries for the SAPT electrostatic, exchange, induction, dispersion and total energies
saptelst = {}
saptexch = {}
saptind = {}
saptdisp = {}
sapttotal = {}
# Loop through all values in rx_values
for rx_value in rx_values:
# Set rx for the waterdimer coordinates from the loop
water_dimer.rx = rx_value
# Calculate h1rx and h2rx using rx
water_dimer.h1rx = water_dimer.rx + 0.00005
water_dimer.h2rx = water_dimer.rx
# Calculate energies using 2nd-order SAPT, traditional definition
energy('sapt2')
# Add the relevant values to their dictionaries, units are given in Hartree energy (atomic units)
saptelst[rx_value] = psi4.get_variable("SAPT ELST ENERGY")
saptexch[rx_value] = psi4.get_variable("SAPT EXCH ENERGY")
saptind[rx_value] = psi4.get_variable("SAPT IND ENERGY")
saptdisp[rx_value] = psi4.get_variable("SAPT DISP ENERGY")
sapttotal[rx_value] = psi4.get_variable("SAPT TOTAL ENERGY")
# Remove scratch files, important to run between independent jobs
# Not sure about this yet
clean()
# Output a header for the table of results
psi4.print_out("\n")
psi4.print_out("-------------------------------------------------------------------------\n")
psi4.print_out(" SAPT2 energies using aug-cc-pVTZ basis \n")
psi4.print_out("-------------------------------------------------------------------------\n")
psi4.print_out(" a.u. | kJ/mol \n")
psi4.print_out("-------------------------------------------------------------------------\n")
psi4.print_out(" Rx | elst exch ind disp total \n")
psi4.print_out("-------------------------------------------------------------------------\n")
# Output the results to the table
for rx_value in rx_values:
# Output each value on a line, converting energy units to kJ/mol and distances to Bohr radius, then ending with a $
# Two spaces are given between each value, nine characters with eight decimal places are used
psi4.print_out(" %11.5f" % (rx_value / psi_bohr2angstroms))
psi4.print_out(" %11.5f" % (saptelst[rx_value] * psi_hartree2kJmol))
psi4.print_out(" %11.5f" % (saptexch[rx_value] * psi_hartree2kJmol))
psi4.print_out(" %11.5f" % (saptind[rx_value] * psi_hartree2kJmol))
psi4.print_out(" %11.5f" % (saptdisp[rx_value] * psi_hartree2kJmol))
psi4.print_out(" %11.5f" % (sapttotal[rx_value] * psi_hartree2kJmol))
psi4.print_out(" \n")
If anybody could help me out at all I would greatly appreciate it!
Kind regards,
Charlie
