<> = Ar$_2$ : DFT calculations of the PES = In this lab we will re-calculate the Ar dimer PES using various density functionals. The idea here is to get a feel for what we may expect of density functional theory. == Structure of the DFT block == The NWCHem input file will be similar to what you have used for HF, but now you include the DFT commands as follows [[attachment:NWChem DFT block]] {{{ DFT ... END Task DFT }}} === Options in the DFT block === The DFT calculations are conceptually similar to the HF ones you've already done, but they typically involve a lot more options. The main [[http://www.nwchem-sw.org/index.php/Release65:Density_Functional_Theory_for_Molecules|NWChem page on DFT is here]], and here are the issues you will need to keep in mind: * '''Functionals''': There are many density functional included in NWChem. See the DFT page for a list. The ones we will use are: PBE, PBE0, LC-PBE and LC-PBE0. * '''Asymptotic correction''': If we are interested in properties that depend on the density tails (high ranking multipole moments, polarizabilities, excitation energies), then we will need some kind of asymptotic correction. This is not needed for the LC functionals as these already should have the correct long-range behaviour. There are two corrections possible in NWChem: '''LB94''' and '''CS00'''. The CS00 correction is more suitable. You normally need to specify a shift, but, if not supplied, it will be estimated. * '''SIC''': Sometimes you need to include a self-interaction correction (SIC) as DFT is not free from this source of error. * '''Dispersion correction''': One of the big failings of commonly used density functionals is the lack of van der Waals interactions. NWChem provides the '''DISP''' and '''XDM''' models to (partially) correct this. The '''DISP''' model simply includes a simple dispersion correction based on pre-computed $C_6$ dispersion coefficients and fitted splicing models. The '''XDM''' model is more sophisticated. A sample DFT block might be: [[attachment:DFT commands with LC-PBE and dispersion correction]] {{{ DFT xc xcampbe96 1.0 cpbe96 1.0 HFexch 1.0 cam 0.30 cam_alpha 0.0 cam_beta 1.0 Direct Iterations 60 DISP END }}} This set of commands defines the '''LC-PBE''' functional with the Grimme dispersion correction. We have set a maximum number of iterations to 60 and requested that NWChem performs the calculation in memory (the '''Direct''' option). You find the commands needed to define functionals on the DFT webpage of the NWChem site. But here are the commands for the four functionals listed above: [[attachment:PBE]] {{{ xc xpbe96 cpbe96 }}} [[attachment:PBE0]] {{{ xc pbe0 }}} [[attachment:LC-PBE]] {{{ xc xcampbe96 1.0 cpbe96 1.0 HFexch 1.0 cam 0.30 cam_alpha 0.0 cam_beta 1.0 }}} [[attachment:LC-PBE0]] {{{ xc xcampbe96 1.0 cpbe96 1.0 HFexch 1.0 cam 0.30 cam_alpha 0.25 cam_beta 0.75 }}} The LC-type functionals are the long-range corrected functionals that use HF-type exchange for the long-range part and DFT-type exchange for the short-range. The ''cam'' scheme is a splicing scheme that allows us to mesh the two parts together. == Interaction energies using DFT == You will encounter some of the biggest problems with DFT when using it to compute interaction energies. These are * Lack of van der Waals: Density functionals are generally local or semi-local. These functionals are fundamentally unable to describe the long-range van der Waals interaction. The ''DISP'' term includes an empirical dispersion correction to alleviate the problem, but it is approximate and doesn't always work. * Self-interaction: Unlike HF, most formulations of DFT suffer from a self-interaction error (they are not even correct for the hydrogen atom!). There are schemes for (partially) correcting this. For interaction energies it is important we use an asymptotic correction to correct for the 1-electron SI error and get more well-behaved density tails. The ''CS00'' correction does this. Another way of enforcing an SI correction is to use an LC-type of functional. === Ar$_2$ === === Tasks === Here you will be investigating the performance of various density functionals on the argon dimer system. Frankly, this is not really representative and the water dimer, with its hydrogen bond, would have been a better example. But this example is simple enough as we can vary the dimer geometry quite easily. Here is what you will be investigating: * '''Convergence with basis''': First of all we need to decide on which basis is appropriate. DFT should behave like HF, but we should check. So, perform the Argon dimer interaction energy calculations with, say, PBE0, using the six basis sets you used in the last exercise. You will need to modify the Python code used to include commands for the DFT calculation. Do not use any dispersion, asymptotic, or self-interaction corrections. * '''Which basis set do you think is suitable?''' * '''Dispersion correction''': We need to decide on which dispersion correction to use. The choices are ''DISP'' and ''XDM''. Repeat the previous calculation with your optimum choice of basis and each of these two dispersion corrections. Compare your results with the reference MP2 aQZ/CP results from last week. Evaluate the dispersion corrections. Can you decide on a winner? * '''Which functional''': There are four you should try: PBE, PBE0, LC-PBE and LC-PBE0. Use your optimal basis and dispersion correction with each of these. Compare against the reference MP2 potential and decide on a suitable functional. {{{#!wiki important Important How will you go about investigating all these aspects of DFT? * You could either compute interaction energies at all $R$ values you used to make the argon dimer interaction energy plot using MP2. But this could take time. * Or you could fix the dimer geometry at the minimum energy separation, and decide on the best basis, functional and dispersion correction, and subsequently compute the interaction energies at all separations. The latter is definitely less work, but is not a good idea as it can be misleading to look at the performance of a method for a single geometry. Instead I suggest you do the following: 1. Fix the argon dimer at its minimum energy geometry. 1. Use the PBE functional to compute the interaction energy using a variety of basis sets. Decide on the best basis set. Keep this fixed for the next steps. 1. Using the basis from above, calculate the interaction energy using the PBE, PBE0, LC-PBE, and LC-PBE0 functionals (you may include more) at '''all''' dimer geometries. There is a script below which will allow you to do this easily. 1. Then repeat this using the two dispersion corrections. Make plots. Don't put too many curves on a graph as it will not be easy to see them. Analyse your data. }}} {{{#!wiki important Important You should end this exercise with a reasonably good idea of the best basis set, XC functional and dispersion correction for use with Ar$_2$. Unfortunately, these conclusions will not necessarily extend to other systems, but we will tackle that problem some other time. }}} === Inputs and codes === Here is a sample template for the Python code. This is the template for the ''PBE0'' functional with the ''DISP'' dispersion correction. You will need to insert this into the Python program. Make sure you make the required changes to the template for each of the tasks you need to perform. [[attachment:Sample template for the Python code]] {{{ input_file_template = """ title "Ar dimer BSSE corrected DFT 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 BqAr 0 0 {R} end basis Ar1 library {basis} Ar2 library {basis} BqAr library Ar {basis} end dft xc pbe0 direct disp iterations 60 end set geometry "Ar+Ar" task dft unset geometry "Ar+Ar" dft; vectors atomic; end set geometry "Ar+ghost" task dft unset geometry "Ar+ghost" """ }}} == Python code to do the batch jobs == This one is from Mark. Thanks, Mark! [[attachment:Ar$_2$ job automation]] {{{ #!/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=2.5, help="minimum Ar separation in A") parser.add_argument("-M", "--max", dest="maxR", type=float, default=5.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"] 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) # This command uses the minR, maxR and deltaR arguments to create an array. positions = numpy.arange(args.minR, args.maxR + 0.00001, args.deltaR) # If you need a non-uniform set of points, you can define them # using the numpy.array command as follows. If you wish to use this approach, # un-comment the following line and comment out the one above. #positions = numpy.array([ 3.0, 3.5, 3.75, 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 DFT interaction energy" scratch_dir /scratch/HDD_2T/esm/{user}/nwchem geometry units angstrom "Ar+Ar" Ar1 0 0 0 Ar2 0 0 {R} end geometry units angstrom "Ar+ghost" Ar1 0 0 0 BqAr 0 0 {R} end basis Ar1 library {basis} Ar2 library {basis} BqAr library Ar {basis} end dft xc pbe0 direct disp iterations 60 end set geometry "Ar+Ar" task dft unset geometry "Ar+Ar" dft; vectors hcore; end set geometry "Ar+ghost" task dft unset geometry "Ar+ghost" """ # Regular expressions to find the right lines in the file dft_re = re.compile("Total DFT 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)) # 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") dft_energies = [] for line in output_file: dft_match = dft_re.search(line) if dft_match: dft_energies.append(float(dft_match.group(1))) output_file.close() if not (len(dft_energies) == 2): print "Apparent error in output file! {0:d} DFT energies found.".format(len(dft_energies)) sys.exit(1) print "dft_energies ", dft_energies energies[ar_position] = dft_energies dft_dimer = energies[ar_position][0] dft_monomer = energies[ar_position][1] dft_interaction = dft_dimer - 2*dft_monomer print "Eint[DFT] = ",dft_interaction*au2cm # prepare plots dft_dimer = numpy.array([energies[p][0] for p in positions]) dft_monomer = numpy.array([energies[p][1] for p in positions]) dft_interaction = dft_dimer - 2*dft_monomer # convert to cm-1 dft_interaction = au2cm*dft_interaction output.write(('%10s '*4 + '\n') % (" R ","DFT mono","DFT dimr","DFT int ")) for i in range(len(positions)): output.write(('%10.6f '*4 + '\n') % (positions[i], dft_monomer[i], dft_dimer[i], dft_interaction[i])) output.close() }}}