ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/OpenMD/branches/development/src/primitives/DirectionalAtom.cpp

Comparing:
trunk/src/primitives/DirectionalAtom.cpp (file contents), Revision 2 by gezelter, Fri Sep 24 04:16:43 2004 UTC vs.
branches/development/src/primitives/DirectionalAtom.cpp (file contents), Revision 1465 by chuckv, Fri Jul 9 23:08:25 2010 UTC

#	Line 1 | Line 1
1	<	#include <math.h>
2	<
3	<	#include "Atom.hpp"
4	<	#include "DirectionalAtom.hpp"
5	<	#include "simError.h"
6	<	#include "MatVec3.h"
7	<
8	<	void DirectionalAtom::zeroForces() {
9	<	if( hasCoords ){
10	<
11	<	Atom::zeroForces();
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. Redistributions of source code must retain the above copyright
10	>	*    notice, this list of conditions and the following disclaimer.
11	>	*
12	>	* 2. Redistributions in binary form must reproduce the above copyright
13	>	*    notice, this list of conditions and the following disclaimer in the
14	>	*    documentation and/or other materials provided with the
15	>	*    distribution.
16	>	*
17	>	* This software is provided "AS IS," without a warranty of any
18	>	* kind. All express or implied conditions, representations and
19	>	* warranties, including any implied warranty of merchantability,
20	>	* fitness for a particular purpose or non-infringement, are hereby
21	>	* excluded.  The University of Notre Dame and its licensors shall not
22	>	* be liable for any damages suffered by licensee as a result of
23	>	* using, modifying or distributing the software or its
24	>	* derivatives. In no event will the University of Notre Dame or its
25	>	* licensors be liable for any lost revenue, profit or data, or for
26	>	* direct, indirect, special, consequential, incidental or punitive
27	>	* damages, however caused and regardless of the theory of liability,
28	>	* arising out of the use of or inability to use software, even if the
29	>	* University of Notre Dame has been advised of the possibility of
30	>	* such damages.
31	>	*
32	>	* SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33	>	* research, please cite the appropriate papers when you publish your
34	>	* work.  Good starting points are:
35	>	*
36	>	* [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37	>	* [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38	>	* [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	>	* [4]  Vardeman & Gezelter, in progress (2009).
40	>	*/
41	>
42	>	#include "primitives/DirectionalAtom.hpp"
43	>	#include "utils/simError.h"
44	>	namespace OpenMD {
45	>
46	>	DirectionalAtom::DirectionalAtom(DirectionalAtomType* dAtomType)
47	>	: Atom(dAtomType){
48	>	objType_= otDAtom;
49	>	if (dAtomType->isMultipole()) {
50	>	electroBodyFrame_ = dAtomType->getElectroBodyFrame();
51	>	}
52
53	<	trq[offsetX] = 0.0;
54	<	trq[offsetY] = 0.0;
55	<	trq[offsetZ] = 0.0;
56	<	}
57	<	else{
58	<
59	<	sprintf( painCave.errMsg,
60	<	"Attempt to zero frc and trq for atom %d before coords set.\n",
61	<	index );
62	<	painCave.isFatal = 1;
63	<	simError();
24	<	}
25	<	}
53	>	// Check if one of the diagonal inertia tensor of this directional
54	>	// atom is zero:
55	>	int nLinearAxis = 0;
56	>	Mat3x3d inertiaTensor = getI();
57	>	for (int i = 0; i < 3; i++) {
58	>	if (fabs(inertiaTensor(i, i)) < OpenMD::epsilon) {
59	>	linear_ = true;
60	>	linearAxis_ = i;
61	>	++ nLinearAxis;
62	>	}
63	>	}
64
65	<	void DirectionalAtom::setCoords(void){
66	<
67	<	if( myConfig->isAllocated() ){
68	<
69	<	myConfig->getAtomPointers( index,
70	<	&pos,
71	<	&vel,
72	<	&frc,
35	<	&trq,
36	<	&Amat,
37	<	&mu,
38	<	&ul);
65	>	if (nLinearAxis > 1) {
66	>	sprintf( painCave.errMsg,
67	>	"Directional Atom warning.\n"
68	>	"\tOpenMD found more than one axis in this directional atom with a vanishing \n"
69	>	"\tmoment of inertia.");
70	>	painCave.isFatal = 0;
71	>	simError();
72	>	}
73		}
74	<	else{
75	<	sprintf( painCave.errMsg,
76	<	"Attempted to set Atom %d  coordinates with an unallocated "
77	<	"SimState object.\n", index );
78	<	painCave.isFatal = 1;
79	<	simError();
74	>
75	>	Mat3x3d DirectionalAtom::getI() {
76	>	return static_cast<DirectionalAtomType*>(getAtomType())->getI();
77	>	}
78	>
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	>	}
84		}
85	<
86	<	hasCoords = true;
87	<
88	<	}
51	<
52	<	void DirectionalAtom::setA( double the_A[3][3] ){
53	<
54	<	if( hasCoords ){
55	<	Amat[Axx] = the_A[0][0]; Amat[Axy] = the_A[0][1]; Amat[Axz] = the_A[0][2];
56	<	Amat[Ayx] = the_A[1][0]; Amat[Ayy] = the_A[1][1]; Amat[Ayz] = the_A[1][2];
57	<	Amat[Azx] = the_A[2][0]; Amat[Azy] = the_A[2][1]; Amat[Azz] = the_A[2][2];
85	>
86	>
87	>	void DirectionalAtom::setA(const RotMat3x3d& a) {
88	>	((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
89
90	<	this->updateU();
91	<	}
92	<	else{
90	>	if (atomType_->isMultipole()) {
91	>	((snapshotMan_->getCurrentSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;
92	>	}
93	>	}
94	>
95	>	void DirectionalAtom::setA(const RotMat3x3d& a, int snapshotNo) {
96	>	((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a;
97
98	<	sprintf( painCave.errMsg,
99	<	"Attempt to set Amat for atom %d before coords set.\n",
100	<	index );
101	<	painCave.isFatal = 1;
102	<	simError();
98	>	if (atomType_->isMultipole()) {
99	>	((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;
100	>	}
101	>	}
102	>
103	>	void DirectionalAtom::rotateBy(const RotMat3x3d& m) {
104	>	setA(m *getA());
105		}
69	–	}
70	–
71	–	void DirectionalAtom::setI( double the_I[3][3] ){
106
107	<	Ixx = the_I[0][0]; Ixy = the_I[0][1]; Ixz = the_I[0][2];
108	<	Iyx = the_I[1][0]; Iyy = the_I[1][1]; Iyz = the_I[1][2];
109	<	Izx = the_I[2][0]; Izy = the_I[2][1]; Izz = the_I[2][2];
110	<	}
111	<
112	<	void DirectionalAtom::setQ( double the_q[4] ){
113	<
114	<	double q0Sqr, q1Sqr, q2Sqr, q3Sqr;
115	<
116	<	if( hasCoords ){
83	<	q0Sqr = the_q[0] * the_q[0];
84	<	q1Sqr = the_q[1] * the_q[1];
85	<	q2Sqr = the_q[2] * the_q[2];
86	<	q3Sqr = the_q[3] * the_q[3];
107	>	std::vector<RealType> DirectionalAtom::getGrad() {
108	>	std::vector<RealType> grad(6, 0.0);
109	>	Vector3d force;
110	>	Vector3d torque;
111	>	Vector3d myEuler;
112	>	RealType phi, theta, psi;
113	>	RealType cphi, sphi, ctheta, stheta;
114	>	Vector3d ephi;
115	>	Vector3d etheta;
116	>	Vector3d epsi;
117
118	+	force = getFrc();
119	+	torque =getTrq();
120	+	myEuler = getA().toEulerAngles();
121
122	<	Amat[Axx] = q0Sqr + q1Sqr - q2Sqr - q3Sqr;
123	<	Amat[Axy] = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] );
124	<	Amat[Axz] = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] );
122	>	phi = myEuler[0];
123	>	theta = myEuler[1];
124	>	psi = myEuler[2];
125
126	<	Amat[Ayx] = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] );
127	<	Amat[Ayy] = q0Sqr - q1Sqr + q2Sqr - q3Sqr;
128	<	Amat[Ayz] = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] );
126	>	cphi = cos(phi);
127	>	sphi = sin(phi);
128	>	ctheta = cos(theta);
129	>	stheta = sin(theta);
130
131	<	Amat[Azx] = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] );
98	<	Amat[Azy] = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] );
99	<	Amat[Azz] = q0Sqr - q1Sqr -q2Sqr +q3Sqr;
131	>	// get unit vectors along the phi, theta and psi rotation axes
132
133	<	this->updateU();
134	<	}
135	<	else{
133	>	ephi[0] = 0.0;
134	>	ephi[1] = 0.0;
135	>	ephi[2] = 1.0;
136
137	<	sprintf( painCave.errMsg,
138	<	"Attempt to set Q for atom %d before coords set.\n",
139	<	index );
108	<	painCave.isFatal = 1;
109	<	simError();
110	<	}
111	<
112	<	}
113	<
114	<	void DirectionalAtom::getA( double the_A[3][3] ){
115	<
116	<	if( hasCoords ){
117	<	the_A[0][0] = Amat[Axx];
118	<	the_A[0][1] = Amat[Axy];
119	<	the_A[0][2] = Amat[Axz];
137	>	//etheta[0] = -sphi;
138	>	//etheta[1] =  cphi;
139	>	//etheta[2] =  0.0;
140
141	<	the_A[1][0] = Amat[Ayx];
142	<	the_A[1][1] = Amat[Ayy];
143	<	the_A[1][2] = Amat[Ayz];
141	>	etheta[0] = cphi;
142	>	etheta[1] = sphi;
143	>	etheta[2] = 0.0;
144
145	<	the_A[2][0] = Amat[Azx];
146	<	the_A[2][1] = Amat[Azy];
147	<	the_A[2][2] = Amat[Azz];
128	<	}
129	<	else{
145	>	epsi[0] = stheta * cphi;
146	>	epsi[1] = stheta * sphi;
147	>	epsi[2] = ctheta;
148
149	<	sprintf( painCave.errMsg,
150	<	"Attempt to get Amat for atom %d before coords set.\n",
151	<	index );
134	<	painCave.isFatal = 1;
135	<	simError();
136	<	}
137	<
138	<	}
139	<
140	<	void DirectionalAtom::printAmatIndex( void ){
141	<
142	<	if( hasCoords ){
143	<	std::cerr << "Atom[" << index << "] index =>\n"
144	<	<< "[ " << Axx << ", " << Axy << ", " << Axz << " ]\n"
145	<	<< "[ " << Ayx << ", " << Ayy << ", " << Ayz << " ]\n"
146	<	<< "[ " << Azx << ", " << Azy << ", " << Azz << " ]\n";
147	<	}
148	<	else{
149	>	//gradient is equal to -force
150	>	for (int j = 0 ; j<3; j++)
151	>	grad[j] = -force[j];
152
153	<	sprintf( painCave.errMsg,
154	<	"Attempt to print Amat indices for atom %d before coords set.\n",
155	<	index );
156	<	painCave.isFatal = 1;
154	<	simError();
155	<	}
156	<	}
157	<
158	<
159	<	void DirectionalAtom::getU( double the_u[3] ){
160	<
161	<	the_u[0] = sU[2][0];
162	<	the_u[1] = sU[2][1];
163	<	the_u[2] = sU[2][2];
164	<
165	<	this->body2Lab( the_u );
166	<	}
167	<
168	<	void DirectionalAtom::getQ( double q[4] ){
169	<
170	<	double t, s;
171	<	double ad1, ad2, ad3;
172	<
173	<	if( hasCoords ){
174	<
175	<	t = Amat[Axx] + Amat[Ayy] + Amat[Azz] + 1.0;
176	<	if( t > 0.0 ){
177	<
178	<	s = 0.5 / sqrt( t );
179	<	q[0] = 0.25 / s;
180	<	q[1] = (Amat[Ayz] - Amat[Azy]) * s;
181	<	q[2] = (Amat[Azx] - Amat[Axz]) * s;
182	<	q[3] = (Amat[Axy] - Amat[Ayx]) * s;
153	>	for (int j = 0; j < 3; j++ ) {
154	>	grad[3] -= torque[j]*ephi[j];
155	>	grad[4] -= torque[j]*etheta[j];
156	>	grad[5] -= torque[j]*epsi[j];
157		}
184	–	else{
185	–
186	–	ad1 = fabs( Amat[Axx] );
187	–	ad2 = fabs( Amat[Ayy] );
188	–	ad3 = fabs( Amat[Azz] );
189	–
190	–	if( ad1 >= ad2 && ad1 >= ad3 ){
191	–
192	–	s = 2.0 * sqrt( 1.0 + Amat[Axx] - Amat[Ayy] - Amat[Azz] );
193	–	q[0] = (Amat[Ayz] + Amat[Azy]) / s;
194	–	q[1] = 0.5 / s;
195	–	q[2] = (Amat[Axy] + Amat[Ayx]) / s;
196	–	q[3] = (Amat[Axz] + Amat[Azx]) / s;
197	–	}
198	–	else if( ad2 >= ad1 && ad2 >= ad3 ){
199	–
200	–	s = sqrt( 1.0 + Amat[Ayy] - Amat[Axx] - Amat[Azz] ) * 2.0;
201	–	q[0] = (Amat[Axz] + Amat[Azx]) / s;
202	–	q[1] = (Amat[Axy] + Amat[Ayx]) / s;
203	–	q[2] = 0.5 / s;
204	–	q[3] = (Amat[Ayz] + Amat[Azy]) / s;
205	–	}
206	–	else{
207	–
208	–	s = sqrt( 1.0 + Amat[Azz] - Amat[Axx] - Amat[Ayy] ) * 2.0;
209	–	q[0] = (Amat[Axy] + Amat[Ayx]) / s;
210	–	q[1] = (Amat[Axz] + Amat[Azx]) / s;
211	–	q[2] = (Amat[Ayz] + Amat[Azy]) / s;
212	–	q[3] = 0.5 / s;
213	–	}
214	–	}
215	–	}
216	–	else{
158
159	<	sprintf( painCave.errMsg,
160	<	"Attempt to get Q for atom %d before coords set.\n",
220	<	index );
221	<	painCave.isFatal = 1;
222	<	simError();
223	<	}
224	<	}
225	<
226	<	void DirectionalAtom::setUnitFrameFromEuler(double phi,
227	<	double theta,
228	<	double psi) {
229	<
230	<	double myA[3][3];
231	<	double uFrame[3][3];
232	<	double len;
233	<	int i, j;
234	<
235	<	myA[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
236	<	myA[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
237	<	myA[0][2] = sin(theta) * sin(psi);
238	<
239	<	myA[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
240	<	myA[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
241	<	myA[1][2] = sin(theta) * cos(psi);
242	<
243	<	myA[2][0] = sin(phi) * sin(theta);
244	<	myA[2][1] = -cos(phi) * sin(theta);
245	<	myA[2][2] = cos(theta);
246	<
247	<	// Make the unit Frame:
248	<
249	<	for (i=0; i < 3; i++)
250	<	for (j=0; j < 3; j++)
251	<	uFrame[i][j] = 0.0;
252	<
253	<	for (i=0; i < 3; i++)
254	<	uFrame[i][i] = 1.0;
255	<
256	<	// rotate by the given rotation matrix:
257	<
258	<	matMul3(myA, uFrame, sU);
259	<
260	<	// renormalize column vectors:
261	<
262	<	for (i=0; i < 3; i++) {
263	<	len = 0.0;
264	<	for (j = 0; j < 3; j++) {
265	<	len += sU[i][j]*sU[i][j];
266	<	}
267	<	len = sqrt(len);
268	<	for (j = 0; j < 3; j++) {
269	<	sU[i][j] /= len;
270	<	}
271	<	}
272	<
273	<	// sU now contains the coordinates of the 'special' frame;
274	<
275	<	}
276	<
277	<	void DirectionalAtom::setEuler( double phi, double theta, double psi ){
159	>	return grad;
160	>	}
161
162	<	if( hasCoords ){
163	<	Amat[Axx] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
281	<	Amat[Axy] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
282	<	Amat[Axz] = sin(theta) * sin(psi);
283	<
284	<	Amat[Ayx] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
285	<	Amat[Ayy] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
286	<	Amat[Ayz] = sin(theta) * cos(psi);
287	<
288	<	Amat[Azx] = sin(phi) * sin(theta);
289	<	Amat[Azy] = -cos(phi) * sin(theta);
290	<	Amat[Azz] = cos(theta);
291	<
292	<	this->updateU();
162	>	void DirectionalAtom::accept(BaseVisitor* v) {
163	>	v->visit(this);
164		}
294	–	else{
295	–
296	–	sprintf( painCave.errMsg,
297	–	"Attempt to set Euler angles for atom %d before coords set.\n",
298	–	index );
299	–	painCave.isFatal = 1;
300	–	simError();
301	–	}
165		}
166
304	–
305	–	void DirectionalAtom::lab2Body( double r[3] ){
306	–
307	–	double rl[3]; // the lab frame vector
308	–
309	–	if( hasCoords ){
310	–	rl[0] = r[0];
311	–	rl[1] = r[1];
312	–	rl[2] = r[2];
313	–
314	–	r[0] = (Amat[Axx] * rl[0]) + (Amat[Axy] * rl[1]) + (Amat[Axz] * rl[2]);
315	–	r[1] = (Amat[Ayx] * rl[0]) + (Amat[Ayy] * rl[1]) + (Amat[Ayz] * rl[2]);
316	–	r[2] = (Amat[Azx] * rl[0]) + (Amat[Azy] * rl[1]) + (Amat[Azz] * rl[2]);
317	–	}
318	–	else{
319	–
320	–	sprintf( painCave.errMsg,
321	–	"Attempt to convert lab2body for atom %d before coords set.\n",
322	–	index );
323	–	painCave.isFatal = 1;
324	–	simError();
325	–	}
326	–
327	–	}
328	–
329	–	void DirectionalAtom::rotateBy( double by_A[3][3]) {
330	–
331	–	// Check this
332	–
333	–	double r00, r01, r02, r10, r11, r12, r20, r21, r22;
334	–
335	–	if( hasCoords ){
336	–
337	–	r00 = by_A[0][0]*Amat[Axx] + by_A[0][1]*Amat[Ayx] + by_A[0][2]*Amat[Azx];
338	–	r01 = by_A[0][0]*Amat[Axy] + by_A[0][1]*Amat[Ayy] + by_A[0][2]*Amat[Azy];
339	–	r02 = by_A[0][0]*Amat[Axz] + by_A[0][1]*Amat[Ayz] + by_A[0][2]*Amat[Azz];
340	–
341	–	r10 = by_A[1][0]*Amat[Axx] + by_A[1][1]*Amat[Ayx] + by_A[1][2]*Amat[Azx];
342	–	r11 = by_A[1][0]*Amat[Axy] + by_A[1][1]*Amat[Ayy] + by_A[1][2]*Amat[Azy];
343	–	r12 = by_A[1][0]*Amat[Axz] + by_A[1][1]*Amat[Ayz] + by_A[1][2]*Amat[Azz];
344	–
345	–	r20 = by_A[2][0]*Amat[Axx] + by_A[2][1]*Amat[Ayx] + by_A[2][2]*Amat[Azx];
346	–	r21 = by_A[2][0]*Amat[Axy] + by_A[2][1]*Amat[Ayy] + by_A[2][2]*Amat[Azy];
347	–	r22 = by_A[2][0]*Amat[Axz] + by_A[2][1]*Amat[Ayz] + by_A[2][2]*Amat[Azz];
348	–
349	–	Amat[Axx] = r00; Amat[Axy] = r01; Amat[Axz] = r02;
350	–	Amat[Ayx] = r10; Amat[Ayy] = r11; Amat[Ayz] = r12;
351	–	Amat[Azx] = r20; Amat[Azy] = r21; Amat[Azz] = r22;
352	–
353	–	}
354	–	else{
355	–
356	–	sprintf( painCave.errMsg,
357	–	"Attempt to rotate frame for atom %d before coords set.\n",
358	–	index );
359	–	painCave.isFatal = 1;
360	–	simError();
361	–	}
362	–
363	–	}
364	–
365	–
366	–	void DirectionalAtom::body2Lab( double r[3] ){
367	–
368	–	double rb[3]; // the body frame vector
369	–
370	–	if( hasCoords ){
371	–	rb[0] = r[0];
372	–	rb[1] = r[1];
373	–	rb[2] = r[2];
374	–
375	–	r[0] = (Amat[Axx] * rb[0]) + (Amat[Ayx] * rb[1]) + (Amat[Azx] * rb[2]);
376	–	r[1] = (Amat[Axy] * rb[0]) + (Amat[Ayy] * rb[1]) + (Amat[Azy] * rb[2]);
377	–	r[2] = (Amat[Axz] * rb[0]) + (Amat[Ayz] * rb[1]) + (Amat[Azz] * rb[2]);
378	–	}
379	–	else{
380	–
381	–	sprintf( painCave.errMsg,
382	–	"Attempt to convert body2lab for atom %d before coords set.\n",
383	–	index );
384	–	painCave.isFatal = 1;
385	–	simError();
386	–	}
387	–	}
388	–
389	–	void DirectionalAtom::updateU( void ){
390	–
391	–	if( hasCoords ){
392	–	ul[offsetX] = (Amat[Axx] * sU[2][0]) +
393	–	(Amat[Ayx] * sU[2][1]) + (Amat[Azx] * sU[2][2]);
394	–	ul[offsetY] = (Amat[Axy] * sU[2][0]) +
395	–	(Amat[Ayy] * sU[2][1]) + (Amat[Azy] * sU[2][2]);
396	–	ul[offsetZ] = (Amat[Axz] * sU[2][0]) +
397	–	(Amat[Ayz] * sU[2][1]) + (Amat[Azz] * sU[2][2]);
398	–	}
399	–	else{
400	–
401	–	sprintf( painCave.errMsg,
402	–	"Attempt to updateU for atom %d before coords set.\n",
403	–	index );
404	–	painCave.isFatal = 1;
405	–	simError();
406	–	}
407	–	}
408	–
409	–	void DirectionalAtom::getJ( double theJ[3] ){
410	–
411	–	theJ[0] = jx;
412	–	theJ[1] = jy;
413	–	theJ[2] = jz;
414	–	}
415	–
416	–	void DirectionalAtom::setJ( double theJ[3] ){
417	–
418	–	jx = theJ[0];
419	–	jy = theJ[1];
420	–	jz = theJ[2];
421	–	}
422	–
423	–	void DirectionalAtom::getTrq( double theT[3] ){
424	–
425	–	if( hasCoords ){
426	–	theT[0] = trq[offsetX];
427	–	theT[1] = trq[offsetY];
428	–	theT[2] = trq[offsetZ];
429	–	}
430	–	else{
431	–
432	–	sprintf( painCave.errMsg,
433	–	"Attempt to get Trq for atom %d before coords set.\n",
434	–	index );
435	–	painCave.isFatal = 1;
436	–	simError();
437	–	}
438	–	}
439	–
440	–	void DirectionalAtom::addTrq( double theT[3] ){
441	–
442	–	if( hasCoords ){
443	–	trq[offsetX] += theT[0];
444	–	trq[offsetY] += theT[1];
445	–	trq[offsetZ] += theT[2];
446	–	}
447	–	else{
448	–
449	–	sprintf( painCave.errMsg,
450	–	"Attempt to add Trq for atom %d before coords set.\n",
451	–	index );
452	–	painCave.isFatal = 1;
453	–	simError();
454	–	}
455	–	}
456	–
457	–
458	–	void DirectionalAtom::getI( double the_I[3][3] ){
459	–
460	–	the_I[0][0] = Ixx;
461	–	the_I[0][1] = Ixy;
462	–	the_I[0][2] = Ixz;
463	–
464	–	the_I[1][0] = Iyx;
465	–	the_I[1][1] = Iyy;
466	–	the_I[1][2] = Iyz;
467	–
468	–	the_I[2][0] = Izx;
469	–	the_I[2][1] = Izy;
470	–	the_I[2][2] = Izz;
471	–	}
472	–
473	–	void DirectionalAtom::getGrad( double grad[6] ) {
474	–
475	–	double myEuler[3];
476	–	double phi, theta, psi;
477	–	double cphi, sphi, ctheta, stheta;
478	–	double ephi[3];
479	–	double etheta[3];
480	–	double epsi[3];
481	–
482	–	this->getEulerAngles(myEuler);
483	–
484	–	phi = myEuler[0];
485	–	theta = myEuler[1];
486	–	psi = myEuler[2];
487	–
488	–	cphi = cos(phi);
489	–	sphi = sin(phi);
490	–	ctheta = cos(theta);
491	–	stheta = sin(theta);
492	–
493	–	// get unit vectors along the phi, theta and psi rotation axes
494	–
495	–	ephi[0] = 0.0;
496	–	ephi[1] = 0.0;
497	–	ephi[2] = 1.0;
498	–
499	–	etheta[0] = cphi;
500	–	etheta[1] = sphi;
501	–	etheta[2] = 0.0;
502	–
503	–	epsi[0] = stheta * cphi;
504	–	epsi[1] = stheta * sphi;
505	–	epsi[2] = ctheta;
506	–
507	–	for (int j = 0 ; j<3; j++)
508	–	grad[j] = frc[j];
509	–
510	–	grad[3] = 0;
511	–	grad[4] = 0;
512	–	grad[5] = 0;
513	–
514	–	for (int j = 0; j < 3; j++ ) {
515	–
516	–	grad[3] += trq[j]*ephi[j];
517	–	grad[4] += trq[j]*etheta[j];
518	–	grad[5] += trq[j]*epsi[j];
519	–
520	–	}
521	–
522	–	}
523	–
524	–	/**
525	–	* getEulerAngles computes a set of Euler angle values consistent
526	–	*  with an input rotation matrix.  They are returned in the following
527	–	* order:
528	–	*  myEuler[0] = phi;
529	–	*  myEuler[1] = theta;
530	–	*  myEuler[2] = psi;
531	–	*/
532	–	void DirectionalAtom::getEulerAngles(double myEuler[3]) {
533	–
534	–	// We use so-called "x-convention", which is the most common definition.
535	–	// In this convention, the rotation given by Euler angles (phi, theta, psi), where the first
536	–	// rotation is by an angle phi about the z-axis, the second is by an angle
537	–	// theta (0 <= theta <= 180)about the x-axis, and thethird is by an angle psi about the
538	–	//z-axis (again).
539	–
540	–
541	–	double phi,theta,psi,eps;
542	–	double ctheta,stheta;
543	–
544	–	// set the tolerance for Euler angles and rotation elements
545	–
546	–	eps = 1.0e-8;
547	–
548	–	theta = acos(min(1.0,max(-1.0,Amat[Azz])));
549	–	ctheta = Amat[Azz];
550	–	stheta = sqrt(1.0 - ctheta * ctheta);
551	–
552	–	// when sin(theta) is close to 0, we need to consider singularity
553	–	// In this case, we can assign an arbitary value to phi (or psi), and then determine
554	–	// the psi (or phi) or vice-versa. We'll assume that phi always gets the rotation, and psi is 0
555	–	// in cases of singularity.
556	–	// we use atan2 instead of atan, since atan2 will give us -Pi to Pi.
557	–	// Since 0 <= theta <= 180, sin(theta) will be always non-negative. Therefore, it never
558	–	// change the sign of both of the parameters passed to atan2.
559	–
560	–	if (fabs(stheta) <= eps){
561	–	psi = 0.0;
562	–	phi = atan2(-Amat[Ayx], Amat[Axx]);
563	–	}
564	–	// we only have one unique solution
565	–	else{
566	–	phi = atan2(Amat[Azx], -Amat[Azy]);
567	–	psi = atan2(Amat[Axz], Amat[Ayz]);
568	–	}
569	–
570	–	//wrap phi and psi, make sure they are in the range from 0 to 2*Pi
571	–	//if (phi < 0)
572	–	//  phi += M_PI;
573	–
574	–	//if (psi < 0)
575	–	//  psi += M_PI;
576	–
577	–	myEuler[0] = phi;
578	–	myEuler[1] = theta;
579	–	myEuler[2] = psi;
580	–
581	–	return;
582	–	}
583	–
584	–	double DirectionalAtom::getZangle( ){
585	–
586	–	if( hasCoords ){
587	–	return zAngle;
588	–	}
589	–	else{
590	–
591	–	sprintf( painCave.errMsg,
592	–	"Attempt to get zAngle for atom %d before coords set.\n",
593	–	index );
594	–	painCave.isFatal = 1;
595	–	simError();
596	–	return 0;
597	–	}
598	–	}
599	–
600	–	void DirectionalAtom::setZangle( double zAng ){
601	–
602	–	if( hasCoords ){
603	–	zAngle = zAng;
604	–	}
605	–	else{
606	–
607	–	sprintf( painCave.errMsg,
608	–	"Attempt to set zAngle for atom %d before coords set.\n",
609	–	index );
610	–	painCave.isFatal = 1;
611	–	simError();
612	–	}
613	–	}
614	–
615	–	void DirectionalAtom::addZangle( double zAng ){
616	–
617	–	if( hasCoords ){
618	–	zAngle += zAng;
619	–	}
620	–	else{
621	–
622	–	sprintf( painCave.errMsg,
623	–	"Attempt to add zAngle to atom %d before coords set.\n",
624	–	index );
625	–	painCave.isFatal = 1;
626	–	simError();
627	–	}
628	–	}
629	–
630	–	double DirectionalAtom::max(double x, double  y) {
631	–	return (x > y) ? x : y;
632	–	}
633	–
634	–	double DirectionalAtom::min(double x, double  y) {
635	–	return (x > y) ? y : x;
636	–	}

Comparing:
trunk/src/primitives/DirectionalAtom.cpp (property svn:keywords), Revision 2 by gezelter, Fri Sep 24 04:16:43 2004 UTC vs.
branches/development/src/primitives/DirectionalAtom.cpp (property svn:keywords), Revision 1465 by chuckv, Fri Jul 9 23:08:25 2010 UTC

#	Line 0 | Line 1
1	+	Author Id Revision Date

Diff Legend

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