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 599 by mmeineke, Mon Jul 14 21:48:43 2003 UTC vs.
Revision 787 by mmeineke, Thu Sep 25 19:27:15 2003 UTC

# Line 1 | Line 1
1   #include <cmath>
2  
3   #include "Atom.hpp"
4 + #include "simError.h"
5  
6 + void DirectionalAtom::zeroForces() {
7 +  if( hasCoords ){
8 +    frc[offsetX] = 0.0;
9 +    frc[offsetY] = 0.0;
10 +    frc[offsetZ] = 0.0;
11 +    
12 +    trq[offsetX] = 0.0;
13 +    trq[offsetY] = 0.0;
14 +    trq[offsetZ] = 0.0;
15 +  }
16 +  else{
17 +    
18 +    sprintf( painCave.errMsg,
19 +             "Attempt to zero frc and trq for atom %d before coords set.\n",
20 +             index );
21 +    painCave.isFatal = 1;
22 +    simError();
23 +  }
24 + }
25  
26 + void DirectionalAtom::setCoords(void){
27  
28 +  if( myConfig->isAllocated() ){
29 +
30 +    myConfig->getAtomPointers( index,
31 +                     &pos,
32 +                     &vel,
33 +                     &frc,
34 +                     &trq,
35 +                     &Amat,
36 +                     &mu,  
37 +                     &ul );
38 +  }
39 +  else{
40 +    sprintf( painCave.errMsg,
41 +             "Attempted to set Atom %d  coordinates with an unallocated "
42 +             "SimState object.\n", index );
43 +    painCave.isFatal = 1;
44 +    simError();
45 +  }
46 +
47 +  hasCoords = true;
48 +
49 +  *mu = myMu;
50 +
51 + }
52 +
53 + double DirectionalAtom::getMu( void ) {
54 +
55 +  if( hasCoords ){
56 +    return *mu;
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 = the_mu;
68 +    myMu = the_mu;
69 +  }
70 +  else{
71 +    myMu = the_mu;
72 +  }
73 + }
74 +
75   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];
76  
77 <  this->updateU();  
77 >  if( hasCoords ){
78 >    Amat[Axx] = the_A[0][0]; Amat[Axy] = the_A[0][1]; Amat[Axz] = the_A[0][2];
79 >    Amat[Ayx] = the_A[1][0]; Amat[Ayy] = the_A[1][1]; Amat[Ayz] = the_A[1][2];
80 >    Amat[Azx] = the_A[2][0]; Amat[Azy] = the_A[2][1]; Amat[Azz] = the_A[2][2];
81 >    
82 >    this->updateU();  
83 >  }
84 >  else{
85 >    
86 >    sprintf( painCave.errMsg,
87 >             "Attempt to set Amat for atom %d before coords set.\n",
88 >             index );
89 >    painCave.isFatal = 1;
90 >    simError();
91 >  }
92   }
93  
94   void DirectionalAtom::setI( double the_I[3][3] ){
# Line 24 | Line 102 | void DirectionalAtom::setQ( double the_q[4] ){
102  
103    double q0Sqr, q1Sqr, q2Sqr, q3Sqr;
104  
105 <  q0Sqr = the_q[0] * the_q[0];
106 <  q1Sqr = the_q[1] * the_q[1];
107 <  q2Sqr = the_q[2] * the_q[2];
108 <  q3Sqr = the_q[3] * the_q[3];
109 <  
105 >  if( hasCoords ){
106 >    q0Sqr = the_q[0] * the_q[0];
107 >    q1Sqr = the_q[1] * the_q[1];
108 >    q2Sqr = the_q[2] * the_q[2];
109 >    q3Sqr = the_q[3] * the_q[3];
110 >    
111 >    
112 >    Amat[Axx] = q0Sqr + q1Sqr - q2Sqr - q3Sqr;
113 >    Amat[Axy] = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] );
114 >    Amat[Axz] = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] );
115 >    
116 >    Amat[Ayx] = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] );
117 >    Amat[Ayy] = q0Sqr - q1Sqr + q2Sqr - q3Sqr;
118 >    Amat[Ayz] = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] );
119 >    
120 >    Amat[Azx] = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] );
121 >    Amat[Azy] = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] );
122 >    Amat[Azz] = q0Sqr - q1Sqr -q2Sqr +q3Sqr;
123 >    
124 >    this->updateU();
125 >  }
126 >  else{
127 >    
128 >    sprintf( painCave.errMsg,
129 >             "Attempt to set Q for atom %d before coords set.\n",
130 >             index );
131 >    painCave.isFatal = 1;
132 >    simError();
133 >  }
134  
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();
135   }
136  
137   void DirectionalAtom::getA( double the_A[3][3] ){
138    
139 <  the_A[0][0] = Amat[Axx];
140 <  the_A[0][1] = Amat[Axy];
141 <  the_A[0][2] = Amat[Axz];
139 >  if( hasCoords ){
140 >    the_A[0][0] = Amat[Axx];
141 >    the_A[0][1] = Amat[Axy];
142 >    the_A[0][2] = Amat[Axz];
143 >    
144 >    the_A[1][0] = Amat[Ayx];
145 >    the_A[1][1] = Amat[Ayy];
146 >    the_A[1][2] = Amat[Ayz];
147 >    
148 >    the_A[2][0] = Amat[Azx];
149 >    the_A[2][1] = Amat[Azy];
150 >    the_A[2][2] = Amat[Azz];
151 >  }
152 >  else{
153 >    
154 >    sprintf( painCave.errMsg,
155 >             "Attempt to get Amat for atom %d before coords set.\n",
156 >             index );
157 >    painCave.isFatal = 1;
158 >    simError();
159 >  }
160  
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];
161   }
162  
163   void DirectionalAtom::printAmatIndex( void ){
164  
165 <  std::cerr << "Atom[" << index << "] index =>\n"
166 <            << "[ " << Axx << ", " << Axy << ", " << Axz << " ]\n"
167 <            << "[ " << Ayx << ", " << Ayy << ", " << Ayz << " ]\n"
168 <            << "[ " << Azx << ", " << Azy << ", " << Azz << " ]\n";
165 >  if( hasCoords ){
166 >    std::cerr << "Atom[" << index << "] index =>\n"
167 >              << "[ " << Axx << ", " << Axy << ", " << Axz << " ]\n"
168 >              << "[ " << Ayx << ", " << Ayy << ", " << Ayz << " ]\n"
169 >              << "[ " << Azx << ", " << Azy << ", " << Azz << " ]\n";
170 >  }
171 >  else{
172 >    
173 >    sprintf( painCave.errMsg,
174 >             "Attempt to print Amat indices for atom %d before coords set.\n",
175 >             index );
176 >    painCave.isFatal = 1;
177 >    simError();
178 >  }
179   }
180  
181  
# Line 83 | Line 193 | void DirectionalAtom::getQ( double q[4] ){
193    double t, s;
194    double ad1, ad2, ad3;
195  
196 <  t = Amat[Axx] + Amat[Ayy] + Amat[Azz] + 1.0;
87 <  if( t > 0.0 ){
196 >  if( hasCoords ){
197      
198 <    s = 0.5 / sqrt( t );
199 <    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 ){
198 >    t = Amat[Axx] + Amat[Ayy] + Amat[Azz] + 1.0;
199 >    if( t > 0.0 ){
200        
201 <      s = 2.0 * sqrt( 1.0 + Amat[Axx] - Amat[Ayy] - Amat[Azz] );
202 <      q[0] = (Amat[Ayz] + Amat[Azy]) / s;
203 <      q[1] = 0.5 / s;
204 <      q[2] = (Amat[Axy] + Amat[Ayx]) / s;
205 <      q[3] = (Amat[Axz] + Amat[Azx]) / s;
201 >      s = 0.5 / sqrt( t );
202 >      q[0] = 0.25 / s;
203 >      q[1] = (Amat[Ayz] - Amat[Azy]) * s;
204 >      q[2] = (Amat[Azx] - Amat[Axz]) * s;
205 >      q[3] = (Amat[Axy] - Amat[Ayx]) * s;
206      }
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 –    }
207      else{
208        
209 <      s = sqrt( 1.0 + Amat[Azz] - Amat[Axx] - Amat[Ayy] ) * 2.0;
210 <      q[0] = (Amat[Axy] + Amat[Ayx]) / s;
211 <      q[1] = (Amat[Axz] + Amat[Azx]) / s;
212 <      q[2] = (Amat[Ayz] + Amat[Azy]) / s;
213 <      q[3] = 0.5 / s;
209 >      ad1 = fabs( Amat[Axx] );
210 >      ad2 = fabs( Amat[Ayy] );
211 >      ad3 = fabs( Amat[Azz] );
212 >      
213 >      if( ad1 >= ad2 && ad1 >= ad3 ){
214 >        
215 >        s = 2.0 * sqrt( 1.0 + Amat[Axx] - Amat[Ayy] - Amat[Azz] );
216 >        q[0] = (Amat[Ayz] + Amat[Azy]) / s;
217 >        q[1] = 0.5 / s;
218 >        q[2] = (Amat[Axy] + Amat[Ayx]) / s;
219 >        q[3] = (Amat[Axz] + Amat[Azx]) / s;
220 >      }
221 >      else if( ad2 >= ad1 && ad2 >= ad3 ){
222 >        
223 >        s = sqrt( 1.0 + Amat[Ayy] - Amat[Axx] - Amat[Azz] ) * 2.0;
224 >        q[0] = (Amat[Axz] + Amat[Azx]) / s;
225 >        q[1] = (Amat[Axy] + Amat[Ayx]) / s;
226 >        q[2] = 0.5 / s;
227 >        q[3] = (Amat[Ayz] + Amat[Azy]) / s;
228 >      }
229 >      else{
230 >        
231 >        s = sqrt( 1.0 + Amat[Azz] - Amat[Axx] - Amat[Ayy] ) * 2.0;
232 >        q[0] = (Amat[Axy] + Amat[Ayx]) / s;
233 >        q[1] = (Amat[Axz] + Amat[Azx]) / s;
234 >        q[2] = (Amat[Ayz] + Amat[Azy]) / s;
235 >        q[3] = 0.5 / s;
236 >      }
237      }
238    }
239 +  else{
240 +    
241 +    sprintf( painCave.errMsg,
242 +             "Attempt to get Q for atom %d before coords set.\n",
243 +             index );
244 +    painCave.isFatal = 1;
245 +    simError();
246 +  }
247   }
248  
249  
250   void DirectionalAtom::setEuler( double phi, double theta, double psi ){
251    
252 <  Amat[Axx] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
253 <  Amat[Axy] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
254 <  Amat[Axz] = sin(theta) * sin(psi);
255 <  
256 <  Amat[Ayx] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
257 <  Amat[Ayy] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
258 <  Amat[Ayz] = sin(theta) * cos(psi);
259 <
260 <  Amat[Azx] = sin(phi) * sin(theta);
261 <  Amat[Azy] = -cos(phi) * sin(theta);
262 <  Amat[Azz] = cos(theta);
263 <
264 <  this->updateU();
252 >  if( hasCoords ){
253 >    Amat[Axx] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
254 >    Amat[Axy] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
255 >    Amat[Axz] = sin(theta) * sin(psi);
256 >    
257 >    Amat[Ayx] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
258 >    Amat[Ayy] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
259 >    Amat[Ayz] = sin(theta) * cos(psi);
260 >    
261 >    Amat[Azx] = sin(phi) * sin(theta);
262 >    Amat[Azy] = -cos(phi) * sin(theta);
263 >    Amat[Azz] = cos(theta);
264 >    
265 >    this->updateU();
266 >  }
267 >  else{
268 >    
269 >    sprintf( painCave.errMsg,
270 >             "Attempt to set Euler angles for atom %d before coords set.\n",
271 >             index );
272 >    painCave.isFatal = 1;
273 >    simError();
274 >  }
275   }
276  
277  
# Line 148 | Line 279 | void DirectionalAtom::lab2Body( double r[3] ){
279  
280    double rl[3]; // the lab frame vector
281    
282 <  rl[0] = r[0];
283 <  rl[1] = r[1];
284 <  rl[2] = r[2];
282 >  if( hasCoords ){
283 >    rl[0] = r[0];
284 >    rl[1] = r[1];
285 >    rl[2] = r[2];
286 >    
287 >    r[0] = (Amat[Axx] * rl[0]) + (Amat[Axy] * rl[1]) + (Amat[Axz] * rl[2]);
288 >    r[1] = (Amat[Ayx] * rl[0]) + (Amat[Ayy] * rl[1]) + (Amat[Ayz] * rl[2]);
289 >    r[2] = (Amat[Azx] * rl[0]) + (Amat[Azy] * rl[1]) + (Amat[Azz] * rl[2]);
290 >  }
291 >  else{
292 >    
293 >    sprintf( painCave.errMsg,
294 >             "Attempt to convert lab2body for atom %d before coords set.\n",
295 >             index );
296 >    painCave.isFatal = 1;
297 >    simError();
298 >  }
299  
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]);
300   }
301  
302   void DirectionalAtom::body2Lab( double r[3] ){
303  
304    double rb[3]; // the body frame vector
305    
306 <  rb[0] = r[0];
307 <  rb[1] = r[1];
308 <  rb[2] = r[2];
309 <
310 <  r[0] = (Amat[Axx] * rb[0]) + (Amat[Ayx] * rb[1]) + (Amat[Azx] * rb[2]);
311 <  r[1] = (Amat[Axy] * rb[0]) + (Amat[Ayy] * rb[1]) + (Amat[Azy] * rb[2]);
312 <  r[2] = (Amat[Axz] * rb[0]) + (Amat[Ayz] * rb[1]) + (Amat[Azz] * rb[2]);
306 >  if( hasCoords ){
307 >    rb[0] = r[0];
308 >    rb[1] = r[1];
309 >    rb[2] = r[2];
310 >    
311 >    r[0] = (Amat[Axx] * rb[0]) + (Amat[Ayx] * rb[1]) + (Amat[Azx] * rb[2]);
312 >    r[1] = (Amat[Axy] * rb[0]) + (Amat[Ayy] * rb[1]) + (Amat[Azy] * rb[2]);
313 >    r[2] = (Amat[Axz] * rb[0]) + (Amat[Ayz] * rb[1]) + (Amat[Azz] * rb[2]);
314 >  }
315 >  else{
316 >    
317 >    sprintf( painCave.errMsg,
318 >             "Attempt to convert body2lab for atom %d before coords set.\n",
319 >             index );
320 >    painCave.isFatal = 1;
321 >    simError();
322 >  }
323   }
324  
325   void DirectionalAtom::updateU( void ){
326  
327 <  ul[offsetX] = (Amat[Axx] * sux) + (Amat[Ayx] * suy) + (Amat[Azx] * suz);
328 <  ul[offsetY] = (Amat[Axy] * sux) + (Amat[Ayy] * suy) + (Amat[Azy] * suz);
329 <  ul[offsetZ] = (Amat[Axz] * sux) + (Amat[Ayz] * suy) + (Amat[Azz] * suz);
327 >  if( hasCoords ){
328 >    ul[offsetX] = (Amat[Axx] * sux) + (Amat[Ayx] * suy) + (Amat[Azx] * suz);
329 >    ul[offsetY] = (Amat[Axy] * sux) + (Amat[Ayy] * suy) + (Amat[Azy] * suz);
330 >    ul[offsetZ] = (Amat[Axz] * sux) + (Amat[Ayz] * suy) + (Amat[Azz] * suz);
331 >  }
332 >  else{
333 >    
334 >    sprintf( painCave.errMsg,
335 >             "Attempt to updateU for atom %d before coords set.\n",
336 >             index );
337 >    painCave.isFatal = 1;
338 >    simError();
339 >  }
340   }
341  
342   void DirectionalAtom::getJ( double theJ[3] ){
# Line 193 | Line 355 | void DirectionalAtom::getTrq( double theT[3] ){
355  
356   void DirectionalAtom::getTrq( double theT[3] ){
357    
358 <  theT[0] = trq[offsetX];
359 <  theT[1] = trq[offsetY];
360 <  theT[2] = trq[offsetZ];
358 >  if( hasCoords ){
359 >    theT[0] = trq[offsetX];
360 >    theT[1] = trq[offsetY];
361 >    theT[2] = trq[offsetZ];
362 >  }
363 >  else{
364 >    
365 >    sprintf( painCave.errMsg,
366 >             "Attempt to get Trq for atom %d before coords set.\n",
367 >             index );
368 >    painCave.isFatal = 1;
369 >    simError();
370 >  }
371   }
372  
373   void DirectionalAtom::addTrq( double theT[3] ){
374    
375 <  trq[offsetX] += theT[0];
376 <  trq[offsetY] += theT[1];
377 <  trq[offsetZ] += theT[2];
375 >  if( hasCoords ){
376 >    trq[offsetX] += theT[0];
377 >    trq[offsetY] += theT[1];
378 >    trq[offsetZ] += theT[2];
379 >  }
380 >  else{
381 >    
382 >    sprintf( painCave.errMsg,
383 >             "Attempt to add Trq for atom %d before coords set.\n",
384 >             index );
385 >    painCave.isFatal = 1;
386 >    simError();
387 >  }
388   }
389  
390  

Diff Legend

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