OpenMD 3.1
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 x = float(L[2])
259 y = float(L[3])
260 z = float(L[4])
261 positions2.append([x, y, z])
262 vx = float(L[5])
263 vy = float(L[6])
264 vz = float(L[7])
265 velocities2.append([vx, vy, vz])
266 if 'pvqj' in L[1]:
267 qw = float(L[8])
268 qx = float(L[9])
269 qy = float(L[10])
270 qz = float(L[11])
271 quaternions2.append([qw, qx, qy, qz])
272 jx = float(L[12])
273 jy = float(L[13])
274 jz = float(L[14])
275 angVels2.append([jx, jy, jz])
276 else:
277 quaternions2.append([0.0, 0.0, 0.0, 0.0])
278 angVels1.append([0.0, 0.0, 0.0])
279
280 line = mdFile.readline()
281 if '</StuntDoubles>' in line:
282 break
283 break
284 line = mdFile.readline()
285 if not line: break
286
287 mdFile.close()
288
289def writeFile(outputFileName):
290 outputFile = open(outputFileName, 'w')
291
292 outputFile.write("<OpenMD version=1>\n")
293
294# for metaline in metaData1:
295# outputFile.write(metaline)
296 outputFile.write(" <MetaData>\n")
297 outputFile.write("\n\n")
298 outputFile.write("component{\n")
299 outputFile.write(soluteTypeLine)
300 outputFile.write(soluteMolLine)
301 outputFile.write("}\n")
302
303 outputFile.write("component{\n")
304 outputFile.write(solventTypeLine)
305 outputFile.write("nMol = %d;\n" % (nSolvents))
306 outputFile.write("}\n")
307 outputFile.write("\n\n")
308 outputFile.write(" </MetaData>\n")
309 outputFile.write(" <Snapshot>\n")
310
311 for frameline in frameData1:
312 outputFile.write(frameline)
313
314 outputFile.write(" <StuntDoubles>\n")
315
316
317 newIndex = 0
318 for i in range(len(indices1)):
319 if (pvqj1[i] == 'pv'):
320 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]))
321 elif(pvqj1[i] == 'pvqj'):
322 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]))
323
324 newIndex = newIndex + 1
325
326 outputFile.write(" </StuntDoubles>\n")
327 outputFile.write(" </Snapshot>\n")
328 outputFile.write("</OpenMD>\n")
329 outputFile.close()
330
331def checkBoxes():
332 boxTolerance = 1.0e-3
333 maxDiff = 0.0
334 for i in range(3):
335 for j in range(3):
336 diff = math.fabs( Hmat1[i][j] - Hmat2[i][j])
337 if (diff > maxDiff):
338 maxDiff = diff
339 if (maxDiff > boxTolerance):
340 print("The solute and solvent boxes have different geometries:")
341 print(" Solute | Solvent")
342 print(" -------------------------------------|------------------------------------")
343 for i in range(3):
344 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])))
345
346 print(" -------------------------------------|------------------------------------")
347 sys.exit()
348
349
350def roundMe(x):
351 if (x >= 0.0):
352 return math.floor(x + 0.5)
353 else:
354 return math.ceil(x - 0.5)
355
356def frange(start,stop,step=1.0):
357 while start < stop:
358 yield start
359 start += step
360
361
362def wrapVector(myVect):
363 scaled = [0.0, 0.0, 0.0]
364 for i in range(3):
365 scaled[i] = myVect[i] * BoxInv1[i]
366 scaled[i] = scaled[i] - roundMe(scaled[i])
367 myVect[i] = scaled[i] * Hmat1[i][i]
368 return myVect
369
370def dot(L1, L2):
371 myDot = 0.0
372 for i in range(len(L1)):
373 myDot = myDot + L1[i]*L2[i]
374 return myDot
375
376def normalize(L1):
377 L2 = []
378 myLength = math.sqrt(dot(L1, L1))
379 for i in range(len(L1)):
380 L2.append(L1[i] / myLength)
381 return L2
382
383def cross(L1, L2):
384 # don't call this with anything other than length 3 lists please
385 # or you'll be sorry
386 L3 = [0.0, 0.0, 0.0]
387 L3[0] = L1[1]*L2[2] - L1[2]*L2[1]
388 L3[1] = L1[2]*L2[0] - L1[0]*L2[2]
389 L3[2] = L1[0]*L2[1] - L1[1]*L2[0]
390 return L3
391
392def removeOverlaps(rcut, nSolventAtoms, nSoluteAtoms):
393
394 rcut2 = rcut*rcut
395 nextMol = 0
396 for i in range(0, len(indices2), nSolventAtoms):
397 keepThisMolecule = 1
398 for atom1 in range (i, (i+nSolventAtoms)):
399
400 iPos = positions2[atom1]
401 for j in range(0, len(indices1)):
402 for atom2 in range (j, (j+nSoluteAtoms), nSoluteAtoms):
403 jPos = positions1[atom2]
404 dpos = [jPos[0]-iPos[0], jPos[1]-iPos[1], jPos[2]-iPos[2]]
405 dpos = wrapVector(dpos)
406 dist2 = dot(dpos, dpos)
407
408
409 if (dist2 < rcut2):
410 keepThisMolecule = 0
411 break
412 if (keepThisMolecule == 0):
413 break
414
415 keepers.append(keepThisMolecule)
416
417
418 global nSolvents
419 myIndex = len(indices2) - 1
420 for i in range(0, len(keepers)):
421
422 if (keepers[i] == 1):
423 nSolvents = nSolvents + 1
424 atomStartIndex = i * nSolventAtoms
425 for j in range (atomStartIndex, (atomStartIndex+nSolventAtoms)):
426 indices1.append(myIndex)
427 pvqj1.append(pvqj2[j])
428 if (pvqj2[j] == 'pv'):
429 positions1.append(positions2[j])
430 velocities1.append(velocities2[j])
431 quaternions1.append([0.0, 0.0, 0.0, 0.0])
432 angVels1.append([0.0, 0.0, 0.0])
433 else:
434 positions1.append(positions2[j])
435 velocities1.append(velocities2[j])
436 quaternions1.append(quaternions2[j])
437 angVels1.append(angVels2[j])
438 # indices1.append(indices2[j])
439 myIndex = myIndex +1
440
441def main(argv):
442 try:
443 opts, args = getopt.getopt(argv, "hu:v:n:p:r:o:", ["help", "solute=", "solvent=", "nSoluteAtoms=", "nSolventAtoms=", "rcut=" "output-file="])
444 except getopt.GetoptError:
445 usage()
446 sys.exit(2)
447 for opt, arg in opts:
448 if opt in ("-h", "--help"):
449 usage()
450 sys.exit()
451 elif opt in ("-u", "--solute"):
452 mdFileName1 = arg
453 global _haveMDFileName1
454 _haveMDFileName1 = 1
455 elif opt in ("-v", "--solvent"):
456 mdFileName2 = arg
457 global _haveMDFileName2
458 _haveMDFileName2 = 1
459 elif opt in ("-n", "--nSoluteAtoms"):
460 nSoluteAtoms = int(arg)
461 global _haveNSoluteAtoms
462 _haveNSoluteAtoms = 1
463 elif opt in ("-p", "--nSolventAtoms"):
464 nSolventAtoms = int(arg)
465 global _haveNSolventAtoms
466 _haveNSolventAtoms = 1
467 elif opt in ("-r", "--rcut"):
468 rcut = float(arg)
469 global _haveRcut
470 _haveRcut = 1
471 elif opt in ("-o", "--output-file"):
472 outputFileName = arg
473 global _haveOutputFileName
474 _haveOutputFileName = 1
475
476 if (_haveMDFileName1 != 1):
477 usage()
478 print("No OpenMD (omd) file was specified for the solute")
479 sys.exit()
480
481 if (_haveMDFileName2 != 1):
482 usage()
483 print("No OpenMD (omd) file was specified for the solvent")
484 sys.exit()
485
486 if (_haveOutputFileName != 1):
487 usage()
488 print("No output file was specified")
489 sys.exit()
490
491 if (_haveRcut != 1):
492 print("No cutoff radius was specified, using 4 angstroms")
493 rcut =4.0
494
495 if (_haveNSoluteAtoms != 1):
496 print("Number of solute atoms was not specified. Using 1 atom.")
497 nSoluteAtoms = 1
498
499 if (_haveNSolventAtoms != 1):
500 print("Number of solute atoms was not specified. Using 1 atom.")
501 nSolventAtoms = 1
502
503 readFile1(mdFileName1)
504 readFile2(mdFileName2)
505 checkBoxes()
506 removeOverlaps(rcut, nSolventAtoms, nSoluteAtoms)
507 writeFile(outputFileName)
508
509if __name__ == "__main__":
510 if len(sys.argv) == 1:
511 usage()
512 sys.exit()
513 main(sys.argv[1:])