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 670 by mmeineke, Thu Aug 7 21:47:18 2003 UTC vs.
Revision 1136 by tim, Tue Apr 27 16:26:44 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 23 | Line 24 | double DirectionalAtom::getMu( void ) {
24    }
25   }
26  
27 < double DirectionalAtom::getMu( void ) {
27 > void DirectionalAtom::setCoords(void){
28  
29 <  if( hasCoords ){
30 <    return mu[index];
29 >  if( myConfig->isAllocated() ){
30 >
31 >    myConfig->getAtomPointers( index,
32 >                     &pos,
33 >                     &vel,
34 >                     &frc,
35 >                     &trq,
36 >                     &Amat,
37 >                     &mu,  
38 >                     &ul,
39 >                 &rc,
40 >                     &massRatio);
41    }
42    else{
32 –    
43      sprintf( painCave.errMsg,
44 <             "Attempt to get Mu for atom %d before coords set.\n",
45 <             index );
44 >             "Attempted to set Atom %d  coordinates with an unallocated "
45 >             "SimState object.\n", index );
46      painCave.isFatal = 1;
47      simError();
48    }
39 –  return 0;
40 – }
49  
50 < void DirectionalAtom::setMu( double the_mu ) {
50 >  hasCoords = true;
51  
44 –  if( hasCoords ){
45 –    mu[index] = the_mu;
46 –  }
47 –  else{
48 –    
49 –    sprintf( painCave.errMsg,
50 –             "Attempt to set Mu for atom %d before coords set.\n",
51 –             index );
52 –    painCave.isFatal = 1;
53 –    simError();
54 –  }
52   }
53  
54   void DirectionalAtom::setA( double the_A[3][3] ){
# Line 73 | Line 70 | void DirectionalAtom::setI( double the_I[3][3] ){
70    }
71   }
72  
73 < void DirectionalAtom::setI( double the_I[3][3] ){
73 > void DirectionalAtom::setI( double the_I[3][3] ){  
74    
75    Ixx = the_I[0][0]; Ixy = the_I[0][1]; Ixz = the_I[0][2];
76    Iyx = the_I[1][0]; Iyy = the_I[1][1]; Iyz = the_I[1][2];
# Line 163 | Line 160 | void DirectionalAtom::getU( double the_u[3] ){
160  
161   void DirectionalAtom::getU( double the_u[3] ){
162    
163 <  the_u[0] = sux;
164 <  the_u[1] = suy;
165 <  the_u[2] = suz;
166 <
163 >  the_u[0] = sU[2][0];
164 >  the_u[1] = sU[2][1];
165 >  the_u[2] = sU[2][2];
166 >  
167    this->body2Lab( the_u );
168   }
169  
# Line 228 | Line 225 | void DirectionalAtom::getQ( double q[4] ){
225    }
226   }
227  
228 + void DirectionalAtom::setUnitFrameFromEuler(double phi,
229 +                                            double theta,
230 +                                            double psi) {
231  
232 +  double myA[3][3];
233 +  double uFrame[3][3];
234 +  double len;
235 +  int i, j;
236 +  
237 +  myA[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
238 +  myA[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
239 +  myA[0][2] = sin(theta) * sin(psi);
240 +  
241 +  myA[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
242 +  myA[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
243 +  myA[1][2] = sin(theta) * cos(psi);
244 +  
245 +  myA[2][0] = sin(phi) * sin(theta);
246 +  myA[2][1] = -cos(phi) * sin(theta);
247 +  myA[2][2] = cos(theta);
248 +  
249 +  // Make the unit Frame:
250 +
251 +  for (i=0; i < 3; i++)
252 +    for (j=0; j < 3; j++)
253 +      uFrame[i][j] = 0.0;
254 +
255 +  for (i=0; i < 3; i++)
256 +    uFrame[i][i] = 1.0;
257 +
258 +  // rotate by the given rotation matrix:
259 +
260 +  matMul3(myA, uFrame, sU);
261 +
262 +  // renormalize column vectors:
263 +
264 +  for (i=0; i < 3; i++) {
265 +    len = 0.0;
266 +    for (j = 0; j < 3; j++) {
267 +      len += sU[i][j]*sU[i][j];
268 +    }
269 +    len = sqrt(len);
270 +    for (j = 0; j < 3; j++) {
271 +      sU[i][j] /= len;    
272 +    }
273 +  }
274 +  
275 +  // sU now contains the coordinates of the 'special' frame;
276 +    
277 + }
278 +
279   void DirectionalAtom::setEuler( double phi, double theta, double psi ){
280    
281    if( hasCoords ){
# Line 274 | Line 321 | void DirectionalAtom::lab2Body( double r[3] ){
321      
322      sprintf( painCave.errMsg,
323               "Attempt to convert lab2body for atom %d before coords set.\n",
324 +             index );
325 +    painCave.isFatal = 1;
326 +    simError();
327 +  }
328 +
329 + }
330 +
331 + void DirectionalAtom::rotateBy( double by_A[3][3]) {
332 +
333 +  // Check this
334 +  
335 +  double r00, r01, r02, r10, r11, r12, r20, r21, r22;
336 +
337 +  if( hasCoords ){
338 +
339 +    r00 = by_A[0][0]*Amat[Axx] + by_A[0][1]*Amat[Ayx] + by_A[0][2]*Amat[Azx];
340 +    r01 = by_A[0][0]*Amat[Axy] + by_A[0][1]*Amat[Ayy] + by_A[0][2]*Amat[Azy];
341 +    r02 = by_A[0][0]*Amat[Axz] + by_A[0][1]*Amat[Ayz] + by_A[0][2]*Amat[Azz];
342 +    
343 +    r10 = by_A[1][0]*Amat[Axx] + by_A[1][1]*Amat[Ayx] + by_A[1][2]*Amat[Azx];
344 +    r11 = by_A[1][0]*Amat[Axy] + by_A[1][1]*Amat[Ayy] + by_A[1][2]*Amat[Azy];
345 +    r12 = by_A[1][0]*Amat[Axz] + by_A[1][1]*Amat[Ayz] + by_A[1][2]*Amat[Azz];
346 +    
347 +    r20 = by_A[2][0]*Amat[Axx] + by_A[2][1]*Amat[Ayx] + by_A[2][2]*Amat[Azx];
348 +    r21 = by_A[2][0]*Amat[Axy] + by_A[2][1]*Amat[Ayy] + by_A[2][2]*Amat[Azy];
349 +    r22 = by_A[2][0]*Amat[Axz] + by_A[2][1]*Amat[Ayz] + by_A[2][2]*Amat[Azz];
350 +    
351 +    Amat[Axx] = r00; Amat[Axy] = r01; Amat[Axz] = r02;
352 +    Amat[Ayx] = r10; Amat[Ayy] = r11; Amat[Ayz] = r12;
353 +    Amat[Azx] = r20; Amat[Azy] = r21; Amat[Azz] = r22;
354 +
355 +  }
356 +  else{
357 +    
358 +    sprintf( painCave.errMsg,
359 +             "Attempt to rotate frame for atom %d before coords set.\n",
360               index );
361      painCave.isFatal = 1;
362      simError();
# Line 281 | Line 364 | void DirectionalAtom::body2Lab( double r[3] ){
364  
365   }
366  
367 +
368   void DirectionalAtom::body2Lab( double r[3] ){
369  
370    double rb[3]; // the body frame vector
# Line 307 | Line 391 | void DirectionalAtom::updateU( void ){
391   void DirectionalAtom::updateU( void ){
392  
393    if( hasCoords ){
394 <    ul[offsetX] = (Amat[Axx] * sux) + (Amat[Ayx] * suy) + (Amat[Azx] * suz);
395 <    ul[offsetY] = (Amat[Axy] * sux) + (Amat[Ayy] * suy) + (Amat[Azy] * suz);
396 <    ul[offsetZ] = (Amat[Axz] * sux) + (Amat[Ayz] * suy) + (Amat[Azz] * suz);
394 >    ul[offsetX] = (Amat[Axx] * sU[2][0]) +
395 >      (Amat[Ayx] * sU[2][1]) + (Amat[Azx] * sU[2][2]);
396 >    ul[offsetY] = (Amat[Axy] * sU[2][0]) +
397 >      (Amat[Ayy] * sU[2][1]) + (Amat[Azy] * sU[2][2]);
398 >    ul[offsetZ] = (Amat[Axz] * sU[2][0]) +
399 >      (Amat[Ayz] * sU[2][1]) + (Amat[Azz] * sU[2][2]);
400    }
401    else{
402      
# Line 384 | Line 471 | void DirectionalAtom::getI( double the_I[3][3] ){
471    the_I[2][1] = Izy;
472    the_I[2][2] = Izz;
473   }
474 +
475 + void DirectionalAtom::getGrad( double grad[6] ) {
476 +
477 +  double myEuler[3];
478 +  double phi, theta, psi;
479 +  double cphi, sphi, ctheta, stheta;
480 +  double ephi[3];
481 +  double etheta[3];
482 +  double epsi[3];
483 +
484 +  this->getEulerAngles(myEuler);
485 +
486 +  phi = myEuler[0];
487 +  theta = myEuler[1];
488 +  psi = myEuler[2];
489 +
490 +  cphi = cos(phi);
491 +  sphi = sin(phi);
492 +  ctheta = cos(theta);
493 +  stheta = sin(theta);
494 +
495 +  // get unit vectors along the phi, theta and psi rotation axes
496 +
497 +  ephi[0] = 0.0;
498 +  ephi[1] = 0.0;
499 +  ephi[2] = 1.0;
500 +
501 +  etheta[0] = cphi;
502 +  etheta[1] = sphi;
503 +  etheta[2] = 0.0;
504 +  
505 +  epsi[0] = stheta * cphi;
506 +  epsi[1] = stheta * sphi;
507 +  epsi[2] = ctheta;
508 +  
509 +  for (int j = 0 ; j<3; j++)
510 +    grad[j] = frc[j];
511 +
512 +  grad[3] = 0;
513 +  grad[4] = 0;
514 +  grad[5] = 0;
515 +
516 +  for (int j = 0; j < 3; j++ ) {
517 +    
518 +    grad[3] += trq[j]*ephi[j];
519 +    grad[4] += trq[j]*etheta[j];
520 +    grad[5] += trq[j]*epsi[j];
521 +    
522 +  }
523 +
524 + }
525 +
526 + /**
527 +  * getEulerAngles computes a set of Euler angle values consistent
528 +  *  with an input rotation matrix.  They are returned in the following
529 +  * order:
530 +  *  myEuler[0] = phi;
531 +  *  myEuler[1] = theta;
532 +  *  myEuler[2] = psi;
533 + */
534 + void DirectionalAtom::getEulerAngles(double myEuler[3]) {
535 +
536 +  // We use so-called "x-convention", which is the most common definition.
537 +  // In this convention, the rotation given by Euler angles (phi, theta, psi), where the first
538 +  // rotation is by an angle phi about the z-axis, the second is by an angle  
539 +  // theta (0 <= theta <= 180)about the x-axis, and thethird is by an angle psi about the
540 +  //z-axis (again).
541 +  
542 +  
543 +  double phi,theta,psi,eps;
544 +  double pi;
545 +  double cphi,ctheta,cpsi;
546 +  double sphi,stheta,spsi;
547 +  double b[3];
548 +  int flip[3];
549 +
550 +  // set the tolerance for Euler angles and rotation elements
551 +  
552 +  eps = 1.0e-8;
553 +
554 +  theta = acos(min(1.0,max(-1.0,Amat[Azz])));
555 +  ctheta = Amat[Azz];
556 +  stheta = sqrt(1.0 - ctheta * ctheta);
557 +
558 +  // when sin(theta) is close to 0, we need to consider singularity
559 +  // In this case, we can assign an arbitary value to phi (or psi), and then determine
560 +  // the psi (or phi) or vice-versa. We'll assume that phi always gets the rotation, and psi is 0
561 +  // in cases of singularity.  
562 +  // we use atan2 instead of atan, since atan2 will give us -Pi to Pi.
563 +  // Since 0 <= theta <= 180, sin(theta) will be always non-negative. Therefore, it never
564 +  // change the sign of both of the parameters passed to atan2.
565 +  
566 +  if (fabs(stheta) <= eps){
567 +    psi = 0.0;
568 +    phi = atan2(-Amat[Ayx], Amat[Axx]);  
569 +  }
570 +  // we only have one unique solution
571 +  else{    
572 +      phi = atan2(Amat[Azx], -Amat[Azy]);
573 +      psi = atan2(Amat[Axz], Amat[Ayz]);
574 +  }
575 +
576 +  //wrap phi and psi, make sure they are in the range from 0 to 2*Pi
577 +  //if (phi < 0)
578 +  //  phi += M_PI;
579 +
580 +  //if (psi < 0)
581 +  //  psi += M_PI;
582 +
583 +  myEuler[0] = phi;
584 +  myEuler[1] = theta;
585 +  myEuler[2] = psi;
586 +  
587 +  return;
588 + }
589 +
590 + double DirectionalAtom::max(double x, double  y) {  
591 +  return (x > y) ? x : y;
592 + }
593 +
594 + double DirectionalAtom::min(double x, double  y) {  
595 +  return (x > y) ? y : x;
596 + }

Diff Legend

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