Interaction energies

The interaction energy is often more important than the energies of the individual components. For example, we don't really care much what the energies of each of the atoms in a water molecule are, as long as they are more strongly bound in the molecule than in isolation. The difference in energies of the bound system from that of the unbound system is the interaction energy. In a chemical bond this energy can be quite large, but in a van der Waals system, it will be small. For more on this topic see my notes on Intermolecular Interaction energies.

Goals of this exercise

We aim to compute the potential energy surface (PES) for the Ar$_2$ system. This is a weakly bound van der Waals system with a rather simple, 1-dimensional PES. That is good as this makes it relatively simple to calculate and subsequently model using an analytical function. Once we have this function, we will, in principle, be able to input it into a simulation code like DL_POLY to perform simulations of the bulk.

But before we can compute this PES, we need to know which methods to use, and how to get accurate interaction energies, so that our PES is accurate. In this exercise, we will aim to understand the Hartree--Fock (HF) and second-order Moller--Plesset (MP2, also called MBPT2) methods. We will also aim to understand something called the basis-set superposition error (BSSE) and how it can be corrected. At the end of this exercise we will have (or should have) reference energies for this dimer.

The Argon dimer

For two Ar atoms separated by distance R this takes the simple form \begin{equation} E_{\rm int}(R) = E[\mbox{Ar}\cdots\mbox{Ar}](R) - 2*E[\mbox{Ar}] \end{equation} where $E[\mbox{Ar}\cdots\mbox{Ar}](R)$ is the energy of the complex of two Ar atoms separated by distance $R$. And $E[\mbox{Ar}]$ is the energy of a single argon atom. So this looks easy to compute.

An example input file to do this calculation is:

title "Ar dimer no BSSE correction : MP2 interaction energy"

scratch_dir /home/{your user name here}/scratch/

geometry units angstrom "Ar+Ar"
  Ar1 0 0 0
  Ar2 0 0 3.76
end

geometry units angstorm "Ar"
  Ar1 0 0 0
end

basis
  Ar1 library aug-cc-pvtz
  Ar2 library aug-cc-pvtz
end

set geometry "Ar+Ar"
task mp2
unset geometry "Ar+Ar"

set geometry "Ar"
task mp2
unset geometry "Ar"

Save this in file Ar2-noCPcorr-atz.nw. The noCPcorr stands for no counterpoise (CP) correction - I will explain what this means soon.

More about the above file:

Is your user name listed there? If so, type:

  $ mkdir /scratch/HDD_2T/esm/${USER}/nwchem

This should create a directory for your NWChem scratch files. If it fails that means we have a problem. Let me know.

From now on we will NOT run any job on /home/${USER}!!! }}}

OK, now to the other parts of this file:

• The file contains two calculations in one.

• We do this by defining two geometries: Ar+Ar and Ar

• The Ar+Ar job has the argon atoms placed 3.76 Angstrom apart. This is close to the minimum energy separation.

• Then use the set command to load the Ar+Ar system into the NWChem database. And calculate the MP2 energy of this system. And then remove it from the database using the unset command.

• Then do the same for the Ar system.

Run the job using

  $ mpirun.mpich -np 1 nwchem Ar2-noCPcorr-atz.nw >& Ar2-noCPcorr-atz.out

Run it to get the energies (in a.u. unless specified):

Basis: aug-cc-pVTZ
Method    E(Ar...Ar)             E(Ar)           Eint (cm-1)
=====
HF       -1053.627091087923  -526.813754874704     +91.88
MP2      -1054.148750447087  -527.074074460254    -132.02
=====

Here I have used the separation of R = 3.76 Angstrom Reference $E_{\rm int}$ at a separation of $R=3.75$ Angstrom is: -99.47 cm-1 (Ref: Patkowski, Murdachaew, Fou, and Szalewicz, Mol. Phys., 103, 10-20 (2005) Have a look at this paper as it contains a lot more data.)

Points to note:

• No binding with HF. Why?
• MP2 appears to give a *larger* binding energy than the reference.

Let's try and understand this. Make a calculation of $E_{\rm int}$ at the HF and MP2 levels of theory using a range of basis sets:

• cc-pVDZ, cc-pVTZ, cc-pVQZ
• and their augmented varieties: aug-cc-pVDZ etc.

Basis Set Superposition Error (BSSE)

There is a problem with these calculations. It's called the Basis set superposition error (BSSE). Here's how it goes:

