--- trunk/src/primitives/Bend.cpp 2004/09/24 16:27:58 3 +++ branches/development/src/primitives/Bend.cpp 2011/11/22 20:38:56 1665 @@ -1,134 +1,105 @@ -#include "primitives/SRI.hpp" -#include "primitives/Atom.hpp" - -#include -#include -#include - -void Bend::set_atoms( Atom &a, Atom &b, Atom &c){ - - c_p_a = &a; - c_p_b = &b; - c_p_c = &c; - - c_potential_E = 0.0; -} - - -void Bend::calc_forces(){ - - double dx,dy,dz,gx,gy,gz,dx2,dy2,dz2,gx2,gy2,gz2; - double rij2, rkj2, riji2, rkji2, dot, denom, cosang, angl; - - double sina2, sinai; - - double comf2, comf3, comf4; - double dcsidx, dcsidy, dcsidz, dcskdx, dcskdy, dcskdz; - // double dcsjdx, dcsjdy, dcsjdz; - double dadxi, dadyi, dadzi; - double dadxk, dadyk, dadzk;//, dadxj, dadyj, dadzj; - double daxi, dayi, dazi, daxk, dayk, dazk, daxj, dayj, dazj; - - double aR[3], bR[3], cR[3]; - double aF[3], bF[3], cF[3]; - - c_p_a->getPos( aR ); - c_p_b->getPos( bR ); - c_p_c->getPos( cR ); - - - dx = aR[0] - bR[0]; - dy = aR[1] - bR[1]; - dz = aR[2] - bR[2]; +/* + * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. + * + * The University of Notre Dame grants you ("Licensee") a + * non-exclusive, royalty free, license to use, modify and + * redistribute this software in source and binary code form, provided + * that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * + * 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. + * + * This software is provided "AS IS," without a warranty of any + * kind. All express or implied conditions, representations and + * warranties, including any implied warranty of merchantability, + * fitness for a particular purpose or non-infringement, are hereby + * excluded. The University of Notre Dame and its licensors shall not + * be liable for any damages suffered by licensee as a result of + * using, modifying or distributing the software or its + * derivatives. In no event will the University of Notre Dame or its + * licensors be liable for any lost revenue, profit or data, or for + * direct, indirect, special, consequential, incidental or punitive + * damages, however caused and regardless of the theory of liability, + * 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] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). + */ - gx = cR[0] - bR[0]; - gy = cR[1] - bR[1]; - gz = cR[2] - bR[2]; - - dx2 = dx * dx; - dy2 = dy * dy; - dz2 = dz * dz; +#include "primitives/Bend.hpp" - gx2 = gx * gx; - gy2 = gy * gy; - gz2 = gz * gz; +namespace OpenMD { - rij2 = dx2 + dy2 + dz2; - rkj2 = gx2 + gy2 + gz2; - - riji2 = 1.0 / rij2; - rkji2 = 1.0 / rkj2; + /**@todo still a lot left to improve*/ + void Bend::calcForce(RealType& angle) { + Vector3d pos1 = atom1_->getPos(); + Vector3d pos2 = atom2_->getPos(); + Vector3d pos3 = atom3_->getPos(); + + Vector3d r21 = pos1 - pos2; + RealType d21 = r21.length(); + + RealType d21inv = 1.0 / d21; + + Vector3d r23 = pos3 - pos2; + RealType d23 = r23.length(); + + RealType d23inv = 1.0 / d23; + + RealType cosTheta = dot(r21, r23) / (d21 * d23); + + //check roundoff + if (cosTheta > 1.0) { + cosTheta = 1.0; + } else if (cosTheta < -1.0) { + cosTheta = -1.0; + } + + RealType theta = acos(cosTheta); + + RealType dVdTheta; - dot = dx * gx + dy * gy + dz * gz; - denom = sqrt((riji2 * rkji2)); - cosang = dot * denom; + bendType_->calcForce(theta, potential_, dVdTheta); + + RealType sinTheta = sqrt(1.0 - cosTheta * cosTheta); + + if (fabs(sinTheta) < 1.0E-6) { + sinTheta = 1.0E-6; + } + + RealType commonFactor1 = dVdTheta / sinTheta * d21inv; + RealType commonFactor2 = dVdTheta / sinTheta * d23inv; + + Vector3d force1 = commonFactor1 * (r23 * d23inv - r21*d21inv*cosTheta); + Vector3d force3 = commonFactor2 * (r21 * d21inv - r23*d23inv*cosTheta); - if(cosang > 1.0)cosang = 1.0; - if(cosang < -1.0) cosang = -1.0; + // Total force in current bend is zero + Vector3d force2 = force1 + force3; + force2 *= -1.0; + + atom1_->addFrc(force1); + atom2_->addFrc(force2); + atom3_->addFrc(force3); - angl = acos(cosang); - angl = angl * 180.0 / M_PI; + atom1_->addParticlePot(potential_); + atom2_->addParticlePot(potential_); + atom3_->addParticlePot(potential_); + + angle = theta /M_PI * 180.0; + } - sina2 = 1.0 - cosang*cosang; - if(fabs(sina2) < 1.0E-12 ) sina2 = 1.0E-12; - sinai = 1.0 / sqrt(sina2); - - comf2 = cosang * riji2; - comf3 = cosang * rkji2; - comf4 = bend_force(angl); - - - dcsidx = gx*denom - comf2*dx; - dcsidy = gy*denom - comf2*dy; - dcsidz = gz*denom - comf2*dz; - - dcskdx = dx*denom - comf3*gx; - dcskdy = dy*denom - comf3*gy; - dcskdz = dz*denom - comf3*gz; - -// dcsjdx = -dcsidx - dcskdx; -// dcsjdy = -dcsidy - dcskdy; -// dcsjdz = -dcsidz - dcskdz; - - dadxi = -sinai*dcsidx; - dadyi = -sinai*dcsidy; - dadzi = -sinai*dcsidz; - - dadxk = -sinai*dcskdx; - dadyk = -sinai*dcskdy; - dadzk = -sinai*dcskdz; - -// dadxj = -dadxi - dadxk; -// dadyj = -dadyi - dadyk; -// dadzj = -dadzi - dadzk; - - daxi = comf4*dadxi; - dayi = comf4*dadyi; - dazi = comf4*dadzi; - - daxk = comf4*dadxk; - dayk = comf4*dadyk; - dazk = comf4*dadzk; - - daxj = -daxi - daxk; - dayj = -dayi - dayk; - dazj = -dazi - dazk; - - aF[0] = daxi; - aF[1] = dayi; - aF[2] = dazi; - - bF[0] = daxj; - bF[1] = dayj; - bF[2] = dazj; - - cF[0] = daxk; - cF[1] = dayk; - cF[2] = dazk; - - c_p_a->addFrc(aF); - c_p_b->addFrc(bF); - c_p_c->addFrc(cF); - - return; -} +} //end namespace OpenMD