--- trunk/src/primitives/Torsion.cpp 2005/01/12 22:41:40 246 +++ trunk/src/primitives/Torsion.cpp 2010/06/17 14:48:02 1446 @@ -1,4 +1,4 @@ - /* +/* * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. * * The University of Notre Dame grants you ("Licensee") a @@ -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,17 +28,27 @@ * 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, 24107 (2008). + * [4] Vardeman & Gezelter, in progress (2009). */ #include "primitives/Torsion.hpp" -namespace oopse { +namespace OpenMD { -Torsion::Torsion(Atom *atom1, Atom *atom2, Atom *atom3, Atom *atom4, - TorsionType *tt) : + Torsion::Torsion(Atom *atom1, Atom *atom2, Atom *atom3, Atom *atom4, + TorsionType *tt) : atom1_(atom1), atom2_(atom2), atom3_(atom3), atom4_(atom4), torsionType_(tt) { } -void Torsion::calcForce() { + void Torsion::calcForce(RealType& angle) { + Vector3d pos1 = atom1_->getPos(); Vector3d pos2 = atom2_->getPos(); Vector3d pos3 = atom3_->getPos(); @@ -59,75 +60,49 @@ void Torsion::calcForce() { // Calculate the cross products and distances Vector3d A = cross(r21, r32); - double rA = A.length(); + RealType rA = A.length(); Vector3d B = cross(r32, r43); - double rB = B.length(); - Vector3d C = cross(r32, A); - double rC = C.length(); + RealType rB = B.length(); + /* + If either of the two cross product vectors is tiny, that means + the three atoms involved are colinear, and the torsion angle is + going to be undefined. The easiest check for this problem is + to use the product of the two lengths. + */ + if (rA * rB < OpenMD::epsilon) return; + A.normalize(); - B.normalize(); - C.normalize(); + B.normalize(); // Calculate the sin and cos - double cos_phi = dot(A, B) ; - double sin_phi = dot(C, B); - - double dVdPhi; - torsionType_->calcForce(cos_phi, sin_phi, potential_, dVdPhi); - - Vector3d f1; - Vector3d f2; - Vector3d f3; - - // Next, we want to calculate the forces. In order - // to do that, we first need to figure out whether the - // sin or cos form will be more stable. For this, - // just look at the value of phi - //if (fabs(sin_phi) > 0.1) { - // use the sin version to avoid 1/cos terms - - Vector3d dcosdA = (cos_phi * A - B) /rA; - Vector3d dcosdB = (cos_phi * B - A) /rB; - - double dVdcosPhi = -dVdPhi / sin_phi; - - f1 = dVdcosPhi * cross(r32, dcosdA); - f2 = dVdcosPhi * ( cross(r43, dcosdB) - cross(r21, dcosdA)); - f3 = dVdcosPhi * cross(dcosdB, r32); - - /** @todo fix below block, must be something wrong with the sign somewhere */ - //} else { - // This angle is closer to 0 or 180 than it is to - // 90, so use the cos version to avoid 1/sin terms - - //double dVdsinPhi = dVdPhi /cos_phi; - //Vector3d dsindB = (sin_phi * B - C) /rB; - //Vector3d dsindC = (sin_phi * C - B) /rC; - - //f1.x() = dVdsinPhi*((r32.y()*r32.y() + r32.z()*r32.z())*dsindC.x() - r32.x()*r32.y()*dsindC.y() - r32.x()*r32.z()*dsindC.z()); - - //f1.y() = dVdsinPhi*((r32.z()*r32.z() + r32.x()*r32.x())*dsindC.y() - r32.y()*r32.z()*dsindC.z() - r32.y()*r32.x()*dsindC.x()); - - //f1.z() = dVdsinPhi*((r32.x()*r32.x() + r32.y()*r32.y())*dsindC.z() - r32.z()*r32.x()*dsindC.x() - r32.z()*r32.y()*dsindC.y()); - - //f2.x() = dVdsinPhi*(-(r32.y()*r21.y() + r32.z()*r21.z())*dsindC.x() + (2.0*r32.x()*r21.y() - r21.x()*r32.y())*dsindC.y() - //+ (2.0*r32.x()*r21.z() - r21.x()*r32.z())*dsindC.z() + dsindB.z()*r43.y() - dsindB.y()*r43.z()); - - //f2.y() = dVdsinPhi*(-(r32.z()*r21.z() + r32.x()*r21.x())*dsindC.y() + (2.0*r32.y()*r21.z() - r21.y()*r32.z())*dsindC.z() - //+ (2.0*r32.y()*r21.x() - r21.y()*r32.x())*dsindC.x() + dsindB.x()*r43.z() - dsindB.z()*r43.x()); - - //f2.z() = dVdsinPhi*(-(r32.x()*r21.x() + r32.y()*r21.y())*dsindC.z() + (2.0*r32.z()*r21.x() - r21.z()*r32.x())*dsindC.x() - //+(2.0*r32.z()*r21.y() - r21.z()*r32.y())*dsindC.y() + dsindB.y()*r43.x() - dsindB.x()*r43.y()); - - //f3 = dVdsinPhi * cross(r32, dsindB); - - //} - + RealType cos_phi = dot(A, B) ; + if (cos_phi > 1.0) cos_phi = 1.0; + if (cos_phi < -1.0) cos_phi = -1.0; + + RealType dVdcosPhi; + torsionType_->calcForce(cos_phi, potential_, dVdcosPhi); + Vector3d f1 ; + Vector3d f2 ; + Vector3d f3 ; + + Vector3d dcosdA = (cos_phi * A - B) /rA; + Vector3d dcosdB = (cos_phi * B - A) /rB; + + f1 = dVdcosPhi * cross(r32, dcosdA); + f2 = dVdcosPhi * ( cross(r43, dcosdB) - cross(r21, dcosdA)); + f3 = dVdcosPhi * cross(dcosdB, r32); + atom1_->addFrc(f1); atom2_->addFrc(f2 - f1); atom3_->addFrc(f3 - f2); atom4_->addFrc(-f3); + + atom1_->addParticlePot(potential_); + atom2_->addParticlePot(potential_); + atom3_->addParticlePot(potential_); + atom4_->addParticlePot(potential_); + + angle = acos(cos_phi) /M_PI * 180.0; + } } - -}