| 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 |
|
} |
