ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DirectionalAtom.cpp
Revision: 689
Committed: Tue Aug 12 19:56:49 2003 UTC (20 years, 11 months ago) by tim
File size: 8936 byte(s)
Log Message:
debugging globals

File Contents

# User Rev Content
1 mmeineke 377 #include <cmath>
2    
3     #include "Atom.hpp"
4 mmeineke 670 #include "simError.h"
5 mmeineke 377
6 mmeineke 670 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 mmeineke 377
26 tim 689 void DirectionalAtom::setCoords(void){
27 mmeineke 414
28 tim 689 if( myConfig->isAllocated() ){
29    
30     myConfig->getAtomPointers( index,
31     &pos,
32     &vel,
33     &frc,
34     &trq,
35     &Amat,
36     &mu,
37     &ul );
38 mmeineke 670 }
39     else{
40     sprintf( painCave.errMsg,
41 tim 689 "Attempted to set Atom %d coordinates with an unallocated "
42     "SimState object.\n" );
43 mmeineke 670 painCave.isFatal = 1;
44     simError();
45     }
46 tim 689
47     hasCoords = true;
48    
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 mmeineke 670 return 0;
62     }
63    
64     void DirectionalAtom::setMu( double the_mu ) {
65    
66     if( hasCoords ){
67     mu[index] = the_mu;
68 tim 689 myMu = the_mu;
69 mmeineke 670 }
70     else{
71 tim 689 myMu = the_mu;
72 mmeineke 670 }
73     }
74    
75 mmeineke 377 void DirectionalAtom::setA( double the_A[3][3] ){
76    
77 mmeineke 670 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 mmeineke 377 }
93    
94     void DirectionalAtom::setI( double the_I[3][3] ){
95    
96     Ixx = the_I[0][0]; Ixy = the_I[0][1]; Ixz = the_I[0][2];
97     Iyx = the_I[1][0]; Iyy = the_I[1][1]; Iyz = the_I[1][2];
98     Izx = the_I[2][0]; Izy = the_I[2][1]; Izz = the_I[2][2];
99     }
100    
101     void DirectionalAtom::setQ( double the_q[4] ){
102    
103     double q0Sqr, q1Sqr, q2Sqr, q3Sqr;
104    
105 mmeineke 670 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 mmeineke 377
135     }
136    
137     void DirectionalAtom::getA( double the_A[3][3] ){
138    
139 mmeineke 670 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 mmeineke 377
161     }
162    
163 mmeineke 597 void DirectionalAtom::printAmatIndex( void ){
164 mmeineke 377
165 mmeineke 670 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 mmeineke 597 }
180    
181    
182 mmeineke 377 void DirectionalAtom::getU( double the_u[3] ){
183    
184     the_u[0] = sux;
185     the_u[1] = suy;
186     the_u[2] = suz;
187    
188     this->body2Lab( the_u );
189     }
190    
191     void DirectionalAtom::getQ( double q[4] ){
192    
193     double t, s;
194     double ad1, ad2, ad3;
195    
196 mmeineke 670 if( hasCoords ){
197 mmeineke 377
198 mmeineke 670 t = Amat[Axx] + Amat[Ayy] + Amat[Azz] + 1.0;
199     if( t > 0.0 ){
200 mmeineke 377
201 mmeineke 670 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 mmeineke 377 }
207     else{
208    
209 mmeineke 670 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 mmeineke 377 }
238     }
239 mmeineke 670 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 mmeineke 377 }
248    
249    
250     void DirectionalAtom::setEuler( double phi, double theta, double psi ){
251    
252 mmeineke 670 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 mmeineke 377 }
276    
277    
278     void DirectionalAtom::lab2Body( double r[3] ){
279    
280     double rl[3]; // the lab frame vector
281    
282 mmeineke 670 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 mmeineke 377
300     }
301    
302     void DirectionalAtom::body2Lab( double r[3] ){
303    
304     double rb[3]; // the body frame vector
305    
306 mmeineke 670 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 mmeineke 377 }
324    
325     void DirectionalAtom::updateU( void ){
326    
327 mmeineke 670 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 mmeineke 377 }
341    
342 mmeineke 599 void DirectionalAtom::getJ( double theJ[3] ){
343    
344     theJ[0] = jx;
345     theJ[1] = jy;
346     theJ[2] = jz;
347     }
348    
349     void DirectionalAtom::setJ( double theJ[3] ){
350    
351     jx = theJ[0];
352     jy = theJ[1];
353     jz = theJ[2];
354     }
355    
356     void DirectionalAtom::getTrq( double theT[3] ){
357    
358 mmeineke 670 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 mmeineke 599 }
372    
373     void DirectionalAtom::addTrq( double theT[3] ){
374    
375 mmeineke 670 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 mmeineke 599 }
389    
390    
391     void DirectionalAtom::getI( double the_I[3][3] ){
392    
393     the_I[0][0] = Ixx;
394     the_I[0][1] = Ixy;
395     the_I[0][2] = Ixz;
396    
397     the_I[1][0] = Iyx;
398     the_I[1][1] = Iyy;
399     the_I[1][2] = Iyz;
400    
401     the_I[2][0] = Izx;
402     the_I[2][1] = Izy;
403     the_I[2][2] = Izz;
404     }