ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/DirectionalAtom.cpp
Revision: 348
Committed: Fri Mar 14 21:33:10 2003 UTC (21 years, 4 months ago) by mmeineke
File size: 4896 byte(s)
Log Message:
mostly compiles. interface twixt c and fortran is broken. (c needs to be brought up to date with fortran.)

File Contents

# User Rev Content
1 mmeineke 270 #include <cmath>
2    
3     #include "Atom.hpp"
4    
5     double* Atom::pos; // the position array
6     double* Atom::vel; // the velocity array
7     double* Atom::frc; // the forc array
8     double* Atom::trq; // the torque vector ( space fixed )
9 mmeineke 296 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 mmeineke 270
13     void DirectionalAtom::setA( double the_A[3][3] ){
14    
15 mmeineke 296 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 mmeineke 299
19     this->updateU();
20 mmeineke 270 }
21    
22     void DirectionalAtom::setI( double the_I[3][3] ){
23    
24     Ixx = the_I[0][0]; Ixy = the_I[0][1]; Ixz = the_I[0][2];
25     Iyx = the_I[1][0]; Iyy = the_I[1][1]; Iyz = the_I[1][2];
26     Izx = the_I[2][0]; Izy = the_I[2][1]; Izz = the_I[2][2];
27     }
28    
29     void DirectionalAtom::setQ( double the_q[4] ){
30    
31     double q0Sqr, q1Sqr, q2Sqr, q3Sqr;
32    
33     q0Sqr = the_q[0] * the_q[0];
34     q1Sqr = the_q[1] * the_q[1];
35     q2Sqr = the_q[2] * the_q[2];
36     q3Sqr = the_q[3] * the_q[3];
37    
38    
39 mmeineke 296 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 mmeineke 270
43 mmeineke 296 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 mmeineke 270
47 mmeineke 296 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 mmeineke 299 this->updateU();
52 mmeineke 270 }
53    
54     void DirectionalAtom::getA( double the_A[3][3] ){
55    
56 mmeineke 296 the_A[0][0] = Amat[Axx];
57     the_A[0][1] = Amat[Axy];
58     the_A[0][2] = Amat[Axz];
59 mmeineke 270
60 mmeineke 296 the_A[1][0] = Amat[Ayx];
61     the_A[1][1] = Amat[Ayy];
62     the_A[1][2] = Amat[Ayz];
63 mmeineke 270
64 mmeineke 296 the_A[2][0] = Amat[Azx];
65     the_A[2][1] = Amat[Azy];
66     the_A[2][2] = Amat[Azz];
67 mmeineke 270 }
68    
69    
70     void DirectionalAtom::getU( double the_u[3] ){
71    
72     the_u[0] = sux;
73     the_u[1] = suy;
74     the_u[2] = suz;
75    
76     this->body2Lab( the_u );
77     }
78    
79     void DirectionalAtom::getQ( double q[4] ){
80    
81     double t, s;
82     double ad1, ad2, ad3;
83    
84 mmeineke 296 t = Amat[Axx] + Amat[Ayy] + Amat[Azz] + 1.0;
85 mmeineke 270 if( t > 0.0 ){
86    
87     s = 0.5 / sqrt( t );
88     q[0] = 0.25 / s;
89 mmeineke 296 q[1] = (Amat[Ayz] - Amat[Azy]) * s;
90     q[2] = (Amat[Azx] - Amat[Axz]) * s;
91     q[3] = (Amat[Axy] - Amat[Ayx]) * s;
92 mmeineke 270 }
93     else{
94    
95 mmeineke 296 ad1 = fabs( Amat[Axx] );
96     ad2 = fabs( Amat[Ayy] );
97     ad3 = fabs( Amat[Azz] );
98 mmeineke 270
99     if( ad1 >= ad2 && ad1 >= ad3 ){
100    
101 mmeineke 296 s = 2.0 * sqrt( 1.0 + Amat[Axx] - Amat[Ayy] - Amat[Azz] );
102     q[0] = (Amat[Ayz] + Amat[Azy]) / s;
103 mmeineke 270 q[1] = 0.5 / s;
104 mmeineke 296 q[2] = (Amat[Axy] + Amat[Ayx]) / s;
105     q[3] = (Amat[Axz] + Amat[Azx]) / s;
106 mmeineke 270 }
107     else if( ad2 >= ad1 && ad2 >= ad3 ){
108    
109 mmeineke 296 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 mmeineke 270 q[2] = 0.5 / s;
113 mmeineke 296 q[3] = (Amat[Ayz] + Amat[Azy]) / s;
114 mmeineke 270 }
115     else{
116    
117 mmeineke 296 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 mmeineke 270 q[3] = 0.5 / s;
122     }
123     }
124     }
125    
126    
127     void DirectionalAtom::setEuler( double phi, double theta, double psi ){
128    
129 mmeineke 296 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 mmeineke 270
133 mmeineke 296 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 mmeineke 270
137 mmeineke 296 Amat[Azx] = sin(phi) * sin(theta);
138     Amat[Azy] = -cos(phi) * sin(theta);
139     Amat[Azz] = cos(theta);
140    
141 mmeineke 299 this->updateU();
142 mmeineke 270 }
143    
144    
145     void DirectionalAtom::lab2Body( double r[3] ){
146    
147     double rl[3]; // the lab frame vector
148    
149     rl[0] = r[0];
150     rl[1] = r[1];
151     rl[2] = r[2];
152    
153 mmeineke 296 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 mmeineke 270 }
157    
158     void DirectionalAtom::body2Lab( double r[3] ){
159    
160     double rb[3]; // the body frame vector
161    
162     rb[0] = r[0];
163     rb[1] = r[1];
164     rb[2] = r[2];
165    
166 mmeineke 296 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 mmeineke 270 }
170    
171 mmeineke 299 void DirectionalAtom::updateU( void ){
172    
173 mmeineke 348 ul[offsetX] = (Amat[Axx] * sux) + (Amat[Ayx] * suy) + (Amat[Azx] * suz);
174     ul[offsetY] = (Amat[Axy] * sux) + (Amat[Ayy] * suy) + (Amat[Azy] * suz);
175     ul[offsetZ] = (Amat[Axz] * sux) + (Amat[Ayz] * suy) + (Amat[Azz] * suz);
176 mmeineke 299 }
177