OpenMD 3.1
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
affineScale
1#!@Python3_EXECUTABLE@
2"""OpenMD affine scaling transform
3
4Takes an OpenMD file and scales both the periodic box and the
5coordinates of all StuntDoubles in the system by the same amount.
6
7You can either specify a new volume scaling for isotropic scaling, or
8specify one (or more) of the coordinates for non-isotropic scaling.
9
10Usage: affineScale
11
12Options:
13 -h, --help show this help
14 -m, --meta-data=... use specified meta-data (.omd, .eor) file
15 -o, --output-file=... use specified output file
16 -x, --newX=... scale the system to a new x dimension
17 -y, --newY=... scale the system to a new y dimension
18 -z, --newZ=... scale the system to a new z dimension
19 -v, --newV=... scale the system to a new volume
20
21Example:
22 affineScale -m lipidSystem.omd -o scaledSystem.omd -v 77000.0
23
24"""
25
26__author__ = "Dan Gezelter (gezelter@nd.edu)"
27__copyright__ = "Copyright (c) 2004-present The University of Notre Dame. All Rights Reserved."
28__license__ = "OpenMD"
29
30import sys
31import getopt
32import string
33import math
34import random
35
36metaData = []
37frameData = []
38indices = []
39pvqj = []
40p = []
41wrap = []
42v = []
43q = []
44j = []
45siteData = []
46Hmat = []
47BoxInv = []
48siteFlag = [False]
49def usage():
50 print(__doc__)
51
52def readFile(mdFileName):
53 mdFile = open(mdFileName, 'r')
54 # Find OpenMD version info first
55 line = mdFile.readline()
56 while True:
57 if '<OOPSE version=' in line or '<OpenMD version=' in line:
58 OpenMDversion = line
59 break
60 line = mdFile.readline()
61
62 # Rewind file and find start of MetaData block
63
64 mdFile.seek(0)
65 line = mdFile.readline()
66 print("Reading MetaData")
67 while True:
68 if '<MetaData>' in line:
69 while 2:
70 metaData.append(line)
71 line = mdFile.readline()
72 if '</MetaData>' in line:
73 metaData.append(line)
74 break
75 break
76 line = mdFile.readline()
77
78 line = mdFile.readline()
79 print("Reading Snapshot")
80 while True:
81 if '<Snapshot>' in line:
82 line = mdFile.readline()
83 while True:
84 print("Reading FrameData")
85 if '<FrameData>' in line:
86 while 2:
87 frameData.append(line)
88 if 'Hmat:' in line:
89 line = line.replace("{","").replace("}","").replace("Hmat:","")
90 L = line.split(",")
91 Hxx = float(L[0])
92 Hxy = float(L[1])
93 Hxz = float(L[2])
94 Hyx = float(L[3])
95 Hyy = float(L[4])
96 Hyz = float(L[5])
97 Hzx = float(L[6])
98 Hzy = float(L[7])
99 Hzz = float(L[8])
100 Hmat.append([Hxx, Hxy, Hxz])
101 Hmat.append([Hyx, Hyy, Hyz])
102 Hmat.append([Hzx, Hzy, Hzz])
103 BoxInv.append(1.0/Hxx)
104 BoxInv.append(1.0/Hyy)
105 BoxInv.append(1.0/Hzz)
106 line = mdFile.readline()
107 if '</FrameData>' in line:
108 frameData.append(line)
109 break
110 break
111
112 line = mdFile.readline()
113 while True:
114 if '<StuntDoubles>' in line:
115 line = mdFile.readline()
116 while 2:
117 L = line.split()
118 myIndex = int(L[0])
119 indices.append(myIndex)
120 pvqj.append(L[1])
121 i = 2
122 if 'p' in L[1]:
123 x = float(L[i])
124 y = float(L[i+1])
125 z = float(L[i+2])
126 p.append([x, y, z])
127 i = i+3
128 else:
129 p.append([0.0, 0.0, 0.0])
130 if 'v' in L[1]:
131 vx = float(L[i])
132 vy = float(L[i+1])
133 vz = float(L[i+2])
134 v.append([vx, vy, vz])
135 i = i+3
136 else:
137 v.append([0.0, 0.0, 0.0])
138 if 'q' in L[1]:
139 qw = float(L[i])
140 qx = float(L[i+1])
141 qy = float(L[i+2])
142 qz = float(L[i+3])
143 q.append([qw, qx, qy, qz])
144 i = i+4
145 else:
146 q.append([0.0, 0.0, 0.0, 0.0])
147 if 'j' in L[1]:
148 jx = float(L[i])
149 jy = float(L[i+1])
150 jz = float(L[i+2])
151 j.append([jx, jy, jz])
152 i = i+3
153 else:
154 j.append([0.0, 0.0, 0.0])
155 wrap.append([0, 0, 0])
156
157 line = mdFile.readline()
158 if '</StuntDoubles>' in line:
159 break
160 break
161 line = mdFile.readline()
162 while True:
163 if '<SiteData>' in line:
164 print("Reading SiteData")
165 siteFlag[0] = True
166 while True:
167 line = mdFile.readline()
168 siteData.append(line)
169 if '</SiteData>' in line:
170 break
171 break
172 break
173
174 line = mdFile.readline()
175 if not line: break
176 mdFile.close()
177
178def writeFile(outputFileName, repeatX, repeatY, repeatZ):
179 outputFile = open(outputFileName, 'w')
180
181 outputFile.write("<OpenMD version=2>\n");
182
183 print("Writing MetaData")
184 for metaline in metaData:
185 if 'nMol' in metaline:
186 metasplit = metaline.split()
187 nMol = float(metasplit[2].strip(';'))
188 newNmol = nMol * repeatX * repeatY * repeatZ
189 outputFile.write(' nMol = %10d;\n' % (newNmol))
190 else:
191 outputFile.write(metaline)
192
193 print("Writing Snapshot")
194 outputFile.write(" <Snapshot>\n")
195
196 print("Writing FrameData")
197 for frameline in frameData:
198 if 'Hmat:' in frameline:
199 myH = []
200 myH.append([repeatX * Hmat[0][0], repeatX * Hmat[0][1], repeatX * Hmat[0][2]])
201 myH.append([repeatY * Hmat[1][0], repeatY * Hmat[1][1], repeatY * Hmat[1][2]])
202 myH.append([repeatZ * Hmat[2][0], repeatZ * Hmat[2][1], repeatZ * Hmat[2][2]])
203 outputFile.write(" Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n" % (myH[0][0], myH[0][1], myH[0][2], myH[1][0], myH[1][1], myH[1][2], myH[2][0], myH[2][1], myH[2][2]))
204 else:
205 outputFile.write(frameline)
206
207 print("Writing StuntDoubles")
208 outputFile.write(" <StuntDoubles>\n")
209
210 #print(repeatX, repeatY, repeatZ)
211
212 deco = sorted([ (index, i) for i, index in enumerate(indices) ])
213 whichSD = 0
214 for index in range(len(deco)):
215 (index, i) = deco[index]
216 for ii in range(repeatX):
217 for jj in range(repeatY):
218 for kk in range(repeatZ):
219 myP = []
220 myP.append(p[i][0] + ii*Hmat[0][0] + jj*Hmat[1][0] + kk*Hmat[2][0])
221 myP.append(p[i][1] + ii*Hmat[0][1] + jj*Hmat[1][1] + kk*Hmat[2][1])
222 myP.append(p[i][2] + ii*Hmat[0][2] + jj*Hmat[1][2] + kk*Hmat[2][2])
223
224 if (pvqj[i] == 'p'):
225 outputFile.write("%10d %7s %18.10g %18.10g %18.10g\n" % (whichSD, pvqj[i], myP[0], myP[1], myP[2]))
226 elif (pvqj[i] == 'pv'):
227 outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e\n" % (whichSD, pvqj[i], myP[0], myP[1], myP[2], v[i][0], v[i][1], v[i][2]))
228 elif (pvqj[i] == 'pq'):
229 outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e %13e\n" % (whichSD, pvqj[i], myP[0], myP[1], myP[2], q[i][0], q[i][1], q[i][2], q[i][3]))
230 elif (pvqj[i] == 'pvqj'):
231 outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e %13e %13e %13e %13e %13e %13e %13e\n" % (whichSD, pvqj[i], myP[0], myP[1], myP[2], v[i][0], v[i][1], v[i][2], q[i][0], q[i][1], q[i][2], q[i][3], j[i][0], j[i][1], j[i][2]))
232 whichSD = whichSD + 1
233 outputFile.write(" </StuntDoubles>\n")
234 if (siteFlag[0]):
235 print("Writing SiteData")
236 outputFile.write(" <SiteData>\n")
237 for line in siteData:
238 outputFile.write(line)
239
240 print("Writing finished")
241 outputFile.write(" </Snapshot>\n")
242 outputFile.write("</OpenMD>\n")
243 outputFile.close()
244
245def roundMe(x):
246 if (x >= 0.0):
247 return math.floor(x + 0.5)
248 else:
249 return math.ceil(x - 0.5)
250
251def wrapVector(myVect):
252 scaled = [0.0, 0.0, 0.0]
253 wrappingNumber = [0, 0, 0]
254 for i in range(3):
255 scaled[i] = myVect[i] * BoxInv[i]
256 wrappingNumber[i] = roundMe(scaled[i])
257 scaled[i] = scaled[i] - wrappingNumber[i]
258 myVect[i] = scaled[i] * Hmat[i][i]
259 return myVect, wrappingNumber
260
261def dot(L1, L2):
262 myDot = 0.0
263 for i in range(len(L1)):
264 myDot = myDot + L1[i]*L2[i]
265 return myDot
266
267def normalize(L1):
268 L2 = []
269 myLength = math.sqrt(dot(L1, L1))
270 for i in range(len(L1)):
271 L2.append(L1[i] / myLength)
272 return L2
273
274def cross(L1, L2):
275 # don't call this with anything other than length 3 lists please
276 # or you'll be sorry
277 L3 = [0.0, 0.0, 0.0]
278 L3[0] = L1[1]*L2[2] - L1[2]*L2[1]
279 L3[1] = L1[2]*L2[0] - L1[0]*L2[2]
280 L3[2] = L1[0]*L2[1] - L1[1]*L2[0]
281 return L3
282
283def mapToBox():
284 for i in range(len(indices)):
285 dpos = []
286 dpos.append(p[i][0])
287 dpos.append(p[i][1])
288 dpos.append(p[i][2])
289 p[i], wrap[i] = wrapVector(dpos)
290
291def scaleBox(scaleX, scaleY, scaleZ):
292 for i in range(3):
293 Hmat[0][i] = scaleX * Hmat[0][i]
294 for i in range(3):
295 Hmat[1][i] = scaleY * Hmat[1][i]
296 for i in range(3):
297 Hmat[2][i] = scaleZ * Hmat[2][i]
298 for i in range(3):
299 BoxInv[i] = 1.0/Hmat[i][i]
300
301def scaleCoordinates(scaleX, scaleY, scaleZ):
302 for i in range(len(indices)):
303 p[i][0] = p[i][0]*scaleX + wrap[i][0] * Hmat[0][0]
304 p[i][1] = p[i][1]*scaleY + wrap[i][1] * Hmat[1][1]
305 p[i][2] = p[i][2]*scaleZ + wrap[i][2] * Hmat[2][2]
306
307def main(argv):
308 _haveMDFileName = 0
309 _haveOutputFileName = 0
310 try:
311 opts, args = getopt.getopt(argv, "hm:o:x:y:z:v:", ["help", "meta-data=", "output-file=", "newX=", "newY=", "newZ=", "newV="])
312 except getopt.GetoptError:
313 usage()
314 sys.exit(2)
315 doV = 0
316 doX = 0
317 doY = 0
318 doZ = 0
319 for opt, arg in opts:
320 if opt in ("-h", "--help"):
321 usage()
322 sys.exit()
323 elif opt in ("-m", "--meta-data"):
324 mdFileName = arg
325 _haveMDFileName = 1
326 elif opt in ("-o", "--output-file"):
327 outputFileName = arg
328 _haveOutputFileName = 1
329 if opt in ("-x", "--newX"):
330 newX = float(arg)
331 doX = 1
332 elif opt in ("-y", "--newY"):
333 newY = float(arg)
334 doY = 1
335 elif opt in ("-z", "--newZ"):
336 newZ = float(arg)
337 doZ = 1
338 elif opt in ("-v", "--newV"):
339 newV = float(arg)
340 doV = 1
341
342 if (_haveMDFileName != 1):
343 usage()
344 print("No meta-data file was specified")
345 sys.exit()
346 if (_haveOutputFileName != 1):
347 usage()
348 print("No output file was specified")
349 sys.exit()
350
351 if not (doV or doX or doY or doZ):
352 usage()
353 print("no scaling options given. Nothing to do!")
354 sys.exit()
355
356 if doV and (doX or doY or doZ):
357 usage()
358 print("-v is mutually exclusive with any of the -x, -y, and -z options")
359 sys.exit()
360
361 readFile(mdFileName)
362
363 scaleX = 1.0
364 scaleY = 1.0
365 scaleZ = 1.0
366
367 if doX:
368 scaleX = newX / Hmat[0][0]
369 if doY:
370 scaleY = newY / Hmat[1][1]
371 if doZ:
372 scaleZ = newZ / Hmat[2][2]
373
374 if doV:
375 oldV = Hmat[0][0] * Hmat[1][1] * Hmat[2][2]
376 scaleX = pow(newV/oldV, 1.0/3.0)
377 scaleY = scaleX
378 scaleZ = scaleX
379
380 mapToBox()
381 scaleBox(scaleX, scaleY, scaleZ)
382 scaleCoordinates(scaleX, scaleY, scaleZ)
383 writeFile(outputFileName, 1, 1, 1)
384
385if __name__ == "__main__":
386 if len(sys.argv) == 1:
387 usage()
388 sys.exit()
389 main(sys.argv[1:])