BSSE: We always use finite basis sets for our calculations of energies. These are incomplete and introduce an inconsistency in the calculations of the energies of the dimer (Ar...Ar) and monomer (Ar). In the former the effective basis space is larger than the latter. If the energy method was variational, this would mean that E(Ar) would not be as low as it should be compared with E(Ar...Ar). The resulting error in Eint is called the basis set superposition error.

Resolution of the problem: Boys and Bernardi proposed using the dimer basis for both calculations. We calculate E(Ar...Ar) as before. But calculate E(Ar) using additionally the basis of the partner. This method of correcting the BSSE is called the counterpoise of CP correction.

In file Ar2/Ar2-CPcorr-atz.nw is an example of how the CP correction can be applied:

title "Ar dimer BSSE corrected MP2 interaction energy" 

scratch_dir /home/{your user name here}/scratch/

geometry units angstrom "Ar+Ar" 
  Ar1 0 0 0 
  Ar2 0 0 3.76
end

geometry units angstrom "Ar+ghost" 
  Ar1 0 0 0 
  Bq2 0 0 3.76
end

basis 
  Ar1 library aug-cc-pvtz 
  Ar2 library aug-cc-pvtz 
  Bq2 library Ar aug-cc-pvtz 
end

set geometry "Ar+Ar"
task mp2
unset geometry "Ar+Ar"

scf; vectors atomic; end

set geometry "Ar+ghost"
task mp2
unset geometry "Ar+ghost"

• Now we have defined a "ghost" atom when calculating the energy of Ar. This ghost atom doesn't contain any electrons. All it contains is the basis set of the missing partner.
• In NWChem these ghost atoms are designated 'Bq'.

• The line scf; vectors atomic; end is needed for the convergence of the Ar+ghost case. I'm not entirely sure why we need it, but we do.

Let's now repeat the above set of calculations for $E_{\rm int}$ with the CP correction. Use all 6 basis sets specified above.

Python code to perform an Ar$_2$ scan

If you wish to calculate the interaction energy for a set of argon separations, the above procedure will be tedious. Instead we can use a script like the one given here:

batch_Ar2_HF_MP2.py

 
import argparse, numpy, subprocess, re, sys
import os
 
# Parse command-line arguments
 
parser = argparse.ArgumentParser(description = """
Create a set of NWChem input files for the Ar2 dimer, run NWChem on them, and
parse the output for interaction energy.
""")
parser.add_argument("-m", "--min", dest="minR", type=float, default=1.,
                    help="minimum Ar separation in A")
parser.add_argument("-M", "--max", dest="maxR", type=float, default=5.,
                    help="maximum Ar separation in A")
parser.add_argument("-d", "--delta", dest="deltaR", type=float, default=0.5,
                    help="step size for Ar separation in A")
parser.add_argument("-b", "--basis", dest="basis", type=str, default="aug-cc-pvtz",
                    help="basis set used in calculation")
parser.add_argument("-o", "--output", dest="output", type=str, default="output.txt",
                    help="output text file")
 
args = parser.parse_args()
 
nwchem_dir = os.environ["NWCHEM_INSTALL"]
username = os.environ["USER"]
scratch = os.environ["SCRATCH"]
nproc = "2"
 
positions = numpy.arange(args.minR, args.maxR + 0.00001, args.deltaR)
print "Basis set ",args.basis
print "positions = ", positions
print "Output to : ",args.output
output = open(args.output,'w')
 
 
# Conversion factor:
au2cm = 219474.631371
 
# Construct input files
 
input_file_template = """
title "Ar dimer BSSE corrected MP2 interaction energy"
 
scratch_dir {scratch}
 
geometry units angstrom "Ar+Ar" noautoz noautosym
 Ar1 0 0 0
 Ar2 0 0 {R}
end
 
geometry units angstrom "Ar+ghost" noautoz noautosym
 Ar1 0 0 0
 Bq2 0 0 {R}
end
 
basis
 Ar1 library    {basis}
 Ar2 library    {basis}
 Bq2 library Ar {basis}
end
 
set geometry "Ar+Ar"
task mp2
unset geometry "Ar+Ar"
 
scf; vectors atomic; end
 
set geometry "Ar+ghost"
task mp2
unset geometry "Ar+ghost"
"""
 
# Regular expressions to find the right lines in the file
 
scf_re = re.compile("SCF energy\s+(-?[0-9]+.[0-9]+)")
mp2_re = re.compile("Total MP2 energy\s+(-?[0-9]+.[0-9]+)")
 
