ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/utilities/affineScale
(Generate patch)

Comparing trunk/src/applications/utilities/affineScale (file contents):
Revision 1064 by chrisfen, Tue Oct 10 14:36:02 2006 UTC vs.
Revision 1383 by gezelter, Thu Oct 22 19:12:14 2009 UTC

# Line 1 | Line 1
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:])

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines