OOPSE/libmdtools/DirectionalAtom.cpp

#include <cmath>

#include "Atom.hpp"


void DirectionalAtom::setA( double the_A[3][3] ){
  
  Amat[Axx] = the_A[0][0]; Amat[Axy] = the_A[0][1]; Amat[Axz] = the_A[0][2];
  Amat[Ayx] = the_A[1][0]; Amat[Ayy] = the_A[1][1]; Amat[Ayz] = the_A[1][2];
  Amat[Azx] = the_A[2][0]; Amat[Azy] = the_A[2][1]; Amat[Azz] = the_A[2][2];

  this->updateU();  
}

void DirectionalAtom::setI( double the_I[3][3] ){
  
  Ixx = the_I[0][0]; Ixy = the_I[0][1]; Ixz = the_I[0][2];
  Iyx = the_I[1][0]; Iyy = the_I[1][1]; Iyz = the_I[1][2];
  Izx = the_I[2][0]; Izy = the_I[2][1]; Izz = the_I[2][2];
}

void DirectionalAtom::setQ( double the_q[4] ){

  double q0Sqr, q1Sqr, q2Sqr, q3Sqr;

  q0Sqr = the_q[0] * the_q[0];
  q1Sqr = the_q[1] * the_q[1];
  q2Sqr = the_q[2] * the_q[2];
  q3Sqr = the_q[3] * the_q[3];
  

  Amat[Axx] = q0Sqr + q1Sqr - q2Sqr - q3Sqr;
  Amat[Axy] = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] );
  Amat[Axz] = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] );
  
  Amat[Ayx] = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] );
  Amat[Ayy] = q0Sqr - q1Sqr + q2Sqr - q3Sqr;
  Amat[Ayz] = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] );

  Amat[Azx] = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] );
  Amat[Azy] = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] );
  Amat[Azz] = q0Sqr - q1Sqr -q2Sqr +q3Sqr;

  this->updateU();
}

void DirectionalAtom::getA( double the_A[3][3] ){
  
  the_A[0][0] = Amat[Axx];
  the_A[0][1] = Amat[Axy];
  the_A[0][2] = Amat[Axz];

  the_A[1][0] = Amat[Ayx];
  the_A[1][1] = Amat[Ayy];
  the_A[1][2] = Amat[Ayz];

  the_A[2][0] = Amat[Azx];
  the_A[2][1] = Amat[Azy];
  the_A[2][2] = Amat[Azz];
}


void DirectionalAtom::getU( double the_u[3] ){
  
  the_u[0] = sux;
  the_u[1] = suy;
  the_u[2] = suz;

  this->body2Lab( the_u );
}

void DirectionalAtom::getQ( double q[4] ){
  
  double t, s;
  double ad1, ad2, ad3;

  t = Amat[Axx] + Amat[Ayy] + Amat[Azz] + 1.0;
  if( t > 0.0 ){
    
    s = 0.5 / sqrt( t );
    q[0] = 0.25 / s;
    q[1] = (Amat[Ayz] - Amat[Azy]) * s;
    q[2] = (Amat[Azx] - Amat[Axz]) * s;
    q[3] = (Amat[Axy] - Amat[Ayx]) * s;
  }
  else{
    
    ad1 = fabs( Amat[Axx] );
    ad2 = fabs( Amat[Ayy] );
    ad3 = fabs( Amat[Azz] );

    if( ad1 >= ad2 && ad1 >= ad3 ){
      
      s = 2.0 * sqrt( 1.0 + Amat[Axx] - Amat[Ayy] - Amat[Azz] );
      q[0] = (Amat[Ayz] + Amat[Azy]) / s;
      q[1] = 0.5 / s;
      q[2] = (Amat[Axy] + Amat[Ayx]) / s;
      q[3] = (Amat[Axz] + Amat[Azx]) / s;
    }
    else if( ad2 >= ad1 && ad2 >= ad3 ){

      s = sqrt( 1.0 + Amat[Ayy] - Amat[Axx] - Amat[Azz] ) * 2.0;
      q[0] = (Amat[Axz] + Amat[Azx]) / s;
      q[1] = (Amat[Axy] + Amat[Ayx]) / s;
      q[2] = 0.5 / s;
      q[3] = (Amat[Ayz] + Amat[Azy]) / s;
    }
    else{
      
      s = sqrt( 1.0 + Amat[Azz] - Amat[Axx] - Amat[Ayy] ) * 2.0;
      q[0] = (Amat[Axy] + Amat[Ayx]) / s;
      q[1] = (Amat[Axz] + Amat[Azx]) / s;
      q[2] = (Amat[Ayz] + Amat[Azy]) / s;
      q[3] = 0.5 / s;
    }
  }
}


