Began converting AMOEBA to the new CUDA platform: AmoebaHarmonicBondForce,...

Began converting AMOEBA to the new CUDA platform: AmoebaHarmonicBondForce, AmoebaHarmonicAngleForce, AmoebaHarmonicInPlaneAngleForce, AmoebaPiTorsionForce

plugins/amoeba/platforms/cuda2/CMakeLists.txt

+#---------------------------------------------------
+# OpenMM CUDA Amoeba Implementation
+#
+# Creates OpenMM library, base name=OpenMMAmoebaCUDA.
+# Default libraries are shared & optimized. Variants
+# are created for debug (_d).
+#
+# Windows:
+#   OpenMMAmoebaCUDA[_d].dll
+#   OpenMMAmoebaCUDA[_d].lib
+# Unix:
+#   libOpenMMAmoebaCUDA[_d].so
+#----------------------------------------------------
+
+# The source is organized into subdirectories, but we handle them all from
+# this CMakeLists file rather than letting CMake visit them as SUBDIRS.
+SET(OPENMM_SOURCE_SUBDIRS .)
+
+
+# Collect up information about the version of the OpenMM library we're building
+# and make it available to the code so it can be built into the binaries.
+
+SET(OPENMMAMOEBACUDA_LIBRARY_NAME OpenMMAmoebaCUDA)
+
+SET(SHARED_TARGET ${OPENMMAMOEBACUDA_LIBRARY_NAME})
+
+
+# Ensure that debug libraries have "_d" appended to their names.
+# CMake gets this right on Windows automatically with this definition.
+IF (${CMAKE_GENERATOR} MATCHES "Visual Studio")
+    SET(CMAKE_DEBUG_POSTFIX "_d" CACHE INTERNAL "" FORCE)
+ENDIF (${CMAKE_GENERATOR} MATCHES "Visual Studio")
+
+# But on Unix or Cygwin we have to add the suffix manually
+IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
+    SET(SHARED_TARGET ${SHARED_TARGET}_d)
+ENDIF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
+
+
+# These are all the places to search for header files which are
+# to be part of the API.
+SET(API_INCLUDE_DIRS) # start empty
+FOREACH(subdir ${OPENMM_SOURCE_SUBDIRS})
+    # append
+    SET(API_INCLUDE_DIRS ${API_INCLUDE_DIRS}
+                         ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/include
+                         ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/include/internal)
+ENDFOREACH(subdir)
+
+# We'll need both *relative* path names, starting with their API_INCLUDE_DIRS,
+# and absolute pathnames.
+SET(API_REL_INCLUDE_FILES)   # start these out empty
+SET(API_ABS_INCLUDE_FILES)
+
+FOREACH(dir ${API_INCLUDE_DIRS})
+    FILE(GLOB fullpaths ${dir}/*.h)	# returns full pathnames
+    SET(API_ABS_INCLUDE_FILES ${API_ABS_INCLUDE_FILES} ${fullpaths})
+
+    FOREACH(pathname ${fullpaths})
+        GET_FILENAME_COMPONENT(filename ${pathname} NAME)
+        SET(API_REL_INCLUDE_FILES ${API_REL_INCLUDE_FILES} ${dir}/${filename})
+    ENDFOREACH(pathname)
+ENDFOREACH(dir)
+
+# collect up source files
+SET(SOURCE_FILES) # empty
+SET(SOURCE_INCLUDE_FILES)
+
+FOREACH(subdir ${OPENMM_SOURCE_SUBDIRS})
+    FILE(GLOB_RECURSE src_files  ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/*.cpp ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/*.c)
+    FILE(GLOB incl_files ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/*.h)
+    SET(SOURCE_FILES         ${SOURCE_FILES}         ${src_files})   #append
+    SET(SOURCE_INCLUDE_FILES ${SOURCE_INCLUDE_FILES} ${incl_files})
+    INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/include)
+ENDFOREACH(subdir)
+
+INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/src)
+INCLUDE_DIRECTORIES(BEFORE ${CMAKE_SOURCE_DIR}/platforms/cuda2/include)
+INCLUDE_DIRECTORIES(BEFORE ${CMAKE_SOURCE_DIR}/platforms/cuda2/src)
+INCLUDE_DIRECTORIES(BEFORE ${CMAKE_BINARY_DIR}/platforms/cuda2/src)
+
+# Set variables needed for encoding kernel sources into a C++ class
+
+SET(CUDA_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src)
+SET(CUDA_SOURCE_CLASS CudaAmoebaKernelSources)
+SET(CUDA_KERNELS_CPP ${CMAKE_CURRENT_BINARY_DIR}/src/${CUDA_SOURCE_CLASS}.cpp)
+SET(CUDA_KERNELS_H ${CMAKE_CURRENT_BINARY_DIR}/src/${CUDA_SOURCE_CLASS}.h)
+SET(SOURCE_FILES ${SOURCE_FILES} ${CUDA_KERNELS_CPP} ${CUDA_KERNELS_H})
+INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}/src)
+
+# Create the library
+
+INCLUDE_DIRECTORIES(${CUDA_TOOLKIT_INCLUDE})
+
+FILE(GLOB CUDA_KERNELS ${CUDA_SOURCE_DIR}/kernels/*.cu)
+ADD_CUSTOM_COMMAND(OUTPUT ${CUDA_KERNELS_CPP} ${CUDA_KERNELS_H}
+    COMMAND ${CMAKE_COMMAND}
+    ARGS -D CUDA_SOURCE_DIR=${CUDA_SOURCE_DIR} -D CUDA_KERNELS_CPP=${CUDA_KERNELS_CPP} -D CUDA_KERNELS_H=${CUDA_KERNELS_H} -D CUDA_SOURCE_CLASS=${CUDA_SOURCE_CLASS} -P ${CMAKE_SOURCE_DIR}/platforms/cuda2/EncodeCUDAFiles.cmake
+    DEPENDS ${CUDA_KERNELS}
+)
+SET_SOURCE_FILES_PROPERTIES(${CUDA_KERNELS_CPP} ${CUDA_KERNELS_H} PROPERTIES GENERATED TRUE)
+ADD_LIBRARY(${SHARED_TARGET} SHARED ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES} ${API_ABS_INCLUDE_FILES})
+
+IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
+    SET(MAIN_OPENMM_LIB ${OPENMM_LIBRARY_NAME}_d)
+ELSE (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
+    SET(MAIN_OPENMM_LIB ${OPENMM_LIBRARY_NAME})
+ENDIF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
+TARGET_LINK_LIBRARIES(${SHARED_TARGET} ${MAIN_OPENMM_LIB}  ${CUDA_LIBRARIES} ${PTHREADS_LIB})
+TARGET_LINK_LIBRARIES(${SHARED_TARGET} debug ${OPENMM_LIBRARY_NAME}CUDA_d optimized ${OPENMM_LIBRARY_NAME}CUDA)
+TARGET_LINK_LIBRARIES(${SHARED_TARGET} debug ${SHARED_AMOEBA_TARGET} optimized ${SHARED_AMOEBA_TARGET})
+SET_TARGET_PROPERTIES(${SHARED_TARGET} PROPERTIES COMPILE_FLAGS "-DOPENMM_BUILDING_SHARED_LIBRARY")
+
+INSTALL(TARGETS ${SHARED_TARGET} DESTINATION ${CMAKE_INSTALL_PREFIX}/lib/plugins)
+# Ensure that links to the main CUDA library will be resolved.
+IF (APPLE)
+    IF (CMAKE_BUILD_TYPE MATCHES Debug)
+        SET(CUDA_LIBRARY libOpenMMCUDA_d.dylib)
+    ELSE (CMAKE_BUILD_TYPE MATCHES Debug)
+        SET(CUDA_LIBRARY libOpenMMCUDA.dylib)
+    ENDIF (CMAKE_BUILD_TYPE MATCHES Debug)
+    INSTALL(CODE "EXECUTE_PROCESS(COMMAND install_name_tool -change ${CUDA_LIBRARY} @loader_path/${CUDA_LIBRARY} ${CMAKE_INSTALL_PREFIX}/lib/plugins/lib${SHARED_TARGET}.dylib)")
+ENDIF (APPLE)
+SUBDIRS (tests)

plugins/amoeba/platforms/cuda2/include/AmoebaCudaKernelFactory.h

+#ifndef AMOEBA_OPENMM_CUDAKERNELFACTORY_H_
+#define AMOEBA_OPENMM_CUDAKERNELFACTORY_H_
+
+/* -------------------------------------------------------------------------- *
+ *                            OpenMMAmoeba                                    *
+ * -------------------------------------------------------------------------- *
+ * 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:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/KernelFactory.h"
+
+namespace OpenMM {
+
+/**
+ * This KernelFactory creates all kernels for AmoebaCudaPlatform.
+ */
+
+class AmoebaCudaKernelFactory : public KernelFactory {
+public:
+    KernelImpl* createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const;
+};
+
+} // namespace OpenMM
+
+#endif /*AMOEBA_OPENMM_CUDAKERNELFACTORY_H_*/

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

