[ViewVC] Annotation of: OpenMD/branches/development/src/applications/utilities/affineScale

#!/usr/bin/env python
"""MetaData affine scaling transform

Takes a MetaData file and scales both the periodic box and the
coordinates of all StuntDoubles in the system by the same amount.

You can either specify a new volume scaling for isotropic scaling,
or specify one (or more) of the coordinates for non-isotropic scaling.

Usage: affineScale

Options:
  -h,  --help             show this help
  -m,  --meta-data=...    use specified meta-data (.md, .eor) file
  -o,  --output-file=...  use specified output file
  -x,  --newX=...         scale the system to a new x dimension
  -y,  --newY=...         scale the system to a new y dimension
  -z,  --newZ=...         scale the system to a new z dimension
  -v,  --newV=...         scale the system to a new volume

Example:
  affineScale -m lipidSystem.md -o scaledSystem.md -v 77000.0

"""

__author__ = "Dan Gezelter (gezelter@nd.edu)"
__version__ = "$Revision$"
__date__ = "$Date$"
__copyright__ = "Copyright (c) 2009 by the University of Notre Dame"
__license__ = "OpenMD"

import sys
import getopt
import string
import math
import random
from sets import *

metaData = []
frameData = []
indices = []
pvqj = []
p = []
v = []
q = []
j = []
Hmat = []
BoxInv = []

def usage():
    print __doc__

def readFile(mdFileName):
    mdFile = open(mdFileName, 'r')        
    # Find OpenMD version info first
    line = mdFile.readline()
    while 1:
        if '<OOPSE version=' in line or '<OpenMD version=' in line:
            OpenMDversion = line
            break
        line = mdFile.readline()
        
        # Rewind file and find start of MetaData block
        
    mdFile.seek(0)
    line = mdFile.readline()
    print "reading MetaData"
    while 1:
        if '<MetaData>' in line:
            while 2:
                metaData.append(line)
                line = mdFile.readline()
                if '</MetaData>' in line:
                    metaData.append(line)
                    break
            break
        line = mdFile.readline()

    mdFile.seek(0)
    print "reading Snapshot"
    line = mdFile.readline()
    while 1:
        if '<Snapshot>' in line:
            line = mdFile.readline()
            while 1:
                print "reading FrameData"
                if '<FrameData>' in line:
                    while 2:
                        frameData.append(line)
                        if 'Hmat:' in line:
                            L = line.split()
                            Hxx = float(L[2].strip(','))
                            Hxy = float(L[3].strip(','))
                            Hxz = float(L[4].strip(','))
                            Hyx = float(L[7].strip(','))
                            Hyy = float(L[8].strip(','))
                            Hyz = float(L[9].strip(','))
                            Hzx = float(L[12].strip(','))
                            Hzy = float(L[13].strip(','))
                            Hzz = float(L[14].strip(','))
                            Hmat.append([Hxx, Hxy, Hxz])
                            Hmat.append([Hyx, Hyy, Hyz])
                            Hmat.append([Hzx, Hzy, Hzz])
                            BoxInv.append(1.0/Hxx)
                            BoxInv.append(1.0/Hyy)
                            BoxInv.append(1.0/Hzz)
                        line = mdFile.readline()
                        if '</FrameData>' in line:
                            frameData.append(line)
                            break
                    break

            line = mdFile.readline()
            while 1:
                if '<StuntDoubles>' in line:
                    line = mdFile.readline()
                    while 2:
                        L = line.split()
                        myIndex = int(L[0])
                        indices.append(myIndex)
                        pvqj.append(L[1])
                        x = float(L[2])
                        y = float(L[3])
                        z = float(L[4])
                        p.append([x, y, z])
                        vx = float(L[5])
                        vy = float(L[6])
                        vz = float(L[7])
                        v.append([vx, vy, vz])
                        if 'pvqj' in L[1]:
                            qw = float(L[8])
                            qx = float(L[9])
                            qy = float(L[10])
                            qz = float(L[11])
                            q.append([qw, qx, qy, qz])
                            jx = float(L[12])
                            jy = float(L[13])
                            jz = float(L[14])
                            j.append([jx, jy, jz])
                        else:
                            q.append([0.0, 0.0, 0.0, 0.0])
                            j.append([0.0, 0.0, 0.0])
                                               
                        line = mdFile.readline()
                        if '</StuntDoubles>' in line:
                            break
                    break
        line = mdFile.readline()
        if not line: break
    
    mdFile.close()

