Reference code for Urey-Bradley

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

@@ -39,6 +39,7 @@ extern "C" void OPENMM_EXPORT registerKernelFactories() {
         if( platform.getName() == "Reference" ){
              AmoebaReferenceKernelFactory* factory = new AmoebaReferenceKernelFactory();
              platform.registerKernelFactory(CalcAmoebaHarmonicBondForceKernel::Name(), factory);
+             platform.registerKernelFactory(CalcAmoebaUreyBradleyForceKernel::Name(), factory);
              platform.registerKernelFactory(CalcAmoebaHarmonicAngleForceKernel::Name(), factory);
              platform.registerKernelFactory(CalcAmoebaHarmonicInPlaneAngleForceKernel::Name(), factory);
              platform.registerKernelFactory(CalcAmoebaTorsionForceKernel::Name(), factory);
@@ -64,6 +65,9 @@ KernelImpl* AmoebaReferenceKernelFactory::createKernelImpl(std::string name, con
     if (name == CalcAmoebaHarmonicBondForceKernel::Name())
         return new ReferenceCalcAmoebaHarmonicBondForceKernel(name, platform, context.getSystem());
 
+    if (name == CalcAmoebaUreyBradleyForceKernel::Name())
+        return new ReferenceCalcAmoebaUreyBradleyForceKernel(name, platform, context.getSystem());
+
     if (name == CalcAmoebaHarmonicAngleForceKernel::Name())
         return new ReferenceCalcAmoebaHarmonicAngleForceKernel(name, platform, context.getSystem());
 

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

@@ -37,6 +37,7 @@
 #include "AmoebaReferenceVdwForce.h"
 #include "AmoebaReferenceWcaDispersionForce.h"
 #include "internal/AmoebaWcaDispersionForceImpl.h"
+#include "AmoebaReferenceUreyBradleyForce.h"
 #include "ReferencePlatform.h"
 #include "openmm/internal/ContextImpl.h"
 #include "AmoebaMultipoleForce.h"
@@ -107,6 +108,43 @@ double ReferenceCalcAmoebaHarmonicBondForceKernel::execute(ContextImpl& context,
     return static_cast<double>(energy);
 }
 
+// ***************************************************************************
+
+ReferenceCalcAmoebaUreyBradleyForceKernel::ReferenceCalcAmoebaUreyBradleyForceKernel(std::string name, const Platform& platform, System& system) : 
+                CalcAmoebaUreyBradleyForceKernel(name, platform), system(system) {
+}
+
+ReferenceCalcAmoebaUreyBradleyForceKernel::~ReferenceCalcAmoebaUreyBradleyForceKernel() {
+}
+
+void ReferenceCalcAmoebaUreyBradleyForceKernel::initialize(const System& system, const AmoebaUreyBradleyForce& force) {
+
+    numIxns = force.getNumInteractions();
+    for( int ii = 0; ii < numIxns; ii++) {
+
+        int particle1Index, particle2Index;
+        double lengthValue, kValue;
+        force.getUreyBradleyParameters(ii, particle1Index, particle2Index, lengthValue, kValue );
+
+        particle1.push_back( particle1Index ); 
+        particle2.push_back( particle2Index ); 
+        length.push_back(    static_cast<RealOpenMM>( lengthValue ) );
+        kQuadratic.push_back( static_cast<RealOpenMM>( kValue ) );
+    } 
+    globalUreyBradleyCubic   = static_cast<RealOpenMM>(force.getAmoebaGlobalUreyBradleyCubic());
+    globalUreyBradleyQuartic = static_cast<RealOpenMM>(force.getAmoebaGlobalUreyBradleyQuartic());
+}
+
+double ReferenceCalcAmoebaUreyBradleyForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    RealOpenMM** posData   = extractPositions(context);
+    RealOpenMM** forceData = extractForces(context);
+    AmoebaReferenceUreyBradleyForce amoebaReferenceUreyBradleyForce;
+    RealOpenMM energy      = amoebaReferenceUreyBradleyForce.calculateForceAndEnergy( numIxns, posData, particle1, particle2, length, kQuadratic,
+                                                                                      globalUreyBradleyCubic, globalUreyBradleyQuartic,
+                                                                                      forceData );
+    return static_cast<double>(energy);
+}
+
 ReferenceCalcAmoebaHarmonicAngleForceKernel::ReferenceCalcAmoebaHarmonicAngleForceKernel(std::string name, const Platform& platform, System& system) :
             CalcAmoebaHarmonicAngleForceKernel(name, platform), system(system) {
 }

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

@@ -69,6 +69,42 @@ private:
     System& system;
 };
 
+/**
+ * This kernel is invoked by AmoebaUreyBradleyForce to calculate the forces acting on the system and the energy of the system.
+ */
+class ReferenceCalcAmoebaUreyBradleyForceKernel : public CalcAmoebaUreyBradleyForceKernel {
+public:
+    ReferenceCalcAmoebaUreyBradleyForceKernel(std::string name, 
+                                               const Platform& platform,
+                                               System& system);
+    ~ReferenceCalcAmoebaUreyBradleyForceKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param force      the AmoebaUreyBradleyForce this kernel will be used for
+     */
+    void initialize(const System& system, const AmoebaUreyBradleyForce& 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 numIxns;
+    std::vector<int>   particle1;
+    std::vector<int>   particle2;
+    std::vector<RealOpenMM> length;
+    std::vector<RealOpenMM> kQuadratic;
+    RealOpenMM globalUreyBradleyCubic;
+    RealOpenMM globalUreyBradleyQuartic;
+    System& system;
+};
+
 /**
  * This kernel is invoked by AmoebaHarmonicAngleForce to calculate the forces acting on the system and the energy of the system.
  */

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceUreyBradleyForce.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 "AmoebaReferenceUreyBradleyForce.h"
+#include "AmoebaReferenceForce.h"
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate Amoeba UB ixn (force and energy)
+
+   @param positionAtomA           Cartesian coordinates of atom A
+   @param positionAtomB           Cartesian coordinates of atom B
+   @param idealLength             ideal length
+   @param kForceConstant          force constant
+   @param cubicK                  cubic force parameter
+   @param quarticK                quartic force parameter
+   @param forces                  force vector
+
+   @return energy
+
+   --------------------------------------------------------------------------------------- */
+
+RealOpenMM AmoebaReferenceUreyBradleyForce::calculateUreyBradleyIxn( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
+                                                                     RealOpenMM idealLength, RealOpenMM kForceConstant,
+                                                                     RealOpenMM cubicK, RealOpenMM quarticK,
+                                                                     RealOpenMM** forces ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName = "AmoebaReferenceUreyBradleyForce::calculateHarmonicForce";
+
+   static const RealOpenMM zero          = 0.0;
+   static const RealOpenMM one           = 1.0;
+   static const RealOpenMM onePt5        = 1.5;
+   static const RealOpenMM two           = 2.0;
+
+   // ---------------------------------------------------------------------------------------
+
+   // get deltaR, R2, and R between 2 atoms
+
+   std::vector<RealOpenMM> deltaR;
+   AmoebaReferenceForce::loadDeltaR( positionAtomA, positionAtomB, deltaR );
+   RealOpenMM r               = AmoebaReferenceForce::getNorm3( deltaR );
+
+   // deltaIdeal = r - r_0
+
+   RealOpenMM deltaIdeal      = r - idealLength;
+   RealOpenMM deltaIdeal2     = deltaIdeal*deltaIdeal;
+
+   RealOpenMM dEdR            = (one + onePt5*cubicK*deltaIdeal + two*quarticK*deltaIdeal2);
+   dEdR                      *= two*kForceConstant*deltaIdeal;
+   dEdR                       = r > zero ? (dEdR/r) : zero;
+
+   forces[0][0]              += dEdR*deltaR[0];
+   forces[0][1]              += dEdR*deltaR[1];
+   forces[0][2]              += dEdR*deltaR[2];
+
+   forces[1][0]              -= dEdR*deltaR[0];
+   forces[1][1]              -= dEdR*deltaR[1];
+   forces[1][2]              -= dEdR*deltaR[2];
+
+   RealOpenMM energy          = kForceConstant*deltaIdeal2*( one + cubicK*deltaIdeal + quarticK*deltaIdeal2 );
+   return energy;
+}
+
+RealOpenMM AmoebaReferenceUreyBradleyForce::calculateForceAndEnergy( int numIxns,
+                                                                      RealOpenMM** particlePositions,
+                                                                      const std::vector<int>&   particle1,
+                                                                      const std::vector<int>&   particle2,
+                                                                      const std::vector<RealOpenMM>& length,
+                                                                      const std::vector<RealOpenMM>& kQuadratic,
+                                                                      RealOpenMM globalUreyBradleyCubic,
+                                                                      RealOpenMM globalUreyBradleyQuartic,
+                                                                      RealOpenMM** forceData ) const {
+    RealOpenMM energy      = 0.0; 
+    for( int ii = 0; ii < numIxns; ii++ ){
+        int particle1Index            = particle1[ii];
+        int particle2Index            = particle2[ii];
+        RealOpenMM idealLength        = length[ii];
+        RealOpenMM kForceConstant     = kQuadratic[ii];
+        RealOpenMM* forces[2];
+        forces[0]                     = forceData[particle1Index];
+        forces[1]                     = forceData[particle2Index];
+        energy                       += calculateUreyBradleyIxn( particlePositions[particle1Index], particlePositions[particle2Index],
+                                                                 idealLength, kForceConstant, globalUreyBradleyCubic, globalUreyBradleyQuartic,
+                                                                 forces );
+    }   
+    return energy;
+}

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceUreyBradleyForce.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 __AmoebaReferenceUreyBradleyForce_H__
+#define __AmoebaReferenceUreyBradleyForce_H__
+
+#include "SimTKUtilities/SimTKOpenMMRealType.h"
+#include <vector>
+
+// ---------------------------------------------------------------------------------------
+
+class AmoebaReferenceUreyBradleyForce {
+
+public:
+ 
+    /**---------------------------------------------------------------------------------------
+       
+        Constructor
+       
+        --------------------------------------------------------------------------------------- */
+ 
+    AmoebaReferenceUreyBradleyForce( ){};
+ 
+    /**---------------------------------------------------------------------------------------
+       
+        Destructor
+       
+        --------------------------------------------------------------------------------------- */
+ 
+    ~AmoebaReferenceUreyBradleyForce( ){};
+ 
+     /**---------------------------------------------------------------------------------------
+     
+        Calculate Amoeba Urey-Bradley ixns (force and energy)
+     
+        @param numIxns                 number of ixns
+        @param posData                 particle positions
+        @param particle1               particle 1 indices
+        @param particle2               particle 2 indices
+        @param length                  ideal length
+        @param kForce                  force constant
+        @param cubic                   cubic force parameter
+        @param quartic                 quartic force parameter
+        @param forceDate               output force vector
+     
+        @return total energy
+     
+        --------------------------------------------------------------------------------------- */
+     
+    RealOpenMM calculateForceAndEnergy( int numIxns, RealOpenMM** posData,
+                                        const std::vector<int>& particle1,
+                                        const std::vector<int>&  particle2,
+                                        const std::vector<RealOpenMM>& length,
+                                        const std::vector<RealOpenMM>& kForce,
+                                        RealOpenMM cubic, RealOpenMM quartic,
+                                        RealOpenMM** forceData ) const;
+
+private:
+
+     /**---------------------------------------------------------------------------------------
+     
+        Calculate Amoeba Urey-Bradley ixns (force and energy)
+     
+        @param positionAtomA           Cartesian coordinates of atom A
+        @param positionAtomB           Cartesian coordinates of atom B
+        @param length                  ideal length
+        @param kForce                  force constant
+        @param cubic                   cubic force parameter
+        @param quartic                 quartic force parameter
+        @param forces                  force vector
+     
+        @return energy
+     
+        --------------------------------------------------------------------------------------- */
+     
+    RealOpenMM calculateUreyBradleyIxn( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
+                                        RealOpenMM idealLength, RealOpenMM kForce,
+                                        RealOpenMM cubic, RealOpenMM quartic,
+                                        RealOpenMM** forces ) const;
+     
+};
+
+// ---------------------------------------------------------------------------------------
+
+#endif // _AmoebaReferenceUreyBradleyForce___

plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaUreyBradleyForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the Cuda implementation of UreyBradleyForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+//#include "AmoebaTinkerParameterFile.h"
+#include "openmm/Context.h"
+#include "AmoebaOpenMM.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+
+const double TOL = 1e-5;
+
+static void computeAmoebaUreyBradleyForce(int bondIndex,  std::vector<Vec3>& positions, AmoebaUreyBradleyForce& amoebaUreyBradleyForce,
+                                           std::vector<Vec3>& forces, double* energy ) {
+
+    int particle1, particle2;
+    double bondLength;
+    double quadraticK;
+    double cubicK    = amoebaUreyBradleyForce.getAmoebaGlobalUreyBradleyCubic();
+    double quarticK  = amoebaUreyBradleyForce.getAmoebaGlobalUreyBradleyQuartic();
+    amoebaUreyBradleyForce.getUreyBradleyParameters(bondIndex, particle1, particle2,  bondLength,  quadraticK );
+
+    double deltaR[3];
+    double r2 = 0.0;
+    for( int ii = 0; ii < 3; ii++ ){
+           deltaR[ii]    = positions[particle2][ii] - positions[particle1][ii];
+           r2           += deltaR[ii]*deltaR[ii];
+    }
+    double r                   = sqrt( r2 );
+
+    double bondDelta           = (r - bondLength);
+    double bondDelta2          = bondDelta*bondDelta;
+    double dEdR                = 1.0 + 1.5*cubicK*bondDelta + 2.0*quarticK*bondDelta2;
+
+           dEdR               *= (r > 0.0) ? (2.0*quadraticK*bondDelta)/r : 0.0;
+
+   forces[particle1][0]       += dEdR*deltaR[0];
+   forces[particle1][1]       += dEdR*deltaR[1];
+   forces[particle1][2]       += dEdR*deltaR[2];
+
+   forces[particle2][0]       -= dEdR*deltaR[0];
+   forces[particle2][1]       -= dEdR*deltaR[1];
+   forces[particle2][2]       -= dEdR*deltaR[2];
+
+   *energy                    += (1.0f + cubicK*bondDelta + quarticK*bondDelta2)*quadraticK*bondDelta2;
+
+}
+
+static void computeAmoebaUreyBradleyForces( Context& context, AmoebaUreyBradleyForce& amoebaUreyBradleyForce,
+                                             std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) {
+
+    // get positions and zero forces
+
+    State state = context.getState(State::Positions);
+    std::vector<Vec3> positions = state.getPositions();
+    expectedForces.resize( positions.size() );
+    
+    for( unsigned int ii = 0; ii < expectedForces.size(); ii++ ){
+        expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
+    }
+
+    // calculates forces/energy
+
+    *expectedEnergy = 0.0;
+    for( int ii = 0; ii < amoebaUreyBradleyForce.getNumInteractions(); ii++ ){
+        computeAmoebaUreyBradleyForce(ii, positions, amoebaUreyBradleyForce, expectedForces, expectedEnergy );
+    }
+
+    if( log ){
+        (void) fprintf( log, "computeAmoebaUreyBradleyForces: expected energy=%15.7e\n", *expectedEnergy );
+        for( unsigned int ii = 0; ii < positions.size(); ii++ ){
+            (void) fprintf( log, "%6u [%15.7e %15.7e %15.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2] );
+        }
+        (void) fflush( log );
+    }
+    return;
+
+}
+
+void compareWithExpectedForceAndEnergy( Context& context, AmoebaUreyBradleyForce& amoebaUreyBradleyForce, double tolerance, const std::string& idString, FILE* log) {
+
+    std::vector<Vec3> expectedForces;
+    double expectedEnergy;
+    computeAmoebaUreyBradleyForces( context, amoebaUreyBradleyForce, expectedForces, &expectedEnergy, NULL );
+   
+    State state                      = context.getState(State::Forces | State::Energy);
+    const std::vector<Vec3> forces   = state.getForces();
+
+    if( log ){
+        (void) fprintf( log, "computeAmoebaUreyBradleyForces: expected energy=%15.7e %15.7e\n", expectedEnergy, state.getPotentialEnergy() );
+        for( unsigned int ii = 0; ii < forces.size(); ii++ ){
+            (void) fprintf( log, "%6u [%15.7e %15.7e %15.7e]   [%15.7e %15.7e %15.7e]\n", ii,
+                            expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
+        }
+        (void) fflush( log );
+    }
+
+    for( unsigned int ii = 0; ii < forces.size(); ii++ ){
+        ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
+    }
+    ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance );
+}
+
+void testOneBond( FILE* log ) {
+
+    System system;
+
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+
+    AmoebaUreyBradleyForce* amoebaUreyBradleyForce = new AmoebaUreyBradleyForce();
+
+    double bondLength = 1.5;
+    double quadraticK = 1.0;
+    double cubicK     = 2.0;
+    double quarticicK = 3.0;
+    amoebaUreyBradleyForce->setAmoebaGlobalUreyBradleyCubic( cubicK );
+    amoebaUreyBradleyForce->setAmoebaGlobalUreyBradleyQuartic( quarticicK );
+    amoebaUreyBradleyForce->addUreyBradley(0, 1, bondLength, quadraticK);
+
+    system.addForce(amoebaUreyBradleyForce);
+    Context context(system, integrator, Platform::getPlatformByName( "Reference"));
+    std::vector<Vec3> positions(2);
+
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+
+    context.setPositions(positions);
+    compareWithExpectedForceAndEnergy( context, *amoebaUreyBradleyForce, TOL, "testOneBond", log );
+}
+
+void testTwoBond( FILE* log ) {
+
+    System system;
+
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+
+    AmoebaUreyBradleyForce* amoebaUreyBradleyForce = new AmoebaUreyBradleyForce();
+
+    double bondLength = 1.5;
+    double quadraticK = 1.0;
+    double cubicK     = 2.0;
+    double quarticicK = 3.0;
+    amoebaUreyBradleyForce->setAmoebaGlobalUreyBradleyCubic( cubicK );
+    amoebaUreyBradleyForce->setAmoebaGlobalUreyBradleyQuartic( quarticicK );
+    amoebaUreyBradleyForce->addUreyBradley(0, 1, bondLength, quadraticK);
+    amoebaUreyBradleyForce->addUreyBradley(1, 2, bondLength, quadraticK);
+
+    system.addForce(amoebaUreyBradleyForce);
+    Context context(system, integrator, Platform::getPlatformByName( "Reference"));
+    std::vector<Vec3> positions(3);
+
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 1);
+
+    context.setPositions(positions);
+    compareWithExpectedForceAndEnergy( context, *amoebaUreyBradleyForce, TOL, "testTwoBond", log );
+}
+
+int main( int numberOfArguments, char* argv[] ) {
+
+    try {
+        std::cout << "TestCudaAmoebaUreyBradleyForce running test..." << std::endl;
+        Platform::loadPluginsFromDirectory( Platform::getDefaultPluginsDirectory() );
+        FILE* log = NULL;
+        //FILE* log = stderr;
+
+        //testOneBond( log );
+        testTwoBond( log );
+
+        if( log && log != stderr )
+            (void) fclose( log );
+
+    }
+    catch(const std::exception& e) {
+        std::cout << "exception: " << e.what() << std::endl;
+        std::cout << "FAIL - ERROR.  Test failed." << std::endl;
+        return 1;
+    }
+    std::cout << "PASS - Test succeeded." << std::endl;
+    return 0;
+}