Psi4: Ar$_2$ interaction energies

In these series of examples we will explore how to best calculate the interaction energies of the Ar dimer system. We will use only MP2 (Moller-Plesset perturbation theory at second-order) and Hartree-Fock (HF) theory. The latter does not lead to binding of rare-gas complexes as these are bound entirely by the van der Waals interaction (dispersion energy) with is the result of electron correlation. MP2 is not a perfect theory but we use it here as it is computationally efficient enough for our purposes.

We will use the Psi4 code here. This is a cutting-edge research-level electronic structure program that is fast, versatile, and which contains a large number of techniques. You will find extensive documentation on the Psi4 website. The code is written in C++, F90 and Python. The input file can use a lot of Python commands - something we will make use of.

Setting up Psi4

In addition to the NWChem set up instructions, do the following.

Include the lines in your .bashrc file:

  # Psi4 paths
  export PSI4_HOME=/opt/ajm/Psi4/current
  export PATH="${PSI4_HOME}/bin:$PATH"
  # Make sure that this command comes after you have defined the SCRATCH variable:
  export PSI_SCRATCH=${SCRATCH}

Then execute it using

  $ source .bashrc

If done correctly, you will be able to find the Psi4 executable using the which command:

  $ which psi4
  /opt/ajm/Psi4/current/bin/psi4

If you get //command not found// then something is wrong and ask your demonstrators for help.

Running Psi4

Psi4 is a parallel code and you need to specify the number of processors it will use using the invocation:

  $ psi4 -n 1 job.psi4 job_basis.out &

where job.psi4 is the Psi4 command file and it will typically use a certain basis. I tend to label my output files using both the job name and the basis like job_basis.out. This makes it clear to me which Psi4 command file I used, and with which basis.

In this class you will always use only 1 processor.

Basis sets

The basis sets we will use are the so called Dunning correlation-consistent basis set. You can see this from the EMSL Basis Set Exchange. Go there and see if you can find the basis sets for Argon.

If you want to learn more on basis sets see this introductory article by a colleague of mine, Dr Mike Towler.

In the examples below you will need to modify the files to use these basis sets:

• cc-pVDZ
• cc-pVTZ
• aug-cc-pVDZ
• aug-cc-pVTZ