def writeFile(outputFileName,repeatX,repeatY,repeatZ):
    outputFile = open(outputFileName, 'w')

    outputFile.write("<OpenMD version=1>\n");

    print "writing MetaData"
    for metaline in metaData:
        if 'nMol' in metaline:
            metasplit = metaline.split()
            nMol = float(metasplit[2].strip(';'))
            newNmol = nMol * repeatX * repeatY * repeatZ
            outputFile.write('  nMol = %10d;\n' % (newNmol))
        else:
            outputFile.write(metaline)

    print "writing Snapshot"
    outputFile.write("  <Snapshot>\n")

    print "writing FrameData"
    for frameline in frameData:
        if 'Hmat:' in frameline:
            myH = []
            myH.append([repeatX * Hmat[0][0], repeatX * Hmat[0][1], repeatX * Hmat[0][2]])
            myH.append([repeatY * Hmat[1][0], repeatY * Hmat[1][1], repeatY * Hmat[1][2]])
            myH.append([repeatZ * Hmat[2][0], repeatZ * Hmat[2][1], repeatZ * Hmat[2][2]])
            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]))
        else:
            outputFile.write(frameline)

    print "writing StuntDoubles"
    outputFile.write("    <StuntDoubles>\n")

    print repeatX, repeatY, repeatZ

    deco = [ (index, i) for i, index in enumerate(indices) ]
    deco.sort()
    whichSD = 0
    for index in range(len(deco)):
        (index,i) = deco[index]
        print i
        for ii in range(repeatX):
            for jj in range(repeatY):
                for kk in range(repeatZ):                    
                    myP = [] 
                    myP.append(p[i][0] + ii*Hmat[0][0] + jj*Hmat[1][0] + kk*Hmat[2][0])
                    myP.append(p[i][1] + ii*Hmat[0][1] + jj*Hmat[1][1] + kk*Hmat[2][1])
                    myP.append(p[i][2] + ii*Hmat[0][2] + jj*Hmat[1][2] + kk*Hmat[2][2])
                        
                    if (pvqj[i] == 'pv'):
                        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]))
                    elif (pvqj[i] == 'pvqj'):
                        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]))
                    whichSD = whichSD + 1

    outputFile.write("    </StuntDoubles>\n")
    outputFile.write("  </Snapshot>\n")
    outputFile.write("</OpenMD>\n")
    outputFile.close()

def roundMe(x):
    if (x >= 0.0):
        return math.floor(x + 0.5)
    else:
        return math.ceil(x - 0.5)

def wrapVector(myVect):
    scaled = [0.0, 0.0, 0.0]
    for i in range(3):
        scaled[i] = myVect[i] * BoxInv[i]
        scaled[i] = scaled[i] - roundMe(scaled[i])
        myVect[i] = scaled[i] * Hmat[i][i]
    return myVect

def dot(L1, L2):
    myDot = 0.0
    for i in range(len(L1)):
        myDot = myDot + L1[i]*L2[i]
    return myDot

def normalize(L1):
    L2 = []
    myLength = math.sqrt(dot(L1, L1))
    for i in range(len(L1)):
        L2.append(L1[i] / myLength)
    return L2

def cross(L1, L2):
    # don't call this with anything other than length 3 lists please
    # or you'll be sorry
    L3 = [0.0, 0.0, 0.0]
    L3[0] = L1[1]*L2[2] - L1[2]*L2[1] 
    L3[1] = L1[2]*L2[0] - L1[0]*L2[2]
    L3[2] = L1[0]*L2[1] - L1[1]*L2[0]
    return L3

def mapToBox():
    for i in range(len(indices)):
        dpos = []
        dpos.append(p[i][0])
        dpos.append(p[i][1])
        dpos.append(p[i][2])
        p[i] = wrapVector(dpos)

def scaleBox(scaleX, scaleY, scaleZ):
    for i in range(3):
        Hmat[0][i] = scaleX * Hmat[0][i]
    for i in range(3):
        Hmat[1][i] = scaleY * Hmat[1][i]
    for i in range(3):
        Hmat[2][i] = scaleZ * Hmat[2][i]
    for i in range(3):            
        BoxInv[i] = 1.0/Hmat[i][i]

def scaleCoordinates(scaleX, scaleY, scaleZ):
    for i in range(len(indices)):
        p[i][0] = p[i][0]*scaleX
        p[i][1] = p[i][1]*scaleY
        p[i][2] = p[i][2]*scaleZ

