QM/MM documentation in user guide
[namd.git] / lib / qmmm / run_qchem.py
1 #!/usr/bin/python
2
3 # by Maximilian Scheurer (mscheurer@ks.uiuc.edu), June 2017
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 import math
14
15
16 ################## Command file parameters
17 qchemNumberOfCores = 2
18
19 qchemConfigLine1 = """\
20 $rem
21 JOBTYPE force
22 METHOD b3lyp
23 BASIS 6-31G*
24 SCF_CONVERGENCE 7
25 $end
26
27 $molecule
28 0 1
29 """
30
31
32
33
34
35 # if you specified an excited state ( CIS, ADC or whatever ) in the input file and want to give its gradient to NAMD,
36 # specify the number of the excited state here,
37 # so that the python script can read the corresponding excited state energy
38 # if no exc. state calculation is requested, set excitedState to 0
39 excitedState=0
40
41 # add the following lines to qchemConfigLine1 to calculate e.g. 4 excited states with TDDFT,
42 # where the gradient will be calculated on the 1st excited state
43 # CIS_N_ROOTS 4
44 # CIS_STATE_DERIV 1
45 # CIS_SINGLETS true
46 # CIS_TRIPLETS false
47
48
49
50 ################## Processing of file names
51
52 inputFilename = sargv[1]
53
54 #print(inputFilename)
55
56 directory = dirname(inputFilename)
57
58 # Prepares the name of the configuration file based on full path to data file provided by NAMD
59 qchemInFileName = directory + "/"
60 qchemInFileName += "qmmm.in"
61
62 # Name of the qchem log file
63 qchemOutFileName = qchemInFileName +".out"
64
65 # Prepares the file name for the file which will be read by NAMD
66 finalResFileName = inputFilename
67 finalResFileName += ".result"
68
69 ################## Reading and parsing NAMD's data ; Writing gaussian's Input
70
71 # Reads NAMD data
72 infile = open(inputFilename,"r")
73
74 line = infile.readline()
75
76 # Gets number of atoms in the Quantum Chemistry region (= QM atoms + Link atoms)
77 numQMatms = int(line.split()[0])
78 # Gets number of point charges
79 numPntChr = int(line.split()[1].replace("\n",""))
80
81 # print("numQMatms:",numQMatms,"; numPntChr",numPntChr)
82
83 # stores all lines written to gaussian's input file
84 outLinesQM = []
85
86 # Identation
87 ident = "  "
88
89 lineIndx = 1
90 pointChargeList = []
91 pointChargeDict = {}
92 charges = []
93 for line in infile:
94
95     posx = line.split()[0]
96     posy = line.split()[1]
97     posz = line.split()[2]
98
99     if lineIndx <= numQMatms:
100
101         # qchem's format requires the fileds to be ordered begining with the
102         # atom's element symbol, and followed by the XYZ coordinates.
103
104         element = line.split()[3].replace("\n","")
105
106         outLinesQM.append(" ".join([element,posx,posy,posz]) + "\n")
107
108     else:
109         #output linebreak to separate atoms from charges
110         if lineIndx == numQMatms+1:
111             outLinesQM.append("\n")
112
113         # qchems's format requires the fields to be ordered begining with the
114         # XYZ, and followed by the charge .
115         pos = " ".join(line.split()[0:3])
116         charge = line.split()[3]
117         charges.append(charge)
118         if pos in pointChargeDict:
119                 print "double occurence: ", pos, pointChargeDict[pos] , " with new charge ", charge
120                 pointChargeDict[pos] += float(charge)
121                 print "new merged charge: ", pointChargeDict[pos]
122         else:
123                 pointChargeDict[pos] = float(charge)
124
125
126     lineIndx += 1
127
128 cnp = np.array(charges,dtype=float)
129 print "Sum of the charges: ", np.sum(cnp)
130
131 outLinesQM.append("$end \n \n $external_charges \n")
132 for k in pointChargeDict:
133         c = pointChargeDict[k]
134         p = k.split()
135         pointChargeList.append([p[0],p[1],p[2],str(c)])
136         outLinesQM.append(" ".join([p[0],p[1],p[2],'{0:.16f}'.format(c)]) + "\n")
137
138 outLinesQM.append("$end")
139 # print len(pointChargeList)
140 infile.close()
141
142 ###
143
144 with open(qchemInFileName,"w") as outQMFile:
145
146     outQMFile.write(qchemConfigLine1)
147
148     for line in outLinesQM:
149         outQMFile.write(line)
150
151     # outQMFile.write(gaussianWhitespace)
152
153 ################## Run qchem
154
155 # We first move the shell to the target directory where calculations are to be
156 # performed
157 cmdline = "cd " + directory + "; "
158 # Then we run qchem with our output file receiving all standard output.
159
160 # we probably need to set some environment variables:
161 import subprocess, os
162 current_env = os.environ.copy()
163 current_env["QCSCRATCH"] = directory
164 # subprocess.Popen(my_command, env=my_env)
165
166 #load qchem from modules
167 cmdline += "qchem -nt %d " % qchemNumberOfCores
168 cmdline += qchemInFileName + " " + qchemOutFileName
169
170 print "command:", cmdline
171 proc = sp.Popen(args=cmdline, shell=True, env=current_env)
172 proc.wait()
173
174 ########### READING QCHEM output
175 excitedStateEnergy=0
176
177 tmpOutFile = open(qchemOutFileName,"r")
178
179 qmCharges = []
180 gradient = []
181
182 # Bohr radius for conversion
183 a0 = 0.52917721067
184 conversionHB = 627.509469/a0
185 selfEnergy = 0.0
186
187 # Iterates ultil we find the section of output that contains atomic partial
188 # charges for QM atoms
189 chargeSection = False
190
191 gradientSection = False
192
193 scfenergyFound = False
194
195 excitedFound = False
196
197 iterate = True
198
199 selfEnergyFound = False
200
201 while iterate:
202
203     line = tmpOutFile.readline()
204
205     if line.find("Ground-State Mulliken Net Atomic Charges") != -1:
206         chargeSection = True
207         # Skips 3 dividing lines
208         for x in range(3):
209             line = tmpOutFile.readline()
210         continue
211
212     if not scfenergyFound:
213         if line.find("SCF   energy in the final basis set") != -1:
214             scfenergy = 627.509469 * float(line.strip().split()[-1])
215             print "SCF energy: ", scfenergy
216             scfenergyFound = True
217
218     if not excitedFound and excitedState != 0:
219         if line.find("Excited State") != -1:
220             line = line.strip().replace(":","").split()
221             if int(line[2]) != excitedState:
222                 continue
223             line =tmpOutFile.readline()
224             # check if we really requested the excited state that the gradient will be calculated for!
225             while line.find("This state for optimization and/or second-order correction.") == -1:
226                 line = tmpOutFile.readline()
227             line = tmpOutFile.readline()
228             line = line.strip()
229             reg = re.search("([-+]?\d*\.\d+|\d+)(.+)",line)
230             excitedStateEnergy = 627.509469 * float(reg.group(1))
231             print "Excited State energy: " , excitedStateEnergy
232             excitedFound = True
233
234     if line.find("Gradient of") != -1:
235         gradientSection = True
236         # line = tmpOutFile.readline()
237
238
239     if gradientSection:
240         gradient = np.zeros(shape=(numQMatms,3))
241
242         ln = 0
243         blockColumnNumber = -1
244         atomLineIndex = 0
245         gradientRead = 0
246         while (gradientRead == 0):
247             line = tmpOutFile.readline()
248             line = line.strip().split()
249             if ln % 4 == 0:
250                 blockColumnNumber = len(line)
251                 # print(line, blockColumnNumber)
252                 atomLineIndex += blockColumnNumber
253             else:
254                 xyzIndex = (ln % 4)-1
255                 # print("current coordinate: "+ str(xyzIndex))
256                 currentColumnNumber = blockColumnNumber
257                 columnCounter = 1
258                 for atomIndex in range(atomLineIndex-currentColumnNumber, atomLineIndex):
259                     gradient[atomIndex,xyzIndex] = float(line[columnCounter])
260                     #print(atomIndex,columnCounter)
261                     if atomIndex == (numQMatms-1) and xyzIndex == 2:
262                         gradientRead = 1
263                     columnCounter += 1
264             ln+=1
265
266         iterate = 0
267
268     if chargeSection:
269         qmCharges.append(line.split()[2].replace("\n","").strip())
270
271     if len(qmCharges) == numQMatms:
272         chargeSection = False
273
274     if selfEnergyFound != True and line.find("Charge-charge energy") != -1:
275         line = line.strip().split()
276         selfEnergy = float(line[-2])
277         print "Self energy of the charges: ", selfEnergy
278         selfEnergyFound = True
279
280 tmpOutFile.close()
281
282 finFile = open(finalResFileName,"w")
283
284 if excitedState == 0:
285         finFile.write(str( scfenergy - selfEnergy*627.509469 ) + "\n")
286         print "Corrected SCF energy: ", scfenergy-selfEnergy*627.509469
287 else:
288         finFile.write(str( excitedStateEnergy - selfEnergy*627.509469 ) + "\n")
289
290 forces=np.multiply(-1.0*conversionHB,gradient)
291 forces=forces.tolist()
292 #print forces
293 for i in range(numQMatms):
294     finFile.write(' '.join('{0:.7f}'.format(c) for c in forces[i]) + " " + qmCharges[i] + "\n")
295
296 finFile.close()
297
298
299 ##########
300
301 exit(0)