These are all Dunning basis sets. The important parts are the 'D' and 'T': these give you an idea of the size of the basis. A double-zeta ('D') basis is smaller than a triple-zeta ('T') basis. Also, the 'aug' means that the basis has been augmented with basis functions that are more slowly decaying (the //diffuse// functions). These tend to be better for van der Waals systems.

No Counterpoise (CP) correction

First of all let's see how the interaction energy can be calculated for a pair of argon atoms at a fixed separation. By definition (see my notes on Intermolecular Interactions) the interaction energy is given by $$

• E_{\rm int} = E[A..B] - E[A] - E[B].

$$ Here $E[A..B]$ is the energy of the complex A..B, and $E[A]$ and $E[B]$ are the energies of the isolated subsystems. Here we have two argon atoms separated by distance $R$, so this becomes $$

• E_{\rm int}(R) = E[{\rm Ar..Ar}] - 2 E[{\rm Ar}]

$$ Therefore all we need to do is calculate the energy of two argon atoms separated by $R$ and subtract twice the energy of one argon. This can be done using the following command file for Psi4. See the comments in the file for explanations of what's going on.

Ar2_MP2_single_noCP.psi4

# Ar2 energy single geometry : MP2 : no counterpoise correction

memory 8 gb

molecule dimer {
  Ar 
  Ar 1  3.75
}

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

molecule mono {
   Ar 
}

molecules=[dimer,mono]

set {
  Print 3
}

# Use basis sets: 
#   cc-pVDZ
#   cc-pVTZ
#   aug-cc-pVDZ
#   aug-cc-pVTZ
#
basis {
   assign aug-cc-pVTZ
}

# This option makes the calculations faster without degrading the
# energies by a significant amount. It causes the inner electrons 
# to be 'frozen', i.e., not considered for the more computationally
# demanding parts of the calculation.
set freeze_core True

# First perform the monomer energy calculation.
 
energy('mp2',molecule=mono)
# Psi4 has variables that store important energies. We extract the
# HF (==SCF) and MP2 energies using these and store them in our own (Python) 
# variables:
e_m_hf  = get_variable('SCF TOTAL ENERGY')
e_m_mp2 = get_variable('MP2 TOTAL ENERGY')

# We need this statement to get Psi4 to clear its settings. 
# This is needed as it currently has various arrays set to the size
# of the previous calculation, but when we perform the Ar2 calculation,
# the system and basis size will change. This will cause the code to 
# fail with a very cryptic error message. If you want to see it, 
# comment out this line and run this code.
clean()


# Now the same for the dimer:
energy('mp2',molecule=dimer)
# Once again we store the HF and MP2 energies in our variables:
e_d_hf  = get_variable('SCF TOTAL ENERGY')
e_d_mp2 = get_variable('MP2 TOTAL ENERGY')

#Remove tmp files etc or else the mono jobs will fail.
clean()

# Psi4 conversion factors
# These constants can be used to convert energy units:
#   psi_hartree2kJmol
#   psi_hartree2wavenumbers
#
# Now define the interaction energy as
#   Eint = E[Ar..Ar] - E[Ar] - E[Ar]
#        = E[Ar..Ar] - 2*E[Ar]
# and convert  units from Hartree (atomic units) to wavenumbers (cm^{-1})
# We do this as the interaction energies of rare gas dimers is usually 
# expressed in these units.

psi4.print_out("\n")
psi4.print_out("no-CP correction HF/MP2 interaction energies\n\n")
psi4.print_out("        R [Bohr]     HF         MP2          [cm^{-1}] \n")
psi4.print_out("-------------------------------------------------------\n")
e_hf  = (e_d_hf  -2.0*e_m_hf ) * psi_hartree2wavenumbers
e_mp2 = (e_d_mp2 -2.0*e_m_mp2) * psi_hartree2wavenumbers
psi4.print_out("        %6.3f     %12.6f     %12.6f       \n" % (R, e_hf,e_mp2))

It looks like a lot, but examine it carefully and you will see that most of it makes sense. We have used the //aug-cc-pVTZ// basis here. Run this using

  $ psi4 -n 1 Ar2_MP2_single_noCP.psi4 Ar2_MP2_single_noCP_aTZ_R3.75.out &

Notice I have included the basis set as well as the separation $R=3.75$ Angstrom in the output file name.

The output file will be quite large. Have a look at it and see if you can spot the SCF (same as Hartree-Fock) energies and MP2 energies. At the very end of the file you should see the following:

no-CP correction HF/MP2 interaction energies

        R [Bohr]     HF         MP2          [cm^{-1}] 
------------------------------------------------------
        8.3        89.551891      -110.712792       

'''* Psi4 exiting successfully. Buy a developer a beer!

The interaction energy at the HF level of theory is $89.55 ~ cm^{-1}$ and at the MP2 level it is $-110.71 ~ cm^{-1}$. This is using the //aug-cc-pVTZ// basis.

There are a few things to note:

• With such an input file we get the interaction energy for a single geometry only.
• Also, we have used a single basis. Try out this calculation with each of the four basis sets listed above.

• There is an error in this calculation that originated from the fact that the basis sets we have used are finite in extent in space. This is called the basis set superposition error, or BSSE. For now ignore this error.

Energy scan: no CP correction : Dimer basis

In the previous example we calculated the interaction energy for a single separation $R$. This is not sufficient and we can gain far more information when we use multiple separations. In this example you will calculate the interaction energy without the counterpoise correction for a range of $R$ values. Here we use a standard basis for the dimer (Ar$_2$) calculation and a monomer basis for Ar. Compute the interaction energies without the CP correction for the set of basis sets listed above.

Ar2_MP2_scan_noCP.psi4

# Ar2 energy scan : MP2 : no counterpoise correction

memory 8 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 
}

# This array contains the list of separations at which the calculations will be performed. 
# Use the longer array only when you are sure the code is working correctly.
# Rvals=[3.0, 3.25, 3.5, 3.65, 3.75, 3.85, 4.0, 4.25, 4.5, 4.75, 5.0, 5.5, 6.0]
Rvals=[3.0, 3.5, 3.75 ]
molecules=[dimer,mono]

set {
  Print 3
}

# Use basis sets: 
#   cc-pVDZ
#   cc-pVTZ
#   aug-cc-pVDZ
#   aug-cc-pVTZ
#
basis {
   assign aug-cc-pVTZ
}

# This option makes the calculations faster without degrading the
# energies by a significant amount. It causes the inner electrons 
# to be 'frozen', i.e., not considered for the more computationally
# demanding parts of the calculation.
set freeze_core True

# Initialize a blank dictionary of counterpoise corrected energies
# (Need this for the syntax below to work)
e_d_hf  = {}
e_d_mp2 = {}
e_d_cc  = {}

# First perform the monomer energy calculation. As we are not using
# the CP correction, this needs to be done only once.
# Also, we do not need to put the energies in arrays.
#
energy('mp2',molecule=mono)
# energy('ccsd(t)',molecule=mono)
e_m_hf  = get_variable('SCF TOTAL ENERGY')
e_m_mp2 = get_variable('MP2 TOTAL ENERGY')
# e_m_cc  = get_variable('CCSD(T) TOTAL ENERGY')

# We need this statement to get Psi4 to clear its settings. 
# This is needed as it currently has various arrays set to the size
# of the previous calculation, but when we perform the Ar2 calculation,
# the system and basis size will change. This will cause the code to 
# fail with a very cryptic error message. If you want to see it, 
# comment out this line and run this code.
clean()


for R in Rvals:
   dimer.R = R
   energy('mp2',molecule=dimer)
   # 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()

# Psi4 conversion factors
# psi_hartree2kJmol
# psi_hartree2wavenumbers

psi4.print_out("\n")
psi4.print_out("no-CP correction HF/MP2/CCSD(T)/aug-cc-pVDZ interaction energies\n\n")
psi4.print_out("        R [Ang]     HF         MP2          CCSD(T)  [cm^{-1}] \n")
psi4.print_out("---------------------------------------------------------------\n")
for R in Rvals:
   e_hf  = (e_d_hf[R]  -2.0*e_m_hf ) * psi_hartree2wavenumbers
   e_mp2 = (e_d_mp2[R] -2.0*e_m_mp2) * psi_hartree2wavenumbers
   # e_cc  = (e_d_cc[R]  -2.0*e_m_cc ) * psi_hartree2wavenumbers
   e_cc  = 0.0
   psi4.print_out("        %6.3f     %12.6f     %12.6f     %12.6f  \n" % (R, e_hf,e_mp2,e_cc))

Run it using

  $ psi4 -n 1 Ar2_MP2_scan_noCP.psi4 Ar2_MP2_scan_noCP_aTZ.out &

This may take a few minutes.

The result summary will be at the end of the output and should look like this

no-CP correction HF/MP2/CCSD(T)/aug-cc-pVDZ interaction energies

        R [Ang]     HF         MP2          CCSD(T)  [cm^{-1}]
---------------------------------------------------------------
        3.0      1421.719436       683.451464         0.000000
        3.5       229.679934       -75.276006         0.000000
        3.8        89.551891      -110.712792         0.000000

'''* Psi4 exiting successfully. Buy a developer a beer!

Ignore the last column. It would normally contant energies calculated using a very high level of theory (CCSD(T)) but we will not use this method as it is quite computationally intensive.

Dimer basis with mid-bonds

Mid-bond functions are basis functions placed between the interacting systems. In our case, this would be mid-way between the two argon atoms. These basis functions are known to significantly accelerate the convergence of the interaction energies. I.e., with the mid-bonds present you can use smaller basis sets and yet achieve results close to those from large basis calculations obtained without the mid-bonds. This is great as the computational cost thereby decreases.

However mid-bond functions pose a massive problem when the counterpoise correction is not applied. And the sole purpose of this part of the exercise is to show you how bad a problem this can be.

Ar2_MP2_scan_noCP_Mb.psi4

# Ar2 energy scan : MP2 : no counterpoise correction
# using Mid-Bonds (MB)
#
# The mid-bond functions accelerate convergence with basis 
# size. But here, when the CP correction is not applied, 
# they wreck havoc with the interaction energies!
# The situation improves somewhat when we move to 
# larger basis sets, so you should see that the aug-cc-pVDZ 
# basis results in nonsense for the interaction energy, but
# aug-cc-pVTZ is slightly better and aug-cc-pVQZ still better.
# But none of these result in interaction energies we can 
# use. 
#

memory 8 gb

molecule dimer {
  Ar 
 @Rn 1 RMb
  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 
}

# Rvals=[3.0, 3.25, 3.5, 3.65, 3.75, 3.85, 4.0, 4.25, 4.5, 4.75, 5.0, 5.5, 6.0]
Rvals=[3.0, 3.5, 3.75 ]
molecules=[dimer,mono]

set {
  Print 3
}

basis {
   assign aug-cc-pVDZ
   assign Rn  mb-set
}

set freeze_core True

# Initialize a blank dictionary of counterpoise corrected energies
# (Need this for the syntax below to work)
e_d_hf  = {}
e_d_mp2 = {}
e_d_cc  = {}

# First perform the monomer energy calculation. As we are not using
# the CP correction, this needs to be done only once.
# Also, we do not need to put the energies in arrays.
#
energy('mp2',molecule=mono)
# energy('ccsd(t)',molecule=mono)
e_m_hf  = get_variable('SCF TOTAL ENERGY')
e_m_mp2 = get_variable('MP2 TOTAL ENERGY')
# e_m_cc  = get_variable('CCSD(T) TOTAL ENERGY')

# We need this statement to get Psi4 to clear its settings. 
# This is needed as it currently has various arrays set to the size
# of the previous calculation, but when we perform the Ar2 calculation,
# the system and basis size will change. This will cause the code to 
# fail with a very cryptic error message. If you want to see it, 
# comment out this line and run this code.
clean()


for R in Rvals:
   dimer.R = R
   dimer.RMb = R/2.0  # position of the midbond
   energy('mp2',molecule=dimer)
   # 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()

# Psi4 conversion factors
# psi_hartree2kJmol
# psi_hartree2wavenumbers

psi4.print_out("\n")
psi4.print_out("no-CP correction HF/MP2/CCSD(T)/aug-cc-pVDZ interaction energies \n\n")
psi4.print_out("        R [Ang]     HF         MP2          CCSD(T)  [cm^{-1}] \n")
psi4.print_out("---------------------------------------------------------------\n")
for R in Rvals:
   e_hf  = (e_d_hf[R]  -2.0*e_m_hf ) * psi_hartree2wavenumbers
   e_mp2 = (e_d_mp2[R] -2.0*e_m_mp2) * psi_hartree2wavenumbers
   # e_cc  = (e_d_cc[R]  -2.0*e_m_cc ) * psi_hartree2wavenumbers
   e_cc  = 0.0
   psi4.print_out("        %6.3f     %12.6f     %12.6f     %12.6f  \n" % (R, e_hf,e_mp2,e_cc))

Counterpoise Corrected interaction energies

To correct for the errors the result with finite basis sets we use what is known as the Boys and Bernadi //counterpoise correction//. The idea is very simple: calculate the energies of the dimer and monomers in the same basis set. That way we avoid, or suppress, any tendency for one of these energies to be inconsistent with the others. So here we do the following:

• For separation //R//,
• compute the dimer energy as before: Ar...Ar
• But compute the monomer energy with basis functions located both on the Argon and //at the position where its partner would be//: Ar...Ar-Ghost. I have indicate the partner with //Ar-Ghost// as it is, in a way, a ghost atom.
• The complication is that now the monomer energy depends on //R// as the Ar-Ghost position varies with //R//.

Here is a command file to perform a CP-corrected calculation of the interaction energy for a fixed geometry. The only change from the uncorrected file (first one on this page) is that the monomer is defined with a ghost atom denoted as //@Ar//:

Ar2_MP2_single_CP.psi4

# Ar2 energy single geometry : MP2 : no counterpoise correction

memory 8 gb
R=3.75

molecule dimer {
  Ar 
  Ar 1  3.75
}

#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  3.75
}

molecules=[dimer,mono]

set {
  Print 3
}

# Use basis sets: 
#   cc-pVDZ
#   cc-pVTZ
#   aug-cc-pVDZ
#   aug-cc-pVTZ
#
basis {
   assign aug-cc-pVTZ
}

# This option makes the calculations faster without degrading the
# energies by a significant amount. It causes the inner electrons 
# to be 'frozen', i.e., not considered for the more computationally
# demanding parts of the calculation.
set freeze_core True

# First perform the monomer energy calculation.
 
energy('mp2',molecule=mono)
# Psi4 has variables that store important energies. We extract the
# HF (==SCF) and MP2 energies using these and store them in our own (Python) 
# variables:
e_m_hf  = get_variable('SCF TOTAL ENERGY')
e_m_mp2 = get_variable('MP2 TOTAL ENERGY')

# We need this statement to get Psi4 to clear its settings. 
# This is needed as it currently has various arrays set to the size
# of the previous calculation, but when we perform the Ar2 calculation,
# the system and basis size will change. This will cause the code to 
# fail with a very cryptic error message. If you want to see it, 
# comment out this line and run this code.
clean()


# Now the same for the dimer:
energy('mp2',molecule=dimer)
# Once again we store the HF and MP2 energies in our variables:
e_d_hf  = get_variable('SCF TOTAL ENERGY')
e_d_mp2 = get_variable('MP2 TOTAL ENERGY')

