attachment:batch_Ar2_PBE0.py of AJMPublic/teaching/electronic-structure/argon2-pes-dft

   1 #!/usr/bin/env python2.7
   2  
   3 import argparse, numpy, subprocess, re, sys
   4 import os
   5  
   6 # Parse command-line arguments
   7  
   8 parser = argparse.ArgumentParser(description = """
   9 Create a set of NWChem input files for the Ar2 dimer, run NWChem on them, and
  10 parse the output for interaction energy.
  11 """)
  12 parser.add_argument("-m", "--min", dest="minR", type=float, default=2.5,
  13                     help="minimum Ar separation in A")
  14 parser.add_argument("-M", "--max", dest="maxR", type=float, default=5.5,
  15                     help="maximum Ar separation in A")
  16 parser.add_argument("-d", "--delta", dest="deltaR", type=float, default=0.5,
  17                     help="step size for Ar separation in A")
  18 parser.add_argument("-b", "--basis", dest="basis", type=str, default="aug-cc-pvtz",
  19                     help="basis set used in calculation")
  20 parser.add_argument("-o", "--output", dest="output", type=str, default="output.txt",
  21                     help="output text file")
  22  
  23 args = parser.parse_args()
  24  
  25 nwchem_dir = os.environ["NWCHEM_INSTALL"]
  26 username = os.environ["USER"]
  27 scratch = os.environ["SCRATCH"]
  28 nproc = "1"
  29  
  30  
  31 # This command uses the minR, maxR and deltaR arguments to create an array.
  32 positions = numpy.arange(args.minR, args.maxR + 0.00001, args.deltaR)
  33 # If you need a non-uniform set of points, you can define them 
  34 # using the numpy.array command as follows. If you wish to use this approach,
  35 # un-comment the following line and comment out the one above.
  36 #positions = numpy.array([ 3.0, 3.5, 3.75, 4.0, 4.25, 4.5, 4.75, 5.0, 5.5, 6.0 ])
  37 print "Basis set ",args.basis
  38 print "positions = ", positions
  39 print "Output to : ",args.output
  40 output = open(args.output,'w')
  41  
  42  
  43 # Conversion factor:
  44 au2cm = 219474.631371
  45  
  46 # Construct input files
  47  
  48 input_file_template = """
  49 title "Ar dimer BSSE corrected DFT interaction energy"
  50  
  51 scratch_dir {scratch}
  52  
  53 geometry units angstrom "Ar+Ar" noautoz
  54  Ar1 0 0 0
  55  Ar2 0 0 {R}
  56 end
  57  
  58 geometry units angstrom "Ar+ghost" noautoz
  59  Ar1  0 0 0
  60  BqAr 0 0 {R}
  61 end
  62  
  63 basis
  64  Ar1   library    {basis}
  65  Ar2   library    {basis}
  66  BqAr  library Ar {basis}
  67 end
  68  
  69 dft
  70   xc pbe0
  71   direct
  72   disp
  73   iterations 60
  74 end
  75  
  76 set geometry "Ar+Ar"
  77 task dft
  78 unset geometry "Ar+Ar"
  79  
  80 dft; vectors hcore; end
  81  
  82 set geometry "Ar+ghost"
  83 task dft
  84 unset geometry "Ar+ghost"
  85 """
  86  
  87 # Regular expressions to find the right lines in the file
  88  
  89 dft_re = re.compile("Total DFT energy =\s+(-?[0-9]+.[0-9]+)")
  90  
  91  
  92 # Dictionary to store the energies we find
  93  
  94 energies = {}
  95  
  96 # Main loop through different position values
  97  
  98 for ar_position in positions:
  99     print "Current position = ",ar_position
 100     input_file_name = "in-R{0:5.3f}.nw".format(ar_position)
 101     output_file_name = "out-R{0:5.3f}.out".format(ar_position)
 102  
 103     # Write input file
 104     with open(input_file_name,"w") as INP:
 105         INP.write(input_file_template.format(R=ar_position, basis=args.basis, 
 106         user=username, scratch=scratch))
 107  
 108     # Run NWChem
 109     output_file = open(output_file_name, "w")
 110     subprocess.call(["mpirun","-np",nproc,"nwchem", input_file_name],
 111                                                 stdout=output_file, stderr=subprocess.STDOUT)
 112     output_file.close()
 113  
 114     # Read output file back in
 115  
 116     output_file = open(output_file_name, "r")
 117     dft_energies = []
 118     for line in output_file:
 119         dft_match = dft_re.search(line)
 120         if dft_match:
 121             dft_energies.append(float(dft_match.group(1)))
 122     output_file.close()
 123  
 124     if not (len(dft_energies) == 2):
 125         print "Apparent error in output file! {0:d} DFT energies found.".format(len(dft_energies))
 126         sys.exit(1)
 127  
 128     print "dft_energies ", dft_energies     
 129     energies[ar_position] = dft_energies
 130     dft_dimer   = energies[ar_position][0]
 131     dft_monomer = energies[ar_position][1]
 132     dft_interaction = dft_dimer - 2*dft_monomer
 133     print "Eint[DFT] = ",dft_interaction*au2cm
 134  
 135 # prepare plots
 136  
 137 dft_dimer = numpy.array([energies[p][0] for p in positions])
 138 dft_monomer = numpy.array([energies[p][1] for p in positions])
 139  
 140 dft_interaction = dft_dimer - 2*dft_monomer
 141  
 142 # convert to cm-1
 143 dft_interaction = au2cm*dft_interaction
 144  
 145 output.write(('%10s '*4 + '\n') % 
 146              ("   R    ","DFT mono","DFT dimr","DFT int "))
 147  
 148 for i in range(len(positions)):
 149     output.write(('%10.6f '*4 + '\n') %
 150     (positions[i],
 151     dft_monomer[i], dft_dimer[i], dft_interaction[i]))
 152  
 153 output.close()