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

Comparing branches/development/src/applications/utilities/principalAxisCalculator (file contents):
Revision 1700 by gezelter, Mon Sep 26 13:30:00 2011 UTC vs.
Revision 1701 by kstocke1, Thu Apr 5 19:37:58 2012 UTC

# Line 1 | Line 1
1 < #!@PYTHON_EXECUTABLE@
1 > #!/usr/bin/env python
2   """principalAxisCalculator
3 <
4 < Opens an XYZ file and computes the moments of inertia and principal axes
5 < for the structure in the XYZ file.  Optionally rotate the structure so that
6 < the long axis (that with the smallest eigenvalue) is pointing along the
7 < z-axis.
8 <
9 < Usage: principalAxisCalculator
10 <
11 < Options:
12 <  -h,  --help              show this help
13 <  -x,  --xyz=...           use specified XYZ (.xyz) file for the structure
14 <  -o,  --out=...           rotate the structure so that the smallest eigenvalue
15 <                           of the rotation matrix points along the z-axis.
16 <
17 < Example:
18 <   principalAxisCalculator -x junk.xyz -o rot.xyz
19 <
20 < """
21 <
22 < __author__ = "Dan Gezelter (gezelter@nd.edu)"
23 < __version__ = "$Revision$"
24 < __date__ = "$Date$"
25 < __copyright__ = "Copyright (c) 2006 by the University of Notre Dame"
26 < __license__ = "OpenMD"
27 <
28 < import sys
29 < import getopt
30 < import string
31 < import math
32 < import random
33 < from sets import *
34 < import numpy
35 <
36 <
37 < _haveXYZFileName = 0
38 < _haveOutFileName = 0
39 <
40 < positions = []
41 < indices = []
42 < atypes = []
43 < Hmass = 1.0079
44 < Cmass = 12.011
45 < Omass = 15.999
46 < Nmass = 14.007
47 < Smass = 32.066
3 >
4 > Opens an XYZ file and computes the moments of inertia and principal axes
5 > for the structure in the XYZ file.  Optionally rotate the structure so that
6 > the long axis (that with the smallest eigenvalue) is pointing along the
7 > z-axis.
8 >
9 > Usage: principalAxisCalculator
10 >
11 > Options:
12 >  -h,  --help              show this help
13 >  -x,  --xyz=...           use specified XYZ (.xyz) file for the structure
14 >  -o,  --out=...           rotate the structure so that the smallest eigenvalue
15 >                           of the rotation matrix points along the z-axis.
16 >
17 > Example:
18 >   principalAxisCalculator -x junk.xyz -o rot.xyz
19 >
20 > """
21 >
22 > __author__ = "Dan Gezelter (gezelter@nd.edu)"
23 > __version__ = "$Revision$"
24 > __date__ = "$Date$"
25 > __copyright__ = "Copyright (c) 2006 by the University of Notre Dame"
26 > __license__ = "OpenMD"
27 >
28 > import sys
29 > import getopt
30 > import string
31 > import math
32 > import random
33 > from sets import *
34 > import numpy
35 >
36 >
37 > _haveXYZFileName = 0
38 > _haveOutFileName = 0
39 >
40 > positions = []
41 > indices = []
42 > atypes = []
43 > Hmass = 1.0079
44 > Cmass = 12.011
45 > Omass = 15.999
46 > Nmass = 14.007
47 > Smass = 32.066
48 > Aumass = 196.466569
49   Au1mass = 196.466569
50   Au2mass =   0.5
51 <
52 < def usage():
53 <    print __doc__
54 <
55 < def add(x,y):
56 <    return x+y
57 <
58 < def sum(seq):
59 <    return reduce(add, seq)
60 <
61 < def readFile(XYZFileName):
62 <    print "reading XYZ file"
63 <
64 <    XYZFile = open(XYZFileName, 'r')        
65 <    # Find number of atoms first
66 <    line = XYZFile.readline()
67 <    L = line.split()
68 <    nAtoms = int(L[0])
69 <    # skip comment line
70 <    line = XYZFile.readline()
71 <    for i in range(nAtoms):
72 <        line = XYZFile.readline()
73 <        L = line.split()
74 <        myIndex = i
75 <        indices.append(myIndex)
76 <        atomType = L[0]
77 <        atypes.append(atomType)
78 <        x = float(L[1])
79 <        y = float(L[2])
80 <        z = float(L[3])
81 <        positions.append([x, y, z])
82 <    XYZFile.close()
83 <
84 <
85 < def findCOM():
86 <    #find center of mass
87 <    Xcom = 0.0
88 <    Ycom = 0.0
89 <    Zcom = 0.0
90 <    totalMass = 0.0
91 <
92 <    for i in range(0,len(indices)):
93 <        if (atypes[i] == "H"):
94 <            myMass = Hmass
95 <        elif (atypes[i] == "C"):
96 <            myMass = Cmass
97 <        elif (atypes[i] == "O"):
98 <            myMass = Omass
99 <        elif (atypes[i] == "N"):
100 <            myMass = Nmass
101 <        elif (atypes[i] == "S"):
102 <            myMass = Smass
103 <        elif (atypes[i] == "Au1"):
104 <            myMass = Au1mass
105 <        elif (atypes[i] == "Au2"):
106 <            myMass = Au2mass
107 <        else:
108 <            print "unknown atom type!"
109 <        
110 <        Xcom = Xcom + myMass * positions[i][0]
111 <        Ycom = Ycom + myMass * positions[i][1]
112 <        Zcom = Zcom + myMass * positions[i][2]
113 <        totalMass = totalMass + myMass
114 <
115 <    Xcom = Xcom / totalMass
116 <    Ycom = Ycom / totalMass
117 <    Zcom = Zcom / totalMass
118 <
119 <    COM = [Xcom, Ycom, Zcom]
120 <
121 <    return COM
122 <
123 < def findMoments():
124 <
125 <    COM = findCOM()
126 <    
127 <    #find inertia tensor matrix elements
128 <
129 <    I = numpy.zeros((3,3), numpy.float)
130 <    
131 <    for i in range(0,len(indices)):
132 <        if (atypes[i] == "H"):
133 <            myMass = Hmass
134 <        elif (atypes[i] == "C"):
135 <            myMass = Cmass
136 <        elif (atypes[i] == "O"):
137 <            myMass = Omass
138 <        elif (atypes[i] == "N"):
139 <            myMass = Nmass
140 <        elif (atypes[i] == "S"):
141 <            myMass = Smass
142 <        elif (atypes[i] == "Au1"):
143 <            myMass = Au1mass
144 <        elif (atypes[i] == "Au2"):
145 <            myMass = Au2mass
146 <        else:
147 <            print "unknown atom type!"
148 <
149 <        dx = positions[i][0] - COM[0]
150 <        dy = positions[i][1] - COM[1]
151 <        dz = positions[i][2] - COM[2]
152 <
153 <        I[0,0] = I[0,0] + myMass * ( dy * dy + dz * dz )
154 <        I[1,1] = I[1,1] + myMass * ( dx * dx + dz * dz )
155 <        I[2,2] = I[2,2] + myMass * ( dx * dx + dy * dy )
156 <
157 <        I[0,1] = I[0,1] - myMass * ( dx * dy )
158 <        I[0,2] = I[0,2] - myMass * ( dx * dz )
159 <
160 <        I[1,2] = I[1,2] - myMass * ( dy * dz )
161 <
162 <        I[1,0] = I[0,1]
163 <        I[2,0] = I[0,2]
164 <        I[2,1] = I[1,2]
165 <
166 <    print "Inertia Tensor:"
167 <    print I
168 <    print
169 <
170 <    (evals, evects) = numpy.linalg.eig(I)
171 <    print "evals:"
172 <    print evals
173 <    print
174 <    print "evects:"
175 <    print evects
176 <    print
177 <
178 <    return (COM, evals, evects)
179 <
180 < def writeFile(OutFileName):
181 <
182 <    (COM, evals, evects) = findMoments()
183 <
184 <
185 <    # we need to re-order the axes so that the smallest moment of inertia
186 <    # (which corresponds to the long axis of the molecule) is along the z-axis
187 <    # we'll just reverse the order of the three axes:
188 <    
189 <    axOrder = numpy.argsort(evals)    
190 <    RotMat = numpy.zeros((3,3), numpy.float)
191 <    RotMat[0] = evects[axOrder[2]]
192 <    RotMat[1] = evects[axOrder[1]]
193 <    RotMat[2] = evects[axOrder[0]]
194 <
195 <    nAtoms = len(indices)
196 <    
197 <    print "writing output XYZ file"
198 <
199 <    OutFile = open(OutFileName, 'w')        
200 <
201 <    OutFile.write('%10d\n' % (nAtoms))
202 <    OutFile.write('\n')
203 <
204 <    for i in range(nAtoms):
205 <
206 <        dx = positions[i][0] - COM[0]
207 <        dy = positions[i][1] - COM[1]
208 <        dz = positions[i][2] - COM[2]
209 <
210 <        r = numpy.array([dx,dy,dz])
211 <        rnew = numpy.dot(RotMat, r)
212 <          
213 <        OutFile.write('%s\t%f\t%f\t%f\t%d\n' % (atypes[i], rnew[0], rnew[1], rnew[2], i))
214 <    OutFile.close()        
215 <
216 < def main(argv):                        
217 <    try:                                
218 <        opts, args = getopt.getopt(argv, "hx:o:", ["help", "xyz=", "out="])
219 <    except getopt.GetoptError:          
220 <        usage()                          
221 <        sys.exit(2)                    
222 <    for opt, arg in opts:                
223 <        if opt in ("-h", "--help"):      
224 <            usage()                    
225 <            sys.exit()                  
226 <        elif opt in ("-x", "--xyz"):
227 <            XYZFileName = arg
228 <            global _haveXYZFileName
229 <            _haveXYZFileName = 1
230 <        elif opt in ("-o", "--out"):
231 <            OutFileName = arg
232 <            global _haveOutFileName
233 <            _haveOutFileName = 1
234 <
235 <
236 <    if (_haveXYZFileName != 1):
237 <        usage()
238 <        print "No xyz file was specified"
239 <        sys.exit()
240 <
241 <    readFile(XYZFileName)
242 <
243 <    if (_haveOutFileName == 1):
244 <        writeFile(OutFileName)
245 <    else:
246 <        findMoments()
247 <        
248 < if __name__ == "__main__":
249 <    if len(sys.argv) == 1:
250 <        usage()
251 <        sys.exit()
252 <    main(sys.argv[1:])
51 >
52 > def usage():
53 >    print __doc__
54 >
55 > def add(x,y):
56 >    return x+y
57 >
58 > def sum(seq):
59 >    return reduce(add, seq)
60 >
61 > def readFile(XYZFileName):
62 >    print "reading XYZ file"
63 >
64 >    XYZFile = open(XYZFileName, 'r')        
65 >    # Find number of atoms first
66 >    line = XYZFile.readline()
67 >    L = line.split()
68 >    nAtoms = int(L[0])
69 >    # skip comment line
70 >    line = XYZFile.readline()
71 >    for i in range(nAtoms):
72 >        line = XYZFile.readline()
73 >        L = line.split()
74 >        myIndex = i
75 >        indices.append(myIndex)
76 >        atomType = L[0]
77 >        atypes.append(atomType)
78 >        x = float(L[1])
79 >        y = float(L[2])
80 >        z = float(L[3])
81 >        positions.append([x, y, z])
82 >    XYZFile.close()
83 >
84 >
85 > def findCOM():
86 >    #find center of mass
87 >    Xcom = 0.0
88 >    Ycom = 0.0
89 >    Zcom = 0.0
90 >    totalMass = 0.0
91 >
92 >    for i in range(0,len(indices)):
93 >        if (atypes[i] == "H"):
94 >            myMass = Hmass
95 >        elif (atypes[i] == "C"):
96 >            myMass = Cmass
97 >        elif (atypes[i] == "O"):
98 >            myMass = Omass
99 >        elif (atypes[i] == "N"):
100 >            myMass = Nmass
101 >        elif (atypes[i] == "S"):
102 >            myMass = Smass
103 >        elif (atypes[i] == "Au1"):
104 >            myMass = Au1mass
105 >        elif (atypes[i] == "Au2"):
106 >            myMass = Au2mass
107 >        elif (atypes[i] == "Au"):
108 >            myMass = Aumass
109 >        else:
110 >            print "unknown atom type! %s" % (atypes[i])
111 >        
112 >        Xcom = Xcom + myMass * positions[i][0]
113 >        Ycom = Ycom + myMass * positions[i][1]
114 >        Zcom = Zcom + myMass * positions[i][2]
115 >        totalMass = totalMass + myMass
116 >
117 >    Xcom = Xcom / totalMass
118 >    Ycom = Ycom / totalMass
119 >    Zcom = Zcom / totalMass
120 >
121 >    COM = [Xcom, Ycom, Zcom]
122 >
123 >    return COM
124 >
125 > def findMoments():
126 >
127 >    COM = findCOM()
128 >    
129 >    #find inertia tensor matrix elements
130 >
131 >    I = numpy.zeros((3,3), numpy.float)
132 >    
133 >    for i in range(0,len(indices)):
134 >        if (atypes[i] == "H"):
135 >            myMass = Hmass
136 >        elif (atypes[i] == "C"):
137 >            myMass = Cmass
138 >        elif (atypes[i] == "O"):
139 >            myMass = Omass
140 >        elif (atypes[i] == "N"):
141 >            myMass = Nmass
142 >        elif (atypes[i] == "S"):
143 >            myMass = Smass
144 >        elif (atypes[i] == "Au1"):
145 >            myMass = Au1mass
146 >        elif (atypes[i] == "Au2"):
147 >            myMass = Au2mass
148 >        elif (atypes[i] == "Au"):
149 >            myMass = Aumass
150 >        else:
151 >            print "unknown atom type!"
152 >
153 >        dx = positions[i][0] - COM[0]
154 >        dy = positions[i][1] - COM[1]
155 >        dz = positions[i][2] - COM[2]
156 >
157 >        I[0,0] = I[0,0] + myMass * ( dy * dy + dz * dz )
158 >        I[1,1] = I[1,1] + myMass * ( dx * dx + dz * dz )
159 >        I[2,2] = I[2,2] + myMass * ( dx * dx + dy * dy )
160 >
161 >        I[0,1] = I[0,1] - myMass * ( dx * dy )
162 >        I[0,2] = I[0,2] - myMass * ( dx * dz )
163 >
164 >        I[1,2] = I[1,2] - myMass * ( dy * dz )
165 >
166 >        I[1,0] = I[0,1]
167 >        I[2,0] = I[0,2]
168 >        I[2,1] = I[1,2]
169 >
170 >    print "Inertia Tensor:"
171 >    print I
172 >    print
173 >
174 >    (evals, evects) = numpy.linalg.eig(I)
175 >    print "evals:"
176 >    print evals
177 >    print
178 >    print "evects:"
179 >    print evects
180 >    print
181 >
182 >    return (COM, evals, evects)
183 >
184 > def writeFile(OutFileName):
185 >
186 >    (COM, evals, evects) = findMoments()
187 >
188 >
189 >    # we need to re-order the axes so that the smallest moment of inertia
190 >    # (which corresponds to the long axis of the molecule) is along the z-axis
191 >    # we'll just reverse the order of the three axes:
192 >    
193 >    axOrder = numpy.argsort(evals)    
194 >    RotMat = numpy.zeros((3,3), numpy.float)
195 >    RotMat[0] = evects[axOrder[2]]
196 >    RotMat[1] = evects[axOrder[1]]
197 >    RotMat[2] = evects[axOrder[0]]
198 >
199 >    nAtoms = len(indices)
200 >    
201 >    print "writing output XYZ file"
202 >
203 >    OutFile = open(OutFileName, 'w')        
204 >
205 >    OutFile.write('%10d\n' % (nAtoms))
206 >    OutFile.write('\n')
207 >
208 >    for i in range(nAtoms):
209 >
210 >        dx = positions[i][0] - COM[0]
211 >        dy = positions[i][1] - COM[1]
212 >        dz = positions[i][2] - COM[2]
213 >
214 >        r = numpy.array([dx,dy,dz])
215 >        rnew = numpy.dot(RotMat, r)
216 >          
217 >        OutFile.write('%s\t%f\t%f\t%f\t%d\n' % (atypes[i], rnew[0], rnew[1], rnew[2], i))
218 >    OutFile.close()        
219 >
220 > def main(argv):                        
221 >    try:                                
222 >        opts, args = getopt.getopt(argv, "hx:o:", ["help", "xyz=", "out="])
223 >    except getopt.GetoptError:          
224 >        usage()                          
225 >        sys.exit(2)                    
226 >    for opt, arg in opts:                
227 >        if opt in ("-h", "--help"):      
228 >            usage()                    
229 >            sys.exit()                  
230 >        elif opt in ("-x", "--xyz"):
231 >            XYZFileName = arg
232 >            global _haveXYZFileName
233 >            _haveXYZFileName = 1
234 >        elif opt in ("-o", "--out"):
235 >            OutFileName = arg
236 >            global _haveOutFileName
237 >            _haveOutFileName = 1
238 >
239 >
240 >    if (_haveXYZFileName != 1):
241 >        usage()
242 >        print "No xyz file was specified"
243 >        sys.exit()
244 >
245 >    readFile(XYZFileName)
246 >
247 >    if (_haveOutFileName == 1):
248 >        writeFile(OutFileName)
249 >    else:
250 >        findMoments()
251 >        
252 > if __name__ == "__main__":
253 >    if len(sys.argv) == 1:
254 >        usage()
255 >        sys.exit()
256 >    main(sys.argv[1:])

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines