## page was renamed from AJMGrpOnly/teaching/electronic-structure/practical3 ## page was renamed from ajm/teaching/electronic-structure/practical3 <> = 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. {{{#!wiki important Important The STO-3G basis is quite small, so let's now converge the geometry with basis by repeating the above calculations with the * cc-pVDZ * cc-pVTZ * aug-cc-pVDZ * aug-cc-pVTZ basis sets. Use separate files for each basis as you will need to keep track of your results. Make a table of your results. Do you see any trend? Does the geometry converge to the reference results? Is there a basis set where you see an '''apparent''' convergence? }}} === 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''. {{{#!wiki important Important Repeat the above exercise with MP2. Include these new geometries in your table. Do you see a better agreement with the reference results? }}} {{{#!wiki important Important By the end of this exercise you should have a table of the form: || Basis set || RHF |||| MP2 |||| || || $R$ (O-H) || $\theta$ (H-O-H) || $R$ (O-H) || $\theta$ (H-O-H) || || STO-3G || || || || || || cc-pVDZ || || || || || || cc-pVTZ || || || || || || aug-cc-pVDZ || || || || || || aug-cc-pVTZ || || || || || || Reference values || || || || || What are the trends you have observed? How well do they correlate with the material covered in the lecture? You may use the data in the paper by Mas and Szalewicz (http://aip.scitation.org/doi/pdf/10.1063/1.471469) for your reference values. They report the experimental geometry of water on page 4. These data have been taken from Ref 30 of that paper. }}} == 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//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: {{{#!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 || || || || || || || cc-pVQZ || || || || || || || aug-cc-pVDZ || || || || || || || aug-cc-pVTZ || || || || || || || aug-cc-pVQZ || || || || || || }}} == 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. {{{#!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: {{{ #!/usr/bin/env python2.7 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''. {{{#!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? }}}