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

Comparing trunk/OOPSE/libmdtools/DirectionalAtom.cpp (file contents):
Revision 878 by gezelter, Fri Dec 12 15:42:13 2003 UTC vs.
Revision 1097 by gezelter, Mon Apr 12 20:32:20 2004 UTC

# Line 1 | Line 1
1   #include <math.h>
2  
3   #include "Atom.hpp"
4 + #include "DirectionalAtom.hpp"
5   #include "simError.h"
6 + #include "MatVec3.h"
7  
6 –
7 –
8   void DirectionalAtom::zeroForces() {
9    if( hasCoords ){
10 <    frc[offsetX] = 0.0;
11 <    frc[offsetY] = 0.0;
12 <    frc[offsetZ] = 0.0;
10 >
11 >    Atom::zeroForces();
12      
13      trq[offsetX] = 0.0;
14      trq[offsetY] = 0.0;
# Line 48 | Line 47 | void DirectionalAtom::setCoords(void){
47  
48    hasCoords = true;
49  
51 –  *mu = myMu;
52 –
50   }
51  
55 – double DirectionalAtom::getMu( void ) {
56 –
57 –  if( hasCoords ){
58 –    return *mu;
59 –  }
60 –  else{
61 –    return myMu;
62 –  }
63 – }
64 –
65 – void DirectionalAtom::setMu( double the_mu ) {
66 –
67 –  if( hasCoords ){
68 –    *mu = the_mu;
69 –    myMu = the_mu;
70 –  }
71 –  else{
72 –    myMu = the_mu;
73 –  }
74 – }
75 –
52   void DirectionalAtom::setA( double the_A[3][3] ){
53  
54    if( hasCoords ){
# Line 92 | Line 68 | void DirectionalAtom::setI( double the_I[3][3] ){
68    }
69   }
70  
71 < void DirectionalAtom::setI( double the_I[3][3] ){
71 > void DirectionalAtom::setI( double the_I[3][3] ){  
72    
73    Ixx = the_I[0][0]; Ixy = the_I[0][1]; Ixz = the_I[0][2];
74    Iyx = the_I[1][0]; Iyy = the_I[1][1]; Iyz = the_I[1][2];
# Line 182 | Line 158 | void DirectionalAtom::getU( double the_u[3] ){
158  
159   void DirectionalAtom::getU( double the_u[3] ){
160    
161 <  the_u[0] = sux;
162 <  the_u[1] = suy;
163 <  the_u[2] = suz;
164 <
161 >  the_u[0] = sU[2][0];
162 >  the_u[1] = sU[2][1];
163 >  the_u[2] = sU[2][2];
164 >  
165    this->body2Lab( the_u );
166   }
167  
# Line 246 | Line 222 | void DirectionalAtom::getQ( double q[4] ){
222      simError();
223    }
224   }
225 +
226 + void DirectionalAtom::setUnitFrameFromEuler(double phi,
227 +                                            double theta,
228 +                                            double psi) {
229 +
230 +  double myA[3][3];
231 +  double uFrame[3][3];
232 +  double len;
233 +  int i, j;
234 +  
235 +  myA[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
236 +  myA[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
237 +  myA[0][2] = sin(theta) * sin(psi);
238 +  
239 +  myA[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
240 +  myA[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
241 +  myA[1][2] = sin(theta) * cos(psi);
242 +  
243 +  myA[2][0] = sin(phi) * sin(theta);
244 +  myA[2][1] = -cos(phi) * sin(theta);
245 +  myA[2][2] = cos(theta);
246 +  
247 +  // Make the unit Frame:
248 +
249 +  for (i=0; i < 3; i++)
250 +    for (j=0; j < 3; j++)
251 +      uFrame[i][j] = 0.0;
252 +
253 +  for (i=0; i < 3; i++)
254 +    uFrame[i][i] = 1.0;
255  
256 +  // rotate by the given rotation matrix:
257  
258 +  matMul3(myA, uFrame, sU);
259 +
260 +  // renormalize column vectors:
261 +
262 +  for (i=0; i < 3; i++) {
263 +    len = 0.0;
264 +    for (j = 0; j < 3; j++) {
265 +      len += sU[i][j]*sU[i][j];
266 +    }
267 +    len = sqrt(len);
268 +    for (j = 0; j < 3; j++) {
269 +      sU[i][j] /= len;    
270 +    }
271 +  }
272 +  
273 +  // sU now contains the coordinates of the 'special' frame;
274 +    
275 + }
276 +
277   void DirectionalAtom::setEuler( double phi, double theta, double psi ){
278    
279    if( hasCoords ){
# Line 300 | Line 326 | void DirectionalAtom::body2Lab( double r[3] ){
326  
327   }
328  
329 + void DirectionalAtom::rotateBy( double by_A[3][3]) {
330 +
331 +  // Check this
332 +  
333 +  double r00, r01, r02, r10, r11, r12, r20, r21, r22;
334 +
335 +  if( hasCoords ){
336 +
337 +    r00 = by_A[0][0]*Amat[Axx] + by_A[0][1]*Amat[Ayx] + by_A[0][2]*Amat[Azx];
338 +    r01 = by_A[0][0]*Amat[Axy] + by_A[0][1]*Amat[Ayy] + by_A[0][2]*Amat[Azy];
339 +    r02 = by_A[0][0]*Amat[Axz] + by_A[0][1]*Amat[Ayz] + by_A[0][2]*Amat[Azz];
340 +    
341 +    r10 = by_A[1][0]*Amat[Axx] + by_A[1][1]*Amat[Ayx] + by_A[1][2]*Amat[Azx];
342 +    r11 = by_A[1][0]*Amat[Axy] + by_A[1][1]*Amat[Ayy] + by_A[1][2]*Amat[Azy];
343 +    r12 = by_A[1][0]*Amat[Axz] + by_A[1][1]*Amat[Ayz] + by_A[1][2]*Amat[Azz];
344 +    
345 +    r20 = by_A[2][0]*Amat[Axx] + by_A[2][1]*Amat[Ayx] + by_A[2][2]*Amat[Azx];
346 +    r21 = by_A[2][0]*Amat[Axy] + by_A[2][1]*Amat[Ayy] + by_A[2][2]*Amat[Azy];
347 +    r22 = by_A[2][0]*Amat[Axz] + by_A[2][1]*Amat[Ayz] + by_A[2][2]*Amat[Azz];
348 +    
349 +    Amat[Axx] = r00; Amat[Axy] = r01; Amat[Axz] = r02;
350 +    Amat[Ayx] = r10; Amat[Ayy] = r11; Amat[Ayz] = r12;
351 +    Amat[Azx] = r20; Amat[Azy] = r21; Amat[Azz] = r22;
352 +
353 +  }
354 +  else{
355 +    
356 +    sprintf( painCave.errMsg,
357 +             "Attempt to rotate frame for atom %d before coords set.\n",
358 +             index );
359 +    painCave.isFatal = 1;
360 +    simError();
361 +  }
362 +
363 + }
364 +
365 +
366   void DirectionalAtom::body2Lab( double r[3] ){
367  
368    double rb[3]; // the body frame vector
# Line 326 | Line 389 | void DirectionalAtom::updateU( void ){
389   void DirectionalAtom::updateU( void ){
390  
391    if( hasCoords ){
392 <    ul[offsetX] = (Amat[Axx] * sux) + (Amat[Ayx] * suy) + (Amat[Azx] * suz);
393 <    ul[offsetY] = (Amat[Axy] * sux) + (Amat[Ayy] * suy) + (Amat[Azy] * suz);
394 <    ul[offsetZ] = (Amat[Axz] * sux) + (Amat[Ayz] * suy) + (Amat[Azz] * suz);
392 >    ul[offsetX] = (Amat[Axx] * sU[2][0]) +
393 >      (Amat[Ayx] * sU[2][1]) + (Amat[Azx] * sU[2][2]);
394 >    ul[offsetY] = (Amat[Axy] * sU[2][0]) +
395 >      (Amat[Ayy] * sU[2][1]) + (Amat[Azy] * sU[2][2]);
396 >    ul[offsetZ] = (Amat[Axz] * sU[2][0]) +
397 >      (Amat[Ayz] * sU[2][1]) + (Amat[Azz] * sU[2][2]);
398    }
399    else{
400      
# Line 429 | Line 495 | void DirectionalAtom::getGrad( double grad[6] ) {
495    ephi[0] = 0.0;
496    ephi[1] = 0.0;
497    ephi[2] = 1.0;
498 <  etheta[0] = -sphi;
499 <  etheta[1] = cphi;
498 >
499 >  etheta[0] = cphi;
500 >  etheta[1] = sphi;
501    etheta[2] = 0.0;
435 –  epsi[0] = ctheta * cphi;
436 –  epsi[1] = ctheta * sphi;
437 –  epsi[2] = -stheta;
502    
503 +  epsi[0] = stheta * cphi;
504 +  epsi[1] = stheta * sphi;
505 +  epsi[2] = ctheta;
506 +  
507    for (int j = 0 ; j<3; j++)
508      grad[j] = frc[j];
509  
510 +  grad[3] = 0;
511 +  grad[4] = 0;
512 +  grad[5] = 0;
513 +
514    for (int j = 0; j < 3; j++ ) {
515      
516      grad[3] += trq[j]*ephi[j];
# Line 449 | Line 521 | void DirectionalAtom::getGrad( double grad[6] ) {
521  
522   }
523  
524 <
524 > /**
525 >  * getEulerAngles computes a set of Euler angle values consistent
526 >  *  with an input rotation matrix.  They are returned in the following
527 >  * order:
528 >  *  myEuler[0] = phi;
529 >  *  myEuler[1] = theta;
530 >  *  myEuler[2] = psi;
531 > */
532   void DirectionalAtom::getEulerAngles(double myEuler[3]) {
533  
534 <  // getEulerAngles computes a set of Euler angle values consistent
535 <  // with an input rotation matrix.  They are returned in the following
536 <  // order:
537 <  //  myEuler[0] = phi;
538 <  //  myEuler[1] = theta;
460 <  //  myEuler[2] = psi;
534 >  // We use so-called "x-convention", which is the most common definition.
535 >  // In this convention, the rotation given by Euler angles (phi, theta, psi), where the first
536 >  // rotation is by an angle phi about the z-axis, the second is by an angle  
537 >  // theta (0 <= theta <= 180)about the x-axis, and thethird is by an angle psi about the
538 >  //z-axis (again).
539    
540 +  
541    double phi,theta,psi,eps;
542    double pi;
543    double cphi,ctheta,cpsi;
# Line 469 | Line 548 | void DirectionalAtom::getEulerAngles(double myEuler[3]
548    // set the tolerance for Euler angles and rotation elements
549    
550    eps = 1.0e-8;
472 –    
473 –  // get a trial value of theta from a single rotation element
474 –  
475 –  theta = asin(min(1.0,max(-1.0,-Amat[Axz])));
476 –  ctheta = cos(theta);
477 –  stheta = -Amat[Axz];
478 –  
479 –  // set the phi/psi difference when theta is either 90 or -90
480 –  
481 –  if (fabs(ctheta) <= eps) {
482 –    phi = 0.0;
483 –    if (fabs(Amat[Azx]) < eps) {
484 –      psi = asin(min(1.0,max(-1.0,-Amat[Ayx]/Amat[Axz])));
485 –    } else {
486 –      if (fabs(Amat[Ayx]) < eps) {
487 –        psi = acos(min(1.0,max(-1.0,-Amat[Azx]/Amat[Axz])));
488 –      } else {
489 –        psi = atan(Amat[Ayx]/Amat[Azx]);
490 –      }    
491 –    }
492 –  }
551  
552 <  // set the phi and psi values for all other theta values
552 >  theta = acos(min(1.0,max(-1.0,Amat[Azz])));
553 >  ctheta = Amat[Azz];
554 >  stheta = sqrt(1.0 - ctheta * ctheta);
555 >
556 >  // when sin(theta) is close to 0, we need to consider singularity
557 >  // In this case, we can assign an arbitary value to phi (or psi), and then determine
558 >  // the psi (or phi) or vice-versa. We'll assume that phi always gets the rotation, and psi is 0
559 >  // in cases of singularity.  
560 >  // we use atan2 instead of atan, since atan2 will give us -Pi to Pi.
561 >  // Since 0 <= theta <= 180, sin(theta) will be always non-negative. Therefore, it never
562 >  // change the sign of both of the parameters passed to atan2.
563    
564 <  else {
565 <    if (fabs(Amat[Axx]) < eps) {
566 <      phi = asin(min(1.0,max(-1.0,Amat[Axy]/ctheta)));
499 <    } else {
500 <      if (fabs(Amat[Axy]) < eps) {
501 <        phi = acos(min(1.0,max(-1.0,Amat[Axx]/ctheta)));
502 <      } else {
503 <        phi = atan(Amat[Axy]/Amat[Axx]);
504 <      }
505 <    }
506 <    if (fabs(Amat[Azz]) < eps) {
507 <      psi = asin(min(1.0,max(-1.0,Amat[Ayz]/ctheta)));
508 <    } else {
509 <      if (fabs(Amat[Ayz]) < eps) {
510 <        psi = acos(min(1.0,max(-1.0,Amat[Azz]/ctheta)));
511 <      }
512 <      psi = atan(Amat[Ayz]/Amat[Azz]);
513 <    }
564 >  if (fabs(stheta) <= eps){
565 >    psi = 0.0;
566 >    phi = atan2(-Amat[Ayx], Amat[Axx]);  
567    }
568 <
569 <  // find sine and cosine of the trial phi and psi values
570 <
571 <  cphi = cos(phi);
519 <  sphi = sin(phi);
520 <  cpsi = cos(psi);
521 <  spsi = sin(psi);
522 <
523 <  // reconstruct the diagonal of the rotation matrix
524 <
525 <  b[0] = ctheta * cphi;
526 <  b[1] = spsi*stheta*sphi + cpsi*cphi;
527 <  b[2] = ctheta * cpsi;
528 <
529 <  // compare the correct matrix diagonal to rebuilt diagonal
530 <
531 <  for (int i = 0; i < 3; i++) {
532 <    flip[i] = 0;
533 <    if (fabs(Amat[3*i + i] - b[i]) > eps)  flip[i] = 1;
568 >  // we only have one unique solution
569 >  else{    
570 >      phi = atan2(Amat[Azx], -Amat[Azy]);
571 >      psi = atan2(Amat[Axz], Amat[Ayz]);
572    }
573  
574 <  // alter Euler angles to get correct rotation matrix values
575 <  
576 <  if (flip[0] && flip[1]) phi = phi - copysign(M_PI,phi);
539 <  if (flip[0] && flip[2]) theta = -theta + copysign(M_PI, theta);
540 <  if (flip[1] && flip[2]) psi = psi - copysign(M_PI, psi);
574 >  //wrap phi and psi, make sure they are in the range from 0 to 2*Pi
575 >  //if (phi < 0)
576 >  //  phi += M_PI;
577  
578 +  //if (psi < 0)
579 +  //  psi += M_PI;
580 +
581    myEuler[0] = phi;
582    myEuler[1] = theta;
583    myEuler[2] = psi;
584 <
584 >  
585    return;
586   }
587  

Diff Legend

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