QM/MM documentation in user guide
[namd.git] / lib / qmmm / runORCA.py
1 #!/usr/bin/python3
2
3 # by Marcelo Melo (melomcr@gmail.com)
4
5 # If you would like to run a Quantum Chemistry software which is not supported by NAMD,
6 # an interface script can be created in order to convert NAMD's input into the format
7 # expected by the QC software, and subsequently convert the software's output inro the
8 # format expected by NAMD.
9
10 # This script creates an interface for ORCA, so that it can be ran using the custom 
11 # QM software option provided by NAMD. It is an example script, since ORCA is supported
12 # natively by NAMD's QM/MM interface, and is intended to be used by those wishing to 
13 # test theyr own software in hybrid MD simulations through NAMD.
14
15 #######  NAMD Files Format
16 # INPUT: This input file will contain on the first line the number of QM atoms (X)
17 # and the number of point charges in the file (Y, which may be 0 or more), separated 
18 # by a space. The following X+Y lines will have four (4) fields: X, Y and Z coordinates,
19 # and a fourth field which will depend on the type of entry. For QM atoms, the 
20 # field will contain the element of the QM atom. For point charge lines, the field
21 # will contain the charge of the point charge.
22
23 # OUTPUT: The expected output file whould be placed in the same directory as the 
24 # input file, and should be named "*inputfile*.result" (meaning it will have the
25 # same path and name of the input file, plus the suffix '.result'). This file
26 # should have, on its first line, the energy of the system, and on the following 
27 # X lines (where X is the number of QM atoms passed in the input file), four (4) 
28 # fields: the x, y and z components of the TOTAL FORCE applied on that atom, and 
29 # on the fourth field, the charge of the atom. If the user indicates that charges
30 # from the QM software should not be used (see "QMChargeMode"), the fourth field
31 # should have zeroes, but should not be empty.
32
33
34 #######  ORCA Files Format
35
36 # ORCA input file is composed of a variable-sized header, where all info pertaining
37 # to the simulation is specified, and a corrdinates section, where all atom positions
38 # are specified, along with some informations regarding the overall system (namely
39 # the charge and multiplicity of the system).
40 #
41 # We will get the header and the initial part of the coordinates section from 
42 # lines specified in this script. The data would be gathered and writen by NAMD,
43 # but in this case we need to previously prepare it from hardcoded values or by 
44 # reading from another file.
45
46 # ORCA expects a separate file with point charge information. This info is provided
47 # by NAMD in the same data file, so we will extract it and paste it in a separate 
48 # file in the same folder. The name and path to the point charge file needs to be 
49 # included in the configuration portion of ORCA's input file.
50
51 ##################  Imports
52
53 from sys import argv as sargv
54 from sys import exit 
55 from os.path import dirname
56 import subprocess as sp
57
58 ################## Hardcoded parameters
59
60 # Here we set up a template for the header and options of the input file.
61
62 # The choice of method and basis-set, as well as parallelization and charge
63 # output are all made here.
64 # For ORCA, the keyword "ENGRAD" is essential for the proper output of gradients
65 # over QM atoms.
66
67 orcaConfigLines1 = """\
68 !  B3LYP 6-31G* Grid4 PAL4 EnGrad TightSCF
69 %output 
70   PrintLevel Mini 
71   Print [ P_Mulliken ] 1
72   Print [ P_AtCharges_M ] 1
73 end
74
75 """
76
77 # We set up the template for the point charge file indicator
78
79 orcaConfigLines2 = "%pointcharges \""
80
81 # And set up the header for the coordinates section of the input file.
82
83 orcaConfigLines3 = """\
84 %coords
85   CTyp xyz
86   Charge 0.000000
87   Mult 1.000000
88   Units Angs
89   coords
90
91 """
92
93 ################## Processing of file names
94
95 inputFilename = sargv[1]
96
97 #print(inputFilename)
98
99 directory = dirname(inputFilename)
100
101 # Prepares the name of the configuration file based on full path to data file provided by NAMD
102 orcaInFileName = directory + "/"
103 orcaInFileName += "qmmm.input"
104
105 # Prepares the name of the point charge file based on full path to data file provided by NAMD
106 #%pointcharges "/dev/shm/NAMD_test/0/qmmm_0.input.pntchrg"
107 pcFileName = orcaInFileName
108 pcFileName += ".pntchrg"
109
110 orcaConfigLines2 += pcFileName + "\"\n"
111
112
113 # Prepares the name of the output file based on full path to data file provided by NAMD.
114 # This is where we will direct all ORCA output, so that we can grab the atom charge
115 # information to pass back to NAMD.
116 orcaOutFileName = orcaInFileName
117 orcaOutFileName += ".TmpOut"
118
119 # Prepares the name of the gradient file based on full path to data file provided by NAMD.
120 # This is where ORCA will write the gradient information on QM atoms.
121 orcaGradFileName = orcaInFileName
122 orcaGradFileName += ".engrad"
123
124 # Prepares the file name for the file which will be read by NAMD
125 finalResFileName = inputFilename
126 finalResFileName += ".result"
127
128 #print("orcaInFileName:",orcaInFileName)
129 #print("pcFileName:",pcFileName)
130 #print("orcaOutFileName:",orcaOutFileName)
131 #print("orcaGradFileName:",orcaGradFileName)
132 #print("finalResFileName:",finalResFileName)
133
134 ################## Reading and parsing NAMD's data ; Writing ORCA's Input
135
136 # Reads NAMD data
137 infile = open(inputFilename,"r") 
138
139 line = infile.readline()
140
141 # Gets number of atoms in the Quantum Chemistry region (= QM atoms + Link atoms)
142 numQMatms = int(line.split()[0])
143 # Gets number of point charges
144 numPntChr = int(line.split()[1].replace("\n",""))
145
146 #print("numQMatms:",numQMatms,"; numPntChr",numPntChr)
147
148 # stores all lines written to ORCA's input file
149 outLinesQM = []
150
151 # stores all lines written to ORCA's point charge file
152 outLinesPC = []
153
154 # The first line in the point charge file is composed only of the total number
155 # of point charges in the file.
156 outLinesPC.append(str(numPntChr) + "\n")
157
158 # Identation
159 ident = "  "
160
161 lineIndx = 1
162 for line in infile:
163     
164     posx = line.split()[0]
165     posy = line.split()[1]
166     posz = line.split()[2]
167     
168     if lineIndx <= numQMatms:
169         
170         # ORCA's format requires the fileds to be ordered begining with the
171         # atom's element symbol, and followed by the XYZ coordinates.
172                 
173         element = line.split()[3].replace("\n","")
174         
175         outLinesQM.append(ident + " ".join([element,posx,posy,posz]) + "\n")
176         
177     else:
178         
179         # ORCA's format requires the fileds to be ordered begining with the
180         # charge, and followed by the XYZ coordinates.
181         
182         charge = line.split()[3]
183         
184         outLinesPC.append(" ".join([charge,posx,posy,posz]) + "\n")
185     
186     lineIndx += 1
187     
188
189 # Finalizes the formating for ORCA's input file, one "end" to terminate the
190 # atomic coordinates, and nother to terminate the section.
191 outLinesQM.append(ident + "end" + "\n")
192 outLinesQM.append("end" + "\n")
193
194
195 infile.close()
196
197 ###
198
199 with open(orcaInFileName,"w") as outQMFile:
200     
201     outQMFile.write(orcaConfigLines1)
202     outQMFile.write(orcaConfigLines2)
203     outQMFile.write(orcaConfigLines3)
204     
205     for line in outLinesQM:
206         outQMFile.write(line)
207
208 with open(pcFileName,"w") as outPCFile:
209     
210     for line in outLinesPC:
211         outPCFile.write(line)
212
213
214
215 ################## Run ORCA
216
217 # We first move the shell to the target directory where calculations are to be
218 # performed
219 cmdline = "cd " + directory + "; "
220 # Then we run orca with our output file receiving all standard output.
221 cmdline += "/data/Programas/orca_3_0_3_linux_x86-64/orca "
222 cmdline += orcaInFileName + " > " + orcaOutFileName
223
224 #print("command:",cmdline)
225
226 proc = sp.Popen(args=cmdline, shell=True)
227 proc.wait()
228
229
230 # Runs a secondary process
231 cmdline2 = "/home/melomcr/Research/NAMD_QMMM/Scripts/saveDipole.sh "
232 cmdline2 += orcaInFileName
233
234 proc2 = sp.Popen(args=cmdline2, shell=True)
235
236 ################## Process ORCA results
237
238 # Gets system energy and gradients from ORCA's output file.
239
240 gradFile = open(orcaGradFileName,"r")
241     
242 # Skips to the line with number of atoms
243 for i in range(3):
244     gradFile.readline()
245
246 orcaNumQMAtms = int(gradFile.readline().replace("\n",""))
247
248 # Runs basic sanity check
249 if orcaNumQMAtms != numQMatms:
250     print("ERROR: Expected",numQMatms,"but found",orcaNumQMAtms,"atoms in engrad file!")
251     exit(1)
252
253 # Skips to the line with final system energy
254 for i in range(3):
255     gradFile.readline()
256
257 finalEnergy = gradFile.readline().replace("\n","").strip()
258
259 print("ORCA energy: ", finalEnergy,"Eh")
260
261 # All energies are given in Eh (Hartree)
262 # NAMD needs energies in kcal/mol
263 # The conversion factor is 627.509469
264 finalEnergy = str( float(finalEnergy) * 627.509469 )
265
266 print("ORCA energy: ", finalEnergy,"kcal/mol")
267
268 # Skips to the lines with gradients
269 for i in range(3):
270     gradFile.readline()
271
272 # Stores all gradients (yes, it would be faster with numpy, but this is just an example)
273 grads = []
274 for i in range(orcaNumQMAtms):
275     
276     grads.append( list() )
277     
278     # All *GRADIENTS* are given in Eh/a0 (Hartree over Bohr radius)
279     # NAMD needs *FORCES* in kcal/mol/angstrons
280     # The conversion factor is -1*627.509469/0.529177 = -1185.82151
281     for j in range(3):
282         gradComp = gradFile.readline().replace("\n","").strip()
283         gradComp = float(gradComp) * -1185.82151
284         grads[i].append( str(gradComp) )
285
286 gradFile.close()
287
288
289 ###
290
291 # Gets atom charges from the temporary output file.
292 tmpOutFile = open(orcaOutFileName,"r")
293
294 qmCharges = []
295
296 # Iterates ultil we find the section of output that contains atomic partial 
297 # charges for QM atoms
298 chargeSection = False
299
300 iterate = True
301
302 while iterate:
303     
304     line = tmpOutFile.readline()
305     
306     if line.find("MULLIKEN ATOMIC CHARGES") != -1:
307         chargeSection = True
308         
309         # Skips a dividing line
310         line = tmpOutFile.readline()
311         
312         continue
313     
314     if chargeSection:
315         qmCharges.append(line.split()[3].replace("\n","").strip())
316         pass
317     
318     if len(qmCharges) == numQMatms:
319         break
320
321 tmpOutFile.close()
322
323 # Writes the final output file that will be read by NAMD.
324
325 finFile = open(finalResFileName,"w")
326
327 # The first line constains the enegy
328 finFile.write(finalEnergy + "\n")
329
330 # And all follwoing lines contain gradients and partial charges
331 for i in range(numQMatms):
332     
333     finFile.write(" ".join(grads[i]) + " " + qmCharges[i] + "\n")
334
335 finFile.close()
336
337
338 ##########
339
340 # Makes sure the secondary process has ended befire we return to NAMD.
341 proc2.wait()
342
343
344 exit(0)