1 |
< |
#!/usr/bin/env perl |
1 |
> |
#!/usr/bin/env python |
2 |
> |
"""MetaData affine scaling transform |
3 |
|
|
4 |
< |
# program that scales an OOPSE .in or .eor file to a new volume |
4 |
> |
Takes a MetaData file and scales both the periodic box and the |
5 |
> |
coordinates of all StuntDoubles in the system by the same amount. |
6 |
|
|
7 |
< |
# author = "Chris Fennell |
8 |
< |
# version = "$Revision: 1.2 $" |
7 |
< |
# date = "$Date: 2006-10-10 14:36:02 $" |
8 |
< |
# copyright = "Copyright (c) 2006 by the University of Notre Dame" |
9 |
< |
# license = "OOPSE" |
7 |
> |
You can either specify a new volume scaling for isotropic scaling, |
8 |
> |
or specify one (or more) of the coordinates for non-isotropic scaling. |
9 |
|
|
10 |
+ |
Usage: affineScale |
11 |
|
|
12 |
< |
use Getopt::Std; |
12 |
> |
Options: |
13 |
> |
-h, --help show this help |
14 |
> |
-m, --meta-data=... use specified meta-data (.md, .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 |
|
|
21 |
< |
# some variables to get things going |
22 |
< |
$startSnap = 0; |
16 |
< |
$startFrame = 0; |
17 |
< |
$startStunts = 0; |
21 |
> |
Example: |
22 |
> |
affineScale -m lipidSystem.md -o scaledSystem.md -v 77000.0 |
23 |
|
|
24 |
< |
# get our options |
20 |
< |
getopts('hv', \%opts); |
24 |
> |
""" |
25 |
|
|
26 |
< |
# if we don't have a filename, drop to -h |
27 |
< |
$opts{h} = 'true' if $#ARGV != 1; |
26 |
> |
__author__ = "Dan Gezelter (gezelter@nd.edu)" |
27 |
> |
__version__ = "$Revision: 1.3 $" |
28 |
> |
__date__ = "$Date: 2009-10-22 19:12:14 $" |
29 |
> |
__copyright__ = "Copyright (c) 2009 by the University of Notre Dame" |
30 |
> |
__license__ = "OpenMD" |
31 |
|
|
32 |
< |
# our option output |
33 |
< |
if ($opts{h}){ |
34 |
< |
print "affineScale: performs an affine transform on an OOPSE .in or .eor\n"; |
35 |
< |
print "\tand writes to a new file named [file name].scale\n\n"; |
36 |
< |
print "usage: affineScale [-hv] [file name] [new volume (Ang^3)]\n\n"; |
37 |
< |
print " -h : show this message\n"; |
31 |
< |
die " -v : more verbose output\n"; |
32 |
< |
} |
32 |
> |
import sys |
33 |
> |
import getopt |
34 |
> |
import string |
35 |
> |
import math |
36 |
> |
import random |
37 |
> |
from sets import * |
38 |
|
|
39 |
< |
# set some variables to be used in the code |
40 |
< |
$fileName = $ARGV[0]; |
41 |
< |
$newVolume = $ARGV[1]; |
39 |
> |
metaData = [] |
40 |
> |
frameData = [] |
41 |
> |
indices = [] |
42 |
> |
pvqj = [] |
43 |
> |
p = [] |
44 |
> |
v = [] |
45 |
> |
q = [] |
46 |
> |
j = [] |
47 |
> |
Hmat = [] |
48 |
> |
BoxInv = [] |
49 |
|
|
50 |
< |
# some crazy input checking |
51 |
< |
if ($newVolume =~ /^[0-9]/){ |
40 |
< |
} else { |
41 |
< |
die "\t[new volume] value ($newVolume) is not a valid number\n\tPlease choose a non-negative numeric value for [new volume]\n"; |
42 |
< |
} |
50 |
> |
def usage(): |
51 |
> |
print __doc__ |
52 |
|
|
53 |
< |
# split up the name to obtain an output filename |
54 |
< |
@names = split('\.', $fileName); |
55 |
< |
$names[$#names] = 'scale'; |
56 |
< |
$outName = join('.',@names); |
53 |
> |
def readFile(mdFileName): |
54 |
> |
mdFile = open(mdFileName, 'r') |
55 |
> |
# Find OOPSE version info first |
56 |
> |
line = mdFile.readline() |
57 |
> |
while 1: |
58 |
> |
if '<OOPSE version=' in line: |
59 |
> |
OOPSEversion = line |
60 |
> |
break |
61 |
> |
line = mdFile.readline() |
62 |
> |
|
63 |
> |
# Rewind file and find start of MetaData block |
64 |
> |
|
65 |
> |
mdFile.seek(0) |
66 |
> |
line = mdFile.readline() |
67 |
> |
print "reading MetaData" |
68 |
> |
while 1: |
69 |
> |
if '<MetaData>' in line: |
70 |
> |
while 2: |
71 |
> |
metaData.append(line) |
72 |
> |
line = mdFile.readline() |
73 |
> |
if '</MetaData>' in line: |
74 |
> |
metaData.append(line) |
75 |
> |
break |
76 |
> |
break |
77 |
> |
line = mdFile.readline() |
78 |
|
|
79 |
< |
open(STATFILE, "./$fileName") || die "\tError: can't find file $fileName\n"; |
80 |
< |
open (SCALEFILE, ">./$outName") || die "\tError: can't open $outName"; |
79 |
> |
mdFile.seek(0) |
80 |
> |
print "reading Snapshot" |
81 |
> |
line = mdFile.readline() |
82 |
> |
while 1: |
83 |
> |
if '<Snapshot>' in line: |
84 |
> |
line = mdFile.readline() |
85 |
> |
while 1: |
86 |
> |
print "reading FrameData" |
87 |
> |
if '<FrameData>' in line: |
88 |
> |
while 2: |
89 |
> |
frameData.append(line) |
90 |
> |
if 'Hmat:' in line: |
91 |
> |
L = line.split() |
92 |
> |
Hxx = float(L[2].strip(',')) |
93 |
> |
Hxy = float(L[3].strip(',')) |
94 |
> |
Hxz = float(L[4].strip(',')) |
95 |
> |
Hyx = float(L[7].strip(',')) |
96 |
> |
Hyy = float(L[8].strip(',')) |
97 |
> |
Hyz = float(L[9].strip(',')) |
98 |
> |
Hzx = float(L[12].strip(',')) |
99 |
> |
Hzy = float(L[13].strip(',')) |
100 |
> |
Hzz = float(L[14].strip(',')) |
101 |
> |
Hmat.append([Hxx, Hxy, Hxz]) |
102 |
> |
Hmat.append([Hyx, Hyy, Hyz]) |
103 |
> |
Hmat.append([Hzx, Hzy, Hzz]) |
104 |
> |
BoxInv.append(1.0/Hxx) |
105 |
> |
BoxInv.append(1.0/Hyy) |
106 |
> |
BoxInv.append(1.0/Hzz) |
107 |
> |
line = mdFile.readline() |
108 |
> |
if '</FrameData>' in line: |
109 |
> |
frameData.append(line) |
110 |
> |
break |
111 |
> |
break |
112 |
|
|
113 |
< |
while (<STATFILE>){ |
114 |
< |
$startSnap = 0 if /\/Snapshot/; |
115 |
< |
$startFrame = 0 if /\/FrameData/; |
116 |
< |
$startStunts = 0 if /\/StuntDoubles/; |
113 |
> |
line = mdFile.readline() |
114 |
> |
while 1: |
115 |
> |
if '<StuntDoubles>' in line: |
116 |
> |
line = mdFile.readline() |
117 |
> |
while 2: |
118 |
> |
L = line.split() |
119 |
> |
myIndex = int(L[0]) |
120 |
> |
indices.append(myIndex) |
121 |
> |
pvqj.append(L[1]) |
122 |
> |
x = float(L[2]) |
123 |
> |
y = float(L[3]) |
124 |
> |
z = float(L[4]) |
125 |
> |
p.append([x, y, z]) |
126 |
> |
vx = float(L[5]) |
127 |
> |
vy = float(L[6]) |
128 |
> |
vz = float(L[7]) |
129 |
> |
v.append([vx, vy, vz]) |
130 |
> |
if 'pvqj' in L[1]: |
131 |
> |
qw = float(L[8]) |
132 |
> |
qx = float(L[9]) |
133 |
> |
qy = float(L[10]) |
134 |
> |
qz = float(L[11]) |
135 |
> |
q.append([qw, qx, qy, qz]) |
136 |
> |
jx = float(L[12]) |
137 |
> |
jy = float(L[13]) |
138 |
> |
jz = float(L[14]) |
139 |
> |
j.append([jx, jy, jz]) |
140 |
> |
else: |
141 |
> |
q.append([0.0, 0.0, 0.0, 0.0]) |
142 |
> |
j.append([0.0, 0.0, 0.0]) |
143 |
> |
|
144 |
> |
line = mdFile.readline() |
145 |
> |
if '</StuntDoubles>' in line: |
146 |
> |
break |
147 |
> |
break |
148 |
> |
line = mdFile.readline() |
149 |
> |
if not line: break |
150 |
> |
|
151 |
> |
mdFile.close() |
152 |
|
|
153 |
< |
if ($startSnap == 1){ |
154 |
< |
if ($startFrame == 1){ |
59 |
< |
if (/Hmat/){ |
60 |
< |
@line = split; |
61 |
< |
|
62 |
< |
chop $line[2]; |
63 |
< |
chop $line[8]; |
64 |
< |
$oldVolume = $line[2]*$line[8]*$line[14]; |
65 |
< |
|
66 |
< |
$scale = ($newVolume/$oldVolume)**0.333333333333333; |
67 |
< |
# scale the hmat vectors (only orthorhombic) |
68 |
< |
$line[2] *= $scale; |
69 |
< |
$line[8] *= $scale; |
70 |
< |
$line[14] *= $scale; |
71 |
< |
$volume = $line[2]*$line[8]*$line[14]; |
72 |
< |
|
73 |
< |
print SCALEFILE "\t$line[0]"; |
74 |
< |
for ($i=1; $i<=$#line; $i++) { |
75 |
< |
print SCALEFILE " $line[$i]"; |
76 |
< |
if ($i == 2 || $i == 8) {print SCALEFILE ",";} |
77 |
< |
} |
78 |
< |
print SCALEFILE "\n"; |
79 |
< |
} else { |
80 |
< |
print SCALEFILE; |
81 |
< |
} |
82 |
< |
} |
83 |
< |
if ($startStunts == 1){ |
84 |
< |
@line = split; |
85 |
< |
$line[2] *= $scale; |
86 |
< |
$line[3] *= $scale; |
87 |
< |
$line[4] *= $scale; |
88 |
< |
print SCALEFILE "$line[0]\t$line[1]\t"; |
89 |
< |
for ($i=2; $i<=$#line; $i++) { |
90 |
< |
print SCALEFILE "$line[$i] "; |
91 |
< |
} |
92 |
< |
print SCALEFILE "\n"; |
93 |
< |
} |
94 |
< |
if ($startFrame == 0 && $startStunts == 0) { |
95 |
< |
print SCALEFILE; |
96 |
< |
} |
97 |
< |
} else { |
98 |
< |
print SCALEFILE; |
99 |
< |
} |
100 |
< |
$startSnap = 1 if /Snapshot/; |
101 |
< |
$startFrame = 1 if /FrameData/; |
102 |
< |
$startStunts = 1 if /StuntDoubles/; |
103 |
< |
# check again since "/value" contains "value" |
104 |
< |
$startSnap = 0 if /\/Snapshot/; |
105 |
< |
$startFrame = 0 if /\/FrameData/; |
106 |
< |
$startStunts = 0 if /\/StuntDoubles/; |
107 |
< |
} |
153 |
> |
def writeFile(outputFileName,repeatX,repeatY,repeatZ): |
154 |
> |
outputFile = open(outputFileName, 'w') |
155 |
|
|
156 |
< |
print "Affine transformed configuration written to $outName\n"; |
156 |
> |
outputFile.write("<OOPSE version=4>\n"); |
157 |
|
|
158 |
< |
if (defined($opts{v})){ |
159 |
< |
print "New Volume : $volume\n"; |
160 |
< |
print "Old Volume : $oldVolume\n"; |
161 |
< |
print "Scale Factor : $scale\n"; |
162 |
< |
} |
158 |
> |
print "writing MetaData" |
159 |
> |
for metaline in metaData: |
160 |
> |
if 'nMol' in metaline: |
161 |
> |
metasplit = metaline.split() |
162 |
> |
nMol = float(metasplit[2].strip(';')) |
163 |
> |
newNmol = nMol * repeatX * repeatY * repeatZ |
164 |
> |
outputFile.write(' nMol = %10d;\n' % (newNmol)) |
165 |
> |
else: |
166 |
> |
outputFile.write(metaline) |
167 |
|
|
168 |
+ |
print "writing Snapshot" |
169 |
+ |
outputFile.write(" <Snapshot>\n") |
170 |
+ |
|
171 |
+ |
print "writing FrameData" |
172 |
+ |
for frameline in frameData: |
173 |
+ |
if 'Hmat:' in frameline: |
174 |
+ |
myH = [] |
175 |
+ |
myH.append([repeatX * Hmat[0][0], repeatX * Hmat[0][1], repeatX * Hmat[0][2]]) |
176 |
+ |
myH.append([repeatY * Hmat[1][0], repeatY * Hmat[1][1], repeatY * Hmat[1][2]]) |
177 |
+ |
myH.append([repeatZ * Hmat[2][0], repeatZ * Hmat[2][1], repeatZ * Hmat[2][2]]) |
178 |
+ |
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])) |
179 |
+ |
else: |
180 |
+ |
outputFile.write(frameline) |
181 |
+ |
|
182 |
+ |
print "writing StuntDoubles" |
183 |
+ |
outputFile.write(" <StuntDoubles>\n") |
184 |
+ |
|
185 |
+ |
print repeatX, repeatY, repeatZ |
186 |
+ |
|
187 |
+ |
deco = [ (index, i) for i, index in enumerate(indices) ] |
188 |
+ |
deco.sort() |
189 |
+ |
whichSD = 0 |
190 |
+ |
for index in range(len(deco)): |
191 |
+ |
(index,i) = deco[index] |
192 |
+ |
print i |
193 |
+ |
for ii in range(repeatX): |
194 |
+ |
for jj in range(repeatY): |
195 |
+ |
for kk in range(repeatZ): |
196 |
+ |
myP = [] |
197 |
+ |
myP.append(p[i][0] + ii*Hmat[0][0] + jj*Hmat[1][0] + kk*Hmat[2][0]) |
198 |
+ |
myP.append(p[i][1] + ii*Hmat[0][1] + jj*Hmat[1][1] + kk*Hmat[2][1]) |
199 |
+ |
myP.append(p[i][2] + ii*Hmat[0][2] + jj*Hmat[1][2] + kk*Hmat[2][2]) |
200 |
+ |
|
201 |
+ |
if (pvqj[i] == 'pv'): |
202 |
+ |
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])) |
203 |
+ |
elif (pvqj[i] == 'pvqj'): |
204 |
+ |
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])) |
205 |
+ |
whichSD = whichSD + 1 |
206 |
+ |
|
207 |
+ |
outputFile.write(" </StuntDoubles>\n") |
208 |
+ |
outputFile.write(" </Snapshot>\n") |
209 |
+ |
outputFile.write("</OOPSE>\n") |
210 |
+ |
outputFile.close() |
211 |
+ |
|
212 |
+ |
def roundMe(x): |
213 |
+ |
if (x >= 0.0): |
214 |
+ |
return math.floor(x + 0.5) |
215 |
+ |
else: |
216 |
+ |
return math.ceil(x - 0.5) |
217 |
+ |
|
218 |
+ |
def wrapVector(myVect): |
219 |
+ |
scaled = [0.0, 0.0, 0.0] |
220 |
+ |
for i in range(3): |
221 |
+ |
scaled[i] = myVect[i] * BoxInv[i] |
222 |
+ |
scaled[i] = scaled[i] - roundMe(scaled[i]) |
223 |
+ |
myVect[i] = scaled[i] * Hmat[i][i] |
224 |
+ |
return myVect |
225 |
+ |
|
226 |
+ |
def dot(L1, L2): |
227 |
+ |
myDot = 0.0 |
228 |
+ |
for i in range(len(L1)): |
229 |
+ |
myDot = myDot + L1[i]*L2[i] |
230 |
+ |
return myDot |
231 |
+ |
|
232 |
+ |
def normalize(L1): |
233 |
+ |
L2 = [] |
234 |
+ |
myLength = math.sqrt(dot(L1, L1)) |
235 |
+ |
for i in range(len(L1)): |
236 |
+ |
L2.append(L1[i] / myLength) |
237 |
+ |
return L2 |
238 |
+ |
|
239 |
+ |
def cross(L1, L2): |
240 |
+ |
# don't call this with anything other than length 3 lists please |
241 |
+ |
# or you'll be sorry |
242 |
+ |
L3 = [0.0, 0.0, 0.0] |
243 |
+ |
L3[0] = L1[1]*L2[2] - L1[2]*L2[1] |
244 |
+ |
L3[1] = L1[2]*L2[0] - L1[0]*L2[2] |
245 |
+ |
L3[2] = L1[0]*L2[1] - L1[1]*L2[0] |
246 |
+ |
return L3 |
247 |
+ |
|
248 |
+ |
def mapToBox(): |
249 |
+ |
for i in range(len(indices)): |
250 |
+ |
dpos = [] |
251 |
+ |
dpos.append(p[i][0]) |
252 |
+ |
dpos.append(p[i][1]) |
253 |
+ |
dpos.append(p[i][2]) |
254 |
+ |
p[i] = wrapVector(dpos) |
255 |
+ |
|
256 |
+ |
def scaleBox(scaleX, scaleY, scaleZ): |
257 |
+ |
for i in range(3): |
258 |
+ |
Hmat[0][i] = scaleX * Hmat[0][i] |
259 |
+ |
for i in range(3): |
260 |
+ |
Hmat[1][i] = scaleY * Hmat[1][i] |
261 |
+ |
for i in range(3): |
262 |
+ |
Hmat[2][i] = scaleZ * Hmat[2][i] |
263 |
+ |
for i in range(3): |
264 |
+ |
BoxInv[i] = 1.0/Hmat[i][i] |
265 |
+ |
|
266 |
+ |
def scaleCoordinates(scaleX, scaleY, scaleZ): |
267 |
+ |
for i in range(len(indices)): |
268 |
+ |
p[i][0] = p[i][0]*scaleX |
269 |
+ |
p[i][1] = p[i][1]*scaleY |
270 |
+ |
p[i][2] = p[i][2]*scaleZ |
271 |
+ |
|
272 |
+ |
def main(argv): |
273 |
+ |
_haveMDFileName = 0 |
274 |
+ |
_haveOutputFileName = 0 |
275 |
+ |
try: |
276 |
+ |
opts, args = getopt.getopt(argv, "hm:o:x:y:z:v:", ["help", "meta-data=", "output-file=", "newX=", "newY=", "newZ=","newV="]) |
277 |
+ |
except getopt.GetoptError: |
278 |
+ |
usage() |
279 |
+ |
sys.exit(2) |
280 |
+ |
doV = 0 |
281 |
+ |
doX = 0 |
282 |
+ |
doY = 0 |
283 |
+ |
doZ = 0 |
284 |
+ |
for opt, arg in opts: |
285 |
+ |
if opt in ("-h", "--help"): |
286 |
+ |
usage() |
287 |
+ |
sys.exit() |
288 |
+ |
elif opt in ("-m", "--meta-data"): |
289 |
+ |
mdFileName = arg |
290 |
+ |
_haveMDFileName = 1 |
291 |
+ |
elif opt in ("-o", "--output-file"): |
292 |
+ |
outputFileName = arg |
293 |
+ |
_haveOutputFileName = 1 |
294 |
+ |
if opt in ("-x", "--newX"): |
295 |
+ |
newX = float(arg) |
296 |
+ |
doX = 1 |
297 |
+ |
elif opt in ("-y", "--newY"): |
298 |
+ |
newY = float(arg) |
299 |
+ |
doY = 1 |
300 |
+ |
elif opt in ("-z", "--newZ"): |
301 |
+ |
newZ = float(arg) |
302 |
+ |
doZ = 1 |
303 |
+ |
elif opt in ("-v", "--newV"): |
304 |
+ |
newV = float(arg) |
305 |
+ |
doV = 1 |
306 |
+ |
|
307 |
+ |
if (_haveMDFileName != 1): |
308 |
+ |
usage() |
309 |
+ |
print "No meta-data file was specified" |
310 |
+ |
sys.exit() |
311 |
+ |
if (_haveOutputFileName != 1): |
312 |
+ |
usage() |
313 |
+ |
print "No output file was specified" |
314 |
+ |
sys.exit() |
315 |
+ |
|
316 |
+ |
if not (doV or doX or doY or doZ): |
317 |
+ |
usage() |
318 |
+ |
print "no scaling options given. Nothing to do!" |
319 |
+ |
sys.exit() |
320 |
+ |
|
321 |
+ |
if doV and (doX or doY or doZ): |
322 |
+ |
usage() |
323 |
+ |
print "-v is mutually exclusive with any of the -x, -y, and -z options" |
324 |
+ |
sys.exit() |
325 |
+ |
|
326 |
+ |
readFile(mdFileName) |
327 |
+ |
|
328 |
+ |
scaleX = 1.0 |
329 |
+ |
scaleY = 1.0 |
330 |
+ |
scaleZ = 1.0 |
331 |
+ |
|
332 |
+ |
if doX: |
333 |
+ |
scaleX = newX / Hmat[0][0] |
334 |
+ |
if doY: |
335 |
+ |
scaleY = newY / Hmat[1][1] |
336 |
+ |
if doZ: |
337 |
+ |
scaleZ = newZ / Hmat[2][2] |
338 |
+ |
|
339 |
+ |
if doV: |
340 |
+ |
oldV = Hmat[0][0] * Hmat[1][1] * Hmat[2][2] |
341 |
+ |
scaleX = pow(newV/oldV, 1.0/3.0) |
342 |
+ |
scaleY = scaleX |
343 |
+ |
scaleZ = scaleX |
344 |
+ |
|
345 |
+ |
mapToBox() |
346 |
+ |
scaleBox(scaleX, scaleY, scaleZ) |
347 |
+ |
scaleCoordinates(scaleX, scaleY, scaleZ) |
348 |
+ |
writeFile(outputFileName, 1, 1, 1) |
349 |
+ |
|
350 |
+ |
if __name__ == "__main__": |
351 |
+ |
if len(sys.argv) == 1: |
352 |
+ |
usage() |
353 |
+ |
sys.exit() |
354 |
+ |
main(sys.argv[1:]) |