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 689 by tim, Tue Aug 12 19:56:49 2003 UTC vs.
Revision 1452 by tim, Mon Aug 23 15:11:36 2004 UTC

# Line 1 | Line 1
1 < #include <cmath>
1 > #include <math.h>
2  
3   #include "Atom.hpp"
4 + #include "DirectionalAtom.hpp"
5   #include "simError.h"
6 + #include "MatVec3.h"
7  
8   void DirectionalAtom::zeroForces() {
9    if( hasCoords ){
10 <    frc[offsetX] = 0.0;
11 <    frc[offsetY] = 0.0;
10 <    frc[offsetZ] = 0.0;
10 >
11 >    Atom::zeroForces();
12      
13      trq[offsetX] = 0.0;
14      trq[offsetY] = 0.0;
# Line 34 | Line 35 | void DirectionalAtom::setCoords(void){
35                       &trq,
36                       &Amat,
37                       &mu,  
38 <                     &ul );
38 >                     &ul,
39 >                     &quat);
40    }
41    else{
42      sprintf( painCave.errMsg,
43               "Attempted to set Atom %d  coordinates with an unallocated "
44 <             "SimState object.\n" );
44 >             "SimState object.\n", index );
45      painCave.isFatal = 1;
46      simError();
47    }
48  
49    hasCoords = true;
50  
49 –  mu[index] = myMu;
50 –
51   }
52  
53 – double DirectionalAtom::getMu( void ) {
54 –
55 –  if( hasCoords ){
56 –    return mu[index];
57 –  }
58 –  else{
59 –    return myMu;
60 –  }
61 –  return 0;
62 – }
63 –
64 – void DirectionalAtom::setMu( double the_mu ) {
65 –
66 –  if( hasCoords ){
67 –    mu[index] = the_mu;
68 –    myMu = the_mu;
69 –  }
70 –  else{
71 –    myMu = the_mu;
72 –  }
73 – }
74 –
53   void DirectionalAtom::setA( double the_A[3][3] ){
54  
55    if( hasCoords ){
# Line 91 | Line 69 | void DirectionalAtom::setA( double the_A[3][3] ){
69    }
70   }
71  
72 < void DirectionalAtom::setI( double the_I[3][3] ){
72 > void DirectionalAtom::setI( double the_I[3][3] ){  
73    
74    Ixx = the_I[0][0]; Ixy = the_I[0][1]; Ixz = the_I[0][2];
75    Iyx = the_I[1][0]; Iyy = the_I[1][1]; Iyz = the_I[1][2];
# Line 181 | Line 159 | void DirectionalAtom::getU( double the_u[3] ){
159  
160   void DirectionalAtom::getU( double the_u[3] ){
161    
162 <  the_u[0] = sux;
163 <  the_u[1] = suy;
164 <  the_u[2] = suz;
165 <
162 >  the_u[0] = sU[2][0];
163 >  the_u[1] = sU[2][1];
164 >  the_u[2] = sU[2][2];
165 >  
166    this->body2Lab( the_u );
167   }
168  
# Line 246 | Line 224 | void DirectionalAtom::getQ( double q[4] ){
224    }
225   }
226  
227 + void DirectionalAtom::setUnitFrameFromEuler(double phi,
228 +                                            double theta,
229 +                                            double psi) {
230  
231 +  double myA[3][3];
232 +  double uFrame[3][3];
233 +  double len;
234 +  int i, j;
235 +  
236 +  myA[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
237 +  myA[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
238 +  myA[0][2] = sin(theta) * sin(psi);
239 +  
240 +  myA[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
241 +  myA[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
242 +  myA[1][2] = sin(theta) * cos(psi);
243 +  
244 +  myA[2][0] = sin(phi) * sin(theta);
245 +  myA[2][1] = -cos(phi) * sin(theta);
246 +  myA[2][2] = cos(theta);
247 +  
248 +  // Make the unit Frame:
249 +
250 +  for (i=0; i < 3; i++)
251 +    for (j=0; j < 3; j++)
252 +      uFrame[i][j] = 0.0;
253 +
254 +  for (i=0; i < 3; i++)
255 +    uFrame[i][i] = 1.0;
256 +
257 +  // rotate by the given rotation matrix:
258 +
259 +  matMul3(myA, uFrame, sU);
260 +
261 +  // renormalize column vectors:
262 +
263 +  for (i=0; i < 3; i++) {
264 +    len = 0.0;
265 +    for (j = 0; j < 3; j++) {
266 +      len += sU[i][j]*sU[i][j];
267 +    }
268 +    len = sqrt(len);
269 +    for (j = 0; j < 3; j++) {
270 +      sU[i][j] /= len;    
271 +    }
272 +  }
273 +  
274 +  // sU now contains the coordinates of the 'special' frame;
275 +    
276 + }
277 +
278   void DirectionalAtom::setEuler( double phi, double theta, double psi ){
279    
280    if( hasCoords ){
# Line 299 | Line 327 | void DirectionalAtom::lab2Body( double r[3] ){
327  
328   }
329  
330 + void DirectionalAtom::rotateBy( double by_A[3][3]) {
331 +
332 +  // Check this
333 +  
334 +  double r00, r01, r02, r10, r11, r12, r20, r21, r22;
335 +
336 +  if( hasCoords ){
337 +
338 +    r00 = by_A[0][0]*Amat[Axx] + by_A[0][1]*Amat[Ayx] + by_A[0][2]*Amat[Azx];
339 +    r01 = by_A[0][0]*Amat[Axy] + by_A[0][1]*Amat[Ayy] + by_A[0][2]*Amat[Azy];
340 +    r02 = by_A[0][0]*Amat[Axz] + by_A[0][1]*Amat[Ayz] + by_A[0][2]*Amat[Azz];
341 +    
342 +    r10 = by_A[1][0]*Amat[Axx] + by_A[1][1]*Amat[Ayx] + by_A[1][2]*Amat[Azx];
343 +    r11 = by_A[1][0]*Amat[Axy] + by_A[1][1]*Amat[Ayy] + by_A[1][2]*Amat[Azy];
344 +    r12 = by_A[1][0]*Amat[Axz] + by_A[1][1]*Amat[Ayz] + by_A[1][2]*Amat[Azz];
345 +    
346 +    r20 = by_A[2][0]*Amat[Axx] + by_A[2][1]*Amat[Ayx] + by_A[2][2]*Amat[Azx];
347 +    r21 = by_A[2][0]*Amat[Axy] + by_A[2][1]*Amat[Ayy] + by_A[2][2]*Amat[Azy];
348 +    r22 = by_A[2][0]*Amat[Axz] + by_A[2][1]*Amat[Ayz] + by_A[2][2]*Amat[Azz];
349 +    
350 +    Amat[Axx] = r00; Amat[Axy] = r01; Amat[Axz] = r02;
351 +    Amat[Ayx] = r10; Amat[Ayy] = r11; Amat[Ayz] = r12;
352 +    Amat[Azx] = r20; Amat[Azy] = r21; Amat[Azz] = r22;
353 +
354 +  }
355 +  else{
356 +    
357 +    sprintf( painCave.errMsg,
358 +             "Attempt to rotate frame for atom %d before coords set.\n",
359 +             index );
360 +    painCave.isFatal = 1;
361 +    simError();
362 +  }
363 +
364 + }
365 +
366 +
367   void DirectionalAtom::body2Lab( double r[3] ){
368  
369    double rb[3]; // the body frame vector
# Line 325 | Line 390 | void DirectionalAtom::updateU( void ){
390   void DirectionalAtom::updateU( void ){
391  
392    if( hasCoords ){
393 <    ul[offsetX] = (Amat[Axx] * sux) + (Amat[Ayx] * suy) + (Amat[Azx] * suz);
394 <    ul[offsetY] = (Amat[Axy] * sux) + (Amat[Ayy] * suy) + (Amat[Azy] * suz);
395 <    ul[offsetZ] = (Amat[Axz] * sux) + (Amat[Ayz] * suy) + (Amat[Azz] * suz);
393 >    ul[offsetX] = (Amat[Axx] * sU[2][0]) +
394 >      (Amat[Ayx] * sU[2][1]) + (Amat[Azx] * sU[2][2]);
395 >    ul[offsetY] = (Amat[Axy] * sU[2][0]) +
396 >      (Amat[Ayy] * sU[2][1]) + (Amat[Azy] * sU[2][2]);
397 >    ul[offsetZ] = (Amat[Axz] * sU[2][0]) +
398 >      (Amat[Ayz] * sU[2][1]) + (Amat[Azz] * sU[2][2]);
399    }
400    else{
401      
# Line 370 | Line 438 | void DirectionalAtom::getTrq( double theT[3] ){
438    }
439   }
440  
441 + void DirectionalAtom::setTrq( double theT[3] ){
442 +  
443 +  if( hasCoords ){
444 +    trq[offsetX] = theT[0];
445 +    trq[offsetY] = theT[1];
446 +    trq[offsetZ] = theT[2];
447 +  }
448 +  else{
449 +    
450 +    sprintf( painCave.errMsg,
451 +             "Attempt to add Trq for atom %d before coords set.\n",
452 +             index );
453 +    painCave.isFatal = 1;
454 +    simError();
455 +  }
456 + }
457 +
458   void DirectionalAtom::addTrq( double theT[3] ){
459    
460    if( hasCoords ){
# Line 401 | Line 486 | void DirectionalAtom::getI( double the_I[3][3] ){
486    the_I[2][0] = Izx;
487    the_I[2][1] = Izy;
488    the_I[2][2] = Izz;
489 + }
490 +
491 + void DirectionalAtom::getGrad( double grad[6] ) {
492 +
493 +  double myEuler[3];
494 +  double phi, theta, psi;
495 +  double cphi, sphi, ctheta, stheta;
496 +  double ephi[3];
497 +  double etheta[3];
498 +  double epsi[3];
499 +
500 +  this->getEulerAngles(myEuler);
501 +
502 +  phi = myEuler[0];
503 +  theta = myEuler[1];
504 +  psi = myEuler[2];
505 +
506 +  cphi = cos(phi);
507 +  sphi = sin(phi);
508 +  ctheta = cos(theta);
509 +  stheta = sin(theta);
510 +
511 +  // get unit vectors along the phi, theta and psi rotation axes
512 +
513 +  ephi[0] = 0.0;
514 +  ephi[1] = 0.0;
515 +  ephi[2] = 1.0;
516 +
517 +  etheta[0] = cphi;
518 +  etheta[1] = sphi;
519 +  etheta[2] = 0.0;
520 +  
521 +  epsi[0] = stheta * cphi;
522 +  epsi[1] = stheta * sphi;
523 +  epsi[2] = ctheta;
524 +  
525 +  for (int j = 0 ; j<3; j++)
526 +    grad[j] = frc[j];
527 +
528 +  grad[3] = 0;
529 +  grad[4] = 0;
530 +  grad[5] = 0;
531 +
532 +  for (int j = 0; j < 3; j++ ) {
533 +    
534 +    grad[3] += trq[j]*ephi[j];
535 +    grad[4] += trq[j]*etheta[j];
536 +    grad[5] += trq[j]*epsi[j];
537 +    
538 +  }
539 +
540   }
541 +
542 + /**
543 +  * getEulerAngles computes a set of Euler angle values consistent
544 +  *  with an input rotation matrix.  They are returned in the following
545 +  * order:
546 +  *  myEuler[0] = phi;
547 +  *  myEuler[1] = theta;
548 +  *  myEuler[2] = psi;
549 + */
550 + void DirectionalAtom::getEulerAngles(double myEuler[3]) {
551 +
552 +  // We use so-called "x-convention", which is the most common definition.
553 +  // In this convention, the rotation given by Euler angles (phi, theta, psi), where the first
554 +  // rotation is by an angle phi about the z-axis, the second is by an angle  
555 +  // theta (0 <= theta <= 180)about the x-axis, and thethird is by an angle psi about the
556 +  //z-axis (again).
557 +  
558 +  
559 +  double phi,theta,psi,eps;
560 +  double ctheta,stheta;
561 +
562 +  // set the tolerance for Euler angles and rotation elements
563 +  
564 +  eps = 1.0e-8;
565 +
566 +  theta = acos(min(1.0,max(-1.0,Amat[Azz])));
567 +  ctheta = Amat[Azz];
568 +  stheta = sqrt(1.0 - ctheta * ctheta);
569 +
570 +  // when sin(theta) is close to 0, we need to consider singularity
571 +  // In this case, we can assign an arbitary value to phi (or psi), and then determine
572 +  // the psi (or phi) or vice-versa. We'll assume that phi always gets the rotation, and psi is 0
573 +  // in cases of singularity.  
574 +  // we use atan2 instead of atan, since atan2 will give us -Pi to Pi.
575 +  // Since 0 <= theta <= 180, sin(theta) will be always non-negative. Therefore, it never
576 +  // change the sign of both of the parameters passed to atan2.
577 +  
578 +  if (fabs(stheta) <= eps){
579 +    psi = 0.0;
580 +    phi = atan2(-Amat[Ayx], Amat[Axx]);  
581 +  }
582 +  // we only have one unique solution
583 +  else{    
584 +      phi = atan2(Amat[Azx], -Amat[Azy]);
585 +      psi = atan2(Amat[Axz], Amat[Ayz]);
586 +  }
587 +
588 +  //wrap phi and psi, make sure they are in the range from 0 to 2*Pi
589 +  //if (phi < 0)
590 +  //  phi += M_PI;
591 +
592 +  //if (psi < 0)
593 +  //  psi += M_PI;
594 +
595 +  myEuler[0] = phi;
596 +  myEuler[1] = theta;
597 +  myEuler[2] = psi;
598 +  
599 +  return;
600 + }
601 +
602 + double DirectionalAtom::getZangle( ){
603 +  
604 +  if( hasCoords ){
605 +    return zAngle;
606 +  }
607 +  else{
608 +    
609 +    sprintf( painCave.errMsg,
610 +             "Attempt to get zAngle for atom %d before coords set.\n",
611 +             index );
612 +    painCave.isFatal = 1;
613 +    simError();
614 +    return 0;
615 +  }
616 + }
617 +
618 + void DirectionalAtom::setZangle( double zAng ){
619 +  
620 +  if( hasCoords ){
621 +    zAngle = zAng;
622 +  }
623 +  else{
624 +    
625 +    sprintf( painCave.errMsg,
626 +             "Attempt to set zAngle for atom %d before coords set.\n",
627 +             index );
628 +    painCave.isFatal = 1;
629 +    simError();
630 +  }
631 + }
632 +
633 + void DirectionalAtom::addZangle( double zAng ){
634 +  
635 +  if( hasCoords ){
636 +    zAngle += zAng;
637 +  }
638 +  else{
639 +    
640 +    sprintf( painCave.errMsg,
641 +             "Attempt to add zAngle to atom %d before coords set.\n",
642 +             index );
643 +    painCave.isFatal = 1;
644 +    simError();
645 +  }
646 + }
647 +
648 + double DirectionalAtom::max(double x, double  y) {  
649 +  return (x > y) ? x : y;
650 + }
651 +
652 + double DirectionalAtom::min(double x, double  y) {  
653 +  return (x > y) ? y : x;
654 + }

Diff Legend

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