QM/MM documentation in user guide
[namd.git] / lib / qmmm / run_gaussian.py
1 #!/usr/bin/python
2
3 # by Maximilian Scheurer (mscheurer@ks.uiuc.edu), Dec 2016
4
5 ##################  Imports
6
7 from sys import argv as sargv
8 from sys import exit
9 from os.path import dirname
10 import subprocess as sp
11 import re
12 import numpy as np
13
14
15 ################## Command file parameters
16
17 gaussianConfigLines1 = """\
18 %NProcShare=4
19 # B3LYP/6-31G* Charge Force
20
21 qmmm
22
23 0 1
24 """
25
26 # if you specified an excited state ( TD(NStates=2,Root=1) ) in the com file and want to give its gradient to NAMD,
27 # specify the number of the excited state here,
28 # so that the python script can read the corresponding excited state energy
29 # if no TD calculation is requested, set excitedState to 0
30 excitedState=0
31
32
33 ################## Processing of file names
34 gaussianWhitespace = "\n\n"
35
36 inputFilename = sargv[1]
37
38 #print(inputFilename)
39
40 directory = dirname(inputFilename)
41
42 # Prepares the name of the configuration file based on full path to data file provided by NAMD
43 gaussianInFileName = directory + "/"
44 gaussianInFileName += "qmmm.com"
45
46 # Name of the Gaussian log file
47 gaussianOutFileName = gaussianInFileName +".log"
48
49 # Prepares the file name for the file which will be read by NAMD
50 finalResFileName = inputFilename
51 finalResFileName += ".result"
52
53 ################## Reading and parsing NAMD's data ; Writing gaussian's Input
54
55 # Reads NAMD data
56 infile = open(inputFilename,"r")
57
58 line = infile.readline()
59
60 # Gets number of atoms in the Quantum Chemistry region (= QM atoms + Link atoms)
61 numQMatms = int(line.split()[0])
62 # Gets number of point charges
63 numPntChr = int(line.split()[1].replace("\n",""))
64
65 # print("numQMatms:",numQMatms,"; numPntChr",numPntChr)
66
67 # stores all lines written to gaussian's input file
68 outLinesQM = []
69
70 # Identation
71 ident = "  "
72
73 lineIndx = 1
74 pointChargeList = []
75 pointChargeDict = {}
76 charges = []
77 for line in infile:
78
79     posx = line.split()[0]
80     posy = line.split()[1]
81     posz = line.split()[2]
82
83     if lineIndx <= numQMatms:
84
85         # gaussian's format requires the fileds to be ordered begining with the
86         # atom's element symbol, and followed by the XYZ coordinates.
87
88         element = line.split()[3].replace("\n","")
89
90         outLinesQM.append(" ".join([element,posx,posy,posz]) + "\n")
91
92     else:
93         #output linebreak to separate atoms from charges
94         if lineIndx == numQMatms+1:
95             outLinesQM.append("\n")
96
97         # gaussian's format requires the fields to be ordered begining with the
98         # XYZ, and followed by the charge .
99         pos = " ".join(line.split()[0:3])
100         charge = line.split()[3]
101         charges.append(charge)
102         if pos in pointChargeDict:
103                 print "double occurence: ", pos, pointChargeDict[pos] , " with new charge ", charge
104                 pointChargeDict[pos] += float(charge)
105                 print "new merged charge: ", pointChargeDict[pos]
106         else:
107                 pointChargeDict[pos] = float(charge)
108
109
110     lineIndx += 1
111
112 cnp = np.array(charges,dtype=float)
113 print "Sum of the charges: ", np.sum(cnp)
114
115 for k in pointChargeDict:
116         c = pointChargeDict[k]
117         p = k.split()
118         pointChargeList.append([p[0],p[1],p[2],str(c)])
119         outLinesQM.append(" ".join([p[0],p[1],p[2],'{0:.16f}'.format(c)]) + "\n")
120
121 # print len(pointChargeList)
122 infile.close()
123
124 ###
125
126 with open(gaussianInFileName,"w") as outQMFile:
127
128     outQMFile.write(gaussianConfigLines1)
129
130     for line in outLinesQM:
131         outQMFile.write(line)
132
133     outQMFile.write(gaussianWhitespace)
134
135 ################## Run gaussian
136
137 # We first move the shell to the target directory where calculations are to be
138 # performed
139 cmdline = "cd " + directory + "; "
140 # Then we run gaussian with our output file receiving all standard output.
141
142 # we probably need to set some environment variables:
143 import subprocess, os
144 current_env = os.environ.copy()
145 current_env["GAUSS_EXEDIR"] = "/path/to/gaussian/exedir"
146 # subprocess.Popen(my_command, env=my_env)
147
148 cmdline += "/path/to/gaussian/exedir/g09 "
149 cmdline += gaussianInFileName + " " + gaussianOutFileName
150
151 # print "command:", cmdline
152 proc = sp.Popen(args=cmdline, shell=True, env=current_env)
153 proc.wait()
154
155 ########### READING Gaussian output
156 excitedStateEnergy=0
157
158 tmpOutFile = open(gaussianOutFileName,"r")
159
160 qmCharges = []
161 gradient = []
162
163 # Bohr radius for conversion
164 a0 = 0.52917721067
165 conversionHB = 627.509469/a0
166 selfEnergyGaussian = 0.0
167
168 # Iterates ultil we find the section of output that contains atomic partial
169 # charges for QM atoms
170 chargeSection = False
171
172 gradientSection = False
173
174 scfenergyFound = False
175
176 excitedFound = False
177
178 iterate = True
179
180 selfEnergyFound = False
181
182 while iterate:
183
184     line = tmpOutFile.readline()
185
186     if line.find("Mulliken charges:") != -1:
187         chargeSection = True
188         # Skips a dividing line
189         line = tmpOutFile.readline()
190         continue
191
192     if not scfenergyFound:
193         if line.find("SCF Done") != -1:
194             reg = re.search("SCF Done:(.+)=(.*)([-+]?\d*\.\d+|\d+)(.+)A.U.",line)
195             scfenergy = 627.509469 * float(reg.group(2).strip())
196             print "SCF energy: ", scfenergy
197             scfenergyFound = True
198
199     if not excitedFound and excitedState != 0:
200         if line.find("Excited State") != -1:
201             line = line.strip().replace(":","").split()
202             if int(line[2]) != excitedState:
203                 continue
204             line =tmpOutFile.readline()
205             # check if we really requested the excited state that the gradient will be calculated for!
206             while line.find("This state for optimization and/or second-order correction.") == -1:
207                 line = tmpOutFile.readline()
208             line = tmpOutFile.readline()
209             line = line.strip()
210             reg = re.search("([-+]?\d*\.\d+|\d+)(.+)",line)
211             excitedStateEnergy = 627.509469 * float(reg.group(1))
212             print "Excited State energy: " , excitedStateEnergy
213             excitedFound = True
214
215     if line.find("Center     Atomic                   Forces (Hartrees/Bohr)") != -1:
216         gradientSection = True
217         #skip 2 lines
218         line = tmpOutFile.readline()
219         line = tmpOutFile.readline()
220         line = tmpOutFile.readline()
221
222     if gradientSection:
223         line = line.strip().split()
224         # don't switch the sign of the number, as Gaussian prints F = -dE/dx, so we can directly pass the numbers to NAMD
225         gradient.append( [ str( float(line[2])*conversionHB ), str( float(line[3])*conversionHB ), str( float(line[4])*conversionHB ) ])
226         if len(gradient) == numQMatms:
227             gradientSection = False
228             iterate = False
229
230     if chargeSection:
231         qmCharges.append(line.split()[2].replace("\n","").strip())
232
233     if len(qmCharges) == numQMatms:
234         chargeSection = False
235
236     if selfEnergyFound != True and line.find("Self energy of the charges") != -1:
237         line = line.strip().split()
238         selfEnergyGaussian = float(line[-2])
239         print "Gaussian self energy of the charges: ", selfEnergyGaussian
240         selfEnergyFound = True
241
242 tmpOutFile.close()
243
244 finFile = open(finalResFileName,"w")
245
246 if excitedState == 0:
247         finFile.write(str( scfenergy - selfEnergyGaussian*627.509469 ) + "\n")
248 else:
249         finFile.write(str( excitedStateEnergy - selfEnergyGaussian*627.509469 ) + "\n")
250
251 for i in range(numQMatms):
252
253     finFile.write(" ".join(gradient[i]) + " " + qmCharges[i] + "\n")
254
255 finFile.close()
256
257
258 ##########
259
260 exit(0)