ReferenceAmoebaStretchBendForce

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernelFactory.cpp

@@ -44,8 +44,8 @@ fprintf( stderr,"In registerKernelFactories AmoebaReferenceKernelFactory\n" ); f
              platform.registerKernelFactory(CalcAmoebaHarmonicInPlaneAngleForceKernel::Name(), factory);
              platform.registerKernelFactory(CalcAmoebaTorsionForceKernel::Name(), factory);
              platform.registerKernelFactory(CalcAmoebaPiTorsionForceKernel::Name(), factory);
-/*
              platform.registerKernelFactory(CalcAmoebaStretchBendForceKernel::Name(), factory);
+/*
              platform.registerKernelFactory(CalcAmoebaOutOfPlaneBendForceKernel::Name(), factory);
              platform.registerKernelFactory(CalcAmoebaTorsionTorsionForceKernel::Name(), factory);
              platform.registerKernelFactory(CalcAmoebaMultipoleForceKernel::Name(), factory);
@@ -77,10 +77,10 @@ KernelImpl* AmoebaReferenceKernelFactory::createKernelImpl(std::string name, con
     if (name == CalcAmoebaPiTorsionForceKernel::Name())
         return new ReferenceCalcAmoebaPiTorsionForceKernel(name, platform, context.getSystem());
 
-/*
     if (name == CalcAmoebaStretchBendForceKernel::Name())
         return new ReferenceCalcAmoebaStretchBendForceKernel(name, platform, context.getSystem());
 
+/*
     if (name == CalcAmoebaOutOfPlaneBendForceKernel::Name())
         return new ReferenceCalcAmoebaOutOfPlaneBendForceKernel(name, platform, context.getSystem());
 

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp

@@ -30,6 +30,7 @@
 #include "AmoebaReferenceHarmonicInPlaneAngleForce.h"
 #include "AmoebaReferenceTorsionForce.h"
 #include "AmoebaReferencePiTorsionForce.h"
+#include "AmoebaReferenceStretchBendForce.h"
 #include "ReferencePlatform.h"
 #include "openmm/internal/ContextImpl.h"
 //#include "internal/AmoebaMultipoleForceImpl.h"
@@ -326,49 +327,53 @@ double ReferenceCalcAmoebaPiTorsionForceKernel::execute(ContextImpl& context, bo
     return static_cast<double>(energy);
 }
 
-//
-//ReferenceCalcAmoebaStretchBendForceKernel::ReferenceCalcAmoebaStretchBendForceKernel(std::string name, const Platform& platform, System& system) :
-//                   CalcAmoebaStretchBendForceKernel(name, platform), system(system) {
-//    data.incrementKernelCount();
-//}
-//
-//ReferenceCalcAmoebaStretchBendForceKernel::~ReferenceCalcAmoebaStretchBendForceKernel() {
-//    data.decrementKernelCount();
-//}
-//
-//void ReferenceCalcAmoebaStretchBendForceKernel::initialize(const System& system, const AmoebaStretchBendForce& force) {
-//
-//    data.setAmoebaLocalForcesKernel( this );
-//    numStretchBends                     = force.getNumStretchBends();
-//
-//    std::vector<int>   particle1(numStretchBends);
-//    std::vector<int>   particle2(numStretchBends);
-//    std::vector<int>   particle3(numStretchBends);
-//    std::vector<RealOpenMM> lengthABParameters(numStretchBends);
-//    std::vector<RealOpenMM> lengthCBParameters(numStretchBends);
-//    std::vector<RealOpenMM> angleParameters(numStretchBends);
-//    std::vector<RealOpenMM> kParameters(numStretchBends);
-//
-//    for (int i = 0; i < numStretchBends; i++) {
-//
-//        double lengthAB, lengthCB, angle, k;
-//
-//        force.getStretchBendParameters(i, particle1[i], particle2[i], particle3[i], lengthAB, lengthCB, angle, k);
-//        lengthABParameters[i] = static_cast<RealOpenMM>(lengthAB);
-//        lengthCBParameters[i] = static_cast<RealOpenMM>(lengthCB);
-//        angleParameters[i]    = static_cast<RealOpenMM>(angle);
-//        kParameters[i]        = static_cast<RealOpenMM>(k);
-//    }
-//    gpuSetAmoebaStretchBendParameters(data.getAmoebaGpu(), particle1, particle2, particle3, lengthABParameters, lengthCBParameters, angleParameters, kParameters);
-//
-//}
-//
-//double ReferenceCalcAmoebaStretchBendForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-//    if( data.getAmoebaLocalForcesKernel() == this ){
-//        computeAmoebaLocalForces( data );
-//    }
-//    return 0.0;
-//}
+ReferenceCalcAmoebaStretchBendForceKernel::ReferenceCalcAmoebaStretchBendForceKernel(std::string name, const Platform& platform, System& system) :
+                   CalcAmoebaStretchBendForceKernel(name, platform), system(system) {
+}
+
+ReferenceCalcAmoebaStretchBendForceKernel::~ReferenceCalcAmoebaStretchBendForceKernel() {
+}
+
+void ReferenceCalcAmoebaStretchBendForceKernel::initialize(const System& system, const AmoebaStretchBendForce& force) {
+
+    numStretchBends = force.getNumStretchBends();
+    for ( int ii = 0; ii < numStretchBends; ii++) {
+        int particle1Index, particle2Index, particle3Index;
+        double lengthAB, lengthCB, angle, k;
+        force.getStretchBendParameters(ii, particle1Index, particle2Index, particle3Index, lengthAB, lengthCB, angle, k);
+        particle1.push_back( particle1Index ); 
+        particle2.push_back( particle2Index ); 
+        particle3.push_back( particle3Index ); 
+        lengthABParameters.push_back( static_cast<RealOpenMM>(lengthAB) );
+        lengthCBParameters.push_back( static_cast<RealOpenMM>(lengthCB) );
+        angleParameters.push_back(    static_cast<RealOpenMM>(angle) );
+        kParameters.push_back(        static_cast<RealOpenMM>(k) );
+    }
+}
+
+double ReferenceCalcAmoebaStretchBendForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    RealOpenMM** posData   = extractPositions(context);
+    RealOpenMM** forceData = extractForces(context);
+    RealOpenMM energy      = 0.0; 
+    for( unsigned int ii = 0; ii < numStretchBends; ii++ ){
+
+        int particle1Index      = particle1[ii];
+        int particle2Index      = particle2[ii];
+        int particle3Index      = particle3[ii];
+
+        RealOpenMM* forces[3];
+        forces[0]               = forceData[particle1Index];
+        forces[1]               = forceData[particle2Index];
+        forces[2]               = forceData[particle3Index];
+
+        energy                 += AmoebaReferenceStretchBendForce::calculateForceAndEnergy( 
+                                       posData[particle1Index], posData[particle2Index], posData[particle3Index],
+                                       lengthABParameters[ii], lengthCBParameters[ii],
+                                       angleParameters[ii], kParameters[ii], forces );
+    }
+    return static_cast<double>(energy);
+}
+
 //ReferenceCalcAmoebaOutOfPlaneBendForceKernel::ReferenceCalcAmoebaOutOfPlaneBendForceKernel(std::string name, const Platform& platform, System& system) :
 //          CalcAmoebaOutOfPlaneBendForceKernel(name, platform), system(system) {
 //    data.incrementKernelCount();

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h

@@ -214,33 +214,40 @@ private:
     System& system;
 };
 
-// /**
-//  * This kernel is invoked by AmoebaStretchBendForce to calculate the forces acting on the system and the energy of the system.
-//  */
-// class ReferenceCalcAmoebaStretchBendForceKernel : public CalcAmoebaStretchBendForceKernel {
-// public:
-//     ReferenceCalcAmoebaStretchBendForceKernel(std::string name, const Platform& platform, System& system);
-//     ~ReferenceCalcAmoebaStretchBendForceKernel();
-//     /**
-//      * Initialize the kernel.
-//      * 
-//      * @param system     the System this kernel will be applied to
-//      * @param force      the AmoebaStretchBendForce this kernel will be used for
-//      */
-//     void initialize(const System& system, const AmoebaStretchBendForce& force);
-//     /**
-//      * Execute the kernel to calculate the forces and/or energy.
-//      *
-//      * @param context        the context in which to execute this kernel
-//      * @param includeForces  true if forces should be calculated
-//      * @param includeEnergy  true if the energy should be calculated
-//      * @return the potential energy due to the force
-//      */
-//     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
-// private:
-//     int numStretchBends;
-//     System& system;
-// };
+/**
+ * This kernel is invoked by AmoebaStretchBendForce to calculate the forces acting on the system and the energy of the system.
+ */
+class ReferenceCalcAmoebaStretchBendForceKernel : public CalcAmoebaStretchBendForceKernel {
+public:
+    ReferenceCalcAmoebaStretchBendForceKernel(std::string name, const Platform& platform, System& system);
+    ~ReferenceCalcAmoebaStretchBendForceKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param force      the AmoebaStretchBendForce this kernel will be used for
+     */
+    void initialize(const System& system, const AmoebaStretchBendForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+private:
+    int numStretchBends;
+    std::vector<int>   particle1;
+    std::vector<int>   particle2;
+    std::vector<int>   particle3;
+    std::vector<RealOpenMM> lengthABParameters;
+    std::vector<RealOpenMM> lengthCBParameters;
+    std::vector<RealOpenMM> angleParameters;
+    std::vector<RealOpenMM> kParameters;
+    System& system;
+};
 // 
 // /**
 //  * This kernel is invoked by AmoebaOutOfPlaneBendForce to calculate the forces acting on the system and the energy of the system.

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceStretchBendForce.cpp

+
+/* Portions copyright (c) 2006 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#include "AmoebaReferenceForce.h"
+#include "AmoebaReferenceStretchBendForce.h"
+#include <vector>
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate Amoeba harmonic angle ixn (force and energy)
+
+   @param positionAtomA           Cartesian coordinates of atom A
+   @param positionAtomB           Cartesian coordinates of atom B
+   @param positionAtomC           Cartesian coordinates of atom C
+   @param positionAtomD           Cartesian coordinates of atom D
+   @param angleLength             angle
+   @param angleK                  quadratic angle force
+   @param angleCubic              cubic angle force parameter
+   @param angleQuartic            quartic angle force parameter
+   @param anglePentic             pentic angle force parameter
+   @param angleSextic             sextic angle force parameter
+   @param forces                  force vector
+
+   @return energy
+
+   --------------------------------------------------------------------------------------- */
+
+RealOpenMM AmoebaReferenceStretchBendForce::calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
+                                                                     const RealOpenMM* positionAtomC,
+                                                                     RealOpenMM lengthAB,      RealOpenMM lengthCB,
+                                                                     RealOpenMM idealAngle,    RealOpenMM kParameter,
+                                                                     RealOpenMM** forces ){
+
+   // ---------------------------------------------------------------------------------------
+
+    //static const std::string methodName = "AmoebaReferenceStretchBendForce::calculateHarmonicForce";
+ 
+    static const RealOpenMM zero          = 0.0;
+    static const RealOpenMM one           = 1.0;
+
+   // ---------------------------------------------------------------------------------------
+
+   enum { A, B, C, LastAtomIndex };
+   enum { AB, CB, CBxAB, ABxP, CBxP, LastDeltaIndex };
+
+   // ---------------------------------------------------------------------------------------
+
+   // get deltaR between various combinations of the 3 atoms
+   // and various intermediate terms
+
+    std::vector<RealOpenMM> deltaR[LastDeltaIndex];
+    for( unsigned int ii = 0; ii < LastDeltaIndex; ii++ ){
+        deltaR[ii].resize(3);
+    }
+    AmoebaReferenceForce::loadDeltaR( positionAtomB, positionAtomA, deltaR[AB] );
+    AmoebaReferenceForce::loadDeltaR( positionAtomB, positionAtomC, deltaR[CB] );
+    RealOpenMM  rAB2 = AmoebaReferenceForce::getNormSquared3( deltaR[AB] );
+    RealOpenMM  rAB  = SQRT( rAB2 );
+    RealOpenMM  rCB2 = AmoebaReferenceForce::getNormSquared3( deltaR[CB] );
+    RealOpenMM  rCB  = SQRT( rCB2 );
+
+    AmoebaReferenceForce::getCrossProduct( deltaR[CB], deltaR[AB], deltaR[CBxAB] );
+    RealOpenMM  rP   = AmoebaReferenceForce::getNorm3( deltaR[CBxAB] );
+    if( rP <= zero ){
+       return zero;
+    }
+    RealOpenMM dot    = AmoebaReferenceForce::getDotProduct3( deltaR[CB], deltaR[AB] );
+    RealOpenMM cosine = dot/(rAB*rCB);
+ 
+    RealOpenMM angle;
+    if( cosine >= one ){
+       angle = zero;
+    } else if( cosine <= -one ){
+       angle = PI_M;
+    } else {
+       angle = RADIAN*ACOS(cosine);
+    }
+ 
+    RealOpenMM termA = -RADIAN/(rAB2*rP);
+    RealOpenMM termC =  RADIAN/(rCB2*rP);
+ 
+    // P = CBxAB
+ 
+    AmoebaReferenceForce::getCrossProduct( deltaR[AB], deltaR[CBxAB], deltaR[ABxP] );
+    AmoebaReferenceForce::getCrossProduct( deltaR[CB], deltaR[CBxAB], deltaR[CBxP] );
+    for( int ii = 0; ii < 3; ii++ ){
+       deltaR[ABxP][ii] *= termA;
+       deltaR[CBxP][ii] *= termC;
+    }
+ 
+    RealOpenMM dr    = rAB - lengthAB + rCB - lengthCB;
+    termA            = one/rAB;
+    termC            = one/rCB;
+ 
+    // ---------------------------------------------------------------------------------------
+ 
+    // forces
+ 
+    // calculate forces for atoms a, b, c
+    // the force for b is then -( a + c)
+ 
+    std::vector<RealOpenMM> subForce[LastAtomIndex];
+    for( int ii = 0; ii < LastAtomIndex; ii++ ){
+        subForce[ii].resize(3);
+    }
+    RealOpenMM dt = angle - idealAngle*RADIAN;
+    for( int jj = 0; jj < 3; jj++ ){
+       subForce[A][jj] = kParameter*(dt*termA*deltaR[AB][jj] + dr*deltaR[ABxP][jj] );
+       subForce[C][jj] = kParameter*(dt*termC*deltaR[CB][jj] + dr*deltaR[CBxP][jj] );
+       subForce[B][jj] = -( subForce[A][jj] + subForce[C][jj] );
+    }
+ 
+    // add in forces
+ 
+    for( int jj = 0; jj < LastAtomIndex; jj++ ){
+        forces[jj][0] -= subForce[jj][0];
+        forces[jj][1] -= subForce[jj][1];
+        forces[jj][2] -= subForce[jj][2];
+    }
+ 
+    // ---------------------------------------------------------------------------------------
+ 
+    // return energy
+    return (kParameter*dt*dr);
+}

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceStretchBendForce.h