#Remove tmp files etc or else the mono jobs will fail.
clean()

# Psi4 conversion factors
# These constants can be used to convert energy units:
#   psi_hartree2kJmol
#   psi_hartree2wavenumbers
#
# Now define the interaction energy as
#   Eint = E[Ar..Ar] - E[Ar] - E[Ar]
#        = E[Ar..Ar] - 2*E[Ar]
# and convert  units from Hartree (atomic units) to wavenumbers (cm^{-1})
# We do this as the interaction energies of rare gas dimers is usually 
# expressed in these units.

psi4.print_out("\n")
psi4.print_out("CP corrected HF/MP2 interaction energies\n\n")
psi4.print_out("        R [Bohr]     HF         MP2          [cm^{-1}] \n")
psi4.print_out("-------------------------------------------------------\n")
e_hf  = (e_d_hf  -2.0*e_m_hf ) * psi_hartree2wavenumbers
e_mp2 = (e_d_mp2 -2.0*e_m_mp2) * psi_hartree2wavenumbers
psi4.print_out("        %6.3f     %12.6f     %12.6f       \n" % (R, e_hf,e_mp2))

As before, run this using

  $ psi4 -n 1 Ar2_MP2_single_CP.psi4 Ar2_MP2_single_CP_aTZ_R3.75.out &

