OpenMD 3.2
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
omd-solvator
1#!@Python3_EXECUTABLE@
2"""OMD Solvator
3
4Opens two omd files, one with a solute structure and one with a
5solvent structure. Deletes any solvent molecules that overlap with
6solute molecules and produces a new combined omd file. The output omd
7file must be edited to run properly in OpenMD. Note that the two
8boxes must have identical box geometries (specified on the Hmat line).
9
10Usage: omd-solvator
11
12Options:
13 -h, --help show this help
14 -u, --solute=... use specified OpenMD (.omd) file as the solute
15 -v, --solvent=... use specified OpenMD (.omd) file as the solvent
16 -r, --rcut=... specify the cutoff radius for deleting solvent
17 -o, --output-file=... use specified output (.omd) file
18 -n, --nSoluteAtoms=... Number of atoms in solute molecule,
19 default is 1 atom.
20 -p, --nSolventAtoms=... Number of atoms in solvent molecule,
21 default is 1 atom.
22
23Example:
24 omd-solvator -u solute.omd -v solvent.omd -n 3 -p 3 -r 4.0 -o combined.omd
25
26"""
27
28__author__ = "Charles Vardeman (cvardema@nd.edu)"
29__version__ = "$Revision$"
30__date__ = "$Date$"
31__copyright__ = "Copyright (c) 2004-present The University of Notre Dame. All Rights Reserved."
32__license__ = "OpenMD"
33
34import sys
35import getopt
36import string
37import math
38import random
39
40_haveMDFileName1 = 0
41_haveMDFileName2 = 0
42_haveRcut = 0
43_haveOutputFileName = 0
44_haveNSoluteAtoms = 0
45_haveNSolventAtoms = 0
46
47metaData1 = []
48frameData1 = []
49positions1 = []
50velocities1 = []
51quaternions1 = []
52angVels1 = []
53indices1 = []
54Hmat1 = []
55BoxInv1 = []
56pvqj1 = []
57
58metaData2 = []
59frameData2 = []
60positions2 = []
61velocities2 = []
62quaternions2 = []
63angVels2 = []
64indices2 = []
65Hmat2 = []
66BoxInv2 = []
67pvqj2 = []
68
69keepers = []
70
71soluteTypeLine = str()
72solventTypeLine = str()
73soluteMolLine = str()
74nSolvents = 0
75
76def usage():
77 print(__doc__)
78
79def readFile1(mdFileName):
80 mdFile = open(mdFileName, 'r')
81 # Find OpenMD version info first
82 line = mdFile.readline()
83 while True:
84 if '<OpenMD version=' in line or '<OOPSE version=' in line:
85 OpenMDversion = line
86 break
87 line = mdFile.readline()
88
89 # Rewind file and find start of MetaData block
90
91 mdFile.seek(0)
92 line = mdFile.readline()
93
94 print("reading solute MetaData")
95 while True:
96 if '<MetaData>' in line:
97 while 2:
98 metaData1.append(line)
99 line = mdFile.readline()
100 if 'type' in line:
101 global soluteTypeLine
102 soluteTypeLine = line
103 if 'nMol' in line:
104 global soluteMolLine
105 soluteMolLine = line
106 if '</MetaData>' in line:
107 metaData1.append(line)
108 break
109 break
110 line = mdFile.readline()
111
112 mdFile.seek(0)
113 print("reading solute Snapshot")
114 line = mdFile.readline()
115 while True:
116 if '<Snapshot>' in line:
117 line = mdFile.readline()
118 while True:
119 print("reading solute FrameData")
120 if '<FrameData>' in line:
121 while 2:
122 frameData1.append(line)
123 if 'Hmat:' in line:
124 L = line.split()
125 Hxx = float(L[2].strip(','))
126 Hxy = float(L[3].strip(','))
127 Hxz = float(L[4].strip(','))
128 Hyx = float(L[7].strip(','))
129 Hyy = float(L[8].strip(','))
130 Hyz = float(L[9].strip(','))
131 Hzx = float(L[12].strip(','))
132 Hzy = float(L[13].strip(','))
133 Hzz = float(L[14].strip(','))
134 Hmat1.append([Hxx, Hxy, Hxz])
135 Hmat1.append([Hyx, Hyy, Hyz])
136 Hmat1.append([Hzx, Hzy, Hzz])
137 BoxInv1.append(1.0/Hxx)
138 BoxInv1.append(1.0/Hyy)
139 BoxInv1.append(1.0/Hzz)
140 line = mdFile.readline()
141 if '</FrameData>' in line:
142 frameData1.append(line)
143 break
144 break
145
146 line = mdFile.readline()
147 while True:
148 if '<StuntDoubles>' in line:
149 line = mdFile.readline()
150 while 2:
151 L = line.split()
152 myIndex = int(L[0])
153 indices1.append(myIndex)
154 pvqj1.append(L[1])
155 x = float(L[2])
156 y = float(L[3])
157 z = float(L[4])
158 positions1.append([x, y, z])
159 vx = float(L[5])
160 vy = float(L[6])
161 vz = float(L[7])
162 velocities1.append([vx, vy, vz])
163 if 'pvqj' in L[1]:
164 qw = float(L[8])
165 qx = float(L[9])
166 qy = float(L[10])
167 qz = float(L[11])
168 quaternions1.append([qw, qx, qy, qz])
169 jx = float(L[12])
170 jy = float(L[13])
171 jz = float(L[14])
172 angVels1.append([jx, jy, jz])
173 else:
174 quaternions1.append([0.0, 0.0, 0.0, 0.0])
175 angVels1.append([0.0, 0.0, 0.0])
176
177 line = mdFile.readline()
178 if '</StuntDoubles>' in line:
179 break
180 break
181 line = mdFile.readline()
182 if not line: break
183
184 mdFile.close()
185
186def readFile2(mdFileName):
187 mdFile = open(mdFileName, 'r')
188 # Find OpenMD version info first
189 line = mdFile.readline()
190 while True:
191 if '<OpenMD version=' in line or '<OOPSE version=':
192 OpenMDversion = line
193 break
194 line = mdFile.readline()
195
196 # Rewind file and find start of MetaData block
197
198 mdFile.seek(0)
199 line = mdFile.readline()
200 print("reading solvent MetaData")
201 while True:
202 if '<MetaData>' in line:
203 while 2:
204 if 'type' in line:
205 global solventTypeLine
206 solventTypeLine = line
207 metaData2.append(line)
208 line = mdFile.readline()
209 if '</MetaData>' in line:
210 metaData2.append(line)
211 break
212 break
213 line = mdFile.readline()
214
215 mdFile.seek(0)
216 print("reading solvent Snapshot")
217 line = mdFile.readline()
218 while True:
219 if '<Snapshot>' in line:
220 line = mdFile.readline()
221 while True:
222 print("reading solvent FrameData")
223 if '<FrameData>' in line:
224 while 2:
225 frameData2.append(line)
226 if 'Hmat:' in line:
227 L = line.split()
228 Hxx = float(L[2].strip(','))
229 Hxy = float(L[3].strip(','))
230 Hxz = float(L[4].strip(','))
231 Hyx = float(L[7].strip(','))
232 Hyy = float(L[8].strip(','))
233 Hyz = float(L[9].strip(','))
234 Hzx = float(L[12].strip(','))
235 Hzy = float(L[13].strip(','))
236 Hzz = float(L[14].strip(','))
237 Hmat2.append([Hxx, Hxy, Hxz])
238 Hmat2.append([Hyx, Hyy, Hyz])
239 Hmat2.append([Hzx, Hzy, Hzz])
240 BoxInv2.append(1.0/Hxx)
241 BoxInv2.append(1.0/Hyy)
242 BoxInv2.append(1.0/Hzz)
243 line = mdFile.readline()
244 if '</FrameData>' in line:
245 frameData2.append(line)
246 break
247 break
248
249 line = mdFile.readline()
250 while True:
251 if '<StuntDoubles>' in line:
252 line = mdFile.readline()
253 while 2:
254 L = line.split()
255 myIndex = int(L[0])
256 indices2.append(myIndex)
257 pvqj2.append(L[1])
258 loc = 2
259 if 'p' in L[1]:
260 x = float(L[loc])
261 y = float(L[loc+1])
262 z = float(L[loc+2])
263 positions2.append([x, y, z])
264 loc += 3
265 else:
266 positions2.append([0.0, 0.0, 0.0])
267 if 'v' in L[1]:
268 vx = float(L[loc])
269 vy = float(L[loc+1])
270 vz = float(L[loc+2])
271 velocities2.append([vx, vy, vz])
272 loc += 3
273 else:
274 velocities2.append([0.0, 0.0, 0.0])
275 if 'q' in L[1]:
276 qw = float(L[loc])
277 qx = float(L[loc+1])
278 qy = float(L[loc+2])
279 qz = float(L[loc+3])
280 quaternions2.append([qw, qx, qy, qz])
281 loc += 4
282 else:
283 quaternions2.append([0.0, 0.0, 0.0, 0.0])
284 if 'j' in L[1]:
285 jx = float(L[loc])
286 jy = float(L[loc+1])
287 jz = float(L[loc+2])
288 angVels2.append([jx, jy, jz])
289 loc += 3
290 else:
291 angVels2.append([0.0, 0.0, 0.0])
292
293 line = mdFile.readline()
294 if '</StuntDoubles>' in line:
295 break
296 break
297 line = mdFile.readline()
298 if not line: break
299
300 mdFile.close()
301
302def writeFile(outputFileName):
303 outputFile = open(outputFileName, 'w')
304
305 outputFile.write("<OpenMD version=1>\n")
306
307# for metaline in metaData1:
308# outputFile.write(metaline)
309 outputFile.write(" <MetaData>\n")
310 outputFile.write("\n\n")
311 outputFile.write("component{\n")
312 outputFile.write(soluteTypeLine)
313 outputFile.write(soluteMolLine)
314 outputFile.write("}\n")
315
316 outputFile.write("component{\n")
317 outputFile.write(solventTypeLine)
318 outputFile.write("nMol = %d;\n" % (nSolvents))
319 outputFile.write("}\n")
320 outputFile.write("\n\n")
321 outputFile.write(" </MetaData>\n")
322 outputFile.write(" <Snapshot>\n")
323
324 for frameline in frameData1:
325 outputFile.write(frameline)
326
327 outputFile.write(" <StuntDoubles>\n")
328
329 newIndex = 0
330 for i in range(len(indices1)):
331 if ('p' in pvqj1[i]):
332 outputFile.write("%10d %7s %18.10g %18.10g %18.10g " % (newIndex, pvqj1[i], positions1[i][0], positions1[i][1], positions1[i][2]))
333 if ('v' in pvqj1[i]):
334 outputFile.write("%13e %13e %13e " % (velocities1[i][0], velocities1[i][1], velocities1[i][2]))
335 if('q' in pvqj1[i]):
336 outputFile.write("%13e %13e %13e %13e " % (quaternions1[i][0], quaternions1[i][1], quaternions1[i][2], quaternions1[i][3]))
337 if('j' in pvqj1[i]):
338 outputFile.write("%13e %13e %13e " % (angVels1[i][0], angVels1[i][1], angVels1[i][2]))
339 outputFile.write("\n")
340
341 newIndex = newIndex + 1
342
343 outputFile.write(" </StuntDoubles>\n")
344 outputFile.write(" </Snapshot>\n")
345 outputFile.write("</OpenMD>\n")
346 outputFile.close()
347
348def checkBoxes():
349 boxTolerance = 1.0e-3
350 maxDiff = 0.0
351 for i in range(3):
352 for j in range(3):
353 diff = math.fabs( Hmat1[i][j] - Hmat2[i][j])
354 if (diff > maxDiff):
355 maxDiff = diff
356 if (maxDiff > boxTolerance):
357 print("The solute and solvent boxes have different geometries:")
358 print(" Solute | Solvent")
359 print(" -------------------------------------|------------------------------------")
360 for i in range(3):
361 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])))
362
363 print(" -------------------------------------|------------------------------------")
364 sys.exit()
365
366
367def roundMe(x):
368 if (x >= 0.0):
369 return math.floor(x + 0.5)
370 else:
371 return math.ceil(x - 0.5)
372
373def frange(start,stop,step=1.0):
374 while start < stop:
375 yield start
376 start += step
377
378
379def wrapVector(myVect):
380 scaled = [0.0, 0.0, 0.0]
381 for i in range(3):
382 scaled[i] = myVect[i] * BoxInv1[i]
383 scaled[i] = scaled[i] - roundMe(scaled[i])
384 myVect[i] = scaled[i] * Hmat1[i][i]
385 return myVect
386
387def dot(L1, L2):
388 myDot = 0.0
389 for i in range(len(L1)):
390 myDot = myDot + L1[i]*L2[i]
391 return myDot
392
393def normalize(L1):
394 L2 = []
395 myLength = math.sqrt(dot(L1, L1))
396 for i in range(len(L1)):
397 L2.append(L1[i] / myLength)
398 return L2
399
400def cross(L1, L2):
401 # don't call this with anything other than length 3 lists please
402 # or you'll be sorry
403 L3 = [0.0, 0.0, 0.0]
404 L3[0] = L1[1]*L2[2] - L1[2]*L2[1]
405 L3[1] = L1[2]*L2[0] - L1[0]*L2[2]
406 L3[2] = L1[0]*L2[1] - L1[1]*L2[0]
407 return L3
408
409def removeOverlaps(rcut, nSolventAtoms, nSoluteAtoms):
410
411 rcut2 = rcut*rcut
412 nextMol = 0
413 for i in range(0, len(indices2), nSolventAtoms):
414 keepThisMolecule = 1
415 for atom1 in range (i, (i+nSolventAtoms)):
416
417 iPos = positions2[atom1]
418 for j in range(0, len(indices1)):
419 for atom2 in range (j, (j+nSoluteAtoms), nSoluteAtoms):
420 jPos = positions1[atom2]
421 dpos = [jPos[0]-iPos[0], jPos[1]-iPos[1], jPos[2]-iPos[2]]
422 dpos = wrapVector(dpos)
423 dist2 = dot(dpos, dpos)
424
425
426 if (dist2 < rcut2):
427 keepThisMolecule = 0
428 break
429 if (keepThisMolecule == 0):
430 break
431
432 keepers.append(keepThisMolecule)
433
434
435 global nSolvents
436 myIndex = len(indices2) - 1
437 for i in range(0, len(keepers)):
438
439 if (keepers[i] == 1):
440 nSolvents = nSolvents + 1
441 atomStartIndex = i * nSolventAtoms
442 for j in range (atomStartIndex, (atomStartIndex+nSolventAtoms)):
443 indices1.append(myIndex)
444 pvqj1.append(pvqj2[j])
445 if ("p" in pvqj2[j]):
446 positions1.append(positions2[j])
447 else:
448 positions1.append([0.0, 0.0, 0.0])
449 if ("v" in pvqj2[j]):
450 velocities1.append(velocities2[j])
451 else:
452 velocities1.append([0.0, 0.0, 0.0])
453 if ("q" in pvqj2[j]):
454 quaternions1.append(quaternions2[j])
455 else:
456 quaternions1.append([0.0, 0.0, 0.0, 0.0])
457 if ("j" in pvqj2[j]):
458 angVels1.append(angVels2[j])
459 else:
460 angVels1.append([0.0, 0.0, 0.0])
461
462 myIndex = myIndex +1
463
464def main(argv):
465 try:
466 opts, args = getopt.getopt(argv, "hu:v:n:p:r:o:", ["help", "solute=", "solvent=", "nSoluteAtoms=", "nSolventAtoms=", "rcut=" "output-file="])
467 except getopt.GetoptError:
468 usage()
469 sys.exit(2)
470 for opt, arg in opts:
471 if opt in ("-h", "--help"):
472 usage()
473 sys.exit()
474 elif opt in ("-u", "--solute"):
475 mdFileName1 = arg
476 global _haveMDFileName1
477 _haveMDFileName1 = 1
478 elif opt in ("-v", "--solvent"):
479 mdFileName2 = arg
480 global _haveMDFileName2
481 _haveMDFileName2 = 1
482 elif opt in ("-n", "--nSoluteAtoms"):
483 nSoluteAtoms = int(arg)
484 global _haveNSoluteAtoms
485 _haveNSoluteAtoms = 1
486 elif opt in ("-p", "--nSolventAtoms"):
487 nSolventAtoms = int(arg)
488 global _haveNSolventAtoms
489 _haveNSolventAtoms = 1
490 elif opt in ("-r", "--rcut"):
491 rcut = float(arg)
492 global _haveRcut
493 _haveRcut = 1
494 elif opt in ("-o", "--output-file"):
495 outputFileName = arg
496 global _haveOutputFileName
497 _haveOutputFileName = 1
498
499 if (_haveMDFileName1 != 1):
500 usage()
501 print("No OpenMD (omd) file was specified for the solute")
502 sys.exit()
503
504 if (_haveMDFileName2 != 1):
505 usage()
506 print("No OpenMD (omd) file was specified for the solvent")
507 sys.exit()
508
509 if (_haveOutputFileName != 1):
510 usage()
511 print("No output file was specified")
512 sys.exit()
513
514 if (_haveRcut != 1):
515 print("No cutoff radius was specified, using 4 angstroms")
516 rcut =4.0
517
518 if (_haveNSoluteAtoms != 1):
519 print("Number of solute atoms was not specified. Using 1 atom.")
520 nSoluteAtoms = 1
521
522 if (_haveNSolventAtoms != 1):
523 print("Number of solute atoms was not specified. Using 1 atom.")
524 nSolventAtoms = 1
525
526 readFile1(mdFileName1)
527 readFile2(mdFileName2)
528 checkBoxes()
529 removeOverlaps(rcut, nSolventAtoms, nSoluteAtoms)
530 writeFile(outputFileName)
531
532if __name__ == "__main__":
533 if len(sys.argv) == 1:
534 usage()
535 sys.exit()
536 main(sys.argv[1:])