OpenMD 3.1
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
principalAxisCalculator
1#!@Python3_EXECUTABLE@
2"""principalAxisCalculator
3
4Opens an XYZ file and computes the moments of inertia and principal
5axes for the structure in the XYZ file. Optionally rotates the
6structure so that the long axis (that with the smallest eigenvalue) is
7pointing along the z-axis.
8
9Usage: principalAxisCalculator
10
11Options:
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
17Example:
18 principalAxisCalculator -x junk.xyz -o rot.xyz
19
20"""
21
22__author__ = "Dan Gezelter (gezelter@nd.edu)"
23__copyright__ = "Copyright (c) 2004-present The University of Notre Dame. All Rights Reserved."
24__license__ = "OpenMD"
25
26import sys
27import getopt
28import string
29import math
30import random
31import numpy
32from functools import reduce
33
34
35_haveXYZFileName = 0
36_haveOutFileName = 0
37
38positions = []
39indices = []
40atypes = []
41
42mass_table = {
43 'H': 1.00794,
44 'C': 12.0107,
45 'Cl': 35.453,
46 'O': 15.999,
47 'N': 14.007,
48 'S': 32.0655,
49 'Au': 196.466569,
50 'CH4': 16.05,
51 'CH3': 15.04,
52 'CH2': 14.03,
53 'CH': 13.02,
54 'CHar': 13.02,
55 'CHa': 13.02,
56 'RCHar': 12.0107,
57 'RCH': 12.0107,
58 'CH3S': 15.04,
59 'CH2S': 14.03,
60 'CHS': 13.02,
61 'CS': 12.0107,
62 'SYZ': 32.0655,
63 'SH': 32.0655,
64 'HS': 1.0079,
65 'S': 32.0655,
66 'SZ': 32.0655,
67 'SS': 32.0655,
68 'SP': 32.0655,
69 'CS': 12.0107,
70 'SAu': 228.9807,
71 'SSD': 18.0153,
72 'SSD1': 18.0153,
73 'SSD_E': 18.0153,
74 'SSD_RF': 18.0153,
75 'O_TIP3P': 15.9994,
76 'O_TIP4P': 15.9994,
77 'O_TIP4P-Ew': 15.9994,
78 'O_TIP5P': 15.9994,
79 'O_TIP5P-E': 15.9994,
80 'O_SPCE': 15.9994,
81 'O_SPC': 15.9994,
82 'H_TIP3P': 1.0079,
83 'H_TIP4P': 1.0079,
84 'H_TIP4P-Ew': 1.0079,
85 'H_TIP5P': 1.0079,
86 'H_SPCE': 1.0079,
87 'H_SPC': 1.0079,
88 'EP_TIP4P': 0.0,
89 'EP_TIP4P-Ew':0.0,
90 'EP_TIP5P': 0.0,
91 'Ni': 58.710,
92 'Cu': 63.550,
93 'Rh': 102.90550,
94 'Pd': 106.42,
95 'Ag': 107.8682,
96 'Ir': 192.217,
97 'Pt': 195.09
98}
99
100
101def usage():
102 print(__doc__)
103
104def add(x, y):
105 return x+y
106
107def sum(seq):
108 return reduce(add, seq)
109
110def readFile(XYZFileName):
111 print("reading XYZ file")
112
113 XYZFile = open(XYZFileName, 'r')
114 # Find number of atoms first
115 line = XYZFile.readline()
116 L = line.split()
117 nAtoms = int(L[0])
118 # skip comment line
119 line = XYZFile.readline()
120 for i in range(nAtoms):
121 line = XYZFile.readline()
122 L = line.split()
123 myIndex = i
124 indices.append(myIndex)
125 atomType = L[0]
126 atypes.append(atomType)
127 x = float(L[1])
128 y = float(L[2])
129 z = float(L[3])
130 positions.append([x, y, z])
131 XYZFile.close()
132
133
134def findCOM():
135 #find center of mass
136 Xcom = 0.0
137 Ycom = 0.0
138 Zcom = 0.0
139 totalMass = 0.0
140
141 for i in range(0, len(indices)):
142 myMass = mass_table[atypes[i]]
143
144 Xcom = Xcom + myMass * positions[i][0]
145 Ycom = Ycom + myMass * positions[i][1]
146 Zcom = Zcom + myMass * positions[i][2]
147 totalMass = totalMass + myMass
148
149 Xcom = Xcom / totalMass
150 Ycom = Ycom / totalMass
151 Zcom = Zcom / totalMass
152
153 COM = [Xcom, Ycom, Zcom]
154
155 return COM
156
157def findMoments():
158
159 COM = findCOM()
160
161 #find inertia tensor matrix elements
162
163 I = numpy.zeros((3, 3), numpy.float)
164
165 for i in range(0, len(indices)):
166 myMass = mass_table[atypes[i]]
167
168 dx = positions[i][0] - COM[0]
169 dy = positions[i][1] - COM[1]
170 dz = positions[i][2] - COM[2]
171
172 I[0, 0] = I[0, 0] + myMass * ( dy * dy + dz * dz )
173 I[1, 1] = I[1, 1] + myMass * ( dx * dx + dz * dz )
174 I[2, 2] = I[2, 2] + myMass * ( dx * dx + dy * dy )
175
176 I[0, 1] = I[0, 1] - myMass * ( dx * dy )
177 I[0, 2] = I[0, 2] - myMass * ( dx * dz )
178
179 I[1, 2] = I[1, 2] - myMass * ( dy * dz )
180
181 I[1, 0] = I[0, 1]
182 I[2, 0] = I[0, 2]
183 I[2, 1] = I[1, 2]
184
185 print("Inertia Tensor:")
186 print(I)
187 print()
188
189 (evals, evects) = numpy.linalg.eig(I)
190 print("evals:")
191 print(evals)
192 print()
193 print("evects:")
194 print(evects)
195 print()
196
197 return (COM, evals, evects)
198
199def writeFile(OutFileName):
200
201 (COM, evals, evects) = findMoments()
202
203 # we need to re-order the axes so that the smallest moment of inertia
204 # (which corresponds to the long axis of the molecule) is along the z-axis
205 # we'll just reverse the order of the three axes:
206
207 axOrder = numpy.argsort(evals)
208 RotMat = numpy.zeros((3, 3), numpy.float)
209 RotMat[0] = evects[axOrder[2]]
210 RotMat[1] = evects[axOrder[1]]
211 RotMat[2] = evects[axOrder[0]]
212
213 q = [0.0, 0.0, 0.0, 0.0]
214 myEuler = [0.0, 0.0, 0.0]
215
216 # RotMat to Quat code is out of OpenMD's SquareMatrix3.hpp code:
217
218 t = RotMat[0][0] + RotMat[1][1] + RotMat[2][2] + 1.0
219
220 if( t > 1e-6 ):
221 s = 0.5 / math.sqrt( t )
222 q[0] = 0.25 / s
223 q[1] = (RotMat[1][2] - RotMat[2][1]) * s
224 q[2] = (RotMat[2][0] - RotMat[0][2]) * s
225 q[3] = (RotMat[0][1] - RotMat[1][0]) * s
226 else:
227 ad1 = RotMat[0][0]
228 ad2 = RotMat[1][1]
229 ad3 = RotMat[2][2]
230
231 if( ad1 >= ad2 and ad1 >= ad3 ):
232 s = 0.5 / math.sqrt( 1.0 + RotMat[0][0] - RotMat[1][1] - RotMat[2][2] )
233 q[0] = (RotMat[1][2] - RotMat[2][1]) * s
234 q[1] = 0.25 / s
235 q[2] = (RotMat[0][1] + RotMat[1][0]) * s
236 q[3] = (RotMat[0][2] + RotMat[2][0]) * s
237 elif ( ad2 >= ad1 and ad2 >= ad3 ):
238 s = 0.5 / math.sqrt( 1.0 + RotMat[1][1] - RotMat[0][0] - RotMat[2][2] )
239 q[0] = (RotMat[2][0] - RotMat[0][2] ) * s
240 q[1] = (RotMat[0][1] + RotMat[1][0]) * s
241 q[2] = 0.25 / s
242 q[3] = (RotMat[1][2] + RotMat[2][1]) * s
243 else:
244 s = 0.5 / math.sqrt( 1.0 + RotMat[2][2] - RotMat[0][0] - RotMat[1][1] )
245 q[0] = (RotMat[0][1] - RotMat[1][0]) * s
246 q[1] = (RotMat[0][2] + RotMat[2][0]) * s
247 q[2] = (RotMat[1][2] + RotMat[2][1]) * s
248
249 print("Quaternions:")
250 print(q)
251
252 theta = math.acos(RotMat[2][2])
253 ctheta = RotMat[2][2]
254 stheta = math.sqrt(1.0 - ctheta * ctheta)
255
256 if (math.fabs(stheta) < 1e-6):
257 psi = 0.0
258 phi = math.atan2(-RotMat[1][0], RotMat[0][0])
259 else:
260 phi = math.atan2(RotMat[2][0], -RotMat[2][1])
261 psi = math.atan2(RotMat[0][2], RotMat[1][2])
262
263 if (phi < 0):
264 phi = phi + 2.0 * math.pi;
265
266 if (psi < 0):
267 psi = psi + 2.0 * math.pi;
268
269 myEuler[0] = phi * 180.0 / math.pi;
270 myEuler[1] = theta * 180.0 / math.pi;
271 myEuler[2] = psi * 180.0 / math.pi;
272
273 print("Euler Angles:")
274 print(myEuler)
275
276 nAtoms = len(indices)
277
278 print("writing output XYZ file")
279
280 OutFile = open(OutFileName, 'w')
281
282 OutFile.write('%10d\n' % (nAtoms))
283 OutFile.write('\n')
284
285 for i in range(nAtoms):
286
287 dx = positions[i][0] - COM[0]
288 dy = positions[i][1] - COM[1]
289 dz = positions[i][2] - COM[2]
290
291 r = numpy.array([dx, dy, dz])
292 rnew = numpy.dot(RotMat, r)
293
294 OutFile.write('%s\t%f\t%f\t%f\t%d\n' % (atypes[i], rnew[0], rnew[1], rnew[2], i))
295 OutFile.close()
296
297def main(argv):
298 try:
299 opts, args = getopt.getopt(argv, "hx:o:", ["help", "xyz=", "out="])
300 except getopt.GetoptError:
301 usage()
302 sys.exit(2)
303 for opt, arg in opts:
304 if opt in ("-h", "--help"):
305 usage()
306 sys.exit()
307 elif opt in ("-x", "--xyz"):
308 XYZFileName = arg
309 global _haveXYZFileName
310 _haveXYZFileName = 1
311 elif opt in ("-o", "--out"):
312 OutFileName = arg
313 global _haveOutFileName
314 _haveOutFileName = 1
315
316
317 if (_haveXYZFileName != 1):
318 usage()
319 print("No xyz file was specified")
320 sys.exit()
321
322 readFile(XYZFileName)
323
324 if (_haveOutFileName == 1):
325 writeFile(OutFileName)
326 else:
327 findMoments()
328
329if __name__ == "__main__":
330 if len(sys.argv) == 1:
331 usage()
332 sys.exit()
333 main(sys.argv[1:])