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

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 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:

PBE

  xc xpbe96 cpbe96

PBE0

  xc pbe0

LC-PBE

    xc xcampbe96 1.0 cpbe96 1.0 HFexch 1.0
    cam 0.30 cam_alpha 0.0 cam_beta 1.0

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.

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!

Ar$_2$ job automation

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