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 1797 by kstocke1, Thu Apr 5 19:37:58 2012 UTC vs.
Revision 1798 by gezelter, Thu Sep 13 14:10:11 2012 UTC

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