ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/applications/utilities/md-solvator
Revision: 1263
Committed: Thu Jun 26 14:23:49 2008 UTC (17 years ago) by chuckv
Original Path: trunk/src/applications/utilities/md-solvator
File size: 14821 byte(s)
Log Message:
More work on python version of md-solvator. Pretty sure mol loop is busted.

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 1263 __version__ = "$Revision: 1.5 $"
25     __date__ = "$Date: 2008-06-26 14:23:49 $"
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 gezelter 1116
42 kfletch2 1262 metaData1 = []
43     frameData1 = []
44     positions1 = []
45     velocities1 = []
46     quaternions1 = []
47     angVels1 = []
48     indices1 = []
49     Hmat1 = []
50     BoxInv1 = []
51 chuckv 1263 pvqj1 = []
52 gezelter 1116
53 kfletch2 1262 metaData2 = []
54     frameData2 = []
55     positions2 = []
56     velocities2 = []
57     quaternions2 = []
58     angVels2 = []
59     indices2 = []
60     Hmat2 = []
61     BoxInv2 = []
62 chuckv 1263 pvqj2 = []
63 gezelter 1116
64 kfletch2 1262 keepers = []
65 gezelter 1116
66 kfletch2 1262 def usage():
67     print __doc__
68 gezelter 1116
69 kfletch2 1262 def readFile1(mdFileName):
70     mdFile = open(mdFileName, 'r')
71     # Find OOPSE version info first
72     line = mdFile.readline()
73     while 1:
74     if '<OOPSE version=' in line:
75     OOPSEversion = line
76     break
77     line = mdFile.readline()
78    
79     # Rewind file and find start of MetaData block
80    
81     mdFile.seek(0)
82     line = mdFile.readline()
83     print "reading MetaData"
84     while 1:
85     if '<MetaData>' in line:
86     while 2:
87     metaData1.append(line)
88     line = mdFile.readline()
89     if '</MetaData>' in line:
90     metaData1.append(line)
91     break
92     break
93     line = mdFile.readline()
94 gezelter 1116
95 kfletch2 1262 mdFile.seek(0)
96     print "reading Snapshot"
97     line = mdFile.readline()
98     while 1:
99     if '<Snapshot>' in line:
100     line = mdFile.readline()
101     while 1:
102     print "reading FrameData"
103     if '<FrameData>' in line:
104     while 2:
105     frameData1.append(line)
106     if 'Hmat:' in line:
107     L = line.split()
108     Hxx = float(L[2].strip(','))
109     Hxy = float(L[3].strip(','))
110     Hxz = float(L[4].strip(','))
111     Hyx = float(L[7].strip(','))
112     Hyy = float(L[8].strip(','))
113     Hyz = float(L[9].strip(','))
114     Hzx = float(L[12].strip(','))
115     Hzy = float(L[13].strip(','))
116     Hzz = float(L[14].strip(','))
117     Hmat1.append([Hxx, Hxy, Hxz])
118     Hmat1.append([Hyx, Hyy, Hyz])
119     Hmat1.append([Hzx, Hzy, Hzz])
120     print Hmat1
121     BoxInv1.append(1.0/Hxx)
122     BoxInv1.append(1.0/Hyy)
123     BoxInv1.append(1.0/Hzz)
124     print BoxInv1
125     line = mdFile.readline()
126     if '</FrameData>' in line:
127     frameData1.append(line)
128     break
129     break
130 gezelter 1116
131 kfletch2 1262 line = mdFile.readline()
132     while 1:
133     if '<StuntDoubles>' in line:
134     line = mdFile.readline()
135     while 2:
136     L = line.split()
137     myIndex = int(L[0])
138     indices1.append(myIndex)
139     pvqj1.append(L[1])
140     x = float(L[2])
141     y = float(L[3])
142     z = float(L[4])
143     positions1.append([x, y, z])
144     vx = float(L[5])
145     vy = float(L[6])
146     vz = float(L[7])
147     velocities1.append([vx, vy, vz])
148 chuckv 1263 if 'pvqj' in L[1]:
149     qw = float(L[8])
150     qx = float(L[9])
151     qy = float(L[10])
152     qz = float(L[11])
153     quaternions1.append([qw, qx, qy, qz])
154     jx = float(L[12])
155     jy = float(L[13])
156     jz = float(L[14])
157     angVels1.append([jx, jy, jz])
158     else:
159     quaternions1.append([0.0, 0.0, 0.0, 0.0])
160     angVels1.append([0.0, 0.0, 0.0])
161    
162 kfletch2 1262 line = mdFile.readline()
163     if '</StuntDoubles>' in line:
164     break
165     break
166     line = mdFile.readline()
167     if not line: break
168    
169     mdFile.close()
170 gezelter 1116
171 kfletch2 1262 def readFile2(mdFileName):
172     mdFile = open(mdFileName, 'r')
173     # Find OOPSE version info first
174     line = mdFile.readline()
175     while 1:
176     if '<OOPSE version=' in line:
177     OOPSEversion = line
178     break
179     line = mdFile.readline()
180    
181     # Rewind file and find start of MetaData block
182    
183     mdFile.seek(0)
184     line = mdFile.readline()
185     print "reading MetaData"
186     while 1:
187     if '<MetaData>' in line:
188     while 2:
189     metaData2.append(line)
190     line = mdFile.readline()
191     if '</MetaData>' in line:
192     metaData2.append(line)
193     break
194     break
195     line = mdFile.readline()
196 gezelter 1116
197 kfletch2 1262 mdFile.seek(0)
198     print "reading Snapshot"
199     line = mdFile.readline()
200     while 1:
201     if '<Snapshot>' in line:
202     line = mdFile.readline()
203     while 1:
204     print "reading FrameData"
205     if '<FrameData>' in line:
206     while 2:
207     frameData2.append(line)
208     if 'Hmat:' in line:
209     L = line.split()
210     Hxx = float(L[2].strip(','))
211     Hxy = float(L[3].strip(','))
212     Hxz = float(L[4].strip(','))
213     Hyx = float(L[7].strip(','))
214     Hyy = float(L[8].strip(','))
215     Hyz = float(L[9].strip(','))
216     Hzx = float(L[12].strip(','))
217     Hzy = float(L[13].strip(','))
218     Hzz = float(L[14].strip(','))
219     Hmat2.append([Hxx, Hxy, Hxz])
220     Hmat2.append([Hyx, Hyy, Hyz])
221     Hmat2.append([Hzx, Hzy, Hzz])
222     print Hmat2
223     BoxInv2.append(1.0/Hxx)
224     BoxInv2.append(1.0/Hyy)
225     BoxInv2.append(1.0/Hzz)
226     print BoxInv2
227     line = mdFile.readline()
228     if '</FrameData>' in line:
229     frameData2.append(line)
230     break
231     break
232 gezelter 1116
233 kfletch2 1262 line = mdFile.readline()
234     while 1:
235     if '<StuntDoubles>' in line:
236     line = mdFile.readline()
237     while 2:
238     L = line.split()
239     myIndex = int(L[0])
240     indices2.append(myIndex)
241     pvqj2.append(L[1])
242     x = float(L[2])
243     y = float(L[3])
244     z = float(L[4])
245     positions2.append([x, y, z])
246     vx = float(L[5])
247     vy = float(L[6])
248     vz = float(L[7])
249     velocities2.append([vx, vy, vz])
250 chuckv 1263 if 'pvqj' in L[1]:
251     qw = float(L[8])
252     qx = float(L[9])
253     qy = float(L[10])
254     qz = float(L[11])
255     quaternions2.append([qw, qx, qy, qz])
256     jx = float(L[12])
257     jy = float(L[13])
258     jz = float(L[14])
259     angVels2.append([jx, jy, jz])
260     else:
261     quaternions1.append([0.0, 0.0, 0.0, 0.0])
262     angVels1.append([0.0, 0.0, 0.0])
263    
264 kfletch2 1262 line = mdFile.readline()
265     if '</StuntDoubles>' in line:
266     break
267     break
268     line = mdFile.readline()
269     if not line: break
270    
271     mdFile.close()
272 xsun 1214
273 kfletch2 1262 def writeFile(outputFileName):
274     outputFile = open(outputFileName, 'w')
275 gezelter 1116
276 kfletch2 1262 outputFile.write("<OOPSE version=4>\n");
277    
278     for metaline in metaData1:
279     outputFile.write(metaline)
280 gezelter 1116
281 kfletch2 1262 outputFile.write(" <Snapshot>\n")
282 gezelter 1116
283 kfletch2 1262 for frameline in frameData1:
284     outputFile.write(frameline)
285    
286     outputFile.write(" <StuntDoubles>\n")
287 gezelter 1116
288 chuckv 1263
289 kfletch2 1262 newIndex = 0
290     for i in range(len(indices1)):
291 chuckv 1263 if (pvqj[i] == 'pv'):
292     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]))
293     elif(pvqj[i] == 'pvqj'):
294     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]))
295    
296 kfletch2 1262 newIndex = newIndex + 1
297 gezelter 1116
298 kfletch2 1262 outputFile.write(" </StuntDoubles>\n")
299     outputFile.write(" </Snapshot>\n")
300     outputFile.write("</OOPSE>\n")
301     outputFile.close()
302    
303     def roundMe(x):
304     if (x >= 0.0):
305     return math.floor(x + 0.5)
306     else:
307     return math.ceil(x - 0.5)
308 gezelter 1116
309 chuckv 1263 def frange(start,stop,step=1.0):
310     while start < stop:
311     yield start
312     start += step
313    
314    
315 kfletch2 1262 def wrapVector(myVect):
316     scaled = [0.0, 0.0, 0.0]
317     for i in range(3):
318     scaled[i] = myVect[i] * BoxInv1[i]
319     scaled[i] = scaled[i] - roundMe(scaled[i])
320     myVect[i] = scaled[i] * Hmat1[i][i]
321     return myVect
322 gezelter 1116
323 kfletch2 1262 def dot(L1, L2):
324     myDot = 0.0
325     for i in range(len(L1)):
326     myDot = myDot + L1[i]*L2[i]
327     return myDot
328 gezelter 1116
329 kfletch2 1262 def normalize(L1):
330     L2 = []
331     myLength = math.sqrt(dot(L1, L1))
332     for i in range(len(L1)):
333     L2.append(L1[i] / myLength)
334     return L2
335 gezelter 1116
336 kfletch2 1262 def cross(L1, L2):
337     # don't call this with anything other than length 3 lists please
338     # or you'll be sorry
339     L3 = [0.0, 0.0, 0.0]
340     L3[0] = L1[1]*L2[2] - L1[2]*L2[1]
341     L3[1] = L1[2]*L2[0] - L1[0]*L2[2]
342     L3[2] = L1[0]*L2[1] - L1[1]*L2[0]
343     return L3
344 gezelter 1116
345 chuckv 1263 def removeOverlaps(rcut,nSolventAtoms,nSoluteAtoms:
346 gezelter 1116
347 kfletch2 1262 rcut2 = rcut*rcut
348 chuckv 1263 nextMol = 0
349    
350    
351     for i in range(0,len(indices2),nSolventAtoms):
352 kfletch2 1262 keepThisMolecule = 1
353 chuckv 1263 for atom1 in range (i, (i+nSolventAtoms-1)):
354     iPos = positions2[atom1]
355 kfletch2 1262
356 chuckv 1263 for j in range(0,len(indices1),nSoluteAtoms):
357     for atom2 in range (j, (j+nSoluteAtoms-1)):
358     jPos = positions1[j]
359     dpos = [jPos[0]-iPos[0], jPos[1]-iPos[1], jPos[2]-iPos[2]]
360     dpos = wrapVector(dpos)
361     dist2 = dot(dpos,dpos)
362 gezelter 1116
363 chuckv 1263 if (dist2 < rcut2):
364     keepThisMolecule = 0
365     break
366    
367 kfletch2 1262 keepers.append(keepThisMolecule)
368 gezelter 1116
369 chuckv 1263 for i in range(0,len(indices2),nSolventAtoms):
370 kfletch2 1262 if (keepers[i] == 1):
371 chuckv 1263 for j in range (i, i+nSolventAtoms):
372     positions1.append(positions2[j])
373     velocities1.append(velocities2[j])
374     quaternions1.append(quaternions2[j])
375     angVels1.append(angVels2[j])
376     indices1.append(indices2[j])
377 gezelter 1116
378 kfletch2 1262 def main(argv):
379     try:
380 chuckv 1263 opts, args = getopt.getopt(argv, "hu:v:n:p:r:o:", ["help", "solute=","solvent=","nSoluteAtoms","nSoluteAtoms", "rcut" "output-file="])
381 kfletch2 1262 except getopt.GetoptError:
382     usage()
383     sys.exit(2)
384     for opt, arg in opts:
385     if opt in ("-h", "--help"):
386     usage()
387     sys.exit()
388     elif opt in ("-u", "--solute"):
389     mdFileName1 = arg
390     global _haveMDFileName1
391     _haveMDFileName1 = 1
392     elif opt in ("-v", "--solvent"):
393     mdFileName2 = arg
394     global _haveMDFileName2
395     _haveMDFileName2 = 1
396 chuckv 1263 elif opt in ("-n", "--nSoluteAtoms"):
397     nSoluteAtoms = arg
398     global _haveNSoluteAtoms
399     _haveNSoluteAtoms = 1
400     elif opt in ("-p", "--nSolventAtoms"):
401     mdFileName2 = arg
402     global _haveNSolventAtoms
403     _haveNSolventAtoms = 1
404 kfletch2 1262 elif opt in ("-r", "--rcut"):
405     rcut = float(arg)
406     global _haveRcut
407     _haveRcut = 1
408     elif opt in ("-o", "--output-file"):
409     outputFileName = arg
410     global _haveOutputFileName
411     _haveOutputFileName = 1
412 gezelter 1116
413 kfletch2 1262 if (_haveMDFileName1 != 1):
414     usage()
415     print "No meta-data file was specified for the solute"
416     sys.exit()
417 gezelter 1116
418 kfletch2 1262 if (_haveMDFileName2 != 1):
419     usage()
420     print "No meta-data file was specified for the solvent"
421     sys.exit()
422 gezelter 1116
423 kfletch2 1262 if (_haveOutputFileName != 1):
424     usage()
425     print "No output file was specified"
426     sys.exit()
427 gezelter 1116
428 kfletch2 1262 if (_haveRcut != 1):
429     print "No cutoff radius was specified, using 4 angstroms"
430     rcut = 4.0
431 gezelter 1116
432 chuckv 1263 if (_haveNSoluteAtoms != 1):
433     nSoluteAtoms = 1
434    
435     if (_haveNSolventAtoms != 1):
436     nSolventAtoms = 1
437    
438 kfletch2 1262 readFile1(mdFileName1)
439     readFile2(mdFileName2)
440 chuckv 1263 removeOverlaps(rcut,nSolventAtoms,nSoluteAtoms)
441 kfletch2 1262 writeFile(outputFileName)
442 gezelter 1116
443 kfletch2 1262 if __name__ == "__main__":
444     if len(sys.argv) == 1:
445     usage()
446     sys.exit()
447     main(sys.argv[1:])

Properties

Name Value
svn:executable *