Added AmoebaReferenceOutOfPlaneBend

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

@@ -45,8 +45,8 @@ fprintf( stderr,"In registerKernelFactories AmoebaReferenceKernelFactory\n" ); f
              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);
              platform.registerKernelFactory(CalcAmoebaGeneralizedKirkwoodForceKernel::Name(), factory);
@@ -80,10 +80,10 @@ KernelImpl* AmoebaReferenceKernelFactory::createKernelImpl(std::string name, con
     if (name == CalcAmoebaStretchBendForceKernel::Name())
         return new ReferenceCalcAmoebaStretchBendForceKernel(name, platform, context.getSystem());
 
-/*
     if (name == CalcAmoebaOutOfPlaneBendForceKernel::Name())
         return new ReferenceCalcAmoebaOutOfPlaneBendForceKernel(name, platform, context.getSystem());
 
+/*
     if (name == CalcAmoebaTorsionTorsionForceKernel::Name())
         return new ReferenceCalcAmoebaTorsionTorsionForceKernel(name, platform, context.getSystem());
 

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

@@ -31,6 +31,7 @@
 #include "AmoebaReferenceTorsionForce.h"
 #include "AmoebaReferencePiTorsionForce.h"
 #include "AmoebaReferenceStretchBendForce.h"
+#include "AmoebaReferenceOutOfPlaneBendForce.h"
 #include "ReferencePlatform.h"
 #include "openmm/internal/ContextImpl.h"
 //#include "internal/AmoebaMultipoleForceImpl.h"
@@ -374,48 +375,58 @@ double ReferenceCalcAmoebaStretchBendForceKernel::execute(ContextImpl& context,
     return static_cast<double>(energy);
 }
 
-//ReferenceCalcAmoebaOutOfPlaneBendForceKernel::ReferenceCalcAmoebaOutOfPlaneBendForceKernel(std::string name, const Platform& platform, System& system) :
-//          CalcAmoebaOutOfPlaneBendForceKernel(name, platform), system(system) {
-//    data.incrementKernelCount();
-//}
-//
-//ReferenceCalcAmoebaOutOfPlaneBendForceKernel::~ReferenceCalcAmoebaOutOfPlaneBendForceKernel() {
-//    data.decrementKernelCount();
-//}
-//
-//void ReferenceCalcAmoebaOutOfPlaneBendForceKernel::initialize(const System& system, const AmoebaOutOfPlaneBendForce& force) {
-//
-//    data.setAmoebaLocalForcesKernel( this );
-//    numOutOfPlaneBends                     = force.getNumOutOfPlaneBends();
-//
-//    std::vector<int>   particle1(numOutOfPlaneBends);
-//    std::vector<int>   particle2(numOutOfPlaneBends);
-//    std::vector<int>   particle3(numOutOfPlaneBends);
-//    std::vector<int>   particle4(numOutOfPlaneBends);
-//    std::vector<RealOpenMM> kParameters(numOutOfPlaneBends);
-//
-//    for (int i = 0; i < numOutOfPlaneBends; i++) {
-//
-//        double k;
-//
-//        force.getOutOfPlaneBendParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], k);
-//        kParameters[i] = static_cast<RealOpenMM>(k);
-//    }
-//    gpuSetAmoebaOutOfPlaneBendParameters(data.getAmoebaGpu(), particle1, particle2, particle3, particle4, kParameters,
-//                                         static_cast<RealOpenMM>( force.getAmoebaGlobalOutOfPlaneBendCubic()),
-//                                         static_cast<RealOpenMM>( force.getAmoebaGlobalOutOfPlaneBendQuartic()),
-//                                         static_cast<RealOpenMM>( force.getAmoebaGlobalOutOfPlaneBendPentic()),
-//                                         static_cast<RealOpenMM>( force.getAmoebaGlobalOutOfPlaneBendSextic() ) );
-//
-//}
-//
-//double ReferenceCalcAmoebaOutOfPlaneBendForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-//    if( data.getAmoebaLocalForcesKernel() == this ){
-//        computeAmoebaLocalForces( data );
-//    }
-//    return 0.0;
-//}
-//
+ReferenceCalcAmoebaOutOfPlaneBendForceKernel::ReferenceCalcAmoebaOutOfPlaneBendForceKernel(std::string name, const Platform& platform, System& system) :
+          CalcAmoebaOutOfPlaneBendForceKernel(name, platform), system(system) {
+}
+
+ReferenceCalcAmoebaOutOfPlaneBendForceKernel::~ReferenceCalcAmoebaOutOfPlaneBendForceKernel() {
+}
+
+void ReferenceCalcAmoebaOutOfPlaneBendForceKernel::initialize(const System& system, const AmoebaOutOfPlaneBendForce& force) {
+
+    numOutOfPlaneBends = force.getNumOutOfPlaneBends();
+    for (int ii = 0; ii < numOutOfPlaneBends; ii++) {
+
+        int particle1Index, particle2Index, particle3Index, particle4Index;
+        double k;
+
+        force.getOutOfPlaneBendParameters(ii, particle1Index, particle2Index, particle3Index, particle4Index, k);
+        particle1.push_back( particle1Index ); 
+        particle2.push_back( particle2Index ); 
+        particle3.push_back( particle3Index ); 
+        particle4.push_back( particle4Index ); 
+        kParameters.push_back( static_cast<RealOpenMM>(k) );
+    }
+    globalOutOfPlaneBendAngleCubic      = static_cast<RealOpenMM>( force.getAmoebaGlobalOutOfPlaneBendCubic());
+    globalOutOfPlaneBendAngleQuartic    = static_cast<RealOpenMM>( force.getAmoebaGlobalOutOfPlaneBendQuartic());
+    globalOutOfPlaneBendAnglePentic     = static_cast<RealOpenMM>( force.getAmoebaGlobalOutOfPlaneBendPentic());
+    globalOutOfPlaneBendAngleSextic     = static_cast<RealOpenMM>( force.getAmoebaGlobalOutOfPlaneBendSextic());
+
+}
+
+double ReferenceCalcAmoebaOutOfPlaneBendForceKernel::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 < numOutOfPlaneBends; ii++) {
+        int particle1Index      = particle1[ii];
+        int particle2Index      = particle2[ii];
+        int particle3Index      = particle3[ii];
+        int particle4Index      = particle4[ii];
+        RealOpenMM angleK       = kParameters[ii];
+        RealOpenMM* forces[4];
+        forces[0]               = forceData[particle1Index];
+        forces[1]               = forceData[particle2Index];
+        forces[2]               = forceData[particle3Index];
+        forces[3]               = forceData[particle4Index];
+        energy                 += AmoebaReferenceOutOfPlaneBendForce::calculateForceAndEnergy( 
+                                       posData[particle1Index], posData[particle2Index], posData[particle3Index], posData[particle4Index],
+                                       angleK, globalOutOfPlaneBendAngleCubic, globalOutOfPlaneBendAngleQuartic,
+                                       globalOutOfPlaneBendAnglePentic, globalOutOfPlaneBendAngleSextic, forces );
+    }
+    return static_cast<double>(energy);
+}
+
 //ReferenceCalcAmoebaTorsionTorsionForceKernel::ReferenceCalcAmoebaTorsionTorsionForceKernel(std::string name, const Platform& platform, System& system) :
 //                CalcAmoebaTorsionTorsionForceKernel(name, platform), system(system) {
 //    data.incrementKernelCount();

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

@@ -248,35 +248,44 @@ private:
     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.
-//  */
-// class ReferenceCalcAmoebaOutOfPlaneBendForceKernel : public CalcAmoebaOutOfPlaneBendForceKernel {
-// public:
-//     ReferenceCalcAmoebaOutOfPlaneBendForceKernel(std::string name, const Platform& platform, System& system);
-//     ~ReferenceCalcAmoebaOutOfPlaneBendForceKernel();
-//     /**
-//      * Initialize the kernel.
-//      * 
-//      * @param system     the System this kernel will be applied to
-//      * @param force      the AmoebaOutOfPlaneBendForce this kernel will be used for
-//      */
-//     void initialize(const System& system, const AmoebaOutOfPlaneBendForce& 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 numOutOfPlaneBends;
-//     System& system;
-// };
-// 
+
+/**
+ * This kernel is invoked by AmoebaOutOfPlaneBendForce to calculate the forces acting on the system and the energy of the system.
+ */
+class ReferenceCalcAmoebaOutOfPlaneBendForceKernel : public CalcAmoebaOutOfPlaneBendForceKernel {
+public:
+    ReferenceCalcAmoebaOutOfPlaneBendForceKernel(std::string name, const Platform& platform, System& system);
+    ~ReferenceCalcAmoebaOutOfPlaneBendForceKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param force      the AmoebaOutOfPlaneBendForce this kernel will be used for
+     */
+    void initialize(const System& system, const AmoebaOutOfPlaneBendForce& 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 numOutOfPlaneBends;
+    std::vector<int>   particle1;
+    std::vector<int>   particle2;
+    std::vector<int>   particle3;
+    std::vector<int>   particle4;
+    std::vector<RealOpenMM> kParameters;
+    RealOpenMM globalOutOfPlaneBendAngleCubic;
+    RealOpenMM globalOutOfPlaneBendAngleQuartic;
+    RealOpenMM globalOutOfPlaneBendAnglePentic;
+    RealOpenMM globalOutOfPlaneBendAngleSextic;
+    System& system;
+};
+
 // /**
 //  * This kernel is invoked by AmoebaTorsionTorsionForce to calculate the forces acting on the system and the energy of the system.
 //  */

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceOutOfPlaneBendForce.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 "AmoebaReferenceOutOfPlaneBendForce.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 AmoebaReferenceOutOfPlaneBendForce::calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
+                                                                              const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD,
+                                                                              RealOpenMM angleK,
+                                                                              RealOpenMM angleCubic,                 RealOpenMM angleQuartic,
+                                                                              RealOpenMM anglePentic,                RealOpenMM angleSextic,
+                                                                              RealOpenMM** forces ){
+
+   // ---------------------------------------------------------------------------------------
+
+    //static const std::string methodName = "AmoebaReferenceOutOfPlaneBendForce::calculateHarmonicForce";
+ 
+    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;
+
+    enum { A, B, C, D, LastAtomIndex };
+    enum { AB, CB, DB, AD, CD, LastDeltaIndex };
+ 
+    // ---------------------------------------------------------------------------------------
+ 
+    // get deltaR between various combinations of the 4 atoms
+    // and various intermediate terms
+ 
+    std::vector<RealOpenMM> deltaR[LastDeltaIndex];
+    for( int ii = 0; ii < LastDeltaIndex; ii++ ){
+        deltaR[ii].resize(3);
+    }   
+    AmoebaReferenceForce::loadDeltaR( positionAtomB, positionAtomA, deltaR[AB] );
+    AmoebaReferenceForce::loadDeltaR( positionAtomB, positionAtomC, deltaR[CB] );
+    AmoebaReferenceForce::loadDeltaR( positionAtomB, positionAtomD, deltaR[DB] );
+    AmoebaReferenceForce::loadDeltaR( positionAtomD, positionAtomA, deltaR[AD] );
+    AmoebaReferenceForce::loadDeltaR( positionAtomD, positionAtomC, deltaR[CD] );
+
+    RealOpenMM rDB2  = AmoebaReferenceForce::getNormSquared3( deltaR[DB] );
+    RealOpenMM rAD2  = AmoebaReferenceForce::getNormSquared3( deltaR[AD] );
+    RealOpenMM rCD2  = AmoebaReferenceForce::getNormSquared3( deltaR[CD] );
+ 
+    std::vector<RealOpenMM> tempVector(3);
+    AmoebaReferenceForce::getCrossProduct( deltaR[CB], deltaR[DB], tempVector );
+    RealOpenMM   eE  = AmoebaReferenceForce::getDotProduct3( deltaR[AB], tempVector  );
+    RealOpenMM  dot  = AmoebaReferenceForce::getDotProduct3( deltaR[AD],  deltaR[CD] );
+    RealOpenMM   cc  = rAD2*rCD2 - dot*dot;
+ 
+    if( rDB2 <= zero || cc == zero ){
+       return zero;
+    }
+    RealOpenMM bkk2   = rDB2 - eE*eE/cc;
+    RealOpenMM cosine = SQRT(bkk2/rDB2);
+    RealOpenMM angle;
+    if( cosine >= one ){
+       angle = zero;
+    } else if( cosine <= -one ){
+       angle = PI_M;
+    } else {
+       angle = RADIAN*ACOS(cosine);
+    }
+ 
+    // chain rule
+ 
+    RealOpenMM dt    = angle;
+    RealOpenMM dt2   = dt*dt;
+    RealOpenMM dt3   = dt2*dt;
+    RealOpenMM dt4   = dt2*dt2;
+ 
+    RealOpenMM dEdDt = two + three*angleCubic*dt + four*angleQuartic*dt2 +
+                       five*anglePentic*dt3 + six*angleSextic*dt4;
+ 
+          dEdDt     *= angleK*dt*RADIAN;
+ 
+    RealOpenMM dEdCos  = dEdDt/SQRT(cc*bkk2);
+    if( eE > zero ){
+       dEdCos *= -one;
+    }
+ 
+    RealOpenMM term = eE/cc;
+ 
+    std::vector<RealOpenMM> dccd[LastAtomIndex];
+    std::vector<RealOpenMM> deed[LastAtomIndex];
+    std::vector<RealOpenMM> subForce[LastAtomIndex];
+    for( int ii = 0; ii < LastAtomIndex; ii++ ){
+        dccd[ii].resize(3);
+        deed[ii].resize(3);
+        subForce[ii].resize(3);
+    }   
+    for( int ii = 0; ii < 3; ii++ ){
+       dccd[A][ii] = (deltaR[AD][ii]*rCD2 - deltaR[CD][ii]*dot)*term;
+       dccd[C][ii] = (deltaR[CD][ii]*rAD2 - deltaR[AD][ii]*dot)*term;
+       dccd[D][ii] = -one*(dccd[A][ii] + dccd[C][ii]);
+    }
+ 
+    AmoebaReferenceForce::getCrossProduct( deltaR[DB], deltaR[CB], deed[A] );
+    AmoebaReferenceForce::getCrossProduct( deltaR[AB], deltaR[DB], deed[C] );
+    AmoebaReferenceForce::getCrossProduct( deltaR[CB], deltaR[AB], deed[D] );
+ 
+    term        = eE/rDB2;
+    deed[D][0] += deltaR[DB][0]*term;
+    deed[D][1] += deltaR[DB][1]*term;
+    deed[D][2] += deltaR[DB][2]*term;
+ 
+    // ---------------------------------------------------------------------------------------
+ 
+    // forces
+ 
+    // calculate forces for atoms a, c, d
+    // the force for b is then -( a+ c + d)
+ 
+ 
+    for( int jj = 0; jj < LastAtomIndex; jj++ ){
+ 
+       // A, C, D
+ 
+       for( int ii = 0; ii < 3; ii++ ){
+          subForce[jj][ii] = dEdCos*( dccd[jj][ii] + deed[jj][ii] );
+       }
+ 
+       if( jj == 0 )jj++; // skip B
+ 
+       // now compute B
+ 
+       if( jj == 3 ){
+          for( int ii = 0; ii < 3; ii++ ){
+             subForce[1][ii] = -one*(subForce[0][ii] + subForce[2][ii] + subForce[3][ii]);
+          }
+       }
+    }
+ 
+    // add in forces
+ 
+    for( int jj = 0; jj < LastAtomIndex; jj++ ){
+       for( int ii = 0; ii < 3; ii++ ){
+          forces[jj][ii] -= subForce[jj][ii];
+       }
+    }
+ 
+    // ---------------------------------------------------------------------------------------
+ 
+    // calculate energy if 'energy' is set
+ 
+    RealOpenMM energy   = one + angleCubic*dt + angleQuartic*dt2 + anglePentic*dt3 + angleSextic*dt4;
+    energy             *= angleK*dt2;
+
+    return energy;
+}

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceOutOfPlaneBendForce.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 __AmoebaReferenceOutOfPlaneBendForce_H__
+#define __AmoebaReferenceOutOfPlaneBendForce_H__
+
+#include "SimTKUtilities/SimTKOpenMMRealType.h"
+
+// ---------------------------------------------------------------------------------------
+
+class AmoebaReferenceOutOfPlaneBendForce {
+
+public:
+ 
+    /**---------------------------------------------------------------------------------------
+       
+       Constructor
+       
+       --------------------------------------------------------------------------------------- */
+ 
+    AmoebaReferenceOutOfPlaneBendForce( );
+ 
+    /**---------------------------------------------------------------------------------------
+       
+       Destructor
+       
+          --------------------------------------------------------------------------------------- */
+ 
+    ~AmoebaReferenceOutOfPlaneBendForce( );
+ 
+    /**---------------------------------------------------------------------------------------
+    
+       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 angleK                  quadratic angle force parameter
+       @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( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
+                                               const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD,
+                                               RealOpenMM angleK,
+                                               RealOpenMM angleCubic,     RealOpenMM angleQuartic,
+                                               RealOpenMM anglePentic,    RealOpenMM angleSextic,
+                                               RealOpenMM** forces );
+         
+};
+
+// ---------------------------------------------------------------------------------------
+
+#endif // _AmoebaReferenceOutOfPlaneBendForce___