Optimizations and the PES

Here we will learn how to use NWChem for geometry optimizations and to generate the PES of a dimer.

Geometry Optimization

This is a common task: you have a system and want to obtain an optimized structure. How do you do it? Here's an example:

Memory 500 mb

charge 0
Geometry units bohr
  O     0.0000000000    0.0000000000    0.0000000000
  H1   -1.4536519600    0.0000000000   -1.1216873200
  H2    1.4536519600    0.0000000000   -1.1216873200
End

Basis "ao basis" spherical
  H  library  STO-3G
  O  library  STO-3G
End

Title "H2O STO-3G Geometry Optimization"

scf
  Singlet
end
task scf optimize

Put these commands in a file titled h2o-sto3g-scf-opt.nw. Run it the usual way:

$ mpirun.mpich -np 1 nwchem h2o-sto3g-scf-opt.nw >& h2o-sto3g-scf-opt.out

Towards the bottom of the output you will find the lines:

                Final and change from initial internal coordinates
                --------------------------------------------------



                                Z-matrix (autoz)
                                --------

 Units are Angstrom for bonds and degrees for angles

      Type          Name      I     J     K     L     M      Value       Change
      ----------- --------  ----- ----- ----- ----- ----- ---------- ----------
    1 Stretch                  1     2                       0.98941    0.01778
    2 Stretch                  1     3                       0.98941    0.01778
    3 Bend                     2     1     3               100.02684   -4.66316

 ======
                                internuclear distances
 ------------------------------------------------------------------------------
       center one      ||      center two      || atomic units ||       a.u.
 ------------------------------------------------------------------------------
    2 H1               ||   1 O                ||     1.86971  ||     1.86971
    3 H2               ||   1 O                ||     1.86971  ||     1.86971
 ------------------------------------------------------------------------------
                         number of included internuclear distances:          2
 ======



 ======
                                 internuclear angles
 ------------------------------------------------------------------------------
        center 1       ||       center 2       ||       center 3       ||  degrees
 ------------------------------------------------------------------------------
    2 H1               ||   1 O                ||   3 H2               ||   100.03
 ------------------------------------------------------------------------------
                            number of included internuclear angles:          1
 ======

This bit contains the final (optimised) bond lengths and angles. The reference O-H bond length is 0.095639 nm (Halkier et al, Chem. Phys. Lett. 274, 235-241 (1997)) and H-O-H bond angle (from the same reference) is 104.54 degrees. These results here are quite different from these.

Optimisation at the MP2 level

It will not, and the reason is that Hartree-Fock is after all, approximate. We will use a higher level of theory to repeat the above task, this time using second-order perturbation theory: MP2. We will learn about this next week, but for now all you need to know is that MP2 is the next level up from HF.

Here's the input file for an MP2 optimization:

Memory 500 mb

charge 0
Geometry units bohr
  O     0.0000000000    0.0000000000    0.0000000000
  H1   -1.4536519600    0.0000000000   -1.1216873200
  H2    1.4536519600    0.0000000000   -1.1216873200
End

Basis "ao basis" spherical
  H  library  STO-3G
  O  library  STO-3G
End

Title "H2O STO-3G MP2 Geometry Optimization"

scf
  Singlet
end
task mp2 optimize

Save this as h2o-sto3g-mp2-opt.nw.

Interaction Energies

The interaction energy is the energy of stabilization when two of more species come together and form some sort of bond.

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 /scratch/HDD_2T/esm/<YOUR USER NAME>/nwchem

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:

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 /scratch/HDD_2T/esm/{YOUR USER NAME}/nwchem

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 Eint 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:

 
import argparse, numpy, subprocess, re, sys
import os
#import matplotlib.pyplot as plt
 
# 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_TOP"]
username = os.environ["USER"]
scratch = os.environ["SCRATCH"]
nproc = "1"

#usage = 'Usage: %s minR maxR dR output' % sys.argv[0]
#try:
#        minR   = float(sys.argv[1])
#        maxR   = float(sys.argv[2])
#        deltaR = float(sys.argv[3])
#        file   = sys.argv[4]
#except:
#        print usage; sys.exit(1)
 
 
positions = numpy.arange(args.minR, args.maxR + 0.00001, args.deltaR)
#positions = numpy.arange(minR, maxR + 0.00001, deltaR)
#positions = numpy.array([ 3.0, 3.3, 3.4, 3.5, 3.6, 3.7, 3.75, 3.8, 3.9, 4.0, 4.25, 4.5, 4.75, 5.0, 5.5, 6.0 ])
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}/nwchem
 
geometry units angstrom "Ar+Ar"
 Ar1 0 0 0
 Ar2 0 0 {R}
end
 
geometry units angstrom "Ar+ghost"
 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.mpich","-np 1",os.path.join(nwchem_dir,"bin/nwchem"), input_file_name],
    subprocess.call(["mpirun.mpich","-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()
 
 
#args.output.write(("{:>10s},  "*6 + "\n").format(
#    "SCF mono", "SCF di", "SCF int",
#    "MP2 mono", "MP2 di", "MP2 int"))
#
#for i in range(len(positions)):
#    args.output.write(("{:10.6f},  "*6 + "\n").format(
#        scf_monomer[i], scf_dimer[i], scf_interaction[i],
#        mp2_monomer[i], mp2_dimer[i], mp2_interaction[i]))
#
#args.output.close()
 
## Plot:
#
#fig = plt.figure(figsize=(9,9))
#ax = fig.add_subplot(111)
#
## edit this to change which lines plotted
#
#lines_to_plot = [scf_interaction, mp2_interaction]
#
#for line in lines_to_plot:
#    ax.plot(positions, line)
#
#plt.show()

This is a rather long code, but it contains a number of python2.7 lines that have been commented out and replaced by python2.6-compatible lines. You need to make a change to the input_file_template and change the scratch_dir to be the one relevant to you. Otherwise no further change should be needed.

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

  $ chmod +x nwchem_batch_comanche.py

To see the usage use:

  $ ./nwchem_batch_comanche.py --help

you should get:

usage: nwchem_batch_comanche.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,

  $ ./nwchem_batch_comanche.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.