OpenMD 3.1
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
omdCombined
1#!@Python3_EXECUTABLE@
2"""OMDCombined
3
4Opens two omd files, one with a solute structure and one with a
5solvent structure. Viable with solvents containing both rigid and
6non-rigid atoms, with the stipulation that no atoms overlap.
7Produces a new combined omd file. The output omd file must be
8edited to run properly in OpenMD. Note that the two boxes must
9have identical box geometries (specified on the Hmat line).
10
11Usage: omdCombined
12
13Options:
14 -h, --help show this help
15 -u, --solute=... use specified OpenMD (.omd) file as the solute
16 -v, --solvent=... use specified OpenMD (.omd) file as the solvent
17 -o, --output-file=... use specified output (.omd) file
18
19
20Example:
21 omdCombined -u solute.omd -v solvent.omd -o combined.omd
22
23"""
24
25__author__ = "Sydney Shavalier (sshavali@nd.edu)"
26__version__ = "$Revision$"
27__date__ = "$Date$"
28__copyright__ = "Copyright (c) 2004-present The University of Notre Dame. All Rights Reserved."
29__license__ = "OpenMD"
30
31import sys
32import getopt
33import string
34import math
35import random
36
37_haveMDFileName1 = 0
38_haveMDFileName2 = 0
39_haveOutputFileName = 0
40
41metaData1 = []
42frameData1 = []
43positions1 = []
44velocities1 = []
45quaternions1 = []
46angVels1 = []
47indices1 = []
48Hmat1 = []
49BoxInv1 = []
50pvqj1 = []
51
52metaData2 = []
53frameData2 = []
54positions2 = []
55velocities2 = []
56quaternions2 = []
57angVels2 = []
58indices2 = []
59Hmat2 = []
60BoxInv2 = []
61pvqj2 = []
62componentLines = []
63ensembleLines = []
64
65def usage():
66 print(__doc__)
67
68def readFile1(mdFileName):
69 mdFile = open(mdFileName, 'r')
70 # Find OpenMD version info first
71 line = mdFile.readline()
72 while True:
73 if '<OpenMD version=' in line or '<OOPSE version=' in line:
74 OpenMDversion = line
75 break
76 line = mdFile.readline()
77
78 mdFile.seek(0)
79
80
81 line = mdFile.readline()
82
83 print("reading solute MetaData")
84 while True:
85 if '<MetaData>' in line:
86 while 2:
87 line = mdFile.readline()
88 if 'molecule' in line:
89 while not 'component' in line:
90 componentLines.append(line)
91 line = mdFile.readline()
92 componentLines.append("\n")
93 if 'component' in line:
94 while not '}' in line:
95 componentLines.append(line)
96 line = mdFile.readline()
97 componentLines.append("}\n\n")
98 if 'ensemble' in line:
99 while not '</MetaData>' in line:
100 ensembleLines.append(line)
101 line = mdFile.readline()
102 metaData1.append(line)
103
104 if '</MetaData>' in line:
105 metaData1.append(line)
106 break
107 break
108 line = mdFile.readline()
109
110 mdFile.seek(0)
111
112 print("reading solute Snapshot")
113 line = mdFile.readline()
114 while True:
115 if '<Snapshot>' in line:
116 line = mdFile.readline()
117 while True:
118 print("reading solute FrameData")
119 if '<FrameData>' in line:
120 while 2:
121 frameData1.append(line)
122 if 'Hmat:' in line:
123 L = line.split()
124 Hxx = float(L[2].strip(','))
125 Hxy = float(L[3].strip(','))
126 Hxz = float(L[4].strip(','))
127 Hyx = float(L[7].strip(','))
128 Hyy = float(L[8].strip(','))
129 Hyz = float(L[9].strip(','))
130 Hzx = float(L[12].strip(','))
131 Hzy = float(L[13].strip(','))
132 Hzz = float(L[14].strip(','))
133 Hmat1.append([Hxx, Hxy, Hxz])
134 Hmat1.append([Hyx, Hyy, Hyz])
135 Hmat1.append([Hzx, Hzy, Hzz])
136 BoxInv1.append(1.0/Hxx)
137 BoxInv1.append(1.0/Hyy)
138 BoxInv1.append(1.0/Hzz)
139 line = mdFile.readline()
140 if '</FrameData>' in line:
141 frameData1.append(line)
142 break
143 break
144
145 line = mdFile.readline()
146 while True:
147 if '<StuntDoubles>' in line:
148 line = mdFile.readline()
149 while 2:
150 L = line.split()
151 myIndex = int(L[0])
152 indices1.append(myIndex)
153 pvqj1.append(L[1])
154 x = float(L[2])
155 y = float(L[3])
156 z = float(L[4])
157 positions1.append([x, y, z])
158 vx = float(L[5])
159 vy = float(L[6])
160 vz = float(L[7])
161 velocities1.append([vx, vy, vz])
162 if 'pvqj' in L[1]:
163 qw = float(L[8])
164 qx = float(L[9])
165 qy = float(L[10])
166 qz = float(L[11])
167 quaternions1.append([qw, qx, qy, qz])
168 jx = float(L[12])
169 jy = float(L[13])
170 jz = float(L[14])
171 angVels1.append([jx, jy, jz])
172 else:
173 quaternions1.append([0.0, 0.0, 0.0, 0.0])
174 angVels1.append([0.0, 0.0, 0.0])
175
176 line = mdFile.readline()
177 if '</StuntDoubles>' in line:
178 break
179 break
180 line = mdFile.readline()
181 if not line: break
182
183 mdFile.close()
184
185def readFile2(mdFileName):
186 mdFile = open(mdFileName, 'r')
187 # Find OpenMD version info first
188 line = mdFile.readline()
189 while True:
190 if '<OpenMD version=' in line or '<OOPSE version=':
191 OpenMDversion = line
192 break
193 line = mdFile.readline()
194
195 mdFile.seek(0)
196
197 line = mdFile.readline()
198
199 print("reading solvent MetaData")
200 while True:
201 if '<MetaData>' in line:
202 while 2:
203 line = mdFile.readline()
204 if 'molecule' in line:
205 while not 'component' in line:
206 componentLines.append(line)
207 line = mdFile.readline()
208 componentLines.append("\n")
209 if 'component' in line:
210 while not '}' in line:
211 componentLines.append(line)
212 line = mdFile.readline()
213 componentLines.append("}\n\n")
214 metaData2.append(line)
215 if '</MetaData>' in line:
216 metaData2.append(line)
217 break
218 break
219 line = mdFile.readline()
220
221 mdFile.seek(0)
222
223 print("reading solvent Snapshot")
224 line = mdFile.readline()
225 while True:
226 if '<Snapshot>' in line:
227 line = mdFile.readline()
228 while True:
229 print("reading solvent FrameData")
230 if '<FrameData>' in line:
231 while 2:
232 frameData2.append(line)
233 if 'Hmat:' in line:
234 L = line.split()
235 Hxx = float(L[2].strip(','))
236 Hxy = float(L[3].strip(','))
237 Hxz = float(L[4].strip(','))
238 Hyx = float(L[7].strip(','))
239 Hyy = float(L[8].strip(','))
240 Hyz = float(L[9].strip(','))
241 Hzx = float(L[12].strip(','))
242 Hzy = float(L[13].strip(','))
243 Hzz = float(L[14].strip(','))
244 Hmat2.append([Hxx, Hxy, Hxz])
245 Hmat2.append([Hyx, Hyy, Hyz])
246 Hmat2.append([Hzx, Hzy, Hzz])
247 BoxInv2.append(1.0/Hxx)
248 BoxInv2.append(1.0/Hyy)
249 BoxInv2.append(1.0/Hzz)
250 line = mdFile.readline()
251 if '</FrameData>' in line:
252 frameData2.append(line)
253 break
254 break
255
256 line = mdFile.readline()
257 while True:
258 if '<StuntDoubles>' in line:
259 line = mdFile.readline()
260 while 2:
261 L = line.split()
262 myIndex = int(L[0])
263 indices2.append(myIndex)
264 pvqj2.append(L[1])
265 x = float(L[2])
266 y = float(L[3])
267 z = float(L[4])
268 positions2.append([x, y, z])
269 vx = float(L[5])
270 vy = float(L[6])
271 vz = float(L[7])
272 velocities2.append([vx, vy, vz])
273 if 'pvqj' in L[1]:
274 qw = float(L[8])
275 qx = float(L[9])
276 qy = float(L[10])
277 qz = float(L[11])
278 quaternions2.append([qw, qx, qy, qz])
279 jx = float(L[12])
280 jy = float(L[13])
281 jz = float(L[14])
282 angVels2.append([jx, jy, jz])
283 else:
284 quaternions2.append([0.0, 0.0, 0.0, 0.0])
285 angVels1.append([0.0, 0.0, 0.0])
286
287 line = mdFile.readline()
288 if '</StuntDoubles>' in line:
289 break
290 break
291 line = mdFile.readline()
292 if not line: break
293
294 mdFile.close()
295
296def writeFile(outputFileName):
297 outputFile = open(outputFileName, 'w')
298
299 outputFile.write("<OpenMD version=1>\n")
300 outputFile.write(" <MetaData>\n")
301 for componentLine in componentLines:
302 outputFile.write(componentLine)
303 for ensembleLine in ensembleLines:
304 outputFile.write(ensembleLine)
305 outputFile.write("\n")
306 outputFile.write(" </MetaData>\n")
307 outputFile.write(" <Snapshot>\n")
308 for frameline in frameData1:
309 outputFile.write(frameline)
310 outputFile.write(" <StuntDoubles>\n")
311
312
313 newIndex = 0
314 for i in range(len(indices1)):
315 if (pvqj1[i] == 'pv'):
316 outputFile.write("%10d %7s %18.10g %18.10g %18.10g %14e %13e %13e\n" % (newIndex, pvqj1[i], positions1[i][0], positions1[i][1], positions1[i][2], velocities1[i][0], velocities1[i][1], velocities1[i][2]))
317 elif(pvqj1[i] == 'pvqj'):
318 outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e %13e %13e %13e %13e %13e %13e %13e\n" % (newIndex, pvqj1[i], positions1[i][0], positions1[i][1], positions1[i][2], velocities1[i][0], velocities1[i][1], velocities1[i][2], quaternions1[i][0], quaternions1[i][1], quaternions1[i][2], quaternions1[i][3], angVels1[i][0], angVels1[i][1], angVels1[i][2]))
319
320 newIndex+=1
321
322 newIndexAlt = newIndex
323 for j in range(len(indices2)):
324 if (pvqj2[j] == 'pv'):
325 outputFile.write("%10d %7s %18.10g %18.10g %18.10g %14e %13e %13e\n" % (newIndexAlt, pvqj2[j], positions2[j][0], positions2[j][1], positions2[j][2], velocities2[j][0], velocities2[j][1], velocities2[j][2]))
326 elif(pvqj2[j] == 'pvqj'):
327 outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e %13e %13e %13e %13e %13e %13e %13e\n" % (newIndexAlt, pvqj2[j], positions2[j][0], positions2[j][1], positions2[j][2], velocities2[j][0], velocities2[j][1], velocities2[j][2], quaternions2[j][0], quaternions2[j][1], quaternions2[j][2], quaternions2[j][3], angVels2[j][0], angVels2[j][1], angVels2[j][2]))
328
329 newIndexAlt+=1
330
331 outputFile.write(" </StuntDoubles>\n")
332 outputFile.write(" </Snapshot>\n")
333 outputFile.write("</OpenMD>\n")
334 outputFile.close()
335
336
337def checkBoxes():
338 boxTolerance = 1.0e-3
339 maxDiff = 0.0
340 for i in range(3):
341 for j in range(3):
342 diff = math.fabs( Hmat1[i][j] - Hmat2[i][j])
343 if (diff > maxDiff):
344 maxDiff = diff
345 if (maxDiff > boxTolerance):
346 print("The solute and solvent boxes have different geometries:")
347 print(" Solute | Solvent")
348 print(" -------------------------------------|------------------------------------")
349 for i in range(3):
350 print(( "| %10.4g %10.4g %10.4g | %10.4g %10.4g %10.4g |" % (Hmat1[i][0], Hmat1[i][1], Hmat1[i][2], Hmat2[i][0], Hmat2[i][1], Hmat2[i][2])))
351
352 print(" -------------------------------------|------------------------------------")
353 sys.exit()
354
355def checkOverlap(sd1_index,sd2_index): #checks to see if solute and solvent atoms overlap (within 1 Angstrom of each other)
356 overlapTolerance = 1.0
357 diff = 0.0
358 delta = []
359 for i in range(3):
360 delta.append(positions1[sd1_index][i] - positions2[sd2_index][i])
361 wrapVector(delta)
362 diff = delta[0]**2 + delta[1]**2 + delta[2]**2
363 diff = math.sqrt(diff)
364 if (diff < overlapTolerance):
365 print("Warning: Solute and solvent atoms overlap! No output file was created.")
366 sys.exit()
367
368
369def roundMe(x):
370 if (x >= 0.0):
371 return math.floor(x + 0.5)
372 else:
373 return math.ceil(x - 0.5)
374
375def frange(start,stop,step=1.0):
376 while start < stop:
377 yield start
378 start += step
379
380
381def wrapVector(myVect):
382 scaled = [0.0, 0.0, 0.0]
383 for i in range(3):
384 scaled[i] = myVect[i] * BoxInv1[i]
385 scaled[i] = scaled[i] - roundMe(scaled[i])
386 myVect[i] = scaled[i] * Hmat1[i][i]
387 return myVect
388
389
390def main(argv):
391 try:
392 opts, args = getopt.getopt(argv, "hu:v:o:", ["help", "solute=", "solvent=", "output-file="])
393 except getopt.GetoptError:
394 usage()
395 sys.exit(2)
396 for opt, arg in opts:
397 if opt in ("-h", "--help"):
398 usage()
399 sys.exit()
400 elif opt in ("-u", "--solute"):
401 mdFileName1 = arg
402 global _haveMDFileName1
403 _haveMDFileName1 = 1
404 elif opt in ("-v", "--solvent"):
405 mdFileName2 = arg
406 global _haveMDFileName2
407 _haveMDFileName2 = 1
408 elif opt in ("-o", "--output-file"):
409 outputFileName = arg
410 global _haveOutputFileName
411 _haveOutputFileName = 1
412
413 if (_haveMDFileName1 != 1):
414 usage()
415 print("No OpenMD (omd) file was specified for the solute")
416 sys.exit()
417
418 if (_haveMDFileName2 != 1):
419 usage()
420 print("No OpenMD (omd) file was specified for the solvent")
421 sys.exit()
422
423 if (_haveOutputFileName != 1):
424 usage()
425 print("No output file was specified")
426 sys.exit()
427
428 readFile1(mdFileName1)
429 readFile2(mdFileName2)
430 checkBoxes()
431 for sd1 in range(len(indices1)):
432 for sd2 in range(len(indices2)):
433 checkOverlap(sd1,sd2)
434 writeFile(outputFileName)
435
436if __name__ == "__main__":
437 if len(sys.argv) == 1:
438 usage()
439 sys.exit()
440 main(sys.argv[1:])