| 1 |
|
#include <math.h> |
| 2 |
+ |
#include <iostream> |
| 3 |
|
#include "RigidBody.hpp" |
| 4 |
|
#include "VDWAtom.hpp" |
| 5 |
|
#include "MatVec3.h" |
| 182 |
|
double Itmp[3][3]; |
| 183 |
|
double pAxisMat[3][3], pAxisRotMat[3][3]; |
| 184 |
|
double evals[3]; |
| 184 |
– |
double prePos[3], rotPos[3]; |
| 185 |
|
double r, r2, len; |
| 186 |
|
double iMat[3][3]; |
| 187 |
< |
|
| 187 |
> |
double test[3]; |
| 188 |
> |
|
| 189 |
|
// First, find the center of mass: |
| 190 |
|
|
| 191 |
|
mass = 0.0; |
| 197 |
|
mass += mtmp; |
| 198 |
|
|
| 199 |
|
apos = refCoords[i]; |
| 199 |
– |
|
| 200 |
|
for(j = 0; j < 3; j++) { |
| 201 |
|
refCOM[j] += apos[j]*mtmp; |
| 202 |
|
} |
| 260 |
|
if (n_linear_coords > 1) { |
| 261 |
|
printf( |
| 262 |
|
"RigidBody error.\n" |
| 263 |
< |
"\tOOPSE found more than one axis in this rigid body with a vanishing \n" |
| 263 |
> |
"\tSHAPES found more than one axis in this rigid body with a vanishing \n" |
| 264 |
|
"\tmoment of inertia. This can happen in one of three ways:\n" |
| 265 |
|
"\t 1) Only one atom was specified, or \n" |
| 266 |
|
"\t 2) All atoms were specified at the same location, or\n" |
| 269 |
|
); |
| 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 |
< |
|
| 272 |
> |
|
| 273 |
|
// renormalize column vectors: |
| 274 |
|
|
| 275 |
|
for (i=0; i < 3; i++) { |
| 282 |
|
sU[i][j] /= len; |
| 283 |
|
} |
| 284 |
|
} |
| 285 |
+ |
|
| 286 |
+ |
//sort and reorder the moment axes |
| 287 |
+ |
|
| 288 |
+ |
// The only problem below is for molecules like C60 with 3 nearly identical |
| 289 |
+ |
// non-zero moments of inertia. In this case it doesn't really matter which is |
| 290 |
+ |
// the principal axis, so they get assigned nearly randomly depending on the |
| 291 |
+ |
// floating point comparison between eigenvalues |
| 292 |
+ |
if (! is_linear) { |
| 293 |
+ |
pAxis = 0; |
| 294 |
+ |
maxAxis = 0; |
| 295 |
+ |
|
| 296 |
+ |
for (i = 0; i < 3; i++) { |
| 297 |
+ |
if (evals[i] < evals[pAxis]) pAxis = i; |
| 298 |
+ |
if (evals[i] > evals[maxAxis]) maxAxis = i; |
| 299 |
+ |
} |
| 300 |
+ |
|
| 301 |
+ |
midAxis = 0; |
| 302 |
+ |
for (i=0; i < 3; i++) { |
| 303 |
+ |
if (pAxis != i && maxAxis != i) midAxis = i; |
| 304 |
+ |
} |
| 305 |
+ |
} else { |
| 306 |
+ |
pAxis = linear_axis; |
| 307 |
+ |
// linear molecules have one zero moment of inertia and two identical |
| 308 |
+ |
// moments of inertia. In this case, it doesn't matter which is chosen |
| 309 |
+ |
// as mid and which is max, so just permute from the pAxis: |
| 310 |
+ |
midAxis = (pAxis + 1)%3; |
| 311 |
+ |
maxAxis = (pAxis + 2)%3; |
| 312 |
+ |
} |
| 313 |
+ |
|
| 314 |
+ |
//let z be the smallest and x be the largest eigenvalue axes |
| 315 |
+ |
for (i=0; i<3; i++){ |
| 316 |
+ |
pAxisMat[i][2] = sU[i][pAxis]; |
| 317 |
+ |
pAxisMat[i][1] = sU[i][midAxis]; |
| 318 |
+ |
pAxisMat[i][0] = sU[i][maxAxis]; |
| 319 |
+ |
} |
| 320 |
+ |
|
| 321 |
+ |
//calculate the proper rotation matrix |
| 322 |
+ |
transposeMat3(pAxisMat, pAxisRotMat); |
| 323 |
+ |
|
| 324 |
+ |
for (i=0; i<myAtoms.size(); i++){ |
| 325 |
+ |
getAtomPos(test, i); |
| 326 |
+ |
printf("%d\t%d\t%d\n",test[0],test[1],test[2]); |
| 327 |
+ |
} |
| 328 |
+ |
|
| 329 |
+ |
//rotate the rigid body to the principle axis frame |
| 330 |
+ |
for (i = 0; i < myAtoms.size(); i++) { |
| 331 |
+ |
matVecMul3(pAxisRotMat, refCoords[i].vec, refCoords[i].vec); |
| 332 |
+ |
myAtoms[i]->setPos(refCoords[i].vec); |
| 333 |
+ |
} |
| 334 |
+ |
|
| 335 |
+ |
for (i=0; i<myAtoms.size(); i++){ |
| 336 |
+ |
getAtomPos(test,i); |
| 337 |
+ |
printf("%d\t%d\t%d\n",test[0],test[1],test[2]); |
| 338 |
+ |
} |
| 339 |
+ |
|
| 340 |
+ |
identityMat3(iMat); |
| 341 |
+ |
setA(iMat); |
| 342 |
|
} |
| 343 |
|
|
| 344 |
|
void RigidBody::doEulerToRotMat(double euler[3], double myA[3][3] ){ |
| 533 |
|
pos[2] = ref[2]; |
| 534 |
|
|
| 535 |
|
} |
| 536 |
+ |
|
| 537 |
+ |
double RigidBody::getAtomRpar(int index){ |
| 538 |
+ |
|
| 539 |
+ |
return myAtoms[index]->getRpar(); |
| 540 |
+ |
|
| 541 |
+ |
} |
| 542 |
+ |
|
| 543 |
+ |
double RigidBody::getAtomEps(int index){ |
| 544 |
+ |
|
| 545 |
+ |
return myAtoms[index]->getEps(); |
| 546 |
+ |
|
| 547 |
+ |
} |
