mdtools/libmdCode/DirectionalAtom.cpp

#include <cmath>

#include "Atom.hpp"

double* Atom::pos; // the position array
double* Atom::vel; // the velocity array
double* Atom::frc; // the forc array
double* Atom::trq; // the torque vector  ( space fixed )
double* Atom::Amat; // the rotation matrix
double* Atom::mu;   // the array of dipole moments
double* Atom::ul;   // the lab frame unit directional vector

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 ){

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

Revision:	299
Committed:	Fri Mar 7 19:31:02 2003 UTC (22 years, 8 months ago) by mmeineke
File size:	4893 byte(s)
Log Message:	implemented the new fortran force interface
#	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			r[offsetX] = (Amat[Axx] * sux) + (Amat[Ayx] * suy) + (Amat[Azx] * suz);
174			r[offsetY] = (Amat[Axy] * sux) + (Amat[Ayy] * suy) + (Amat[Azy] * suz);
175			r[offsetZ] = (Amat[Axz] * sux) + (Amat[Ayz] * suy) + (Amat[Azz] * suz);
176			}
177