ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/applications/utilities/principalAxisCalculator
Revision: 3520
Committed: Mon Sep 7 16:31:51 2009 UTC (15 years, 10 months ago) by cli2
File size: 6201 byte(s)
Log Message:
Added new restraint infrastructure
Added MolecularRestraints
Added ObjectRestraints
Added RestraintStamp
Updated thermodynamic integration to use ObjectRestraints
Added Quaternion mathematics for twist swing decompositions
Significantly updated RestWriter and RestReader to use dump-like files
Added selections for x, y, and z coordinates of atoms
Removed monolithic Restraints class
Fixed a few bugs in gradients of Euler angles in DirectionalAtom and RigidBody
Added some rotational capabilities to prinicpalAxisCalculator

File Contents

# Content
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: 1.2 $"
24 __date__ = "$Date: 2009-09-07 16:31:49 $"
25 __copyright__ = "Copyright (c) 2006 by the University of Notre Dame"
26 __license__ = "OOPSE"
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
49 def usage():
50 print __doc__
51
52 def add(x,y):
53 return x+y
54
55 def sum(seq):
56 return reduce(add, seq)
57
58 def readFile(XYZFileName):
59 print "reading XYZ file"
60
61 XYZFile = open(XYZFileName, 'r')
62 # Find number of atoms first
63 line = XYZFile.readline()
64 L = line.split()
65 nAtoms = int(L[0])
66 # skip comment line
67 line = XYZFile.readline()
68 for i in range(nAtoms):
69 line = XYZFile.readline()
70 L = line.split()
71 myIndex = i
72 indices.append(myIndex)
73 atomType = L[0]
74 atypes.append(atomType)
75 x = float(L[1])
76 y = float(L[2])
77 z = float(L[3])
78 positions.append([x, y, z])
79 XYZFile.close()
80
81
82 def findCOM():
83 #find center of mass
84 Xcom = 0.0
85 Ycom = 0.0
86 Zcom = 0.0
87 totalMass = 0.0
88
89 for i in range(0,len(indices)):
90 if (atypes[i] == "H"):
91 myMass = Hmass
92 elif (atypes[i] == "C"):
93 myMass = Cmass
94 elif (atypes[i] == "O"):
95 myMass = Omass
96 elif (atypes[i] == "N"):
97 myMass = Nmass
98 elif (atypes[i] == "S"):
99 myMass = Smass
100 else:
101 print "unknown atom type!"
102
103 Xcom = Xcom + myMass * positions[i][0]
104 Ycom = Ycom + myMass * positions[i][1]
105 Zcom = Zcom + myMass * positions[i][2]
106 totalMass = totalMass + myMass
107
108 Xcom = Xcom / totalMass
109 Ycom = Ycom / totalMass
110 Zcom = Zcom / totalMass
111
112 COM = [Xcom, Ycom, Zcom]
113
114 return COM
115
116 def findMoments():
117
118 COM = findCOM()
119
120 #find inertia tensor matrix elements
121
122 I = numpy.zeros((3,3), numpy.float)
123
124 for i in range(0,len(indices)):
125 if (atypes[i] == "H"):
126 myMass = Hmass
127 elif (atypes[i] == "C"):
128 myMass = Cmass
129 elif (atypes[i] == "O"):
130 myMass = Omass
131 elif (atypes[i] == "N"):
132 myMass = Nmass
133 elif (atypes[i] == "S"):
134 myMass = Smass
135 else:
136 print "unknown atom type!"
137
138 dx = positions[i][0] - COM[0]
139 dy = positions[i][1] - COM[1]
140 dz = positions[i][2] - COM[2]
141
142 I[0,0] = I[0,0] + myMass * ( dy * dy + dz * dz )
143 I[1,1] = I[1,1] + myMass * ( dx * dx + dz * dz )
144 I[2,2] = I[2,2] + myMass * ( dx * dx + dy * dy )
145
146 I[0,1] = I[0,1] - myMass * ( dx * dy )
147 I[0,2] = I[0,2] - myMass * ( dx * dz )
148
149 I[1,2] = I[1,2] - myMass * ( dy * dz )
150
151 I[1,0] = I[0,1]
152 I[2,0] = I[0,2]
153 I[2,1] = I[1,2]
154
155 print "Inertia Tensor:"
156 print I
157 print
158
159 (evals, evects) = numpy.linalg.eig(I)
160 print "evals:"
161 print evals
162 print
163 print "evects:"
164 print evects
165 print
166
167 return (COM, evals, evects)
168
169 def writeFile(OutFileName):
170
171 (COM, evals, evects) = findMoments()
172
173
174 # we need to re-order the axes so that the smallest moment of inertia
175 # (which corresponds to the long axis of the molecule) is along the z-axis
176 # we'll just reverse the order of the three axes:
177
178 axOrder = numpy.argsort(evals)
179 RotMat = numpy.zeros((3,3), numpy.float)
180 RotMat[0] = evects[axOrder[2]]
181 RotMat[1] = evects[axOrder[1]]
182 RotMat[2] = evects[axOrder[0]]
183
184 nAtoms = len(indices)
185
186 print "writing output XYZ file"
187
188 OutFile = open(OutFileName, 'w')
189
190 OutFile.write('%10d\n' % (nAtoms))
191 OutFile.write('\n')
192
193 for i in range(nAtoms):
194
195 dx = positions[i][0] - COM[0]
196 dy = positions[i][1] - COM[1]
197 dz = positions[i][2] - COM[2]
198
199 r = numpy.array([dx,dy,dz])
200 rnew = numpy.dot(RotMat, r)
201
202 OutFile.write('%s\t%f\t%f\t%f\t%d\n' % (atypes[i], rnew[0], rnew[1], rnew[2], i))
203 OutFile.close()
204
205 def main(argv):
206 try:
207 opts, args = getopt.getopt(argv, "hx:o:", ["help", "xyz=", "out="])
208 except getopt.GetoptError:
209 usage()
210 sys.exit(2)
211 for opt, arg in opts:
212 if opt in ("-h", "--help"):
213 usage()
214 sys.exit()
215 elif opt in ("-x", "--xyz"):
216 XYZFileName = arg
217 global _haveXYZFileName
218 _haveXYZFileName = 1
219 elif opt in ("-o", "--out"):
220 OutFileName = arg
221 global _haveOutFileName
222 _haveOutFileName = 1
223
224
225 if (_haveXYZFileName != 1):
226 usage()
227 print "No xyz file was specified"
228 sys.exit()
229
230 readFile(XYZFileName)
231
232 if (_haveOutFileName == 1):
233 writeFile(OutFileName)
234 else:
235 findMoments()
236
237 if __name__ == "__main__":
238 if len(sys.argv) == 1:
239 usage()
240 sys.exit()
241 main(sys.argv[1:])