1 |
gezelter |
3345 |
#!/usr/bin/env python |
2 |
|
|
"""MetaData file remapper |
3 |
|
|
|
4 |
|
|
Takes a MetaData file and maps all StuntDoubles back to the periodic box. |
5 |
|
|
Will optionally replicate the system in the three box directions, and |
6 |
|
|
then writes out a new MetaData file. |
7 |
|
|
|
8 |
|
|
Usage: md2md |
9 |
|
|
|
10 |
|
|
Options: |
11 |
|
|
-h, --help show this help |
12 |
|
|
-m, --meta-data=... use specified meta-data (.md, .eor) file |
13 |
|
|
-o, --output-file=... use specified output file |
14 |
|
|
-x, --repeatX=... make the system repeat in the x direction |
15 |
|
|
-y, --repeatY=... make the system repeat in the y direction |
16 |
|
|
-z, --repeatZ=... make the system repeat in the z direction |
17 |
|
|
|
18 |
|
|
|
19 |
|
|
Example: |
20 |
|
|
md2md -m lipidSystem.md -o bigLipidSystem.md -x 2 -y 2 -z 1 |
21 |
|
|
|
22 |
|
|
""" |
23 |
|
|
|
24 |
|
|
__author__ = "Dan Gezelter (gezelter@nd.edu)" |
25 |
xsun |
3393 |
__version__ = "$Revision: 1.2 $" |
26 |
|
|
__date__ = "$Date: 2008-05-14 20:49:17 $" |
27 |
gezelter |
3345 |
__copyright__ = "Copyright (c) 2008 by the University of Notre Dame" |
28 |
|
|
__license__ = "OOPSE" |
29 |
|
|
|
30 |
|
|
import sys |
31 |
|
|
import getopt |
32 |
|
|
import string |
33 |
|
|
import math |
34 |
|
|
import random |
35 |
|
|
from sets import * |
36 |
|
|
|
37 |
|
|
_haveMDFileName = 0 |
38 |
|
|
_haveOutputFileName = 0 |
39 |
|
|
|
40 |
|
|
repeatX = 1 |
41 |
|
|
repeatY = 1 |
42 |
|
|
repeatZ = 1 |
43 |
|
|
|
44 |
|
|
metaData = [] |
45 |
|
|
frameData = [] |
46 |
|
|
indices = [] |
47 |
|
|
pvqj = [] |
48 |
|
|
p = [] |
49 |
|
|
v = [] |
50 |
|
|
q = [] |
51 |
|
|
j = [] |
52 |
|
|
Hmat = [] |
53 |
|
|
BoxInv = [] |
54 |
|
|
|
55 |
|
|
def usage(): |
56 |
|
|
print __doc__ |
57 |
|
|
|
58 |
|
|
def readFile(mdFileName): |
59 |
|
|
mdFile = open(mdFileName, 'r') |
60 |
|
|
# Find OOPSE version info first |
61 |
|
|
line = mdFile.readline() |
62 |
|
|
while 1: |
63 |
|
|
if '<OOPSE version=' in line: |
64 |
|
|
OOPSEversion = line |
65 |
|
|
break |
66 |
|
|
line = mdFile.readline() |
67 |
|
|
|
68 |
|
|
# Rewind file and find start of MetaData block |
69 |
|
|
|
70 |
|
|
mdFile.seek(0) |
71 |
|
|
line = mdFile.readline() |
72 |
|
|
print "reading MetaData" |
73 |
|
|
while 1: |
74 |
|
|
if '<MetaData>' in line: |
75 |
|
|
while 2: |
76 |
|
|
metaData.append(line) |
77 |
|
|
line = mdFile.readline() |
78 |
|
|
if '</MetaData>' in line: |
79 |
|
|
metaData.append(line) |
80 |
|
|
break |
81 |
|
|
break |
82 |
|
|
line = mdFile.readline() |
83 |
|
|
|
84 |
|
|
mdFile.seek(0) |
85 |
|
|
print "reading Snapshot" |
86 |
|
|
line = mdFile.readline() |
87 |
|
|
while 1: |
88 |
|
|
if '<Snapshot>' in line: |
89 |
|
|
line = mdFile.readline() |
90 |
|
|
while 1: |
91 |
|
|
print "reading FrameData" |
92 |
|
|
if '<FrameData>' in line: |
93 |
|
|
while 2: |
94 |
|
|
frameData.append(line) |
95 |
|
|
if 'Hmat:' in line: |
96 |
|
|
L = line.split() |
97 |
|
|
Hxx = float(L[2].strip(',')) |
98 |
|
|
Hxy = float(L[3].strip(',')) |
99 |
|
|
Hxz = float(L[4].strip(',')) |
100 |
|
|
Hyx = float(L[7].strip(',')) |
101 |
|
|
Hyy = float(L[8].strip(',')) |
102 |
|
|
Hyz = float(L[9].strip(',')) |
103 |
|
|
Hzx = float(L[12].strip(',')) |
104 |
|
|
Hzy = float(L[13].strip(',')) |
105 |
|
|
Hzz = float(L[14].strip(',')) |
106 |
|
|
Hmat.append([Hxx, Hxy, Hxz]) |
107 |
|
|
Hmat.append([Hyx, Hyy, Hyz]) |
108 |
|
|
Hmat.append([Hzx, Hzy, Hzz]) |
109 |
|
|
BoxInv.append(1.0/Hxx) |
110 |
|
|
BoxInv.append(1.0/Hyy) |
111 |
|
|
BoxInv.append(1.0/Hzz) |
112 |
|
|
line = mdFile.readline() |
113 |
|
|
if '</FrameData>' in line: |
114 |
|
|
frameData.append(line) |
115 |
|
|
break |
116 |
|
|
break |
117 |
|
|
|
118 |
|
|
line = mdFile.readline() |
119 |
|
|
while 1: |
120 |
|
|
if '<StuntDoubles>' in line: |
121 |
|
|
line = mdFile.readline() |
122 |
|
|
while 2: |
123 |
|
|
L = line.split() |
124 |
|
|
myIndex = int(L[0]) |
125 |
|
|
indices.append(myIndex) |
126 |
|
|
pvqj.append(L[1]) |
127 |
|
|
x = float(L[2]) |
128 |
|
|
y = float(L[3]) |
129 |
|
|
z = float(L[4]) |
130 |
|
|
p.append([x, y, z]) |
131 |
|
|
vx = float(L[5]) |
132 |
|
|
vy = float(L[6]) |
133 |
|
|
vz = float(L[7]) |
134 |
|
|
v.append([vx, vy, vz]) |
135 |
|
|
if 'pvqj' in L[1]: |
136 |
|
|
qw = float(L[8]) |
137 |
|
|
qx = float(L[9]) |
138 |
|
|
qy = float(L[10]) |
139 |
|
|
qz = float(L[11]) |
140 |
|
|
q.append([qw, qx, qy, qz]) |
141 |
|
|
jx = float(L[12]) |
142 |
|
|
jy = float(L[13]) |
143 |
|
|
jz = float(L[14]) |
144 |
|
|
j.append([jx, jy, jz]) |
145 |
|
|
else: |
146 |
|
|
q.append([0.0, 0.0, 0.0, 0.0]) |
147 |
|
|
j.append([0.0, 0.0, 0.0]) |
148 |
|
|
|
149 |
|
|
line = mdFile.readline() |
150 |
|
|
if '</StuntDoubles>' in line: |
151 |
|
|
break |
152 |
|
|
break |
153 |
|
|
line = mdFile.readline() |
154 |
|
|
if not line: break |
155 |
|
|
|
156 |
|
|
mdFile.close() |
157 |
|
|
|
158 |
|
|
def writeFile(outputFileName): |
159 |
|
|
outputFile = open(outputFileName, 'w') |
160 |
|
|
|
161 |
|
|
outputFile.write("<OOPSE version=4>\n"); |
162 |
|
|
|
163 |
|
|
print "writing MetaData" |
164 |
|
|
for metaline in metaData: |
165 |
|
|
if 'nMol' in metaline: |
166 |
|
|
metasplit = metaline.split() |
167 |
|
|
nMol = float(metasplit[2].strip(';')) |
168 |
|
|
newNmol = nMol * repeatX * repeatY * repeatZ |
169 |
|
|
outputFile.write(' nMol = %10d;\n' % (newNmol)) |
170 |
|
|
else: |
171 |
|
|
outputFile.write(metaline) |
172 |
|
|
|
173 |
|
|
print "writing Snapshot" |
174 |
|
|
outputFile.write(" <Snapshot>\n") |
175 |
|
|
|
176 |
|
|
print "writing FrameData" |
177 |
|
|
for frameline in frameData: |
178 |
|
|
if 'Hmat:' in frameline: |
179 |
|
|
myH = [] |
180 |
|
|
myH.append([repeatX * Hmat[0][0], repeatX * Hmat[0][1], repeatX * Hmat[0][2]]) |
181 |
|
|
myH.append([repeatY * Hmat[1][0], repeatY * Hmat[1][1], repeatY * Hmat[1][2]]) |
182 |
|
|
myH.append([repeatZ * Hmat[2][0], repeatZ * Hmat[2][1], repeatZ * Hmat[2][2]]) |
183 |
|
|
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])) |
184 |
|
|
else: |
185 |
|
|
outputFile.write(frameline) |
186 |
|
|
|
187 |
|
|
print "writing StuntDoubles" |
188 |
|
|
outputFile.write(" <StuntDoubles>\n") |
189 |
|
|
whichSD = 0 |
190 |
|
|
|
191 |
xsun |
3393 |
print repeatX, repeatY, repeatZ |
192 |
|
|
|
193 |
gezelter |
3345 |
for i in range(len(indices)): |
194 |
|
|
for ii in range(repeatX): |
195 |
|
|
for jj in range(repeatY): |
196 |
|
|
for kk in range(repeatZ): |
197 |
xsun |
3393 |
|
198 |
|
|
myP = [] |
199 |
|
|
myP.append(p[i][0] + ii*Hmat[0][0] + jj*Hmat[1][0] + kk*Hmat[2][0]) |
200 |
|
|
myP.append(p[i][1] + ii*Hmat[0][1] + jj*Hmat[1][1] + kk*Hmat[2][1]) |
201 |
|
|
myP.append(p[i][2] + ii*Hmat[0][2] + jj*Hmat[1][2] + kk*Hmat[2][2]) |
202 |
gezelter |
3345 |
|
203 |
|
|
if (pvqj[i] == 'pv'): |
204 |
|
|
outputFile.write("%10d %7s %18.10g %18.10g %18.10g %14e %13e %13e\n" % (whichSD, pvqj[i], myP[0], myP[1], myP[2], v[i][0], v[i][1], v[i][2])) |
205 |
xsun |
3393 |
elif (pvqj[i] == 'pvqj'): |
206 |
gezelter |
3345 |
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])) |
207 |
|
|
|
208 |
|
|
whichSD = whichSD + 1 |
209 |
|
|
|
210 |
|
|
outputFile.write(" </StuntDoubles>\n") |
211 |
|
|
outputFile.write(" </Snapshot>\n") |
212 |
|
|
outputFile.write("</OOPSE>\n") |
213 |
|
|
outputFile.close() |
214 |
|
|
|
215 |
|
|
def roundMe(x): |
216 |
|
|
if (x >= 0.0): |
217 |
|
|
return math.floor(x + 0.5) |
218 |
|
|
else: |
219 |
|
|
return math.ceil(x - 0.5) |
220 |
|
|
|
221 |
|
|
def wrapVector(myVect): |
222 |
|
|
scaled = [0.0, 0.0, 0.0] |
223 |
|
|
for i in range(3): |
224 |
|
|
scaled[i] = myVect[i] * BoxInv[i] |
225 |
|
|
scaled[i] = scaled[i] - roundMe(scaled[i]) |
226 |
|
|
myVect[i] = scaled[i] * Hmat[i][i] |
227 |
|
|
return myVect |
228 |
|
|
|
229 |
|
|
def dot(L1, L2): |
230 |
|
|
myDot = 0.0 |
231 |
|
|
for i in range(len(L1)): |
232 |
|
|
myDot = myDot + L1[i]*L2[i] |
233 |
|
|
return myDot |
234 |
|
|
|
235 |
|
|
def normalize(L1): |
236 |
|
|
L2 = [] |
237 |
|
|
myLength = math.sqrt(dot(L1, L1)) |
238 |
|
|
for i in range(len(L1)): |
239 |
|
|
L2.append(L1[i] / myLength) |
240 |
|
|
return L2 |
241 |
|
|
|
242 |
|
|
def cross(L1, L2): |
243 |
|
|
# don't call this with anything other than length 3 lists please |
244 |
|
|
# or you'll be sorry |
245 |
|
|
L3 = [0.0, 0.0, 0.0] |
246 |
|
|
L3[0] = L1[1]*L2[2] - L1[2]*L2[1] |
247 |
|
|
L3[1] = L1[2]*L2[0] - L1[0]*L2[2] |
248 |
|
|
L3[2] = L1[0]*L2[1] - L1[1]*L2[0] |
249 |
|
|
return L3 |
250 |
|
|
|
251 |
|
|
def mapToBox(): |
252 |
|
|
for i in range(len(indices)): |
253 |
|
|
dpos = p[i] |
254 |
|
|
p[i] = wrapVector(dpos) |
255 |
|
|
|
256 |
|
|
|
257 |
|
|
def main(argv): |
258 |
|
|
try: |
259 |
|
|
opts, args = getopt.getopt(argv, "hm:o:x:y:z:", ["help", "meta-data=", "output-file=", "repeatX=", "repeatY=", "repeatZ="]) |
260 |
|
|
except getopt.GetoptError: |
261 |
|
|
usage() |
262 |
|
|
sys.exit(2) |
263 |
|
|
for opt, arg in opts: |
264 |
|
|
if opt in ("-h", "--help"): |
265 |
|
|
usage() |
266 |
|
|
sys.exit() |
267 |
|
|
elif opt in ("-m", "--meta-data"): |
268 |
|
|
mdFileName = arg |
269 |
|
|
global _haveMDFileName |
270 |
|
|
_haveMDFileName = 1 |
271 |
|
|
elif opt in ("-o", "--output-file"): |
272 |
|
|
outputFileName = arg |
273 |
|
|
global _haveOutputFileName |
274 |
|
|
_haveOutputFileName = 1 |
275 |
xsun |
3393 |
elif opt in ("-x", "--repeatX"): |
276 |
|
|
global repeatX |
277 |
gezelter |
3345 |
repeatX = int(arg) |
278 |
|
|
elif opt in ("-y", "--repeatY"): |
279 |
xsun |
3393 |
global repeatY |
280 |
gezelter |
3345 |
repeatY = int(arg) |
281 |
|
|
elif opt in ("-z", "--repeatZ"): |
282 |
xsun |
3393 |
global repeatZ |
283 |
gezelter |
3345 |
repeatZ = int(arg) |
284 |
|
|
if (_haveMDFileName != 1): |
285 |
|
|
usage() |
286 |
|
|
print "No meta-data file was specified" |
287 |
|
|
sys.exit() |
288 |
|
|
if (_haveOutputFileName != 1): |
289 |
|
|
usage() |
290 |
|
|
print "No output file was specified" |
291 |
|
|
sys.exit() |
292 |
|
|
readFile(mdFileName) |
293 |
|
|
mapToBox() |
294 |
|
|
writeFile(outputFileName) |
295 |
|
|
|
296 |
|
|
if __name__ == "__main__": |
297 |
|
|
if len(sys.argv) == 1: |
298 |
|
|
usage() |
299 |
|
|
sys.exit() |
300 |
|
|
main(sys.argv[1:]) |