Added Urey-Bradley term

plugins/amoeba/openmmapi/include/AmoebaOpenMM.h

@@ -33,6 +33,7 @@
  * -------------------------------------------------------------------------- */
 
 #include "AmoebaHarmonicBondForce.h"
+#include "AmoebaUreyBradleyForce.h"
 #include "AmoebaHarmonicAngleForce.h"
 #include "AmoebaHarmonicInPlaneAngleForce.h"
 #include "AmoebaTorsionForce.h"

plugins/amoeba/openmmapi/include/AmoebaUreyBradleyForce.h

+#ifndef OPENMM_AMOEBA_UREY_BRADLEY_FORCE_H_
+#define OPENMM_AMOEBA_UREY_BRADLEY_FORCE_H_
+
+/* -------------------------------------------------------------------------- *
+ *                              AmoebaOpenMM                                  *
+ * -------------------------------------------------------------------------- *
+ * 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-2009 Stanford University and the Authors.      *
+ * Authors:                                                                   *
+ * 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/Force.h"
+#include "openmm/Vec3.h"
+#include <map>
+#include <vector>
+#include "openmm/internal/windowsExport.h"
+
+namespace OpenMM {
+
+/**
+ * This class implements an interaction between pairs of particles that varies harmonically with the distance
+ * between them.  To use it, create a AmoebaUreyBradleyForce object then call addUreyBradley() once for each bond.  After
+ * a bond has been added, you can modify its force field parameters by calling setUreyBradleyParameters().
+ */
+
+ class OPENMM_EXPORT AmoebaUreyBradleyForce : public Force {
+public:
+    /**
+     * Create a Amoeba UreyBradleyForce.
+     */
+    AmoebaUreyBradleyForce();
+    /**
+     * Get the number of UB terms in the potential function
+     */
+    int getNumInteractions() const {
+        return bonds.size();
+    }
+
+    /**
+     * Set the global cubic term
+     * 
+     * @param cubicK      the cubic force constant
+     */
+    void setAmoebaGlobalUreyBradleyCubic( double cubicK );
+
+    /**
+     * Get the global cubic term
+     * 
+     * @return global cubicK term
+     */
+    double getAmoebaGlobalUreyBradleyCubic( void ) const;
+
+    /**
+     * Set the global cubic term
+     * 
+     * @param quarticK  the quartic force constant
+     */
+    void setAmoebaGlobalUreyBradleyQuartic( double quarticK );
+
+    /**
+     * Get the global quartic term
+     * 
+     * @return global  quartic term
+     */
+    double getAmoebaGlobalUreyBradleyQuartic( void ) const;
+
+    /**
+     * Add a UB term to the force field.
+     *
+     * @param particle1     the index of the first particle 
+     * @param particle2     the index of the second particle
+     * @param length        the equilibrium length, measured in nm
+     * @param k             the quadratic harmonic force constant 
+     * @return the index of the bond that was added
+     */
+
+    int addUreyBradley(int particle1, int particle2, double length, double quadraticK );
+
+    /**
+     * Get the force field parameters for a bond term.
+     * 
+     * @param index         the index of the ixn for which to get parameters
+     * @param particle1     the index of the first particle
+     * @param particle2     the index of the second particle
+     * @param length        the equilibrium distance, measured in nm
+     * @param quadratic k   the quadratic harmonic force constant
+     */
+
+    void getUreyBradleyParameters(int index, int& particle1, int& particle2, double& length, double& quadraticK ) const;
+
+    /**
+     * Set the force field parameters for a UB term.
+     * 
+     * @param index     the index of the ixn for which to set parameters
+     * @param particle1 the index of the first particle
+     * @param particle2 the index of the second particle
+     * @param length    the equilibrium distance, measured in nm
+     * @param k         the quadratic harmonic force constant for the bond
+     */
+    void setUreyBradleyParameters(int index, int particle1, int particle2, double length, double quadraticK );
+
+protected:
+    double _globalQuarticK, _globalCubicK;
+    ForceImpl* createImpl();
+private:
+
+    class UreyBradleyInfo;
+
+// Retarded visual studio compiler complains about being unable to 
+// export private stl class members.
+// This stanza explains that it should temporarily shut up.
+#if defined(_MSC_VER)
+#pragma warning(push)
+#pragma warning(disable:4251)
+#endif
+    std::vector<UreyBradleyInfo> bonds;
+#if defined(_MSC_VER)
+#pragma warning(pop)
+#endif
+
+};
+
+class AmoebaUreyBradleyForce::UreyBradleyInfo {
+public:
+    int particle1, particle2;
+    double length, quadraticK;
+    UreyBradleyInfo() {
+        particle1 = particle2    = -1;
+        length    = quadraticK   = 0.0;
+    }
+    UreyBradleyInfo(int particle1, int particle2, double length, double  quadraticK ) :
+        particle1(particle1), particle2(particle2), length(length), quadraticK(quadraticK) {
+    }
+};
+
+} // namespace AmoebaOpenMM
+
+#endif /*OPENMM_AMOEBA_UREY_BRADLEY_FORCE_H_*/

plugins/amoeba/openmmapi/include/amoebaKernels.h

@@ -71,6 +71,32 @@ public:
     virtual double execute(ContextImpl& context, bool includeForces, bool includeEnergy) = 0;
 };
 
+class CalcAmoebaUreyBradleyForceKernel : public KernelImpl {
+public:
+    static std::string Name() {
+        return "CalcAmoebaUreyBradleyForce";
+    }
+
+    CalcAmoebaUreyBradleyForceKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
+    }
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param force      the AmoebaUreyBradleyForce this kernel will be used for
+     */
+    virtual void initialize(const System& system, const AmoebaUreyBradleyForce& force) = 0;
+    /**
+     * 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
+     */
+    virtual double execute(ContextImpl& context, bool includeForces, bool includeEnergy) = 0;
+};
+
 /**
  * This kernel is invoked by AmoebaHarmonicAngleForce to calculate the forces acting on the system and the energy of the system.
  */

plugins/amoeba/openmmapi/include/internal/AmoebaUreyBradleyForceImpl.h

+#ifndef OPENMM_AMOEBA_UREY_BRADLEY_FORCE_IMPL_H_
+#define OPENMM_AMOEBA_UREY_BRADLEY_FORCE_IMPL_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                AmoebaOpenMM                                *
+ * -------------------------------------------------------------------------- *
+ * 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:                                                                   *
+ * 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/internal/ForceImpl.h"
+#include "AmoebaUreyBradleyForce.h"
+#include "openmm/Kernel.h"
+#include <utility>
+#include <set>
+#include <string>
+
+namespace OpenMM {
+
+/**
+ * This is the internal implementation of AmoebaUreyBradleyForce.
+ */
+
+class AmoebaUreyBradleyForceImpl : public ForceImpl {
+public:
+    AmoebaUreyBradleyForceImpl(AmoebaUreyBradleyForce& owner);
+    ~AmoebaUreyBradleyForceImpl();
+    void initialize(ContextImpl& context);
+    AmoebaUreyBradleyForce& getOwner() {
+        return owner;
+    }
+    void updateContextState(ContextImpl& context) {
+        // This force field doesn't update the state directly.
+    }
+    double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy);
+    std::map<std::string, double> getDefaultParameters() {
+        return std::map<std::string, double>(); // This force field doesn't define any parameters.
+    }
+    std::vector<std::string> getKernelNames();
+private:
+    AmoebaUreyBradleyForce& owner;
+    Kernel kernel;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_AMOEBA_UREY_BRADLEY_FORCE_IMPL_H_*/

plugins/amoeba/openmmapi/src/AmoebaUreyBradleyForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                AmoebaOpenMM                                *
+ * -------------------------------------------------------------------------- *
+ * 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-2009 Stanford University and the Authors.      *
+ * Authors:                                                                   *
+ * 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/Force.h"
+#include "openmm/OpenMMException.h"
+#include "AmoebaUreyBradleyForce.h"
+#include "internal/AmoebaUreyBradleyForceImpl.h"
+
+using namespace OpenMM;
+
+AmoebaUreyBradleyForce::AmoebaUreyBradleyForce() {
+   _globalCubicK = _globalQuarticK = 0.0;
+}
+
+int AmoebaUreyBradleyForce::addUreyBradley(int particle1, int particle2, double length, double quadraticK) {
+    bonds.push_back(UreyBradleyInfo(particle1, particle2, length, quadraticK ));
+    return bonds.size()-1;
+}
+
+void AmoebaUreyBradleyForce::getUreyBradleyParameters(int index, int& particle1, int& particle2, double& length, double&  quadraticK ) const {
+    particle1       = bonds[index].particle1;
+    particle2       = bonds[index].particle2;
+    length          = bonds[index].length;
+    quadraticK      = bonds[index].quadraticK;
+}
+
+void AmoebaUreyBradleyForce::setUreyBradleyParameters(int index, int particle1, int particle2, double length, double quadraticK ) {
+    bonds[index].particle1  = particle1;
+    bonds[index].particle2  = particle2;
+    bonds[index].length     = length;
+    bonds[index].quadraticK = quadraticK;
+}
+
+void AmoebaUreyBradleyForce::setAmoebaGlobalUreyBradleyCubic(double cubicK ) {
+    _globalCubicK           = cubicK;
+}
+
+void AmoebaUreyBradleyForce::setAmoebaGlobalUreyBradleyQuartic(double quarticK ) {
+    _globalQuarticK         = quarticK;
+}
+
+double AmoebaUreyBradleyForce::getAmoebaGlobalUreyBradleyCubic( void ) const {
+    return _globalCubicK;
+}
+
+double AmoebaUreyBradleyForce::getAmoebaGlobalUreyBradleyQuartic( void ) const {
+    return _globalQuarticK;
+}
+
+ForceImpl* AmoebaUreyBradleyForce::createImpl() {
+    return new AmoebaUreyBradleyForceImpl(*this);
+}

plugins/amoeba/openmmapi/src/AmoebaUreyBradleyForceImpl.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:                                                                   *
+ * 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/internal/ContextImpl.h"
+#include "internal/AmoebaUreyBradleyForceImpl.h"
+#include "amoebaKernels.h"
+#include "openmm/Platform.h"
+
+using namespace OpenMM;
+
+using std::pair;
+using std::vector;
+using std::set;
+
+AmoebaUreyBradleyForceImpl::AmoebaUreyBradleyForceImpl(AmoebaUreyBradleyForce& owner) : owner(owner) {
+}
+
+AmoebaUreyBradleyForceImpl::~AmoebaUreyBradleyForceImpl() {
+}
+
+void AmoebaUreyBradleyForceImpl::initialize(ContextImpl& context) {
+    kernel = context.getPlatform().createKernel(CalcAmoebaUreyBradleyForceKernel::Name(), context);
+    dynamic_cast<CalcAmoebaUreyBradleyForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
+}
+
+double AmoebaUreyBradleyForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    return dynamic_cast<CalcAmoebaUreyBradleyForceKernel&>(kernel.getImpl()).execute(context, includeForces, includeEnergy);
+}
+
+std::vector<std::string> AmoebaUreyBradleyForceImpl::getKernelNames() {
+    std::vector<std::string> names;
+    names.push_back(CalcAmoebaUreyBradleyForceKernel::Name());
+    return names;
+}

