ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/DirectionalAtom.cpp
(Generate patch)

Comparing trunk/OOPSE_old/src/mdtools/libmdCode/DirectionalAtom.cpp (file contents):
Revision 270 by mmeineke, Fri Feb 14 17:08:46 2003 UTC vs.
Revision 296 by mmeineke, Thu Mar 6 20:05:39 2003 UTC

# Line 6 | Line 6 | double* Atom::trq; // the torque vector  ( space fixed
6   double* Atom::vel; // the velocity array
7   double* Atom::frc; // the forc array
8   double* Atom::trq; // the torque vector  ( space fixed )
9 + double* Atom::Amat; // the rotation matrix
10 + double* Atom::mu;   // the array of dipole moments
11 + double* Atom::ul;   // the lab frame unit directional vector
12  
13   void DirectionalAtom::setA( double the_A[3][3] ){
14    
15 <  Axx = the_A[0][0]; Axy = the_A[0][1]; Axz = the_A[0][2];
16 <  Ayx = the_A[1][0]; Ayy = the_A[1][1]; Ayz = the_A[1][2];
17 <  Azx = the_A[2][0]; Azy = the_A[2][1]; Azz = the_A[2][2];
15 >  Amat[Axx] = the_A[0][0]; Amat[Axy] = the_A[0][1]; Amat[Axz] = the_A[0][2];
16 >  Amat[Ayx] = the_A[1][0]; Amat[Ayy] = the_A[1][1]; Amat[Ayz] = the_A[1][2];
17 >  Amat[Azx] = the_A[2][0]; Amat[Azy] = the_A[2][1]; Amat[Azz] = the_A[2][2];
18 >  
19 >  this->body2Lab( &ul[offsetX] );
20   }
21  
22   void DirectionalAtom::setI( double the_I[3][3] ){
# Line 31 | Line 36 | void DirectionalAtom::setQ( double the_q[4] ){
36    q3Sqr = the_q[3] * the_q[3];
37    
38  
39 <  Axx = q0Sqr + q1Sqr - q2Sqr - q3Sqr;
40 <  Axy = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] );
41 <  Axz = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] );
39 >  Amat[Axx] = q0Sqr + q1Sqr - q2Sqr - q3Sqr;
40 >  Amat[Axy] = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] );
41 >  Amat[Axz] = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] );
42    
43 <  Ayx = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] );
44 <  Ayy = q0Sqr - q1Sqr + q2Sqr - q3Sqr;
45 <  Ayz = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] );
43 >  Amat[Ayx] = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] );
44 >  Amat[Ayy] = q0Sqr - q1Sqr + q2Sqr - q3Sqr;
45 >  Amat[Ayz] = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] );
46  
47 <  Azx = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] );
48 <  Azy = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] );
49 <  Azz = q0Sqr - q1Sqr -q2Sqr +q3Sqr;
47 >  Amat[Azx] = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] );
48 >  Amat[Azy] = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] );
49 >  Amat[Azz] = q0Sqr - q1Sqr -q2Sqr +q3Sqr;
50 >
51 >  this->body2Lab( &ul[offsetX] );
52   }
53  
54   void DirectionalAtom::getA( double the_A[3][3] ){
55    
56 <  the_A[0][0] = Axx;
57 <  the_A[0][1] = Axy;
58 <  the_A[0][2] = Axz;
56 >  the_A[0][0] = Amat[Axx];
57 >  the_A[0][1] = Amat[Axy];
58 >  the_A[0][2] = Amat[Axz];
59  
60 <  the_A[1][0] = Ayx;
61 <  the_A[1][1] = Ayy;
62 <  the_A[1][2] = Ayz;
60 >  the_A[1][0] = Amat[Ayx];
61 >  the_A[1][1] = Amat[Ayy];
62 >  the_A[1][2] = Amat[Ayz];
63  
64 <  the_A[2][0] = Azx;
65 <  the_A[2][1] = Azy;
66 <  the_A[2][2] = Azz;
64 >  the_A[2][0] = Amat[Azx];
65 >  the_A[2][1] = Amat[Azy];
66 >  the_A[2][2] = Amat[Azz];
67   }
68  
69  
# Line 74 | Line 81 | void DirectionalAtom::getQ( double q[4] ){
81    double t, s;
82    double ad1, ad2, ad3;
83  
84 <  t = Axx + Ayy + Azz + 1.0;
84 >  t = Amat[Axx] + Amat[Ayy] + Amat[Azz] + 1.0;
85    if( t > 0.0 ){
86      
87      s = 0.5 / sqrt( t );
88      q[0] = 0.25 / s;
89 <    q[1] = (Ayz - Azy) * s;
90 <    q[2] = (Azx - Axz) * s;
91 <    q[3] = (Axy - Ayx) * s;
89 >    q[1] = (Amat[Ayz] - Amat[Azy]) * s;
90 >    q[2] = (Amat[Azx] - Amat[Axz]) * s;
91 >    q[3] = (Amat[Axy] - Amat[Ayx]) * s;
92    }
93    else{
94      
95 <    ad1 = fabs( Axx );
96 <    ad2 = fabs( Ayy );
97 <    ad3 = fabs( Azz );
95 >    ad1 = fabs( Amat[Axx] );
96 >    ad2 = fabs( Amat[Ayy] );
97 >    ad3 = fabs( Amat[Azz] );
98  
99      if( ad1 >= ad2 && ad1 >= ad3 ){
100        
101 <      s = 2.0 * sqrt( 1.0 + Axx - Ayy - Azz );
102 <      q[0] = (Ayz + Azy) / s;
101 >      s = 2.0 * sqrt( 1.0 + Amat[Axx] - Amat[Ayy] - Amat[Azz] );
102 >      q[0] = (Amat[Ayz] + Amat[Azy]) / s;
103        q[1] = 0.5 / s;
104 <      q[2] = (Axy + Ayx) / s;
105 <      q[3] = (Axz + Azx) / s;
104 >      q[2] = (Amat[Axy] + Amat[Ayx]) / s;
105 >      q[3] = (Amat[Axz] + Amat[Azx]) / s;
106      }
107      else if( ad2 >= ad1 && ad2 >= ad3 ){
108  
109 <      s = sqrt( 1.0 + Ayy - Axx - Azz ) * 2.0;
110 <      q[0] = (Axz + Azx) / s;
111 <      q[1] = (Axy + Ayx) / s;
109 >      s = sqrt( 1.0 + Amat[Ayy] - Amat[Axx] - Amat[Azz] ) * 2.0;
110 >      q[0] = (Amat[Axz] + Amat[Azx]) / s;
111 >      q[1] = (Amat[Axy] + Amat[Ayx]) / s;
112        q[2] = 0.5 / s;
113 <      q[3] = (Ayz + Azy) / s;
113 >      q[3] = (Amat[Ayz] + Amat[Azy]) / s;
114      }
115      else{
116        
117 <      s = sqrt( 1.0 + Azz - Axx - Ayy ) * 2.0;
118 <      q[0] = (Axy + Ayx) / s;
119 <      q[1] = (Axz + Azx) / s;
120 <      q[2] = (Ayz + Azy) / s;
117 >      s = sqrt( 1.0 + Amat[Azz] - Amat[Axx] - Amat[Ayy] ) * 2.0;
118 >      q[0] = (Amat[Axy] + Amat[Ayx]) / s;
119 >      q[1] = (Amat[Axz] + Amat[Azx]) / s;
120 >      q[2] = (Amat[Ayz] + Amat[Azy]) / s;
121        q[3] = 0.5 / s;
122      }
123    }
# Line 119 | Line 126 | void DirectionalAtom::setEuler( double phi, double the
126  
127   void DirectionalAtom::setEuler( double phi, double theta, double psi ){
128    
129 <  Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
130 <  Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
131 <  Axz = sin(theta) * sin(psi);
129 >  Amat[Axx] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
130 >  Amat[Axy] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
131 >  Amat[Axz] = sin(theta) * sin(psi);
132    
133 <  Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
134 <  Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
135 <  Ayz = sin(theta) * cos(psi);
133 >  Amat[Ayx] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
134 >  Amat[Ayy] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
135 >  Amat[Ayz] = sin(theta) * cos(psi);
136  
137 <  Azx = sin(phi) * sin(theta);
138 <  Azy = -cos(phi) * sin(theta);
139 <  Azz = cos(theta);
137 >  Amat[Azx] = sin(phi) * sin(theta);
138 >  Amat[Azy] = -cos(phi) * sin(theta);
139 >  Amat[Azz] = cos(theta);
140 >
141 >  this->body2Lab( &ul[offsetX] );
142   }
143  
144  
# Line 141 | Line 150 | void DirectionalAtom::lab2Body( double r[3] ){
150    rl[1] = r[1];
151    rl[2] = r[2];
152  
153 <  r[0] = (Axx * rl[0]) + (Axy * rl[1]) + (Axz * rl[2]);
154 <  r[1] = (Ayx * rl[0]) + (Ayy * rl[1]) + (Ayz * rl[2]);
155 <  r[2] = (Azx * rl[0]) + (Azy * rl[1]) + (Azz * rl[2]);
153 >  r[0] = (Amat[Axx] * rl[0]) + (Amat[Axy] * rl[1]) + (Amat[Axz] * rl[2]);
154 >  r[1] = (Amat[Ayx] * rl[0]) + (Amat[Ayy] * rl[1]) + (Amat[Ayz] * rl[2]);
155 >  r[2] = (Amat[Azx] * rl[0]) + (Amat[Azy] * rl[1]) + (Amat[Azz] * rl[2]);
156   }
157  
158   void DirectionalAtom::body2Lab( double r[3] ){
# Line 154 | Line 163 | void DirectionalAtom::body2Lab( double r[3] ){
163    rb[1] = r[1];
164    rb[2] = r[2];
165  
166 <  r[0] = (Axx * rb[0]) + (Ayx * rb[1]) + (Azx * rb[2]);
167 <  r[1] = (Axy * rb[0]) + (Ayy * rb[1]) + (Azy * rb[2]);
168 <  r[2] = (Axz * rb[0]) + (Ayz * rb[1]) + (Azz * rb[2]);
166 >  r[0] = (Amat[Axx] * rb[0]) + (Amat[Ayx] * rb[1]) + (Amat[Azx] * rb[2]);
167 >  r[1] = (Amat[Axy] * rb[0]) + (Amat[Ayy] * rb[1]) + (Amat[Azy] * rb[2]);
168 >  r[2] = (Amat[Axz] * rb[0]) + (Amat[Ayz] * rb[1]) + (Amat[Azz] * rb[2]);
169   }
170  

Diff Legend

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