## page was renamed from AJMGrpOnly/teaching/electronic-structure/argon2-pes-hf-mp2 ## page was renamed from ajm/teaching/electronic-structure/argon2-pes-hf-mp2 <> = 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 [[ajm/teaching/intermolecular-interactions|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: {{{#!wiki warning Warning First of all, we are now performing significant calculations so we cannot use the home directory to run our jobs. Instead we use a scratch directory. This is done by defining ''scratch_dir'' in the input file. Please change this directory path to include your own user name!!! I have asked our administrators to create scratch directories for each of you. to check to see if you have one type: {{{ $ ls /scratch/HDD_2T/esm/ }}} 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. {{{#!wiki tip Tip For more on the '''set''' and '''unset''' command see the [[http://www.nwchem-sw.org/index.php/Top-level#SET|NWChem Wiki SET/UNSET]] }}} 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: [[attachment:ar2_accurate_potential_patkowskimfs05_molp103.pdf|Patkowski, Murdachaew, Fou, and Szalewicz, Mol. Phys., 103, 10-20 (2005)]] Have a look at this paper as it contains a lot more data.) {{{#!wiki tip Tip Conversion factor: 1 Hartree = 219474.631 cm-1 }}} Points to note: * No binding with HF. Why? * MP2 appears to give a *larger* binding energy than the reference. {{{#!wiki important Important Repeat the calculation using the aug-cc-pVDZ basis. This time you should get * $E_{\rm int}({\rm HF}) = 82.575$ cm$^{-1}$, and * $E_{\rm int}({\rm MP2}) = -103.24$ cm$^{-1}$ Amazing! MP2 is now in better agreement with the reference energy. But how can this be when we have reduced the basis size? }}} 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. {{{#!wiki important Important Plot these energies (both the HF and MP2) as a function of basis set. Use the basis size on the $x$-axis - you will find the size printed in the NWChem output. Q: Is there a pattern that will allow you to make a reasonable guess at the complete basis limit? Put your data in a table of the form. Your interaction energies (in cm$^{-1}$) go in the columns labeled 'NoCP' = no counterpoise correction. We will fill in the ''CP'' columns next. || Basis set || Basis size || RHF || MP2 || RHF || MP2 || || || || NoCP || NoCP || CP || CP || || cc-pVDZ || || || || || || || cc-pVTZ || || || || || || || aug-cc-pVDZ || || || || || || || aug-cc-pVTZ || || || || || || }}} == 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. {{{#!wiki important Important Include these new energies on your plot of interaction energy versus basis size. Q: Can you see a pattern this time? Is it clearer or less clear than last time? Put your data in the table we began in the previous step. }}} = 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: [[attachment:batch_Ar2_HF_MP2.py]] {{{ #!/usr/bin/env python2.7 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''. {{{#!wiki important Important Tasks: * Calculate the interaction energies using the aug-cc-pVDZ, aug-cc-pVTZ, and aug-cc-pVQZ basis sets. * Plot these using Gnuplot. Make sure you include the data points and use ''smooth csplines'' to draw a smooth line between the data points. * What pattern do you see? * Include data from the paper linked above. These data will be your reference values. * How close are your results to these reference values? }}} {{{#!wiki important Important Sometimes it is good to do a calculation another way to see what happens: * Modify the python code to calculate the interaction energies without the CP correction. * Calculate the interaction energies without this correction and plot them as above. * What do you conclude? }}}