| 185 |
|
|
| 186 |
|
(COM, evals, evects) = findMoments() |
| 187 |
|
|
| 188 |
– |
|
| 188 |
|
# we need to re-order the axes so that the smallest moment of inertia |
| 189 |
|
# (which corresponds to the long axis of the molecule) is along the z-axis |
| 190 |
|
# we'll just reverse the order of the three axes: |
| 194 |
|
RotMat[0] = evects[axOrder[2]] |
| 195 |
|
RotMat[1] = evects[axOrder[1]] |
| 196 |
|
RotMat[2] = evects[axOrder[0]] |
| 197 |
+ |
|
| 198 |
+ |
q = [0.0, 0.0, 0.0, 0.0] |
| 199 |
+ |
myEuler = [0.0, 0.0, 0.0] |
| 200 |
+ |
|
| 201 |
+ |
# RotMat to Quat code is out of OpenMD's SquareMatrix3.hpp code: |
| 202 |
+ |
|
| 203 |
+ |
t = RotMat[0][0] + RotMat[1][1] + RotMat[2][2] + 1.0 |
| 204 |
+ |
|
| 205 |
+ |
if( t > 1e-6 ): |
| 206 |
+ |
s = 0.5 / math.sqrt( t ) |
| 207 |
+ |
q[0] = 0.25 / s |
| 208 |
+ |
q[1] = (RotMat[1][2] - RotMat[2][1]) * s |
| 209 |
+ |
q[2] = (RotMat[2][0] - RotMat[0][2]) * s |
| 210 |
+ |
q[3] = (RotMat[0][1] - RotMat[1][0]) * s |
| 211 |
+ |
else: |
| 212 |
+ |
ad1 = RotMat[0][0] |
| 213 |
+ |
ad2 = RotMat[1][1] |
| 214 |
+ |
ad3 = RotMat[2][2] |
| 215 |
+ |
|
| 216 |
+ |
if( ad1 >= ad2 and ad1 >= ad3 ): |
| 217 |
+ |
s = 0.5 / math.sqrt( 1.0 + RotMat[0][0] - RotMat[1][1] - RotMat[2][2] ) |
| 218 |
+ |
q[0] = (RotMat[1][2] - RotMat[2][1]) * s |
| 219 |
+ |
q[1] = 0.25 / s |
| 220 |
+ |
q[2] = (RotMat[0][1] + RotMat[1][0]) * s |
| 221 |
+ |
q[3] = (RotMat[0][2] + RotMat[2][0]) * s |
| 222 |
+ |
elif ( ad2 >= ad1 and ad2 >= ad3 ): |
| 223 |
+ |
s = 0.5 / math.sqrt( 1.0 + RotMat[1][1] - RotMat[0][0] - RotMat[2][2] ) |
| 224 |
+ |
q[0] = (RotMat[2][0] - RotMat[0][2] ) * s |
| 225 |
+ |
q[1] = (RotMat[0][1] + RotMat[1][0]) * s |
| 226 |
+ |
q[2] = 0.25 / s |
| 227 |
+ |
q[3] = (RotMat[1][2] + RotMat[2][1]) * s |
| 228 |
+ |
else: |
| 229 |
+ |
s = 0.5 / math.sqrt( 1.0 + RotMat[2][2] - RotMat[0][0] - RotMat[1][1] ) |
| 230 |
+ |
q[0] = (RotMat[0][1] - RotMat[1][0]) * s |
| 231 |
+ |
q[1] = (RotMat[0][2] + RotMat[2][0]) * s |
| 232 |
+ |
q[2] = (RotMat[1][2] + RotMat[2][1]) * s |
| 233 |
|
|
| 234 |
+ |
print "Quaternions:" |
| 235 |
+ |
print q |
| 236 |
+ |
|
| 237 |
+ |
theta = math.acos(RotMat[2][2]) |
| 238 |
+ |
ctheta = RotMat[2][2] |
| 239 |
+ |
stheta = math.sqrt(1.0 - ctheta * ctheta) |
| 240 |
+ |
|
| 241 |
+ |
if (math.fabs(stheta) < 1e-6): |
| 242 |
+ |
psi = 0.0 |
| 243 |
+ |
phi = math.atan2(-RotMat[1][0], RotMat[0][0]) |
| 244 |
+ |
else: |
| 245 |
+ |
phi = math.atan2(RotMat[2][0], -RotMat[2][1]) |
| 246 |
+ |
psi = math.atan2(RotMat[0][2], RotMat[1][2]) |
| 247 |
+ |
|
| 248 |
+ |
if (phi < 0): |
| 249 |
+ |
phi = phi + 2.0 * math.pi; |
| 250 |
+ |
|
| 251 |
+ |
if (psi < 0): |
| 252 |
+ |
psi = psi + 2.0 * math.pi; |
| 253 |
+ |
|
| 254 |
+ |
myEuler[0] = phi * 180.0 / math.pi; |
| 255 |
+ |
myEuler[1] = theta * 180.0 / math.pi; |
| 256 |
+ |
myEuler[2] = psi * 180.0 / math.pi; |
| 257 |
+ |
|
| 258 |
+ |
print "Euler Angles:" |
| 259 |
+ |
print myEuler |
| 260 |
+ |
|
| 261 |
|
nAtoms = len(indices) |
| 262 |
|
|
| 263 |
|
print "writing output XYZ file" |
