AMH: DFT calculations with NWChem

We now re-visit the AMH dimers (both of them) using density functional theory.

We already know from the Ar$_2$ examples we have done that most density functionals are rather poor for van der Waals dimers. The AMH system is not quite a van der Waals system, but it has a strong dispersion contribution, so any density functional we use will need to be corrected for the missing dispersion energy.

Further, according to Griffiths et al., there is a problem with PBE: it overbinds the //P4/nmm// phase due to a charge-transfer problem that manifests in an overbinding at long range.

You will need to test PBE for both of these issues and propose fixes.

You will need to decide on

• Which functional to use: PBE or PBE0 or some other? Remember you need to justify your choice.
• Dispersion correction or no dispersion correction?
• Explain the observations highlighted on the previous lab page using these calculations.
• The goal of this report is to use your calculations to best explain the main results of the paper - particularly those I have highlighted in the previous laboratory page.

Scripts

The script on the AMH laboratory page can be modified to do a DFT calculation. I've included an example here.

Here's a script for the molecular dimer. This one does a PBE+D calculation. It's provided only as a template. You will need to decide just which calculation you need.

AMH PBE+D

import argparse, numpy, subprocess, re, sys, os
import numpy, re, sys, subprocess
#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("-o", "--output", dest="output", type=str, default="output.txt",
                    help="output text file")
parser.add_argument("-np", "--nproc", dest="nproc", type=int, default=1,
                    help="Number of processors to be used.")
parser.add_argument("-b", "--basis", dest="basis", type=str, default="aug-cc-pvtz",
                    help="basis set used in calculation")

args = parser.parse_args()

nwchem_dir = os.environ["NWCHEM_TOP"]
username = os.environ["USER"]
nproc = args.nproc


# I have commented out the uniform grid:
#positions = numpy.arange(args.minR, args.maxR + 0.00001, args.deltaR)
# Instead, use a non-uniform grid specific to your system:
positions = numpy.array([1.7,2.1,2.3,2.5,2.7,3.1,4.0,5.0,10.0,20.0])

print "Basis set ",args.basis
print "positions = ", positions
print "Output to : ",args.output
print "Number of processors : ",nproc
output = open(args.output,'w')


# Conversion factor:
au2cm = 219474.631371
au2kJ = 2625.499

# Construct input files

input_file_template = """
title "NH3..H2O BSSE corrected DFT interaction energy"

scratch_dir /scratch/HDD_2T/{user}/nwchem

geometry units angstrom "NH3+H2O"
  N       0.00000     0.00000     0.00000
  H1      0.98601     0.00000    -0.30465
  H2     -0.46058     0.84108    -0.38142
  H3     -0.46058    -0.84108    -0.38142
  O       0.00000     0.00000     {RzO}
  H4      0.00000     0.00000     {RzH4}
  H5     -0.96568    -0.00000     {RzH5}
end

geometry units angstrom "NH3+ghost"
  N       0.00000     0.00000     0.00000
  H1      0.98601     0.00000    -0.30465
  H2     -0.46058     0.84108    -0.38142
  H3     -0.46058    -0.84108    -0.38142
  BqO     0.00000     0.00000     {RzO}
  BqH     0.00000     0.00000     {RzH4}
  BqH    -0.96568    -0.00000     {RzH5}
end

geometry units angstrom "ghost+H2O"
  BqN       0.00000     0.00000     0.00000
  BqH       0.98601     0.00000    -0.30465
  BqH      -0.46058     0.84108    -0.38142
  BqH      -0.46058    -0.84108    -0.38142
  O         0.00000     0.00000     {RzO}
  H4        0.00000     0.00000     {RzH4}
  H5       -0.96568    -0.00000     {RzH5}
end

basis
 N   library    {basis}
 H1  library    {basis}
 H2  library    {basis}
 H3  library    {basis}
 O   library    {basis}
 H4  library    {basis}
 H5  library    {basis}
 BqN   library  N  {basis}
 BqH   library  H  {basis}
 BqO   library  O  {basis}
end

dft 
  xc xpbe96 cpbe96
  iterations 200
  disp
  direct
end

set geometry "NH3+H2O"
charge 0
task dft
unset geometry "NH3+H2O"

dft; vectors hcore; end

set geometry "NH3+ghost"
charge 0
task dft
unset geometry "NH3+ghost"

dft; vectors hcore; end

set geometry "ghost+H2O"
charge 0
task dft
unset geometry "ghost+H2O"

"""

# z coordinates for O, H4 and H5 : H2O will be translated along the z-axis
# Units: Angstrom
zO  =  0.0
zH4 = -1.03900 
zH5 =  0.27836

# 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 delta_z in positions:
    print "Current position = ",delta_z
    input_file_name = "in-R{0:5.3f}.nw".format(delta_z)
    output_file_name = "out-R{0:5.3f}.out".format(delta_z)

    # Write input file
    Rz_O  = zO + delta_z  
    Rz_H4 = zH4 + delta_z
    Rz_H5 = zH5 + delta_z
    RzO  = """{0:10.6f}""".format(Rz_O)
    RzH4 = """{0:10.6f}""".format(Rz_H4)
    RzH5 = """{0:10.6f}""".format(Rz_H5)
    input_file = open(input_file_name, "w")
    input_file.write(input_file_template.format(RzO=RzO,RzH4=RzH4,RzH5=RzH5,user=username,basis=args.basis))
    input_file.close()

    # Run NWChem

    output_file = open(output_file_name, "w")
    subprocess.call(["mpirun.mpich","-np",str(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) == 3):
        print "Apparent error in output file! {0:d} DFT energies found.".format(len(dft_energies))
        sys.exit(1)

    energies[delta_z] = dft_energies
    dft_dimer    = energies[delta_z][0]
    dft_monomerA = energies[delta_z][1]
    dft_monomerB = energies[delta_z][2]
    dft_interaction = dft_dimer - dft_monomerA - dft_monomerB
    print "Eint[DFT] = ",dft_interaction*au2kJ

# prepare plots

dft_dimer = numpy.array([energies[p][0] for p in positions])
dft_monomerA = numpy.array([energies[p][1] for p in positions])
dft_monomerB = numpy.array([energies[p][2] for p in positions])

dft_interaction = dft_dimer - dft_monomerA - dft_monomerB

# convert to kJ/mol
dft_interaction = au2kJ*dft_interaction

output.write(('%10s '*5 + '\n') % 
             ("   R    ","DFT monoA","DFT monoB","DFT dimr","DFT int "))

for i in range(len(positions)):
    output.write(('%10.6f '*5 + '\n') %
    (positions[i],
    dft_monomerA[i], dft_monomerB[i], dft_dimer[i], dft_interaction[i]))

output.close()