ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/primitives/DirectionalAtom.cpp
(Generate patch)

Comparing trunk/src/primitives/DirectionalAtom.cpp (file contents):
Revision 2 by gezelter, Fri Sep 24 04:16:43 2004 UTC vs.
Revision 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC

# Line 1 | Line 1
1 < #include <math.h>
1 > /*
2 > * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 > *
4 > * The University of Notre Dame grants you ("Licensee") a
5 > * non-exclusive, royalty free, license to use, modify and
6 > * redistribute this software in source and binary code form, provided
7 > * that the following conditions are met:
8 > *
9 > * 1. Acknowledgement of the program authors must be made in any
10 > *    publication of scientific results based in part on use of the
11 > *    program.  An acceptable form of acknowledgement is citation of
12 > *    the article in which the program was described (Matthew
13 > *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 > *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 > *    Parallel Simulation Engine for Molecular Dynamics,"
16 > *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 > *
18 > * 2. Redistributions of source code must retain the above copyright
19 > *    notice, this list of conditions and the following disclaimer.
20 > *
21 > * 3. Redistributions in binary form must reproduce the above copyright
22 > *    notice, this list of conditions and the following disclaimer in the
23 > *    documentation and/or other materials provided with the
24 > *    distribution.
25 > *
26 > * This software is provided "AS IS," without a warranty of any
27 > * kind. All express or implied conditions, representations and
28 > * warranties, including any implied warranty of merchantability,
29 > * fitness for a particular purpose or non-infringement, are hereby
30 > * excluded.  The University of Notre Dame and its licensors shall not
31 > * be liable for any damages suffered by licensee as a result of
32 > * using, modifying or distributing the software or its
33 > * derivatives. In no event will the University of Notre Dame or its
34 > * licensors be liable for any lost revenue, profit or data, or for
35 > * direct, indirect, special, consequential, incidental or punitive
36 > * damages, however caused and regardless of the theory of liability,
37 > * arising out of the use of or inability to use software, even if the
38 > * University of Notre Dame has been advised of the possibility of
39 > * such damages.
40 > */
41 >
42 > #include "primitives/DirectionalAtom.hpp"
43  
44 < #include "Atom.hpp"
4 < #include "DirectionalAtom.hpp"
5 < #include "simError.h"
6 < #include "MatVec3.h"
44 > namespace oopse {
45  
46 < void DirectionalAtom::zeroForces() {
47 <  if( hasCoords ){
48 <
49 <    Atom::zeroForces();
50 <    
51 <    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 <  }
46 > DirectionalAtom::DirectionalAtom(DirectionalAtomType* dAtomType)
47 >                         : Atom(dAtomType){
48 >    objType_= otDAtom;
49 >    if (dAtomType->isMultipole()) {
50 >        electroBodyFrame_ = dAtomType->getElectroBodyFrame();
51 >    }
52   }
53  
54 < void DirectionalAtom::setCoords(void){
54 > Mat3x3d DirectionalAtom::getI() {
55 >    return static_cast<DirectionalAtomType*>(getAtomType())->getI();
56 > }    
57  
58 <  if( myConfig->isAllocated() ){
59 <
60 <    myConfig->getAtomPointers( index,
61 <                     &pos,
62 <                     &vel,
34 <                     &frc,
35 <                     &trq,
36 <                     &Amat,
37 <                     &mu,  
38 <                     &ul);
39 <  }
40 <  else{
41 <    sprintf( painCave.errMsg,
42 <             "Attempted to set Atom %d  coordinates with an unallocated "
43 <             "SimState object.\n", index );
44 <    painCave.isFatal = 1;
45 <    simError();
46 <  }
47 <
48 <  hasCoords = true;
49 <
58 > void DirectionalAtom::setPrevA(const RotMat3x3d& a) {
59 >    ((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
60 >    if (atomType_->isMultipole()) {
61 >        ((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;
62 >    }
63   }
64  
65 < void DirectionalAtom::setA( double the_A[3][3] ){
65 >      
66 > void DirectionalAtom::setA(const RotMat3x3d& a) {
67 >    ((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
68  
69 <  if( hasCoords ){
70 <    Amat[Axx] = the_A[0][0]; Amat[Axy] = the_A[0][1]; Amat[Axz] = the_A[0][2];
71 <    Amat[Ayx] = the_A[1][0]; Amat[Ayy] = the_A[1][1]; Amat[Ayz] = the_A[1][2];
72 <    Amat[Azx] = the_A[2][0]; Amat[Azy] = the_A[2][1]; Amat[Azz] = the_A[2][2];
69 >    if (atomType_->isMultipole()) {
70 >        ((snapshotMan_->getCurrentSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;
71 >    }
72 > }    
73      
74 <    this->updateU();  
75 <  }
61 <  else{
62 <    
63 <    sprintf( painCave.errMsg,
64 <             "Attempt to set Amat for atom %d before coords set.\n",
65 <             index );
66 <    painCave.isFatal = 1;
67 <    simError();
68 <  }
69 < }
74 > void DirectionalAtom::setA(const RotMat3x3d& a, int snapshotNo) {
75 >    ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a;
76  
77 < void DirectionalAtom::setI( double the_I[3][3] ){  
78 <  
79 <  Ixx = the_I[0][0]; Ixy = the_I[0][1]; Ixz = the_I[0][2];
80 <  Iyx = the_I[1][0]; Iyy = the_I[1][1]; Iyz = the_I[1][2];
75 <  Izx = the_I[2][0]; Izy = the_I[2][1]; Izz = the_I[2][2];
76 < }
77 >    if (atomType_->isMultipole()) {
78 >        ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;    
79 >    }
80 > }    
81  
82 < void DirectionalAtom::setQ( double the_q[4] ){
83 <
80 <  double q0Sqr, q1Sqr, q2Sqr, q3Sqr;
81 <
82 <  if( hasCoords ){
83 <    q0Sqr = the_q[0] * the_q[0];
84 <    q1Sqr = the_q[1] * the_q[1];
85 <    q2Sqr = the_q[2] * the_q[2];
86 <    q3Sqr = the_q[3] * the_q[3];
87 <    
88 <    
89 <    Amat[Axx] = q0Sqr + q1Sqr - q2Sqr - q3Sqr;
90 <    Amat[Axy] = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] );
91 <    Amat[Axz] = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] );
92 <    
93 <    Amat[Ayx] = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] );
94 <    Amat[Ayy] = q0Sqr - q1Sqr + q2Sqr - q3Sqr;
95 <    Amat[Ayz] = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] );
96 <    
97 <    Amat[Azx] = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] );
98 <    Amat[Azy] = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] );
99 <    Amat[Azz] = q0Sqr - q1Sqr -q2Sqr +q3Sqr;
100 <    
101 <    this->updateU();
102 <  }
103 <  else{
104 <    
105 <    sprintf( painCave.errMsg,
106 <             "Attempt to set Q for atom %d before coords set.\n",
107 <             index );
108 <    painCave.isFatal = 1;
109 <    simError();
110 <  }
111 <
82 > void DirectionalAtom::rotateBy(const RotMat3x3d& m) {
83 >    setA(m *getA());
84   }
85  
86 < void DirectionalAtom::getA( double the_A[3][3] ){
87 <  
88 <  if( hasCoords ){
89 <    the_A[0][0] = Amat[Axx];
90 <    the_A[0][1] = Amat[Axy];
91 <    the_A[0][2] = Amat[Axz];
92 <    
93 <    the_A[1][0] = Amat[Ayx];
94 <    the_A[1][1] = Amat[Ayy];
95 <    the_A[1][2] = Amat[Ayz];
124 <    
125 <    the_A[2][0] = Amat[Azx];
126 <    the_A[2][1] = Amat[Azy];
127 <    the_A[2][2] = Amat[Azz];
128 <  }
129 <  else{
130 <    
131 <    sprintf( painCave.errMsg,
132 <             "Attempt to get Amat for atom %d before coords set.\n",
133 <             index );
134 <    painCave.isFatal = 1;
135 <    simError();
136 <  }
86 > std::vector<double> DirectionalAtom::getGrad() {
87 >    std::vector<double> grad(6, 0.0);
88 >    Vector3d force;
89 >    Vector3d torque;
90 >    Vector3d myEuler;
91 >    double phi, theta, psi;
92 >    double cphi, sphi, ctheta, stheta;
93 >    Vector3d ephi;
94 >    Vector3d etheta;
95 >    Vector3d epsi;
96  
97 < }
97 >    force = getFrc();
98 >    torque =getTrq();
99 >    myEuler = getA().toEulerAngles();
100  
101 < void DirectionalAtom::printAmatIndex( void ){
101 >    phi = myEuler[0];
102 >    theta = myEuler[1];
103 >    psi = myEuler[2];
104  
105 <  if( hasCoords ){
106 <    std::cerr << "Atom[" << index << "] index =>\n"
107 <              << "[ " << Axx << ", " << Axy << ", " << Axz << " ]\n"
108 <              << "[ " << Ayx << ", " << Ayy << ", " << Ayz << " ]\n"
146 <              << "[ " << Azx << ", " << Azy << ", " << Azz << " ]\n";
147 <  }
148 <  else{
149 <    
150 <    sprintf( painCave.errMsg,
151 <             "Attempt to print Amat indices for atom %d before coords set.\n",
152 <             index );
153 <    painCave.isFatal = 1;
154 <    simError();
155 <  }
156 < }
105 >    cphi = cos(phi);
106 >    sphi = sin(phi);
107 >    ctheta = cos(theta);
108 >    stheta = sin(theta);
109  
110 +    // get unit vectors along the phi, theta and psi rotation axes
111  
112 < void DirectionalAtom::getU( double the_u[3] ){
113 <  
114 <  the_u[0] = sU[2][0];
162 <  the_u[1] = sU[2][1];
163 <  the_u[2] = sU[2][2];
164 <  
165 <  this->body2Lab( the_u );
166 < }
112 >    ephi[0] = 0.0;
113 >    ephi[1] = 0.0;
114 >    ephi[2] = 1.0;
115  
116 < void DirectionalAtom::getQ( double q[4] ){
117 <  
118 <  double t, s;
171 <  double ad1, ad2, ad3;
116 >    etheta[0] = cphi;
117 >    etheta[1] = sphi;
118 >    etheta[2] = 0.0;
119  
120 <  if( hasCoords ){
121 <    
122 <    t = Amat[Axx] + Amat[Ayy] + Amat[Azz] + 1.0;
176 <    if( t > 0.0 ){
177 <      
178 <      s = 0.5 / sqrt( t );
179 <      q[0] = 0.25 / s;
180 <      q[1] = (Amat[Ayz] - Amat[Azy]) * s;
181 <      q[2] = (Amat[Azx] - Amat[Axz]) * s;
182 <      q[3] = (Amat[Axy] - Amat[Ayx]) * s;
183 <    }
184 <    else{
185 <      
186 <      ad1 = fabs( Amat[Axx] );
187 <      ad2 = fabs( Amat[Ayy] );
188 <      ad3 = fabs( Amat[Azz] );
189 <      
190 <      if( ad1 >= ad2 && ad1 >= ad3 ){
191 <        
192 <        s = 2.0 * sqrt( 1.0 + Amat[Axx] - Amat[Ayy] - Amat[Azz] );
193 <        q[0] = (Amat[Ayz] + Amat[Azy]) / s;
194 <        q[1] = 0.5 / s;
195 <        q[2] = (Amat[Axy] + Amat[Ayx]) / s;
196 <        q[3] = (Amat[Axz] + Amat[Azx]) / s;
197 <      }
198 <      else if( ad2 >= ad1 && ad2 >= ad3 ){
199 <        
200 <        s = sqrt( 1.0 + Amat[Ayy] - Amat[Axx] - Amat[Azz] ) * 2.0;
201 <        q[0] = (Amat[Axz] + Amat[Azx]) / s;
202 <        q[1] = (Amat[Axy] + Amat[Ayx]) / s;
203 <        q[2] = 0.5 / s;
204 <        q[3] = (Amat[Ayz] + Amat[Azy]) / s;
205 <      }
206 <      else{
207 <        
208 <        s = sqrt( 1.0 + Amat[Azz] - Amat[Axx] - Amat[Ayy] ) * 2.0;
209 <        q[0] = (Amat[Axy] + Amat[Ayx]) / s;
210 <        q[1] = (Amat[Axz] + Amat[Azx]) / s;
211 <        q[2] = (Amat[Ayz] + Amat[Azy]) / s;
212 <        q[3] = 0.5 / s;
213 <      }
214 <    }
215 <  }
216 <  else{
217 <    
218 <    sprintf( painCave.errMsg,
219 <             "Attempt to get Q for atom %d before coords set.\n",
220 <             index );
221 <    painCave.isFatal = 1;
222 <    simError();
223 <  }
224 < }
120 >    epsi[0] = stheta * cphi;
121 >    epsi[1] = stheta * sphi;
122 >    epsi[2] = ctheta;
123  
124 < void DirectionalAtom::setUnitFrameFromEuler(double phi,
125 <                                            double theta,
126 <                                            double psi) {
124 >    //gradient is equal to -force
125 >    for (int j = 0 ; j<3; j++)
126 >        grad[j] = -force[j];
127  
128 <  double myA[3][3];
231 <  double uFrame[3][3];
232 <  double len;
233 <  int i, j;
234 <  
235 <  myA[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
236 <  myA[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
237 <  myA[0][2] = sin(theta) * sin(psi);
238 <  
239 <  myA[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
240 <  myA[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
241 <  myA[1][2] = sin(theta) * cos(psi);
242 <  
243 <  myA[2][0] = sin(phi) * sin(theta);
244 <  myA[2][1] = -cos(phi) * sin(theta);
245 <  myA[2][2] = cos(theta);
246 <  
247 <  // Make the unit Frame:
128 >    for (int j = 0; j < 3; j++ ) {
129  
130 <  for (i=0; i < 3; i++)
131 <    for (j=0; j < 3; j++)
132 <      uFrame[i][j] = 0.0;
130 >        grad[3] += torque[j]*ephi[j];
131 >        grad[4] += torque[j]*etheta[j];
132 >        grad[5] += torque[j]*epsi[j];
133  
253 –  for (i=0; i < 3; i++)
254 –    uFrame[i][i] = 1.0;
255 –
256 –  // rotate by the given rotation matrix:
257 –
258 –  matMul3(myA, uFrame, sU);
259 –
260 –  // renormalize column vectors:
261 –
262 –  for (i=0; i < 3; i++) {
263 –    len = 0.0;
264 –    for (j = 0; j < 3; j++) {
265 –      len += sU[i][j]*sU[i][j];
134      }
267 –    len = sqrt(len);
268 –    for (j = 0; j < 3; j++) {
269 –      sU[i][j] /= len;    
270 –    }
271 –  }
272 –  
273 –  // sU now contains the coordinates of the 'special' frame;
274 –    
275 – }
276 –
277 – void DirectionalAtom::setEuler( double phi, double theta, double psi ){
278 –  
279 –  if( hasCoords ){
280 –    Amat[Axx] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
281 –    Amat[Axy] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
282 –    Amat[Axz] = sin(theta) * sin(psi);
135      
136 <    Amat[Ayx] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
137 <    Amat[Ayy] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
286 <    Amat[Ayz] = sin(theta) * cos(psi);
287 <    
288 <    Amat[Azx] = sin(phi) * sin(theta);
289 <    Amat[Azy] = -cos(phi) * sin(theta);
290 <    Amat[Azz] = cos(theta);
291 <    
292 <    this->updateU();
293 <  }
294 <  else{
295 <    
296 <    sprintf( painCave.errMsg,
297 <             "Attempt to set Euler angles for atom %d before coords set.\n",
298 <             index );
299 <    painCave.isFatal = 1;
300 <    simError();
301 <  }
302 < }
136 >    return grad;
137 > }    
138  
139 + void DirectionalAtom::accept(BaseVisitor* v) {
140 +    v->visit(this);
141 + }    
142  
305 – void DirectionalAtom::lab2Body( double r[3] ){
306 –
307 –  double rl[3]; // the lab frame vector
308 –  
309 –  if( hasCoords ){
310 –    rl[0] = r[0];
311 –    rl[1] = r[1];
312 –    rl[2] = r[2];
313 –    
314 –    r[0] = (Amat[Axx] * rl[0]) + (Amat[Axy] * rl[1]) + (Amat[Axz] * rl[2]);
315 –    r[1] = (Amat[Ayx] * rl[0]) + (Amat[Ayy] * rl[1]) + (Amat[Ayz] * rl[2]);
316 –    r[2] = (Amat[Azx] * rl[0]) + (Amat[Azy] * rl[1]) + (Amat[Azz] * rl[2]);
317 –  }
318 –  else{
319 –    
320 –    sprintf( painCave.errMsg,
321 –             "Attempt to convert lab2body for atom %d before coords set.\n",
322 –             index );
323 –    painCave.isFatal = 1;
324 –    simError();
325 –  }
326 –
143   }
144  
329 – void DirectionalAtom::rotateBy( double by_A[3][3]) {
330 –
331 –  // Check this
332 –  
333 –  double r00, r01, r02, r10, r11, r12, r20, r21, r22;
334 –
335 –  if( hasCoords ){
336 –
337 –    r00 = by_A[0][0]*Amat[Axx] + by_A[0][1]*Amat[Ayx] + by_A[0][2]*Amat[Azx];
338 –    r01 = by_A[0][0]*Amat[Axy] + by_A[0][1]*Amat[Ayy] + by_A[0][2]*Amat[Azy];
339 –    r02 = by_A[0][0]*Amat[Axz] + by_A[0][1]*Amat[Ayz] + by_A[0][2]*Amat[Azz];
340 –    
341 –    r10 = by_A[1][0]*Amat[Axx] + by_A[1][1]*Amat[Ayx] + by_A[1][2]*Amat[Azx];
342 –    r11 = by_A[1][0]*Amat[Axy] + by_A[1][1]*Amat[Ayy] + by_A[1][2]*Amat[Azy];
343 –    r12 = by_A[1][0]*Amat[Axz] + by_A[1][1]*Amat[Ayz] + by_A[1][2]*Amat[Azz];
344 –    
345 –    r20 = by_A[2][0]*Amat[Axx] + by_A[2][1]*Amat[Ayx] + by_A[2][2]*Amat[Azx];
346 –    r21 = by_A[2][0]*Amat[Axy] + by_A[2][1]*Amat[Ayy] + by_A[2][2]*Amat[Azy];
347 –    r22 = by_A[2][0]*Amat[Axz] + by_A[2][1]*Amat[Ayz] + by_A[2][2]*Amat[Azz];
348 –    
349 –    Amat[Axx] = r00; Amat[Axy] = r01; Amat[Axz] = r02;
350 –    Amat[Ayx] = r10; Amat[Ayy] = r11; Amat[Ayz] = r12;
351 –    Amat[Azx] = r20; Amat[Azy] = r21; Amat[Azz] = r22;
352 –
353 –  }
354 –  else{
355 –    
356 –    sprintf( painCave.errMsg,
357 –             "Attempt to rotate frame for atom %d before coords set.\n",
358 –             index );
359 –    painCave.isFatal = 1;
360 –    simError();
361 –  }
362 –
363 – }
364 –
365 –
366 – void DirectionalAtom::body2Lab( double r[3] ){
367 –
368 –  double rb[3]; // the body frame vector
369 –  
370 –  if( hasCoords ){
371 –    rb[0] = r[0];
372 –    rb[1] = r[1];
373 –    rb[2] = r[2];
374 –    
375 –    r[0] = (Amat[Axx] * rb[0]) + (Amat[Ayx] * rb[1]) + (Amat[Azx] * rb[2]);
376 –    r[1] = (Amat[Axy] * rb[0]) + (Amat[Ayy] * rb[1]) + (Amat[Azy] * rb[2]);
377 –    r[2] = (Amat[Axz] * rb[0]) + (Amat[Ayz] * rb[1]) + (Amat[Azz] * rb[2]);
378 –  }
379 –  else{
380 –    
381 –    sprintf( painCave.errMsg,
382 –             "Attempt to convert body2lab for atom %d before coords set.\n",
383 –             index );
384 –    painCave.isFatal = 1;
385 –    simError();
386 –  }
387 – }
388 –
389 – void DirectionalAtom::updateU( void ){
390 –
391 –  if( hasCoords ){
392 –    ul[offsetX] = (Amat[Axx] * sU[2][0]) +
393 –      (Amat[Ayx] * sU[2][1]) + (Amat[Azx] * sU[2][2]);
394 –    ul[offsetY] = (Amat[Axy] * sU[2][0]) +
395 –      (Amat[Ayy] * sU[2][1]) + (Amat[Azy] * sU[2][2]);
396 –    ul[offsetZ] = (Amat[Axz] * sU[2][0]) +
397 –      (Amat[Ayz] * sU[2][1]) + (Amat[Azz] * sU[2][2]);
398 –  }
399 –  else{
400 –    
401 –    sprintf( painCave.errMsg,
402 –             "Attempt to updateU for atom %d before coords set.\n",
403 –             index );
404 –    painCave.isFatal = 1;
405 –    simError();
406 –  }
407 – }
408 –
409 – void DirectionalAtom::getJ( double theJ[3] ){
410 –  
411 –  theJ[0] = jx;
412 –  theJ[1] = jy;
413 –  theJ[2] = jz;
414 – }
415 –
416 – void DirectionalAtom::setJ( double theJ[3] ){
417 –  
418 –  jx = theJ[0];
419 –  jy = theJ[1];
420 –  jz = theJ[2];
421 – }
422 –
423 – void DirectionalAtom::getTrq( double theT[3] ){
424 –  
425 –  if( hasCoords ){
426 –    theT[0] = trq[offsetX];
427 –    theT[1] = trq[offsetY];
428 –    theT[2] = trq[offsetZ];
429 –  }
430 –  else{
431 –    
432 –    sprintf( painCave.errMsg,
433 –             "Attempt to get Trq for atom %d before coords set.\n",
434 –             index );
435 –    painCave.isFatal = 1;
436 –    simError();
437 –  }
438 – }
439 –
440 – void DirectionalAtom::addTrq( double theT[3] ){
441 –  
442 –  if( hasCoords ){
443 –    trq[offsetX] += theT[0];
444 –    trq[offsetY] += theT[1];
445 –    trq[offsetZ] += theT[2];
446 –  }
447 –  else{
448 –    
449 –    sprintf( painCave.errMsg,
450 –             "Attempt to add Trq for atom %d before coords set.\n",
451 –             index );
452 –    painCave.isFatal = 1;
453 –    simError();
454 –  }
455 – }
456 –
457 –
458 – void DirectionalAtom::getI( double the_I[3][3] ){
459 –  
460 –  the_I[0][0] = Ixx;
461 –  the_I[0][1] = Ixy;
462 –  the_I[0][2] = Ixz;
463 –
464 –  the_I[1][0] = Iyx;
465 –  the_I[1][1] = Iyy;
466 –  the_I[1][2] = Iyz;
467 –
468 –  the_I[2][0] = Izx;
469 –  the_I[2][1] = Izy;
470 –  the_I[2][2] = Izz;
471 – }
472 –
473 – void DirectionalAtom::getGrad( double grad[6] ) {
474 –
475 –  double myEuler[3];
476 –  double phi, theta, psi;
477 –  double cphi, sphi, ctheta, stheta;
478 –  double ephi[3];
479 –  double etheta[3];
480 –  double epsi[3];
481 –
482 –  this->getEulerAngles(myEuler);
483 –
484 –  phi = myEuler[0];
485 –  theta = myEuler[1];
486 –  psi = myEuler[2];
487 –
488 –  cphi = cos(phi);
489 –  sphi = sin(phi);
490 –  ctheta = cos(theta);
491 –  stheta = sin(theta);
492 –
493 –  // get unit vectors along the phi, theta and psi rotation axes
494 –
495 –  ephi[0] = 0.0;
496 –  ephi[1] = 0.0;
497 –  ephi[2] = 1.0;
498 –
499 –  etheta[0] = cphi;
500 –  etheta[1] = sphi;
501 –  etheta[2] = 0.0;
502 –  
503 –  epsi[0] = stheta * cphi;
504 –  epsi[1] = stheta * sphi;
505 –  epsi[2] = ctheta;
506 –  
507 –  for (int j = 0 ; j<3; j++)
508 –    grad[j] = frc[j];
509 –
510 –  grad[3] = 0;
511 –  grad[4] = 0;
512 –  grad[5] = 0;
513 –
514 –  for (int j = 0; j < 3; j++ ) {
515 –    
516 –    grad[3] += trq[j]*ephi[j];
517 –    grad[4] += trq[j]*etheta[j];
518 –    grad[5] += trq[j]*epsi[j];
519 –    
520 –  }
521 –
522 – }
523 –
524 – /**
525 –  * getEulerAngles computes a set of Euler angle values consistent
526 –  *  with an input rotation matrix.  They are returned in the following
527 –  * order:
528 –  *  myEuler[0] = phi;
529 –  *  myEuler[1] = theta;
530 –  *  myEuler[2] = psi;
531 – */
532 – void DirectionalAtom::getEulerAngles(double myEuler[3]) {
533 –
534 –  // We use so-called "x-convention", which is the most common definition.
535 –  // In this convention, the rotation given by Euler angles (phi, theta, psi), where the first
536 –  // rotation is by an angle phi about the z-axis, the second is by an angle  
537 –  // theta (0 <= theta <= 180)about the x-axis, and thethird is by an angle psi about the
538 –  //z-axis (again).
539 –  
540 –  
541 –  double phi,theta,psi,eps;
542 –  double ctheta,stheta;
543 –
544 –  // set the tolerance for Euler angles and rotation elements
545 –  
546 –  eps = 1.0e-8;
547 –
548 –  theta = acos(min(1.0,max(-1.0,Amat[Azz])));
549 –  ctheta = Amat[Azz];
550 –  stheta = sqrt(1.0 - ctheta * ctheta);
551 –
552 –  // when sin(theta) is close to 0, we need to consider singularity
553 –  // In this case, we can assign an arbitary value to phi (or psi), and then determine
554 –  // the psi (or phi) or vice-versa. We'll assume that phi always gets the rotation, and psi is 0
555 –  // in cases of singularity.  
556 –  // we use atan2 instead of atan, since atan2 will give us -Pi to Pi.
557 –  // Since 0 <= theta <= 180, sin(theta) will be always non-negative. Therefore, it never
558 –  // change the sign of both of the parameters passed to atan2.
559 –  
560 –  if (fabs(stheta) <= eps){
561 –    psi = 0.0;
562 –    phi = atan2(-Amat[Ayx], Amat[Axx]);  
563 –  }
564 –  // we only have one unique solution
565 –  else{    
566 –      phi = atan2(Amat[Azx], -Amat[Azy]);
567 –      psi = atan2(Amat[Axz], Amat[Ayz]);
568 –  }
569 –
570 –  //wrap phi and psi, make sure they are in the range from 0 to 2*Pi
571 –  //if (phi < 0)
572 –  //  phi += M_PI;
573 –
574 –  //if (psi < 0)
575 –  //  psi += M_PI;
576 –
577 –  myEuler[0] = phi;
578 –  myEuler[1] = theta;
579 –  myEuler[2] = psi;
580 –  
581 –  return;
582 – }
583 –
584 – double DirectionalAtom::getZangle( ){
585 –  
586 –  if( hasCoords ){
587 –    return zAngle;
588 –  }
589 –  else{
590 –    
591 –    sprintf( painCave.errMsg,
592 –             "Attempt to get zAngle for atom %d before coords set.\n",
593 –             index );
594 –    painCave.isFatal = 1;
595 –    simError();
596 –    return 0;
597 –  }
598 – }
599 –
600 – void DirectionalAtom::setZangle( double zAng ){
601 –  
602 –  if( hasCoords ){
603 –    zAngle = zAng;
604 –  }
605 –  else{
606 –    
607 –    sprintf( painCave.errMsg,
608 –             "Attempt to set zAngle for atom %d before coords set.\n",
609 –             index );
610 –    painCave.isFatal = 1;
611 –    simError();
612 –  }
613 – }
614 –
615 – void DirectionalAtom::addZangle( double zAng ){
616 –  
617 –  if( hasCoords ){
618 –    zAngle += zAng;
619 –  }
620 –  else{
621 –    
622 –    sprintf( painCave.errMsg,
623 –             "Attempt to add zAngle to atom %d before coords set.\n",
624 –             index );
625 –    painCave.isFatal = 1;
626 –    simError();
627 –  }
628 – }
629 –
630 – double DirectionalAtom::max(double x, double  y) {  
631 –  return (x > y) ? x : y;
632 – }
633 –
634 – double DirectionalAtom::min(double x, double  y) {  
635 –  return (x > y) ? y : x;
636 – }

Diff Legend

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