# Dictionary to store the energies we find
 
energies = {}
 
# Main loop through different position values
 
for ar_position in positions:
    print "Current position = ",ar_position
    input_file_name = "in-R{0:5.3f}.nw".format(ar_position)
    output_file_name = "out-R{0:5.3f}.out".format(ar_position)
 
    # Write input file
    with open(input_file_name,"w") as INP:
        INP.write(input_file_template.format(R=ar_position, basis=args.basis,
                    user=username, scratch=scratch))
 
    # Run NWChem
    output_file = open(output_file_name, "w")
    subprocess.call(["mpirun","-np",nproc,"nwchem", input_file_name],
                                                stdout=output_file, stderr=subprocess.STDOUT)
    output_file.close()
 
    # Read output file back in
 
    output_file = open(output_file_name, "r")
    scf_energies = []
    mp2_energies = []
    for line in output_file:
        scf_match = scf_re.search(line)
        if scf_match:
            scf_energies.append(float(scf_match.group(1)))
        else:
            mp2_match = mp2_re.search(line)
            if mp2_match:
                mp2_energies.append(float(mp2_match.group(1)))
    output_file.close()
 
    if not (len(scf_energies) == 2 and len(mp2_energies) == 2):
        print "Apparent error in output file! {:d} SCF energies and {:d} MP2 energies found.".format(
            len(scf_energies), len(mp2_energies))
        sys.exit(1)
 
    energies[ar_position] = scf_energies + mp2_energies
    scf_dimer   = energies[ar_position][0]
    scf_monomer = energies[ar_position][1]
    mp2_dimer   = energies[ar_position][2]
    mp2_monomer = energies[ar_position][3]
    scf_interaction = scf_dimer - 2*scf_monomer
    mp2_interaction = mp2_dimer - 2*mp2_monomer
    print "Eint[SCF] = ",scf_interaction*au2cm
    print "Eint[MP2] = ",mp2_interaction*au2cm
 
# prepare plots
 
scf_dimer = numpy.array([energies[p][0] for p in positions])
scf_monomer = numpy.array([energies[p][1] for p in positions])
mp2_dimer = numpy.array([energies[p][2] for p in positions])
mp2_monomer = numpy.array([energies[p][3] for p in positions])
 
scf_interaction = scf_dimer - 2*scf_monomer
mp2_interaction = mp2_dimer - 2*mp2_monomer
 
# convert to cm-1
scf_interaction = au2cm*scf_interaction
mp2_interaction = au2cm*mp2_interaction
 
output.write(('%10s '*7 + '\n') % 
             ("   R    ","SCF mono","SCF dimr",
              "SCF int ","MP2 mono","MP2 dimr","MP2 int "))
 
for i in range(len(positions)):
    output.write(('%10.6f '*7 + '\n') %
    (positions[i],
    scf_monomer[i], scf_dimer[i], scf_interaction[i],
    mp2_monomer[i], mp2_dimer[i], mp2_interaction[i]))
 
output.close()

Save the above commands in a file called batch_Ar2_HF_MP2.py and make it executable using

  $ chmod +x batch_Ar2_HF_MP2.py

To see the usage use:

  $ ./batch_Ar2_HF_MP2 --help  

you should get:

usage: batch_Ar2_HF_MP2.py [-h] [-m MINR] [-M MAXR] [-d DELTAR] [-b BASIS]
                           [-o OUTPUT]

Create a set of NWChem input files for the Ar2 dimer, run NWChem on them, and
parse the output for interaction energy.

optional arguments:
  -h, --help            show this help message and exit
  -m MINR, --min MINR   minimum Ar separation in A
  -M MAXR, --max MAXR   maximum Ar separation in A
  -d DELTAR, --delta DELTAR
                        step size for Ar separation in A
  -b BASIS, --basis BASIS
                        basis set used in calculation
  -o OUTPUT, --output OUTPUT
                        output text file

Now you should be able to run it using, for example,

  $ ./batch_Ar2_HF_MP2.py -m 3.0 -M 6.0 -d 0.5 -b aug-cc-pVDZ -o Ar2_adz_hf_mp2_scan.txt

This will calculate interaction energies at the HF and MP2 levels of theory for the argon dimer at separations 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.5, 5.0, 5.5 and 6.0 Angstrom and save the results in file Ar2_adz_hf_mp2_scan.txt.