plugins/amoeba/platforms/cuda/src/AmoebaCudaKernelFactory.cpp

@@ -42,6 +42,7 @@ extern "C" void OPENMMCUDA_EXPORT registerKernelFactories() {
              AmoebaCudaKernelFactory* factory = new AmoebaCudaKernelFactory();
 
              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);
@@ -131,5 +132,8 @@ KernelImpl* AmoebaCudaKernelFactory::createKernelImpl(std::string name, const Pl
     if (name == CalcAmoebaWcaDispersionForceKernel::Name())
         return new CudaCalcAmoebaWcaDispersionForceKernel(name, platform, *amoebaCudaData, context.getSystem());
 
+    if (name == CalcAmoebaUreyBradleyForceKernel::Name())
+        return new CudaCalcAmoebaUreyBradleyForceKernel(name, platform, *amoebaCudaData, context.getSystem());
+
     throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '")+name+"'").c_str());
 }

plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp

@@ -133,6 +133,79 @@ double CudaCalcAmoebaHarmonicBondForceKernel::execute(ContextImpl& context, bool
     return 0.0;
 }
 
+/* -------------------------------------------------------------------------- *
+ *                           AmoebaUreyBradley                               *
+ * -------------------------------------------------------------------------- */
+
+class CudaCalcAmoebaUreyBradleyForceKernel::ForceInfo : public CudaForceInfo {
+public:
+    ForceInfo(const AmoebaUreyBradleyForce& force) : force(force) {
+    }
+    int getNumParticleGroups() {
+        return force.getNumInteractions();
+    }
+    void getParticlesInGroup(int index, std::vector<int>& particles) {
+        int particle1, particle2;
+        double length, k;
+        force.getUreyBradleyParameters(index, particle1, particle2, length, k);
+        particles.resize(2);
+        particles[0] = particle1;
+        particles[1] = particle2;
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        int particle1, particle2;
+        double length1, length2, k1, k2;
+        force.getUreyBradleyParameters(group1, particle1, particle2, length1, k1);
+        force.getUreyBradleyParameters(group2, particle1, particle2, length2, k2);
+        return (length1 == length2 && k1 == k2);
+    }
+private:
+    const AmoebaUreyBradleyForce& force;
+};
+
+CudaCalcAmoebaUreyBradleyForceKernel::CudaCalcAmoebaUreyBradleyForceKernel(std::string name, const Platform& platform, AmoebaCudaData& data, System& system) : 
+                CalcAmoebaUreyBradleyForceKernel(name, platform), data(data), system(system) {
+    data.incrementKernelCount();
+}
+
+CudaCalcAmoebaUreyBradleyForceKernel::~CudaCalcAmoebaUreyBradleyForceKernel() {
+    data.decrementKernelCount();
+}
+
+void CudaCalcAmoebaUreyBradleyForceKernel::initialize(const System& system, const AmoebaUreyBradleyForce& force) {
+
+    data.setAmoebaLocalForcesKernel( this );
+
+    numInteractions = force.getNumInteractions();
+    std::vector<int>   particle1(numInteractions);
+    std::vector<int>   particle2(numInteractions);
+    std::vector<float> length(numInteractions);
+    std::vector<float> quadratic(numInteractions);
+
+    for (int i = 0; i < numInteractions; i++) {
+
+        int particle1Index, particle2Index;
+        double lengthValue, kValue;
+        force.getUreyBradleyParameters(i, particle1Index, particle2Index, lengthValue, kValue );
+
+        particle1[i]     = particle1Index; 
+        particle2[i]     = particle2Index; 
+        length[i]        = static_cast<float>( lengthValue );
+        quadratic[i]     = static_cast<float>( kValue );
+    } 
+    gpuSetAmoebaUreyBradleyParameters( data.getAmoebaGpu(), particle1, particle2, length, quadratic, 
+                                       static_cast<float>(force.getAmoebaGlobalUreyBradleyCubic()),
+                                       static_cast<float>(force.getAmoebaGlobalUreyBradleyQuartic()) );
+    data.getAmoebaGpu()->gpuContext->forces.push_back(new ForceInfo(force));
+}
+
+double CudaCalcAmoebaUreyBradleyForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    if( data.getAmoebaLocalForcesKernel() == this ){
+        computeAmoebaLocalForces( data );
+    }
+    return 0.0;
+}
+
 /* -------------------------------------------------------------------------- *
  *                           AmoebaHarmonicAngle                              *
  * -------------------------------------------------------------------------- */

plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.h

@@ -66,6 +66,39 @@ private:
     System& system;
 };
 
+/**
+ * This kernel is invoked by AmoebaUreyBradleyForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CudaCalcAmoebaUreyBradleyForceKernel : public CalcAmoebaUreyBradleyForceKernel {
+public:
+    CudaCalcAmoebaUreyBradleyForceKernel(std::string name, 
+                                          const Platform& platform,
+                                          AmoebaCudaData& data,
+                                          System& system);
+    ~CudaCalcAmoebaUreyBradleyForceKernel();
+    /**
+     * 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:
+    class ForceInfo;
+    int numInteractions;
+    AmoebaCudaData& data;
+    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/cuda/src/kernels/amoebaCudaGpu.cpp

@@ -328,6 +328,16 @@ void gpuPrintCudaAmoebaGmxSimulation(amoebaGpuContext amoebaGpu, FILE* log )
 
     (void) fprintf( log, "     pOutputBufferCounter               %p\n",      amoebaGpu->gpuContext->pOutputBufferCounter );
 
+    if( amoebaGpu->psAmoebaUreyBradleyParameter)(void) fprintf( log, "\n" );
+    gpuPrintCudaStreamInt4( amoebaGpu->psAmoebaUreyBradleyID, log );
+    gpuPrintCudaStreamFloat2( amoebaGpu->psAmoebaUreyBradleyParameter, log );
+    (void) fprintf( log, "     amoebaUreyBradleys                 %u\n",      amoebaGpu->amoebaSim.amoebaUreyBradleys );
+    (void) fprintf( log, "     amoebaUreyBradley_offset           %u\n",      amoebaGpu->amoebaSim.amoebaUreyBradley_offset );
+    (void) fprintf( log, "     cubic                              %15.7e\n",  amoebaGpu->amoebaSim.amoebaUreyBradleyCubicParameter);
+    (void) fprintf( log, "     quartic                            %15.7e\n",  amoebaGpu->amoebaSim.amoebaUreyBradleyQuarticicParameter);
+    (void) fprintf( log, "     pAmoebaUreyBradleyID               %p\n",      amoebaGpu->amoebaSim.pAmoebaUreyBradleyID );
+    (void) fprintf( log, "     pAmoebaUreyBradleyParameter        %p\n",      amoebaGpu->amoebaSim.pAmoebaUreyBradleyParameter );
+    
     if( amoebaGpu->psRotationMatrix)(void) fprintf( log, "\n" );
     gpuPrintCudaStreamFloat( amoebaGpu->psRotationMatrix, log );
     gpuPrintCudaStreamInt4( amoebaGpu->psMultipoleParticlesIdsAndAxisType, log );
@@ -480,6 +490,50 @@ if( amoebaGpu->log && (i < DUMP_PARAMETERS || i > bonds - (DUMP_PARAMETERS + 1)
     psBondParameter->Upload();
 }
 
+extern "C"
+void gpuSetAmoebaUreyBradleyParameters(amoebaGpuContext amoebaGpu, const std::vector<int>& particles1, const std::vector<int>& particles2, 
+                                       const std::vector<float>& length, const std::vector<float>& k, float cubic, float quartic)
+{
+    _gpuContext* gpu                                         = amoebaGpu->gpuContext;
+    int bonds                                                = particles1.size();
+    amoebaGpu->amoebaSim.amoebaUreyBradleys                  = bonds;
+
+    CUDAStream<int4>* psUreyBradleyID                        = new CUDAStream<int4>(bonds, 1, "AmoebaUreyBradleyID");
+    amoebaGpu->psAmoebaUreyBradleyID                         = psUreyBradleyID;
+    amoebaGpu->amoebaSim.pAmoebaUreyBradleyID                = psUreyBradleyID->_pDevStream[0];
+
+    CUDAStream<float2>* psUreyBradleyParameter               = new CUDAStream<float2>(bonds, 1, "AmoebaUreyBradleyParameter");
+    amoebaGpu->psAmoebaUreyBradleyParameter                  = psUreyBradleyParameter;
+    amoebaGpu->amoebaSim.pAmoebaUreyBradleyParameter         = psUreyBradleyParameter->_pDevStream[0];
+
+    amoebaGpu->amoebaSim.amoebaUreyBradleyCubicParameter     = cubic;
+    amoebaGpu->amoebaSim.amoebaUreyBradleyQuarticicParameter = quartic;
+    for (int i = 0; i < bonds; i++)
+    {
+        (*psUreyBradleyID)[i].x         = particles1[i];
+        (*psUreyBradleyID)[i].y         = particles2[i];
+        (*psUreyBradleyID)[i].z         = gpu->pOutputBufferCounter[(*psUreyBradleyID)[i].x]++;
+        (*psUreyBradleyID)[i].w         = gpu->pOutputBufferCounter[(*psUreyBradleyID)[i].y]++;
+        (*psUreyBradleyParameter)[i].x  = length[i];
+        (*psUreyBradleyParameter)[i].y  = k[i];
+
+#undef DUMP_PARAMETERS
+#define DUMP_PARAMETERS 5
+#if (DUMP_PARAMETERS > 0 )                
+if( amoebaGpu->log && (i < DUMP_PARAMETERS || i > bonds - (DUMP_PARAMETERS + 1)  ) )
+        fprintf( amoebaGpu->log, "UreyBradleys: %5d [%5d %5d %5d %5d] L=%15.7e k[%15.7e %15.7e %15.7e] [%5d %5d]\n",
+            i, (*psUreyBradleyID)[i].x, (*psUreyBradleyID)[i].y, (*psUreyBradleyID)[i].z, (*psUreyBradleyID)[i].w, 
+            (*psUreyBradleyParameter)[i].x, (*psUreyBradleyParameter)[i].y, cubic, quartic, 
+            gpu->pOutputBufferCounter[(*psUreyBradleyID)[i].x],
+            gpu->pOutputBufferCounter[(*psUreyBradleyID)[i].y] );
+#endif
+#undef DUMP_PARAMETERS
+
+    }
+    psUreyBradleyID->Upload();
+    psUreyBradleyParameter->Upload();
+}
+
 extern "C"
 void gpuSetAmoebaAngleParameters(amoebaGpuContext amoebaGpu, const std::vector<int>& particles1, const std::vector<int>& particles2, const std::vector<int>& particles3,
                                  const std::vector<float>& angle, const std::vector<float>& k,
@@ -579,7 +633,7 @@ void gpuSetAmoebaInPlaneAngleParameters(amoebaGpuContext amoebaGpu, const std::v
         psAngleID2->_pSysData[i].w = gpu->pOutputBufferCounter[psAngleID1->_pSysData[i].w]++;
 
 #undef DUMP_PARAMETERS
-#define DUMP_PARAMETERS 50000
+#define DUMP_PARAMETERS 10
 #if (DUMP_PARAMETERS > 0 )
 if( (i < DUMP_PARAMETERS || i > bond_angles - (DUMP_PARAMETERS + 1)) && amoebaGpu->log )
          fprintf( amoebaGpu->log, "InPlaneAngles: %5d [%5d %5d %5d %5d] [%5d %5d %5d %5d] A=%15.7e k=%15.7e [%5d %5d %5d %5d]\n", i, 
@@ -822,7 +876,7 @@ void gpuSetAmoebaOutOfPlaneBendParameters(amoebaGpuContext amoebaGpu, const std:
     amoebaGpu->amoebaSim.amoebaOutOfPlaneBendSexticK     = sexticK;
 
 #undef DUMP_PARAMETERS
-#define DUMP_PARAMETERS 50000
+#define DUMP_PARAMETERS 10
 #if (DUMP_PARAMETERS > 0 )
     if( amoebaGpu->log )
         fprintf( amoebaGpu->log, "OutOfPlaneBends: global ks[%15.7e %15.7e %15.7e %15.7e]\n", cubicK, quarticK, penticK, sexticK );
@@ -1065,7 +1119,6 @@ void gpuSetAmoebaBondOffsets(amoebaGpuContext amoebaGpu )
     flipped = 1;
     _gpuContext* gpu                                           = amoebaGpu->gpuContext;
 
-
     amoebaGpu->amoebaSim.amoebaBond_offset                     = amoebaGpu->psAmoebaBondParameter ? amoebaGpu->psAmoebaBondParameter->_stride : 0;
 
     amoebaGpu->amoebaSim.amoebaAngle_offset                    = amoebaGpu->amoebaSim.amoebaBond_offset            + 
@@ -1089,7 +1142,12 @@ void gpuSetAmoebaBondOffsets(amoebaGpuContext amoebaGpu )
     amoebaGpu->amoebaSim.amoebaTorsionTorsion_offset           = amoebaGpu->amoebaSim.amoebaOutOfPlaneBend_offset  +
                                                                  (amoebaGpu->psAmoebaTorsionTorsionID1 ?  amoebaGpu->psAmoebaTorsionTorsionID1->_stride : 0);
 
+    amoebaGpu->amoebaSim.amoebaUreyBradley_offset              = amoebaGpu->amoebaSim.amoebaTorsionTorsion_offset +
+                                                                 (amoebaGpu->psAmoebaUreyBradleyParameter ?  amoebaGpu->psAmoebaUreyBradleyParameter->_stride : 0);
+
     //gpu->sim.localForces_threads_per_block  = (std::max(amoebaGpu->amoebaSim.amoebaTorsionTorsion_offset, gpu->sim.customBonds) / gpu->sim.blocks + 15) & 0xfffffff0;
+    //fprintf( amoebaGpu->log, "Final (UreyBradley) offset : %5d \n", amoebaGpu->amoebaSim.amoebaUreyBradley_offset );
+
     unsigned int maxI                                         = (amoebaGpu->amoebaSim.amoebaTorsionTorsion_offset > gpu->sim.customBonds) ? amoebaGpu->amoebaSim.amoebaTorsionTorsion_offset : gpu->sim.customBonds;
     gpu->sim.localForces_threads_per_block                    = (maxI/gpu->sim.blocks + 15) & 0xfffffff0;
     if (gpu->sim.localForces_threads_per_block > gpu->sim.max_localForces_threads_per_block)
@@ -1189,7 +1247,14 @@ void gpuSetAmoebaBondOffsets(amoebaGpuContext amoebaGpu )
             amoebaGpu->psAmoebaTorsionTorsionID2->Upload();
             amoebaGpu->psAmoebaTorsionTorsionID3->Upload();
         }
-
+        for (unsigned int i = 0; i < amoebaGpu->amoebaSim.amoebaUreyBradleys; i++) 
+        {    
+            (*amoebaGpu->psAmoebaUreyBradleyID)[i].z = flip - (*amoebaGpu->psAmoebaUreyBradleyID)[i].z;
+            (*amoebaGpu->psAmoebaUreyBradleyID)[i].w = flip - (*amoebaGpu->psAmoebaUreyBradleyID)[i].w;
+        }    
+        if( amoebaGpu->psAmoebaUreyBradleyID ){
+            amoebaGpu->psAmoebaUreyBradleyID->Upload();
+        }
     }    
     
 
@@ -2606,6 +2671,9 @@ void amoebaGpuShutDown(amoebaGpuContext gpu)
     delete gpu->psAmoebaBondID;
     delete gpu->psAmoebaBondParameter;
 
+    delete gpu->psAmoebaUreyBradleyID;
+    delete gpu->psAmoebaUreyBradleyParameter;
+
     delete gpu->psAmoebaAngleID1;
     delete gpu->psAmoebaAngleID2;
     delete gpu->psAmoebaAngleParameter;

plugins/amoeba/platforms/cuda/src/kernels/amoebaCudaTypes.h

@@ -59,6 +59,7 @@ struct cudaAmoebaGmxSimulation {
     float           amoebaBondQuarticicParameter;   // quartic bond parameters
     unsigned int    amoebaBond_offset;              // Offset to end of bonds
 
+
     unsigned int    amoebaAngles;                   // Number of bond angles
     int4*           pAmoebaAngleID1;                // Bond angle atom and first output buffer IDs
     int2*           pAmoebaAngleID2;                // Bond angle output buffer IDs
@@ -124,6 +125,13 @@ struct cudaAmoebaGmxSimulation {
     float amoebaTorTorGridDelta[4];                 // 15.0
     float4*          pAmoebaTorsionTorsionGrids;    // torsion torsion grids
 
+    unsigned int    amoebaUreyBradleys;             // Number of UB ixns
+    int4*           pAmoebaUreyBradleyID;           // UreyBradley atom and output buffer IDs
+    float2*         pAmoebaUreyBradleyParameter;    // UreyBradley parameters
+    float           amoebaUreyBradleyCubicParameter;// cubic parameter
+    float           amoebaUreyBradleyQuarticicParameter; // quartic parameter
+    unsigned int    amoebaUreyBradley_offset;       // Offset to end of bonds
+
     unsigned int numberOfAtoms;                     // number of atoms
     unsigned int paddedNumberOfAtoms;               // padded number of atoms
     //float cutoffDistance2;                          // cutoff distance squared for PME

plugins/amoeba/platforms/cuda/src/kernels/amoebaGpuTypes.h

@@ -86,6 +86,9 @@ struct _amoebaGpuContext {
     CUDAStream<int4>*   psAmoebaBondID;
     CUDAStream<float2>* psAmoebaBondParameter;
 
+    CUDAStream<int4>*   psAmoebaUreyBradleyID;
+    CUDAStream<float2>* psAmoebaUreyBradleyParameter;
+
     CUDAStream<int4>*   psAmoebaAngleID1;
     CUDAStream<int2>*   psAmoebaAngleID2;
     CUDAStream<float2>* psAmoebaAngleParameter;
@@ -249,6 +252,10 @@ extern "C"
 void gpuSetAmoebaBondParameters(amoebaGpuContext gpu, const std::vector<int>& atom1, const std::vector<int>& atom2, 
                                 const std::vector<float>& length, const std::vector<float>& k, float cubic, float quartic);
 
+extern "C"
+void gpuSetAmoebaUreyBradleyParameters(amoebaGpuContext gpu, const std::vector<int>& atom1, const std::vector<int>& atom2, 
+                                       const std::vector<float>& length, const std::vector<float>& k, float cubic, float quartic);
+
 extern "C"
 void gpuSetAmoebaAngleParameters(amoebaGpuContext gpu, const std::vector<int>& atom1, const std::vector<int>& atom2, const std::vector<int>& atom3,
                                  const std::vector<float>& angle, const std::vector<float>& k, float cubicK,

plugins/amoeba/platforms/cuda/src/kernels/kCalculateAmoebaCudaLocalForces.cu

@@ -1571,6 +1571,49 @@ void kCalculateAmoebaLocalForces_kernel()
         }
         pos += blockDim.x * gridDim.x;
     }
+
+    while (pos < cAmoebaSim.amoebaUreyBradley_offset)
+    {
+        unsigned int pos1   = pos - cAmoebaSim.amoebaTorsionTorsion_offset;
+        if (pos1 < cAmoebaSim.amoebaUreyBradleys)
+        {
+            int4   atom         = cAmoebaSim.pAmoebaUreyBradleyID[pos1];
+            float4 atomA        = cSim.pPosq[atom.x];
+            float4 atomB        = cSim.pPosq[atom.y];
+            float2 bond         = cAmoebaSim.pAmoebaUreyBradleyParameter[pos];
+            float dx            = atomB.x - atomA.x;
+            float dy            = atomB.y - atomA.y;
+            float dz            = atomB.z - atomA.z;
+            float r2            = dx * dx + dy * dy + dz * dz;
+            float r             = sqrt(r2);
+            float deltaIdeal    = r - bond.x;
+            float deltaIdeal2   = deltaIdeal*deltaIdeal;
+            energy             += bond.y * deltaIdeal2*( 1.0f + cAmoebaSim.amoebaUreyBradleyCubicParameter*deltaIdeal +
+                                                                cAmoebaSim.amoebaUreyBradleyQuarticicParameter*deltaIdeal2 );
+
+            float dEdR          = 2.0f*bond.y * deltaIdeal*( 1.0f + 1.5f*cAmoebaSim.amoebaUreyBradleyCubicParameter*deltaIdeal +
+                                                                    2.0f*cAmoebaSim.amoebaUreyBradleyQuarticicParameter*deltaIdeal2 );
+
+            dEdR                = (r > 0.0f) ? (dEdR / r) : 0.0f;
+
+            dx                 *= dEdR;
+            dy                 *= dEdR;
+            dz                 *= dEdR;
+            unsigned int offsetA                = atom.x + atom.z * cSim.stride;
+            unsigned int offsetB                = atom.y + atom.w * cSim.stride;
+            float4 forceA                       = cSim.pForce4[offsetA];
+            float4 forceB                       = cSim.pForce4[offsetB];
+            forceA.x                           += dx;
+            forceA.y                           += dy;
+            forceA.z                           += dz;
+            forceB.x                           -= dx;
+            forceB.y                           -= dy;
+            forceB.z                           -= dz;
+            cSim.pForce4[offsetA]               = forceA;
+            cSim.pForce4[offsetB]               = forceB;    
+        }
+        pos += blockDim.x * gridDim.x;
+    }
     cSim.pEnergy[blockIdx.x * blockDim.x + threadIdx.x] += energy;
 }
 
@@ -1579,7 +1622,7 @@ void kCalculateAmoebaLocalForces(amoebaGpuContext gpu)
 {
    
     if( gpu->log ){
-        static int call = 0;
+        static int call = 1;
         if( call == 0 ){
             (void) fprintf( gpu->log,"kCalculateAmoebaLocalForces: blks=%u thrds/blk=%u\n",
                             gpu->gpuContext->sim.blocks, gpu->gpuContext->sim.localForces_threads_per_block); fflush( gpu->log );
@@ -1588,17 +1631,8 @@ void kCalculateAmoebaLocalForces(amoebaGpuContext gpu)
     }
 
     kCalculateAmoebaLocalForces_kernel<<<gpu->gpuContext->sim.blocks, gpu->gpuContext->sim.localForces_threads_per_block>>>();
+
     LAUNCHERROR("kCalculateAmoebaLocalForces");
 
-/*    
-    fprintf( stderr, "Done kCalculateAmoebaLocalForces %p \n", gpu->gpuContext->psForce4); fflush( stderr );
-    gpu->gpuContext->psForce4->Print( stderr ); fflush( stderr );
-    gpu->gpuContext->psEnergy->Print( stderr ); fflush( stderr );
-    GetCalculateAmoebaLocalForcesSim( gpu );
-    gpu->gpuContext->psForce4->Print( stderr ); fflush( stderr );
-    gpu->gpuContext->psEnergy->Print( stderr ); fflush( stderr );
-    fprintf( stderr, "Ez %p\n", (float* ) gpu->gpuContext->sim.pEnergy );
-*/
-    
 }
 

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaUreyBradleyForce.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( "Cuda"));
+    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( "Cuda"));
+    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;
+}