ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/SHAPES/RigidBody.cpp
(Generate patch)

Comparing trunk/SHAPES/RigidBody.cpp (file contents):
Revision 1271 by gezelter, Tue Jun 15 20:20:36 2004 UTC vs.
Revision 1283 by gezelter, Mon Jun 21 15:54:27 2004 UTC

# Line 1 | Line 1
1   #include <math.h>
2 + #include <iostream>
3   #include "RigidBody.hpp"
4   #include "VDWAtom.hpp"
5   #include "MatVec3.h"
# Line 15 | Line 16 | void RigidBody::addAtom(VDWAtom* at) {
16   void RigidBody::addAtom(VDWAtom* at) {
17  
18    vec3 coords;
18 –  vec3 euler;
19 –  mat3x3 Atmp;
19  
20    myAtoms.push_back(at);
21    
# Line 49 | Line 48 | void RigidBody::setEuler( double phi, double theta, do
48      A[2][1] = -cos(phi) * sin(theta);
49      A[2][2] = cos(theta);
50  
51 +    printf("A[2][x] = %lf\t%lf\t%lf\n", A[2][0], A[2][1], A[2][2]);
52 +
53   }
54  
55   void RigidBody::getQ( double q[4] ){
# Line 175 | Line 176 | void RigidBody::calcRefCoords( ) {
176  
177   void RigidBody::calcRefCoords( ) {
178  
179 <  int i,j,k, it, n_linear_coords;
179 >  int i, j, it, n_linear_coords, pAxis, maxAxis, midAxis;
180    double mtmp;
181    vec3 apos;
182    double refCOM[3];
183    vec3 ptmp;
184    double Itmp[3][3];
185 <  double evals[3];
186 <  double evects[3][3];
185 >  double pAxisMat[3][3], pAxisRotMat[3][3];
186 >  double evals[3];
187    double r, r2, len;
188 <
188 >  double iMat[3][3];
189 >  double test[3];
190 >  
191    // First, find the center of mass:
192    
193    mass = 0.0;
# Line 196 | Line 199 | void RigidBody::calcRefCoords( ) {
199      mass += mtmp;
200  
201      apos = refCoords[i];
199 –    
202      for(j = 0; j < 3; j++) {
203        refCOM[j] += apos[j]*mtmp;    
204      }    
# Line 260 | Line 262 | void RigidBody::calcRefCoords( ) {
262    if (n_linear_coords > 1) {
263            printf(
264                 "RigidBody error.\n"
265 <               "\tOOPSE found more than one axis in this rigid body with a vanishing \n"
265 >               "\tSHAPES found more than one axis in this rigid body with a vanishing \n"
266                 "\tmoment of inertia.  This can happen in one of three ways:\n"
267                 "\t 1) Only one atom was specified, or \n"
268                 "\t 2) All atoms were specified at the same location, or\n"
# Line 269 | Line 271 | void RigidBody::calcRefCoords( ) {
271                 );
272                 exit(-1);    
273    }
274 <  
274 >
275    // renormalize column vectors:
276    
277    for (i=0; i < 3; i++) {
# Line 280 | Line 282 | void RigidBody::calcRefCoords( ) {
282      len = sqrt(len);
283      for (j = 0; j < 3; j++) {
284        sU[i][j] /= len;
285 +    }
286 +  }
287 +  
288 +  //sort and reorder the moment axes
289 +
290 +  // The only problem below is for molecules like C60 with 3 nearly identical
291 +  // non-zero moments of inertia.  In this case it doesn't really matter which is
292 +  // the principal axis, so they get assigned nearly randomly depending on the
293 +  // floating point comparison between eigenvalues
294 +  if (! is_linear) {
295 +    pAxis = 0;
296 +    maxAxis = 0;
297 +  
298 +    for (i = 0; i < 3; i++) {
299 +      if (evals[i] < evals[pAxis]) pAxis = i;
300 +      if (evals[i] > evals[maxAxis]) maxAxis = i;
301      }
302 +  
303 +    midAxis = 0;
304 +    for (i=0; i < 3; i++) {
305 +      if (pAxis != i && maxAxis != i) midAxis = i;
306 +    }
307 +  } else {
308 +    pAxis = linear_axis;
309 +    // linear molecules have one zero moment of inertia and two identical
310 +    // moments of inertia.  In this case, it doesn't matter which is chosen
311 +    // as mid and which is max, so just permute from the pAxis:
312 +    midAxis = (pAxis + 1)%3;
313 +    maxAxis = (pAxis + 2)%3;
314    }
315 +      
316 +  //let z be the smallest and x be the largest eigenvalue axes
317 +  for (i=0; i<3; i++){
318 +    pAxisMat[i][2] = sU[i][pAxis];
319 +    pAxisMat[i][1] = sU[i][midAxis];
320 +    pAxisMat[i][0] = sU[i][maxAxis];
321 +  }
322 +    
323 +  //calculate the proper rotation matrix
324 +  transposeMat3(pAxisMat, pAxisRotMat);
325 +  
326 +  
327 +  for (i=0; i<myAtoms.size(); i++){
328 +    apos = refCoords[i];
329 +    printf("%f\t%f\t%f\n",apos[0],apos[1],apos[2]);
330 +  }
331 +    
332 +  //rotate the rigid body to the principle axis frame
333 +  for (i = 0; i < myAtoms.size(); i++) {
334 +    matVecMul3(pAxisRotMat, refCoords[i].vec, refCoords[i].vec);
335 +    myAtoms[i]->setPos(refCoords[i].vec);
336 +  }
337 +  
338 +  for (i=0; i<myAtoms.size(); i++){
339 +    apos = refCoords[i];
340 +    printf("%f\t%f\t%f\n",apos[0],apos[1],apos[2]);
341 +  }
342 +  
343 +  identityMat3(iMat);
344 +  setA(iMat);
345   }
346  
347 < void RigidBody::doEulerToRotMat(vec3 &euler, mat3x3 &myA ){
347 > void RigidBody::doEulerToRotMat(double euler[3], double myA[3][3] ){
348  
349    double phi, theta, psi;
350    
# Line 344 | Line 404 | void RigidBody::getEulerAngles(double myEuler[3]) {
404    
405    
406    double phi,theta,psi,eps;
407 <  double pi;
408 <  double cphi,ctheta,cpsi;
409 <  double sphi,stheta,spsi;
350 <  double b[3];
351 <  int flip[3];
352 <  
407 >  double ctheta;
408 >  double stheta;
409 >    
410    // set the tolerance for Euler angles and rotation elements
411    
412    eps = 1.0e-8;
# Line 400 | Line 457 | double RigidBody::min(double x, double  y) {  
457    return (x > y) ? y : x;
458   }
459  
460 + double RigidBody::findMaxExtent(){
461 +        int i;
462 +        double refAtomPos[3];
463 +        double maxExtent;
464 +        double tempExtent;
465 +        
466 +        //zero the extent variables
467 +        maxExtent = 0.0;
468 +        tempExtent = 0.0;
469 +        for (i=0; i<3; i++)
470 +                refAtomPos[i] = 0.0;
471 +        
472 +        //loop over all atoms
473 +        for (i=0; i<myAtoms.size(); i++){
474 +                getAtomRefCoor(refAtomPos, i);
475 +                tempExtent = sqrt(refAtomPos[0]*refAtomPos[0] + refAtomPos[1]*refAtomPos[1]
476 +                                                  + refAtomPos[2]*refAtomPos[2]);
477 +                if (tempExtent > maxExtent)
478 +                        maxExtent = tempExtent;
479 +        }
480 +        return maxExtent;
481 + }
482 +
483   void RigidBody::findCOM() {
484    
485    size_t i;
486    int j;
487    double mtmp;
488    double ptmp[3];
489 <  double vtmp[3];
410 <  
489 >  
490    for(j = 0; j < 3; j++) {
491      pos[j] = 0.0;
492    }
# Line 457 | Line 536 | void RigidBody::getAtomRefCoor(double pos[3], int inde
536    pos[2] = ref[2];
537    
538   }
539 +
540 + double RigidBody::getAtomRpar(int index){
541 +  
542 +  return myAtoms[index]->getRpar();
543 +  
544 + }
545 +
546 + double RigidBody::getAtomEps(int index){
547 +  
548 +  return myAtoms[index]->getEps();
549 +  
550 + }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines