# | 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 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |