OpenMD 3.0
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(wrapVector([x, y, z])) #wraps positions back into periodic box
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(wrapVector([x, y, z])) #wraps positions back into periodic box
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 for i in range(3):
359 diff += (math.fabs(positions1[sd1_index][i]-positions2[sd2_index][i]))**2
360 diff = math.sqrt(diff)
361 if (diff < overlapTolerance):
362 print("Warning: Solute and solvent atoms overlap! No output file was created.")
363 sys.exit()
364
365
366def roundMe(x):
367 if (x >= 0.0):
368 return math.floor(x + 0.5)
369 else:
370 return math.ceil(x - 0.5)
371
372def frange(start,stop,step=1.0):
373 while start < stop:
374 yield start
375 start += step
376
377
378def wrapVector(myVect):
379 scaled = [0.0, 0.0, 0.0]
380 for i in range(3):
381 scaled[i] = myVect[i] * BoxInv1[i]
382 scaled[i] = scaled[i] - roundMe(scaled[i])
383 myVect[i] = scaled[i] * Hmat1[i][i]
384 return myVect
385
386
387def main(argv):
388 try:
389 opts, args = getopt.getopt(argv, "hu:v:o:", ["help", "solute=", "solvent=", "output-file="])
390 except getopt.GetoptError:
391 usage()
392 sys.exit(2)
393 for opt, arg in opts:
394 if opt in ("-h", "--help"):
395 usage()
396 sys.exit()
397 elif opt in ("-u", "--solute"):
398 mdFileName1 = arg
399 global _haveMDFileName1
400 _haveMDFileName1 = 1
401 elif opt in ("-v", "--solvent"):
402 mdFileName2 = arg
403 global _haveMDFileName2
404 _haveMDFileName2 = 1
405 elif opt in ("-o", "--output-file"):
406 outputFileName = arg
407 global _haveOutputFileName
408 _haveOutputFileName = 1
409
410 if (_haveMDFileName1 != 1):
411 usage()
412 print("No OpenMD (omd) file was specified for the solute")
413 sys.exit()
414
415 if (_haveMDFileName2 != 1):
416 usage()
417 print("No OpenMD (omd) file was specified for the solvent")
418 sys.exit()
419
420 if (_haveOutputFileName != 1):
421 usage()
422 print("No output file was specified")
423 sys.exit()
424
425 readFile1(mdFileName1)
426 readFile2(mdFileName2)
427 checkBoxes()
428 for sd1 in range(len(indices1)):
429 for sd2 in range(len(indices2)):
430 checkOverlap(sd1,sd2)
431 writeFile(outputFileName)
432
433if __name__ == "__main__":
434 if len(sys.argv) == 1:
435 usage()
436 sys.exit()
437 main(sys.argv[1:])