Code for AmoebaReferenceHarmonicInPlaneAngleForce

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

@@ -41,8 +41,8 @@ fprintf( stderr,"In registerKernelFactories AmoebaReferenceKernelFactory\n" ); f
              AmoebaReferenceKernelFactory* factory = new AmoebaReferenceKernelFactory();
              platform.registerKernelFactory(CalcAmoebaHarmonicBondForceKernel::Name(), factory);
              platform.registerKernelFactory(CalcAmoebaHarmonicAngleForceKernel::Name(), factory);
-/*
              platform.registerKernelFactory(CalcAmoebaHarmonicInPlaneAngleForceKernel::Name(), factory);
+/*
              platform.registerKernelFactory(CalcAmoebaTorsionForceKernel::Name(), factory);
              platform.registerKernelFactory(CalcAmoebaPiTorsionForceKernel::Name(), factory);
              platform.registerKernelFactory(CalcAmoebaStretchBendForceKernel::Name(), factory);
@@ -68,10 +68,10 @@ KernelImpl* AmoebaReferenceKernelFactory::createKernelImpl(std::string name, con
     if (name == CalcAmoebaHarmonicAngleForceKernel::Name())
         return new ReferenceCalcAmoebaHarmonicAngleForceKernel(name, platform, context.getSystem());
 
-/*
     if (name == CalcAmoebaHarmonicInPlaneAngleForceKernel::Name())
         return new ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel(name, platform, context.getSystem());
 
+/*
     if (name == CalcAmoebaTorsionForceKernel::Name())
         return new ReferenceCalcAmoebaTorsionForceKernel(name, platform, context.getSystem());
 

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

@@ -27,6 +27,7 @@
 #include "AmoebaReferenceKernels.h"
 #include "AmoebaReferenceHarmonicBondForce.h"
 #include "AmoebaReferenceHarmonicAngleForce.h"
+#include "AmoebaReferenceHarmonicInPlaneAngleForce.h"
 #include "ReferencePlatform.h"
 #include "openmm/internal/ContextImpl.h"
 //#include "internal/AmoebaMultipoleForceImpl.h"
@@ -155,48 +156,57 @@ double ReferenceCalcAmoebaHarmonicAngleForceKernel::execute(ContextImpl& context
     return static_cast<double>(energy);
 }
 
-//ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel::ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel(std::string name, const Platform& platform, System& system) : 
-//          CalcAmoebaHarmonicInPlaneAngleForceKernel(name, platform), system(system) {
-//}
-//
-//ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel::~ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel() {
-//}
-//
-//void ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel::initialize(const System& system, const AmoebaHarmonicInPlaneAngleForce& force) {
-//
-//    numAngles = force.getNumAngles();
-//
-//    std::vector<int> particle1(numAngles);
-//    std::vector<int> particle2(numAngles);
-//    std::vector<int> particle3(numAngles);
-//    std::vector<int> particle4(numAngles);
-//    std::vector<RealOpenMM> angle(numAngles);
-//    std::vector<RealOpenMM> k(numAngles);
-//
-//    for (int i = 0; i < numAngles; i++) {
-//        double angleValue, kQuadratic;
-//        force.getAngleParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], angleValue, kQuadratic);
-//        //angle[i]            = static_cast<RealOpenMM>( (angleValue*RadiansToDegrees) );
-//        angle[i]            = static_cast<RealOpenMM>( angleValue );
-//        k[i]                = static_cast<RealOpenMM>( kQuadratic );
-//    }
-///*
-//    gpuSetAmoebaInPlaneAngleParameters(data.getAmoebaGpu(), particle1, particle2, particle3, particle4, angle, k,
-//                                       static_cast<RealOpenMM>( force.getAmoebaGlobalHarmonicInPlaneAngleCubic()),
-//                                       static_cast<RealOpenMM>( force.getAmoebaGlobalHarmonicInPlaneAngleQuartic()),
-//                                       static_cast<RealOpenMM>( force.getAmoebaGlobalHarmonicInPlaneAnglePentic()),
-//                                       static_cast<RealOpenMM>( force.getAmoebaGlobalHarmonicInPlaneAngleSextic() ) );
-//*/
-//
-//}
-//
-//double ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-//    if( data.getAmoebaLocalForcesKernel() == this ){
-//        computeAmoebaLocalForces( data );
-//    }
-//    return 0.0;
-//}
-//
+ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel::ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel(std::string name, const Platform& platform, System& system) : 
+          CalcAmoebaHarmonicInPlaneAngleForceKernel(name, platform), system(system) {
+}
+
+ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel::~ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel() {
+}
+
+void ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel::initialize(const System& system, const AmoebaHarmonicInPlaneAngleForce& force) {
+
+    numAngles = force.getNumAngles();
+    for (int ii = 0; ii < numAngles; ii++) {
+        int particle1Index, particle2Index, particle3Index, particle4Index;
+        double angleValue, k;
+        force.getAngleParameters(ii, particle1Index, particle2Index, particle3Index, particle4Index, angleValue, k);
+        particle1.push_back( particle1Index ); 
+        particle2.push_back( particle2Index ); 
+        particle3.push_back( particle3Index ); 
+        particle4.push_back( particle4Index ); 
+        angle.push_back(       static_cast<RealOpenMM>( angleValue ) );
+        kQuadratic.push_back(  static_cast<RealOpenMM>( k ) );
+    }
+    globalHarmonicInPlaneAngleCubic    = static_cast<RealOpenMM>(force.getAmoebaGlobalHarmonicInPlaneAngleCubic());
+    globalHarmonicInPlaneAngleQuartic  = static_cast<RealOpenMM>(force.getAmoebaGlobalHarmonicInPlaneAngleQuartic());
+    globalHarmonicInPlaneAnglePentic   = static_cast<RealOpenMM>(force.getAmoebaGlobalHarmonicInPlaneAnglePentic());
+    globalHarmonicInPlaneAngleSextic   = static_cast<RealOpenMM>(force.getAmoebaGlobalHarmonicInPlaneAngleSextic());
+}
+
+double ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel::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 < numAngles; ii++) {
+        int particle1Index      = particle1[ii];
+        int particle2Index      = particle2[ii];
+        int particle3Index      = particle3[ii];
+        int particle4Index      = particle4[ii];
+        RealOpenMM idealAngle   = angle[ii];
+        RealOpenMM angleK       = kQuadratic[ii];
+        RealOpenMM* forces[4];
+        forces[0]               = forceData[particle1Index];
+        forces[1]               = forceData[particle2Index];
+        forces[2]               = forceData[particle3Index];
+        forces[3]               = forceData[particle4Index];
+        energy                 += AmoebaReferenceHarmonicInPlaneAngleForce::calculateForceAndEnergy( 
+                                       posData[particle1Index], posData[particle2Index], posData[particle3Index], posData[particle4Index],
+                                       idealAngle, angleK, globalHarmonicInPlaneAngleCubic, globalHarmonicInPlaneAngleQuartic,
+                                       globalHarmonicInPlaneAnglePentic, globalHarmonicInPlaneAngleSextic, forces );
+    }
+    return static_cast<double>(energy);
+}
+
 //ReferenceCalcAmoebaTorsionForceKernel::ReferenceCalcAmoebaTorsionForceKernel(std::string name, const Platform& platform, System& system) :
 //             CalcAmoebaTorsionForceKernel(name, platform), system(system) {
 //    data.incrementKernelCount();

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