The output will contain

CP corrected HF/MP2 interaction energies

        R [Bohr]     HF         MP2          [cm^{-1}]
-------------------------------------------------------
        8.3        96.102802       -82.055480

Energy scan with the CP correction : dimer basis

Here is a Psi4 command file that is very similar to the ones above, but this one performs an energy scan with the counterpoise correction applied.

Ar2_MP2_scan_CP.psi4

# Ar2 : MP2 energy scan : Counterpoise correction applied
# Using mid-bonds (MB) and density-fitting (DF)

#! 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 10 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.25, 3.5, 3.65, 3.75, 3.85, 4.0, 4.25, 4.5, 4.75, 5.0, 5.5, 6.0]
Rvals=[3.0, 3.5, 3.75 ]
molecules=[dimer,mono]


# This is one way of defining a basis set:
# In the other examples we use a different, and more flexible way.
set basis aug-cc-pVTZ

set freeze_core True

# Initialize a blank dictionary of counterpoise corrected energies
# (Need this for the syntax below to work)
e_d_hf  = {}
e_d_mp2 = {}
e_d_cc  = {}
e_m_hf  = {}
e_m_mp2 = {}
e_m_cc  = {}

for R in Rvals:
   dimer.R = R
   energy('mp2',molecule=dimer)
   # 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('mp2',molecule=mono)
   # 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')

   clean()

