15 |
|
void RigidBody::addAtom(VDWAtom* at) { |
16 |
|
|
17 |
|
vec3 coords; |
18 |
– |
vec3 euler; |
19 |
– |
mat3x3 Atmp; |
18 |
|
|
19 |
|
myAtoms.push_back(at); |
20 |
|
|
173 |
|
|
174 |
|
void RigidBody::calcRefCoords( ) { |
175 |
|
|
176 |
< |
int i,j,k, it, n_linear_coords; |
176 |
> |
int i, j, it, n_linear_coords, pAxis, maxAxis, midAxis; |
177 |
|
double mtmp; |
178 |
|
vec3 apos; |
179 |
|
double refCOM[3]; |
180 |
|
vec3 ptmp; |
181 |
|
double Itmp[3][3]; |
182 |
+ |
double pAxisMat[3][3], pAxisRotMat[3][3]; |
183 |
|
double evals[3]; |
184 |
< |
double evects[3][3]; |
184 |
> |
double prePos[3], rotPos[3]; |
185 |
|
double r, r2, len; |
186 |
+ |
double iMat[3][3]; |
187 |
|
|
188 |
|
// First, find the center of mass: |
189 |
|
|
270 |
|
exit(-1); |
271 |
|
} |
272 |
|
|
273 |
+ |
//sort and reorder the moment axes |
274 |
+ |
if (evals[0] < evals[1] && evals[0] < evals[2]) |
275 |
+ |
pAxis = 0; |
276 |
+ |
else if (evals[1] < evals[2]) |
277 |
+ |
pAxis = 1; |
278 |
+ |
else |
279 |
+ |
pAxis = 2; |
280 |
+ |
|
281 |
+ |
if (evals[0] > evals[1] && evals[0] > evals[2]) |
282 |
+ |
maxAxis = 0; |
283 |
+ |
else if (evals[1] > evals[2]) |
284 |
+ |
maxAxis = 1; |
285 |
+ |
else |
286 |
+ |
maxAxis = 2; |
287 |
+ |
|
288 |
+ |
midAxis = 0; |
289 |
+ |
if (midAxis == pAxis || midAxis == pAxis) |
290 |
+ |
midAxis = 1; |
291 |
+ |
if (midAxis == pAxis || midAxis == pAxis) |
292 |
+ |
midAxis = 2; |
293 |
+ |
|
294 |
+ |
if (pAxis != maxAxis){ |
295 |
+ |
//zero out our matrices |
296 |
+ |
for (i=0; i<3; i++){ |
297 |
+ |
for (j=0; j<3; j++) { |
298 |
+ |
pAxisMat[i][j] = 0.0; |
299 |
+ |
pAxisRotMat[i][j] = 0.0; |
300 |
+ |
} |
301 |
+ |
} |
302 |
+ |
|
303 |
+ |
//let z be the smallest and x be the largest eigenvalue axes |
304 |
+ |
for (i=0; i<3; i++){ |
305 |
+ |
pAxisMat[i][2] = I[i][pAxis]; |
306 |
+ |
pAxisMat[i][1] = I[i][midAxis]; |
307 |
+ |
pAxisMat[i][0] = I[i][maxAxis]; |
308 |
+ |
} |
309 |
+ |
|
310 |
+ |
//calculate the proper rotation matrix |
311 |
+ |
transposeMat3(pAxisMat, pAxisRotMat); |
312 |
+ |
|
313 |
+ |
//rotate the rigid body to the principle axis frame |
314 |
+ |
for (i = 0; i < myAtoms.size(); i++) { |
315 |
+ |
apos = refCoords[i]; |
316 |
+ |
for (j=0; j<3; j++) |
317 |
+ |
prePos[j] = apos[j]; |
318 |
+ |
|
319 |
+ |
matVecMul3(pAxisRotMat, prePos, rotPos); |
320 |
+ |
|
321 |
+ |
for (j=0; j < 3; j++) |
322 |
+ |
apos[j] = rotPos[j]; |
323 |
+ |
|
324 |
+ |
refCoords[i] = apos; |
325 |
+ |
} |
326 |
+ |
|
327 |
+ |
//the lab and the body frame match up at this point, so A = Identity Matrix |
328 |
+ |
for (i=0; i<3; i++){ |
329 |
+ |
for (j=0; j<3; j++){ |
330 |
+ |
if (i == j) |
331 |
+ |
iMat[i][j] = 1.0; |
332 |
+ |
else |
333 |
+ |
iMat[i][j] = 0.0; |
334 |
+ |
} |
335 |
+ |
} |
336 |
+ |
setA(iMat); |
337 |
+ |
} |
338 |
+ |
|
339 |
|
// renormalize column vectors: |
340 |
|
|
341 |
|
for (i=0; i < 3; i++) { |
350 |
|
} |
351 |
|
} |
352 |
|
|
353 |
< |
void RigidBody::doEulerToRotMat(vec3 &euler, mat3x3 &myA ){ |
353 |
> |
void RigidBody::doEulerToRotMat(double euler[3], double myA[3][3] ){ |
354 |
|
|
355 |
|
double phi, theta, psi; |
356 |
|
|
410 |
|
|
411 |
|
|
412 |
|
double phi,theta,psi,eps; |
413 |
< |
double pi; |
414 |
< |
double cphi,ctheta,cpsi; |
415 |
< |
double sphi,stheta,spsi; |
350 |
< |
double b[3]; |
351 |
< |
int flip[3]; |
352 |
< |
|
413 |
> |
double ctheta; |
414 |
> |
double stheta; |
415 |
> |
|
416 |
|
// set the tolerance for Euler angles and rotation elements |
417 |
|
|
418 |
|
eps = 1.0e-8; |
463 |
|
return (x > y) ? y : x; |
464 |
|
} |
465 |
|
|
466 |
+ |
double RigidBody::findMaxExtent(){ |
467 |
+ |
int i; |
468 |
+ |
double refAtomPos[3]; |
469 |
+ |
double maxExtent; |
470 |
+ |
double tempExtent; |
471 |
+ |
|
472 |
+ |
//zero the extent variables |
473 |
+ |
maxExtent = 0.0; |
474 |
+ |
tempExtent = 0.0; |
475 |
+ |
for (i=0; i<3; i++) |
476 |
+ |
refAtomPos[i] = 0.0; |
477 |
+ |
|
478 |
+ |
//loop over all atoms |
479 |
+ |
for (i=0; i<myAtoms.size(); i++){ |
480 |
+ |
getAtomRefCoor(refAtomPos, i); |
481 |
+ |
tempExtent = sqrt(refAtomPos[0]*refAtomPos[0] + refAtomPos[1]*refAtomPos[1] |
482 |
+ |
+ refAtomPos[2]*refAtomPos[2]); |
483 |
+ |
if (tempExtent > maxExtent) |
484 |
+ |
maxExtent = tempExtent; |
485 |
+ |
} |
486 |
+ |
return maxExtent; |
487 |
+ |
} |
488 |
+ |
|
489 |
|
void RigidBody::findCOM() { |
490 |
|
|
491 |
|
size_t i; |
492 |
|
int j; |
493 |
|
double mtmp; |
494 |
|
double ptmp[3]; |
495 |
< |
double vtmp[3]; |
410 |
< |
|
495 |
> |
|
496 |
|
for(j = 0; j < 3; j++) { |
497 |
|
pos[j] = 0.0; |
498 |
|
} |