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 597 by mmeineke, Mon Jul 14 21:28:54 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  
11 +    Atom::zeroForces();
12 +    
13 +    trq[offsetX] = 0.0;
14 +    trq[offsetY] = 0.0;
15 +    trq[offsetZ] = 0.0;
16 +  }
17 +  else{
18 +    
19 +    sprintf( painCave.errMsg,
20 +             "Attempt to zero frc and trq for atom %d before coords set.\n",
21 +             index );
22 +    painCave.isFatal = 1;
23 +    simError();
24 +  }
25 + }
26  
27 + void DirectionalAtom::setCoords(void){
28 +
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{
43 +    sprintf( painCave.errMsg,
44 +             "Attempted to set Atom %d  coordinates with an unallocated "
45 +             "SimState object.\n", index );
46 +    painCave.isFatal = 1;
47 +    simError();
48 +  }
49 +
50 +  hasCoords = true;
51 +
52 + }
53 +
54   void DirectionalAtom::setA( double the_A[3][3] ){
8 –  
9 –  Amat[Axx] = the_A[0][0]; Amat[Axy] = the_A[0][1]; Amat[Axz] = the_A[0][2];
10 –  Amat[Ayx] = the_A[1][0]; Amat[Ayy] = the_A[1][1]; Amat[Ayz] = the_A[1][2];
11 –  Amat[Azx] = the_A[2][0]; Amat[Azy] = the_A[2][1]; Amat[Azz] = the_A[2][2];
55  
56 <  this->updateU();  
56 >  if( hasCoords ){
57 >    Amat[Axx] = the_A[0][0]; Amat[Axy] = the_A[0][1]; Amat[Axz] = the_A[0][2];
58 >    Amat[Ayx] = the_A[1][0]; Amat[Ayy] = the_A[1][1]; Amat[Ayz] = the_A[1][2];
59 >    Amat[Azx] = the_A[2][0]; Amat[Azy] = the_A[2][1]; Amat[Azz] = the_A[2][2];
60 >    
61 >    this->updateU();  
62 >  }
63 >  else{
64 >    
65 >    sprintf( painCave.errMsg,
66 >             "Attempt to set Amat for atom %d before coords set.\n",
67 >             index );
68 >    painCave.isFatal = 1;
69 >    simError();
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 24 | Line 81 | void DirectionalAtom::setQ( double the_q[4] ){
81  
82    double q0Sqr, q1Sqr, q2Sqr, q3Sqr;
83  
84 <  q0Sqr = the_q[0] * the_q[0];
85 <  q1Sqr = the_q[1] * the_q[1];
86 <  q2Sqr = the_q[2] * the_q[2];
87 <  q3Sqr = the_q[3] * the_q[3];
88 <  
84 >  if( hasCoords ){
85 >    q0Sqr = the_q[0] * the_q[0];
86 >    q1Sqr = the_q[1] * the_q[1];
87 >    q2Sqr = the_q[2] * the_q[2];
88 >    q3Sqr = the_q[3] * the_q[3];
89 >    
90 >    
91 >    Amat[Axx] = q0Sqr + q1Sqr - q2Sqr - q3Sqr;
92 >    Amat[Axy] = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] );
93 >    Amat[Axz] = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] );
94 >    
95 >    Amat[Ayx] = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] );
96 >    Amat[Ayy] = q0Sqr - q1Sqr + q2Sqr - q3Sqr;
97 >    Amat[Ayz] = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] );
98 >    
99 >    Amat[Azx] = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] );
100 >    Amat[Azy] = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] );
101 >    Amat[Azz] = q0Sqr - q1Sqr -q2Sqr +q3Sqr;
102 >    
103 >    this->updateU();
104 >  }
105 >  else{
106 >    
107 >    sprintf( painCave.errMsg,
108 >             "Attempt to set Q for atom %d before coords set.\n",
109 >             index );
110 >    painCave.isFatal = 1;
111 >    simError();
112 >  }
113  
33 –  Amat[Axx] = q0Sqr + q1Sqr - q2Sqr - q3Sqr;
34 –  Amat[Axy] = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] );
35 –  Amat[Axz] = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] );
36 –  
37 –  Amat[Ayx] = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] );
38 –  Amat[Ayy] = q0Sqr - q1Sqr + q2Sqr - q3Sqr;
39 –  Amat[Ayz] = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] );
40 –
41 –  Amat[Azx] = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] );
42 –  Amat[Azy] = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] );
43 –  Amat[Azz] = q0Sqr - q1Sqr -q2Sqr +q3Sqr;
44 –
45 –  this->updateU();
114   }
115  
116   void DirectionalAtom::getA( double the_A[3][3] ){
117    
118 <  the_A[0][0] = Amat[Axx];
119 <  the_A[0][1] = Amat[Axy];
120 <  the_A[0][2] = Amat[Axz];
118 >  if( hasCoords ){
119 >    the_A[0][0] = Amat[Axx];
120 >    the_A[0][1] = Amat[Axy];
121 >    the_A[0][2] = Amat[Axz];
122 >    
123 >    the_A[1][0] = Amat[Ayx];
124 >    the_A[1][1] = Amat[Ayy];
125 >    the_A[1][2] = Amat[Ayz];
126 >    
127 >    the_A[2][0] = Amat[Azx];
128 >    the_A[2][1] = Amat[Azy];
129 >    the_A[2][2] = Amat[Azz];
130 >  }
131 >  else{
132 >    
133 >    sprintf( painCave.errMsg,
134 >             "Attempt to get Amat for atom %d before coords set.\n",
135 >             index );
136 >    painCave.isFatal = 1;
137 >    simError();
138 >  }
139  
54 –  the_A[1][0] = Amat[Ayx];
55 –  the_A[1][1] = Amat[Ayy];
56 –  the_A[1][2] = Amat[Ayz];
57 –
58 –  the_A[2][0] = Amat[Azx];
59 –  the_A[2][1] = Amat[Azy];
60 –  the_A[2][2] = Amat[Azz];
140   }
141  
142   void DirectionalAtom::printAmatIndex( void ){
143  
144 <  std::cerr << "Atom[" << index << "] index =>\n"
145 <            << "[ " << Axx << ", " << Axy << ", " << Axz << " ]\n"
146 <            << "[ " << Ayx << ", " << Ayy << ", " << Ayz << " ]\n"
147 <            << "[ " << Azx << ", " << Azy << ", " << Azz << " ]\n";
144 >  if( hasCoords ){
145 >    std::cerr << "Atom[" << index << "] index =>\n"
146 >              << "[ " << Axx << ", " << Axy << ", " << Axz << " ]\n"
147 >              << "[ " << Ayx << ", " << Ayy << ", " << Ayz << " ]\n"
148 >              << "[ " << Azx << ", " << Azy << ", " << Azz << " ]\n";
149 >  }
150 >  else{
151 >    
152 >    sprintf( painCave.errMsg,
153 >             "Attempt to print Amat indices for atom %d before coords set.\n",
154 >             index );
155 >    painCave.isFatal = 1;
156 >    simError();
157 >  }
158   }
159  
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 83 | Line 172 | void DirectionalAtom::getQ( double q[4] ){
172    double t, s;
173    double ad1, ad2, ad3;
174  
175 <  t = Amat[Axx] + Amat[Ayy] + Amat[Azz] + 1.0;
87 <  if( t > 0.0 ){
175 >  if( hasCoords ){
176      
177 <    s = 0.5 / sqrt( t );
178 <    q[0] = 0.25 / s;
91 <    q[1] = (Amat[Ayz] - Amat[Azy]) * s;
92 <    q[2] = (Amat[Azx] - Amat[Axz]) * s;
93 <    q[3] = (Amat[Axy] - Amat[Ayx]) * s;
94 <  }
95 <  else{
96 <    
97 <    ad1 = fabs( Amat[Axx] );
98 <    ad2 = fabs( Amat[Ayy] );
99 <    ad3 = fabs( Amat[Azz] );
100 <
101 <    if( ad1 >= ad2 && ad1 >= ad3 ){
177 >    t = Amat[Axx] + Amat[Ayy] + Amat[Azz] + 1.0;
178 >    if( t > 0.0 ){
179        
180 <      s = 2.0 * sqrt( 1.0 + Amat[Axx] - Amat[Ayy] - Amat[Azz] );
181 <      q[0] = (Amat[Ayz] + Amat[Azy]) / s;
182 <      q[1] = 0.5 / s;
183 <      q[2] = (Amat[Axy] + Amat[Ayx]) / s;
184 <      q[3] = (Amat[Axz] + Amat[Azx]) / s;
180 >      s = 0.5 / sqrt( t );
181 >      q[0] = 0.25 / s;
182 >      q[1] = (Amat[Ayz] - Amat[Azy]) * s;
183 >      q[2] = (Amat[Azx] - Amat[Axz]) * s;
184 >      q[3] = (Amat[Axy] - Amat[Ayx]) * s;
185      }
109 –    else if( ad2 >= ad1 && ad2 >= ad3 ){
110 –
111 –      s = sqrt( 1.0 + Amat[Ayy] - Amat[Axx] - Amat[Azz] ) * 2.0;
112 –      q[0] = (Amat[Axz] + Amat[Azx]) / s;
113 –      q[1] = (Amat[Axy] + Amat[Ayx]) / s;
114 –      q[2] = 0.5 / s;
115 –      q[3] = (Amat[Ayz] + Amat[Azy]) / s;
116 –    }
186      else{
187        
188 <      s = sqrt( 1.0 + Amat[Azz] - Amat[Axx] - Amat[Ayy] ) * 2.0;
189 <      q[0] = (Amat[Axy] + Amat[Ayx]) / s;
190 <      q[1] = (Amat[Axz] + Amat[Azx]) / s;
191 <      q[2] = (Amat[Ayz] + Amat[Azy]) / s;
192 <      q[3] = 0.5 / s;
188 >      ad1 = fabs( Amat[Axx] );
189 >      ad2 = fabs( Amat[Ayy] );
190 >      ad3 = fabs( Amat[Azz] );
191 >      
192 >      if( ad1 >= ad2 && ad1 >= ad3 ){
193 >        
194 >        s = 2.0 * sqrt( 1.0 + Amat[Axx] - Amat[Ayy] - Amat[Azz] );
195 >        q[0] = (Amat[Ayz] + Amat[Azy]) / s;
196 >        q[1] = 0.5 / s;
197 >        q[2] = (Amat[Axy] + Amat[Ayx]) / s;
198 >        q[3] = (Amat[Axz] + Amat[Azx]) / s;
199 >      }
200 >      else if( ad2 >= ad1 && ad2 >= ad3 ){
201 >        
202 >        s = sqrt( 1.0 + Amat[Ayy] - Amat[Axx] - Amat[Azz] ) * 2.0;
203 >        q[0] = (Amat[Axz] + Amat[Azx]) / s;
204 >        q[1] = (Amat[Axy] + Amat[Ayx]) / s;
205 >        q[2] = 0.5 / s;
206 >        q[3] = (Amat[Ayz] + Amat[Azy]) / s;
207 >      }
208 >      else{
209 >        
210 >        s = sqrt( 1.0 + Amat[Azz] - Amat[Axx] - Amat[Ayy] ) * 2.0;
211 >        q[0] = (Amat[Axy] + Amat[Ayx]) / s;
212 >        q[1] = (Amat[Axz] + Amat[Azx]) / s;
213 >        q[2] = (Amat[Ayz] + Amat[Azy]) / s;
214 >        q[3] = 0.5 / s;
215 >      }
216      }
217    }
218 +  else{
219 +    
220 +    sprintf( painCave.errMsg,
221 +             "Attempt to get Q for atom %d before coords set.\n",
222 +             index );
223 +    painCave.isFatal = 1;
224 +    simError();
225 +  }
226   }
227  
228 + void DirectionalAtom::setUnitFrameFromEuler(double phi,
229 +                                            double theta,
230 +                                            double psi) {
231  
232 < void DirectionalAtom::setEuler( double phi, double theta, double psi ){
232 >  double myA[3][3];
233 >  double uFrame[3][3];
234 >  double len;
235 >  int i, j;
236    
237 <  Amat[Axx] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
238 <  Amat[Axy] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
239 <  Amat[Axz] = sin(theta) * sin(psi);
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 <  Amat[Ayx] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
242 <  Amat[Ayy] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
243 <  Amat[Ayz] = sin(theta) * cos(psi);
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 <  Amat[Azx] = sin(phi) * sin(theta);
252 <  Amat[Azy] = -cos(phi) * sin(theta);
253 <  Amat[Azz] = cos(theta);
251 >  for (i=0; i < 3; i++)
252 >    for (j=0; j < 3; j++)
253 >      uFrame[i][j] = 0.0;
254  
255 <  this->updateU();
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 ){
282 +    Amat[Axx] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
283 +    Amat[Axy] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
284 +    Amat[Axz] = sin(theta) * sin(psi);
285 +    
286 +    Amat[Ayx] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
287 +    Amat[Ayy] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
288 +    Amat[Ayz] = sin(theta) * cos(psi);
289 +    
290 +    Amat[Azx] = sin(phi) * sin(theta);
291 +    Amat[Azy] = -cos(phi) * sin(theta);
292 +    Amat[Azz] = cos(theta);
293 +    
294 +    this->updateU();
295 +  }
296 +  else{
297 +    
298 +    sprintf( painCave.errMsg,
299 +             "Attempt to set Euler angles for atom %d before coords set.\n",
300 +             index );
301 +    painCave.isFatal = 1;
302 +    simError();
303 +  }
304 + }
305  
306 +
307   void DirectionalAtom::lab2Body( double r[3] ){
308  
309    double rl[3]; // the lab frame vector
310    
311 <  rl[0] = r[0];
312 <  rl[1] = r[1];
313 <  rl[2] = r[2];
311 >  if( hasCoords ){
312 >    rl[0] = r[0];
313 >    rl[1] = r[1];
314 >    rl[2] = r[2];
315 >    
316 >    r[0] = (Amat[Axx] * rl[0]) + (Amat[Axy] * rl[1]) + (Amat[Axz] * rl[2]);
317 >    r[1] = (Amat[Ayx] * rl[0]) + (Amat[Ayy] * rl[1]) + (Amat[Ayz] * rl[2]);
318 >    r[2] = (Amat[Azx] * rl[0]) + (Amat[Azy] * rl[1]) + (Amat[Azz] * rl[2]);
319 >  }
320 >  else{
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  
155 –  r[0] = (Amat[Axx] * rl[0]) + (Amat[Axy] * rl[1]) + (Amat[Axz] * rl[2]);
156 –  r[1] = (Amat[Ayx] * rl[0]) + (Amat[Ayy] * rl[1]) + (Amat[Ayz] * rl[2]);
157 –  r[2] = (Amat[Azx] * rl[0]) + (Amat[Azy] * rl[1]) + (Amat[Azz] * rl[2]);
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();
363 +  }
364 +
365 + }
366 +
367 +
368   void DirectionalAtom::body2Lab( double r[3] ){
369  
370    double rb[3]; // the body frame vector
371    
372 <  rb[0] = r[0];
373 <  rb[1] = r[1];
374 <  rb[2] = r[2];
375 <
376 <  r[0] = (Amat[Axx] * rb[0]) + (Amat[Ayx] * rb[1]) + (Amat[Azx] * rb[2]);
377 <  r[1] = (Amat[Axy] * rb[0]) + (Amat[Ayy] * rb[1]) + (Amat[Azy] * rb[2]);
378 <  r[2] = (Amat[Axz] * rb[0]) + (Amat[Ayz] * rb[1]) + (Amat[Azz] * rb[2]);
372 >  if( hasCoords ){
373 >    rb[0] = r[0];
374 >    rb[1] = r[1];
375 >    rb[2] = r[2];
376 >    
377 >    r[0] = (Amat[Axx] * rb[0]) + (Amat[Ayx] * rb[1]) + (Amat[Azx] * rb[2]);
378 >    r[1] = (Amat[Axy] * rb[0]) + (Amat[Ayy] * rb[1]) + (Amat[Azy] * rb[2]);
379 >    r[2] = (Amat[Axz] * rb[0]) + (Amat[Ayz] * rb[1]) + (Amat[Azz] * rb[2]);
380 >  }
381 >  else{
382 >    
383 >    sprintf( painCave.errMsg,
384 >             "Attempt to convert body2lab for atom %d before coords set.\n",
385 >             index );
386 >    painCave.isFatal = 1;
387 >    simError();
388 >  }
389   }
390  
391   void DirectionalAtom::updateU( void ){
392  
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 >  if( hasCoords ){
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 >    
403 >    sprintf( painCave.errMsg,
404 >             "Attempt to updateU for atom %d before coords set.\n",
405 >             index );
406 >    painCave.isFatal = 1;
407 >    simError();
408 >  }
409   }
410  
411 + void DirectionalAtom::getJ( double theJ[3] ){
412 +  
413 +  theJ[0] = jx;
414 +  theJ[1] = jy;
415 +  theJ[2] = jz;
416 + }
417 +
418 + void DirectionalAtom::setJ( double theJ[3] ){
419 +  
420 +  jx = theJ[0];
421 +  jy = theJ[1];
422 +  jz = theJ[2];
423 + }
424 +
425 + void DirectionalAtom::getTrq( double theT[3] ){
426 +  
427 +  if( hasCoords ){
428 +    theT[0] = trq[offsetX];
429 +    theT[1] = trq[offsetY];
430 +    theT[2] = trq[offsetZ];
431 +  }
432 +  else{
433 +    
434 +    sprintf( painCave.errMsg,
435 +             "Attempt to get Trq for atom %d before coords set.\n",
436 +             index );
437 +    painCave.isFatal = 1;
438 +    simError();
439 +  }
440 + }
441 +
442 + void DirectionalAtom::addTrq( double theT[3] ){
443 +  
444 +  if( hasCoords ){
445 +    trq[offsetX] += theT[0];
446 +    trq[offsetY] += theT[1];
447 +    trq[offsetZ] += theT[2];
448 +  }
449 +  else{
450 +    
451 +    sprintf( painCave.errMsg,
452 +             "Attempt to add Trq for atom %d before coords set.\n",
453 +             index );
454 +    painCave.isFatal = 1;
455 +    simError();
456 +  }
457 + }
458 +
459 +
460 + void DirectionalAtom::getI( double the_I[3][3] ){
461 +  
462 +  the_I[0][0] = Ixx;
463 +  the_I[0][1] = Ixy;
464 +  the_I[0][2] = Ixz;
465 +
466 +  the_I[1][0] = Iyx;
467 +  the_I[1][1] = Iyy;
468 +  the_I[1][2] = Iyz;
469 +
470 +  the_I[2][0] = Izx;
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