Contents
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
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 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:
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:
xc xpbe96 cpbe96
xc pbe0
xc xcampbe96 1.0 cpbe96 1.0 HFexch 1.0
cam 0.30 cam_alpha 0.0 cam_beta 1.0xc 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.
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:
- Fix the argon dimer at its minimum energy geometry.
- 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.
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.
- 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.
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.
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!
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()