void DirectionalAtom::setEuler( double phi, double theta, double psi ){
  
  Amat[Axx] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
  Amat[Axy] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
  Amat[Axz] = sin(theta) * sin(psi);
  
  Amat[Ayx] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
  Amat[Ayy] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
  Amat[Ayz] = sin(theta) * cos(psi);

  Amat[Azx] = sin(phi) * sin(theta);
  Amat[Azy] = -cos(phi) * sin(theta);
  Amat[Azz] = cos(theta);

  this->updateU();
}


void DirectionalAtom::lab2Body( double r[3] ){

  double rl[3]; // the lab frame vector 
  
  rl[0] = r[0];
  rl[1] = r[1];
  rl[2] = r[2];

  r[0] = (Amat[Axx] * rl[0]) + (Amat[Axy] * rl[1]) + (Amat[Axz] * rl[2]);
  r[1] = (Amat[Ayx] * rl[0]) + (Amat[Ayy] * rl[1]) + (Amat[Ayz] * rl[2]);
  r[2] = (Amat[Azx] * rl[0]) + (Amat[Azy] * rl[1]) + (Amat[Azz] * rl[2]);
}

void DirectionalAtom::body2Lab( double r[3] ){

  double rb[3]; // the body frame vector 
  
  rb[0] = r[0];
  rb[1] = r[1];
  rb[2] = r[2];

  r[0] = (Amat[Axx] * rb[0]) + (Amat[Ayx] * rb[1]) + (Amat[Azx] * rb[2]);
  r[1] = (Amat[Axy] * rb[0]) + (Amat[Ayy] * rb[1]) + (Amat[Azy] * rb[2]);
  r[2] = (Amat[Axz] * rb[0]) + (Amat[Ayz] * rb[1]) + (Amat[Azz] * rb[2]);
}

void DirectionalAtom::updateU( void ){

  ul[offsetX] = (Amat[Axx] * sux) + (Amat[Ayx] * suy) + (Amat[Azx] * suz);
  ul[offsetY] = (Amat[Axy] * sux) + (Amat[Ayy] * suy) + (Amat[Azy] * suz);
  ul[offsetZ] = (Amat[Axz] * sux) + (Amat[Ayz] * suy) + (Amat[Azz] * suz);
}

