ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/applications/utilities/md-solvator
Revision: 1266
Committed: Fri Jun 27 15:47:06 2008 UTC (17 years ago) by chuckv
Original Path: trunk/src/applications/utilities/md-solvator
File size: 16511 byte(s)
Log Message:
Md solvator writes out md files now. Index problem.

File Contents

# User Rev Content
1 kfletch2 1262 #!/usr/bin/env python
2     """Ice Cube Solvator
3 gezelter 1116
4 kfletch2 1262 Opens two md files, one with water in an ice structure,
5     and one with water in a liquid phase. Deletes any overlapping
6     liquid molecules and merges the two md files.
7 gezelter 1116
8 kfletch2 1262 Usage: waterRotator
9 gezelter 1116
10 kfletch2 1262 Options:
11     -h, --help show this help
12     -u, --solute=... use specified meta-data (.md) file as the solute
13     -v, --solvent=... use specified meta-data (.md) file as the solvent
14     -r, --rcut=... specify the cutoff radius for deleting solvent
15     -o, --output-file=... use specified output (.md) file
16 gezelter 1116
17    
18 kfletch2 1262 Example:
19     iceCubeSolvator -u frosty.md -v tepid.md -r 4.0 -o lukewarm.md
20 gezelter 1116
21 kfletch2 1262 """
22 gezelter 1116
23 kfletch2 1262 __author__ = "Dan Gezelter (gezelter@nd.edu)"
24 chuckv 1266 __version__ = "$Revision: 1.8 $"
25     __date__ = "$Date: 2008-06-27 15:47:06 $"
26 kfletch2 1262 __copyright__ = "Copyright (c) 2006 by the University of Notre Dame"
27     __license__ = "OOPSE"
28 gezelter 1116
29 kfletch2 1262 import sys
30     import getopt
31     import string
32     import math
33     import random
34     from sets import *
35     #from Numeric import *
36 gezelter 1116
37 kfletch2 1262 _haveMDFileName1 = 0
38     _haveMDFileName2 = 0
39     _haveRcut = 0
40     _haveOutputFileName = 0
41 chuckv 1265 _haveNSoluteAtoms = 1
42     _haveNSolventAtoms = 1
43 gezelter 1116
44 kfletch2 1262 metaData1 = []
45     frameData1 = []
46     positions1 = []
47     velocities1 = []
48     quaternions1 = []
49     angVels1 = []
50     indices1 = []
51     Hmat1 = []
52     BoxInv1 = []
53 chuckv 1263 pvqj1 = []
54 gezelter 1116
55 kfletch2 1262 metaData2 = []
56     frameData2 = []
57     positions2 = []
58     velocities2 = []
59     quaternions2 = []
60     angVels2 = []
61     indices2 = []
62     Hmat2 = []
63     BoxInv2 = []
64 chuckv 1263 pvqj2 = []
65 gezelter 1116
66 kfletch2 1262 keepers = []
67 gezelter 1116
68 chuckv 1266 soluteTypeLine = str()
69     solventTypeLine = str()
70     soluteMolLine = str()
71     nSolvents = 0
72    
73 kfletch2 1262 def usage():
74     print __doc__
75 gezelter 1116
76 kfletch2 1262 def readFile1(mdFileName):
77     mdFile = open(mdFileName, 'r')
78     # Find OOPSE version info first
79     line = mdFile.readline()
80     while 1:
81     if '<OOPSE version=' in line:
82     OOPSEversion = line
83     break
84     line = mdFile.readline()
85    
86     # Rewind file and find start of MetaData block
87    
88     mdFile.seek(0)
89     line = mdFile.readline()
90 chuckv 1266
91 kfletch2 1262 print "reading MetaData"
92     while 1:
93     if '<MetaData>' in line:
94     while 2:
95     metaData1.append(line)
96     line = mdFile.readline()
97 chuckv 1266 if 'type' in line:
98     global soluteTypeLine
99     soluteTypeLine = line
100     if 'nMol' in line:
101     global soluteMolLine
102     soluteMolLine = line
103 kfletch2 1262 if '</MetaData>' in line:
104     metaData1.append(line)
105     break
106     break
107     line = mdFile.readline()
108 gezelter 1116
109 kfletch2 1262 mdFile.seek(0)
110     print "reading Snapshot"
111     line = mdFile.readline()
112     while 1:
113     if '<Snapshot>' in line:
114     line = mdFile.readline()
115     while 1:
116     print "reading FrameData"
117     if '<FrameData>' in line:
118     while 2:
119     frameData1.append(line)
120     if 'Hmat:' in line:
121     L = line.split()
122     Hxx = float(L[2].strip(','))
123     Hxy = float(L[3].strip(','))
124     Hxz = float(L[4].strip(','))
125     Hyx = float(L[7].strip(','))
126     Hyy = float(L[8].strip(','))
127     Hyz = float(L[9].strip(','))
128     Hzx = float(L[12].strip(','))
129     Hzy = float(L[13].strip(','))
130     Hzz = float(L[14].strip(','))
131     Hmat1.append([Hxx, Hxy, Hxz])
132     Hmat1.append([Hyx, Hyy, Hyz])
133     Hmat1.append([Hzx, Hzy, Hzz])
134     print Hmat1
135     BoxInv1.append(1.0/Hxx)
136     BoxInv1.append(1.0/Hyy)
137     BoxInv1.append(1.0/Hzz)
138     print BoxInv1
139     line = mdFile.readline()
140     if '</FrameData>' in line:
141     frameData1.append(line)
142     break
143     break
144 gezelter 1116
145 kfletch2 1262 line = mdFile.readline()
146     while 1:
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 chuckv 1263 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 kfletch2 1262 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 gezelter 1116
185 kfletch2 1262 def readFile2(mdFileName):
186     mdFile = open(mdFileName, 'r')
187     # Find OOPSE version info first
188     line = mdFile.readline()
189     while 1:
190     if '<OOPSE version=' in line:
191     OOPSEversion = line
192     break
193     line = mdFile.readline()
194    
195     # Rewind file and find start of MetaData block
196    
197     mdFile.seek(0)
198     line = mdFile.readline()
199     print "reading MetaData"
200     while 1:
201 chuckv 1266 if '<MetaData>' in line:
202 kfletch2 1262 while 2:
203 chuckv 1266 if 'type' in line:
204     global solventTypeLine
205     solventTypeLine = line
206 kfletch2 1262 metaData2.append(line)
207     line = mdFile.readline()
208     if '</MetaData>' in line:
209     metaData2.append(line)
210     break
211     break
212     line = mdFile.readline()
213 gezelter 1116
214 kfletch2 1262 mdFile.seek(0)
215     print "reading Snapshot"
216     line = mdFile.readline()
217     while 1:
218     if '<Snapshot>' in line:
219     line = mdFile.readline()
220     while 1:
221     print "reading FrameData"
222     if '<FrameData>' in line:
223     while 2:
224     frameData2.append(line)
225     if 'Hmat:' in line:
226     L = line.split()
227     Hxx = float(L[2].strip(','))
228     Hxy = float(L[3].strip(','))
229     Hxz = float(L[4].strip(','))
230     Hyx = float(L[7].strip(','))
231     Hyy = float(L[8].strip(','))
232     Hyz = float(L[9].strip(','))
233     Hzx = float(L[12].strip(','))
234     Hzy = float(L[13].strip(','))
235     Hzz = float(L[14].strip(','))
236     Hmat2.append([Hxx, Hxy, Hxz])
237     Hmat2.append([Hyx, Hyy, Hyz])
238     Hmat2.append([Hzx, Hzy, Hzz])
239     print Hmat2
240     BoxInv2.append(1.0/Hxx)
241     BoxInv2.append(1.0/Hyy)
242     BoxInv2.append(1.0/Hzz)
243     print BoxInv2
244     line = mdFile.readline()
245     if '</FrameData>' in line:
246     frameData2.append(line)
247     break
248     break
249 gezelter 1116
250 kfletch2 1262 line = mdFile.readline()
251     while 1:
252     if '<StuntDoubles>' in line:
253     line = mdFile.readline()
254     while 2:
255     L = line.split()
256     myIndex = int(L[0])
257     indices2.append(myIndex)
258     pvqj2.append(L[1])
259     x = float(L[2])
260     y = float(L[3])
261     z = float(L[4])
262     positions2.append([x, y, z])
263     vx = float(L[5])
264     vy = float(L[6])
265     vz = float(L[7])
266     velocities2.append([vx, vy, vz])
267 chuckv 1263 if 'pvqj' in L[1]:
268     qw = float(L[8])
269     qx = float(L[9])
270     qy = float(L[10])
271     qz = float(L[11])
272     quaternions2.append([qw, qx, qy, qz])
273     jx = float(L[12])
274     jy = float(L[13])
275     jz = float(L[14])
276     angVels2.append([jx, jy, jz])
277     else:
278     quaternions1.append([0.0, 0.0, 0.0, 0.0])
279     angVels1.append([0.0, 0.0, 0.0])
280    
281 kfletch2 1262 line = mdFile.readline()
282     if '</StuntDoubles>' in line:
283     break
284     break
285     line = mdFile.readline()
286     if not line: break
287    
288     mdFile.close()
289 xsun 1214
290 kfletch2 1262 def writeFile(outputFileName):
291     outputFile = open(outputFileName, 'w')
292 gezelter 1116
293 chuckv 1266 outputFile.write("<OOPSE version=4>\n")
294 kfletch2 1262
295 chuckv 1266 # for metaline in metaData1:
296     # outputFile.write(metaline)
297     outputFile.write(" <MetaData>\n")
298     outputFile.write("\n\n")
299     outputFile.write("#include \"thiols.md\"\n")
300     outputFile.write("#include \"metals.md\"\n")
301     outputFile.write("\n\n")
302     outputFile.write("component{\n")
303     outputFile.write(soluteTypeLine)
304     outputFile.write(soluteMolLine)
305     outputFile.write("}\n")
306 gezelter 1116
307 chuckv 1266 outputFile.write("component{\n")
308     outputFile.write(solventTypeLine)
309     outputFile.write("nMol = %d;" % (nSolvents))
310     outputFile.write("}\n")
311     outputFile.write("\n\n")
312     outputFile.write("forceField = \"MnM\";\n\n")
313     outputFile.write("\n\n")
314     outputFile.write(" </MetaData>\n")
315 kfletch2 1262 outputFile.write(" <Snapshot>\n")
316 gezelter 1116
317 kfletch2 1262 for frameline in frameData1:
318     outputFile.write(frameline)
319    
320     outputFile.write(" <StuntDoubles>\n")
321 gezelter 1116
322 chuckv 1263
323 kfletch2 1262 newIndex = 0
324     for i in range(len(indices1)):
325 chuckv 1266 if (pvqj1[i] == 'pv'):
326 chuckv 1263 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]))
327 chuckv 1266 elif(pvqj1[i] == 'pvqj'):
328 chuckv 1263 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]))
329    
330 kfletch2 1262 newIndex = newIndex + 1
331 gezelter 1116
332 kfletch2 1262 outputFile.write(" </StuntDoubles>\n")
333     outputFile.write(" </Snapshot>\n")
334     outputFile.write("</OOPSE>\n")
335     outputFile.close()
336    
337     def roundMe(x):
338     if (x >= 0.0):
339     return math.floor(x + 0.5)
340     else:
341     return math.ceil(x - 0.5)
342 gezelter 1116
343 chuckv 1263 def frange(start,stop,step=1.0):
344     while start < stop:
345     yield start
346     start += step
347    
348    
349 kfletch2 1262 def wrapVector(myVect):
350     scaled = [0.0, 0.0, 0.0]
351     for i in range(3):
352     scaled[i] = myVect[i] * BoxInv1[i]
353     scaled[i] = scaled[i] - roundMe(scaled[i])
354     myVect[i] = scaled[i] * Hmat1[i][i]
355     return myVect
356 gezelter 1116
357 kfletch2 1262 def dot(L1, L2):
358     myDot = 0.0
359     for i in range(len(L1)):
360     myDot = myDot + L1[i]*L2[i]
361     return myDot
362 gezelter 1116
363 kfletch2 1262 def normalize(L1):
364     L2 = []
365     myLength = math.sqrt(dot(L1, L1))
366     for i in range(len(L1)):
367     L2.append(L1[i] / myLength)
368     return L2
369 gezelter 1116
370 kfletch2 1262 def cross(L1, L2):
371     # don't call this with anything other than length 3 lists please
372     # or you'll be sorry
373     L3 = [0.0, 0.0, 0.0]
374     L3[0] = L1[1]*L2[2] - L1[2]*L2[1]
375     L3[1] = L1[2]*L2[0] - L1[0]*L2[2]
376     L3[2] = L1[0]*L2[1] - L1[1]*L2[0]
377     return L3
378 gezelter 1116
379 chuckv 1265 def removeOverlaps(rcut,nSolventAtoms,nSoluteAtoms):
380 gezelter 1116
381 kfletch2 1262 rcut2 = rcut*rcut
382 chuckv 1263 nextMol = 0
383    
384 chuckv 1266
385 chuckv 1263 for i in range(0,len(indices2),nSolventAtoms):
386 kfletch2 1262 keepThisMolecule = 1
387 chuckv 1264 for atom1 in range (i, (i+nSolventAtoms)):
388 chuckv 1263 iPos = positions2[atom1]
389 kfletch2 1262
390 chuckv 1263 for j in range(0,len(indices1),nSoluteAtoms):
391 chuckv 1264 for atom2 in range (j, (j+nSoluteAtoms)):
392 chuckv 1263 jPos = positions1[j]
393     dpos = [jPos[0]-iPos[0], jPos[1]-iPos[1], jPos[2]-iPos[2]]
394     dpos = wrapVector(dpos)
395     dist2 = dot(dpos,dpos)
396 chuckv 1266
397 chuckv 1263 if (dist2 < rcut2):
398     keepThisMolecule = 0
399 chuckv 1266 # break
400    
401 kfletch2 1262 keepers.append(keepThisMolecule)
402 gezelter 1116
403 chuckv 1266 global nSolvents
404     myIndex = len(indices2) - 1
405     for i in range(0,len(keepers)):
406    
407 kfletch2 1262 if (keepers[i] == 1):
408 chuckv 1266 nSolvents = nSolvents + 1
409 chuckv 1263 for j in range (i, i+nSolventAtoms):
410 chuckv 1266 indices1.append(myIndex)
411     pvqj1.append(pvqj2[j])
412 chuckv 1265 if (pvqj2[j] == 'pv'):
413     positions1.append(positions2[j])
414     velocities1.append(velocities2[j])
415 chuckv 1266 quaternions1.append([0.0, 0.0, 0.0, 0.0])
416     angVels1.append([0.0, 0.0, 0.0])
417     else:
418     positions1.append(positions2[j])
419     velocities1.append(velocities2[j])
420 chuckv 1265 quaternions1.append(quaternions2[j])
421     angVels1.append(angVels2[j])
422 chuckv 1266 # indices1.append(indices2[j])
423     myIndex = myIndex +1
424    
425 kfletch2 1262 def main(argv):
426     try:
427 chuckv 1265 opts, args = getopt.getopt(argv, "hu:v:n:p:r:o:", ["help", "solute=","solvent=","nSoluteAtoms=","nSolventAtoms=", "rcut=" "output-file="])
428 kfletch2 1262 except getopt.GetoptError:
429     usage()
430     sys.exit(2)
431     for opt, arg in opts:
432     if opt in ("-h", "--help"):
433     usage()
434     sys.exit()
435     elif opt in ("-u", "--solute"):
436     mdFileName1 = arg
437     global _haveMDFileName1
438     _haveMDFileName1 = 1
439     elif opt in ("-v", "--solvent"):
440     mdFileName2 = arg
441     global _haveMDFileName2
442     _haveMDFileName2 = 1
443 chuckv 1263 elif opt in ("-n", "--nSoluteAtoms"):
444 chuckv 1265 nSoluteAtoms = int(arg)
445 chuckv 1263 global _haveNSoluteAtoms
446     _haveNSoluteAtoms = 1
447     elif opt in ("-p", "--nSolventAtoms"):
448 chuckv 1265 nSolventAtoms = int(arg)
449 chuckv 1263 global _haveNSolventAtoms
450     _haveNSolventAtoms = 1
451 kfletch2 1262 elif opt in ("-r", "--rcut"):
452     rcut = float(arg)
453     global _haveRcut
454     _haveRcut = 1
455     elif opt in ("-o", "--output-file"):
456     outputFileName = arg
457     global _haveOutputFileName
458     _haveOutputFileName = 1
459 gezelter 1116
460 kfletch2 1262 if (_haveMDFileName1 != 1):
461     usage()
462     print "No meta-data file was specified for the solute"
463     sys.exit()
464 gezelter 1116
465 kfletch2 1262 if (_haveMDFileName2 != 1):
466     usage()
467     print "No meta-data file was specified for the solvent"
468     sys.exit()
469 gezelter 1116
470 kfletch2 1262 if (_haveOutputFileName != 1):
471     usage()
472     print "No output file was specified"
473     sys.exit()
474 gezelter 1116
475 kfletch2 1262 if (_haveRcut != 1):
476     print "No cutoff radius was specified, using 4 angstroms"
477 chuckv 1266 rcut =4.0
478 gezelter 1116
479 chuckv 1263 if (_haveNSoluteAtoms != 1):
480     nSoluteAtoms = 1
481    
482     if (_haveNSolventAtoms != 1):
483     nSolventAtoms = 1
484    
485 kfletch2 1262 readFile1(mdFileName1)
486     readFile2(mdFileName2)
487 chuckv 1263 removeOverlaps(rcut,nSolventAtoms,nSoluteAtoms)
488 kfletch2 1262 writeFile(outputFileName)
489 gezelter 1116
490 kfletch2 1262 if __name__ == "__main__":
491     if len(sys.argv) == 1:
492     usage()
493     sys.exit()
494     main(sys.argv[1:])

Properties

Name Value
svn:executable *