ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/applications/utilities/md-solvator
Revision: 1262
Committed: Wed Jun 25 18:22:08 2008 UTC (17 years ago) by kfletch2
Original Path: trunk/src/applications/utilities/md-solvator
File size: 13121 byte(s)
Log Message:
*** empty log message ***

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

Properties

Name Value
svn:executable *