QM/MM documentation in user guide
[namd.git] / lib / qmmm / run_terachem.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 numpy as np
12 from scipy import spatial
13
14 # TeraChem does not print the point-charge self energy, so we need to
15 # calculate it ourselves. This may take a lot of time if many point charges
16 # are included due to O(N^2) scaling!
17 def computePointChargeSelfEnergy(pclist):
18       energy = 0.0
19       length = len(pclist)
20       pc = np.array(pclist,dtype=np.float64)
21       a = np.arange(0,length,1)
22       coords = pc[:,:3]
23       r = spatial.distance.cdist(coords,coords)/0.52917721067
24       idc = np.triu_indices(length,1)
25       energy = np.sum(pc[:,3][idc[0]]*pc[:,3][idc[1]]/(r[idc]))
26       return energy
27
28
29
30 ################## Hardcoded parameters
31
32 # Here we set up a template for the header and options of the input file.
33
34 excitedState=0
35
36 tcConfigLines1 = """\
37 basis 6-31G*
38 coordinates qmmm.xyz
39 pointcharges point_charges
40 charge 0
41 spinmult 1
42 poptype mulliken
43 method b3lyp
44 dftd no
45 run gradient
46 end
47 """
48
49 inputFilename = sargv[1]
50
51 directory = dirname(inputFilename)
52
53 # Prepares the name of the configuration file based on full path to data file provided by NAMD
54 tcInFileName = directory + "/"
55 tcInFileName += "qmmm_tc.in"
56
57 coordFile = directory + "/qmmm.xyz"
58
59 # Prepares the name of the point charge file based on full path to data file provided by NAMD
60 #%pointcharges "/dev/shm/NAMD_test/0/qmmm_0.input.pntchrg"
61 pcFileName = directory + "/point_charges"
62
63
64 # Prepares the name of the output file based on full path to data file provided by NAMD.
65 # This is where we will direct all tc output, so that we can grab the atom charge
66 # information to pass back to NAMD.
67 tcOutFileName = tcInFileName
68 tcOutFileName += ".out"
69
70 # Prepares the file name for the file which will be read by NAMD
71 finalResFileName = inputFilename
72 finalResFileName += ".result"
73
74 # print("tcInFileName:",tcInFileName)
75 # print("pcFileName:",pcFileName)
76 # print("tcOutFileName:",tcOutFileName)
77 #print("tcGradFileName:",tcGradFileName)
78 # print("finalResFileName:",finalResFileName)
79
80 ################## Reading and parsing NAMD's data ; Writing tc's Input
81 # Reads NAMD data
82 infile = open(inputFilename,"r")
83
84 line = infile.readline()
85
86 # Gets number of atoms in the Quantum Chemistry region (= QM atoms + Link atoms)
87 numQMatms = int(line.split()[0])
88 # Gets number of point charges
89 numPntChr = int(line.split()[1].replace("\n",""))
90
91 #print("numQMatms:",numQMatms,"; numPntChr",numPntChr)
92
93 # stores all lines written to tc's input file
94 outLinesQM = []
95
96 # stores all lines written to tc's point charge file
97 outLinesPC = []
98
99 # The first line in the point charge file is composed only of the total number
100 # of point charges in the file.
101
102 # Identation
103 ident = "  "
104
105 lineIndx = 1
106 lineIndx = 1
107 pointChargeList = []
108 pointChargeDict = {}
109 charges = []
110 for line in infile:
111
112     posx = line.split()[0]
113     posy = line.split()[1]
114     posz = line.split()[2]
115
116     if lineIndx <= numQMatms:
117
118
119         element = line.split()[3].replace("\n","").capitalize()
120
121         outLinesQM.append(ident + " ".join([element,posx,posy,posz]) + "\n")
122
123     else:
124
125
126         charge = line.split()[3]
127         pos = " ".join(line.split()[0:3])
128         charges.append(charge)
129         if pos in pointChargeDict:
130                 print("double occurence: ", pos, pointChargeDict[pos] , " with new charge ", charge)
131                 pointChargeDict[pos] += float(charge)
132                 print("new merged charge: ", pointChargeDict[pos])
133         else:
134                 pointChargeDict[pos] = float(charge)
135
136
137     lineIndx += 1
138
139
140 cnp = np.array(charges,dtype=float)
141 print("Sum of the charges: ", np.sum(cnp))
142
143 outLinesPC.append(str(len(pointChargeDict)) + "\n\n")
144 for k in pointChargeDict:
145         c = pointChargeDict[k]
146         p = k.split()
147         pointChargeList.append([p[0],p[1],p[2],str(c)])
148         outLinesPC.append(" ".join(['{0:.4f}'.format(c),p[0],p[1],p[2]]) + "\n")
149
150 pcSelfEnergy = computePointChargeSelfEnergy(pointChargeList)
151 print("self-energy: " , str(pcSelfEnergy) , " a.u.")
152
153 infile.close()
154
155 ###
156
157 with open(tcInFileName,"w") as outQMFile:
158
159     outQMFile.write(tcConfigLines1)
160
161
162 with open(pcFileName,"w") as outPCFile:
163
164     for line in outLinesPC:
165         outPCFile.write(line)
166
167 with open(coordFile,"w") as cFile:
168         cFile.write(str(numQMatms))
169         cFile.write("\n\n")
170         for line in outLinesQM:
171                 cFile.write(line)
172
173 ################## Run tc
174
175 # We first move the shell to the target directory where calculations are to be
176 # performed
177 import os
178 current_env = os.environ.copy()
179
180 cmdline = "cd " + directory + "; "
181 # Then we run tc with our output file receiving all standard output.
182
183 # here, scipy is not regularly in the python path :(
184 cmdline += "$TeraChem/bin/terachem "
185 cmdline += tcInFileName + " > " + tcOutFileName
186
187 #print("command:",cmdline)
188
189 proc = sp.Popen(args=cmdline, shell=True,env=current_env)
190 proc.wait()
191
192
193 print("terachem finished running.")
194 ################## Process tc results
195
196 scfEnergy = 0
197 iterate = True
198
199 gradientSection = False
200 excSection = False
201 grads = []
202 a0 = 0.52917721067
203 conversionHB = 627.509469/a0
204 eEnergy = 0
205 resultFile = open("scr/results.dat","r")
206 mainFile = open(tcOutFileName,"r")
207
208 if not excitedState:
209         while iterate:
210                 line = mainFile.readline()
211
212                 if line.find("dE/dX") != -1:
213                         gradientSection = True
214                         line = mainFile.readline()
215
216                 if gradientSection == True:
217                         g = line.strip().split()[0:3]
218                         grads.append([-1.0*float(x)*conversionHB for x in g])
219                         if len(grads) == numQMatms:
220                                 gradientSection = False
221                                 break
222
223
224 iterate = True
225 while iterate:
226         line = resultFile.readline()
227         if line == None:
228                 break
229         if line.find("Ground state energy (a.u.):") != -1:
230                 line = resultFile.readline()
231                 scfEnergy = (float(line.strip()) - pcSelfEnergy)*627.509469
232                 if excitedState == 0:
233                         break
234
235         if line.find("dE/dX") != -1:
236                 gradientSection = True
237                 line = resultFile.readline()
238
239         if gradientSection == True:
240                 g = line.strip().split()[1:4]
241                 grads.append([-1.0*float(x)*conversionHB for x in g])
242                 if len(grads) == numQMatms:
243                         gradientSection = False
244                         if excitedState == 0:
245                                 break
246
247         if line.find("Root      Energy") != -1:
248                 line = resultFile.readline()
249                 excSection = True
250
251         if excSection == True:
252                 state, energy = line.strip().split()
253                 if int(state) == excitedState:
254                         # eV to kcal/mol
255                         eEnergy = float(energy) * 23.06054819
256                         break
257
258
259 resultFile.close()
260 finalEnergy = scfEnergy + eEnergy
261 print(scfEnergy,eEnergy,finalEnergy)
262 ###
263
264 # Gets atom charges from the charge_mull.xls output file.
265 mullikenFile = open("scr/charge_mull.xls","r")
266
267 qmCharges = []
268
269 iterate = True
270
271 while iterate:
272     line = mullikenFile.readline()
273     qmCharges.append(line.split()[2].replace("\n","").strip())
274
275     if len(qmCharges) == numQMatms:
276         break
277
278 mullikenFile.close()
279 # Writes the final output file that will be read by NAMD.
280 finFile = open(finalResFileName,"w")
281
282 # The first line constains the enegy
283 finFile.write(str(finalEnergy) + "\n")
284
285 # And all follwoing lines contain gradients and partial charges
286 for i in range(numQMatms):
287
288     finFile.write(" ".join(str(x) for x in grads[i]) + " " + qmCharges[i] + "\n")
289
290 finFile.close()
291
292
293 ##########
294
295 exit(0)