# 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
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.0
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$

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:

1. Fix the argon dimer at its minimum energy geometry.
2. 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.
3. 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.

4. 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.

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"
unset geometry "Ar+Ar"

dft; vectors atomic; end

set geometry "Ar+ghost"
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.
""")
help="minimum Ar separation in A")
help="maximum Ar separation in A")
help="step size for Ar separation in A")
help="basis set used in calculation")
help="output text file")

args = parser.parse_args()

nwchem_dir = os.environ["NWCHEM_TOP"]
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"
unset geometry "Ar+Ar"

dft; vectors hcore; end

set geometry "Ar+ghost"
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:

# 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()

AJMPublic/teaching/electronic-structure/practical-ar2-dft (last edited 2021-04-14 13:53:56 by apw109)