# Psi4 conversion factors
# psi_hartree2kJmol
# psi_hartree2wavenumbers

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)  [cm^{-1}] \n")
psi4.print_out("---------------------------------------------------------------\n")
for R in Rvals:
   e_hf  = (e_d_hf[R]  -2.0*e_m_hf[R] ) * psi_hartree2wavenumbers
   e_mp2 = (e_d_mp2[R] -2.0*e_m_mp2[R]) * psi_hartree2wavenumbers
   # e_cc  = (e_d_cc[R]  -2.0*e_m_cc[R] ) * psi_hartree2wavenumbers
   e_cc  = 0.0
   psi4.print_out("        %6.3f     %12.6f     %12.6f     %12.6f  \n" % (R, e_hf,e_mp2,e_cc))

Energy scan with CP correction : dimer basis with mid-bonds

The above Psi4 command file works fine, but you may have noticed that there was no clear indication that the interaction energies had converged with basis set. By converged we mean that a large basis would not change the energies by much (or at all). Convergence can be vastly improved if we include the so-called //mid-bond basis sets//. These are basis sets placed in the middle of the interacting species. So in our case we'd have the following

  Ar....Mb....Ar

where Mb indicates the mid-bond set. In the following example we use mid-bond functions. But we also use a numerical device called //density-fitting// to significantly speed up our calculations. This can allow us to use even larger basis sets. The complication is that we now have to define two kinds of basis sets:

basis {
   assign aug-cc-pVDZ
   # Do not alter this one:
   assign Rn  mb-set
}

df_basis_scf {
   assign aug-cc-pVDZ-jkfit
   # Do not alter this one:
   assign Rn mb-set-ri
}