+/* -------------------------------------------------------------------------- *
+ *                              OpenMMAmoeba                                  *
+ * -------------------------------------------------------------------------- *
+ * 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-2012 Stanford University and the Authors.      *
+ * Authors: Mark Friedrichs, Peter Eastman                                    *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "AmoebaCudaKernelFactory.h"
+#include "AmoebaCudaKernels.h"
+#include "CudaPlatform.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/OpenMMException.h"
+#include "openmm/internal/windowsExport.h"
+
+using namespace OpenMM;
+
+extern "C" void registerPlatforms() {
+}
+
+extern "C" OPENMM_EXPORT void registerKernelFactories() {
+    try {
+        Platform& platform = Platform::getPlatformByName("CUDA");
+        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);
+        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);
+        platform.registerKernelFactory(CalcAmoebaVdwForceKernel::Name(), factory);
+        platform.registerKernelFactory(CalcAmoebaWcaDispersionForceKernel::Name(), factory);
+    }
+    catch (...) {
+        // Ignore.  The CUDA platform isn't available.
+    }
+}
+
+extern "C" OPENMM_EXPORT void registerAmoebaCudaKernelFactories() {
+    try {
+        Platform::getPlatformByName("CUDA");
+    }
+    catch (...) {
+        Platform::registerPlatform(new CudaPlatform());
+    }
+    registerKernelFactories();
+}
+
+KernelImpl* AmoebaCudaKernelFactory::createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const {
+    CudaPlatform::PlatformData& data = *static_cast<CudaPlatform::PlatformData*>(context.getPlatformData());
+    CudaContext& cu = *data.contexts[0];
+
+    if (name == CalcAmoebaHarmonicBondForceKernel::Name())
+        return new CudaCalcAmoebaHarmonicBondForceKernel(name, platform, cu, context.getSystem());
+
+    if (name == CalcAmoebaHarmonicAngleForceKernel::Name())
+        return new CudaCalcAmoebaHarmonicAngleForceKernel(name, platform, cu, context.getSystem());
+
+    if (name == CalcAmoebaHarmonicInPlaneAngleForceKernel::Name())
+        return new CudaCalcAmoebaHarmonicInPlaneAngleForceKernel(name, platform, cu, context.getSystem());
+
+//    if (name == CalcAmoebaTorsionForceKernel::Name())
+//        return new CudaCalcAmoebaTorsionForceKernel(name, platform, cu, context.getSystem());
+
+    if (name == CalcAmoebaPiTorsionForceKernel::Name())
+        return new CudaCalcAmoebaPiTorsionForceKernel(name, platform, cu, context.getSystem());
+
+//    if (name == CalcAmoebaStretchBendForceKernel::Name())
+//        return new CudaCalcAmoebaStretchBendForceKernel(name, platform, cu, context.getSystem());
+//
+//    if (name == CalcAmoebaOutOfPlaneBendForceKernel::Name())
+//        return new CudaCalcAmoebaOutOfPlaneBendForceKernel(name, platform, cu, context.getSystem());
+//
+//    if (name == CalcAmoebaTorsionTorsionForceKernel::Name())
+//        return new CudaCalcAmoebaTorsionTorsionForceKernel(name, platform, cu, context.getSystem());
+//
+//    if (name == CalcAmoebaMultipoleForceKernel::Name())
+//        return new CudaCalcAmoebaMultipoleForceKernel(name, platform, cu, context.getSystem());
+//
+//    if (name == CalcAmoebaGeneralizedKirkwoodForceKernel::Name())
+//        return new CudaCalcAmoebaGeneralizedKirkwoodForceKernel(name, platform, cu, context.getSystem());
+//
+//    if (name == CalcAmoebaVdwForceKernel::Name())
+//        return new CudaCalcAmoebaVdwForceKernel(name, platform, cu, context.getSystem());
+//
+//    if (name == CalcAmoebaWcaDispersionForceKernel::Name())
+//        return new CudaCalcAmoebaWcaDispersionForceKernel(name, platform, cu, context.getSystem());
+//
+//    if (name == CalcAmoebaUreyBradleyForceKernel::Name())
+//        return new CudaCalcAmoebaUreyBradleyForceKernel(name, platform, cu, context.getSystem());
+
+    throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '")+name+"'").c_str());
+}

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

+#ifndef AMOEBA_OPENMM_CUDAKERNELS_H_
+#define AMOEBA_OPENMM_CUDAKERNELS_H_
+
+/* -------------------------------------------------------------------------- *
+ *                              OpenMMAmoeba                                  *
+ * -------------------------------------------------------------------------- *
+ * 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-2012 Stanford University and the Authors.      *
+ * Authors: Mark Friedrichs, Peter Eastman                                    *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/amoebaKernels.h"
+#include "openmm/kernels.h"
+#include "openmm/System.h"
+#include "CudaContext.h"
+#include "CudaArray.h"
+
+namespace OpenMM {
+
+/**
+ * This kernel is invoked by AmoebaHarmonicBondForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CudaCalcAmoebaHarmonicBondForceKernel : public CalcAmoebaHarmonicBondForceKernel {
+public:
+    CudaCalcAmoebaHarmonicBondForceKernel(std::string name, 
+                                          const Platform& platform,
+                                          CudaContext& cu,
+                                          System& system);
+    ~CudaCalcAmoebaHarmonicBondForceKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param force      the AmoebaHarmonicBondForce this kernel will be used for
+     */
+    void initialize(const System& system, const AmoebaHarmonicBondForce& 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 numBonds;
+    CudaContext& cu;
+    System& system;
+    CudaArray* params;
+};
+
+/**
+ * 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,
+                                         CudaContext& cu,
+                                         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;
+    CudaContext& cu;
+    System& system;
+};
+
+/**
+ * This kernel is invoked by AmoebaHarmonicAngleForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CudaCalcAmoebaHarmonicAngleForceKernel : public CalcAmoebaHarmonicAngleForceKernel {
+public:
+    CudaCalcAmoebaHarmonicAngleForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system);
+    ~CudaCalcAmoebaHarmonicAngleForceKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param force      the AmoebaHarmonicAngleForce this kernel will be used for
+     */
+    void initialize(const System& system, const AmoebaHarmonicAngleForce& 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 numAngles;
+    CudaContext& cu;
+    System& system;
+    CudaArray* params;
+};
+
+/**
+ * This kernel is invoked by AmoebaHarmonicInPlaneAngleForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CudaCalcAmoebaHarmonicInPlaneAngleForceKernel : public CalcAmoebaHarmonicInPlaneAngleForceKernel {
+public:
+    CudaCalcAmoebaHarmonicInPlaneAngleForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system);
+    ~CudaCalcAmoebaHarmonicInPlaneAngleForceKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param force      the AmoebaHarmonicInPlaneAngleForce this kernel will be used for
+     */
+    void initialize(const System& system, const AmoebaHarmonicInPlaneAngleForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+private:
+    class ForceInfo;
+    int numAngles;
+    CudaContext& cu;
+    System& system;
+    CudaArray* params;
+};
+
+/**
+ * This kernel is invoked by AmoebaTorsionForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CudaCalcAmoebaTorsionForceKernel : public CalcAmoebaTorsionForceKernel {
+public:
+    CudaCalcAmoebaTorsionForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system);
+    ~CudaCalcAmoebaTorsionForceKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param force      the AmoebaTorsionForce this kernel will be used for
+     */
+    void initialize(const System& system, const AmoebaTorsionForce& 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 numTorsions;
+    CudaContext& cu;
+    System& system;
+};
+
+/**
+ * This kernel is invoked by AmoebaPiTorsionForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CudaCalcAmoebaPiTorsionForceKernel : public CalcAmoebaPiTorsionForceKernel {
+public:
+    CudaCalcAmoebaPiTorsionForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system);
+    ~CudaCalcAmoebaPiTorsionForceKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param force      the AmoebaPiTorsionForce this kernel will be used for
+     */
+    void initialize(const System& system, const AmoebaPiTorsionForce& 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 numPiTorsions;
+    CudaContext& cu;
+    System& system;
+    CudaArray* params;
+};
+
+/**
+ * This kernel is invoked by AmoebaStretchBendForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CudaCalcAmoebaStretchBendForceKernel : public CalcAmoebaStretchBendForceKernel {
+public:
+    CudaCalcAmoebaStretchBendForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system);
+    ~CudaCalcAmoebaStretchBendForceKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param force      the AmoebaStretchBendForce this kernel will be used for
+     */
+    void initialize(const System& system, const AmoebaStretchBendForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+private:
+    class ForceInfo;
+    int numStretchBends;
+    CudaContext& cu;
+    System& system;
+};
+
+/**
+ * This kernel is invoked by AmoebaOutOfPlaneBendForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CudaCalcAmoebaOutOfPlaneBendForceKernel : public CalcAmoebaOutOfPlaneBendForceKernel {
+public:
+    CudaCalcAmoebaOutOfPlaneBendForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system);
+    ~CudaCalcAmoebaOutOfPlaneBendForceKernel();
+    /**
+     * 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:
+    class ForceInfo;
+    int numOutOfPlaneBends;
+    CudaContext& cu;
+    System& system;
+};
+
+/**
+ * This kernel is invoked by AmoebaTorsionTorsionForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CudaCalcAmoebaTorsionTorsionForceKernel : public CalcAmoebaTorsionTorsionForceKernel {
+public:
+    CudaCalcAmoebaTorsionTorsionForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system);
+    ~CudaCalcAmoebaTorsionTorsionForceKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param force      the AmoebaTorsionTorsionForce this kernel will be used for
+     */
+    void initialize(const System& system, const AmoebaTorsionTorsionForce& 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 numTorsionTorsions;
+    int numTorsionTorsionGrids;
+    CudaContext& cu;
+    System& system;
+};
+
+/**
+ * This kernel is invoked by AmoebaMultipoleForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CudaCalcAmoebaMultipoleForceKernel : public CalcAmoebaMultipoleForceKernel {
+public:
+    CudaCalcAmoebaMultipoleForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system);
+    ~CudaCalcAmoebaMultipoleForceKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param force      the AmoebaMultipoleForce this kernel will be used for
+     */
+    void initialize(const System& system, const AmoebaMultipoleForce& 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);
+    /**
+     * Execute the kernel to calculate the electrostatic potential
+     *
+     * @param context        the context in which to execute this kernel
+     * @param inputGrid      input grid coordinates
+     * @param outputElectrostaticPotential output potential 
+     */
+    void getElectrostaticPotential(ContextImpl& context, const std::vector< Vec3 >& inputGrid,
+                                   std::vector< double >& outputElectrostaticPotential );
+
+   /** 
+     * Get the system multipole moments
+     *
+     * @param origin       origin
+     * @param context      context
+     * @param outputMultipoleMonents (charge,
+                                      dipole_x, dipole_y, dipole_z,
+                                      quadrupole_xx, quadrupole_xy, quadrupole_xz,
+                                      quadrupole_yx, quadrupole_yy, quadrupole_yz,
+                                      quadrupole_zx, quadrupole_zy, quadrupole_zz )
+     */
+
+    void getSystemMultipoleMoments( ContextImpl& context, const Vec3& origin, std::vector< double >& outputMultipoleMonents );
+
+
+private:
+    class ForceInfo;
+    int numMultipoles;
+    CudaContext& cu;
+    System& system;
+};
+
+/**
+ * This kernel is invoked by AmoebaMultipoleForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CudaCalcAmoebaGeneralizedKirkwoodForceKernel : public CalcAmoebaGeneralizedKirkwoodForceKernel {
+public:
+    CudaCalcAmoebaGeneralizedKirkwoodForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system);
+    ~CudaCalcAmoebaGeneralizedKirkwoodForceKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param force      the AmoebaMultipoleForce this kernel will be used for
+     */
+    void initialize(const System& system, const AmoebaGeneralizedKirkwoodForce& 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;
+    CudaContext& cu;
+    System& system;
+};
+
+/**
+ * This kernel is invoked to calculate the vdw forces acting on the system and the energy of the system.
+ */
+class CudaCalcAmoebaVdwForceKernel : public CalcAmoebaVdwForceKernel {
+public:
+    CudaCalcAmoebaVdwForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system);
+    ~CudaCalcAmoebaVdwForceKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param force      the AmoebaMultipoleForce this kernel will be used for
+     */
+    void initialize(const System& system, const AmoebaVdwForce& 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;
+    CudaContext& cu;
+    System& system;
+};
+
+/**
+ * This kernel is invoked to calculate the WCA dispersion forces acting on the system and the energy of the system.
+ */
+class CudaCalcAmoebaWcaDispersionForceKernel : public CalcAmoebaWcaDispersionForceKernel {
+public:
+    CudaCalcAmoebaWcaDispersionForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system);
+    ~CudaCalcAmoebaWcaDispersionForceKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param force      the AmoebaMultipoleForce this kernel will be used for
+     */
+    void initialize(const System& system, const AmoebaWcaDispersionForce& 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;
+    CudaContext& cu;
+    System& system;
+};
+
+} // namespace OpenMM
+
+#endif /*AMOEBA_OPENMM_CUDAKERNELS_H*/

plugins/amoeba/platforms/cuda2/src/CudaAmoebaKernelSources.cpp.in