def main(argv):
    _haveMDFileName = 0
    _haveOutputFileName = 0
    try:                                
        opts, args = getopt.getopt(argv, "hm:o:x:y:z:v:", ["help", "meta-data=", "output-file=", "newX=", "newY=", "newZ=","newV="]) 
    except getopt.GetoptError:           
        usage()                          
        sys.exit(2)                     
    doV = 0
    doX = 0
    doY = 0
    doZ = 0
    for opt, arg in opts:                
        if opt in ("-h", "--help"):      
            usage()                     
            sys.exit()                  
        elif opt in ("-m", "--meta-data"): 
            mdFileName = arg
            _haveMDFileName = 1
        elif opt in ("-o", "--output-file"): 
            outputFileName = arg
            _haveOutputFileName = 1
        if opt in ("-x", "--newX"):
            newX = float(arg)
            doX = 1
        elif opt in ("-y", "--newY"): 
            newY = float(arg)
            doY = 1
        elif opt in ("-z", "--newZ"): 
            newZ = float(arg)
            doZ = 1
        elif opt in ("-v", "--newV"):        
            newV = float(arg)
            doV = 1

    if (_haveMDFileName != 1):
        usage() 
        print "No meta-data file was specified"
        sys.exit()
    if (_haveOutputFileName != 1):
        usage()
        print "No output file was specified"
        sys.exit()

    if not (doV or doX or doY or doZ):
        usage()
        print "no scaling options given.  Nothing to do!"
        sys.exit()

    if doV and (doX or doY or doZ):
        usage()
        print "-v is mutually exclusive with any of the -x, -y, and -z options"
        sys.exit()

    readFile(mdFileName)

    scaleX = 1.0
    scaleY = 1.0
    scaleZ = 1.0    

    if doX:
        scaleX = newX / Hmat[0][0]
    if doY:
        scaleY = newY / Hmat[1][1]
    if doZ:
        scaleZ = newZ / Hmat[2][2]

    if doV:
        oldV = Hmat[0][0] * Hmat[1][1] * Hmat[2][2]
        scaleX = pow(newV/oldV, 1.0/3.0)
        scaleY = scaleX
        scaleZ = scaleX

    mapToBox()
    scaleBox(scaleX, scaleY, scaleZ)
    scaleCoordinates(scaleX, scaleY, scaleZ)
    writeFile(outputFileName, 1, 1, 1)

if __name__ == "__main__":
    if len(sys.argv) == 1:
        usage()
        sys.exit()
    main(sys.argv[1:])
Revision:	1465
Committed:	Fri Jul 9 23:08:25 2010 UTC (14 years, 9 months ago) by chuckv
File size:	11631 byte(s)
Log Message:	Creating busticated version of OpenMD
#	User	Rev	Content
1	gezelter	1383	#!/usr/bin/env python
2			"""MetaData affine scaling transform
3	chrisfen	1063
4	gezelter	1383	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	chrisfen	1063
7	gezelter	1383	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	chrisfen	1063
10	gezelter	1383	Usage: affineScale
11	chrisfen	1063
12	gezelter	1383	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	chrisfen	1063
21	gezelter	1383	Example:
22			affineScale -m lipidSystem.md -o scaledSystem.md -v 77000.0
23	chrisfen	1063
24	gezelter	1383	"""
25	chrisfen	1063
26	gezelter	1383	__author__ = "Dan Gezelter (gezelter@nd.edu)"
27	gezelter	1442	__version__ = "$Revision$"
28			__date__ = "$Date$"
29	gezelter	1383	__copyright__ = "Copyright (c) 2009 by the University of Notre Dame"
30			__license__ = "OpenMD"
31	chrisfen	1063
32	gezelter	1383	import sys
33			import getopt
34			import string
35			import math
36			import random
37			from sets import *
38	chrisfen	1063
39	gezelter	1383	metaData = []
40			frameData = []
41			indices = []
42			pvqj = []
43			p = []
44			v = []
45			q = []
46			j = []
47			Hmat = []
48			BoxInv = []
49	chrisfen	1063
50	gezelter	1383	def usage():
51			print __doc__
52	chrisfen	1063
53	gezelter	1383	def readFile(mdFileName):
54			mdFile = open(mdFileName, 'r')
55	gezelter	1390	# Find OpenMD version info first
56	gezelter	1383	line = mdFile.readline()
57			while 1:
58	gezelter	1390	if '<OOPSE version=' in line or '<OpenMD version=' in line:
59			OpenMDversion = line
60	gezelter	1383	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	chrisfen	1063
79	gezelter	1383	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	chrisfen	1063
113	gezelter	1383	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	chrisfen	1063
153	gezelter	1383	def writeFile(outputFileName,repeatX,repeatY,repeatZ):
154			outputFile = open(outputFileName, 'w')
155	chrisfen	1063
156	gezelter	1390	outputFile.write("<OpenMD version=1>\n");
157	chrisfen	1063
158	gezelter	1383	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	chrisfen	1063
168	gezelter	1383	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] + iiHmat[0][0] + jjHmat[1][0] + kk*Hmat[2][0])
198			myP.append(p[i][1] + iiHmat[0][1] + jjHmat[1][1] + kk*Hmat[2][1])
199			myP.append(p[i][2] + iiHmat[0][2] + jjHmat[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	gezelter	1390	outputFile.write("</OpenMD>\n")
210	gezelter	1383	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:])