ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/primitives/DirectionalAtom.cpp
(Generate patch)

Comparing trunk/OOPSE-4/src/primitives/DirectionalAtom.cpp (file contents):
Revision 1709 by gezelter, Thu Nov 4 16:22:03 2004 UTC vs.
Revision 1957 by tim, Tue Jan 25 17:45:23 2005 UTC

# Line 1 | Line 1
1 < #include <math.h>
2 <
3 < #include "primitives/Atom.hpp"
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   #include "utils/simError.h"
44 < #include "math/MatVec3.h"
44 > namespace oopse {
45  
46 < void DirectionalAtom::zeroForces() {
47 <  if( hasCoords ){
46 > DirectionalAtom::DirectionalAtom(DirectionalAtomType* dAtomType)
47 >                         : Atom(dAtomType){
48 >    objType_= otDAtom;
49 >    if (dAtomType->isMultipole()) {
50 >        electroBodyFrame_ = dAtomType->getElectroBodyFrame();
51 >    }
52  
53 <    Atom::zeroForces();
54 <    
55 <    trq[offsetX] = 0.0;
56 <    trq[offsetY] = 0.0;
57 <    trq[offsetZ] = 0.0;
58 <  }
59 <  else{
60 <    
61 <    sprintf( painCave.errMsg,
62 <             "Attempt to zero frc and trq for atom %d before coords set.\n",
21 <             index );
22 <    painCave.isFatal = 1;
23 <    simError();
24 <  }
25 < }
53 >    //check if one of the diagonal inertia tensor of this directional atom  is zero
54 >    int nLinearAxis = 0;
55 >    Mat3x3d inertiaTensor = getI();
56 >    for (int i = 0; i < 3; i++) {    
57 >        if (fabs(inertiaTensor(i, i)) < oopse::epsilon) {
58 >            linear_ = true;
59 >            linearAxis_ = i;
60 >            ++ nLinearAxis;
61 >        }
62 >    }
63  
64 < void DirectionalAtom::setCoords(void){
65 <
66 <  if( myConfig->isAllocated() ){
67 <
68 <    myConfig->getAtomPointers( index,
69 <                     &pos,
70 <                     &vel,
71 <                     &frc,
72 <                     &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 <
64 >    if (nLinearAxis > 1) {
65 >        sprintf( painCave.errMsg,
66 >            "Directional Atom error.\n"
67 >            "\tOOPSE found more than one axis in this directional atom with a vanishing \n"
68 >            "\tmoment of inertia.");
69 >        painCave.isFatal = 1;
70 >        simError();
71 >    }
72 >      
73   }
74  
75 < void DirectionalAtom::setA( double the_A[3][3] ){
75 > Mat3x3d DirectionalAtom::getI() {
76 >    return static_cast<DirectionalAtomType*>(getAtomType())->getI();
77 > }    
78  
79 <  if( hasCoords ){
80 <    Amat[Axx] = the_A[0][0]; Amat[Axy] = the_A[0][1]; Amat[Axz] = the_A[0][2];
81 <    Amat[Ayx] = the_A[1][0]; Amat[Ayy] = the_A[1][1]; Amat[Ayz] = the_A[1][2];
82 <    Amat[Azx] = the_A[2][0]; Amat[Azy] = the_A[2][1]; Amat[Azz] = the_A[2][2];
58 <    
59 <    this->updateU();  
60 <  }
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 < }
70 <
71 < void DirectionalAtom::setI( double the_I[3][3] ){  
72 <
73 <  int n_linear_coords, i, j;
74 <  
75 <  Ixx = the_I[0][0]; Ixy = the_I[0][1]; Ixz = the_I[0][2];
76 <  Iyx = the_I[1][0]; Iyy = the_I[1][1]; Iyz = the_I[1][2];
77 <  Izx = the_I[2][0]; Izy = the_I[2][1]; Izz = the_I[2][2];
78 <  
79 <  n_linear_coords = 0;
80 <
81 <  for (i = 0; i<3; i++) {
82 <    if (fabs(the_I[i][i]) < momIntTol) {
83 <      is_linear = true;
84 <      n_linear_coords++;
85 <      linear_axis = i;
79 > void DirectionalAtom::setPrevA(const RotMat3x3d& a) {
80 >    ((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
81 >    if (atomType_->isMultipole()) {
82 >        ((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;
83      }
87 –  }
88 –  
89 –  if (n_linear_coords > 1) {
90 –    sprintf( painCave.errMsg,
91 –             "DirectionalAtom error.\n"
92 –             "\tOOPSE was told to set more than one axis in this\n"
93 –             "\tDirectionalAtom to a vanishing moment of inertia.\n"
94 –             "\tThis should not be a DirectionalAtom.  Use an Atom.\n"
95 –             );
96 –      painCave.isFatal = 1;
97 –      simError();
98 –  }
99 –
100 –
84   }
85  
86 < void DirectionalAtom::setQ( double the_q[4] ){
86 >      
87 > void DirectionalAtom::setA(const RotMat3x3d& a) {
88 >    ((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
89  
90 <  double q0Sqr, q1Sqr, q2Sqr, q3Sqr;
91 <
92 <  if( hasCoords ){
93 <    q0Sqr = the_q[0] * the_q[0];
109 <    q1Sqr = the_q[1] * the_q[1];
110 <    q2Sqr = the_q[2] * the_q[2];
111 <    q3Sqr = the_q[3] * the_q[3];
90 >    if (atomType_->isMultipole()) {
91 >        ((snapshotMan_->getCurrentSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;
92 >    }
93 > }    
94      
95 <    
96 <    Amat[Axx] = q0Sqr + q1Sqr - q2Sqr - q3Sqr;
115 <    Amat[Axy] = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] );
116 <    Amat[Axz] = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] );
117 <    
118 <    Amat[Ayx] = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] );
119 <    Amat[Ayy] = q0Sqr - q1Sqr + q2Sqr - q3Sqr;
120 <    Amat[Ayz] = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] );
121 <    
122 <    Amat[Azx] = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] );
123 <    Amat[Azy] = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] );
124 <    Amat[Azz] = q0Sqr - q1Sqr -q2Sqr +q3Sqr;
125 <    
126 <    this->updateU();
127 <  }
128 <  else{
129 <    
130 <    sprintf( painCave.errMsg,
131 <             "Attempt to set Q for atom %d before coords set.\n",
132 <             index );
133 <    painCave.isFatal = 1;
134 <    simError();
135 <  }
95 > void DirectionalAtom::setA(const RotMat3x3d& a, int snapshotNo) {
96 >    ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a;
97  
98 < }
98 >    if (atomType_->isMultipole()) {
99 >        ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;    
100 >    }
101 > }    
102  
103 < void DirectionalAtom::getA( double the_A[3][3] ){
104 <  
141 <  if( hasCoords ){
142 <    the_A[0][0] = Amat[Axx];
143 <    the_A[0][1] = Amat[Axy];
144 <    the_A[0][2] = Amat[Axz];
145 <    
146 <    the_A[1][0] = Amat[Ayx];
147 <    the_A[1][1] = Amat[Ayy];
148 <    the_A[1][2] = Amat[Ayz];
149 <    
150 <    the_A[2][0] = Amat[Azx];
151 <    the_A[2][1] = Amat[Azy];
152 <    the_A[2][2] = Amat[Azz];
153 <  }
154 <  else{
155 <    
156 <    sprintf( painCave.errMsg,
157 <             "Attempt to get Amat for atom %d before coords set.\n",
158 <             index );
159 <    painCave.isFatal = 1;
160 <    simError();
161 <  }
162 <
103 > void DirectionalAtom::rotateBy(const RotMat3x3d& m) {
104 >    setA(m *getA());
105   }
106  
107 < void DirectionalAtom::printAmatIndex( void ){
107 > std::vector<double> DirectionalAtom::getGrad() {
108 >    std::vector<double> grad(6, 0.0);
109 >    Vector3d force;
110 >    Vector3d torque;
111 >    Vector3d myEuler;
112 >    double phi, theta, psi;
113 >    double cphi, sphi, ctheta, stheta;
114 >    Vector3d ephi;
115 >    Vector3d etheta;
116 >    Vector3d epsi;
117  
118 <  if( hasCoords ){
119 <    std::cerr << "Atom[" << index << "] index =>\n"
120 <              << "[ " << Axx << ", " << Axy << ", " << Axz << " ]\n"
170 <              << "[ " << Ayx << ", " << Ayy << ", " << Ayz << " ]\n"
171 <              << "[ " << Azx << ", " << Azy << ", " << Azz << " ]\n";
172 <  }
173 <  else{
174 <    
175 <    sprintf( painCave.errMsg,
176 <             "Attempt to print Amat indices for atom %d before coords set.\n",
177 <             index );
178 <    painCave.isFatal = 1;
179 <    simError();
180 <  }
181 < }
118 >    force = getFrc();
119 >    torque =getTrq();
120 >    myEuler = getA().toEulerAngles();
121  
122 +    phi = myEuler[0];
123 +    theta = myEuler[1];
124 +    psi = myEuler[2];
125  
126 < void DirectionalAtom::getU( double the_u[3] ){
127 <  
128 <  the_u[0] = sU[2][0];
129 <  the_u[1] = sU[2][1];
188 <  the_u[2] = sU[2][2];
189 <  
190 <  this->body2Lab( the_u );
191 < }
126 >    cphi = cos(phi);
127 >    sphi = sin(phi);
128 >    ctheta = cos(theta);
129 >    stheta = sin(theta);
130  
131 < void DirectionalAtom::getQ( double q[4] ){
194 <  
195 <  double t, s;
196 <  double ad1, ad2, ad3;
131 >    // get unit vectors along the phi, theta and psi rotation axes
132  
133 <  if( hasCoords ){
134 <    
135 <    t = Amat[Axx] + Amat[Ayy] + Amat[Azz] + 1.0;
201 <    if( t > 0.0 ){
202 <      
203 <      s = 0.5 / sqrt( t );
204 <      q[0] = 0.25 / s;
205 <      q[1] = (Amat[Ayz] - Amat[Azy]) * s;
206 <      q[2] = (Amat[Azx] - Amat[Axz]) * s;
207 <      q[3] = (Amat[Axy] - Amat[Ayx]) * s;
208 <    }
209 <    else{
210 <      
211 <      ad1 = fabs( Amat[Axx] );
212 <      ad2 = fabs( Amat[Ayy] );
213 <      ad3 = fabs( Amat[Azz] );
214 <      
215 <      if( ad1 >= ad2 && ad1 >= ad3 ){
216 <        
217 <        s = 2.0 * sqrt( 1.0 + Amat[Axx] - Amat[Ayy] - Amat[Azz] );
218 <        q[0] = (Amat[Ayz] + Amat[Azy]) / s;
219 <        q[1] = 0.5 / s;
220 <        q[2] = (Amat[Axy] + Amat[Ayx]) / s;
221 <        q[3] = (Amat[Axz] + Amat[Azx]) / s;
222 <      }
223 <      else if( ad2 >= ad1 && ad2 >= ad3 ){
224 <        
225 <        s = sqrt( 1.0 + Amat[Ayy] - Amat[Axx] - Amat[Azz] ) * 2.0;
226 <        q[0] = (Amat[Axz] + Amat[Azx]) / s;
227 <        q[1] = (Amat[Axy] + Amat[Ayx]) / s;
228 <        q[2] = 0.5 / s;
229 <        q[3] = (Amat[Ayz] + Amat[Azy]) / s;
230 <      }
231 <      else{
232 <        
233 <        s = sqrt( 1.0 + Amat[Azz] - Amat[Axx] - Amat[Ayy] ) * 2.0;
234 <        q[0] = (Amat[Axy] + Amat[Ayx]) / s;
235 <        q[1] = (Amat[Axz] + Amat[Azx]) / s;
236 <        q[2] = (Amat[Ayz] + Amat[Azy]) / s;
237 <        q[3] = 0.5 / s;
238 <      }
239 <    }
240 <  }
241 <  else{
242 <    
243 <    sprintf( painCave.errMsg,
244 <             "Attempt to get Q for atom %d before coords set.\n",
245 <             index );
246 <    painCave.isFatal = 1;
247 <    simError();
248 <  }
249 < }
133 >    ephi[0] = 0.0;
134 >    ephi[1] = 0.0;
135 >    ephi[2] = 1.0;
136  
137 < void DirectionalAtom::setUnitFrameFromEuler(double phi,
138 <                                            double theta,
139 <                                            double psi) {
137 >    etheta[0] = cphi;
138 >    etheta[1] = sphi;
139 >    etheta[2] = 0.0;
140  
141 <  double myA[3][3];
142 <  double uFrame[3][3];
143 <  double len;
258 <  int i, j;
259 <  
260 <  myA[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
261 <  myA[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
262 <  myA[0][2] = sin(theta) * sin(psi);
263 <  
264 <  myA[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
265 <  myA[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
266 <  myA[1][2] = sin(theta) * cos(psi);
267 <  
268 <  myA[2][0] = sin(phi) * sin(theta);
269 <  myA[2][1] = -cos(phi) * sin(theta);
270 <  myA[2][2] = cos(theta);
271 <  
272 <  // Make the unit Frame:
141 >    epsi[0] = stheta * cphi;
142 >    epsi[1] = stheta * sphi;
143 >    epsi[2] = ctheta;
144  
145 <  for (i=0; i < 3; i++)
146 <    for (j=0; j < 3; j++)
147 <      uFrame[i][j] = 0.0;
145 >    //gradient is equal to -force
146 >    for (int j = 0 ; j<3; j++)
147 >        grad[j] = -force[j];
148  
149 <  for (i=0; i < 3; i++)
279 <    uFrame[i][i] = 1.0;
149 >    for (int j = 0; j < 3; j++ ) {
150  
151 <  // rotate by the given rotation matrix:
152 <
153 <  matMul3(myA, uFrame, sU);
151 >        grad[3] += torque[j]*ephi[j];
152 >        grad[4] += torque[j]*etheta[j];
153 >        grad[5] += torque[j]*epsi[j];
154  
285 –  // renormalize column vectors:
286 –
287 –  for (i=0; i < 3; i++) {
288 –    len = 0.0;
289 –    for (j = 0; j < 3; j++) {
290 –      len += sU[i][j]*sU[i][j];
155      }
292 –    len = sqrt(len);
293 –    for (j = 0; j < 3; j++) {
294 –      sU[i][j] /= len;    
295 –    }
296 –  }
297 –  
298 –  // sU now contains the coordinates of the 'special' frame;
156      
157 < }
157 >    return grad;
158 > }    
159  
160 < void DirectionalAtom::setEuler( double phi, double theta, double psi ){
161 <  
162 <  if( hasCoords ){
305 <    Amat[Axx] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
306 <    Amat[Axy] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
307 <    Amat[Axz] = sin(theta) * sin(psi);
308 <    
309 <    Amat[Ayx] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
310 <    Amat[Ayy] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
311 <    Amat[Ayz] = sin(theta) * cos(psi);
312 <    
313 <    Amat[Azx] = sin(phi) * sin(theta);
314 <    Amat[Azy] = -cos(phi) * sin(theta);
315 <    Amat[Azz] = cos(theta);
316 <    
317 <    this->updateU();
318 <  }
319 <  else{
320 <    
321 <    sprintf( painCave.errMsg,
322 <             "Attempt to set Euler angles for atom %d before coords set.\n",
323 <             index );
324 <    painCave.isFatal = 1;
325 <    simError();
326 <  }
327 < }
160 > void DirectionalAtom::accept(BaseVisitor* v) {
161 >    v->visit(this);
162 > }    
163  
329 –
330 – void DirectionalAtom::lab2Body( double r[3] ){
331 –
332 –  double rl[3]; // the lab frame vector
333 –  
334 –  if( hasCoords ){
335 –    rl[0] = r[0];
336 –    rl[1] = r[1];
337 –    rl[2] = r[2];
338 –    
339 –    r[0] = (Amat[Axx] * rl[0]) + (Amat[Axy] * rl[1]) + (Amat[Axz] * rl[2]);
340 –    r[1] = (Amat[Ayx] * rl[0]) + (Amat[Ayy] * rl[1]) + (Amat[Ayz] * rl[2]);
341 –    r[2] = (Amat[Azx] * rl[0]) + (Amat[Azy] * rl[1]) + (Amat[Azz] * rl[2]);
342 –  }
343 –  else{
344 –    
345 –    sprintf( painCave.errMsg,
346 –             "Attempt to convert lab2body for atom %d before coords set.\n",
347 –             index );
348 –    painCave.isFatal = 1;
349 –    simError();
350 –  }
351 –
164   }
165  
354 – void DirectionalAtom::rotateBy( double by_A[3][3]) {
355 –
356 –  // Check this
357 –  
358 –  double r00, r01, r02, r10, r11, r12, r20, r21, r22;
359 –
360 –  if( hasCoords ){
361 –
362 –    r00 = by_A[0][0]*Amat[Axx] + by_A[0][1]*Amat[Ayx] + by_A[0][2]*Amat[Azx];
363 –    r01 = by_A[0][0]*Amat[Axy] + by_A[0][1]*Amat[Ayy] + by_A[0][2]*Amat[Azy];
364 –    r02 = by_A[0][0]*Amat[Axz] + by_A[0][1]*Amat[Ayz] + by_A[0][2]*Amat[Azz];
365 –    
366 –    r10 = by_A[1][0]*Amat[Axx] + by_A[1][1]*Amat[Ayx] + by_A[1][2]*Amat[Azx];
367 –    r11 = by_A[1][0]*Amat[Axy] + by_A[1][1]*Amat[Ayy] + by_A[1][2]*Amat[Azy];
368 –    r12 = by_A[1][0]*Amat[Axz] + by_A[1][1]*Amat[Ayz] + by_A[1][2]*Amat[Azz];
369 –    
370 –    r20 = by_A[2][0]*Amat[Axx] + by_A[2][1]*Amat[Ayx] + by_A[2][2]*Amat[Azx];
371 –    r21 = by_A[2][0]*Amat[Axy] + by_A[2][1]*Amat[Ayy] + by_A[2][2]*Amat[Azy];
372 –    r22 = by_A[2][0]*Amat[Axz] + by_A[2][1]*Amat[Ayz] + by_A[2][2]*Amat[Azz];
373 –    
374 –    Amat[Axx] = r00; Amat[Axy] = r01; Amat[Axz] = r02;
375 –    Amat[Ayx] = r10; Amat[Ayy] = r11; Amat[Ayz] = r12;
376 –    Amat[Azx] = r20; Amat[Azy] = r21; Amat[Azz] = r22;
377 –
378 –  }
379 –  else{
380 –    
381 –    sprintf( painCave.errMsg,
382 –             "Attempt to rotate frame for atom %d before coords set.\n",
383 –             index );
384 –    painCave.isFatal = 1;
385 –    simError();
386 –  }
387 –
388 – }
389 –
390 –
391 – void DirectionalAtom::body2Lab( double r[3] ){
392 –
393 –  double rb[3]; // the body frame vector
394 –  
395 –  if( hasCoords ){
396 –    rb[0] = r[0];
397 –    rb[1] = r[1];
398 –    rb[2] = r[2];
399 –    
400 –    r[0] = (Amat[Axx] * rb[0]) + (Amat[Ayx] * rb[1]) + (Amat[Azx] * rb[2]);
401 –    r[1] = (Amat[Axy] * rb[0]) + (Amat[Ayy] * rb[1]) + (Amat[Azy] * rb[2]);
402 –    r[2] = (Amat[Axz] * rb[0]) + (Amat[Ayz] * rb[1]) + (Amat[Azz] * rb[2]);
403 –  }
404 –  else{
405 –    
406 –    sprintf( painCave.errMsg,
407 –             "Attempt to convert body2lab for atom %d before coords set.\n",
408 –             index );
409 –    painCave.isFatal = 1;
410 –    simError();
411 –  }
412 – }
413 –
414 – void DirectionalAtom::updateU( void ){
415 –
416 –  if( hasCoords ){
417 –    ul[offsetX] = (Amat[Axx] * sU[2][0]) +
418 –      (Amat[Ayx] * sU[2][1]) + (Amat[Azx] * sU[2][2]);
419 –    ul[offsetY] = (Amat[Axy] * sU[2][0]) +
420 –      (Amat[Ayy] * sU[2][1]) + (Amat[Azy] * sU[2][2]);
421 –    ul[offsetZ] = (Amat[Axz] * sU[2][0]) +
422 –      (Amat[Ayz] * sU[2][1]) + (Amat[Azz] * sU[2][2]);
423 –  }
424 –  else{
425 –    
426 –    sprintf( painCave.errMsg,
427 –             "Attempt to updateU for atom %d before coords set.\n",
428 –             index );
429 –    painCave.isFatal = 1;
430 –    simError();
431 –  }
432 – }
433 –
434 – void DirectionalAtom::getJ( double theJ[3] ){
435 –  
436 –  theJ[0] = jx;
437 –  theJ[1] = jy;
438 –  theJ[2] = jz;
439 – }
440 –
441 – void DirectionalAtom::setJ( double theJ[3] ){
442 –  
443 –  jx = theJ[0];
444 –  jy = theJ[1];
445 –  jz = theJ[2];
446 – }
447 –
448 – void DirectionalAtom::getTrq( double theT[3] ){
449 –  
450 –  if( hasCoords ){
451 –    theT[0] = trq[offsetX];
452 –    theT[1] = trq[offsetY];
453 –    theT[2] = trq[offsetZ];
454 –  }
455 –  else{
456 –    
457 –    sprintf( painCave.errMsg,
458 –             "Attempt to get Trq for atom %d before coords set.\n",
459 –             index );
460 –    painCave.isFatal = 1;
461 –    simError();
462 –  }
463 – }
464 –
465 – void DirectionalAtom::addTrq( double theT[3] ){
466 –  
467 –  if( hasCoords ){
468 –    trq[offsetX] += theT[0];
469 –    trq[offsetY] += theT[1];
470 –    trq[offsetZ] += theT[2];
471 –  }
472 –  else{
473 –    
474 –    sprintf( painCave.errMsg,
475 –             "Attempt to add Trq for atom %d before coords set.\n",
476 –             index );
477 –    painCave.isFatal = 1;
478 –    simError();
479 –  }
480 – }
481 –
482 –
483 – void DirectionalAtom::getI( double the_I[3][3] ){
484 –  
485 –  the_I[0][0] = Ixx;
486 –  the_I[0][1] = Ixy;
487 –  the_I[0][2] = Ixz;
488 –
489 –  the_I[1][0] = Iyx;
490 –  the_I[1][1] = Iyy;
491 –  the_I[1][2] = Iyz;
492 –
493 –  the_I[2][0] = Izx;
494 –  the_I[2][1] = Izy;
495 –  the_I[2][2] = Izz;
496 – }
497 –
498 – void DirectionalAtom::getGrad( double grad[6] ) {
499 –
500 –  double myEuler[3];
501 –  double phi, theta, psi;
502 –  double cphi, sphi, ctheta, stheta;
503 –  double ephi[3];
504 –  double etheta[3];
505 –  double epsi[3];
506 –
507 –  this->getEulerAngles(myEuler);
508 –
509 –  phi = myEuler[0];
510 –  theta = myEuler[1];
511 –  psi = myEuler[2];
512 –
513 –  cphi = cos(phi);
514 –  sphi = sin(phi);
515 –  ctheta = cos(theta);
516 –  stheta = sin(theta);
517 –
518 –  // get unit vectors along the phi, theta and psi rotation axes
519 –
520 –  ephi[0] = 0.0;
521 –  ephi[1] = 0.0;
522 –  ephi[2] = 1.0;
523 –
524 –  etheta[0] = cphi;
525 –  etheta[1] = sphi;
526 –  etheta[2] = 0.0;
527 –  
528 –  epsi[0] = stheta * cphi;
529 –  epsi[1] = stheta * sphi;
530 –  epsi[2] = ctheta;
531 –  
532 –  for (int j = 0 ; j<3; j++)
533 –    grad[j] = frc[j];
534 –
535 –  grad[3] = 0;
536 –  grad[4] = 0;
537 –  grad[5] = 0;
538 –
539 –  for (int j = 0; j < 3; j++ ) {
540 –    
541 –    grad[3] += trq[j]*ephi[j];
542 –    grad[4] += trq[j]*etheta[j];
543 –    grad[5] += trq[j]*epsi[j];
544 –    
545 –  }
546 –
547 – }
548 –
549 – /**
550 –  * getEulerAngles computes a set of Euler angle values consistent
551 –  *  with an input rotation matrix.  They are returned in the following
552 –  * order:
553 –  *  myEuler[0] = phi;
554 –  *  myEuler[1] = theta;
555 –  *  myEuler[2] = psi;
556 – */
557 – void DirectionalAtom::getEulerAngles(double myEuler[3]) {
558 –
559 –  // We use so-called "x-convention", which is the most common definition.
560 –  // In this convention, the rotation given by Euler angles (phi, theta, psi), where the first
561 –  // rotation is by an angle phi about the z-axis, the second is by an angle  
562 –  // theta (0 <= theta <= 180)about the x-axis, and thethird is by an angle psi about the
563 –  //z-axis (again).
564 –  
565 –  
566 –  double phi,theta,psi,eps;
567 –  double ctheta,stheta;
568 –
569 –  // set the tolerance for Euler angles and rotation elements
570 –  
571 –  eps = 1.0e-8;
572 –
573 –  theta = acos(min(1.0,max(-1.0,Amat[Azz])));
574 –  ctheta = Amat[Azz];
575 –  stheta = sqrt(1.0 - ctheta * ctheta);
576 –
577 –  // when sin(theta) is close to 0, we need to consider singularity
578 –  // In this case, we can assign an arbitary value to phi (or psi), and then determine
579 –  // the psi (or phi) or vice-versa. We'll assume that phi always gets the rotation, and psi is 0
580 –  // in cases of singularity.  
581 –  // we use atan2 instead of atan, since atan2 will give us -Pi to Pi.
582 –  // Since 0 <= theta <= 180, sin(theta) will be always non-negative. Therefore, it never
583 –  // change the sign of both of the parameters passed to atan2.
584 –  
585 –  if (fabs(stheta) <= eps){
586 –    psi = 0.0;
587 –    phi = atan2(-Amat[Ayx], Amat[Axx]);  
588 –  }
589 –  // we only have one unique solution
590 –  else{    
591 –      phi = atan2(Amat[Azx], -Amat[Azy]);
592 –      psi = atan2(Amat[Axz], Amat[Ayz]);
593 –  }
594 –
595 –  //wrap phi and psi, make sure they are in the range from 0 to 2*Pi
596 –  //if (phi < 0)
597 –  //  phi += M_PI;
598 –
599 –  //if (psi < 0)
600 –  //  psi += M_PI;
601 –
602 –  myEuler[0] = phi;
603 –  myEuler[1] = theta;
604 –  myEuler[2] = psi;
605 –  
606 –  return;
607 – }
608 –
609 – double DirectionalAtom::getZangle( ){
610 –  
611 –  if( hasCoords ){
612 –    return zAngle;
613 –  }
614 –  else{
615 –    
616 –    sprintf( painCave.errMsg,
617 –             "Attempt to get zAngle for atom %d before coords set.\n",
618 –             index );
619 –    painCave.isFatal = 1;
620 –    simError();
621 –    return 0;
622 –  }
623 – }
624 –
625 – void DirectionalAtom::setZangle( double zAng ){
626 –  
627 –  if( hasCoords ){
628 –    zAngle = zAng;
629 –  }
630 –  else{
631 –    
632 –    sprintf( painCave.errMsg,
633 –             "Attempt to set zAngle for atom %d before coords set.\n",
634 –             index );
635 –    painCave.isFatal = 1;
636 –    simError();
637 –  }
638 – }
639 –
640 – void DirectionalAtom::addZangle( double zAng ){
641 –  
642 –  if( hasCoords ){
643 –    zAngle += zAng;
644 –  }
645 –  else{
646 –    
647 –    sprintf( painCave.errMsg,
648 –             "Attempt to add zAngle to atom %d before coords set.\n",
649 –             index );
650 –    painCave.isFatal = 1;
651 –    simError();
652 –  }
653 – }
654 –
655 – double DirectionalAtom::max(double x, double  y) {  
656 –  return (x > y) ? x : y;
657 – }
658 –
659 – double DirectionalAtom::min(double x, double  y) {  
660 –  return (x > y) ? y : x;
661 – }

Diff Legend

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