+/* -------------------------------------------------------------------------- *
+ *                                   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) 2012 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "CudaAmoebaKernelSources.h"
+
+using namespace OpenMM;
+using namespace std;
+

plugins/amoeba/platforms/cuda2/src/CudaAmoebaKernelSources.h.in

+#ifndef OPENMM_CUDAAMOEBAKERNELSOURCES_H_
+#define OPENMM_CUDAAMOEBAKERNELSOURCES_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   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) 2010-2012 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/internal/windowsExport.h"
+#include <string>
+
+namespace OpenMM {
+
+/**
+ * This class is a central holding place for the source code of CUDA kernels.
+ * The CMake build script inserts declarations into it based on the .cu files in the
+ * kernels subfolder.
+ */
+
+class OPENMM_EXPORT CudaAmoebaKernelSources {
+public:
+@CUDA_FILE_DECLARATIONS@
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CUDAAMOEBAKERNELSOURCES_H_*/

plugins/amoeba/platforms/cuda2/src/kernels/amoebaAngleForce.cu

+float2 angleParams = PARAMS[index];
+real deltaIdeal = theta*RAD_TO_DEG-angleParams.x;
+real deltaIdeal2 = deltaIdeal*deltaIdeal;
+real deltaIdeal3 = deltaIdeal*deltaIdeal2;
+real deltaIdeal4 = deltaIdeal2*deltaIdeal2;
+energy += angleParams.y*deltaIdeal2*(1.0f + CUBIC_K*deltaIdeal + QUARTIC_K*deltaIdeal2 + PENTIC_K*deltaIdeal3 + SEXTIC_K*deltaIdeal4);
+real dEdAngle = angleParams.y*deltaIdeal*(2.0f + 3.0f*CUBIC_K*deltaIdeal + 4.0f*QUARTIC_K*deltaIdeal2 + 5.0f*PENTIC_K*deltaIdeal3 + 6.0f*SEXTIC_K*deltaIdeal4);
+dEdAngle *= RAD_TO_DEG;
\ No newline at end of file

plugins/amoeba/platforms/cuda2/src/kernels/amoebaBondForce.cu

+float2 bondParams = PARAMS[index];
+real deltaIdeal = r-bondParams.x;
+real deltaIdeal2 = deltaIdeal*deltaIdeal;
+energy += bondParams.y*deltaIdeal2*(1.0f + CUBIC_K*deltaIdeal + QUARTIC_K*deltaIdeal2);
+real dEdR = 2.0f*bondParams.y*deltaIdeal*(1.0f + 1.5f*CUBIC_K*deltaIdeal + 2.0f*QUARTIC_K*deltaIdeal2);

plugins/amoeba/platforms/cuda2/src/kernels/amoebaInPlaneForce.cu

+float2 angleParams = PARAMS[index];
+real xad = pos1.x - pos4.x;
+real yad = pos1.y - pos4.y;
+real zad = pos1.z - pos4.z;
+
+real xbd = pos2.x - pos4.x;
+real ybd = pos2.y - pos4.y;
+real zbd = pos2.z - pos4.z;
+
+real xcd = pos3.x - pos4.x;
+real ycd = pos3.y - pos4.y;
+real zcd = pos3.z - pos4.z;
+
+real xt = yad*zcd - zad*ycd;
+real yt = zad*xcd - xad*zcd;
+real zt = xad*ycd - yad*xcd;
+
+real rt2 = xt*xt + yt*yt + zt*zt;
+
+real delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
+
+real xip = pos2.x + xt*delta;
+real yip = pos2.y + yt*delta;
+real zip = pos2.z + zt*delta;
+
+real xap = pos1.x - xip;
+real yap = pos1.y - yip;
+real zap = pos1.z - zip;
+
+real xcp = pos3.x - xip;
+real ycp = pos3.y - yip;
+real zcp = pos3.z - zip;
+
+real rap2 = xap*xap + yap*yap + zap*zap;
+real rcp2 = xcp*xcp + ycp*ycp + zcp*zcp;
+
+real xm = ycp*zap - zcp*yap;
+real ym = zcp*xap - xcp*zap;
+real zm = xcp*yap - ycp*xap;
+
+real rm = max(SQRT(xm*xm + ym*ym + zm*zm), (real) 1e-6f);
+real dot = xap*xcp + yap*ycp + zap*zcp;
+real product = SQRT(rap2*rcp2);
+real cosine = (product > 0 ? (dot/product) : 0);
+cosine = max(min(cosine, (real) 1), (real) -1);
+real angle = ACOS(cosine);
+
+// if product == 0, set force/energy to 0
+
+real deltaIdeal = (product > 0 ? (angle*RAD_TO_DEG - angleParams.x) : 0);
+real deltaIdeal2 = deltaIdeal*deltaIdeal;
+real deltaIdeal3 = deltaIdeal*deltaIdeal2;
+real deltaIdeal4 = deltaIdeal2*deltaIdeal2;
+
+energy += angleParams.y*deltaIdeal2*(1.0f + CUBIC_K*deltaIdeal + QUARTIC_K*deltaIdeal2 + PENTIC_K*deltaIdeal3 + SEXTIC_K*deltaIdeal4);
+real dEdAngle = angleParams.y*deltaIdeal*(2.0f + 3.0f*CUBIC_K*deltaIdeal + 4.0f*QUARTIC_K*deltaIdeal2 + 5.0f*PENTIC_K*deltaIdeal3 + 6.0f*SEXTIC_K*deltaIdeal4);
+dEdAngle *= RAD_TO_DEG;
+
+real terma = -dEdAngle/(rap2*rm);
+real termc = dEdAngle/(rcp2*rm);
+
+real dedxia = terma * (yap*zm-zap*ym);
+real dedyia = terma * (zap*xm-xap*zm);
+real dedzia = terma * (xap*ym-yap*xm);
+
+real dedxic = termc * (ycp*zm-zcp*ym);
+real dedyic = termc * (zcp*xm-xcp*zm);
+real dedzic = termc * (xcp*ym-ycp*xm);
+
+real dedxip = -dedxia - dedxic;
+real dedyip = -dedyia - dedyic;
+real dedzip = -dedzia - dedzic;
+
+real delta2 = 2.0f*delta;
+real ptrt2 = (dedxip*xt + dedyip*yt + dedzip*zt) / rt2;
+
+real term = (zcd*ybd-ycd*zbd) + delta2*(yt*zcd-zt*ycd);
+real dpdxia = delta*(ycd*dedzip-zcd*dedyip) + term*ptrt2;
+
+term = (xcd*zbd-zcd*xbd) + delta2*(zt*xcd-xt*zcd);
+real dpdyia = delta*(zcd*dedxip-xcd*dedzip) + term*ptrt2;
+
+term = (ycd*xbd-xcd*ybd) + delta2*(xt*ycd-yt*xcd);
+real dpdzia = delta*(xcd*dedyip-ycd*dedxip) + term*ptrt2;
+
+term = (yad*zbd-zad*ybd) + delta2*(zt*yad-yt*zad);
+real dpdxic = delta*(zad*dedyip-yad*dedzip) + term*ptrt2;
+
+term = (zad*xbd-xad*zbd) + delta2*(xt*zad-zt*xad);
+real dpdyic = delta*(xad*dedzip-zad*dedxip) + term*ptrt2;
+
+term = (xad*ybd-yad*xbd) + delta2*(yt*xad-xt*yad);
+real dpdzic = delta*(yad*dedxip-xad*dedyip) + term*ptrt2;
+
+dedxia = dedxia + dpdxia;
+dedyia = dedyia + dpdyia;
+dedzia = dedzia + dpdzia;
+
+real dedxib = dedxip;
+real dedyib = dedyip;
+real dedzib = dedzip;
+
+dedxic = dedxic + dpdxic;
+dedyic = dedyic + dpdyic;
+dedzic = dedzic + dpdzic;
+
+real dedxid = -dedxia - dedxib - dedxic;
+real dedyid = -dedyia - dedyib - dedyic;
+real dedzid = -dedzia - dedzib - dedzic;
+
+real3 force1 = make_real3(-dedxia, -dedyia, -dedzia);
+real3 force2 = make_real3(-dedxib, -dedyib, -dedzib);
+real3 force3 = make_real3(-dedxic, -dedyic, -dedzic);
+real3 force4 = make_real3(-dedxid, -dedyid, -dedzid);
\ No newline at end of file

plugins/amoeba/platforms/cuda2/src/kernels/amoebaPiTorsionForce.cu

+// compute the value of the pi-orbital torsion angle
+
+real xad = pos1.x - pos4.x;
+real yad = pos1.y - pos4.y;
+real zad = pos1.z - pos4.z;
+
+real xbd = pos2.x - pos4.x;
+real ybd = pos2.y - pos4.y;
+real zbd = pos2.z - pos4.z;
+
+real xec = pos5.x - pos3.x;
+real yec = pos5.y - pos3.y;
+real zec = pos5.z - pos3.z;
+
+real xgc = pos6.x - pos3.x;
+real ygc = pos6.y - pos3.y;
+real zgc = pos6.z - pos3.z;
+
+real xip = yad*zbd - ybd*zad + pos3.x;
+real yip = zad*xbd - zbd*xad + pos3.y;
+real zip = xad*ybd - xbd*yad + pos3.z;
+
+real xiq = yec*zgc - ygc*zec + pos4.x;
+real yiq = zec*xgc - zgc*xec + pos4.y;
+real ziq = xec*ygc - xgc*yec + pos4.z;
+
+real xcp = pos3.x - xip;
+real ycp = pos3.y - yip;
+real zcp = pos3.z - zip;
+
+real xdc = pos4.x - pos3.x;
+real ydc = pos4.y - pos3.y;
+real zdc = pos4.z - pos3.z;
+
+real xqd = xiq - pos4.x;
+real yqd = yiq - pos4.y;
+real zqd = ziq - pos4.z;
+
+real xt = ycp*zdc - ydc*zcp;
+real yt = zcp*xdc - zdc*xcp;
+real zt = xcp*ydc - xdc*ycp;
+
+real xu = ydc*zqd - yqd*zdc;
+real yu = zdc*xqd - zqd*xdc;
+real zu = xdc*yqd - xqd*ydc;
+
+real xtu = yt*zu - yu*zt;
+real ytu = zt*xu - zu*xt;
+real ztu = xt*yu - xu*yt;
+
+real rt2 = xt*xt + yt*yt + zt*zt;
+real ru2 = xu*xu + yu*yu + zu*zu;
+
+real rtru = sqrtf(rt2 * ru2);
+real rdc = sqrtf(xdc*xdc + ydc*ydc + zdc*zdc);
+
+real cosine = rtru > 0.0f ? (xt*xu + yt*yu + zt*zu) / rtru : 0.0f;
+real sine = (rtru*rdc) > 0.0f ? (xdc*xtu + ydc*ytu + zdc*ztu) / (rdc*rtru) : 0.0f;
+
+// zero energy/force if rtru == 0
+
+float v2 = PARAMS[index];
+v2 = (rtru > 0 ? v2 : 0.0f);
+
+// compute the multiple angle trigonometry and the phase terms
+
+real cosine2 = cosine*cosine - sine*sine;
+real sine2 = 2.0f * cosine * sine;
+real phi2 = 1.0f - cosine2;
+real dphi2 = 2.0f * sine2;
+
+// calculate pi-orbital torsion energy and master chain rule term
+
+energy += v2 * phi2;
+real dedphi = v2 * dphi2;
+
+// chain rule terms for first derivative components
+
+real xdp = pos4.x - xip;
+real ydp = pos4.y - yip;
+real zdp = pos4.z - zip;
+
+real xqc = xiq - pos3.x;
+real yqc = yiq - pos3.y;
+real zqc = ziq - pos3.z;
+
+real dedxt = dedphi * (yt*zdc - ydc*zt) / (rt2*rdc);
+real dedyt = dedphi * (zt*xdc - zdc*xt) / (rt2*rdc);
+real dedzt = dedphi * (xt*ydc - xdc*yt) / (rt2*rdc);
+
+real dedxu = -dedphi * (yu*zdc - ydc*zu) / (ru2*rdc);
+real dedyu = -dedphi * (zu*xdc - zdc*xu) / (ru2*rdc);
+real dedzu = -dedphi * (xu*ydc - xdc*yu) / (ru2*rdc);
+
+// compute first derivative components for pi-orbital angle
+
+real dedxip = zdc*dedyt - ydc*dedzt;
+real dedyip = xdc*dedzt - zdc*dedxt;
+real dedzip = ydc*dedxt - xdc*dedyt;
+
+real dedxic = ydp*dedzt - zdp*dedyt + zqd*dedyu - yqd*dedzu;
+real dedyic = zdp*dedxt - xdp*dedzt + xqd*dedzu - zqd*dedxu;
+real dedzic = xdp*dedyt - ydp*dedxt + yqd*dedxu - xqd*dedyu;
+
+real dedxid = zcp*dedyt - ycp*dedzt + yqc*dedzu - zqc*dedyu;
+real dedyid = xcp*dedzt - zcp*dedxt + zqc*dedxu - xqc*dedzu;
+real dedzid = ycp*dedxt - xcp*dedyt + xqc*dedyu - yqc*dedxu;
+
+real dedxiq = zdc*dedyu - ydc*dedzu;
+real dedyiq = xdc*dedzu - zdc*dedxu;
+real dedziq = ydc*dedxu - xdc*dedyu;
+
+// compute first derivative components for individual atoms
+
+real dedxia = ybd*dedzip - zbd*dedyip;
+real dedyia = zbd*dedxip - xbd*dedzip;
+real dedzia = xbd*dedyip - ybd*dedxip;
+
+real dedxib = zad*dedyip - yad*dedzip;
+real dedyib = xad*dedzip - zad*dedxip;
+real dedzib = yad*dedxip - xad*dedyip;
+
+real dedxie = ygc*dedziq - zgc*dedyiq;
+real dedyie = zgc*dedxiq - xgc*dedziq;
+real dedzie = xgc*dedyiq - ygc*dedxiq;
+
+real dedxig = zec*dedyiq - yec*dedziq;
+real dedyig = xec*dedziq - zec*dedxiq;
+real dedzig = yec*dedxiq - xec*dedyiq;
+
+dedxic = dedxic + dedxip - dedxie - dedxig;
+dedyic = dedyic + dedyip - dedyie - dedyig;
+dedzic = dedzic + dedzip - dedzie - dedzig;
+dedxid = dedxid + dedxiq - dedxia - dedxib;
+dedyid = dedyid + dedyiq - dedyia - dedyib;
+dedzid = dedzid + dedziq - dedzia - dedzib;
+
+real3 force1 = make_real3(-dedxia, -dedyia, -dedzia);
+real3 force2 = make_real3(-dedxib, -dedyib, -dedzib);
+real3 force3 = make_real3(-dedxic, -dedyic, -dedzic);
+real3 force4 = make_real3(-dedxid, -dedyid, -dedzid);
+real3 force5 = make_real3(-dedxie, -dedyie, -dedzie);
+real3 force6 = make_real3(-dedxig, -dedyig, -dedzig);
\ No newline at end of file

plugins/amoeba/platforms/cuda2/tests/AmoebaTinkerParameterFile.h

+/* -------------------------------------------------------------------------- *
+ *                                   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) 2009 Stanford University and the Authors.           *
+ * Authors: Peter Eastman, Mark Friedrichs                                    *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/internal/AssertionUtilities.h"
+
+#include "openmm/CMMotionRemover.h"
+#include "openmm/System.h"
+#include "openmm/Context.h"
+#include "openmm/LangevinIntegrator.h"
+#include "openmm/VariableLangevinIntegrator.h"
+#include "openmm/VerletIntegrator.h"
+#include "openmm/VariableVerletIntegrator.h"
+#include "openmm/BrownianIntegrator.h"
+
+#include "openmm/AmoebaHarmonicBondForce.h"
+#include "openmm/AmoebaHarmonicAngleForce.h"
+#include "openmm/AmoebaHarmonicInPlaneAngleForce.h"
+#include "openmm/AmoebaTorsionForce.h"
+#include "openmm/AmoebaPiTorsionForce.h"
+#include "openmm/AmoebaStretchBendForce.h"
+#include "openmm/AmoebaOutOfPlaneBendForce.h"
+#include "openmm/AmoebaTorsionTorsionForce.h"
+#include "openmm/AmoebaMultipoleForce.h"
+#include "openmm/AmoebaGeneralizedKirkwoodForce.h"
+#include "openmm/AmoebaVdwForce.h"
+#include "openmm/AmoebaWcaDispersionForce.h"
+#include "openmm/AmoebaUreyBradleyForce.h"
+#include "openmm/internal/windowsExport.h"
+#include "openmm/internal/AmoebaWcaDispersionForceImpl.h"
+
+#include <ctime>
+#include <vector>
+#include <algorithm>
+#include <map>
+#include <cfloat>
+#include <cstring>
+#include <cstdlib>
+#include <typeinfo>
+#include <time.h>
+
+extern "C" void registerAmoebaCudaKernelFactories();
+
+// force enums
+
+#define MAX_PRINT 5
+
+static std::string AMOEBA_HARMONIC_BOND_FORCE                         = "AmoebaHarmonicBond";
+static std::string AMOEBA_HARMONIC_ANGLE_FORCE                        = "AmoebaHarmonicAngle"; 
+static std::string AMOEBA_HARMONIC_IN_PLANE_ANGLE_FORCE               = "AmoebaHarmonicInPlaneAngle"; 
+static std::string AMOEBA_TORSION_FORCE                               = "AmoebaTorsion"; 
+static std::string AMOEBA_PI_TORSION_FORCE                            = "AmoebaPiTorsion"; 
+static std::string AMOEBA_STRETCH_BEND_FORCE                          = "AmoebaStretchBend";
+static std::string AMOEBA_OUT_OF_PLANE_BEND_FORCE                     = "AmoebaOutOfPlaneBend";
+static std::string AMOEBA_TORSION_TORSION_FORCE                       = "AmoebaTorsionTorsion";
+static std::string AMOEBA_MULTIPOLE_FORCE                             = "AmoebaMultipole";
+static std::string AMOEBA_GK_FORCE                                    = "AmoebaGk";
+static std::string AMOEBA_GK_CAVITY_FORCE                             = "AmoebaGkAndCavity";
+static std::string AMOEBA_VDW_FORCE                                   = "AmoebaVdw";
+static std::string AMOEBA_WCA_DISPERSION_FORCE                        = "AmoebaWcaDispersion";
+static std::string AMOEBA_UREY_BRADLEY_FORCE                          = "AmoebaUreyBradley";
+static std::string ALL_FORCES                                         = "AllForces";
+
+static std::string AMOEBA_MULTIPOLE_ROTATION_MATRICES                 = "AmoebaMultipoleRotationMatrices";
+static std::string AMOEBA_MULTIPOLE_ROTATED                           = "AmoebaMultipolesRotated";
+static std::string AMOEBA_FIXED_E                                     = "AmoebaFixedE";
+static std::string AMOEBA_FIXED_E_GK                                  = "AmoebaFixedE_GK";
+static std::string AMOEBA_INDUCDED_DIPOLES                            = "AmoebaInducedDipoles";
+static std::string AMOEBA_INDUCDED_DIPOLES_GK                         = "AmoebaInducedDipoles_GK";
+
+static std::string INCLUDE_OBC_CAVITY_TERM                            = "includeObcCavityTerm";
+static std::string MUTUAL_INDUCED_MAX_ITERATIONS                      = "mutualInducedMaxIterations";
+static std::string MUTUAL_INDUCED_TARGET_EPSILON                      = "mutualInducedTargetEpsilon";
+static std::string APPLY_N2                                           = "applyN2";
+static std::string ZERO_HARMONIC_BOND_IXN                             = "zeroHarmonicBondIxn";
+
+#define AmoebaHarmonicBondIndex                            0
+#define AmoebaHarmonicAngleIndex                           1
+#define AmoebaHarmonicInPlaneAngleIndex                    2
+#define AmoebaTorsionIndex                                 3
+#define AmoebaPiTorsionIndex                               4
+#define AmoebaStretchBendIndex                             5
+#define AmoebaOutOfPlaneBendIndex                          6
+#define AmoebaTorsionTorsionIndex                          7
+#define AmoebaMultipoleIndex                               8
+#define AmoebaVdwIndex                                     9
+#define AmoebaWcaDispersionIndex                          10
+#define AmoebaObcIndex                                    11
+#define SumIndex                                          12
+#define UreyBradleyIndex                                  13
+#define AmoebaLastIndex                                   14
+
+#define AngstromToNm 0.1
+#define CalToJoule   4.184
+
+const double DegreesToRadians = 3.14159265/180.0;
+const double RadiansToDegrees = 180/3.14159265;
+
+using namespace OpenMM;
+using namespace std;
+
+// the following are used in parsing parameter file
+
+typedef std::vector<std::string> StringVector;
+typedef StringVector::iterator StringVectorI;
+typedef StringVector::const_iterator StringVectorCI;
+
+typedef std::vector<StringVector> StringVectorVector;
+
+typedef std::vector<std::vector<double> > VectorOfVectors;
+typedef VectorOfVectors::iterator VectorOfVectorsI;
+typedef VectorOfVectors::const_iterator VectorOfVectorsCI;
+
+typedef std::map< std::string, VectorOfVectors > MapStringVectorOfVectors;
+typedef MapStringVectorOfVectors::iterator MapStringVectorOfVectorsI;
+typedef MapStringVectorOfVectors::const_iterator MapStringVectorOfVectorsCI;
+
+typedef std::map< std::string, std::string > MapStringString;
+typedef MapStringString::iterator MapStringStringI;
+typedef MapStringString::const_iterator MapStringStringCI;
+
+typedef std::map< std::string, int > MapStringInt;
+typedef MapStringInt::iterator MapStringIntI;
+typedef MapStringInt::const_iterator MapStringIntCI;
+
+typedef std::map< std::string,  std::vector<Vec3> > MapStringVec3;
+typedef MapStringVec3::iterator MapStringVec3I;
+typedef MapStringVec3::const_iterator MapStringVec3CI;
+
+typedef std::map< std::string, double > MapStringDouble;
+typedef MapStringDouble::iterator MapStringDoubleI;
+typedef MapStringDouble::const_iterator MapStringDoubleCI;
+
+typedef std::map< std::string, Force*> MapStringForce;
+typedef MapStringForce::iterator MapStringForceI;
+typedef MapStringForce::const_iterator MapStringForceCI;
+
+// default return value from methods
+
+static const int DefaultReturnValue               = 0;
+
+
+static const int LengthUnit                       = 0;
+static const int EnergyUnit                       = 1;
+static const int ForceUnit                        = 2;
+static const int LastUnits                        = ForceUnit + 1;
+
+static const int NoUnitsConversion                = 0;
+static const int KcalA_To_kJNm                    = 1;
+
+/**---------------------------------------------------------------------------------------
+ * Initialize units
+ *
+ * @param unitType        has w/ force name as key and int as value
+ * @param units array
+ *
+ *
+    --------------------------------------------------------------------------------------- */
+
+void setUnits( int unitType, double* units );
+
+/**---------------------------------------------------------------------------------------
+
+    Read parameter file
+
+    @param inputParameterFile   input parameter file name
+    @param system               system to which forces based on parameters are to be added
+    @param coordinates          Vec3 array containing coordinates on output
+    @param velocities           Vec3 array containing velocities on output
+    @param inputLog             log file pointer -- may be NULL
+
+    @return number of lines read
+
+    --------------------------------------------------------------------------------------- */
+
+Integrator* readAmoebaParameterFile( const std::string& inputParameterFile, MapStringInt& forceMap, System& system,
+                                     std::vector<Vec3>& coordinates, 
+                                     std::vector<Vec3>& velocities,
+                                     MapStringVec3& forces, MapStringDouble& potentialEnergy,
+                                     MapStringVectorOfVectors& supplementary, FILE* inputLog );
+
+
+/**---------------------------------------------------------------------------------------
+ * Get integrator
+ * 
+ * @param  integratorName       integratorName (VerletIntegrator, BrownianIntegrator, LangevinIntegrator, ...)
+ * @param  timeStep             time step
+ * @param  friction (ps)        friction
+ * @param  temperature          temperature
+ * @param  shakeTolerance       Shake tolerance
+ * @param  errorTolerance       Error tolerance
+ * @param  randomNumberSeed     seed
+ *
+ * @return DefaultReturnValue or ErrorReturnValue
+ *
+    --------------------------------------------------------------------------------------- */
+
+Integrator* getIntegrator( std::string& integratorName, double timeStep,
+                             double friction, double temperature,
+                             double shakeTolerance, double errorTolerance,
+                             int randomNumberSeed, FILE* log );
+
+/**---------------------------------------------------------------------------------------
+ * Initialize forceMap
+ *
+ * @param forceMap        has w/ force name as key and int as value
+ * @param initialValue    initial value
+ *
+ *
+    --------------------------------------------------------------------------------------- */
+
+void initializeForceMap( MapStringInt& forceMap, int initialValue );
+
+void testUsingAmoebaTinkerParameterFile( const std::string& amoebaTinkerParameterFileName, MapStringInt& forceMap,
+                                         double tolerance, FILE* summaryFile, FILE* log );
+
+int OPENMM_EXPORT runTestsUsingAmoebaTinkerParameterFile( MapStringString& argumentMap );
+void OPENMM_EXPORT appendInputArgumentsToArgumentMap( int numberOfArguments, char* argv[], MapStringString& argumentMap );

plugins/amoeba/platforms/cuda2/tests/CMakeLists.txt

+#
+# Testing
+#
+
+ENABLE_TESTING()
+
+INCLUDE_DIRECTORIES(${CUDA_INCLUDE_DIR})
+
+# Automatically create tests using files named "Test*.cpp"
+FILE(GLOB TEST_PROGS "*Test*.cpp")
+FOREACH(TEST_PROG ${TEST_PROGS})
+    GET_FILENAME_COMPONENT(TEST_ROOT ${TEST_PROG} NAME_WE)
+
+    # Link with shared library
+    ADD_EXECUTABLE(${TEST_ROOT} ${TEST_PROG})
+    TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_AMOEBA_TARGET} ${SHARED_TARGET})
+    ADD_TEST(${TEST_ROOT} ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT})
+
+ENDFOREACH(TEST_PROG ${TEST_PROGS})

