attachment:batch_Ar2_HF_MP2.py of AJMPublic/teaching/electronic-structure/argon2-pes-hf-mp2

   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=1.,
  13                     help="minimum Ar separation in A")
  14 parser.add_argument("-M", "--max", dest="maxR", type=float, default=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 = "2"
  29  
  30 positions = numpy.arange(args.minR, args.maxR + 0.00001, args.deltaR)
  31 print "Basis set ",args.basis
  32 print "positions = ", positions
  33 print "Output to : ",args.output
  34 output = open(args.output,'w')
  35  
  36  
  37 # Conversion factor:
  38 au2cm = 219474.631371
  39  
  40 # Construct input files
  41  
  42 input_file_template = """
  43 title "Ar dimer BSSE corrected MP2 interaction energy"
  44  
  45 scratch_dir {scratch}
  46  
  47 geometry units angstrom "Ar+Ar" noautoz noautosym
  48  Ar1 0 0 0
  49  Ar2 0 0 {R}
  50 end
  51  
  52 geometry units angstrom "Ar+ghost" noautoz noautosym
  53  Ar1 0 0 0
  54  Bq2 0 0 {R}
  55 end
  56  
  57 basis
  58  Ar1 library    {basis}
  59  Ar2 library    {basis}
  60  Bq2 library Ar {basis}
  61 end
  62  
  63 set geometry "Ar+Ar"
  64 task mp2
  65 unset geometry "Ar+Ar"
  66  
  67 scf; vectors atomic; end
  68  
  69 set geometry "Ar+ghost"
  70 task mp2
  71 unset geometry "Ar+ghost"
  72 """
  73  
  74 # Regular expressions to find the right lines in the file
  75  
  76 scf_re = re.compile("SCF energy\s+(-?[0-9]+.[0-9]+)")
  77 mp2_re = re.compile("Total MP2 energy\s+(-?[0-9]+.[0-9]+)")
  78  
  79 # Dictionary to store the energies we find
  80  
  81 energies = {}
  82  
  83 # Main loop through different position values
  84  
  85 for ar_position in positions:
  86     print "Current position = ",ar_position
  87     input_file_name = "in-R{0:5.3f}.nw".format(ar_position)
  88     output_file_name = "out-R{0:5.3f}.out".format(ar_position)
  89  
  90     # Write input file
  91     with open(input_file_name,"w") as INP:
  92         INP.write(input_file_template.format(R=ar_position, basis=args.basis,
  93                     user=username, scratch=scratch))
  94  
  95     # Run NWChem
  96     output_file = open(output_file_name, "w")
  97     subprocess.call(["mpirun","-np",nproc,"nwchem", input_file_name],
  98                                                 stdout=output_file, stderr=subprocess.STDOUT)
  99     output_file.close()
 100  
 101     # Read output file back in
 102  
 103     output_file = open(output_file_name, "r")
 104     scf_energies = []
 105     mp2_energies = []
 106     for line in output_file:
 107         scf_match = scf_re.search(line)
 108         if scf_match:
 109             scf_energies.append(float(scf_match.group(1)))
 110         else:
 111             mp2_match = mp2_re.search(line)
 112             if mp2_match:
 113                 mp2_energies.append(float(mp2_match.group(1)))
 114     output_file.close()
 115  
 116     if not (len(scf_energies) == 2 and len(mp2_energies) == 2):
 117         print "Apparent error in output file! {:d} SCF energies and {:d} MP2 energies found.".format(
 118             len(scf_energies), len(mp2_energies))
 119         sys.exit(1)
 120  
 121     energies[ar_position] = scf_energies + mp2_energies
 122     scf_dimer   = energies[ar_position][0]
 123     scf_monomer = energies[ar_position][1]
 124     mp2_dimer   = energies[ar_position][2]
 125     mp2_monomer = energies[ar_position][3]
 126     scf_interaction = scf_dimer - 2*scf_monomer
 127     mp2_interaction = mp2_dimer - 2*mp2_monomer
 128     print "Eint[SCF] = ",scf_interaction*au2cm
 129     print "Eint[MP2] = ",mp2_interaction*au2cm
 130  
 131 # prepare plots
 132  
 133 scf_dimer = numpy.array([energies[p][0] for p in positions])
 134 scf_monomer = numpy.array([energies[p][1] for p in positions])
 135 mp2_dimer = numpy.array([energies[p][2] for p in positions])
 136 mp2_monomer = numpy.array([energies[p][3] for p in positions])
 137  
 138 scf_interaction = scf_dimer - 2*scf_monomer
 139 mp2_interaction = mp2_dimer - 2*mp2_monomer
 140  
 141 # convert to cm-1
 142 scf_interaction = au2cm*scf_interaction
 143 mp2_interaction = au2cm*mp2_interaction
 144  
 145 output.write(('%10s '*7 + '\n') % 
 146              ("   R    ","SCF mono","SCF dimr",
 147               "SCF int ","MP2 mono","MP2 dimr","MP2 int "))
 148  
 149 for i in range(len(positions)):
 150     output.write(('%10.6f '*7 + '\n') %
 151     (positions[i],
 152     scf_monomer[i], scf_dimer[i], scf_interaction[i],
 153     mp2_monomer[i], mp2_dimer[i], mp2_interaction[i]))
 154  
 155 output.close()