Here the first basis (indicated simply as //basis//) is the main basis set. For technical reasons the Rn is not radon but is actually used for the mid-bond basis. The second basis is //df_basis_scf// and is the density-fitting basis. This needs to be matched with the //basis// definition. So if you alter the main basis to be aug-cc-pVTZ then also change the //df_basis_scf// to use aug-cc-pVTZ-jkfit.

Here is the full command file:

Ar2_MP2_scan_CP_DF_Mb.psi4

# Ar2 : CCSD(T) energy scan : Counterpoise correction applied
# Using mid-bonds (MB) and density-fitting (DF)
#
# Mid-Bonds
# ---------
# These are a small set of basis functions placed between the 
# interacting atoms. Here they are at the site labeled 'Rn'.
# 
# Mid-bonds significantly accelerate the convergence of energies
# with basis set. So you should find interaction energies calculated
# using the mid-bond sets to be closer to the reference energies.
#
# Density-Fitting
# ---------------
# This is a numerical device used to significantly lessen the 
# computational cost of a calculation. It does not usually incurr
# a significant loss in accuracy.

memory 8 gb

molecule dimer {
  Ar 
 @Rn 1 RMb
  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 
  @Rn 1 RMb
  @Ar 1 R
}

# Rvals=[3.0, 3.25, 3.5, 3.65, 3.75, 3.85, 4.0, 4.25, 4.5, 4.75, 5.0, 5.5, 6.0]
Rvals=[3.0, 3.5, 3.75 ]
molecules=[dimer,mono]

set {
  Print 3
}

# We have more basis sets defined here as we have mid-bonds and
# are using density-fitting. 
# Alter the basis sets consistently. So if you change the first
# basis to aug-cc-pVTZ, then also change the df_basis_set to 
# aug-cc-pVTZ-jkfit 

basis {
   assign aug-cc-pVDZ
   # Do not alter this one:
   assign Rn  mb-set
}

df_basis_scf {
   assign aug-cc-pVDZ-jkfit
   # Do not alter this one:
   assign Rn mb-set-ri
}

#df_basis_cc {
#   assign aug-cc-pVDZ-RI
#   # Do not alter this one:
#   assign Rn mb-set-ri
#}

set scf_type df
set mp2_type df
set cc_type  df
set freeze_core True

# Initialize a blank dictionary of counterpoise corrected energies
# (Need this for the syntax below to work)
e_d_hf  = {}
e_d_mp2 = {}
e_d_cc  = {}
e_m_hf  = {}
e_m_mp2 = {}
e_m_cc  = {}

for R in Rvals:
   dimer.R = R
   dimer.RMb = R/2.0  # position of the midbond
   energy('mp2',molecule=dimer)
   # 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
   mono.RMb = R/2.0  # position of the midbond
   energy('mp2',molecule=mono)
   # 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')

   clean()

# Psi4 conversion factors
# psi_hartree2kJmol
# psi_hartree2wavenumbers

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)  [cm^{-1}] \n")
psi4.print_out("---------------------------------------------------------------\n")
for R in Rvals:
   e_hf  = (e_d_hf[R]  -2.0*e_m_hf[R] ) * psi_hartree2wavenumbers
   e_mp2 = (e_d_mp2[R] -2.0*e_m_mp2[R]) * psi_hartree2wavenumbers
   # e_cc  = (e_d_cc[R]  -2.0*e_m_cc[R] ) * psi_hartree2wavenumbers
   e_cc  = 0.0
   psi4.print_out("        %6.3f     %12.6f     %12.6f     %12.6f  \n" % (R, e_hf,e_mp2,e_cc))

Final results of the practical

At the end of this lab session you should have consolidated your data into a single table saved as a text file. The table may take the form:

Ar2_MP2.txt

# Ar2 MP2 interaction energies
# Units: R : Angstrom
#        E : cm^{-1}
#       
#                                                   --------Using mid-bonds---             ----No mid-bonds ---
#   R     aDZ-noCP  aTZ-noCP   DZ-noCP    TZ-noCP   aDZ-CP-Mb aTZ-CP-Mb DZ-CP-Mb TZ-CP-Mb  aDZ-CP   aTZ-CP ... 
  3.000    849.203  683.451    ...        ...       693.597    646.098  ...      ...     ...
  3.250    205.203  106.749                         114.204     94.645  
  3.500    -14.944  -75.276                         -66.813    -76.199 
  3.650    -63.415 -105.688                         -97.067   -104.068 
  ...
  ...

This table will contain at least 12 data sets:

• no-CP: 4 basis sets
• CP: 4 basis sets
• CP with mid-bonds : 4 data sets

Use column headers that are short and descriptive. Only the MP2 data are important as HF yields no binding.

You should visualise the data by plotting it. I suggest you choose $x$- and $y$-axis limits to emphasise the bowl of the potential well. I used all $R$ values, but restricted the energy range to $-150$ to $100$ cm$^{-1}$. I suggest you use Gnuplot to graph the data. This code is available for all platforms so you will be able to install it on your laptop. It is not present on our workstations.