@@ -106,41 +106,50 @@ private:
     System& system;
 };
 
-// /**
-//  * This kernel is invoked by AmoebaHarmonicInPlaneAngleForce to calculate the forces acting on the system and the energy of the system.
-//  */
-// class ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel : public CalcAmoebaHarmonicInPlaneAngleForceKernel {
-// public:
-//     ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel(std::string name, const Platform& platform, AmoebaReferenceData& data, System& system);
-//     ~ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel();
-//     /**
-//      * Initialize the kernel.
-//      * 
-//      * @param system     the System this kernel will be applied to
-//      * @param force      the AmoebaHarmonicInPlaneAngleForce this kernel will be used for
-//      */
-//     void initialize(const System& system, const AmoebaHarmonicInPlaneAngleForce& 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 numAngles;
-//     AmoebaReferenceData& data;
-//     System& system;
-// };
-// 
+/**
+ * This kernel is invoked by AmoebaHarmonicInPlaneAngleForce to calculate the forces acting on the system and the energy of the system.
+ */
+class ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel : public CalcAmoebaHarmonicInPlaneAngleForceKernel {
+public:
+    ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel(std::string name, const Platform& platform, System& system);
+    ~ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param force      the AmoebaHarmonicInPlaneAngleForce this kernel will be used for
+     */
+    void initialize(const System& system, const AmoebaHarmonicInPlaneAngleForce& 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 numAngles;
+    std::vector<int>   particle1;
+    std::vector<int>   particle2;
+    std::vector<int>   particle3;
+    std::vector<int>   particle4;
+    std::vector<RealOpenMM> angle;
+    std::vector<RealOpenMM> kQuadratic;
+    RealOpenMM globalHarmonicInPlaneAngleCubic;
+    RealOpenMM globalHarmonicInPlaneAngleQuartic;
+    RealOpenMM globalHarmonicInPlaneAnglePentic;
+    RealOpenMM globalHarmonicInPlaneAngleSextic;
+    System& system;
+};
+
 // /**
 //  * This kernel is invoked by AmoebaTorsionForce to calculate the forces acting on the system and the energy of the system.
 //  */
 // class ReferenceCalcAmoebaTorsionForceKernel : public CalcAmoebaTorsionForceKernel {
 // public:
-//     ReferenceCalcAmoebaTorsionForceKernel(std::string name, const Platform& platform, AmoebaReferenceData& data, System& system);
+//     ReferenceCalcAmoebaTorsionForceKernel(std::string name, const Platform& platform, System& system);
 //     ~ReferenceCalcAmoebaTorsionForceKernel();
 //     /**
 //      * Initialize the kernel.
@@ -160,7 +169,6 @@ private:
 //     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 // private:
 //     int numTorsions;
-//     AmoebaReferenceData& data;
 //     System& system;
 // };
 // 
@@ -169,7 +177,7 @@ private:
 //  */
 // class ReferenceCalcAmoebaPiTorsionForceKernel : public CalcAmoebaPiTorsionForceKernel {
 // public:
-//     ReferenceCalcAmoebaPiTorsionForceKernel(std::string name, const Platform& platform, AmoebaReferenceData& data, System& system);
+//     ReferenceCalcAmoebaPiTorsionForceKernel(std::string name, const Platform& platform, System& system);
 //     ~ReferenceCalcAmoebaPiTorsionForceKernel();
 //     /**
 //      * Initialize the kernel.
@@ -189,7 +197,6 @@ private:
 //     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 // private:
 //     int numPiTorsions;
-//     AmoebaReferenceData& data;
 //     System& system;
 // };
 // 
@@ -198,7 +205,7 @@ private:
 //  */
 // class ReferenceCalcAmoebaStretchBendForceKernel : public CalcAmoebaStretchBendForceKernel {
 // public:
-//     ReferenceCalcAmoebaStretchBendForceKernel(std::string name, const Platform& platform, AmoebaReferenceData& data, System& system);
+//     ReferenceCalcAmoebaStretchBendForceKernel(std::string name, const Platform& platform, System& system);
 //     ~ReferenceCalcAmoebaStretchBendForceKernel();
 //     /**
 //      * Initialize the kernel.
@@ -218,7 +225,6 @@ private:
 //     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 // private:
 //     int numStretchBends;
-//     AmoebaReferenceData& data;
 //     System& system;
 // };
 // 
@@ -227,7 +233,7 @@ private:
 //  */
 // class ReferenceCalcAmoebaOutOfPlaneBendForceKernel : public CalcAmoebaOutOfPlaneBendForceKernel {
 // public:
-//     ReferenceCalcAmoebaOutOfPlaneBendForceKernel(std::string name, const Platform& platform, AmoebaReferenceData& data, System& system);
+//     ReferenceCalcAmoebaOutOfPlaneBendForceKernel(std::string name, const Platform& platform, System& system);
 //     ~ReferenceCalcAmoebaOutOfPlaneBendForceKernel();
 //     /**
 //      * Initialize the kernel.
@@ -247,7 +253,6 @@ private:
 //     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 // private:
 //     int numOutOfPlaneBends;
-//     AmoebaReferenceData& data;
 //     System& system;
 // };
 // 
@@ -256,7 +261,7 @@ private:
 //  */
 // class ReferenceCalcAmoebaTorsionTorsionForceKernel : public CalcAmoebaTorsionTorsionForceKernel {
 // public:
-//     ReferenceCalcAmoebaTorsionTorsionForceKernel(std::string name, const Platform& platform, AmoebaReferenceData& data, System& system);
+//     ReferenceCalcAmoebaTorsionTorsionForceKernel(std::string name, const Platform& platform, System& system);
 //     ~ReferenceCalcAmoebaTorsionTorsionForceKernel();
 //     /**
 //      * Initialize the kernel.
@@ -277,7 +282,6 @@ private:
 // private:
 //     int numTorsionTorsions;
 //     int numTorsionTorsionGrids;
-//     AmoebaReferenceData& data;
 //     System& system;
 // };
 // 
@@ -286,7 +290,7 @@ private:
 //  */
 // class ReferenceCalcAmoebaMultipoleForceKernel : public CalcAmoebaMultipoleForceKernel {
 // public:
-//     ReferenceCalcAmoebaMultipoleForceKernel(std::string name, const Platform& platform, AmoebaReferenceData& data, System& system);
+//     ReferenceCalcAmoebaMultipoleForceKernel(std::string name, const Platform& platform, System& system);
 //     ~ReferenceCalcAmoebaMultipoleForceKernel();
 //     /**
 //      * Initialize the kernel.
@@ -306,7 +310,6 @@ private:
 //     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 // private:
 //     int numMultipoles;
-//     AmoebaReferenceData& data;
 //     System& system;
 // };
 // 
@@ -315,7 +318,7 @@ private:
 //  */
 // class ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel : public CalcAmoebaGeneralizedKirkwoodForceKernel {
 // public:
-//     ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel(std::string name, const Platform& platform, AmoebaReferenceData& data, System& system);
+//     ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel(std::string name, const Platform& platform, System& system);
 //     ~ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel();
 //     /**
 //      * Initialize the kernel.
@@ -334,7 +337,6 @@ private:
 //      */
 //     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 // private:
-//     AmoebaReferenceData& data;
 //     System& system;
 // };
 // 
@@ -343,7 +345,7 @@ private:
 //  */
 // class ReferenceCalcAmoebaVdwForceKernel : public CalcAmoebaVdwForceKernel {
 // public:
-//     ReferenceCalcAmoebaVdwForceKernel(std::string name, const Platform& platform, AmoebaReferenceData& data, System& system);
+//     ReferenceCalcAmoebaVdwForceKernel(std::string name, const Platform& platform, System& system);
 //     ~ReferenceCalcAmoebaVdwForceKernel();
 //     /**
 //      * Initialize the kernel.
@@ -362,7 +364,6 @@ private:
 //      */
 //     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 // private:
-//     AmoebaReferenceData& data;
 //     System& system;
 // };
 // 
@@ -371,7 +372,7 @@ private:
 //  */
 // class ReferenceCalcAmoebaWcaDispersionForceKernel : public CalcAmoebaWcaDispersionForceKernel {
 // public:
-//     ReferenceCalcAmoebaWcaDispersionForceKernel(std::string name, const Platform& platform, AmoebaReferenceData& data, System& system);
+//     ReferenceCalcAmoebaWcaDispersionForceKernel(std::string name, const Platform& platform, System& system);
 //     ~ReferenceCalcAmoebaWcaDispersionForceKernel();
 //     /**
 //      * Initialize the kernel.
@@ -390,7 +391,6 @@ private:
 //      */
 //     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 // private:
-//     AmoebaReferenceData& data;
 //     System& system;
 // };
 

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicInPlaneAngleForce.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 "AmoebaReferenceHarmonicInPlaneAngleForce.h"
+#include <vector>
+
+/**---------------------------------------------------------------------------------------
+
+   Get dEdT and energy prefactor given cosine of angle :: the calculation for different
+   the force types is identical 
+
+   @param  cosine               cosine of angle
+   @param  idealAngle           ideal angle
+   @param  angleK               angle k (quadratic prefactor)
+   @param  angleCubic           cubic prefactor
+   @param  angleQuartic         quartiic prefactor
+   @param  anglePentic          pentic prefactor
+   @param  angleSextic          sextic prefactor
+   @param  dEdR                 dEdR
+
+   @return energy
+
+   --------------------------------------------------------------------------------------- */
+
+RealOpenMM AmoebaReferenceHarmonicInPlaneAngleForce::getPrefactorsGivenAngleCosine( RealOpenMM cosine,
+                                                                                    RealOpenMM idealAngle,     RealOpenMM angleK,
+                                                                                    RealOpenMM angleCubic,     RealOpenMM angleQuartic,
+                                                                                    RealOpenMM anglePentic,    RealOpenMM angleSextic,
+                                                                                    RealOpenMM* dEdR ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const RealOpenMM zero          = 0.0;
+   static const RealOpenMM one           = 1.0;
+   static const RealOpenMM two           = 2.0;
+   static const RealOpenMM three         = 3.0;
+   static const RealOpenMM four          = 4.0;
+   static const RealOpenMM five          = 5.0;
+   static const RealOpenMM six           = 6.0;
+
+   // static const std::string methodName = "AmoebaReferenceHarmonicInPlaneAngleForce::getPrefactorsGivenAngleCosine";
+
+   // ---------------------------------------------------------------------------------------
+
+   RealOpenMM angle;
+   if( cosine >= one ){
+      angle = zero;
+   } else if( cosine <= -one ){
+      angle = RADIAN*PI_M;
+   } else {
+      angle = RADIAN*ACOS(cosine);
+   }   
+   RealOpenMM deltaIdeal         = angle - idealAngle;
+   RealOpenMM deltaIdeal2        = deltaIdeal*deltaIdeal;
+   RealOpenMM deltaIdeal3        = deltaIdeal*deltaIdeal2;
+   RealOpenMM deltaIdeal4        = deltaIdeal2*deltaIdeal2;
+
+   *dEdR                         = ( two + three*angleCubic*deltaIdeal +
+                                     four*angleQuartic*deltaIdeal2     + 
+                                     five*anglePentic*deltaIdeal3      +
+                                     six*angleSextic*deltaIdeal4 );
+
+   *dEdR                        *= RADIAN*angleK*deltaIdeal;
+
+   RealOpenMM energy             = 1.0f + angleCubic*deltaIdeal + angleQuartic*deltaIdeal2 +
+                                   anglePentic*deltaIdeal3 + angleSextic*deltaIdeal4;
+   energy                       *= angleK*deltaIdeal2;
+
+   return energy;
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   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 AmoebaReferenceHarmonicInPlaneAngleForce::calculateForceAndEnergy( RealOpenMM* positionAtomA, RealOpenMM* positionAtomB,
+                                                                              RealOpenMM* positionAtomC, RealOpenMM* positionAtomD,
+                                                                              RealOpenMM angle,          RealOpenMM angleK,
+                                                                              RealOpenMM angleCubic,     RealOpenMM angleQuartic,
+                                                                              RealOpenMM anglePentic,    RealOpenMM angleSextic,
+                                                                              RealOpenMM** forces ){
+
+   // ---------------------------------------------------------------------------------------
+
+    //static const std::string methodName = "AmoebaReferenceHarmonicInPlaneAngleForce::calculateHarmonicForce";
+ 
+    static const RealOpenMM zero          = 0.0;
+    static const RealOpenMM one           = 1.0;
+    static const RealOpenMM two           = 2.0;
+
+   // ---------------------------------------------------------------------------------------
+
+    // T   = AD x CD
+    // P   = B + T*delta
+    // AP  = A - P
+    // CP  = A - P
+    // M   = CP x AP
+    enum { AD, BD, CD, T, AP, P, CP, M, APxM, CPxM, ADxBD, BDxCD, TxCD, ADxT, dBxAD, CDxdB, LastDeltaAtomIndex };
+
+    std::vector<RealOpenMM> deltaR[LastDeltaAtomIndex];
+    for( int ii = 0; ii < LastDeltaAtomIndex; ii++ ){
+        deltaR[ii].resize(3);
+    }
+    AmoebaReferenceForce::loadDeltaR( positionAtomD, positionAtomA, deltaR[AD] );
+    AmoebaReferenceForce::loadDeltaR( positionAtomD, positionAtomB, deltaR[BD] );
+    AmoebaReferenceForce::loadDeltaR( positionAtomD, positionAtomC, deltaR[CD] );
+ 
+    AmoebaReferenceForce::getCrossProduct( deltaR[AD], deltaR[CD], deltaR[T] );
+
+    RealOpenMM rT2     = AmoebaReferenceForce::getNormSquared3( deltaR[T] );
+    RealOpenMM delta   = AmoebaReferenceForce::getDotProduct3( deltaR[T], deltaR[BD] )/rT2;
+         delta        *= -one;
+ 
+    for( int ii = 0; ii < 3; ii++ ){
+       deltaR[P][ii]  = positionAtomB[ii] + deltaR[T][ii]*delta;
+       deltaR[AP][ii] = positionAtomA[ii] - deltaR[P][ii];
+       deltaR[CP][ii] = positionAtomC[ii] - deltaR[P][ii];
+    }
+ 
+    RealOpenMM rAp2 = AmoebaReferenceForce::getNormSquared3( deltaR[AP] );
+    RealOpenMM rCp2 = AmoebaReferenceForce::getNormSquared3( deltaR[CP] );
+    if( rAp2 <= zero && rCp2 <= zero ){
+       return zero;
+    }
+ 
+    AmoebaReferenceForce::getCrossProduct( deltaR[CP], deltaR[AP], deltaR[M] );
+ 
+    RealOpenMM rm = AmoebaReferenceForce::getNorm3( deltaR[M] );
+    if( rm < 1.0e-06 ){
+       rm = 1.0e-06;
+    }
+ 
+    RealOpenMM dot     = AmoebaReferenceForce::getDotProduct3( deltaR[AP], deltaR[CP] );
+    RealOpenMM cosine  = dot/SQRT( rAp2*rCp2 );
+ 
+    RealOpenMM dEdR;
+    RealOpenMM energy = getPrefactorsGivenAngleCosine( cosine, angle, angleK, angleCubic, angleQuartic,
+                                                       anglePentic, angleSextic, &dEdR );
+
+    RealOpenMM termA   = -dEdR/(rAp2*rm);
+    RealOpenMM termC   =  dEdR/(rCp2*rm);
+ 
+    AmoebaReferenceForce::getCrossProduct( deltaR[AP], deltaR[M], deltaR[APxM] );
+    AmoebaReferenceForce::getCrossProduct( deltaR[CP], deltaR[M], deltaR[CPxM] );
+ 
+    // forces will be gathered here
+ 
+    enum { dA, dB, dC, dD, LastDIndex };
+    std::vector<RealOpenMM> forceTerm[LastDIndex];
+    for( int ii = 0; ii < LastDIndex; ii++ ){
+        forceTerm[ii].resize(3);
+    }
+ 
+    for( int ii = 0; ii < 3; ii++ ){
+       forceTerm[dA][ii] = deltaR[APxM][ii]*termA;
+       forceTerm[dC][ii] = deltaR[CPxM][ii]*termC;
+       forceTerm[dB][ii] = -one*( forceTerm[dA][ii] + forceTerm[dC][ii] );
+    }
+ 
+    RealOpenMM pTrT2  = AmoebaReferenceForce::getDotProduct3( forceTerm[dB], deltaR[T] );
+               pTrT2 /= rT2;
+ 
+    AmoebaReferenceForce::getCrossProduct( deltaR[CD], forceTerm[dB], deltaR[CDxdB] );
+    AmoebaReferenceForce::getCrossProduct( forceTerm[dB], deltaR[AD], deltaR[dBxAD] );
+ 
+    if( FABS( pTrT2 ) > 1.0e-08 ){
+ 
+       RealOpenMM delta2 = delta*two;
+ 
+       AmoebaReferenceForce::getCrossProduct( deltaR[BD], deltaR[CD], deltaR[BDxCD] );
+       AmoebaReferenceForce::getCrossProduct( deltaR[T],  deltaR[CD], deltaR[TxCD]  );
+       AmoebaReferenceForce::getCrossProduct( deltaR[AD], deltaR[BD], deltaR[ADxBD] );
+       AmoebaReferenceForce::getCrossProduct( deltaR[AD], deltaR[T],  deltaR[ADxT]  );
+       for( int ii = 0; ii < 3; ii++ ){
+    
+          RealOpenMM term     = deltaR[BDxCD][ii] + delta2*deltaR[TxCD][ii];
+          forceTerm[dA][ii]  += delta*deltaR[CDxdB][ii] + term*pTrT2;
+    
+               term           = deltaR[ADxBD][ii] + delta2*deltaR[ADxT][ii];
+          forceTerm[dC][ii]  += delta*deltaR[dBxAD][ii] + term*pTrT2;
+    
+          forceTerm[dD][ii]  = -( forceTerm[dA][ii] + forceTerm[dB][ii] + forceTerm[dC][ii] );
+       }
+    } else {
+       for( int ii = 0; ii < 3; ii++ ){
+    
+          forceTerm[dA][ii] += delta*deltaR[CDxdB][ii];
+          forceTerm[dC][ii] += delta*deltaR[dBxAD][ii];
+ 
+          forceTerm[dD][ii]  = -( forceTerm[dA][ii] + forceTerm[dB][ii] + forceTerm[dC][ii] );
+       }
+    }
+ 
+    // accumulate forces
+ 
+    for( int jj = 0; jj < 4; jj++ ){
+        forces[jj][0] -= forceTerm[jj][0];
+        forces[jj][1] -= forceTerm[jj][1];
+        forces[jj][2] -= forceTerm[jj][2];
+    }
+ 
+    return energy;
+
+}

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicInPlaneAngleForce.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 __AmoebaReferenceHarmonicInPlaneAngleForce_H__
+#define __AmoebaReferenceHarmonicInPlaneAngleForce_H__
+
+#include "SimTKUtilities/SimTKOpenMMRealType.h"
+
+// ---------------------------------------------------------------------------------------
+
+class AmoebaReferenceHarmonicInPlaneAngleForce {
+
+public:
+ 
+    /**---------------------------------------------------------------------------------------
+       
+       Constructor
+       
+       --------------------------------------------------------------------------------------- */
+ 
+    AmoebaReferenceHarmonicInPlaneAngleForce( );
+ 
+    /**---------------------------------------------------------------------------------------
+       
+       Destructor
+       
+          --------------------------------------------------------------------------------------- */
+ 
+    ~AmoebaReferenceHarmonicInPlaneAngleForce( );
+ 
+    /**---------------------------------------------------------------------------------------
+    
+       Get dEdT and energy prefactor given cosine of angle :: the calculation for different
+       the force types is identical 
+    
+       @param  cosine               cosine of angle
+       @param  idealAngle           ideal angle
+       @param  angleK               angle k (quadratic prefactor)
+       @param  angleCubic           cubic prefactor
+       @param  angleQuartic         quartiic prefactor
+       @param  anglePentic          pentic prefactor
+       @param  angleSextic          sextic prefactor
+       @param  dEdR                 dEdR
+    
+       @return energy
+    
+       --------------------------------------------------------------------------------------- */
+    
+    static RealOpenMM getPrefactorsGivenAngleCosine( RealOpenMM cosine, RealOpenMM idealAngle, RealOpenMM angleK,
+                                                     RealOpenMM angleCubic,     RealOpenMM angleQuartic,
+                                                     RealOpenMM anglePentic,    RealOpenMM angleSextic,
+                                                     RealOpenMM* dEdR );
+    
+    /**---------------------------------------------------------------------------------------
+    
+       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
+    
+       --------------------------------------------------------------------------------------- */
+    
+    static RealOpenMM calculateForceAndEnergy( RealOpenMM* positionAtomA, RealOpenMM* positionAtomB,
+                                               RealOpenMM* positionAtomC, RealOpenMM* positionAtomD,
+                                               RealOpenMM angle,          RealOpenMM angleK,
+                                               RealOpenMM angleCubic,     RealOpenMM angleQuartic,
+                                               RealOpenMM anglePentic,    RealOpenMM angleSextic,
+                                               RealOpenMM** forces );
+         
+};
+
+// ---------------------------------------------------------------------------------------
+
+#endif // _AmoebaReferenceHarmonicInPlaneAngleForce___