+
+/* Portions copyright (c) 2006 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#ifndef __AmoebaReferenceStretchBendForce_H__
+#define __AmoebaReferenceStretchBendForce_H__
+
+#include "SimTKUtilities/SimTKOpenMMRealType.h"
+
+// ---------------------------------------------------------------------------------------
+
+class AmoebaReferenceStretchBendForce {
+
+public:
+ 
+    /**---------------------------------------------------------------------------------------
+       
+       Constructor
+       
+       --------------------------------------------------------------------------------------- */
+ 
+    AmoebaReferenceStretchBendForce( );
+ 
+    /**---------------------------------------------------------------------------------------
+       
+       Destructor
+       
+          --------------------------------------------------------------------------------------- */
+ 
+    ~AmoebaReferenceStretchBendForce( );
+
+    /**---------------------------------------------------------------------------------------
+    
+       Calculate Amoeba harmonic angle ixn (force and energy)
+    
+       @param positionAtomA           Cartesian coordinates of atom A
+       @param positionAtomB           Cartesian coordinates of atom B
+       @param positionAtomC           Cartesian coordinates of atom C
+       @param lengthAB                ideal AB bondlength
+       @param lengthCB                ideal CB bondlength
+       @param idealAngle              ideal angle
+       @param kParameter              k
+       @param forces                  force vector
+    
+       @return energy
+    
+       --------------------------------------------------------------------------------------- */
+    static RealOpenMM calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
+                                               const RealOpenMM* positionAtomC,
+                                               RealOpenMM lengthAB,      RealOpenMM lengthCB,
+                                               RealOpenMM idealAngle,    RealOpenMM kParameter,
+                                               RealOpenMM** forces );
+ 
+};
+
+// ---------------------------------------------------------------------------------------
+
+#endif // _AmoebaReferenceStretchBendForce___