ViewVC Help

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

root/OpenMD/trunk/src/math/SquareMatrix3.hpp

Comparing trunk/src/math/SquareMatrix3.hpp (file contents):
Revision 633 by tim, Tue Sep 27 20:02:57 2005 UTC vs.
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC

#	Line 6 | Line 6
6		* redistribute this software in source and binary code form, provided
7		* that the following conditions are met:
8		*
9	<	* 1. Acknowledgement of the program authors must be made in any
10	<	*    publication of scientific results based in part on use of the
11	<	*    program.  An acceptable form of acknowledgement is citation of
12	<	*    the article in which the program was described (Matthew
13	<	*    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	<	*    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	<	*    Parallel Simulation Engine for Molecular Dynamics,"
16	<	*    J. Comput. Chem. 26, pp. 252-271 (2005))
17	<	*
18	<	* 2. Redistributions of source code must retain the above copyright
9	>	* 1. Redistributions of source code must retain the above copyright
10		*    notice, this list of conditions and the following disclaimer.
11		*
12	<	* 3. Redistributions in binary form must reproduce the above copyright
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.
#	Line 37 | Line 28
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, 234107 (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		/**
#	Line 47 | Line 48
48		*/
49		#ifndef MATH_SQUAREMATRIX3_HPP
50		#define  MATH_SQUAREMATRIX3_HPP
51	<
51	>	#include "config.h"
52	>	#include <cmath>
53	>	#include <vector>
54		#include "Quaternion.hpp"
55		#include "SquareMatrix.hpp"
56		#include "Vector3.hpp"
57		#include "utils/NumericConstant.hpp"
58	<	namespace oopse {
58	>	namespace OpenMD {
59
60		template<typename Real>
61		class SquareMatrix3 : public SquareMatrix<Real, 3> {
#	Line 109 | Line 112 | namespace oopse {
112		return *this;
113		}
114
115	+
116		/**
117		* Sets this matrix to a rotation matrix by three euler angles
118		* @ param euler
#	Line 121 | Line 125 | namespace oopse {
125		* Sets this matrix to a rotation matrix by three euler angles
126		* @param phi
127		* @param theta
128	<	* @psi theta
128	>	* @param psi
129		*/
130		void setupRotMat(Real phi, Real theta, Real psi) {
131		Real sphi, stheta, spsi;
#	Line 168 | Line 172 | namespace oopse {
172		*this = q.toRotationMatrix3();
173		}
174
175	+	void setupSkewMat(Vector3<Real> v) {
176	+	setupSkewMat(v[0], v[1], v[2]);
177	+	}
178	+
179	+	void setupSkewMat(Real v1, Real v2, Real v3) {
180	+	this->data_[0][0] = 0;
181	+	this->data_[0][1] = -v3;
182	+	this->data_[0][2] = v2;
183	+	this->data_[1][0] = v3;
184	+	this->data_[1][1] = 0;
185	+	this->data_[1][2] = -v1;
186	+	this->data_[2][0] = -v2;
187	+	this->data_[2][1] = v1;
188	+	this->data_[2][2] = 0;
189	+
190	+
191	+	}
192	+
193	+
194		/**
195		* Returns the quaternion from this rotation matrix
196		* @return the quaternion from this rotation matrix
#	Line 179 | Line 202 | namespace oopse {
202		Real ad1, ad2, ad3;
203		t = this->data_[0][0] + this->data_[1][1] + this->data_[2][2] + 1.0;
204
205	<	if( t > 0.0 ){
205	>	if( t > NumericConstant::epsilon ){
206
207		s = 0.5 / sqrt( t );
208		q[0] = 0.25 / s;
#	Line 224 | Line 247 | namespace oopse {
247		* @return the euler angles in a vector
248		* @exception invalid rotation matrix
249		* We use so-called "x-convention", which is the most common definition.
250	<	* In this convention, the rotation given by Euler angles (phi, theta, psi), where the first
251	<	* rotation is by an angle phi about the z-axis, the second is by an angle
252	<	* theta (0 <= theta <= 180)about the x-axis, and thethird is by an angle psi about the
253	<	* z-axis (again).
250	>	* In this convention, the rotation given by Euler angles (phi, theta,
251	>	* psi), where the first rotation is by an angle phi about the z-axis,
252	>	* the second is by an angle theta (0 <= theta <= 180) about the x-axis,
253	>	* and the third is by an angle psi about the z-axis (again).
254		*/
255		Vector3<Real> toEulerAngles() {
256		Vector3<Real> myEuler;
#	Line 239 | Line 262 | namespace oopse {
262
263		// set the tolerance for Euler angles and rotation elements
264
265	<	theta = acos(std::min(1.0, std::max(-1.0,this->data_[2][2])));
265	>	theta = acos(std::min((RealType)1.0, std::max((RealType)-1.0,this->data_[2][2])));
266		ctheta = this->data_[2][2];
267		stheta = sqrt(1.0 - ctheta * ctheta);
268
269	<	// when sin(theta) is close to 0, we need to consider singularity
270	<	// In this case, we can assign an arbitary value to phi (or psi), and then determine
271	<	// the psi (or phi) or vice-versa. We'll assume that phi always gets the rotation, and psi is 0
272	<	// in cases of singularity.
269	>	// when sin(theta) is close to 0, we need to consider
270	>	// singularity In this case, we can assign an arbitary value to
271	>	// phi (or psi), and then determine the psi (or phi) or
272	>	// vice-versa. We'll assume that phi always gets the rotation,
273	>	// and psi is 0 in cases of singularity.
274		// we use atan2 instead of atan, since atan2 will give us -Pi to Pi.
275	<	// Since 0 <= theta <= 180, sin(theta) will be always non-negative. Therefore, it never
276	<	// change the sign of both of the parameters passed to atan2.
275	>	// Since 0 <= theta <= 180, sin(theta) will be always
276	>	// non-negative. Therefore, it will never change the sign of both of
277	>	// the parameters passed to atan2.
278
279	<	if (fabs(stheta) <= oopse::epsilon){
279	>	if (fabs(stheta) < 1e-6){
280		psi = 0.0;
281		phi = atan2(-this->data_[1][0], this->data_[0][0]);
282		}
#	Line 263 | Line 288 | namespace oopse {
288
289		//wrap phi and psi, make sure they are in the range from 0 to 2*Pi
290		if (phi < 0)
291	<	phi += M_PI;
291	>	phi += 2.0 * M_PI;
292
293		if (psi < 0)
294	<	psi += M_PI;
294	>	psi += 2.0 * M_PI;
295
296		myEuler[0] = phi;
297		myEuler[1] = theta;
#	Line 298 | Line 323 | namespace oopse {
323		*/
324		SquareMatrix3<Real>  inverse() const {
325		SquareMatrix3<Real> m;
326	<	double det = determinant();
327	<	if (fabs(det) <= oopse::epsilon) {
326	>	RealType det = determinant();
327	>	if (fabs(det) <= OpenMD::epsilon) {
328		//"The method was called on a matrix with |determinant| <= 1e-6.",
329		//"This is a runtime or a programming error in your application.");
330	<	}
330	>	std::vector<int> zeroDiagElementIndex;
331	>	for (int i =0; i < 3; ++i) {
332	>	if (fabs(this->data_[i][i]) <= OpenMD::epsilon) {
333	>	zeroDiagElementIndex.push_back(i);
334	>	}
335	>	}
336
337	<	m(0, 0) = this->data_[1][1]*this->data_[2][2] - this->data_[1][2]*this->data_[2][1];
338	<	m(1, 0) = this->data_[1][2]*this->data_[2][0] - this->data_[1][0]*this->data_[2][2];
339	<	m(2, 0) = this->data_[1][0]*this->data_[2][1] - this->data_[1][1]*this->data_[2][0];
340	<	m(0, 1) = this->data_[2][1]*this->data_[0][2] - this->data_[2][2]*this->data_[0][1];
311	<	m(1, 1) = this->data_[2][2]*this->data_[0][0] - this->data_[2][0]*this->data_[0][2];
312	<	m(2, 1) = this->data_[2][0]*this->data_[0][1] - this->data_[2][1]*this->data_[0][0];
313	<	m(0, 2) = this->data_[0][1]*this->data_[1][2] - this->data_[0][2]*this->data_[1][1];
314	<	m(1, 2) = this->data_[0][2]*this->data_[1][0] - this->data_[0][0]*this->data_[1][2];
315	<	m(2, 2) = this->data_[0][0]*this->data_[1][1] - this->data_[0][1]*this->data_[1][0];
337	>	if (zeroDiagElementIndex.size() == 2) {
338	>	int index = zeroDiagElementIndex[0];
339	>	m(index, index) = 1.0 / this->data_[index][index];
340	>	}else if (zeroDiagElementIndex.size() == 1) {
341
342	<	m /= det;
342	>	int a = (zeroDiagElementIndex[0] + 1) % 3;
343	>	int b = (zeroDiagElementIndex[0] + 2) %3;
344	>	RealType denom = this->data_[a][a] * this->data_[b][b] - this->data_[b][a]*this->data_[a][b];
345	>	m(a, a) = this->data_[b][b] /denom;
346	>	m(b, a) = -this->data_[b][a]/denom;
347	>
348	>	m(a,b) = -this->data_[a][b]/denom;
349	>	m(b, b) = this->data_[a][a]/denom;
350	>
351	>	}
352	>
353	>	/*
354	>	for(std::vector<int>::iterator iter = zeroDiagElementIndex.begin(); iter != zeroDiagElementIndex.end() ++iter) {
355	>	if (this->data_[*iter][0] > OpenMD::epsilon || this->data_[*iter][1] ||this->data_[*iter][2] ||
356	>	this->data_[0][*iter] > OpenMD::epsilon || this->data_[1][*iter] ||this->data_[2][*iter] ) {
357	>	std::cout << "can not inverse matrix" << std::endl;
358	>	}
359	>	}
360	>	*/
361	>	} else {
362	>
363	>	m(0, 0) = this->data_[1][1]*this->data_[2][2] - this->data_[1][2]*this->data_[2][1];
364	>	m(1, 0) = this->data_[1][2]*this->data_[2][0] - this->data_[1][0]*this->data_[2][2];
365	>	m(2, 0) = this->data_[1][0]*this->data_[2][1] - this->data_[1][1]*this->data_[2][0];
366	>	m(0, 1) = this->data_[2][1]*this->data_[0][2] - this->data_[2][2]*this->data_[0][1];
367	>	m(1, 1) = this->data_[2][2]*this->data_[0][0] - this->data_[2][0]*this->data_[0][2];
368	>	m(2, 1) = this->data_[2][0]*this->data_[0][1] - this->data_[2][1]*this->data_[0][0];
369	>	m(0, 2) = this->data_[0][1]*this->data_[1][2] - this->data_[0][2]*this->data_[1][1];
370	>	m(1, 2) = this->data_[0][2]*this->data_[1][0] - this->data_[0][0]*this->data_[1][2];
371	>	m(2, 2) = this->data_[0][0]*this->data_[1][1] - this->data_[0][1]*this->data_[1][0];
372	>
373	>	m /= det;
374	>	}
375		return m;
376		}
377	+
378	+	SquareMatrix3<Real> transpose() const{
379	+	SquareMatrix3<Real> result;
380	+
381	+	for (unsigned int i = 0; i < 3; i++)
382	+	for (unsigned int j = 0; j < 3; j++)
383	+	result(j, i) = this->data_[i][j];
384	+
385	+	return result;
386	+	}
387		/**
388		* Extract the eigenvalues and eigenvectors from a 3x3 matrix.
389		* The eigenvectors (the columns of V) will be normalized.
#	Line 353 | Line 420 | namespace oopse {
420		Vector3<Real> v_maxI, v_k, v_j;
421
422		// diagonalize using Jacobi
423	<	jacobi(a, w, v);
423	>	SquareMatrix3<Real>::jacobi(a, w, v);
424		// if all the eigenvalues are the same, return identity matrix
425		if (w[0] == w[1] && w[0] == w[2] ) {
426		v = SquareMatrix3<Real>::identity();
#	Line 499 | Line 566 | namespace oopse {
566		}
567
568
569	<	typedef SquareMatrix3<double> Mat3x3d;
570	<	typedef SquareMatrix3<double> RotMat3x3d;
569	>	typedef SquareMatrix3<RealType> Mat3x3d;
570	>	typedef SquareMatrix3<RealType> RotMat3x3d;
571
572	<	} //namespace oopse
572	>	} //namespace OpenMD
573		#endif // MATH_SQUAREMATRIX_HPP
574

Comparing trunk/src/math/SquareMatrix3.hpp (property svn:keywords):
Revision 633 by tim, Tue Sep 27 20:02:57 2005 UTC vs.
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC

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

Diff Legend

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