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 670 by mmeineke, Thu Aug 7 21:47:18 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 + double DirectionalAtom::getMu( void ) {
27  
28 +  if( hasCoords ){
29 +    return mu[index];
30 +  }
31 +  else{
32 +    
33 +    sprintf( painCave.errMsg,
34 +             "Attempt to get Mu for atom %d before coords set.\n",
35 +             index );
36 +    painCave.isFatal = 1;
37 +    simError();
38 +  }
39 +  return 0;
40 + }
41 +
42 + void DirectionalAtom::setMu( double the_mu ) {
43 +
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 +  }
55 + }
56 +
57   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];
58  
59 <  this->updateU();  
59 >  if( hasCoords ){
60 >    Amat[Axx] = the_A[0][0]; Amat[Axy] = the_A[0][1]; Amat[Axz] = the_A[0][2];
61 >    Amat[Ayx] = the_A[1][0]; Amat[Ayy] = the_A[1][1]; Amat[Ayz] = the_A[1][2];
62 >    Amat[Azx] = the_A[2][0]; Amat[Azy] = the_A[2][1]; Amat[Azz] = the_A[2][2];
63 >    
64 >    this->updateU();  
65 >  }
66 >  else{
67 >    
68 >    sprintf( painCave.errMsg,
69 >             "Attempt to set Amat for atom %d before coords set.\n",
70 >             index );
71 >    painCave.isFatal = 1;
72 >    simError();
73 >  }
74   }
75  
76   void DirectionalAtom::setI( double the_I[3][3] ){
# Line 24 | Line 84 | void DirectionalAtom::setQ( double the_q[4] ){
84  
85    double q0Sqr, q1Sqr, q2Sqr, q3Sqr;
86  
87 <  q0Sqr = the_q[0] * the_q[0];
88 <  q1Sqr = the_q[1] * the_q[1];
89 <  q2Sqr = the_q[2] * the_q[2];
90 <  q3Sqr = the_q[3] * the_q[3];
91 <  
87 >  if( hasCoords ){
88 >    q0Sqr = the_q[0] * the_q[0];
89 >    q1Sqr = the_q[1] * the_q[1];
90 >    q2Sqr = the_q[2] * the_q[2];
91 >    q3Sqr = the_q[3] * the_q[3];
92 >    
93 >    
94 >    Amat[Axx] = q0Sqr + q1Sqr - q2Sqr - q3Sqr;
95 >    Amat[Axy] = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] );
96 >    Amat[Axz] = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] );
97 >    
98 >    Amat[Ayx] = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] );
99 >    Amat[Ayy] = q0Sqr - q1Sqr + q2Sqr - q3Sqr;
100 >    Amat[Ayz] = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] );
101 >    
102 >    Amat[Azx] = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] );
103 >    Amat[Azy] = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] );
104 >    Amat[Azz] = q0Sqr - q1Sqr -q2Sqr +q3Sqr;
105 >    
106 >    this->updateU();
107 >  }
108 >  else{
109 >    
110 >    sprintf( painCave.errMsg,
111 >             "Attempt to set Q for atom %d before coords set.\n",
112 >             index );
113 >    painCave.isFatal = 1;
114 >    simError();
115 >  }
116  
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();
117   }
118  
119   void DirectionalAtom::getA( double the_A[3][3] ){
120    
121 <  the_A[0][0] = Amat[Axx];
122 <  the_A[0][1] = Amat[Axy];
123 <  the_A[0][2] = Amat[Axz];
121 >  if( hasCoords ){
122 >    the_A[0][0] = Amat[Axx];
123 >    the_A[0][1] = Amat[Axy];
124 >    the_A[0][2] = Amat[Axz];
125 >    
126 >    the_A[1][0] = Amat[Ayx];
127 >    the_A[1][1] = Amat[Ayy];
128 >    the_A[1][2] = Amat[Ayz];
129 >    
130 >    the_A[2][0] = Amat[Azx];
131 >    the_A[2][1] = Amat[Azy];
132 >    the_A[2][2] = Amat[Azz];
133 >  }
134 >  else{
135 >    
136 >    sprintf( painCave.errMsg,
137 >             "Attempt to get Amat for atom %d before coords set.\n",
138 >             index );
139 >    painCave.isFatal = 1;
140 >    simError();
141 >  }
142  
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];
143   }
144  
145   void DirectionalAtom::printAmatIndex( void ){
146  
147 <  std::cerr << "Atom[" << index << "] index =>\n"
148 <            << "[ " << Axx << ", " << Axy << ", " << Axz << " ]\n"
149 <            << "[ " << Ayx << ", " << Ayy << ", " << Ayz << " ]\n"
150 <            << "[ " << Azx << ", " << Azy << ", " << Azz << " ]\n";
147 >  if( hasCoords ){
148 >    std::cerr << "Atom[" << index << "] index =>\n"
149 >              << "[ " << Axx << ", " << Axy << ", " << Axz << " ]\n"
150 >              << "[ " << Ayx << ", " << Ayy << ", " << Ayz << " ]\n"
151 >              << "[ " << Azx << ", " << Azy << ", " << Azz << " ]\n";
152 >  }
153 >  else{
154 >    
155 >    sprintf( painCave.errMsg,
156 >             "Attempt to print Amat indices for atom %d before coords set.\n",
157 >             index );
158 >    painCave.isFatal = 1;
159 >    simError();
160 >  }
161   }
162  
163  
# Line 83 | Line 175 | void DirectionalAtom::getQ( double q[4] ){
175    double t, s;
176    double ad1, ad2, ad3;
177  
178 <  t = Amat[Axx] + Amat[Ayy] + Amat[Azz] + 1.0;
87 <  if( t > 0.0 ){
178 >  if( hasCoords ){
179      
180 <    s = 0.5 / sqrt( t );
181 <    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 ){
180 >    t = Amat[Axx] + Amat[Ayy] + Amat[Azz] + 1.0;
181 >    if( t > 0.0 ){
182        
183 <      s = 2.0 * sqrt( 1.0 + Amat[Axx] - Amat[Ayy] - Amat[Azz] );
184 <      q[0] = (Amat[Ayz] + Amat[Azy]) / s;
185 <      q[1] = 0.5 / s;
186 <      q[2] = (Amat[Axy] + Amat[Ayx]) / s;
187 <      q[3] = (Amat[Axz] + Amat[Azx]) / s;
183 >      s = 0.5 / sqrt( t );
184 >      q[0] = 0.25 / s;
185 >      q[1] = (Amat[Ayz] - Amat[Azy]) * s;
186 >      q[2] = (Amat[Azx] - Amat[Axz]) * s;
187 >      q[3] = (Amat[Axy] - Amat[Ayx]) * s;
188      }
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 –    }
189      else{
190        
191 <      s = sqrt( 1.0 + Amat[Azz] - Amat[Axx] - Amat[Ayy] ) * 2.0;
192 <      q[0] = (Amat[Axy] + Amat[Ayx]) / s;
193 <      q[1] = (Amat[Axz] + Amat[Azx]) / s;
194 <      q[2] = (Amat[Ayz] + Amat[Azy]) / s;
195 <      q[3] = 0.5 / s;
191 >      ad1 = fabs( Amat[Axx] );
192 >      ad2 = fabs( Amat[Ayy] );
193 >      ad3 = fabs( Amat[Azz] );
194 >      
195 >      if( ad1 >= ad2 && ad1 >= ad3 ){
196 >        
197 >        s = 2.0 * sqrt( 1.0 + Amat[Axx] - Amat[Ayy] - Amat[Azz] );
198 >        q[0] = (Amat[Ayz] + Amat[Azy]) / s;
199 >        q[1] = 0.5 / s;
200 >        q[2] = (Amat[Axy] + Amat[Ayx]) / s;
201 >        q[3] = (Amat[Axz] + Amat[Azx]) / s;
202 >      }
203 >      else if( ad2 >= ad1 && ad2 >= ad3 ){
204 >        
205 >        s = sqrt( 1.0 + Amat[Ayy] - Amat[Axx] - Amat[Azz] ) * 2.0;
206 >        q[0] = (Amat[Axz] + Amat[Azx]) / s;
207 >        q[1] = (Amat[Axy] + Amat[Ayx]) / s;
208 >        q[2] = 0.5 / s;
209 >        q[3] = (Amat[Ayz] + Amat[Azy]) / s;
210 >      }
211 >      else{
212 >        
213 >        s = sqrt( 1.0 + Amat[Azz] - Amat[Axx] - Amat[Ayy] ) * 2.0;
214 >        q[0] = (Amat[Axy] + Amat[Ayx]) / s;
215 >        q[1] = (Amat[Axz] + Amat[Azx]) / s;
216 >        q[2] = (Amat[Ayz] + Amat[Azy]) / s;
217 >        q[3] = 0.5 / s;
218 >      }
219      }
220    }
221 +  else{
222 +    
223 +    sprintf( painCave.errMsg,
224 +             "Attempt to get Q for atom %d before coords set.\n",
225 +             index );
226 +    painCave.isFatal = 1;
227 +    simError();
228 +  }
229   }
230  
231  
232   void DirectionalAtom::setEuler( double phi, double theta, double psi ){
233    
234 <  Amat[Axx] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
235 <  Amat[Axy] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
236 <  Amat[Axz] = sin(theta) * sin(psi);
237 <  
238 <  Amat[Ayx] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
239 <  Amat[Ayy] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
240 <  Amat[Ayz] = sin(theta) * cos(psi);
241 <
242 <  Amat[Azx] = sin(phi) * sin(theta);
243 <  Amat[Azy] = -cos(phi) * sin(theta);
244 <  Amat[Azz] = cos(theta);
245 <
246 <  this->updateU();
234 >  if( hasCoords ){
235 >    Amat[Axx] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
236 >    Amat[Axy] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
237 >    Amat[Axz] = sin(theta) * sin(psi);
238 >    
239 >    Amat[Ayx] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
240 >    Amat[Ayy] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
241 >    Amat[Ayz] = sin(theta) * cos(psi);
242 >    
243 >    Amat[Azx] = sin(phi) * sin(theta);
244 >    Amat[Azy] = -cos(phi) * sin(theta);
245 >    Amat[Azz] = cos(theta);
246 >    
247 >    this->updateU();
248 >  }
249 >  else{
250 >    
251 >    sprintf( painCave.errMsg,
252 >             "Attempt to set Euler angles for atom %d before coords set.\n",
253 >             index );
254 >    painCave.isFatal = 1;
255 >    simError();
256 >  }
257   }
258  
259  
# Line 148 | Line 261 | void DirectionalAtom::lab2Body( double r[3] ){
261  
262    double rl[3]; // the lab frame vector
263    
264 <  rl[0] = r[0];
265 <  rl[1] = r[1];
266 <  rl[2] = r[2];
264 >  if( hasCoords ){
265 >    rl[0] = r[0];
266 >    rl[1] = r[1];
267 >    rl[2] = r[2];
268 >    
269 >    r[0] = (Amat[Axx] * rl[0]) + (Amat[Axy] * rl[1]) + (Amat[Axz] * rl[2]);
270 >    r[1] = (Amat[Ayx] * rl[0]) + (Amat[Ayy] * rl[1]) + (Amat[Ayz] * rl[2]);
271 >    r[2] = (Amat[Azx] * rl[0]) + (Amat[Azy] * rl[1]) + (Amat[Azz] * rl[2]);
272 >  }
273 >  else{
274 >    
275 >    sprintf( painCave.errMsg,
276 >             "Attempt to convert lab2body for atom %d before coords set.\n",
277 >             index );
278 >    painCave.isFatal = 1;
279 >    simError();
280 >  }
281  
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]);
282   }
283  
284   void DirectionalAtom::body2Lab( double r[3] ){
285  
286    double rb[3]; // the body frame vector
287    
288 <  rb[0] = r[0];
289 <  rb[1] = r[1];
290 <  rb[2] = r[2];
291 <
292 <  r[0] = (Amat[Axx] * rb[0]) + (Amat[Ayx] * rb[1]) + (Amat[Azx] * rb[2]);
293 <  r[1] = (Amat[Axy] * rb[0]) + (Amat[Ayy] * rb[1]) + (Amat[Azy] * rb[2]);
294 <  r[2] = (Amat[Axz] * rb[0]) + (Amat[Ayz] * rb[1]) + (Amat[Azz] * rb[2]);
288 >  if( hasCoords ){
289 >    rb[0] = r[0];
290 >    rb[1] = r[1];
291 >    rb[2] = r[2];
292 >    
293 >    r[0] = (Amat[Axx] * rb[0]) + (Amat[Ayx] * rb[1]) + (Amat[Azx] * rb[2]);
294 >    r[1] = (Amat[Axy] * rb[0]) + (Amat[Ayy] * rb[1]) + (Amat[Azy] * rb[2]);
295 >    r[2] = (Amat[Axz] * rb[0]) + (Amat[Ayz] * rb[1]) + (Amat[Azz] * rb[2]);
296 >  }
297 >  else{
298 >    
299 >    sprintf( painCave.errMsg,
300 >             "Attempt to convert body2lab for atom %d before coords set.\n",
301 >             index );
302 >    painCave.isFatal = 1;
303 >    simError();
304 >  }
305   }
306  
307   void DirectionalAtom::updateU( void ){
308  
309 <  ul[offsetX] = (Amat[Axx] * sux) + (Amat[Ayx] * suy) + (Amat[Azx] * suz);
310 <  ul[offsetY] = (Amat[Axy] * sux) + (Amat[Ayy] * suy) + (Amat[Azy] * suz);
311 <  ul[offsetZ] = (Amat[Axz] * sux) + (Amat[Ayz] * suy) + (Amat[Azz] * suz);
309 >  if( hasCoords ){
310 >    ul[offsetX] = (Amat[Axx] * sux) + (Amat[Ayx] * suy) + (Amat[Azx] * suz);
311 >    ul[offsetY] = (Amat[Axy] * sux) + (Amat[Ayy] * suy) + (Amat[Azy] * suz);
312 >    ul[offsetZ] = (Amat[Axz] * sux) + (Amat[Ayz] * suy) + (Amat[Azz] * suz);
313 >  }
314 >  else{
315 >    
316 >    sprintf( painCave.errMsg,
317 >             "Attempt to updateU for atom %d before coords set.\n",
318 >             index );
319 >    painCave.isFatal = 1;
320 >    simError();
321 >  }
322   }
323  
324   void DirectionalAtom::getJ( double theJ[3] ){
# Line 193 | Line 337 | void DirectionalAtom::getTrq( double theT[3] ){
337  
338   void DirectionalAtom::getTrq( double theT[3] ){
339    
340 <  theT[0] = trq[offsetX];
341 <  theT[1] = trq[offsetY];
342 <  theT[2] = trq[offsetZ];
340 >  if( hasCoords ){
341 >    theT[0] = trq[offsetX];
342 >    theT[1] = trq[offsetY];
343 >    theT[2] = trq[offsetZ];
344 >  }
345 >  else{
346 >    
347 >    sprintf( painCave.errMsg,
348 >             "Attempt to get Trq for atom %d before coords set.\n",
349 >             index );
350 >    painCave.isFatal = 1;
351 >    simError();
352 >  }
353   }
354  
355   void DirectionalAtom::addTrq( double theT[3] ){
356    
357 <  trq[offsetX] += theT[0];
358 <  trq[offsetY] += theT[1];
359 <  trq[offsetZ] += theT[2];
357 >  if( hasCoords ){
358 >    trq[offsetX] += theT[0];
359 >    trq[offsetY] += theT[1];
360 >    trq[offsetZ] += theT[2];
361 >  }
362 >  else{
363 >    
364 >    sprintf( painCave.errMsg,
365 >             "Attempt to add Trq for atom %d before coords set.\n",
366 >             index );
367 >    painCave.isFatal = 1;
368 >    simError();
369 >  }
370   }
371  
372  

Diff Legend

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