plugins/amoeba/platforms/cuda2/tests/TestCudaAmoebaHarmonicAngleForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMMAmoeba                             *
+ * -------------------------------------------------------------------------- *
+ * 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: Mark Friedrichs                                                   *
+ * 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 CudaAmoebaHarmonicAngleForce.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "AmoebaTinkerParameterFile.h"
+#include "openmm/Context.h"
+#include "OpenMMAmoeba.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+
+const double TOL = 1e-5;
+#define PI_M               3.141592653589
+#define RADIAN            57.29577951308
+#define RADIAN_TO_DEGREE  57.29577951308
+#define DEGREE_TO_RADIAN   0.01745329252
+#define RADIAN_INVERSE     0.01745329252
+
+/* ---------------------------------------------------------------------------------------
+
+   Compute cross product of two 3-vectors and place in 3rd vector
+
+   vectorZ = vectorX x vectorY
+
+   @param vectorX             x-vector
+   @param vectorY             y-vector
+   @param vectorZ             z-vector
+
+   @return vector is vectorZ
+
+   --------------------------------------------------------------------------------------- */
+     
+static void crossProductVector3( double* vectorX, double* vectorY, double* vectorZ ){
+
+    vectorZ[0]  = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1];
+    vectorZ[1]  = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2];
+    vectorZ[2]  = vectorX[0]*vectorY[1] - vectorX[1]*vectorY[0];
+
+    return;
+}
+
+static void getPrefactorsGivenAngleCosine( double cosine, double idealAngle, double quadraticK, double cubicK,
+                                           double quarticK, double penticK, double sexticK,
+                                           double* dEdR, double* energyTerm, FILE* log ) {
+
+    double angle;
+    if( cosine >= 1.0 ){
+        angle = 0.0f;
+    } else if( cosine <= -1.0 ){
+        angle = RADIAN*PI_M;
+    } else {
+        angle = RADIAN*acos(cosine);
+    }
+
+#ifdef AMOEBA_DEBUG
+    if( log ){
+        (void) fprintf( log, "getPrefactorsGivenAngleCosine: cosine=%10.3e angle=%10.3e ideal=%10.3e\n", cosine, angle, idealAngle ); 
+        (void) fflush( log );
+    }
+#endif
+
+    double deltaIdeal         = angle - idealAngle;
+    double deltaIdeal2        = deltaIdeal*deltaIdeal;
+    double deltaIdeal3        = deltaIdeal*deltaIdeal2;
+    double deltaIdeal4        = deltaIdeal2*deltaIdeal2;
+ 
+    // deltaIdeal = r - r_0
+ 
+    *dEdR        = ( 2.0                        +
+                     3.0*cubicK*  deltaIdeal    +
+                     4.0*quarticK*deltaIdeal2   +
+                     5.0*penticK* deltaIdeal3   +
+                     6.0*sexticK* deltaIdeal4     );
+ 
+    *dEdR       *= RADIAN*quadraticK*deltaIdeal;
+ 
+
+    *energyTerm  = 1.0f + cubicK* deltaIdeal    +
+                          quarticK*deltaIdeal2   +
+                          penticK* deltaIdeal3   +
+                          sexticK* deltaIdeal4;
+    *energyTerm *= quadraticK*deltaIdeal2;
+
+    return;
+}
+
+static void computeAmoebaHarmonicAngleForce(int bondIndex,  std::vector<Vec3>& positions, AmoebaHarmonicAngleForce& amoebaHarmonicAngleForce,
+                                             std::vector<Vec3>& forces, double* energy, FILE* log ) {
+
+    int particle1, particle2, particle3;
+    double idealAngle;
+    double quadraticK;
+    amoebaHarmonicAngleForce.getAngleParameters(bondIndex, particle1, particle2, particle3, idealAngle, quadraticK );
+
+    double cubicK         = amoebaHarmonicAngleForce.getAmoebaGlobalHarmonicAngleCubic();
+    double quarticK       = amoebaHarmonicAngleForce.getAmoebaGlobalHarmonicAngleQuartic();
+    double penticK        = amoebaHarmonicAngleForce.getAmoebaGlobalHarmonicAnglePentic();
+    double sexticK        = amoebaHarmonicAngleForce.getAmoebaGlobalHarmonicAngleSextic();
+
+#ifdef AMOEBA_DEBUG
+    if( log ){
+        (void) fprintf( log, "computeAmoebaHarmonicAngleForce: bond %d [%d %d %d] ang=%10.3f k=%10.3f [%10.3e %10.3e %10.3e %10.3e]\n", 
+                             bondIndex, particle1, particle2, particle3, idealAngle, quadraticK, cubicK, quarticK, penticK, sexticK );
+        (void) fflush( log );
+    }
+#endif
+
+    double deltaR[2][3];
+    double r2_0 = 0.0;
+    double r2_1 = 0.0;
+    for( int ii = 0; ii < 3; ii++ ){
+
+           deltaR[0][ii]    = positions[particle1][ii] - positions[particle2][ii];
+           r2_0            += deltaR[0][ii]*deltaR[0][ii];
+
+           deltaR[1][ii]    = positions[particle3][ii] - positions[particle2][ii];
+           r2_1            += deltaR[1][ii]*deltaR[1][ii];
+
+    }
+
+    double pVector[3];
+    crossProductVector3( deltaR[0], deltaR[1], pVector );
+    double rp      = sqrt( pVector[0]*pVector[0] + pVector[1]*pVector[1] + pVector[2]*pVector[2] );
+    if( rp < 1.0e-06 ){
+       rp = 1.0e-06;
+    }   
+    double dot    = deltaR[0][0]*deltaR[1][0] + deltaR[0][1]*deltaR[1][1] + deltaR[0][2]*deltaR[1][2];
+    double cosine = dot/sqrt(r2_0*r2_1);
+
+#ifdef AMOEBA_DEBUG
+    if( log ){
+        (void) fprintf( log, "dot=%10.3e r2_0=%10.3e r2_1=%10.3e\n", dot, r2_0, r2_1 ); 
+        (void) fflush( log );
+    }
+#endif
+
+    double dEdR;
+    double energyTerm;
+    getPrefactorsGivenAngleCosine( cosine, idealAngle, quadraticK, cubicK,
+                                   quarticK, penticK, sexticK, &dEdR, &energyTerm, log );
+
+    double termA  = -dEdR/(r2_0*rp);
+    double termC  =  dEdR/(r2_1*rp);
+
+    double deltaCrossP[3][3];
+    crossProductVector3( deltaR[0], pVector, deltaCrossP[0] );
+    crossProductVector3( deltaR[1], pVector, deltaCrossP[2] );
+    for( int ii = 0; ii < 3; ii++ ){
+        deltaCrossP[0][ii] *= termA;
+        deltaCrossP[2][ii] *= termC;
+        deltaCrossP[1][ii]  = -1.0*(deltaCrossP[0][ii] + deltaCrossP[2][ii]);
+    }
+
+    forces[particle1][0]       += deltaCrossP[0][0];
+    forces[particle1][1]       += deltaCrossP[0][1];
+    forces[particle1][2]       += deltaCrossP[0][2];
+
+    forces[particle2][0]       += deltaCrossP[1][0];
+    forces[particle2][1]       += deltaCrossP[1][1];
+    forces[particle2][2]       += deltaCrossP[1][2];
+
+    forces[particle3][0]       += deltaCrossP[2][0];
+    forces[particle3][1]       += deltaCrossP[2][1];
+    forces[particle3][2]       += deltaCrossP[2][2];
+
+    *energy                    += energyTerm;
+}
+
+static void computeAmoebaHarmonicAngleForces( Context& context, AmoebaHarmonicAngleForce& amoebaHarmonicAngleForce,
+                                             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 < amoebaHarmonicAngleForce.getNumAngles(); ii++ ){
+        computeAmoebaHarmonicAngleForce(ii, positions, amoebaHarmonicAngleForce, expectedForces, expectedEnergy, log );
+    }
+
+#ifdef AMOEBA_DEBUG
+    if( log ){
+        (void) fprintf( log, "computeAmoebaHarmonicAngleForces: expected energy=%14.7e\n", *expectedEnergy );
+        for( unsigned int ii = 0; ii < positions.size(); ii++ ){
+            (void) fprintf( log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2] );
+        }
+        (void) fflush( log );
+    }
+#endif
+
+    return;
+
+}
+
+void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicAngleForce& amoebaHarmonicAngleForce,
+                                        double tolerance, const std::string& idString, FILE* log) {
+
+    std::vector<Vec3> expectedForces;
+    double expectedEnergy;
+    computeAmoebaHarmonicAngleForces( context, amoebaHarmonicAngleForce, expectedForces, &expectedEnergy, log );
+   
+    State state                      = context.getState(State::Forces | State::Energy);
+    const std::vector<Vec3> forces   = state.getForces();
+
+#ifdef AMOEBA_DEBUG
+    if( log ){
+        (void) fprintf( log, "computeAmoebaHarmonicAngleForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
+        for( unsigned int ii = 0; ii < forces.size(); ii++ ){
+            (void) fprintf( log, "%6u [%14.7e %14.7e %14.7e]   [%14.7e %14.7e %14.7e]\n", ii,
+                            expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
+        }
+        (void) fflush( log );
+    }
+#endif
+
+    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 testOneAngle( FILE* log ) {
+
+    System system;
+    int numberOfParticles = 3;
+    for( int ii = 0; ii < numberOfParticles; ii++ ){
+        system.addParticle(1.0);
+    }
+
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+
+    AmoebaHarmonicAngleForce* amoebaHarmonicAngleForce = new AmoebaHarmonicAngleForce();
+
+    double angle      = 100.0;
+    double quadraticK = 1.0;
+    double cubicK     = 1.0e-01;
+    double quarticK   = 1.0e-02;
+    double penticK    = 1.0e-03;
+    double sexticK    = 1.0e-04;
+    amoebaHarmonicAngleForce->addAngle(0, 1, 2, angle, quadraticK);
+
+    amoebaHarmonicAngleForce->setAmoebaGlobalHarmonicAngleCubic(cubicK);
+    amoebaHarmonicAngleForce->setAmoebaGlobalHarmonicAngleQuartic(quarticK);
+    amoebaHarmonicAngleForce->setAmoebaGlobalHarmonicAnglePentic(penticK);
+    amoebaHarmonicAngleForce->setAmoebaGlobalHarmonicAngleSextic(sexticK);
+
+    system.addForce(amoebaHarmonicAngleForce);
+    Context context(system, integrator, Platform::getPlatformByName( "CUDA"));
+
+    std::vector<Vec3> positions(numberOfParticles);
+
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(0, 0, 1);
+
+    context.setPositions(positions);
+    compareWithExpectedForceAndEnergy( context, *amoebaHarmonicAngleForce, TOL, "testOneAngle", log );
+
+}
+
+int main( int numberOfArguments, char* argv[] ) {
+
+    try {
+        std::cout << "TestCudaAmoebaHarmonicAngleForce running test..." << std::endl;
+        registerAmoebaCudaKernelFactories();
+        FILE* log = NULL;
+        testOneAngle( 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 << "Done" << std::endl;
+    return 0;
+}

plugins/amoeba/platforms/cuda2/tests/TestCudaAmoebaHarmonicBondForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMMAmoeba                             *
+ * -------------------------------------------------------------------------- *
+ * 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: Mark Friedrichs                                                   *
+ * 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 HarmonicBondForce.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "CudaPlatform.h"
+#include "AmoebaTinkerParameterFile.h"
+#include "openmm/Context.h"
+#include "OpenMMAmoeba.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+
+const double TOL = 1e-5;
+static void computeAmoebaHarmonicBondForce(int bondIndex,  std::vector<Vec3>& positions, AmoebaHarmonicBondForce& amoebaHarmonicBondForce,
+                                           std::vector<Vec3>& forces, double* energy ) {
+
+    int particle1, particle2;
+    double bondLength;
+    double quadraticK;
+    double cubicK    = amoebaHarmonicBondForce.getAmoebaGlobalHarmonicBondCubic();
+    double quarticK  = amoebaHarmonicBondForce.getAmoebaGlobalHarmonicBondQuartic();
+    amoebaHarmonicBondForce.getBondParameters(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 computeAmoebaHarmonicBondForces( Context& context, AmoebaHarmonicBondForce& amoebaHarmonicBondForce,
+                                             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 < amoebaHarmonicBondForce.getNumBonds(); ii++ ){
+        computeAmoebaHarmonicBondForce(ii, positions, amoebaHarmonicBondForce, expectedForces, expectedEnergy );
+    }
+
+#ifdef AMOEBA_DEBUG
+    if( log ){
+        (void) fprintf( log, "computeAmoebaHarmonicBondForces: 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 );
+    }
+#endif
+    return;
+
+}
+
+void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicBondForce& amoebaHarmonicBondForce, double tolerance, const std::string& idString, FILE* log) {
+
+    std::vector<Vec3> expectedForces;
+    double expectedEnergy;
+    computeAmoebaHarmonicBondForces( context, amoebaHarmonicBondForce, expectedForces, &expectedEnergy, NULL );
+   
+    State state                      = context.getState(State::Forces | State::Energy);
+    const std::vector<Vec3> forces   = state.getForces();
+
+#ifdef AMOEBA_DEBUG
+    if( log ){
+        (void) fprintf( log, "computeAmoebaHarmonicBondForces: 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 );
+    }
+#endif
+
+    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);
+
+    AmoebaHarmonicBondForce* amoebaHarmonicBondForce = new AmoebaHarmonicBondForce();
+
+    double bondLength = 1.5;
+    double quadraticK = 1.0;
+    double cubicK     = 2.0;
+    double quarticicK = 3.0;
+    amoebaHarmonicBondForce->setAmoebaGlobalHarmonicBondCubic( cubicK );
+    amoebaHarmonicBondForce->setAmoebaGlobalHarmonicBondQuartic( quarticicK );
+    amoebaHarmonicBondForce->addBond(0, 1, bondLength, quadraticK);
+
+    system.addForce(amoebaHarmonicBondForce);
+    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, *amoebaHarmonicBondForce, 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);
+
+    AmoebaHarmonicBondForce* amoebaHarmonicBondForce = new AmoebaHarmonicBondForce();
+
+    double bondLength = 1.5;
+    double quadraticK = 1.0;
+    double cubicK     = 2.0;
+    double quarticicK = 3.0;
+    amoebaHarmonicBondForce->setAmoebaGlobalHarmonicBondCubic( cubicK );
+    amoebaHarmonicBondForce->setAmoebaGlobalHarmonicBondQuartic( quarticicK );
+    amoebaHarmonicBondForce->addBond(0, 1, bondLength, quadraticK);
+    amoebaHarmonicBondForce->addBond(1, 2, bondLength, quadraticK);
+
+    system.addForce(amoebaHarmonicBondForce);
+    Context context(system, integrator, Platform::getPlatformByName( "CUDA"));
+    //Context context(system, integrator, platform );
+    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, *amoebaHarmonicBondForce, TOL, "testTwoBond", log );
+}
+
+int main( int numberOfArguments, char* argv[] ) {
+
+    try {
+        std::cout << "TestCudaAmoebaHarmonicBondForce running test..." << std::endl;
+        registerAmoebaCudaKernelFactories();
+        FILE* log = NULL;
+        testTwoBond( 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 << "Done" << std::endl;
+    return 0;
+}

plugins/amoeba/platforms/cuda2/tests/TestCudaAmoebaHarmonicInPlaneAngleForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMMAmoeba                             *
+ * -------------------------------------------------------------------------- *
+ * 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: Mark Friedrichs                                                   *
+ * 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 CudaAmoebaHarmonicInPlaneAngleForce.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "AmoebaTinkerParameterFile.h"
+#include "openmm/Context.h"
+#include "OpenMMAmoeba.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+
+const double TOL = 1e-5;
+#define PI_M               3.141592653589
+#define RADIAN            57.29577951308
+
+/* ---------------------------------------------------------------------------------------
+
+   Compute cross product of two 3-vectors and place in 3rd vector
+
+   vectorZ = vectorX x vectorY
+
+   @param vectorX             x-vector
+   @param vectorY             y-vector
+   @param vectorZ             z-vector
+
+   @return vector is vectorZ
+
+   --------------------------------------------------------------------------------------- */
+     
+static void crossProductVector3( double* vectorX, double* vectorY, double* vectorZ ){
+
+    vectorZ[0]  = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1];
+    vectorZ[1]  = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2];
+    vectorZ[2]  = vectorX[0]*vectorY[1] - vectorX[1]*vectorY[0];
+
+    return;
+}
+
+static double dotVector3( double* vectorX, double* vectorY ){
+    return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2];
+}
+
+static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInPlaneAngle, double quadraticK, double cubicK,
+                                                  double quarticK, double penticK, double sexticK,
+                                                  double* dEdR, double* energyTerm, FILE* log ) {
+
+    double angle;
+    if( cosine >= 1.0 ){
+        angle = 0.0f;
+    } else if( cosine <= -1.0 ){
+        angle = RADIAN*PI_M;
+    } else {
+        angle = RADIAN*acos(cosine);
+    }
+
+#ifdef AMOEBA_DEBUG
+    if( log ){
+        (void) fprintf( log, "getPrefactorsGivenInPlaneAngleCosine: cosine=%10.3e angle=%10.3e ideal=%10.3e\n", cosine, angle, idealInPlaneAngle ); 
+        (void) fflush( log );
+    }
+#endif
+
+    double deltaIdeal         = angle - idealInPlaneAngle;
+    double deltaIdeal2        = deltaIdeal*deltaIdeal;
+    double deltaIdeal3        = deltaIdeal*deltaIdeal2;
+    double deltaIdeal4        = deltaIdeal2*deltaIdeal2;
+ 
+    // deltaIdeal = r - r_0
+ 
+    *dEdR        = ( 2.0                        +
+                     3.0*cubicK*  deltaIdeal    +
+                     4.0*quarticK*deltaIdeal2   +
+                     5.0*penticK* deltaIdeal3   +
+                     6.0*sexticK* deltaIdeal4     );
+ 
+    *dEdR       *= RADIAN*quadraticK*deltaIdeal;
+ 
+
+    *energyTerm  = 1.0f + cubicK* deltaIdeal    +
+                          quarticK*deltaIdeal2   +
+                          penticK* deltaIdeal3   +
+                          sexticK* deltaIdeal4;
+    *energyTerm *= quadraticK*deltaIdeal2;
+
+    return;
+}
+
+static void computeAmoebaHarmonicInPlaneAngleForce(int bondIndex,  std::vector<Vec3>& positions, AmoebaHarmonicInPlaneAngleForce& amoebaHarmonicInPlaneAngleForce,
+                                                   std::vector<Vec3>& forces, double* energy, FILE* log ) {
+
+    int particle1, particle2, particle3, particle4;
+    double idealInPlaneAngle;
+    double quadraticK;
+    amoebaHarmonicInPlaneAngleForce.getAngleParameters(bondIndex, particle1, particle2, particle3, particle4, idealInPlaneAngle, quadraticK );
+
+    double cubicK         = amoebaHarmonicInPlaneAngleForce.getAmoebaGlobalHarmonicInPlaneAngleCubic();
+    double quarticK       = amoebaHarmonicInPlaneAngleForce.getAmoebaGlobalHarmonicInPlaneAngleQuartic();
+    double penticK        = amoebaHarmonicInPlaneAngleForce.getAmoebaGlobalHarmonicInPlaneAnglePentic();
+    double sexticK        = amoebaHarmonicInPlaneAngleForce.getAmoebaGlobalHarmonicInPlaneAngleSextic();
+
+#ifdef AMOEBA_DEBUG
+    if( log ){
+        (void) fprintf( log, "computeAmoebaHarmonicInPlaneAngleForce: bond %d [%d %d %d %d] ang=%10.3f k=%10.3f [%10.3e %10.3e %10.3e %10.3e]\n", 
+                             bondIndex, particle1, particle2, particle3, particle4, idealInPlaneAngle, quadraticK, cubicK, quarticK, penticK, sexticK );
+        (void) fflush( log );
+    }
+#endif
+
+    // T   = AD x CD
+    // P   = B + T*delta
+    // AP  = A - P
+    // CP  = A - P
+    // M   = CP x AP
+
+    enum { AD, BD, CD, T, AP, P, CP, M, APxM, CPxM, ADxBD, BDxCD, TxCD, ADxT, dBxAD, CDxdB, LastDeltaAtomIndex };
+ 
+    // AD   0
+    // BD   1
+    // CD   2 
+    //  T   3
+    // AP   4
+    //  P   5
+    // CP   6
+    // M    7
+    // APxM, CPxM, ADxBD, BDxCD, TxCD, ADxT, dBxAD, CDxdB, LastDeltaAtomIndex
+
+    double deltaR[LastDeltaAtomIndex][3];
+    for( int ii = 0; ii < 3; ii++ ){
+        deltaR[AD][ii] = positions[particle1][ii] - positions[particle4][ii];
+        deltaR[BD][ii] = positions[particle2][ii] - positions[particle4][ii];
+        deltaR[CD][ii] = positions[particle3][ii] - positions[particle4][ii];
+    }
+    crossProductVector3( deltaR[AD], deltaR[CD], deltaR[T] );
+ 
+    double rT2     = dotVector3( deltaR[T], deltaR[T] );
+    double delta   = dotVector3( deltaR[T], deltaR[BD] );
+         delta    *= -1.0/rT2;
+ 
+    for( int ii = 0; ii < 3; ii++ ){
+       deltaR[P][ii]  = positions[particle2][ii] + deltaR[T][ii]*delta;
+       deltaR[AP][ii] = positions[particle1][ii] - deltaR[P][ii];
+       deltaR[CP][ii] = positions[particle3][ii] - deltaR[P][ii];
+    }   
+ 
+    double rAp2 = dotVector3( deltaR[AP],  deltaR[AP] );
+    double rCp2 = dotVector3( deltaR[CP],  deltaR[CP] );
+    if( rAp2 <= 0.0 && rCp2 <= 0.0 ){
+#ifdef AMOEBA_DEBUG
+        if( log ){
+            (void) fprintf( log, "computeAmoebaHarmonicInPlaneAngleForce:  rAp2 or rCp2 <= 0.0\n" );
+            (void) fflush( log );
+        }
+#endif
+
+        return;
+    }
+
+    crossProductVector3( deltaR[CP], deltaR[AP], deltaR[M] );
+ 
+    double rm = dotVector3( deltaR[M], deltaR[M] );
+         rm   = sqrt( rm );
+    if( rm < 0.000001 ){
+       rm = 0.000001;
+    }
+ 
+    double dot     = dotVector3( deltaR[AP], deltaR[CP] );
+    double cosine  = dot/sqrt( rAp2*rCp2 );
+ 
+    double dEdR;
+    double energyTerm;
+    getPrefactorsGivenInPlaneAngleCosine( cosine, idealInPlaneAngle, quadraticK, cubicK,
+                                          quarticK, penticK, sexticK, &dEdR,  &energyTerm, log );
+ 
+    double termA   = -dEdR/(rAp2*rm);
+    double termC   =  dEdR/(rCp2*rm);
+ 
+    crossProductVector3( deltaR[AP], deltaR[M], deltaR[APxM] );
+    crossProductVector3( deltaR[CP], deltaR[M], deltaR[CPxM] );
+ 
+    // forces will be gathered here
+ 
+    enum { dA, dB, dC, dD, LastDIndex };
+    double forceTerm[LastDIndex][3];
+ 
+    for( int ii = 0; ii < 3; ii++ ){
+       forceTerm[dA][ii] = deltaR[APxM][ii]*termA;
+       forceTerm[dC][ii] = deltaR[CPxM][ii]*termC;
+       forceTerm[dB][ii] = -1.0*( forceTerm[dA][ii] + forceTerm[dC][ii] );
+    }
+ 
+    double pTrT2  = dotVector3( forceTerm[dB], deltaR[T] );
+         pTrT2   /= rT2;
+ 
+    crossProductVector3( deltaR[CD], forceTerm[dB], deltaR[CDxdB] );
+    crossProductVector3( forceTerm[dB], deltaR[AD], deltaR[dBxAD] );
+ 
+    if( fabs( pTrT2 ) > 1.0e-08 ){
+       double delta2 = delta*2.0;
+ 
+       crossProductVector3( deltaR[BD], deltaR[CD], deltaR[BDxCD] );
+       crossProductVector3( deltaR[T],  deltaR[CD], deltaR[TxCD]  );
+       crossProductVector3( deltaR[AD], deltaR[BD], deltaR[ADxBD] );
+       crossProductVector3( deltaR[AD], deltaR[T],  deltaR[ADxT]  );
+       for( int ii = 0; ii < 3; ii++ ){
+ 
+          double term           = deltaR[BDxCD][ii] + delta2*deltaR[TxCD][ii];
+          forceTerm[dA][ii]  += delta*deltaR[CDxdB][ii] + term*pTrT2;
+ 
+               term           = deltaR[ADxBD][ii] + delta2*deltaR[ADxT][ii];
+          forceTerm[dC][ii]  += delta*deltaR[dBxAD][ii] + term*pTrT2;
+ 
+          forceTerm[dD][ii]  = -( forceTerm[dA][ii] + forceTerm[dB][ii] + forceTerm[dC][ii] );
+       }
+    } else {
+       for( int ii = 0; ii < 3; ii++ ){
+ 
+          forceTerm[dA][ii] += delta*deltaR[CDxdB][ii];
+          forceTerm[dC][ii] += delta*deltaR[dBxAD][ii];
+ 
+          forceTerm[dD][ii]  = -( forceTerm[dA][ii] + forceTerm[dB][ii] + forceTerm[dC][ii] );
+       }
+    }
+ 
+    // accumulate forces and energy
+ 
+    *energy                    += energyTerm;
+ 
+    forces[particle1][0]       -= forceTerm[0][0];
+    forces[particle1][1]       -= forceTerm[0][1];
+    forces[particle1][2]       -= forceTerm[0][2];
+
+    forces[particle2][0]       -= forceTerm[1][0];
+    forces[particle2][1]       -= forceTerm[1][1];
+    forces[particle2][2]       -= forceTerm[1][2];
+
+    forces[particle3][0]       -= forceTerm[2][0];
+    forces[particle3][1]       -= forceTerm[2][1];
+    forces[particle3][2]       -= forceTerm[2][2];
+
+    forces[particle4][0]       -= forceTerm[3][0];
+    forces[particle4][1]       -= forceTerm[3][1];
+    forces[particle4][2]       -= forceTerm[3][2];
+
+}
+
+static void computeAmoebaHarmonicInPlaneAngleForces( Context& context, AmoebaHarmonicInPlaneAngleForce& amoebaHarmonicInPlaneAngleForce,
+                                                     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 < amoebaHarmonicInPlaneAngleForce.getNumAngles(); ii++ ){
+        computeAmoebaHarmonicInPlaneAngleForce(ii, positions, amoebaHarmonicInPlaneAngleForce, expectedForces, expectedEnergy, log );
+    }
+
+#ifdef AMOEBA_DEBUG
+    if( log ){
+        (void) fprintf( log, "computeAmoebaHarmonicInPlaneAngleForces: expected energy=%14.7e\n", *expectedEnergy );
+        for( unsigned int ii = 0; ii < positions.size(); ii++ ){
+            (void) fprintf( log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2] );
+        }
+        (void) fflush( log );
+    }
+#endif
+    return;
+
+}
+
+void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicInPlaneAngleForce& amoebaHarmonicInPlaneAngleForce,
+                                        double tolerance, const std::string& idString, FILE* log) {
+
+    std::vector<Vec3> expectedForces;
+    double expectedEnergy;
+    computeAmoebaHarmonicInPlaneAngleForces( context, amoebaHarmonicInPlaneAngleForce, expectedForces, &expectedEnergy, log );
+   
+    State state                      = context.getState(State::Forces | State::Energy);
+    const std::vector<Vec3> forces   = state.getForces();
+
+#ifdef AMOEBA_DEBUG
+    if( log ){
+        (void) fprintf( log, "computeAmoebaHarmonicInPlaneAngleForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
+        for( unsigned int ii = 0; ii < forces.size(); ii++ ){
+            (void) fprintf( log, "%6u [%14.7e %14.7e %14.7e]   [%14.7e %14.7e %14.7e]\n", ii,
+                            expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
+        }
+        (void) fflush( log );
+    }
+#endif
+
+    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 testOneAngle( FILE* log ) {
+
+    System system;
+    int numberOfParticles = 4;
+    for( int ii = 0; ii < numberOfParticles; ii++ ){
+        system.addParticle(1.0);
+    }
+
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+
+    AmoebaHarmonicInPlaneAngleForce* amoebaHarmonicInPlaneAngleForce = new AmoebaHarmonicInPlaneAngleForce();
+
+    double angle      = 65.0;
+    double quadraticK = 1.0;
+    double cubicK     = 0.0e-01;
+    double quarticK   = 0.0e-02;
+    double penticK    = 0.0e-03;
+    double sexticK    = 0.0e-04;
+    amoebaHarmonicInPlaneAngleForce->addAngle(0, 1, 2, 3, angle, quadraticK);
+
+    amoebaHarmonicInPlaneAngleForce->setAmoebaGlobalHarmonicInPlaneAngleCubic(cubicK);
+    amoebaHarmonicInPlaneAngleForce->setAmoebaGlobalHarmonicInPlaneAngleQuartic(quarticK);
+    amoebaHarmonicInPlaneAngleForce->setAmoebaGlobalHarmonicInPlaneAnglePentic(penticK);
+    amoebaHarmonicInPlaneAngleForce->setAmoebaGlobalHarmonicInPlaneAngleSextic(sexticK);
+
+    system.addForce(amoebaHarmonicInPlaneAngleForce);
+    Context context(system, integrator, Platform::getPlatformByName( "CUDA"));
+
+    std::vector<Vec3> positions(numberOfParticles);
+
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(0, 0, 1);
+    positions[3] = Vec3(1, 1, 1);
+
+    context.setPositions(positions);
+    compareWithExpectedForceAndEnergy( context, *amoebaHarmonicInPlaneAngleForce, TOL, "testOneInPlaneAngle", log );
+
+}
+
+int main( int numberOfArguments, char* argv[] ) {
+
+    try {
+        std::cout << "TestCudaAmoebaHarmonicInPlaneAngleForce running test..." << std::endl;
+        registerAmoebaCudaKernelFactories();
+        FILE* log = NULL;
+        testOneAngle( NULL );
+    } catch(const std::exception& e) {
+        std::cout << "exception: " << e.what() << std::endl;
+        std::cout << "FAIL - ERROR.  Test failed." << std::endl;
+        return 1;
+    }
+    std::cout << "Done" << std::endl;
+    return 0;
+}

plugins/amoeba/platforms/cuda2/tests/TestCudaAmoebaPiTorsionForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMMAmoeba                             *
+ * -------------------------------------------------------------------------- *
+ * 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: Mark Friedrichs                                                   *
+ * 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 CudaAmoebaPiTorsionForce.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "AmoebaTinkerParameterFile.h"
+#include "openmm/Context.h"
+#include "OpenMMAmoeba.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+
+const double TOL = 1e-5;
+#define PI_M               3.141592653589
+#define RADIAN            57.29577951308
+
+/* ---------------------------------------------------------------------------------------
+
+   Compute cross product of two 3-vectors and place in 3rd vector
+
+   vectorZ = vectorX x vectorY
+
+   @param vectorX             x-vector
+   @param vectorY             y-vector
+   @param vectorZ             z-vector
+
+   @return vector is vectorZ
+
+   --------------------------------------------------------------------------------------- */
+     
+static void crossProductVector3( double* vectorX, double* vectorY, double* vectorZ ){
+
+    vectorZ[0]  = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1];
+    vectorZ[1]  = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2];
+    vectorZ[2]  = vectorX[0]*vectorY[1] - vectorX[1]*vectorY[0];
+
+    return;
+}
+
+static double dotVector3( double* vectorX, double* vectorY ){
+    return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2];
+}
+
+
+static void computeAmoebaPiTorsionForce(int bondIndex,  std::vector<Vec3>& positions, AmoebaPiTorsionForce& amoebaPiTorsionForce,
+                                        std::vector<Vec3>& forces, double* energy, FILE* log ) {
+
+    int particle1, particle2, particle3, particle4, particle5, particle6;
+    double kTorsion;
+
+    amoebaPiTorsionForce.getPiTorsionParameters(bondIndex, particle1, particle2, particle3, particle4,  particle5, particle6, kTorsion);
+
+#ifdef AMOEBA_DEBUG
+    if( log ){
+        (void) fprintf( log, "computeAmoebaPiTorsionForce: bond %d [%d %d %d %d %d %d] k=%10.3e\n", 
+                             bondIndex, particle1, particle2, particle3, particle4, particle5, particle6, kTorsion );
+        (void) fflush( log );
+    }
+#endif
+
+    enum { AD, BD, EC, FC, P, Q, CP, DC, QD, T, U, TU, DP, QC, dT, dU, dP, dQ, dC1, dC2, dD1, dD2, LastDeltaIndex };
+    double deltaR[LastDeltaIndex][3];
+
+    enum { A, B, C, D, E, F, LastAtomIndex };
+    double d[LastAtomIndex][3];
+ 
+    for( int ii = 0; ii < 3; ii++ ){
+         deltaR[AD][ii] = positions[particle1][ii] - positions[particle4][ii];
+         deltaR[BD][ii] = positions[particle2][ii] - positions[particle4][ii];
+         deltaR[EC][ii] = positions[particle5][ii] - positions[particle3][ii];
+         deltaR[FC][ii] = positions[particle6][ii] - positions[particle3][ii];
+    }
+ 
+    crossProductVector3( deltaR[AD], deltaR[BD], deltaR[P] );
+    crossProductVector3( deltaR[EC], deltaR[FC], deltaR[Q] );
+    for( int ii = 0; ii < 3; ii++ ){
+        deltaR[CP][ii]  = -deltaR[P][ii];
+        deltaR[DC][ii]  =  positions[particle4][ii] - positions[particle3][ii];
+        deltaR[QD][ii]  =  deltaR[Q][ii];
+  
+        deltaR[P][ii]  += positions[particle3][ii];
+        deltaR[Q][ii]  += positions[particle4][ii];
+    }
+    crossProductVector3( deltaR[CP], deltaR[DC], deltaR[T]  );
+    crossProductVector3( deltaR[DC], deltaR[QD], deltaR[U]  );
+    crossProductVector3( deltaR[T],  deltaR[U],  deltaR[TU] );
+ 
+    double rT2  = dotVector3( deltaR[T], deltaR[T] );
+    double rU2  = dotVector3( deltaR[U], deltaR[U] );
+    double rTrU = sqrt( rT2*rU2 );
+    if( rTrU <= 0.0 ){
+        return;
+    }
+ 
+    double rDC     = dotVector3( deltaR[DC], deltaR[DC] );
+         rDC       = sqrt( rDC );
+  
+    double cosine  = dotVector3( deltaR[T], deltaR[U] );
+         cosine   /= rTrU;
+ 
+    double sine   = dotVector3( deltaR[DC], deltaR[TU] );
+         sine    /= ( rDC*rTrU );
+ 
+    double cosine2 = cosine*cosine - sine*sine;
+    double sine2   = 2.0*cosine*sine;
+ 
+    double phi2    = 1.0 - cosine2;
+    double dphi2   = 2.0*sine2;
+ 
+    double dedphi = kTorsion*dphi2; 
+ 
+    for( int ii = 0; ii < 3; ii++ ){
+        deltaR[DP][ii] = positions[particle4][ii] -    deltaR[P][ii];
+        deltaR[QC][ii] = deltaR[Q][ii]    - positions[particle3][ii];
+    }
+ 
+    double factorT =  dedphi/( rDC*rT2 );
+    double factorU = -dedphi/( rDC*rU2 );
+ 
+    crossProductVector3( deltaR[T], deltaR[DC], deltaR[dT]  );
+    crossProductVector3( deltaR[U], deltaR[DC], deltaR[dU]  );
+    for( int ii = 0; ii < 3; ii++ ){
+        deltaR[dT][ii] *= factorT;
+        deltaR[dU][ii] *= factorU;
+    }
+ 
+    crossProductVector3( deltaR[dT], deltaR[DC], deltaR[dP]   );
+    crossProductVector3( deltaR[dU], deltaR[DC], deltaR[dQ]   );
+ 
+    crossProductVector3( deltaR[DP], deltaR[dT], deltaR[dC1]  );
+    crossProductVector3( deltaR[dU], deltaR[QD], deltaR[dC2]  );
+ 
+    crossProductVector3( deltaR[dT], deltaR[CP], deltaR[dD1]  );
+    crossProductVector3( deltaR[QC], deltaR[dU], deltaR[dD2]  );
+ 
+    crossProductVector3( deltaR[BD], deltaR[dP], d[A]         );
+    crossProductVector3( deltaR[dP], deltaR[AD], d[B]         );
+ 
+    crossProductVector3( deltaR[FC], deltaR[dQ], d[E]         );
+    crossProductVector3( deltaR[dQ], deltaR[EC], d[F]         );
+ 
+    for( int ii = 0; ii < 3; ii++ ){
+        d[C][ii] = deltaR[dC1][ii] + deltaR[dC2][ii] + deltaR[dP][ii] - d[E][ii] - d[F][ii];
+        d[D][ii] = deltaR[dD1][ii] + deltaR[dD2][ii] + deltaR[dQ][ii] - d[A][ii] - d[B][ii];
+    }
+ 
+    // ---------------------------------------------------------------------------------------
+ 
+    // accumulate forces and energy
+ 
+    forces[particle1][0]       -= d[0][0];
+    forces[particle1][1]       -= d[0][1];
+    forces[particle1][2]       -= d[0][2];
+
+    forces[particle2][0]       -= d[1][0];
+    forces[particle2][1]       -= d[1][1];
+    forces[particle2][2]       -= d[1][2];
+
+    forces[particle3][0]       -= d[2][0];
+    forces[particle3][1]       -= d[2][1];
+    forces[particle3][2]       -= d[2][2];
+
+    forces[particle4][0]       -= d[3][0];
+    forces[particle4][1]       -= d[3][1];
+    forces[particle4][2]       -= d[3][2];
+
+    forces[particle5][0]       -= d[4][0];
+    forces[particle5][1]       -= d[4][1];
+    forces[particle5][2]       -= d[4][2];
+
+    forces[particle6][0]       -= d[5][0];
+    forces[particle6][1]       -= d[5][1];
+    forces[particle6][2]       -= d[5][2];
+
+    *energy                    += kTorsion*phi2;
+ 
+    return;
+}
+ 
+static void computeAmoebaPiTorsionForces( Context& context, AmoebaPiTorsionForce& amoebaPiTorsionForce,
+                                          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 < amoebaPiTorsionForce.getNumPiTorsions(); ii++ ){
+        computeAmoebaPiTorsionForce(ii, positions, amoebaPiTorsionForce, expectedForces, expectedEnergy, log );
+    }
+
+#ifdef AMOEBA_DEBUG
+    if( log ){
+        (void) fprintf( log, "computeAmoebaPiTorsionForces: expected energy=%14.7e\n", *expectedEnergy );
+        for( unsigned int ii = 0; ii < positions.size(); ii++ ){
+            (void) fprintf( log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2] );
+        }
+        (void) fflush( log );
+    }
+#endif
+    return;
+
+}
+
+void compareWithExpectedForceAndEnergy( Context& context, AmoebaPiTorsionForce& amoebaPiTorsionForce,
+                                        double tolerance, const std::string& idString, FILE* log) {
+
+    std::vector<Vec3> expectedForces;
+    double expectedEnergy;
+    computeAmoebaPiTorsionForces( context, amoebaPiTorsionForce, expectedForces, &expectedEnergy, log );
+   
+    State state                      = context.getState(State::Forces | State::Energy);
+    const std::vector<Vec3> forces   = state.getForces();
+
+#ifdef AMOEBA_DEBUG
+    if( log ){
+        (void) fprintf( log, "computeAmoebaPiTorsionForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
+        for( unsigned int ii = 0; ii < forces.size(); ii++ ){
+            (void) fprintf( log, "%6u [%14.7e %14.7e %14.7e]   [%14.7e %14.7e %14.7e]\n", ii,
+                            expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
+        }
+        (void) fflush( log );
+    }
+#endif
+
+    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 testOnePiTorsion( FILE* log ) {
+
+    System system;
+    int numberOfParticles = 6;
+    for( int ii = 0; ii < numberOfParticles; ii++ ){
+        system.addParticle(1.0);
+    }
+
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+
+    AmoebaPiTorsionForce* amoebaPiTorsionForce = new AmoebaPiTorsionForce();
+
+    double kTorsion = 6.85;
+    amoebaPiTorsionForce->addPiTorsion(0, 1, 2, 3, 4, 5, kTorsion );
+
+    system.addForce(amoebaPiTorsionForce);
+    Context context(system, integrator, Platform::getPlatformByName( "CUDA"));
+
+    std::vector<Vec3> positions(numberOfParticles);
+
+    positions[0] = Vec3( 0.262660000E+02,  0.254130000E+02,  0.284200000E+01 );
+    positions[1] = Vec3( 0.278860000E+02,  0.264630000E+02,  0.426300000E+01 );
+    positions[2] = Vec3( 0.269130000E+02,  0.266390000E+02,  0.353100000E+01 );
+
+    positions[3] = Vec3( 0.245568230E+02,  0.250215290E+02,  0.796852800E+01 );
+    positions[4] = Vec3( 0.261000000E+02,  0.292530000E+02,  0.520200000E+01 );
+    positions[5] = Vec3( 0.254124630E+02,  0.234691880E+02,  0.773335400E+01 );
+
+    context.setPositions(positions);
+    compareWithExpectedForceAndEnergy( context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion", log );
+
+}
+
+int main( int numberOfArguments, char* argv[] ) {
+
+    try {
+        std::cout << "TestCudaAmoebaPiTorsionForce running test..." << std::endl;
+        registerAmoebaCudaKernelFactories();
+        FILE* log = NULL;
+        testOnePiTorsion( 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 << "Done" << std::endl;
+    return 0;
+}