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 1292 by chrisfen, Wed Jun 23 23:19:43 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 175 | Line 174 | void RigidBody::calcRefCoords( ) {
174  
175   void RigidBody::calcRefCoords( ) {
176  
177 <  int i,j,k, it, n_linear_coords;
177 >  int i, j, it, n_linear_coords, pAxis, maxAxis, midAxis;
178    double mtmp;
179    vec3 apos;
180    double refCOM[3];
181    vec3 ptmp;
182    double Itmp[3][3];
183 +  double pAxisMat[3][3], pAxisRotMat[3][3];
184    double evals[3];
185 –  double evects[3][3];
185    double r, r2, len;
186 <
186 >  double iMat[3][3];
187 >  double test[3];
188 >  
189    // First, find the center of mass:
190    
191    mass = 0.0;
# Line 196 | Line 197 | void RigidBody::calcRefCoords( ) {
197      mass += mtmp;
198  
199      apos = refCoords[i];
199 –    
200      for(j = 0; j < 3; j++) {
201        refCOM[j] += apos[j]*mtmp;    
202      }    
# Line 260 | Line 260 | void RigidBody::calcRefCoords( ) {
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"
# Line 269 | Line 269 | void RigidBody::calcRefCoords( ) {
269                 );
270                 exit(-1);    
271    }
272 <  
272 >
273    // renormalize column vectors:
274    
275    for (i=0; i < 3; i++) {
# Line 280 | Line 280 | void RigidBody::calcRefCoords( ) {
280      len = sqrt(len);
281      for (j = 0; j < 3; j++) {
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 +  
322 +  //calculate the proper rotation matrix
323 +  transposeMat3(pAxisMat, pAxisRotMat);
324 +  
325 +  
326 +  //rotate the rigid body to the principle axis frame
327 +  for (i = 0; i < myAtoms.size(); i++) {
328 +    matVecMul3(pAxisRotMat, refCoords[i].vec, refCoords[i].vec);
329 +    myAtoms[i]->setPos(refCoords[i].vec);
330 +  }
331 +  
332 +  identityMat3(iMat);
333 +  setA(iMat);
334 +  
335 +  //and resort the moments of intertia to match the new orientation
336 +  for (i=0; i<3; i++)
337 +    if (evals[i]<momIntTol)
338 +      evals[i] = 0.0;
339 +  I[0][0] = evals[maxAxis];
340 +  I[1][1] = evals[midAxis];
341 +  I[2][2] = evals[pAxis];
342   }
343  
344 < void RigidBody::doEulerToRotMat(vec3 &euler, mat3x3 &myA ){
344 > void RigidBody::doEulerToRotMat(double euler[3], double myA[3][3] ){
345  
346    double phi, theta, psi;
347    
# Line 344 | Line 401 | void RigidBody::getEulerAngles(double myEuler[3]) {
401    
402    
403    double phi,theta,psi,eps;
404 <  double pi;
405 <  double cphi,ctheta,cpsi;
406 <  double sphi,stheta,spsi;
350 <  double b[3];
351 <  int flip[3];
352 <  
404 >  double ctheta;
405 >  double stheta;
406 >    
407    // set the tolerance for Euler angles and rotation elements
408    
409    eps = 1.0e-8;
# Line 400 | Line 454 | double RigidBody::min(double x, double  y) {  
454    return (x > y) ? y : x;
455   }
456  
457 + double RigidBody::findMaxExtent(){
458 +        int i;
459 +        double refAtomPos[3];
460 +        double maxExtent;
461 +        double tempExtent;
462 +        
463 +        //zero the extent variables
464 +        maxExtent = 0.0;
465 +        tempExtent = 0.0;
466 +        for (i=0; i<3; i++)
467 +                refAtomPos[i] = 0.0;
468 +        
469 +        //loop over all atoms
470 +        for (i=0; i<myAtoms.size(); i++){
471 +                getAtomRefCoor(refAtomPos, i);
472 +                tempExtent = sqrt(refAtomPos[0]*refAtomPos[0] + refAtomPos[1]*refAtomPos[1]
473 +                                                  + refAtomPos[2]*refAtomPos[2]);
474 +                if (tempExtent > maxExtent)
475 +                        maxExtent = tempExtent;
476 +        }
477 +        return maxExtent;
478 + }
479 +
480   void RigidBody::findCOM() {
481    
482    size_t i;
483    int j;
484    double mtmp;
485    double ptmp[3];
486 <  double vtmp[3];
410 <  
486 >  
487    for(j = 0; j < 3; j++) {
488      pos[j] = 0.0;
489    }
# Line 457 | Line 533 | void RigidBody::getAtomRefCoor(double pos[3], int inde
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 + }
548 +
549 + char *RigidBody::getAtomBase(int index){
550 +  
551 +  return myAtoms[index]->getBase();
552 +  
553 + }

Diff Legend

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