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 205 by gezelter, Thu Nov 4 16:22:03 2004 UTC vs.
branches/development/src/primitives/DirectionalAtom.cpp (file contents), Revision 1787 by gezelter, Wed Aug 29 18:13:11 2012 UTC

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

Comparing:
trunk/src/primitives/DirectionalAtom.cpp (property svn:keywords), Revision 205 by gezelter, Thu Nov 4 16:22:03 2004 UTC vs.
branches/development/src/primitives/DirectionalAtom.cpp (property svn:keywords), Revision 1787 by gezelter, Wed Aug 29 18:13:11 2012 UTC

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

Diff Legend

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