Revision:	414
Committed:	Wed Mar 26 22:02:36 2003 UTC (21 years, 3 months ago) by mmeineke
File size:	4566 byte(s)
Log Message:	the skeleton for making the molecules is in place. ForceField needs to be updated next.
#	User	Rev	Content
1	mmeineke	377	#include <cmath>
2
3			#include "Atom.hpp"
4
5
6	mmeineke	414
7	mmeineke	377	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];
12
13			this->updateU();
14			}
15
16			void DirectionalAtom::setI( double the_I[3][3] ){
17
18			Ixx = the_I[0][0]; Ixy = the_I[0][1]; Ixz = the_I[0][2];
19			Iyx = the_I[1][0]; Iyy = the_I[1][1]; Iyz = the_I[1][2];
20			Izx = the_I[2][0]; Izy = the_I[2][1]; Izz = the_I[2][2];
21			}
22
23			void DirectionalAtom::setQ( double the_q[4] ){
24
25			double q0Sqr, q1Sqr, q2Sqr, q3Sqr;
26
27			q0Sqr = the_q[0] * the_q[0];
28			q1Sqr = the_q[1] * the_q[1];
29			q2Sqr = the_q[2] * the_q[2];
30			q3Sqr = the_q[3] * the_q[3];
31
32
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();
46			}
47
48			void DirectionalAtom::getA( double the_A[3][3] ){
49
50			the_A[0][0] = Amat[Axx];
51			the_A[0][1] = Amat[Axy];
52			the_A[0][2] = Amat[Axz];
53
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];
61			}
62
63
64			void DirectionalAtom::getU( double the_u[3] ){
65
66			the_u[0] = sux;
67			the_u[1] = suy;
68			the_u[2] = suz;
69
70			this->body2Lab( the_u );
71			}
72
73			void DirectionalAtom::getQ( double q[4] ){
74
75			double t, s;
76			double ad1, ad2, ad3;
77
78			t = Amat[Axx] + Amat[Ayy] + Amat[Azz] + 1.0;
79			if( t > 0.0 ){
80
81			s = 0.5 / sqrt( t );
82			q[0] = 0.25 / s;
83			q[1] = (Amat[Ayz] - Amat[Azy]) * s;
84			q[2] = (Amat[Azx] - Amat[Axz]) * s;
85			q[3] = (Amat[Axy] - Amat[Ayx]) * s;
86			}
87			else{
88
89			ad1 = fabs( Amat[Axx] );
90			ad2 = fabs( Amat[Ayy] );
91			ad3 = fabs( Amat[Azz] );
92
93			if( ad1 >= ad2 && ad1 >= ad3 ){
94
95			s = 2.0 * sqrt( 1.0 + Amat[Axx] - Amat[Ayy] - Amat[Azz] );
96			q[0] = (Amat[Ayz] + Amat[Azy]) / s;
97			q[1] = 0.5 / s;
98			q[2] = (Amat[Axy] + Amat[Ayx]) / s;
99			q[3] = (Amat[Axz] + Amat[Azx]) / s;
100			}
101			else if( ad2 >= ad1 && ad2 >= ad3 ){
102
103			s = sqrt( 1.0 + Amat[Ayy] - Amat[Axx] - Amat[Azz] ) * 2.0;
104			q[0] = (Amat[Axz] + Amat[Azx]) / s;
105			q[1] = (Amat[Axy] + Amat[Ayx]) / s;
106			q[2] = 0.5 / s;
107			q[3] = (Amat[Ayz] + Amat[Azy]) / s;
108			}
109			else{
110
111			s = sqrt( 1.0 + Amat[Azz] - Amat[Axx] - Amat[Ayy] ) * 2.0;
112			q[0] = (Amat[Axy] + Amat[Ayx]) / s;
113			q[1] = (Amat[Axz] + Amat[Azx]) / s;
114			q[2] = (Amat[Ayz] + Amat[Azy]) / s;
115			q[3] = 0.5 / s;
116			}
117			}
118			}
119
120
121			void DirectionalAtom::setEuler( double phi, double theta, double psi ){
122
123			Amat[Axx] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
124			Amat[Axy] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
125			Amat[Axz] = sin(theta) * sin(psi);
126
127			Amat[Ayx] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
128			Amat[Ayy] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
129			Amat[Ayz] = sin(theta) * cos(psi);
130
131			Amat[Azx] = sin(phi) * sin(theta);
132			Amat[Azy] = -cos(phi) * sin(theta);
133			Amat[Azz] = cos(theta);
134
135			this->updateU();
136			}
137
138
139			void DirectionalAtom::lab2Body( double r[3] ){
140
141			double rl[3]; // the lab frame vector
142
143			rl[0] = r[0];
144			rl[1] = r[1];
145			rl[2] = r[2];
146
147			r[0] = (Amat[Axx] * rl[0]) + (Amat[Axy] * rl[1]) + (Amat[Axz] * rl[2]);
148			r[1] = (Amat[Ayx] * rl[0]) + (Amat[Ayy] * rl[1]) + (Amat[Ayz] * rl[2]);
149			r[2] = (Amat[Azx] * rl[0]) + (Amat[Azy] * rl[1]) + (Amat[Azz] * rl[2]);
150			}
151
152			void DirectionalAtom::body2Lab( double r[3] ){
153
154			double rb[3]; // the body frame vector
155
156			rb[0] = r[0];
157			rb[1] = r[1];
158			rb[2] = r[2];
159
160			r[0] = (Amat[Axx] * rb[0]) + (Amat[Ayx] * rb[1]) + (Amat[Azx] * rb[2]);
161			r[1] = (Amat[Axy] * rb[0]) + (Amat[Ayy] * rb[1]) + (Amat[Azy] * rb[2]);
162			r[2] = (Amat[Axz] * rb[0]) + (Amat[Ayz] * rb[1]) + (Amat[Azz] * rb[2]);
163			}
164
165			void DirectionalAtom::updateU( void ){
166
167			ul[offsetX] = (Amat[Axx] * sux) + (Amat[Ayx] * suy) + (Amat[Azx] * suz);
168			ul[offsetY] = (Amat[Axy] * sux) + (Amat[Ayy] * suy) + (Amat[Azy] * suz);
169			ul[offsetZ] = (Amat[Axz] * sux) + (Amat[Ayz] * suy) + (Amat[Azz] * suz);
170			}
171