--- trunk/src/math/SquareMatrix3.hpp 2006/05/17 21:51:42 963 +++ trunk/src/math/SquareMatrix3.hpp 2014/05/31 22:35:05 2000 @@ -6,19 +6,10 @@ * redistribute this software in source and binary code form, provided * that the following conditions are met: * - * 1. Acknowledgement of the program authors must be made in any - * publication of scientific results based in part on use of the - * program. An acceptable form of acknowledgement is citation of - * the article in which the program was described (Matthew - * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher - * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented - * Parallel Simulation Engine for Molecular Dynamics," - * J. Comput. Chem. 26, pp. 252-271 (2005)) - * - * 2. Redistributions of source code must retain the above copyright + * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * - * 3. Redistributions in binary form must reproduce the above copyright + * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the * distribution. @@ -37,6 +28,16 @@ * arising out of the use of or inability to use software, even if the * University of Notre Dame has been advised of the possibility of * such damages. + * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * research, please cite the appropriate papers when you publish your + * work. Good starting points are: + * + * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). + * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ /** @@ -47,12 +48,14 @@ */ #ifndef MATH_SQUAREMATRIX3_HPP #define MATH_SQUAREMATRIX3_HPP +#include "config.h" +#include #include #include "Quaternion.hpp" #include "SquareMatrix.hpp" #include "Vector3.hpp" #include "utils/NumericConstant.hpp" -namespace oopse { +namespace OpenMD { template class SquareMatrix3 : public SquareMatrix { @@ -109,6 +112,7 @@ namespace oopse { return *this; } + /** * Sets this matrix to a rotation matrix by three euler angles * @ param euler @@ -121,7 +125,7 @@ namespace oopse { * Sets this matrix to a rotation matrix by three euler angles * @param phi * @param theta - * @psi theta + * @param psi */ void setupRotMat(Real phi, Real theta, Real psi) { Real sphi, stheta, spsi; @@ -187,7 +191,6 @@ namespace oopse { } - /** * Returns the quaternion from this rotation matrix * @return the quaternion from this rotation matrix @@ -244,10 +247,10 @@ namespace oopse { * @return the euler angles in a vector * @exception invalid rotation matrix * We use so-called "x-convention", which is the most common definition. - * In this convention, the rotation given by Euler angles (phi, theta, psi), where the first - * rotation is by an angle phi about the z-axis, the second is by an angle - * theta (0 <= theta <= 180)about the x-axis, and thethird is by an angle psi about the - * z-axis (again). + * In this convention, the rotation given by Euler angles (phi, theta, + * psi), where the first rotation is by an angle phi about the z-axis, + * the second is by an angle theta (0 <= theta <= 180) about the x-axis, + * and the third is by an angle psi about the z-axis (again). */ Vector3 toEulerAngles() { Vector3 myEuler; @@ -263,15 +266,17 @@ namespace oopse { ctheta = this->data_[2][2]; stheta = sqrt(1.0 - ctheta * ctheta); - // when sin(theta) is close to 0, we need to consider singularity - // In this case, we can assign an arbitary value to phi (or psi), and then determine - // the psi (or phi) or vice-versa. We'll assume that phi always gets the rotation, and psi is 0 - // in cases of singularity. + // when sin(theta) is close to 0, we need to consider + // singularity In this case, we can assign an arbitary value to + // phi (or psi), and then determine the psi (or phi) or + // vice-versa. We'll assume that phi always gets the rotation, + // and psi is 0 in cases of singularity. // we use atan2 instead of atan, since atan2 will give us -Pi to Pi. - // Since 0 <= theta <= 180, sin(theta) will be always non-negative. Therefore, it never - // change the sign of both of the parameters passed to atan2. + // Since 0 <= theta <= 180, sin(theta) will be always + // non-negative. Therefore, it will never change the sign of both of + // the parameters passed to atan2. - if (fabs(stheta) <= oopse::epsilon){ + if (fabs(stheta) < 1e-6){ psi = 0.0; phi = atan2(-this->data_[1][0], this->data_[0][0]); } @@ -283,10 +288,10 @@ namespace oopse { //wrap phi and psi, make sure they are in the range from 0 to 2*Pi if (phi < 0) - phi += M_PI; + phi += 2.0 * M_PI; if (psi < 0) - psi += M_PI; + psi += 2.0 * M_PI; myEuler[0] = phi; myEuler[1] = theta; @@ -319,12 +324,12 @@ namespace oopse { SquareMatrix3 inverse() const { SquareMatrix3 m; RealType det = determinant(); - if (fabs(det) <= oopse::epsilon) { + if (fabs(det) <= OpenMD::epsilon) { //"The method was called on a matrix with |determinant| <= 1e-6.", //"This is a runtime or a programming error in your application."); std::vector zeroDiagElementIndex; for (int i =0; i < 3; ++i) { - if (fabs(this->data_[i][i]) <= oopse::epsilon) { + if (fabs(this->data_[i][i]) <= OpenMD::epsilon) { zeroDiagElementIndex.push_back(i); } } @@ -347,8 +352,8 @@ namespace oopse { /* for(std::vector::iterator iter = zeroDiagElementIndex.begin(); iter != zeroDiagElementIndex.end() ++iter) { - if (this->data_[*iter][0] > oopse::epsilon || this->data_[*iter][1] ||this->data_[*iter][2] || - this->data_[0][*iter] > oopse::epsilon || this->data_[1][*iter] ||this->data_[2][*iter] ) { + if (this->data_[*iter][0] > OpenMD::epsilon || this->data_[*iter][1] ||this->data_[*iter][2] || + this->data_[0][*iter] > OpenMD::epsilon || this->data_[1][*iter] ||this->data_[2][*iter] ) { std::cout << "can not inverse matrix" << std::endl; } } @@ -415,7 +420,7 @@ namespace oopse { Vector3 v_maxI, v_k, v_j; // diagonalize using Jacobi - jacobi(a, w, v); + SquareMatrix3::jacobi(a, w, v); // if all the eigenvalues are the same, return identity matrix if (w[0] == w[1] && w[0] == w[2] ) { v = SquareMatrix3::identity(); @@ -564,6 +569,9 @@ namespace oopse { typedef SquareMatrix3 Mat3x3d; typedef SquareMatrix3 RotMat3x3d; -} //namespace oopse + const Mat3x3d M3Zero(0.0); + + +} //namespace OpenMD #endif // MATH_SQUAREMATRIX_HPP