--- trunk/src/math/SquareMatrix3.hpp 2005/09/27 20:02:57 633 +++ 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; @@ -168,6 +172,25 @@ namespace oopse { *this = q.toRotationMatrix3(); } + void setupSkewMat(Vector3 v) { + setupSkewMat(v[0], v[1], v[2]); + } + + void setupSkewMat(Real v1, Real v2, Real v3) { + this->data_[0][0] = 0; + this->data_[0][1] = -v3; + this->data_[0][2] = v2; + this->data_[1][0] = v3; + this->data_[1][1] = 0; + this->data_[1][2] = -v1; + this->data_[2][0] = -v2; + this->data_[2][1] = v1; + this->data_[2][2] = 0; + + + } + + /** * Returns the quaternion from this rotation matrix * @return the quaternion from this rotation matrix @@ -179,7 +202,7 @@ namespace oopse { Real ad1, ad2, ad3; t = this->data_[0][0] + this->data_[1][1] + this->data_[2][2] + 1.0; - if( t > 0.0 ){ + if( t > NumericConstant::epsilon ){ s = 0.5 / sqrt( t ); q[0] = 0.25 / s; @@ -224,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; @@ -239,19 +262,21 @@ namespace oopse { // set the tolerance for Euler angles and rotation elements - theta = acos(std::min(1.0, std::max(-1.0,this->data_[2][2]))); + theta = acos(std::min((RealType)1.0, std::max((RealType)-1.0,this->data_[2][2]))); 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]); } @@ -263,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; @@ -298,25 +323,67 @@ namespace oopse { */ SquareMatrix3 inverse() const { SquareMatrix3 m; - double det = determinant(); - if (fabs(det) <= oopse::epsilon) { + RealType det = determinant(); + 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]) <= OpenMD::epsilon) { + zeroDiagElementIndex.push_back(i); + } + } - m(0, 0) = this->data_[1][1]*this->data_[2][2] - this->data_[1][2]*this->data_[2][1]; - m(1, 0) = this->data_[1][2]*this->data_[2][0] - this->data_[1][0]*this->data_[2][2]; - m(2, 0) = this->data_[1][0]*this->data_[2][1] - this->data_[1][1]*this->data_[2][0]; - m(0, 1) = this->data_[2][1]*this->data_[0][2] - this->data_[2][2]*this->data_[0][1]; - m(1, 1) = this->data_[2][2]*this->data_[0][0] - this->data_[2][0]*this->data_[0][2]; - m(2, 1) = this->data_[2][0]*this->data_[0][1] - this->data_[2][1]*this->data_[0][0]; - m(0, 2) = this->data_[0][1]*this->data_[1][2] - this->data_[0][2]*this->data_[1][1]; - m(1, 2) = this->data_[0][2]*this->data_[1][0] - this->data_[0][0]*this->data_[1][2]; - m(2, 2) = this->data_[0][0]*this->data_[1][1] - this->data_[0][1]*this->data_[1][0]; + if (zeroDiagElementIndex.size() == 2) { + int index = zeroDiagElementIndex[0]; + m(index, index) = 1.0 / this->data_[index][index]; + }else if (zeroDiagElementIndex.size() == 1) { - m /= det; + int a = (zeroDiagElementIndex[0] + 1) % 3; + int b = (zeroDiagElementIndex[0] + 2) %3; + RealType denom = this->data_[a][a] * this->data_[b][b] - this->data_[b][a]*this->data_[a][b]; + m(a, a) = this->data_[b][b] /denom; + m(b, a) = -this->data_[b][a]/denom; + + m(a,b) = -this->data_[a][b]/denom; + m(b, b) = this->data_[a][a]/denom; + + } + +/* + for(std::vector::iterator iter = zeroDiagElementIndex.begin(); iter != zeroDiagElementIndex.end() ++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; + } + } +*/ + } else { + + m(0, 0) = this->data_[1][1]*this->data_[2][2] - this->data_[1][2]*this->data_[2][1]; + m(1, 0) = this->data_[1][2]*this->data_[2][0] - this->data_[1][0]*this->data_[2][2]; + m(2, 0) = this->data_[1][0]*this->data_[2][1] - this->data_[1][1]*this->data_[2][0]; + m(0, 1) = this->data_[2][1]*this->data_[0][2] - this->data_[2][2]*this->data_[0][1]; + m(1, 1) = this->data_[2][2]*this->data_[0][0] - this->data_[2][0]*this->data_[0][2]; + m(2, 1) = this->data_[2][0]*this->data_[0][1] - this->data_[2][1]*this->data_[0][0]; + m(0, 2) = this->data_[0][1]*this->data_[1][2] - this->data_[0][2]*this->data_[1][1]; + m(1, 2) = this->data_[0][2]*this->data_[1][0] - this->data_[0][0]*this->data_[1][2]; + m(2, 2) = this->data_[0][0]*this->data_[1][1] - this->data_[0][1]*this->data_[1][0]; + + m /= det; + } return m; } + + SquareMatrix3 transpose() const{ + SquareMatrix3 result; + + for (unsigned int i = 0; i < 3; i++) + for (unsigned int j = 0; j < 3; j++) + result(j, i) = this->data_[i][j]; + + return result; + } /** * Extract the eigenvalues and eigenvectors from a 3x3 matrix. * The eigenvectors (the columns of V) will be normalized. @@ -353,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(); @@ -499,9 +566,12 @@ namespace oopse { } - typedef SquareMatrix3 Mat3x3d; - typedef SquareMatrix3 RotMat3x3d; + typedef SquareMatrix3 Mat3x3d; + typedef SquareMatrix3 RotMat3x3d; -} //namespace oopse + const Mat3x3d M3Zero(0.0); + + +} //namespace OpenMD #endif // MATH_SQUAREMATRIX_HPP