--- trunk/src/primitives/GhostBend.cpp 2004/09/24 16:27:58 3 +++ trunk/src/primitives/GhostBend.cpp 2009/11/25 20:02:06 1390 @@ -1,176 +1,101 @@ +/* + * 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] Vardeman & Gezelter, in progress (2009). + */ + +#include "primitives/GhostBend.hpp" +#include "primitives/DirectionalAtom.hpp" +namespace OpenMD { -#include -#include -#include - -#include "utils/simError.h" -#include "primitives/SRI.hpp" -#include "primitives/Atom.hpp" - - - -GhostBend::GhostBend( Atom &a, Atom &b ){ - - c_p_a = &a; - - if( !b.isDirectional() ){ + /**@todo still a lot left to improve*/ + void GhostBend::calcForce(RealType& angle) { + DirectionalAtom* ghostAtom = static_cast(atom2_); - // if atom b is not directional, then bad things will happen + Vector3d pos1 = atom1_->getPos(); + Vector3d pos2 = ghostAtom->getPos(); - sprintf( painCave.errMsg, - " Ghost Bend error: Atom # %d of type \"%s\" is not " - "directional.\n", - b.getIndex(), - b.getType() ); - painCave.isFatal = 1; - simError(); - } + Vector3d r12 = pos1 - pos2; + RealType d12 = r12.length(); + + RealType d12inv = 1.0 / d12; + + Vector3d r32 = ghostAtom->getElectroFrame().getColumn(2); + RealType d32 = r32.length(); + + RealType d32inv = 1.0 / d32; + + RealType cosTheta = dot(r12, r32) / (d12 * d32); + + //check roundoff + if (cosTheta > 1.0) { + cosTheta = 1.0; + } else if (cosTheta < -1.0) { + cosTheta = -1.0; + } + + RealType theta = acos(cosTheta); + + RealType firstDerivative; + + bendType_->calcForce(theta, potential_, firstDerivative); + + RealType sinTheta = sqrt(1.0 - cosTheta * cosTheta); + + if (fabs(sinTheta) < 1.0E-12) { + sinTheta = 1.0E-12; + } + + RealType commonFactor1 = -firstDerivative / sinTheta * d12inv; + RealType commonFactor2 = -firstDerivative / sinTheta * d32inv; + + Vector3d force1 = commonFactor1*(r12*(d12inv*cosTheta) - r32*d32inv); + Vector3d force3 = commonFactor2*(r32*(d32inv*cosTheta) - r12*d12inv); + atom1_->addFrc(force1); + ghostAtom->addFrc(-force1); + /**@todo test correctness */ + ghostAtom->addTrq(cross(r32, force3) ); - atomB = ( DirectionalAtom* ) &b; - - c_potential_E = 0.0; + atom1_->addParticlePot(potential_); + ghostAtom->addParticlePot(potential_); -} + angle = theta /M_PI * 180.0; + + } +} //end namespace OpenMD - -void GhostBend::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 u[3]; - - double aR[3], bR[3]; - double aF[3], bF[3], bTrq[3]; - - c_p_a->getPos( aR ); - atomB->getPos( bR ); - - - dx = aR[0] - bR[0]; - dy = aR[1] - bR[1]; - dz = aR[2] - bR[2]; - - atomB->getU(u); - - gx = u[0]; - gy = u[1]; - gz = u[2]; - - dx2 = dx * dx; - dy2 = dy * dy; - dz2 = dz * dz; - - gx2 = gx * gx; - gy2 = gy * gy; - gz2 = gz * gz; - - rij2 = dx2 + dy2 + dz2; - rkj2 = gx2 + gy2 + gz2; - - riji2 = 1.0 / rij2; - rkji2 = 1.0 / rkj2; - - dot = dx * gx + dy * gy + dz * gz; - denom = sqrt((riji2 * rkji2)); - cosang = dot * denom; - - if(cosang > 1.0)cosang = 1.0; - if(cosang < -1.0) cosang = -1.0; - - angl = acos(cosang); - angl = angl * 180.0 / M_PI; - - 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 + daxk; - bF[1] = dayj + dayk; - bF[2] = dazj + dazk; - - bTrq[0] = gy*dazk - gz*dayk; - bTrq[1] = gz*daxk - gx*dazk; - bTrq[2] = gx*dayk - gy*daxk; - - - c_p_a->addFrc( aF ); - atomB->addFrc( bF ); - atomB->addTrq( bTrq ); - - return; -} - -void GhostBend::setConstants( double the_c1, double the_c2, double the_c3, - double the_Th0 ){ - c1 = the_c1; - c2 = the_c2; - c3 = the_c3; - theta0 = the_Th0; -} - - -double GhostBend::bend_force( double theta ){ - - double dt, dt2; - double force; - - dt = ( theta - theta0 ) * M_PI / 180.0; - dt2 = dt * dt; - - c_potential_E = ( c1 * dt2 ) + ( c2 * dt ) + c3; - force = -( ( 2.0 * c1 * dt ) + c2 ); - return force; -}