Initial Amoeba

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

+#include "amoebaCudaKernels.h"
+
+//#define AMOEBA_DEBUG
+#define BLOCK_SIZE 128
+
+using namespace std; 
+
+static __constant__ cudaGmxSimulation cSim;
+static __constant__ cudaAmoebaGmxSimulation cAmoebaSim;
+
+void SetCalculateAmoebaCudaMapTorquesSim(amoebaGpuContext amoebaGpu)
+{
+    cudaError_t status;
+    gpuContext gpu = amoebaGpu->gpuContext;
+    status         = cudaMemcpyToSymbol(cSim, &gpu->sim, sizeof(cudaGmxSimulation));    
+    RTERROR(status, "SetCalculateAmoebaCudaMapTorquesSim: cudaMemcpyToSymbol: SetSim copy to cSim failed");
+    status         = cudaMemcpyToSymbol(cAmoebaSim, &amoebaGpu->amoebaSim, sizeof(cudaAmoebaGmxSimulation));    
+    RTERROR(status, "SetCalculateAmoebaCudaMapTorquesSim: cudaMemcpyToSymbol: SetSim copy to cAmoebaSim failed");
+}
+
+void GetCalculateAmoebaCudaMapTorquesSim(amoebaGpuContext amoebaGpu)
+{
+    cudaError_t status;
+    gpuContext gpu = amoebaGpu->gpuContext;
+    status = cudaMemcpyFromSymbol(&gpu->sim, cSim, sizeof(cudaGmxSimulation));    
+    RTERROR(status, "GetCalculateAmoebaCudaMapTorquesSim: cudaMemcpyFromSymbol: SetSim copy from cSim failed");
+    status = cudaMemcpyFromSymbol(&amoebaGpu->amoebaSim, cAmoebaSim, sizeof(cudaAmoebaGmxSimulation));    
+    RTERROR(status, "GetCalculateAmoebaCudaMapTorquesSim: cudaMemcpyFromSymbol: SetSim copy from cAmoebaSim failed");
+}
+
+__device__ float normVector3( float* vector ) 
+{
+
+    float norm                    = DOT3( vector, vector );
+    float returnNorm              = SQRT( norm );
+    norm                          = returnNorm > 0.0f ? 1.0f/returnNorm : 0.0f;
+
+    vector[0]                    *= norm;
+    vector[1]                    *= norm;
+    vector[2]                    *= norm;
+
+    return returnNorm;
+}
+
+__global__
+#if (__CUDA_ARCH__ >= 200)
+__launch_bounds__(GF1XX_THREADS_PER_BLOCK, 1)
+#elif (__CUDA_ARCH__ >= 130)
+__launch_bounds__(GT2XX_THREADS_PER_BLOCK, 1)
+#else
+__launch_bounds__(G8X_THREADS_PER_BLOCK, 1)
+#endif
+void amoebaMapTorqueToForce_kernel( 
+				    int numOfAtoms,
+				    float4* atomCoord,
+				    float* torque,
+				    int4* multiPoleAtoms,
+				    int maxDiff,
+				    float* tempElecForce
+				    ){
+
+  // ---------------------------------------------------------------------------------------
+  
+  int threadId = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
+  if(threadId >= numOfAtoms)return;
+  
+  int U     = 0;
+  int V     = 1;
+  int W     = 2;
+  
+  int X     = 0;
+  int Y     = 1;
+  int Z     = 2;
+  
+  float forces[3][3];
+  float norms[3];
+  float vector[3][3];
+  
+  // ---------------------------------------------------------------------------------------
+  
+    int axisAtom                  = multiPoleAtoms[threadId].x;
+
+    vector[U][0]                  = atomCoord[axisAtom].x - atomCoord[threadId].x;
+    vector[U][1]                  = atomCoord[axisAtom].y - atomCoord[threadId].y;
+    vector[U][2]                  = atomCoord[axisAtom].z - atomCoord[threadId].z;
+
+    norms[U]                      = normVector3( vector[U] );
+
+    axisAtom                      = multiPoleAtoms[threadId].y;
+
+    vector[V][0]                  = atomCoord[axisAtom].x - atomCoord[threadId].x;
+    vector[V][1]                  = atomCoord[axisAtom].y - atomCoord[threadId].y;
+    vector[V][2]                  = atomCoord[axisAtom].z - atomCoord[threadId].z;
+
+    norms[V]                      = normVector3( vector[V] );
+
+    // W = UxV
+
+    vector[W][0]                  = vector[U][1]*vector[V][2] - vector[U][2]*vector[V][1];
+    vector[W][1]                  = vector[U][2]*vector[V][0] - vector[U][0]*vector[V][2];
+    vector[W][2]                  = vector[U][0]*vector[V][1] - vector[U][1]*vector[V][0];
+
+    norms[W]                      = normVector3( vector[W] );
+
+    float diff[3];
+    diff[0]                       = vector[V][0] - vector[U][0];
+    diff[1]                       = vector[V][1] - vector[U][1];
+    diff[2]                       = vector[V][2] - vector[U][2];
+
+    float dotDu                   = DOT3( vector[U], diff );
+    float dotDv                   = DOT3( vector[V], diff );
+
+    float up[3], vp[3];
+
+    up[0]                         = diff[0] - dotDu*vector[U][0];
+    vp[0]                         = diff[0] - dotDv*vector[V][0];
+
+    up[1]                         = diff[1] - dotDu*vector[U][1];
+    vp[1]                         = diff[1] - dotDv*vector[V][1];
+
+    up[2]                         = diff[2] - dotDu*vector[U][2];
+    vp[2]                         = diff[2] - dotDv*vector[V][2];
+
+    float norm                    = normVector3( up );
+    norm                          = normVector3( vp );
+
+    float dphi[3];
+    dphi[0]                       = vector[0][0]*torque[threadId*3] + vector[0][1]*torque[threadId*3+1] + vector[0][2]*torque[threadId*3+2];
+    dphi[1]                       = vector[1][0]*torque[threadId*3] + vector[1][1]*torque[threadId*3+1] + vector[1][2]*torque[threadId*3+2];
+    dphi[2]                       = vector[2][0]*torque[threadId*3] + vector[2][1]*torque[threadId*3+1] + vector[2][2]*torque[threadId*3+2];
+
+    // clamp c to interval [-1,1]
+
+    float c                       = DOT3( vector[U], vector[V] );
+          c                       = c >  1.0f ?  1.0f : c;
+          c                       = c < -1.0f ? -1.0f : c;
+
+    float s                       = SQRT( 1.0f - (c*c) );
+    float uvdis                   = norms[U]*s;
+    float vudis                   = norms[V]*s;
+
+    float factorX;
+    float factorZ;
+    if( multiPoleAtoms[threadId].w == 1 ){       
+        factorX = 0.5f;
+        factorZ = 0.5f;
+    } else {
+        factorX = 1.0f;
+        factorZ = 0.0f;
+    }
+
+    forces[X][0]                  = -vector[W][0]*dphi[V]/uvdis + factorX*up[0]*dphi[W]/norms[U];
+    forces[Z][0]                  =  vector[W][0]*dphi[U]/vudis + factorZ*vp[0]*dphi[W]/norms[V];
+
+    forces[X][1]                  = -vector[W][1]*dphi[V]/uvdis + factorX*up[1]*dphi[W]/norms[U];
+    forces[Z][1]                  =  vector[W][1]*dphi[U]/vudis + factorZ*vp[1]*dphi[W]/norms[V];
+
+    forces[X][2]                  = -vector[W][2]*dphi[V]/uvdis + factorX*up[2]*dphi[W]/norms[U];
+    forces[Z][2]                  =  vector[W][2]*dphi[U]/vudis + factorZ*vp[2]*dphi[W]/norms[V];
+
+    forces[Y][0]                  = -(forces[X][0] + forces[Z][0]);
+    forces[Y][1]                  = -(forces[X][1] + forces[Z][1]);
+    forces[Y][2]                  = -(forces[X][2] + forces[Z][2]);
+
+    int temp                      = multiPoleAtoms[threadId].x;
+    int min                       = multiPoleAtoms[temp].z;
+    int offset                    = 3*(temp*maxDiff + threadId-min);
+    tempElecForce[offset  ]       = forces[X][0];
+    tempElecForce[offset+1]       = forces[X][1];
+    tempElecForce[offset+2]       = forces[X][2];
+   
+    temp                          = multiPoleAtoms[threadId].y;
+    min                           = multiPoleAtoms[temp].z;
+    offset                        = 3*(temp*maxDiff + threadId-min);
+    tempElecForce[offset  ]       = forces[Z][0];
+    tempElecForce[offset+1]       = forces[Z][1];
+    tempElecForce[offset+2]       = forces[Z][2];
+ 
+    min                           = multiPoleAtoms[threadId].z;
+    offset                        = 3*(threadId*(maxDiff + 1)-min); 
+    tempElecForce[offset  ]       = forces[Y][0];
+    tempElecForce[offset+1]       = forces[Y][1];
+    tempElecForce[offset+2]       = forces[Y][2];
+   
+}
+
+__global__
+#if (__CUDA_ARCH__ >= 200)
+__launch_bounds__(GF1XX_THREADS_PER_BLOCK, 1)
+#elif (__CUDA_ARCH__ >= 130)
+__launch_bounds__(GT2XX_THREADS_PER_BLOCK, 1)
+#else
+__launch_bounds__(G8X_THREADS_PER_BLOCK, 1)
+#endif
+void amoebaMapTorqueReduce_kernel( 
+					   int numThreads,
+					   int numOfAtoms,
+					   int maxDiff,
+					   float* tempElecForce,
+					   float* elecForce
+					   ){
+  
+    unsigned int tid = threadIdx.x;
+    
+    __shared__ float sfx[BLOCK_SIZE];
+    __shared__ float sfy[BLOCK_SIZE];
+    __shared__ float sfz[BLOCK_SIZE];
+  
+    // load values then sum and add results to elecForce
+
+    if( tid < maxDiff ){
+  
+        sfx[tid] = tempElecForce[blockIdx.x*3*maxDiff + tid*3  ];
+        sfy[tid] = tempElecForce[blockIdx.x*3*maxDiff + tid*3+1];
+        sfz[tid] = tempElecForce[blockIdx.x*3*maxDiff + tid*3+2];
+  
+    } else {
+        sfx[tid] = sfy[tid] = sfz[tid] = 0.0f;
+    }
+  
+    __syncthreads();
+  
+    for( unsigned int s = (blockDim.x)/2; s != 0; s >>= 1 ){
+        if(tid<s){
+            sfx[tid] += sfx[tid+s];
+            sfy[tid] += sfy[tid+s];
+            sfz[tid] += sfz[tid+s];
+        }
+        __syncthreads();
+    }
+  
+    if( tid == 0 ){
+        elecForce[blockIdx.x*3]    += sfx[0];
+        elecForce[blockIdx.x*3+1]  += sfy[0];
+        elecForce[blockIdx.x*3+2]  += sfz[0];
+    }  
+
+}
+
+__global__
+#if (__CUDA_ARCH__ >= 200)
+__launch_bounds__(GF1XX_THREADS_PER_BLOCK, 1)
+#elif (__CUDA_ARCH__ >= 130)
+__launch_bounds__(GT2XX_THREADS_PER_BLOCK, 1)
+#else
+__launch_bounds__(G8X_THREADS_PER_BLOCK, 1)
+#endif
+void amoebaMapTorqueReduce_kernel2( 
+					   int numThreads,
+					   int numOfAtoms,
+					   int maxDiff,
+					   float* tempElecForce,
+					   float* elecForce,
+					   float4* outputForce
+					   ){
+  
+    unsigned int tid = threadIdx.x;
+    
+    __shared__ float sfx[BLOCK_SIZE];
+    __shared__ float sfy[BLOCK_SIZE];
+    __shared__ float sfz[BLOCK_SIZE];
+  
+    // load values then sum and add results to elecForce
+
+    if( tid < maxDiff ){
+  
+        sfx[tid] = tempElecForce[blockIdx.x*3*maxDiff + tid*3  ];
+        sfy[tid] = tempElecForce[blockIdx.x*3*maxDiff + tid*3+1];
+        sfz[tid] = tempElecForce[blockIdx.x*3*maxDiff + tid*3+2];
+  
+    } else {
+        sfx[tid] = sfy[tid] = sfz[tid] = 0.0f;
+    }
+  
+    __syncthreads();
+  
+    for( unsigned int s = (blockDim.x)/2; s != 0; s >>= 1 ){
+        if(tid<s){
+            sfx[tid] += sfx[tid+s];
+            sfy[tid] += sfy[tid+s];
+            sfz[tid] += sfz[tid+s];
+        }
+        __syncthreads();
+    }
+  
+    if( tid == 0 ){
+        outputForce[blockIdx.x].x  += elecForce[blockIdx.x*3  ] + sfx[0];
+        outputForce[blockIdx.x].y  += elecForce[blockIdx.x*3+1] + sfy[0];
+        outputForce[blockIdx.x].z  += elecForce[blockIdx.x*3+2] + sfz[0];
+    }  
+
+}
+
+void cudaComputeAmoebaMapTorques( amoebaGpuContext amoebaGpu, CUDAStream<float>* psTorque, CUDAStream<float>* psForce )
+{
+  
+   // ---------------------------------------------------------------------------------------
+
+//#define AMOEBA_DEBUG
+#ifdef AMOEBA_DEBUG
+    static const char* methodName = "cudaMapAmoebaTorqueToForce";
+    static int timestep = 0; 
+    std::vector<int> fileId;
+    timestep++;
+    fileId.resize( 2 ); 
+    fileId[0] = timestep;
+    fileId[1] = 1; 
+#endif
+
+    // check that BLOCK_SIZE >= amoebaGpu->maxMapTorqueDifference
+
+    if( amoebaGpu->maxMapTorqueDifference > BLOCK_SIZE ){
+        (void) fprintf( amoebaGpu->log, "block size (%d) in amoebaMapTorqueReduce_kernel is too small ( > %d)! -- aborting.\n",
+                        BLOCK_SIZE, amoebaGpu->maxMapTorqueDifference );
+        exit(-1);  
+    }
+
+   // ---------------------------------------------------------------------------------------
+    
+    gpuContext gpu    = amoebaGpu->gpuContext;
+    int numThreads    = min(256, (gpu->natoms));
+    int numBlocks     =  1 + (gpu->natoms/numThreads);
+
+//#ifdef AMOEBA_DEBUG
+#if 0
+    if( amoebaGpu->log ){
+        (void) fprintf( amoebaGpu->log, "%s: numBlocks=%d numThreads=%d\n", methodName, numBlocks, numThreads ); (void) fflush( amoebaGpu->log );
+    }
+    amoebaGpu->psForce->Download();
+    amoebaGpu->psTorque->Download();
+    int maxPrint        = 20;
+    (void) fprintf( amoebaGpu->log,"Pre torqueMap\n" );
+    for( int ii = 0; ii < gpu->natoms; ii++ ){
+        (void) fprintf( amoebaGpu->log, "%5d ", ii);
+
+        int indexOffset     = ii*3;
+
+        (void) fprintf( amoebaGpu->log,"E[%16.9e %16.9e %16.9e] ",
+                        amoebaGpu->psForce->_pSysStream[0][indexOffset],
+                        amoebaGpu->psForce->_pSysStream[0][indexOffset+1],
+                        amoebaGpu->psForce->_pSysStream[0][indexOffset+2] );
+        (void) fprintf( amoebaGpu->log,"T[%16.9e %16.9e %16.9e]\n",
+                        amoebaGpu->psTorque->_pSysStream[0][indexOffset],
+                        amoebaGpu->psTorque->_pSysStream[0][indexOffset+1],
+                        amoebaGpu->psTorque->_pSysStream[0][indexOffset+2] );
+        if( ii == maxPrint && (gpu->natoms - maxPrint) > ii )ii = gpu->natoms - maxPrint;
+    }
+    int nansDetected   = checkForNansAndInfinities( gpu->natoms*3, amoebaGpu->psForce );
+    nansDetected      += checkForNansAndInfinities( gpu->natoms*3, amoebaGpu->psTorque );
+    if( nansDetected ){
+        (void) fprintf( amoebaGpu->log,"WARNING: %d nans/infinities detected force/torques.\n", nansDetected );
+    } else {
+        (void) fprintf( amoebaGpu->log,"No nans/infinities detected in force/torques.\n" );
+    }
+
+// zero forces
+#if 0
+    for( int ii = 0; ii < 3*gpu->natoms; ii++ ){
+        amoebaGpu->psForce->_pSysStream[0][ii] = 0.0f;
+    }
+    amoebaGpu->psForce->Upload();
+#endif
+
+#endif
+
+    // torqueMapForce is zeroed when initialized; should not need to be reinitialized 
+
+/*
+    AmoebaTorqueMapZeroKernel<<< numBlocks, numThreads >>>( 
+       gpu->natoms, amoebaGpu->torqueMapForce->_pDevStream[0] );
+    LAUNCHERROR("AmoebaMapTrqZeroKernel");  
+*/
+
+
+    amoebaMapTorqueToForce_kernel<<< numBlocks, numThreads>>> (
+                gpu->natoms,
+                gpu->psPosq4->_pDevStream[0],
+                psTorque->_pDevStream[0],
+                amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pDevStream[0],
+                amoebaGpu->maxMapTorqueDifference,
+                amoebaGpu->torqueMapForce->_pDevStream[0] );
+    LAUNCHERROR("AmoebaMapTrqKernel");
+  
+//#ifdef AMOEBA_DEBUG
+#if 0
+    amoebaGpu->torqueMapForce->Download();
+    //int maxPrint        = 10;
+    (void) fprintf( amoebaGpu->log,"Post AmoebaMapTrqKernel maxMapTorqueDifference=%d\n", amoebaGpu->maxMapTorqueDifference );
+    for( int ii = 0; ii < gpu->natoms; ii++ ){
+        (void) fprintf( amoebaGpu->log, "\n%5d multi[%d %d %d %d]\n", ii, 
+                        amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].x,
+                        amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].y,
+                        amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].z,
+                        amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].w );
+
+        int indexOffset     = ii*3*amoebaGpu->maxMapTorqueDifference;
+        float sum[3]        = { 0.0f,  0.0f,  0.0f };
+        for( int jj = 0; jj < amoebaGpu->maxMapTorqueDifference; jj++ ){
+            if( amoebaGpu->torqueMapForce->_pSysStream[0][indexOffset] != 0.0f ){ 
+                (void) fprintf( amoebaGpu->log,"    %4d %4d Temp[%16.9e %16.9e %16.9e] %d\n",
+                                ii, jj + amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].z,
+                                amoebaGpu->torqueMapForce->_pSysStream[0][indexOffset],
+                                amoebaGpu->torqueMapForce->_pSysStream[0][indexOffset+1],
+                                amoebaGpu->torqueMapForce->_pSysStream[0][indexOffset+2], indexOffset );
+                sum[0] += amoebaGpu->torqueMapForce->_pSysStream[0][indexOffset];
+                sum[1] += amoebaGpu->torqueMapForce->_pSysStream[0][indexOffset+1];
+                sum[2] += amoebaGpu->torqueMapForce->_pSysStream[0][indexOffset+2];
+            }
+            indexOffset += 3;
+        }
+        (void) fprintf( amoebaGpu->log,"               Sum[%16.9e %16.9e %16.9e]\n", sum[0], sum[1], sum[2] );
+        if( ii == maxPrint && (gpu->natoms - maxPrint) > ii )ii = gpu->natoms - maxPrint;
+    }
+#endif
+
+    numBlocks  = gpu->natoms;
+    numThreads = amoebaGpu->maxMapTorqueDifferencePow2;
+
+    amoebaMapTorqueReduce_kernel<<< numBlocks, numThreads>>>(
+            numThreads, gpu->natoms,
+            amoebaGpu->maxMapTorqueDifference,
+            amoebaGpu->torqueMapForce->_pDevStream[0],
+            psForce->_pDevStream[0] );
+    LAUNCHERROR("amoebaMapTorqueReduce_kernel");
+
+#ifdef AMOEBA_DEBUG
+    if( 0 && amoebaGpu->log ){
+        (void) fprintf( amoebaGpu->log, "%s: numBlocks=%d numThreads=%d %d\n", methodName, numBlocks, numThreads, amoebaGpu->maxMapTorqueDifferencePow2); (void) fflush( amoebaGpu->log );
+        amoebaGpu->psForce->Download();
+        amoebaGpu->psTorque->Download();
+        int maxPrint        = 10;
+        (void) fprintf( amoebaGpu->log,"Post torqueMap\n" );
+        for( int ii = 0; ii < gpu->natoms; ii++ ){
+            (void) fprintf( amoebaGpu->log, "%5d ", ii);
+    
+            int indexOffset     = ii*3;
+    
+            (void) fprintf( amoebaGpu->log,"E[%16.9e %16.9e %16.9e] ",
+                            amoebaGpu->psForce->_pSysStream[0][indexOffset],
+                            amoebaGpu->psForce->_pSysStream[0][indexOffset+1],
+                            amoebaGpu->psForce->_pSysStream[0][indexOffset+2] );
+            (void) fprintf( amoebaGpu->log,"T[%16.9e %16.9e %16.9e]\n",
+                            amoebaGpu->psTorque->_pSysStream[0][indexOffset],
+                            amoebaGpu->psTorque->_pSysStream[0][indexOffset+1],
+                            amoebaGpu->psTorque->_pSysStream[0][indexOffset+2] );
+            if( ii == maxPrint && (gpu->natoms - maxPrint) > ii )ii = gpu->natoms - maxPrint;
+        }
+        (void) fflush( amoebaGpu->log );
+    }
+    if( 1 ){
+        //std::vector<int> fileId;
+        //fileId.push_back( 0 );
+        VectorOfDoubleVectors outputVector;
+        cudaLoadCudaFloat4Array( gpu->natoms, 3, gpu->psPosq4,              outputVector );
+        cudaLoadCudaFloatArray( gpu->natoms,  3, amoebaGpu->psForce,      outputVector );
+        cudaLoadCudaFloatArray( gpu->natoms,  3, amoebaGpu->psTorque, outputVector);
+        cudaWriteVectorOfDoubleVectorsToFile( "CudaVacuumElecForce", fileId, outputVector );
+    }
+#endif
+
+}
+
+void cudaComputeAmoebaMapTorquesAndAddTotalForce( amoebaGpuContext amoebaGpu,
+                                                  CUDAStream<float>* psTorque,
+                                                  CUDAStream<float>* psForce,
+                                                  CUDAStream<float4>* psCudaForce4 )
+{
+  
+   // ---------------------------------------------------------------------------------------
+
+//#define AMOEBA_DEBUG
+#ifdef AMOEBA_DEBUG
+    static const char* methodName = "cudaComputeAmoebaMapTorquesAndAddTotalForce";
+    static int timestep = 0; 
+    std::vector<int> fileId;
+    timestep++;
+    fileId.resize( 2 ); 
+    fileId[0] = timestep;
+    fileId[1] = 1; 
+#endif
+
+    // check that BLOCK_SIZE >= amoebaGpu->maxMapTorqueDifference
+
+    if( amoebaGpu->maxMapTorqueDifference > BLOCK_SIZE ){
+        (void) fprintf( amoebaGpu->log, "block size (%d) in amoebaMapTorqueReduce_kernel is too small ( > %d)! -- aborting.\n",
+                        BLOCK_SIZE, amoebaGpu->maxMapTorqueDifference );
+        exit(-1);  
+    }
+
+   // ---------------------------------------------------------------------------------------
+    
+    gpuContext gpu    = amoebaGpu->gpuContext;
+    int numThreads    = min(256, (gpu->natoms));
+    int numBlocks     =  1 + (gpu->natoms/numThreads);
+
+//#ifdef AMOEBA_DEBUG
+#if 0
+    if( amoebaGpu->log ){
+        (void) fprintf( amoebaGpu->log, "%s: numBlocks=%d numThreads=%d\n", methodName, numBlocks, numThreads ); (void) fflush( amoebaGpu->log );
+    }
+    amoebaGpu->psForce->Download();
+    amoebaGpu->psTorque->Download();
+    int maxPrint        = 20;
+    (void) fprintf( amoebaGpu->log,"Pre torqueMap\n" );
+    for( int ii = 0; ii < gpu->natoms; ii++ ){
+        (void) fprintf( amoebaGpu->log, "%5d ", ii);
+
+        int indexOffset     = ii*3;
+
+        (void) fprintf( amoebaGpu->log,"E[%16.9e %16.9e %16.9e] ",
+                        amoebaGpu->psForce->_pSysStream[0][indexOffset],
+                        amoebaGpu->psForce->_pSysStream[0][indexOffset+1],
+                        amoebaGpu->psForce->_pSysStream[0][indexOffset+2] );
+        (void) fprintf( amoebaGpu->log,"T[%16.9e %16.9e %16.9e]\n",
+                        amoebaGpu->psTorque->_pSysStream[0][indexOffset],
+                        amoebaGpu->psTorque->_pSysStream[0][indexOffset+1],
+                        amoebaGpu->psTorque->_pSysStream[0][indexOffset+2] );
+        if( ii == maxPrint && (gpu->natoms - maxPrint) > ii )ii = gpu->natoms - maxPrint;
+    }
+    int nansDetected   = checkForNansAndInfinities( gpu->natoms*3, amoebaGpu->psForce );
+    nansDetected      += checkForNansAndInfinities( gpu->natoms*3, amoebaGpu->psTorque );
+    if( nansDetected ){
+        (void) fprintf( amoebaGpu->log,"WARNING: %d nans/infinities detected force/torques.\n", nansDetected );
+    } else {
+        (void) fprintf( amoebaGpu->log,"No nans/infinities detected in force/torques.\n" );
+    }
+
+// zero forces
+#if 0
+    for( int ii = 0; ii < 3*gpu->natoms; ii++ ){
+        amoebaGpu->psForce->_pSysStream[0][ii] = 0.0f;
+    }
+    amoebaGpu->psForce->Upload();
+#endif
+
+#endif
+
+    // torqueMapForce is zeroed when initialized; should not need to be reinitialized 
+
+/*
+    AmoebaTorqueMapZeroKernel<<< numBlocks, numThreads >>>( 
+       gpu->natoms, amoebaGpu->torqueMapForce->_pDevStream[0] );
+    LAUNCHERROR("AmoebaMapTrqZeroKernel");  
+*/
+
+
+    amoebaMapTorqueToForce_kernel<<< numBlocks, numThreads>>> (
+                gpu->natoms,
+                gpu->psPosq4->_pDevStream[0],
+                psTorque->_pDevStream[0],
+                amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pDevStream[0],
+                amoebaGpu->maxMapTorqueDifference,
+                amoebaGpu->torqueMapForce->_pDevStream[0] );
+    LAUNCHERROR("AmoebaMapTrqKernel");
+  
+//#ifdef AMOEBA_DEBUG
+#if 0
+    amoebaGpu->torqueMapForce->Download();
+    //int maxPrint        = 10;
+    (void) fprintf( amoebaGpu->log,"Post AmoebaMapTrqKernel maxMapTorqueDifference=%d\n", amoebaGpu->maxMapTorqueDifference );
+    for( int ii = 0; ii < gpu->natoms; ii++ ){
+        (void) fprintf( amoebaGpu->log, "\n%5d multi[%d %d %d %d]\n", ii, 
+                        amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].x,
+                        amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].y,
+                        amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].z,
+                        amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].w );
+
+        int indexOffset     = ii*3*amoebaGpu->maxMapTorqueDifference;
+        float sum[3]        = { 0.0f,  0.0f,  0.0f };
+        for( int jj = 0; jj < amoebaGpu->maxMapTorqueDifference; jj++ ){
+            if( amoebaGpu->torqueMapForce->_pSysStream[0][indexOffset] != 0.0f ){ 
+                (void) fprintf( amoebaGpu->log,"    %4d %4d Temp[%16.9e %16.9e %16.9e] %d\n",
+                                ii, jj + amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].z,
+                                amoebaGpu->torqueMapForce->_pSysStream[0][indexOffset],
+                                amoebaGpu->torqueMapForce->_pSysStream[0][indexOffset+1],
+                                amoebaGpu->torqueMapForce->_pSysStream[0][indexOffset+2], indexOffset );
+                sum[0] += amoebaGpu->torqueMapForce->_pSysStream[0][indexOffset];
+                sum[1] += amoebaGpu->torqueMapForce->_pSysStream[0][indexOffset+1];
+                sum[2] += amoebaGpu->torqueMapForce->_pSysStream[0][indexOffset+2];
+            }
+            indexOffset += 3;
+        }
+        (void) fprintf( amoebaGpu->log,"               Sum[%16.9e %16.9e %16.9e]\n", sum[0], sum[1], sum[2] );
+        if( ii == maxPrint && (gpu->natoms - maxPrint) > ii )ii = gpu->natoms - maxPrint;
+    }
+#endif
+
+    numBlocks  = gpu->natoms;
+    numThreads = amoebaGpu->maxMapTorqueDifferencePow2;
+
+    amoebaMapTorqueReduce_kernel2<<< numBlocks, numThreads>>>(
+            numThreads, gpu->natoms,
+            amoebaGpu->maxMapTorqueDifference,
+            amoebaGpu->torqueMapForce->_pDevStream[0],
+            psForce->_pDevStream[0], psCudaForce4->_pDevStream[0] );
+    LAUNCHERROR("amoebaMapTorqueReduce_kernel2");
+
+#ifdef AMOEBA_DEBUG
+    if( amoebaGpu->log ){
+        (void) fprintf( amoebaGpu->log, "%s: numBlocks=%d numThreads=%d %d\n", methodName, numBlocks, numThreads, amoebaGpu->maxMapTorqueDifferencePow2); (void) fflush( amoebaGpu->log );
+        amoebaGpu->psForce->Download();
+        psCudaForce4->Download();
+ 
+        amoebaGpu->psTorque->Download();
+        int maxPrint        = 10;
+        (void) fprintf( amoebaGpu->log,"Post torqueMap\n" );
+        for( int ii = 0; ii < gpu->natoms; ii++ ){
+            (void) fprintf( amoebaGpu->log, "%5d ", ii);
+    
+            int indexOffset     = ii*3;
+    
+            (void) fprintf( amoebaGpu->log,"FTtl[%16.9e %16.9e %16.9e] ",
+                            psCudaForce4->_pSysStream[0][ii].x,
+                            psCudaForce4->_pSysStream[0][ii].y,
+                            psCudaForce4->_pSysStream[0][ii].z );
+            (void) fprintf( amoebaGpu->log,"F[%16.9e %16.9e %16.9e] ",
+                            amoebaGpu->psForce->_pSysStream[0][indexOffset],
+                            amoebaGpu->psForce->_pSysStream[0][indexOffset+1],
+                            amoebaGpu->psForce->_pSysStream[0][indexOffset+2] );
+            (void) fprintf( amoebaGpu->log,"T[%16.9e %16.9e %16.9e]\n",
+                            amoebaGpu->psTorque->_pSysStream[0][indexOffset],
+                            amoebaGpu->psTorque->_pSysStream[0][indexOffset+1],
+                            amoebaGpu->psTorque->_pSysStream[0][indexOffset+2] );
+            if( ii == maxPrint && (gpu->natoms - maxPrint) > ii )ii = gpu->natoms - maxPrint;
+        }
+        (void) fflush( amoebaGpu->log );
+    }
+    if( 1 ){
+        //std::vector<int> fileId;
+        //fileId.push_back( 0 );
+        VectorOfDoubleVectors outputVector;
+        //cudaLoadCudaFloat4Array( gpu->natoms, 3, gpu->psPosq4,              outputVector );
+        cudaLoadCudaFloat4Array( gpu->natoms, 4, gpu->psForce4,             outputVector );
+        cudaLoadCudaFloatArray( gpu->natoms,  3, amoebaGpu->psForce,        outputVector );
+        cudaLoadCudaFloatArray( gpu->natoms,  3, amoebaGpu->psTorque,       outputVector);
+        cudaWriteVectorOfDoubleVectorsToFile( "CudaVacuumElecForce", fileId, outputVector );
+    }
+#endif
+
+}

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

+//-----------------------------------------------------------------------------------------
+
+//-----------------------------------------------------------------------------------------
+
+#include "amoebaCudaKernels.h"
+
+#include <stdio.h>
+#include <cuda.h>
+#include <cstdlib>
+using namespace std; 
+
+#define SQRT sqrtf
+
+static __constant__ cudaGmxSimulation cSim;
+static __constant__ cudaAmoebaGmxSimulation cAmoebaSim;
+
+void SetCalculateAmoebaMultipoleForcesSim(amoebaGpuContext amoebaGpu)
+{
+    cudaError_t status;
+    gpuContext gpu = amoebaGpu->gpuContext;
+    status         = cudaMemcpyToSymbol(cSim, &gpu->sim, sizeof(cudaGmxSimulation));     
+    RTERROR(status, "SetCalculateAmoebaMultipoleForcesSim: cudaMemcpyToSymbol: SetSim copy to cSim failed");
+    status         = cudaMemcpyToSymbol(cAmoebaSim, &amoebaGpu->amoebaSim, sizeof(cudaAmoebaGmxSimulation));     
+    RTERROR(status, "SetCalculateAmoebaMultipoleForcesSim: cudaMemcpyToSymbol: SetSim copy to cAmoebaSim failed");
+}
+
+void GetCalculateAmoebaMultipoleForcesSim(amoebaGpuContext amoebaGpu)
+{
+    cudaError_t status;
+    gpuContext gpu = amoebaGpu->gpuContext;
+    status = cudaMemcpyFromSymbol(&gpu->sim, cSim, sizeof(cudaGmxSimulation));     
+    RTERROR(status, "GetCalculateAmoebaMultipoleForcesSim: cudaMemcpyFromSymbol: SetSim copy from cSim failed");
+    status = cudaMemcpyFromSymbol(&amoebaGpu->amoebaSim, cAmoebaSim, sizeof(cudaAmoebaGmxSimulation));     
+    RTERROR(status, "GetCalculateAmoebaMultipoleForcesSim: cudaMemcpyFromSymbol: SetSim copy from cAmoebaSim failed");
+}
+
+__device__ float normVector3( float* vector ) 
+{
+
+    float norm                    = DOT3( vector, vector );
+    float returnNorm              = SQRT( norm );
+    norm                          = returnNorm > 0.0f ? 1.0f/returnNorm : 0.0f;
+
+    vector[0]                    *= norm;
+    vector[1]                    *= norm;
+    vector[2]                    *= norm;
+
+    return returnNorm;
+}
+
+__global__
+#if (__CUDA_ARCH__ >= 200)
+__launch_bounds__(GF1XX_THREADS_PER_BLOCK, 1)
+#elif (__CUDA_ARCH__ >= 130)
+__launch_bounds__(GT2XX_THREADS_PER_BLOCK, 1)
+#else
+__launch_bounds__(G8X_THREADS_PER_BLOCK, 1)
+#endif
+void kCudaComputeLabFrameMoments_kernel(
+				   int numOfAtoms,
+				   float *rotationMatrix,
+				   float4 *atomCoord,
+				   int4 *multiPoleAtoms,
+				   float *molecularDipole, float *molecularQuadrupole,
+				   float *labFrameDipole,  float *labFrameQuadrupole )
+{
+
+   float* vectorX;
+   float* vectorY;
+   float* vectorZ;
+
+   // ---------------------------------------------------------------------------------------
+
+   int atomIndex = blockIdx.x;//__mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
+
+   // ---------------------------------------------------------------------------------------
+
+   // get coordinates of this atom and the z & x axis atoms
+   // compute the vector between the atoms and 1/sqrt(d2), d2 is distance between
+   // this atom and the axis atom
+
+   // this atom is referred to as the k-atom in notes below
+
+   // code common to ZThenX and Bisector
+   
+   vectorX                          = &(rotationMatrix[atomIndex*9]);
+   vectorY                          = &(rotationMatrix[atomIndex*9+ 3]);
+   vectorZ                          = &(rotationMatrix[atomIndex*9+ 6]);
+
+   float4 coordinatesThisAtom       = atomCoord[atomIndex];
+
+   int multipoleAtomIndex           = multiPoleAtoms[atomIndex].x;
+   float4 coordinatesAxisAtom       = atomCoord[multipoleAtomIndex];
+
+   vectorZ[0]                       = coordinatesAxisAtom.x - coordinatesThisAtom.x;
+   vectorZ[1]                       = coordinatesAxisAtom.y - coordinatesThisAtom.y;
+   vectorZ[2]                       = coordinatesAxisAtom.z - coordinatesThisAtom.z;
+     
+   multipoleAtomIndex               = multiPoleAtoms[atomIndex].y; 
+   coordinatesAxisAtom              = atomCoord[multipoleAtomIndex];
+
+   vectorX[0]                       = coordinatesAxisAtom.x - coordinatesThisAtom.x;
+   vectorX[1]                       = coordinatesAxisAtom.y - coordinatesThisAtom.y;
+   vectorX[2]                       = coordinatesAxisAtom.z - coordinatesThisAtom.z;
+
+   int axisType                     = multiPoleAtoms[atomIndex].w; 
+     
+   float sum                        = normVector3( vectorZ );
+   
+   // branch based on axis type
+    
+   if( axisType == 1 ){
+
+     // bisector
+     // dx = dx1 + dx2 (in Tinker code)
+     
+     sum               = normVector3( vectorX );
+     
+     vectorZ[0]       += vectorX[0];
+     vectorZ[1]       += vectorX[1];
+     vectorZ[2]       += vectorX[2];
+
+     sum               = normVector3( vectorZ );
+
+   }
+   
+   float dot         = vectorZ[0]*vectorX[0] + vectorZ[1]*vectorX[1] + vectorZ[2]*vectorX[2];
+   
+   vectorX[0]       -= dot*vectorZ[0];
+   vectorX[1]       -= dot*vectorZ[1];
+   vectorX[2]       -= dot*vectorZ[2];
+
+   sum               = normVector3( vectorX );
+   
+   vectorY[0]        = (vectorZ[1]*vectorX[2]) - (vectorZ[2]*vectorX[1]);
+   vectorY[1]        = (vectorZ[2]*vectorX[0]) - (vectorZ[0]*vectorX[2]);
+   vectorY[2]        = (vectorZ[0]*vectorX[1]) - (vectorZ[1]*vectorX[0]);
+
+   float* molDipole  = &(molecularDipole[atomIndex*3]);
+   float* labDipole  = &(labFrameDipole[atomIndex*3]);
+   
+   // set out-of-range elements to 0.0f
+
+   labDipole[0]      = atomIndex >= numOfAtoms ? 0.0f : molDipole[0]*vectorX[0] + molDipole[1]*vectorY[0] + molDipole[2]*vectorZ[0];
+   labDipole[1]      = atomIndex >= numOfAtoms ? 0.0f : molDipole[0]*vectorX[1] + molDipole[1]*vectorY[1] + molDipole[2]*vectorZ[1];
+   labDipole[2]      = atomIndex >= numOfAtoms ? 0.0f : molDipole[0]*vectorX[2] + molDipole[1]*vectorY[2] + molDipole[2]*vectorZ[2];
+   
+   // ---------------------------------------------------------------------------------------
+   
+   const float * mPole[3];
+   float* rPole[3];
+   
+   float* molQuadrupole       = &(molecularQuadrupole[atomIndex*9]);
+   float* labQuadrupole       = &(labFrameQuadrupole[atomIndex*9]);
+   
+   for( int ii = 0; ii < 3; ii++ ){
+      mPole[ii]    = molQuadrupole + ii*3;
+      rPole[ii]    = labQuadrupole + ii*3;
+      rPole[ii][0] = rPole[ii][1] = rPole[ii][2] = 0.0f;
+   }
+   
+   int ii = threadIdx.x;
+   if( ii < 3 ){
+      for( int jj = ii; jj < 3; jj++ ){
+
+         rPole[ii][jj] += vectorX[ii]*vectorX[jj]*mPole[0][0];
+      	rPole[ii][jj] += vectorX[ii]*vectorY[jj]*mPole[0][1];
+      	rPole[ii][jj] += vectorX[ii]*vectorZ[jj]*mPole[0][2];
+      	
+      	rPole[ii][jj] += vectorY[ii]*vectorX[jj]*mPole[1][0];
+      	rPole[ii][jj] += vectorY[ii]*vectorY[jj]*mPole[1][1];
+      	rPole[ii][jj] += vectorY[ii]*vectorZ[jj]*mPole[1][2];
+      	
+      	rPole[ii][jj] += vectorZ[ii]*vectorX[jj]*mPole[2][0];
+      	rPole[ii][jj] += vectorZ[ii]*vectorY[jj]*mPole[2][1];
+      	rPole[ii][jj] += vectorZ[ii]*vectorZ[jj]*mPole[2][2];
+      }
+   }
+
+   __syncthreads();
+
+
+   rPole[1][0] = rPole[0][1];
+   rPole[2][0] = rPole[0][2];
+   rPole[2][1] = rPole[1][2];
+
+   // set out-of-range elements to 0.0f
+
+   labQuadrupole[0]   = atomIndex >= numOfAtoms ? 0.0f : labQuadrupole[0];
+   labQuadrupole[1]   = atomIndex >= numOfAtoms ? 0.0f : labQuadrupole[1];
+   labQuadrupole[2]   = atomIndex >= numOfAtoms ? 0.0f : labQuadrupole[2];
+   labQuadrupole[3]   = atomIndex >= numOfAtoms ? 0.0f : labQuadrupole[3];
+   labQuadrupole[4]   = atomIndex >= numOfAtoms ? 0.0f : labQuadrupole[4];
+   labQuadrupole[5]   = atomIndex >= numOfAtoms ? 0.0f : labQuadrupole[5];
+   labQuadrupole[6]   = atomIndex >= numOfAtoms ? 0.0f : labQuadrupole[6];
+   labQuadrupole[7]   = atomIndex >= numOfAtoms ? 0.0f : labQuadrupole[7];
+   labQuadrupole[8]   = atomIndex >= numOfAtoms ? 0.0f : labQuadrupole[8];
+}
+
+void cudaComputeAmoebaLabFrameMoments( amoebaGpuContext amoebaGpu )
+{
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName = "computeCudaAmoebaLabFrameMoments";
+
+   // ---------------------------------------------------------------------------------------
+
+    gpuContext gpu    = amoebaGpu->gpuContext;
+
+    int numBlocks     =  amoebaGpu->paddedNumberOfAtoms;
+    int numThreads    =  20;
+
+//#define AMOEBA_DEBUG  
+#ifdef AMOEBA_DEBUG
+    if( 0 && amoebaGpu->log ){
+        (void) fprintf( amoebaGpu->log, "%s: numBlocks/atoms=%d\n", methodName, numBlocks ); (void) fflush( amoebaGpu->log );
+        amoebaGpu->psMultipoleParticlesIdsAndAxisType->Download();
+        amoebaGpu->psMolecularDipole->Download();
+        gpu->psPosq4->Download();
+        for( int ii = 0; ii < gpu->natoms; ii++ ){
+            int mIndex = 3*ii;
+             (void) fprintf( amoebaGpu->log,"%6d [%6d %6d %6d] x[%16.9e %16.9e %16.9e] dpl[%16.9e %16.9e %16.9e]\nRot[%16.9e %16.9e %16.9e] [%16.9e %16.9e %16.9e] [%16.9e %16.9e %16.9e]\n\n", ii,
+                             amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].x,
+                             amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].y,
+                             amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].w,
+                             gpu->psPosq4->_pSysStream[0][ii].x,
+                             gpu->psPosq4->_pSysStream[0][ii].y,
+                             gpu->psPosq4->_pSysStream[0][ii].z,
+                             amoebaGpu->psMolecularDipole->_pSysStream[0][mIndex],
+                             amoebaGpu->psMolecularDipole->_pSysStream[0][mIndex+1],
+                             amoebaGpu->psMolecularDipole->_pSysStream[0][mIndex+2] );
+        }
+    }
+//    int64 kernelTime = AmoebaTiming::getTimeOfDay();
+    double kernelTime = 0.0;
+#endif
+
+    kCudaComputeLabFrameMoments_kernel<<< numBlocks, numThreads>>> (
+       gpu->natoms,
+       amoebaGpu->psRotationMatrix->_pDevStream[0],
+       gpu->psPosq4->_pDevStream[0],
+       amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pDevStream[0],
+       amoebaGpu->psMolecularDipole->_pDevStream[0],
+       amoebaGpu->psMolecularQuadrupole->_pDevStream[0],
+       amoebaGpu->psLabFrameDipole->_pDevStream[0],
+       amoebaGpu->psLabFrameQuadrupole->_pDevStream[0] 
+       );
+    LAUNCHERROR(methodName);
+
+#ifdef AMOEBA_DEBUG
+    if( amoebaGpu->log ){
+//        kernelTime          = AmoebaTiming::getTimeOfDay() - kernelTime;
+        static int timestep = 0;
+        timestep++;
+        (void) fprintf( amoebaGpu->log, "Finished rotation kernel execution in %lf us\n", kernelTime ); (void) fflush( amoebaGpu->log );
+        (void) fprintf( amoebaGpu->log, "psLabFrameDipole=%p _pSysStream=%p _pSysStream[0]=%p _pDevStream=%p _pDevStream[0]=%p\n",
+                        amoebaGpu->psLabFrameDipole,  amoebaGpu->psLabFrameDipole->_pSysStream, 
+                        amoebaGpu->psLabFrameDipole->_pSysStream[0], amoebaGpu->psLabFrameDipole->_pDevStream, amoebaGpu->psLabFrameDipole->_pDevStream[0] );
+        fflush( amoebaGpu->log );
+
+        amoebaGpu->psRotationMatrix->Download();
+        amoebaGpu->psLabFrameDipole->Download();
+        (void) fprintf( amoebaGpu->log, "psLabFrameDipole completed\n" );  (void) fflush( amoebaGpu->log );
+
+        amoebaGpu->psLabFrameQuadrupole->Download();
+        (void) fprintf( amoebaGpu->log, "psLabFrameQpole completed\n" );  (void) fflush( amoebaGpu->log );
+
+        int maxPrint = 10;
+        for( int ii = 0; ii < amoebaGpu->paddedNumberOfAtoms; ii++ ){
+
+             int dipoleOffset     = 3*ii;
+             int quadrupoleOffset = 9*ii;
+
+             (void) fprintf( amoebaGpu->log,"\n%6d [%6d %6d %6d] ", ii,
+                             amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].x,
+                             amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].y,
+                             amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].w );
+             // coords
+
+             (void) fprintf( amoebaGpu->log,"x[%16.9e %16.9e %16.9e]\n",
+                             gpu->psPosq4->_pSysStream[0][ii].x,
+                             gpu->psPosq4->_pSysStream[0][ii].y,
+                             gpu->psPosq4->_pSysStream[0][ii].z);
+
+             (void) fprintf( amoebaGpu->log,"   R[%16.9e %16.9e %16.9e] [%16.9e %16.9e %16.9e] [%16.9e %16.9e %16.9e]\n",
+                             amoebaGpu->psRotationMatrix->_pSysStream[0][quadrupoleOffset],
+                             amoebaGpu->psRotationMatrix->_pSysStream[0][quadrupoleOffset+1],
+                             amoebaGpu->psRotationMatrix->_pSysStream[0][quadrupoleOffset+2],
+                             amoebaGpu->psRotationMatrix->_pSysStream[0][quadrupoleOffset+3],
+                             amoebaGpu->psRotationMatrix->_pSysStream[0][quadrupoleOffset+4],
+                             amoebaGpu->psRotationMatrix->_pSysStream[0][quadrupoleOffset+5],
+                             amoebaGpu->psRotationMatrix->_pSysStream[0][quadrupoleOffset+6],
+                             amoebaGpu->psRotationMatrix->_pSysStream[0][quadrupoleOffset+7],
+                             amoebaGpu->psRotationMatrix->_pSysStream[0][quadrupoleOffset+8] );
+
+             // dipole
+
+             (void) fprintf( amoebaGpu->log,"   D[%16.9e %16.9e %16.9e]\n",
+                             amoebaGpu->psLabFrameDipole->_pSysStream[0][dipoleOffset],
+                             amoebaGpu->psLabFrameDipole->_pSysStream[0][dipoleOffset+1],
+                             amoebaGpu->psLabFrameDipole->_pSysStream[0][dipoleOffset+2] );
+    
+             // quadrupole
+
+             (void) fprintf( amoebaGpu->log,"   Q[%16.9e %16.9e %16.9e] [%16.9e %16.9e %16.9e] [%16.9e %16.9e %16.9e]\n",
+                             amoebaGpu->psLabFrameQuadrupole->_pSysStream[0][quadrupoleOffset],
+                             amoebaGpu->psLabFrameQuadrupole->_pSysStream[0][quadrupoleOffset+1],
+                             amoebaGpu->psLabFrameQuadrupole->_pSysStream[0][quadrupoleOffset+2],
+                             amoebaGpu->psLabFrameQuadrupole->_pSysStream[0][quadrupoleOffset+3],
+                             amoebaGpu->psLabFrameQuadrupole->_pSysStream[0][quadrupoleOffset+4],
+                             amoebaGpu->psLabFrameQuadrupole->_pSysStream[0][quadrupoleOffset+5],
+                             amoebaGpu->psLabFrameQuadrupole->_pSysStream[0][quadrupoleOffset+6],
+                             amoebaGpu->psLabFrameQuadrupole->_pSysStream[0][quadrupoleOffset+7],
+                             amoebaGpu->psLabFrameQuadrupole->_pSysStream[0][quadrupoleOffset+8] );
+
+            if( ii == maxPrint && (ii < (gpu->natoms - maxPrint)) ){
+                ii = gpu->natoms - maxPrint;
+            }
+        }
+        int nansDetected   = checkForNansAndInfinities( amoebaGpu->paddedNumberOfAtoms*3, amoebaGpu->psLabFrameDipole );
+            nansDetected  += checkForNansAndInfinities( amoebaGpu->paddedNumberOfAtoms*9, amoebaGpu->psLabFrameQuadrupole );
+        if( nansDetected ){
+             (void) fprintf( amoebaGpu->log,"Nans detected in dipole/quadrupoles.\n" );
+             exit(0);
+        }
+        (void) fflush( amoebaGpu->log );
+    }
+#endif
+
+    if( 0 ){
+//        int particles = particles;
+        int particles = amoebaGpu->paddedNumberOfAtoms;
+        std::vector<int> fileId;
+        //fileId.push_back( 0 );
+        VectorOfDoubleVectors outputVector;
+        cudaLoadCudaFloat4Array( particles, 3, gpu->psPosq4,                     outputVector );
+        cudaLoadCudaFloatArray( particles,  9, amoebaGpu->psRotationMatrix,      outputVector );
+        cudaWriteVectorOfDoubleVectorsToFile( "CudaRotationMatrices", fileId, outputVector );
+    }
+    if( 0 ){
+
+        int particles = amoebaGpu->paddedNumberOfAtoms;
+        std::vector<int> fileId;
+        //fileId.push_back( 0 );
+
+        VectorOfDoubleVectors outputVector;
+        cudaLoadCudaFloat4Array( particles, 3, gpu->psPosq4,                     outputVector );
+        cudaLoadCudaFloatArray( particles,  3, amoebaGpu->psLabFrameDipole,      outputVector );
+        cudaLoadCudaFloatArray( particles,  9, amoebaGpu->psLabFrameQuadrupole,  outputVector );
+        cudaWriteVectorOfDoubleVectorsToFile( "CudaRotatedMoments", fileId, outputVector );
+    }
+  
+}
+
+void kCalculateAmoebaMultipoleForces(amoebaGpuContext amoebaGpu, bool hasAmoebaGeneralizedKirkwood ) 
+{
+    std::string methodName = "kCalculateAmoebaMultipoleForces";
+    //printf("%s \n", methodName.c_str() ); fflush( stdout );
+
+    // compute lab frame moments
+
+    cudaComputeAmoebaLabFrameMoments( amoebaGpu );
+
+    // compute fixed E-field and mutual induced field 
+
+    if( hasAmoebaGeneralizedKirkwood ){
+        cudaComputeAmoebaFixedEAndGkFields( amoebaGpu );
+        cudaComputeAmoebaMutualInducedAndGkField( amoebaGpu );
+    } else {
+        cudaComputeAmoebaFixedEField( amoebaGpu );
+        cudaComputeAmoebaMutualInducedField( amoebaGpu );
+    }
+
+    // check if induce dipole calculation converged -- abort if it did not
+
+    if( amoebaGpu->mutualInducedDone ){
+       //cudaComputeAmoebaElectrostatic( amoebaGpuContextGlobal );
+    } else {
+       (void) fprintf( amoebaGpu->log, "%s induced dipole calculation did not converge -- aborting!\n", methodName.c_str() );
+       (void) fflush( amoebaGpu->log );
+       exit(-1);
+    }
+
+    // calculate electrostatic forces
+
+    cudaComputeAmoebaElectrostatic( amoebaGpu );
+
+    // map torques to forces
+
+    cudaComputeAmoebaMapTorquesAndAddTotalForce( amoebaGpu, amoebaGpu->psTorque, amoebaGpu->psForce, amoebaGpu->gpuContext->psForce4 );
+   
+    if( 0 && amoebaGpu->log ){
+        (void) fprintf( amoebaGpu->log, "Done mapping torques -> forces%s\n", methodName.c_str() ); fflush( NULL );
+        (void) fflush( NULL );
+    }
+}
+
+#undef AMOEBA_DEBUG

plugins/amoeba/platforms/cuda/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 "../../../tests/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 "AmoebaHarmonicBondForce.h"
+#include "AmoebaHarmonicAngleForce.h"
+#include "AmoebaHarmonicInPlaneAngleForce.h"
+#include "AmoebaTorsionForce.h"
+#include "AmoebaPiTorsionForce.h"
+#include "AmoebaStretchBendForce.h"
+#include "AmoebaOutOfPlaneBendForce.h"
+#include "AmoebaTorsionTorsionForce.h"
+#include "AmoebaMultipoleForce.h"
+#include "AmoebaGeneralizedKirkwoodForce.h"
+#include "AmoebaVdwForce.h"
+#include "AmoebaWcaDispersionForce.h"
+#include "AmoebaSASAForce.h"
+
+#include <ctime>
+#include <vector>
+#include <algorithm>
+#include <map>
+#include <cfloat>
+#include <cstring>
+#include <cstdlib>
+#include <typeinfo>
+#include <time.h>
+
+// 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_VDW_FORCE                                   = "AmoebaVdw";
+static std::string AMOEBA_WCA_DISPERSION_FORCE                        = "AmoebaWcaDispersion";
+static std::string AMOEBA_SASA_FORCE                                  = "AmoebaSASA";
+
+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                            = "INCLUDE_OBC_CAVITY_TERM";
+
+#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 AmoebaLastIndex                                   13
+
+#define BOLTZMANN  (1.380658e-23)                         /* (J/K) */
+#define AVOGADRO   (6.0221367e23)                         /* ()    */
+#define RGAS       (BOLTZMANN*AVOGADRO)                   /* (J/(mol K))  */
+#define BOLTZ      (RGAS/1.0e+03)                         /* (kJ/(mol K)) */
+
+#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;
+
+// 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 runTestsUsingAmoebaTinkerParameterFile( MapStringString& argumentMap );
+void appendInputArgumentsToArgumentMap( int numberOfArguments, char* argv[], MapStringString& argumentMap );

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

+#
+# Testing
+#
+ENABLE_TESTING()
+# INCLUDE(${CMAKE_SOURCE_DIR}/platforms/cuda/cuda-cmake/FindCuda.cmake)
+INCLUDE_DIRECTORIES(${CUDA_INCLUDE})
+INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/cuda/include)
+INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/cuda/src)
+INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/cuda/src/kernels)
+
+SET(SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET   "AmoebaTinkerParameterFile" )
+SET(AMOEBA_TINKER_PARAMETER_FILE_SOURCE_FILES    "AmoebaTinkerParameterFile.cpp" )
+SET(AMOEBA_TINKER_PARAMETER_FILE_INCLUDE_FILES   "AmoebaTinkerParameterFile.h" )
+ADD_LIBRARY(${SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET} SHARED ${AMOEBA_TINKER_PARAMETER_FILE_SOURCE_FILES} ${AMOEBA_TINKER_PARAMETER_FILE_INCLUDE_FILES} )
+SET_TARGET_PROPERTIES(${SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET} PROPERTIES COMPILE_FLAGS "-DOPENMM_BUILDING_SHARED_LIBRARY -DLEPTON_BUILDING_SHARED_LIBRARY -DOPENMM_VALIDATE_BUILDING_SHARED_LIBRARY")
+
+Set( SHARED_OPENMM__AMOEBA_TARGET OpenMMAmoeba)
+Set( STATIC_OPENMM_TARGET OpenMMAmoeba_static)
+Set( SHARED_CUDA_TARGET OpenMMAmoebaCuda )
+Set( STATIC_CUDA_TARGET OpenMMCuda_static OpenMMAmoebaCuda_static)
+
+IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
+    SET(SHARED_CUDA_TARGET   ${SHARED_CUDA_TARGET}_d)
+    SET(SHARED_OPENMM__AMOEBA_TARGET ${SHARED_OPENMM__AMOEBA_TARGET}_d)
+    #SET(STATIC_CUDA_TARGET ${STATIC_CUDA_TARGET}_d)
+    #Set(STATIC_OPENMM_TARGET ${STATIC_OPENMM_TARGET}_d)
+ENDIF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
+
+#LINK_DIRECTORIES
+
+# 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
+
+    CUDA_ADD_EXECUTABLE(${TEST_ROOT} ${TEST_PROG})
+    TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET} ${SHARED_OPENMM_TARGET} ${SHARED_CUDA_TARGET} ${SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET})
+    ADD_TEST(${TEST_ROOT} ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT})
+
+    # Link with static library
+#     SET(TEST_STATIC ${TEST_ROOT}Static)
+#     CUDA_ADD_EXECUTABLE(${TEST_STATIC} ${TEST_PROG})
+#     SET_TARGET_PROPERTIES(${TEST_STATIC}
+#                 PROPERTIES
+#                 COMPILE_FLAGS "-DOPENMM_USE_STATIC_LIBRARIES"
+#                 )
+#     TARGET_LINK_LIBRARIES(${TEST_STATIC} ${STATIC_TARGET} ${STATIC_OPENMM_TARGET} ${STATIC_CUDA_TARGET})
+#     ADD_TEST(${TEST_STATIC} ${EXECUTABLE_OUTPUT_PATH}/${TEST_STATIC})
+
+ENDFOREACH(TEST_PROG ${TEST_PROGS})
+
+# TestCudaUsingParameterFile customized w/ command-line argument (input file name used in test) 
+
+#ADD_EXECUTABLE(TestAmoebaCudaUsingParameterFile TstAmoebaCudaUsingParameterFile.cpp)
+#TARGET_LINK_LIBRARIES(TestAmoebaCudaUsingParameterFile ${SHARED_TARGET} ${SHARED_OPENMM_TARGET} ${SHARED_CUDA_TARGET})
+#ADD_TEST(TestCudaUsingParameterFile "${EXECUTABLE_OUTPUT_PATH}/TestCudaUsingParameterFile" "-parameterFileName" "${CMAKE_CURRENT_SOURCE_DIR}/lambdaSdObcParameters.txt")
+#ADD_TEST(TestCudaUsingParameterFile "${EXECUTABLE_OUTPUT_PATH}/TestCudaUsingParameterFile" "-parameterFileName" "${CMAKE_CURRENT_SOURCE_DIR}/bptiMdRfNoPbcParameters.txt")
+#
+#SET(TEST_ROOT TestCudaUsingParameterFile)
+#SET(TEST_PROG TstCudaUsingParameterFile.cpp)
+#SET(TEST_STATIC ${TEST_ROOT}Static)
+#SET(INCLUDE_CUDA_STATIC 1)
+#IF(INCLUDE_CUDA_STATIC)
+#   ADD_EXECUTABLE(${TEST_STATIC} ${TEST_PROG})
+#   SET_TARGET_PROPERTIES(${TEST_STATIC}
+#                         PROPERTIES
+#                         COMPILE_FLAGS "-DOPENMM_USE_STATIC_LIBRARIES"
+#                        )
+#   TARGET_LINK_LIBRARIES(${TEST_STATIC} ${STATIC_TARGET} ${STATIC_BROOK_TARGET})
+#   ADD_TEST(${TEST_STATIC} "${EXECUTABLE_OUTPUT_PATH}/TestCudaUsingParameterFileStatic" "-parameterFileName" "${CMAKE_CURRENT_SOURCE_DIR}/lambdaSdObcParameters.txt")
+#   ADD_TEST(${TEST_STATIC} "${EXECUTABLE_OUTPUT_PATH}/TestCudaUsingParameterFileStatic" "-parameterFileName" "${CMAKE_CURRENT_SOURCE_DIR}/bptiMdRfNoPbcParameters.txt")
+#  ADD_TEST(${TEST_STATIC} "${EXECUTABLE_OUTPUT_PATH}/TestCudaUsingParameterFileStatic" "-parameterFileName" "${CMAKE_CURRENT_SOURCE_DIR}/bptiMdRfPbcParameters.txt" " +checkEnergyForceConsistent -checkForces" )
+#ENDIF(INCLUDE_CUDA_STATIC)

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaForceFromParameterFile.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the Cuda implementation of HarmonicBondForce.
+ */
+
+#include "AmoebaTinkerParameterFile.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+
+const double TOL = 1e-5;
+
+int main( int numberOfArguments, char* argv[] ) {
+
+    try {
+        std::cout << "Running test..." << std::endl;
+        std::string openmmPluginDirectory = "/home/friedrim/src/openmm/trunk/OpenMM/bin";
+        Platform::loadPluginsFromDirectory( openmmPluginDirectory );
+
+        if( numberOfArguments > 1 ){
+            MapStringString argumentMap;
+            appendInputArgumentsToArgumentMap( numberOfArguments, argv, argumentMap );
+            argumentMap[INCLUDE_OBC_CAVITY_TERM] = "0";
+            runTestsUsingAmoebaTinkerParameterFile( argumentMap );
+        }        
+    }
+    catch(const std::exception& e) {
+        std::cout << "exception: " << e.what() << std::endl;
+        std::cout << "FAIL - ERROR.  Test failed." << std::endl;
+        return 1;
+    }
+    std::cout << "PASS - Test succeeded." << std::endl;
+    return 0;
+}

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

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the Cuda implementation of CudaAmoebaHarmonicAngleForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "AmoebaTinkerParameterFile.h"
+#include "openmm/Context.h"
+#include "AmoebaOpenMM.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+
+const double TOL = 1e-5;
+#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);
+    }
+
+    if( log ){
+        (void) fprintf( log, "getPrefactorsGivenAngleCosine: cosine=%10.3e angle=%10.3e ideal=%10.3e\n", cosine, angle, idealAngle ); 
+        (void) fflush( log );
+    }
+
+    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();
+
+    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 );
+    }
+
+    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);
+
+    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 );
+    }
+
+    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 );
+    }
+
+    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 );
+    }
+    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();
+
+    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 );
+    }
+
+    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 << "Running test..." << std::endl;
+        std::string openmmPluginDirectory = "/home/friedrim/src/openmm/trunk/OpenMM/bin";
+        Platform::loadPluginsFromDirectory( openmmPluginDirectory );
+        FILE* log = fopen( "AmoebaHarmonicAngleForce.log", "w" );;
+
+        testOneAngle( log );
+        fclose( log );
+
+    }
+    catch(const std::exception& e) {
+        std::cout << "exception: " << e.what() << std::endl;
+        std::cout << "FAIL - ERROR.  Test failed." << std::endl;
+        return 1;
+    }
+    std::cout << "PASS - Test succeeded." << std::endl;
+    return 0;
+}

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

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the Cuda implementation of HarmonicBondForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "AmoebaTinkerParameterFile.h"
+#include "openmm/Context.h"
+#include "AmoebaOpenMM.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+
+const double TOL = 1e-5;
+
+static void computeAmoebaHarmonicBondForce(int bondIndex,  std::vector<Vec3>& positions, AmoebaHarmonicBondForce& amoebaHarmonicBondForce,
+                                           std::vector<Vec3>& forces, double* energy ) {
+
+    int particle1, particle2;
+    double bondLength;
+    double quadraticK;
+    double cubicK;
+    double quarticK;
+    amoebaHarmonicBondForce.getBondParameters(bondIndex, particle1, particle2,  bondLength,  quadraticK,  cubicK,  quarticK );
+
+    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 );
+    }
+
+    if( log ){
+        (void) fprintf( log, "computeAmoebaHarmonicBondForces: 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 );
+    }
+    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();
+
+    if( log ){
+        (void) fprintf( log, "computeAmoebaHarmonicBondForces: 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 );
+    }
+
+    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->addBond(0, 1, bondLength, quadraticK, cubicK,quarticicK);
+
+    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->addBond(0, 1, bondLength, quadraticK, cubicK, quarticicK);
+    amoebaHarmonicBondForce->addBond(1, 2, bondLength, quadraticK, cubicK, quarticicK);
+
+    system.addForce(amoebaHarmonicBondForce);
+    Context context(system, integrator, Platform::getPlatformByName( "Cuda"));
+    std::vector<Vec3> positions(3);
+
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 1);
+
+    context.setPositions(positions);
+    compareWithExpectedForceAndEnergy( context, *amoebaHarmonicBondForce, TOL, "testTwoBond", log );
+}
+
+int main( int numberOfArguments, char* argv[] ) {
+
+    try {
+        std::cout << "Running test..." << std::endl;
+        std::string openmmPluginDirectory = "/home/friedrim/src/openmm/trunk/OpenMM/bin";
+        Platform::loadPluginsFromDirectory( openmmPluginDirectory );
+        FILE* log = stderr;
+
+        testOneBond( log );
+        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 << "PASS - Test succeeded." << std::endl;
+    return 0;
+}

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

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the Cuda implementation of CudaAmoebaHarmonicInPlaneAngleForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "AmoebaTinkerParameterFile.h"
+#include "openmm/Context.h"
+#include "AmoebaOpenMM.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+
+const double TOL = 1e-5;
+#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);
+    }
+
+    if( log ){
+        (void) fprintf( log, "getPrefactorsGivenInPlaneAngleCosine: cosine=%10.3e angle=%10.3e ideal=%10.3e\n", cosine, angle, idealInPlaneAngle ); 
+        (void) fflush( log );
+    }
+
+    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();
+
+    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 );
+    }
+
+    // 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 ){
+        if( log ){
+            (void) fprintf( log, "computeAmoebaHarmonicInPlaneAngleForce:  rAp2 or rCp2 <= 0.0\n" );
+            (void) fflush( log );
+        }
+        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 );
+    }
+
+    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 );
+    }
+    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();
+
+    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 );
+    }
+
+    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 << "Running test..." << std::endl;
+        std::string openmmPluginDirectory = "/home/friedrim/src/openmm/trunk/OpenMM/bin";
+        Platform::loadPluginsFromDirectory( openmmPluginDirectory );
+        FILE* log = fopen( "AmoebaHarmonicInPlaneAngleForce.log", "w" );;
+
+        testOneAngle( NULL );
+        (void) fclose( log );
+    }
+    catch(const std::exception& e) {
+        std::cout << "exception: " << e.what() << std::endl;
+        std::cout << "FAIL - ERROR.  Test failed." << std::endl;
+        return 1;
+    }
+    std::cout << "PASS - Test succeeded." << std::endl;
+    return 0;
+}

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaOutOfPlaneBendForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the Cuda implementation of CudaAmoebaOutOfPlaneBendForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "AmoebaTinkerParameterFile.h"
+#include "openmm/Context.h"
+#include "AmoebaOpenMM.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+
+const double TOL = 1e-5;
+#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 computeAmoebaOutOfPlaneBendForce(int bondIndex,  std::vector<Vec3>& positions, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce,
+                                             std::vector<Vec3>& forces, double* energy, FILE* log ) {
+
+
+    double kAngleCubic     = amoebaOutOfPlaneBendForce.getAmoebaGlobalOutOfPlaneBendCubic();
+    double kAngleQuartic   = amoebaOutOfPlaneBendForce.getAmoebaGlobalOutOfPlaneBendQuartic();
+    double kAnglePentic    = amoebaOutOfPlaneBendForce.getAmoebaGlobalOutOfPlaneBendPentic();
+    double kAngleSextic    = amoebaOutOfPlaneBendForce.getAmoebaGlobalOutOfPlaneBendSextic();
+
+    int particle1, particle2, particle3, particle4;
+    double kAngleQuadratic;
+    amoebaOutOfPlaneBendForce.getOutOfPlaneBendParameters(bondIndex, particle1, particle2, particle3, particle4, kAngleQuadratic );
+
+    if( log ){
+        (void) fprintf( log, "computeAmoebaOutOfPlaneBendForce: bond %d [%d %d %d %d] k=[%10.3e %10.3e %10.3e %10.3e %10.3e]\n", 
+                             bondIndex, particle1, particle2, particle3, particle4, kAngleQuadratic, kAngleCubic, kAngleQuartic, kAnglePentic, kAngleSextic );
+        (void) fflush( log );
+    }
+
+    enum { A, B, C, D, LastAtomIndex };
+    enum { AB, CB, DB, AD, CD, LastDeltaIndex };
+ 
+    // ---------------------------------------------------------------------------------------
+ 
+    // get deltaR between various combinations of the 4 atoms
+    // and various intermediate terms
+ 
+    double deltaR[LastDeltaIndex][6];
+    for( int ii = 0; ii < 3; ii++ ){
+         deltaR[AB][ii] = positions[particle1][ii] - positions[particle2][ii];
+         deltaR[CB][ii] = positions[particle3][ii] - positions[particle2][ii];
+         deltaR[DB][ii] = positions[particle4][ii] - positions[particle2][ii];
+         deltaR[AD][ii] = positions[particle1][ii] - positions[particle4][ii];
+         deltaR[CD][ii] = positions[particle3][ii] - positions[particle4][ii];
+    }   
+
+    double rDB2  = dotVector3( deltaR[DB], deltaR[DB] );
+    double rAD2  = dotVector3( deltaR[AD], deltaR[AD] );
+    double rCD2  = dotVector3( deltaR[CD], deltaR[CD] );
+ 
+    double tempVector[3];
+    crossProductVector3( deltaR[CB], deltaR[DB], tempVector );
+    double   eE  = dotVector3( deltaR[AB], tempVector  );
+    double  dot  = dotVector3( deltaR[AD],  deltaR[CD] );
+    double   cc  = rAD2*rCD2 - dot*dot;
+ 
+    if( rDB2 <= 0.0 || cc == 0.0 ){
+       return;
+    }
+    double bkk2   = rDB2 - eE*eE/cc;
+    double cosine = sqrt(bkk2/rDB2);
+    double angle;
+    if( cosine >= 1.0 ){
+        angle = 0.0;
+    } else if( cosine <= -1.0 ){
+        angle = PI_M;
+    } else {
+        angle = RADIAN*acos(cosine);
+    }
+ 
+    if( log ){
+        (void) fprintf( log, "computeAmoebaOutOfPlaneBendForce: bkk2=%14.7e rDB2=%14.7e cos=%14.7e dt=%14.7e]\n", 
+                             bkk2, rDB2, cosine, angle );
+        (void) fflush( log );
+    }
+
+    // chain rule
+ 
+    double dt    = angle;
+    double dt2   = dt*dt;
+    double dt3   = dt2*dt;
+    double dt4   = dt2*dt2;
+ 
+    double dEdDt = 2.0 + 3.0*kAngleCubic*dt  + 4.0*kAngleQuartic*dt2 + 5.0*kAnglePentic *dt3 + 6.0*kAngleSextic *dt4;
+     dEdDt      *= kAngleQuadratic*dt*RADIAN;
+ 
+    double dEdCos;
+    dEdCos       = dEdDt/sqrt(cc*bkk2);
+    if( eE > 0.0 ){
+        dEdCos *= -1.0;
+    }
+ 
+    double term = eE/cc;
+ 
+    double dccd[LastAtomIndex][3];
+    for( int ii = 0; ii < 3; ii++ ){
+        dccd[A][ii] = (deltaR[AD][ii]*rCD2 - deltaR[CD][ii]*dot)*term;
+        dccd[C][ii] = (deltaR[CD][ii]*rAD2 - deltaR[AD][ii]*dot)*term;
+        dccd[D][ii] = -1.0*(dccd[A][ii] + dccd[C][ii]);
+    }
+ 
+    double deed[LastAtomIndex][3];
+    crossProductVector3( deltaR[DB], deltaR[CB], deed[A] );
+    crossProductVector3( deltaR[AB], deltaR[DB], deed[C] );
+    crossProductVector3( deltaR[CB], deltaR[AB], deed[D] );
+ 
+    term        = eE/rDB2;
+    deed[D][0] += deltaR[DB][0]*term;
+    deed[D][1] += deltaR[DB][1]*term;
+    deed[D][2] += deltaR[DB][2]*term;
+ 
+    // ---------------------------------------------------------------------------------------
+ 
+    // forces
+ 
+    // calculate forces for atoms a, c, d
+    // the force for b is then -( a+ c + d)
+ 
+    double subForce[LastAtomIndex][3];
+ 
+    for( int jj = 0; jj < LastAtomIndex; jj++ ){
+ 
+        // A, C, D
+  
+        for( int ii = 0; ii < 3; ii++ ){
+            subForce[jj][ii] = dEdCos*( dccd[jj][ii] + deed[jj][ii] );
+        }
+  
+        if( jj == 0 )jj++; // skip B
+  
+        // now compute B
+  
+        if( jj == 3 ){
+           for( int ii = 0; ii < 3; ii++ ){
+               subForce[1][ii] = -1.0*(subForce[0][ii] + subForce[2][ii] + subForce[3][ii]);
+           }
+        }
+    }
+ 
+    // accumulate forces and energy
+
+    forces[particle1][0]       -= subForce[0][0];
+    forces[particle1][1]       -= subForce[0][1];
+    forces[particle1][2]       -= subForce[0][2];
+
+    forces[particle2][0]       -= subForce[1][0];
+    forces[particle2][1]       -= subForce[1][1];
+    forces[particle2][2]       -= subForce[1][2];
+
+    forces[particle3][0]       -= subForce[2][0];
+    forces[particle3][1]       -= subForce[2][1];
+    forces[particle3][2]       -= subForce[2][2];
+    
+    forces[particle4][0]       -= subForce[3][0];
+    forces[particle4][1]       -= subForce[3][1];
+    forces[particle4][2]       -= subForce[3][2];
+    
+    // ---------------------------------------------------------------------------------------
+ 
+    // calculate energy if 'energy' is set
+ 
+    double energyTerm  = 1.0 + kAngleCubic  *dt  +
+                               kAngleQuartic*dt2 +
+                               kAnglePentic *dt3 +
+                               kAngleSextic *dt4;
+    energyTerm        *= kAngleQuadratic*dt2;
+    *energy           += energyTerm;
+    return;
+}
+ 
+static void computeAmoebaOutOfPlaneBendForces( Context& context, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce,
+                                               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 < amoebaOutOfPlaneBendForce.getNumOutOfPlaneBends(); ii++ ){
+        computeAmoebaOutOfPlaneBendForce(ii, positions, amoebaOutOfPlaneBendForce, expectedForces, expectedEnergy, log );
+    }
+
+    if( log ){
+        (void) fprintf( log, "computeAmoebaOutOfPlaneBendForces: 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 );
+    }
+    return;
+
+}
+
+void compareWithExpectedForceAndEnergy( Context& context, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce,
+                                        double tolerance, const std::string& idString, FILE* log) {
+
+    std::vector<Vec3> expectedForces;
+    double expectedEnergy;
+    computeAmoebaOutOfPlaneBendForces( context, amoebaOutOfPlaneBendForce, expectedForces, &expectedEnergy, log );
+   
+    State state                      = context.getState(State::Forces | State::Energy);
+    const std::vector<Vec3> forces   = state.getForces();
+
+    if( log ){
+        (void) fprintf( log, "computeAmoebaOutOfPlaneBendForces: 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 );
+    }
+
+    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 testOneOutOfPlaneBend( 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);
+
+    AmoebaOutOfPlaneBendForce* amoebaOutOfPlaneBendForce = new AmoebaOutOfPlaneBendForce();
+
+    amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendCubic(   -0.1400000E-01 );
+    amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendQuartic(  0.5600000E-04 );
+    amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendPentic(  -0.7000000E-06 );
+    amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendSextic(   0.2200000E-07 );
+
+    double kOutOfPlaneBend = 0.328682196E-01;
+    amoebaOutOfPlaneBendForce->addOutOfPlaneBend(0, 1, 2, 3, kOutOfPlaneBend );
+
+    system.addForce(amoebaOutOfPlaneBendForce);
+    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.269130000E+02,  0.266390000E+02,  0.353100000E+01 );
+
+    positions[2] = Vec3( 0.278860000E+02,  0.264630000E+02,  0.426300000E+01 );
+    positions[3] = Vec3( 0.245568230E+02,  0.250215290E+02,  0.796852800E+01 );
+
+    context.setPositions(positions);
+    compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log );
+
+}
+
+void testOneOutOfPlaneBend2( FILE* log, int setId ) {
+
+    System system;
+    int numberOfParticles = 4;
+    for( int ii = 0; ii < numberOfParticles; ii++ ){
+        system.addParticle(1.0);
+    }
+
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+
+    AmoebaOutOfPlaneBendForce* amoebaOutOfPlaneBendForce = new AmoebaOutOfPlaneBendForce();
+
+    amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendCubic(   -0.1400000E-01 );
+    amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendQuartic(  0.5600000E-04 );
+    amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendPentic(  -0.7000000E-06 );
+    amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendSextic(   0.2200000E-07 );
+/*
+   285    441    442    443    444  0.328682196E-01
+   286    441    442    444    443  0.164493407E-01
+   287    443    442    444    441  0.636650407E-02
+   288    442    444    447    448  0.392956472E-02
+   289    442    444    448    447  0.392956472E-02
+   290    447    444    448    442  0.214755281E-01
+  441  0.893800000E+01  0.439800000E+01  0.343100000E+01
+  442  0.779100000E+01  0.614600000E+01  0.390100000E+01
+  443  0.915400000E+01  0.683900000E+01  0.389400000E+01
+  444  0.101770000E+02  0.619000000E+01  0.379900000E+01
+  445  0.921000000E+01  0.813800000E+01  0.398600000E+01
+  446  0.708500000E+01  0.672900000E+01  0.332700000E+01
+  447  0.744300000E+01  0.605200000E+01  0.491900000E+01
+  448  0.100820000E+02  0.859300000E+01  0.398200000E+01
+  449  0.838000000E+01  0.866100000E+01  0.406000000E+01
+*/
+
+    std::map<int,Vec3> coordinates;
+    coordinates[440] = Vec3(  0.893800000E+01,  0.439800000E+01,  0.343100000E+01 );
+    coordinates[441] = Vec3(  0.779100000E+01,  0.614600000E+01,  0.390100000E+01 );
+    coordinates[442] = Vec3(  0.915400000E+01,  0.683900000E+01,  0.389400000E+01 );
+    coordinates[443] = Vec3(  0.101770000E+02,  0.619000000E+01,  0.379900000E+01 );
+    coordinates[444] = Vec3(  0.921000000E+01,  0.813800000E+01,  0.398600000E+01 );
+    coordinates[445] = Vec3(  0.708500000E+01,  0.672900000E+01,  0.332700000E+01 );
+    coordinates[446] = Vec3(  0.744300000E+01,  0.605200000E+01,  0.491900000E+01 );
+    coordinates[447] = Vec3(  0.100820000E+02,  0.859300000E+01,  0.398200000E+01 );
+    coordinates[448] = Vec3(  0.838000000E+01,  0.866100000E+01,  0.406000000E+01 );
+    
+    double kOutOfPlaneBend = 0.328682196E-01;
+    std::vector<int> particleIndices;
+    if( setId == 1 ){
+        particleIndices.push_back( 441 ); 
+        particleIndices.push_back( 442 ); 
+        particleIndices.push_back( 443 ); 
+        particleIndices.push_back( 444 ); 
+        kOutOfPlaneBend = 0.328682196E-01;
+    } else if( setId == 2 ){
+        particleIndices.push_back( 441 ); 
+        particleIndices.push_back( 442 ); 
+        particleIndices.push_back( 444 ); 
+        particleIndices.push_back( 443 ); 
+        kOutOfPlaneBend = 0.164493407E-01;
+    } else if( setId == 3 ){
+        particleIndices.push_back( 443 ); 
+        particleIndices.push_back( 442 ); 
+        particleIndices.push_back( 444 ); 
+        particleIndices.push_back( 441 ); 
+        kOutOfPlaneBend = 0.636650407E-02;
+    } else if( setId == 4 ){
+        particleIndices.push_back( 442 ); 
+        particleIndices.push_back( 444 ); 
+        particleIndices.push_back( 447 ); 
+        particleIndices.push_back( 448 ); 
+        kOutOfPlaneBend = 0.392956472E-02;
+    } else if( setId == 5 ){
+        particleIndices.push_back( 442 ); 
+        particleIndices.push_back( 444 ); 
+        particleIndices.push_back( 448 ); 
+        particleIndices.push_back( 447 ); 
+        kOutOfPlaneBend = 0.392956472E-02;
+    } else if( setId == 6 ){
+        particleIndices.push_back( 447 ); 
+        particleIndices.push_back( 444 ); 
+        particleIndices.push_back( 448 ); 
+        particleIndices.push_back( 442 ); 
+        kOutOfPlaneBend = 0.214755281E-01;
+    } else {
+        if( log ){
+            (void) fprintf( log, "Set id %d not recognized.\n", setId );
+        }
+        char buffer[1024];
+        (void) sprintf( buffer, "Set id %d not recognized.", setId );
+        throw OpenMMException( buffer );
+    }
+    if( log ){
+        (void) fprintf( log, "Set id %d.\n", setId );
+    }
+
+    amoebaOutOfPlaneBendForce->addOutOfPlaneBend(0, 1, 2, 3, kOutOfPlaneBend );
+
+    system.addForce(amoebaOutOfPlaneBendForce);
+    Context context(system, integrator, Platform::getPlatformByName( "Cuda"));
+    std::vector<Vec3> positions(numberOfParticles);
+
+    for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){
+        if( coordinates.find( particleIndices[ii] ) == coordinates.end() ){
+            if( log ){
+                (void) fprintf( log, "Coordinates %d not loaded.", particleIndices[ii] );
+            }
+            char buffer[1024];
+            (void) sprintf( buffer, "Coordinates %d not loaded.", particleIndices[ii] );
+            throw OpenMMException( buffer );
+        }
+        positions[ii] = coordinates[particleIndices[ii]];
+    }
+
+    context.setPositions(positions);
+    compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log );
+
+    static int iter = 0;
+    static std::map<int,Vec3> totalForces;
+    static double totalEnergy;
+    if( iter == 0 ){
+
+        totalForces[441] = Vec3(  0.0, 0.0, 0.0 );
+        totalForces[442] = Vec3(  0.0, 0.0, 0.0 );
+        totalForces[443] = Vec3(  0.0, 0.0, 0.0 );
+        totalForces[444] = Vec3(  0.0, 0.0, 0.0 );
+        totalForces[445] = Vec3(  0.0, 0.0, 0.0 );
+        totalForces[446] = Vec3(  0.0, 0.0, 0.0 );
+        totalForces[447] = Vec3(  0.0, 0.0, 0.0 );
+        totalForces[448] = Vec3(  0.0, 0.0, 0.0 );
+        totalForces[449] = Vec3(  0.0, 0.0, 0.0 );
+        totalEnergy      = 0.0;
+    }
+    iter++;
+
+    std::vector<Vec3> forces;
+    forces.resize( numberOfParticles );
+    double energy;
+    computeAmoebaOutOfPlaneBendForce( 0, positions, *amoebaOutOfPlaneBendForce, forces, &energy, log );
+
+    totalEnergy += energy;
+    for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){
+        for( unsigned int jj = 0; jj < 3; jj++ ){
+            totalForces[particleIndices[ii]][jj] += forces[ii][jj]; 
+        }
+    }
+
+    if( iter == 6 ){
+    if( log ){
+        (void) fprintf( log, "computeAmoebaOutOfPlaneBendForces: energy=%14.7e\n", totalEnergy);
+        for( std::map<int,Vec3>::iterator ii = totalForces.begin(); ii != totalForces.end(); ii++ ){
+            int particleIndex =  ii->first;
+            Vec3 forces       =  ii->second;
+            (void) fprintf( log, "%6d [%14.7e %14.7e %14.7e] \n", particleIndex,
+                            forces[0], forces[1], forces[2] );
+        }
+        (void) fflush( log );
+    }
+    }
+
+}
+
+int main( int numberOfArguments, char* argv[] ) {
+
+    try {
+        std::cout << "Running test..." << std::endl;
+        std::string openmmPluginDirectory = "/home/friedrim/src/openmm/trunk/OpenMM/bin";
+        Platform::loadPluginsFromDirectory( openmmPluginDirectory );
+        //FILE* log = stderr;
+
+        FILE* log = fopen( "AmoebaOutOfPlaneBendForce1.log", "w" );;
+        //testOneOutOfPlaneBend( log );
+        //testOneOutOfPlaneBend2( log, atoi( argv[1] ) );
+        for( int ii = 1; ii <= 6; ii++ ){
+            testOneOutOfPlaneBend2( log, ii );
+        }
+        (void) fclose( log );
+
+    }
+    catch(const std::exception& e) {
+        std::cout << "exception: " << e.what() << std::endl;
+        std::cout << "FAIL - ERROR.  Test failed." << std::endl;
+        return 1;
+    }
+    std::cout << "PASS - Test succeeded." << std::endl;
+    return 0;
+}

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

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the Cuda implementation of CudaAmoebaPiTorsionForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "AmoebaTinkerParameterFile.h"
+#include "openmm/Context.h"
+#include "AmoebaOpenMM.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+
+const double TOL = 1e-5;
+#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);
+
+    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 );
+    }
+
+    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 );
+    }
+
+    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 );
+    }
+    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();
+
+    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 );
+    }
+
+    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 << "Running test..." << std::endl;
+        std::string openmmPluginDirectory = "/home/friedrim/src/openmm/trunk/OpenMM/bin";
+        Platform::loadPluginsFromDirectory( openmmPluginDirectory );
+        //FILE* log = stderr;
+
+        FILE* log = fopen( "AmoebaPiTorsionForce1.log", "w" );;
+        testOnePiTorsion( log );
+        (void) fclose( log );
+
+    }
+    catch(const std::exception& e) {
+        std::cout << "exception: " << e.what() << std::endl;
+        std::cout << "FAIL - ERROR.  Test failed." << std::endl;
+        return 1;
+    }
+    std::cout << "PASS - Test succeeded." << std::endl;
+    return 0;
+}

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaStretchBendForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the Cuda implementation of CudaAmoebaStretchBendForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "AmoebaTinkerParameterFile.h"
+#include "openmm/Context.h"
+#include "AmoebaOpenMM.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+
+const double TOL = 1e-5;
+#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 computeAmoebaStretchBendForce(int bondIndex,  std::vector<Vec3>& positions, AmoebaStretchBendForce& amoebaStretchBendForce,
+                                          std::vector<Vec3>& forces, double* energy, FILE* log ) {
+
+    int particle1, particle2, particle3;
+    double abBondLength, cbBondLength, angleStretchBend, kStretchBend;
+
+    amoebaStretchBendForce.getStretchBendParameters(bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend);
+    angleStretchBend *= RADIAN;
+    if( log ){
+        (void) fprintf( log, "computeAmoebaStretchBendForce: bond %d [%d %d %d] ab=%10.3e cb=%10.3e angle=%10.3e k=%10.3e\n", 
+                             bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend );
+        (void) fflush( log );
+    }
+
+    enum { A, B, C, LastAtomIndex };
+    enum { AB, CB, CBxAB, ABxP, CBxP, LastDeltaIndex };
+ 
+    // ---------------------------------------------------------------------------------------
+ 
+    // get deltaR between various combinations of the 3 atoms
+    // and various intermediate terms
+ 
+    double deltaR[LastDeltaIndex][3];
+    double rAB2 = 0.0;
+    double rCB2 = 0.0;
+    for( int ii = 0; ii < 3; ii++ ){
+         deltaR[AB][ii]  = positions[particle1][ii] - positions[particle2][ii];
+         rAB2           += deltaR[AB][ii]*deltaR[AB][ii];
+
+         deltaR[CB][ii]  = positions[particle3][ii] - positions[particle2][ii];
+         rCB2           += deltaR[CB][ii]*deltaR[CB][ii];
+    }
+    double rAB   = sqrt( rAB2 );
+    double rCB   = sqrt( rCB2 );
+
+    crossProductVector3( deltaR[CB], deltaR[AB], deltaR[CBxAB] );
+    double  rP   = dotVector3( deltaR[CBxAB], deltaR[CBxAB] );
+            rP   = sqrt( rP );
+ 
+    if( rP <= 0.0 ){
+       return;
+    }
+    double dot    = dotVector3( deltaR[CB], deltaR[AB] );
+    double cosine = dot/(rAB*rCB);
+ 
+    double angle;
+    if( cosine >= 1.0 ){
+       angle = 0.0;
+    } else if( cosine <= -1.0 ){
+       angle = PI_M;
+    } else {
+       angle = RADIAN*acos(cosine);
+    }
+ 
+    double termA = -RADIAN/(rAB2*rP);
+    double termC =  RADIAN/(rCB2*rP);
+ 
+    // P = CBxAB
+ 
+    crossProductVector3( deltaR[AB], deltaR[CBxAB], deltaR[ABxP] );
+    crossProductVector3( deltaR[CB], deltaR[CBxAB], deltaR[CBxP] );
+    for( int ii = 0; ii < 3; ii++ ){
+       deltaR[ABxP][ii] *= termA;
+       deltaR[CBxP][ii] *= termC;
+    }
+ 
+    double dr    = rAB - abBondLength + rCB - cbBondLength;
+ 
+    termA        = 1.0/rAB;
+    termC        = 1.0/rCB;
+ 
+    double term  = kStretchBend;
+ 
+    // ---------------------------------------------------------------------------------------
+ 
+    // forces
+ 
+    // calculate forces for atoms a, b, c
+    // the force for b is then -( a + c)
+ 
+    double subForce[LastAtomIndex][3];
+    double dt = angle - angleStretchBend;
+    for( int jj = 0; jj < 3; jj++ ){
+        subForce[A][jj] = term*(dt*termA*deltaR[AB][jj] + dr*deltaR[ABxP][jj] );
+        subForce[C][jj] = term*(dt*termC*deltaR[CB][jj] + dr*deltaR[CBxP][jj] );
+        subForce[B][jj] = -( subForce[A][jj] + subForce[C][jj] );
+    }
+ 
+    // ---------------------------------------------------------------------------------------
+ 
+    // accumulate forces and energy
+ 
+    forces[particle1][0]       -= subForce[0][0];
+    forces[particle1][1]       -= subForce[0][1];
+    forces[particle1][2]       -= subForce[0][2];
+
+    forces[particle2][0]       -= subForce[1][0];
+    forces[particle2][1]       -= subForce[1][1];
+    forces[particle2][2]       -= subForce[1][2];
+
+    forces[particle3][0]       -= subForce[2][0];
+    forces[particle3][1]       -= subForce[2][1];
+    forces[particle3][2]       -= subForce[2][2];
+
+    *energy                    += term*dt*dr;
+
+    if( log ){
+        (void) fprintf( log, "computeAmoebaStretchBendForce: angle=%10.3e dt=%10.3e dr=%10.3e\n", angle, dt, dr ); 
+        (void) fflush( log );
+    }
+
+    return;
+}
+ 
+static void computeAmoebaStretchBendForces( Context& context, AmoebaStretchBendForce& amoebaStretchBendForce,
+                                          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 < amoebaStretchBendForce.getNumStretchBends(); ii++ ){
+        computeAmoebaStretchBendForce(ii, positions, amoebaStretchBendForce, expectedForces, expectedEnergy, log );
+    }
+
+    if( log ){
+        (void) fprintf( log, "computeAmoebaStretchBendForces: 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 );
+    }
+    return;
+
+}
+
+void compareWithExpectedForceAndEnergy( Context& context, AmoebaStretchBendForce& amoebaStretchBendForce,
+                                        double tolerance, const std::string& idString, FILE* log) {
+
+    std::vector<Vec3> expectedForces;
+    double expectedEnergy;
+    computeAmoebaStretchBendForces( context, amoebaStretchBendForce, expectedForces, &expectedEnergy, log );
+   
+    State state                      = context.getState(State::Forces | State::Energy);
+    const std::vector<Vec3> forces   = state.getForces();
+
+    if( log ){
+        (void) fprintf( log, "computeAmoebaStretchBendForces: 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 );
+    }
+
+    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 testOneStretchBend( 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);
+
+    AmoebaStretchBendForce* amoebaStretchBendForce = new AmoebaStretchBendForce();
+
+    double abLength         = 0.144800000E+01;
+    double cbLength         = 0.101500000E+01;
+    double angleStretchBend = 0.108500000E+03*DegreesToRadians;
+    //double kStretchBend     = 0.750491578E-01;
+    double kStretchBend     = 1.0;
+
+    amoebaStretchBendForce->addStretchBend(0, 1, 2, abLength, cbLength, angleStretchBend, kStretchBend );
+
+    system.addForce(amoebaStretchBendForce);
+    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.273400000E+02,  0.244300000E+02,  0.261400000E+01 );
+    positions[2] = Vec3( 0.269573220E+02,  0.236108860E+02,  0.216376800E+01 );
+
+    context.setPositions(positions);
+    compareWithExpectedForceAndEnergy( context, *amoebaStretchBendForce, TOL, "testOneStretchBend", log );
+
+}
+
+int main( int numberOfArguments, char* argv[] ) {
+
+    try {
+        std::cout << "Running test..." << std::endl;
+        std::string openmmPluginDirectory = "/home/friedrim/src/openmm/trunk/OpenMM/bin";
+        Platform::loadPluginsFromDirectory( openmmPluginDirectory );
+        //FILE* log = stderr;
+
+        FILE* log = fopen( "AmoebaStretchBendForce1.log", "w" );;
+        testOneStretchBend( log );
+        (void) fclose( log );
+
+    }
+    catch(const std::exception& e) {
+        std::cout << "exception: " << e.what() << std::endl;
+        std::cout << "FAIL - ERROR.  Test failed." << std::endl;
+        return 1;
+    }
+    std::cout << "PASS - Test succeeded." << std::endl;
+    return 0;
+}
+
+
+

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaTorsionForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the Cuda implementation of CudaAmoebaTorsionForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "AmoebaTinkerParameterFile.h"
+#include "openmm/Context.h"
+#include "AmoebaOpenMM.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+
+const double TOL = 1e-5;
+#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 computeAmoebaTorsionForce(int bondIndex,  std::vector<Vec3>& positions, AmoebaTorsionForce& amoebaTorsionForce,
+                                      std::vector<Vec3>& forces, double* energy, FILE* log ) {
+
+    int particle1, particle2, particle3, particle4;
+
+    std::vector<double> torsion1;
+    std::vector<double> torsion2;
+    std::vector<double> torsion3;
+
+    torsion1.resize(3);
+    torsion2.resize(3);
+    torsion3.resize(3);
+
+    amoebaTorsionForce.getTorsionParameters(bondIndex, particle1, particle2, particle3, particle4, torsion1, torsion2, torsion3);
+
+    std::vector< std::vector<double> > torsions;
+    torsions.push_back( torsion1 );
+    torsions.push_back( torsion2 );
+    torsions.push_back( torsion3 );
+
+    if( log ){
+        (void) fprintf( log, "computeAmoebaTorsionForce: bond %d [%d %d %d %d]\n", 
+                             bondIndex, particle1, particle2, particle3, particle4 );
+        for( unsigned int ii = 0; ii < 3; ii++ ){
+            (void) fprintf( log, "    [%10.3e %10.3e %10.3e]\n", torsions[ii][0], torsions[ii][1], torsions[ii][2] ); 
+        }
+        (void) fflush( log );
+    }
+
+    enum { BA, CB, DC, CA, DB, LastDeltaIndex };
+    double deltaR[LastDeltaIndex][3];
+    for( int ii = 0; ii < 3; ii++ ){
+        deltaR[BA][ii] = positions[particle2][ii] - positions[particle1][ii];
+        deltaR[CB][ii] = positions[particle3][ii] - positions[particle2][ii];
+        deltaR[DC][ii] = positions[particle4][ii] - positions[particle3][ii];
+        deltaR[CA][ii] = positions[particle3][ii] - positions[particle1][ii];
+        deltaR[DB][ii] = positions[particle4][ii] - positions[particle2][ii];
+    }   
+
+    enum { Xt, Xu, Xtu, LastXtIndex };
+    double crossProducts[LastXtIndex][3];
+    crossProductVector3( deltaR[BA], deltaR[CB], crossProducts[Xt] );
+    crossProductVector3( deltaR[CB], deltaR[DC], crossProducts[Xu] );
+    crossProductVector3( crossProducts[Xt], crossProducts[Xu], crossProducts[Xtu] );
+ 
+    double rT2   = dotVector3( crossProducts[Xt], crossProducts[Xt] );
+    double rU2   = dotVector3( crossProducts[Xu], crossProducts[Xu] );
+    double rTrU  = sqrt( rT2*rU2 );
+    if( rTrU <= 0.0 ){
+        return;
+    }
+ 
+    double rCB   = dotVector3( deltaR[CB], deltaR[CB] );
+         rCB     = sqrt( rCB );
+ 
+    // ---------------------------------------------------------------------------------------
+ 
+    // cos(w), cos(2w), cos(3w), ... 
+    // sin(w), sin(2w), sin(3w), ... 
+  
+    double cosine[6], sine[6];
+   
+    cosine[0]  = dotVector3( crossProducts[Xt], crossProducts[Xu] );
+    cosine[0] /= rTrU;
+ 
+    sine[0]    = dotVector3( deltaR[CB], crossProducts[Xtu] );
+    sine[0]   /= (rCB*rTrU);
+ 
+    for( int ii = 1; ii < 3; ii++ ){
+        cosine[ii] = cosine[0]*cosine[ii-1] - sine[0]*  sine[ii-1];
+          sine[ii] = cosine[0]*  sine[ii-1] + sine[0]*cosine[ii-1];
+    }
+ 
+    // ---------------------------------------------------------------------------------------
+ 
+    // dEdPhi prefactor
+  
+    double dEdPhi = 0.0;
+    for( int ii = 0; ii < 3; ii++ ){
+        dEdPhi += torsions[ii][0]*((double) (ii+1))*( cosine[ii]*sin( torsions[ii][1] ) - sine[ii]*cos( torsions[ii][1] ) );
+    }
+ 
+    // ---------------------------------------------------------------------------------------
+  
+    // dEdtu[0]      = dEdT
+    // dEdtu[1]      = dEdU
+  
+    // tempVector[0] == dEdA: dEdT x CB
+    // tempVector[1] == dEdB: (CA x dEdT) + (dEdU x DC)
+    // tempVector[2] == dEdC: (dEdT x BA) + (DB x dEdU)
+    // tempVector[3] == dEdD: (dEdU x CB)
+   
+    double dEdtu[2][3];
+    double tempVector[6][3];
+ 
+    // dEdT & dEdU
+  
+    crossProductVector3( crossProducts[Xt], deltaR[CB], tempVector[0] );
+    crossProductVector3( crossProducts[Xu], deltaR[CB], tempVector[1] );
+    double norm[2] = { dEdPhi/(rT2*rCB ), -dEdPhi/(rU2*rCB ) };
+    for( int jj = 0; jj < 2; jj++ ){
+        for( int ii = 0; ii < 3; ii++ ){
+            dEdtu[jj][ii] = norm[jj]*tempVector[jj][ii];
+        }
+    }
+ 
+    // dEdA
+   
+    crossProductVector3( dEdtu[0], deltaR[CB], tempVector[0] );
+ 
+    // dEdB
+  
+    crossProductVector3( deltaR[CA], dEdtu[0], tempVector[4] );
+    crossProductVector3( dEdtu[1], deltaR[DC], tempVector[1] );
+ 
+    // dEdC
+ 
+    crossProductVector3( dEdtu[0], deltaR[BA], tempVector[5] );
+    crossProductVector3( deltaR[DB], dEdtu[1], tempVector[2] );
+    for( int jj = 0; jj < 3; jj++ ){
+        tempVector[1][jj] += tempVector[4][jj];
+        tempVector[2][jj] += tempVector[5][jj];
+    }
+ 
+    // dEdD
+  
+    crossProductVector3( dEdtu[1], deltaR[CB], tempVector[3] );
+ 
+    // ---------------------------------------------------------------------------------------
+     
+    // accumulate forces and energy
+  
+    forces[particle1][0]       += tempVector[0][0];
+    forces[particle1][1]       += tempVector[0][1];
+    forces[particle1][2]       += tempVector[0][2];
+
+    forces[particle2][0]       += tempVector[1][0];
+    forces[particle2][1]       += tempVector[1][1];
+    forces[particle2][2]       += tempVector[1][2];
+
+    forces[particle3][0]       += tempVector[2][0];
+    forces[particle3][1]       += tempVector[2][1];
+    forces[particle3][2]       += tempVector[2][2];
+
+    forces[particle4][0]       += tempVector[3][0];
+    forces[particle4][1]       += tempVector[3][1];
+    forces[particle4][2]       += tempVector[3][2];
+
+    double energyTerm = 0.0;
+    for( int ii = 0; ii < 3; ii++ ){
+        energyTerm += torsions[ii][0]*( 1.0 + cosine[ii]*cos( torsions[ii][1] ) + sine[ii]*sin( torsions[ii][1] ) );
+    }
+    *energy    += energyTerm;
+
+} 
+
+static void computeAmoebaTorsionForces( Context& context, AmoebaTorsionForce& amoebaTorsionForce,
+                                        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 < amoebaTorsionForce.getNumTorsions(); ii++ ){
+        computeAmoebaTorsionForce(ii, positions, amoebaTorsionForce, expectedForces, expectedEnergy, log );
+    }
+
+    if( log ){
+        (void) fprintf( log, "computeAmoebaTorsionForces: 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 );
+    }
+    return;
+
+}
+
+void compareWithExpectedForceAndEnergy( Context& context, AmoebaTorsionForce& amoebaTorsionForce,
+                                        double tolerance, const std::string& idString, FILE* log) {
+
+    std::vector<Vec3> expectedForces;
+    double expectedEnergy;
+    computeAmoebaTorsionForces( context, amoebaTorsionForce, expectedForces, &expectedEnergy, log );
+   
+    State state                      = context.getState(State::Forces | State::Energy);
+    const std::vector<Vec3> forces   = state.getForces();
+
+    if( log ){
+        (void) fprintf( log, "computeAmoebaTorsionForces: 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 );
+    }
+
+    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 testOneTorsion( 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);
+
+    AmoebaTorsionForce* amoebaTorsionForce = new AmoebaTorsionForce();
+
+    std::vector<double> torsion1;
+    torsion1.push_back(  0.619500000E+00 );
+    torsion1.push_back(  0.000000000E+00 );
+
+    std::vector<double> torsion2;
+    torsion2.push_back( -0.202500000E+00 );
+    torsion2.push_back(  0.180000000E+03 );
+
+    std::vector<double> torsion3;
+    torsion3.push_back(  0.175000000E-01 );
+    torsion3.push_back(  0.000000000E+00 );
+    amoebaTorsionForce->addTorsion(0, 1, 2, 3, torsion1, torsion2, torsion3 );
+
+    system.addForce(amoebaTorsionForce);
+    Context context(system, integrator, Platform::getPlatformByName( "Cuda"));
+
+    std::vector<Vec3> positions(numberOfParticles);
+
+    positions[0] = Vec3( 0.278860000E+02,  0.264630000E+02,  0.426300000E+01 );
+    positions[1] = Vec3( 0.273400000E+02,  0.244300000E+02,  0.261400000E+01 );
+    positions[2] = Vec3( 0.262660000E+02,  0.254130000E+02,  0.284200000E+01 );
+    positions[3] = Vec3( 0.269130000E+02,  0.266390000E+02,  0.353100000E+01 );
+
+    context.setPositions(positions);
+    compareWithExpectedForceAndEnergy( context, *amoebaTorsionForce, TOL, "testOneTorsion", log );
+
+}
+
+int main( int numberOfArguments, char* argv[] ) {
+
+    try {
+        std::cout << "Running test..." << std::endl;
+        std::string openmmPluginDirectory = "/home/friedrim/src/openmm/trunk/OpenMM/bin";
+        Platform::loadPluginsFromDirectory( openmmPluginDirectory );
+        //FILE* log = stderr;
+        FILE* log = fopen( "AmoebaTorsionForce1.log", "w" );;
+        testOneTorsion( log );
+        (void) fclose( log );
+
+    }
+    catch(const std::exception& e) {
+        std::cout << "exception: " << e.what() << std::endl;
+        std::cout << "FAIL - ERROR.  Test failed." << std::endl;
+        return 1;
+    }
+    std::cout << "PASS - Test succeeded." << std::endl;
+    return 0;
+}

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaTorsionTorsionForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the Cuda implementation of CudaAmoebaTorsionTorsionForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "AmoebaTinkerParameterFile.h"
+#include "openmm/Context.h"
+#include "AmoebaOpenMM.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+
+const double TOL = 1e-5;
+#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 computeAmoebaTorsionTorsionForce(int bondIndex,  std::vector<Vec3>& positions, AmoebaTorsionTorsionForce& amoebaTorsionTorsionForce,
+                                        std::vector<Vec3>& forces, double* energy, FILE* log ) {
+/*
+
+    int particle1, particle2, particle3, particle4, particle5, particle6;
+    double kTorsion;
+
+    amoebaTorsionTorsionForce.getTorsionTorsionParameters(bondIndex, particle1, particle2, particle3, particle4,  particle5, particle6, kTorsion);
+
+    if( log ){
+        (void) fprintf( log, "computeAmoebaTorsionTorsionForce: bond %d [%d %d %d %d %d %d] k=%10.3e\n", 
+                             bondIndex, particle1, particle2, particle3, particle4, particle5, particle6, kTorsion );
+        (void) fflush( log );
+    }
+
+    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 computeAmoebaTorsionTorsionForces( Context& context, AmoebaTorsionTorsionForce& amoebaTorsionTorsionForce,
+                                          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 < amoebaTorsionTorsionForce.getNumTorsionTorsions(); ii++ ){
+        computeAmoebaTorsionTorsionForce(ii, positions, amoebaTorsionTorsionForce, expectedForces, expectedEnergy, log );
+    }
+
+    if( log ){
+        (void) fprintf( log, "computeAmoebaTorsionTorsionForces: 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 );
+    }
+    return;
+
+}
+
+void compareWithExpectedForceAndEnergy( Context& context, AmoebaTorsionTorsionForce& amoebaTorsionTorsionForce,
+                                        double tolerance, const std::string& idString, FILE* log) {
+
+    std::vector<Vec3> expectedForces;
+    double expectedEnergy;
+    computeAmoebaTorsionTorsionForces( context, amoebaTorsionTorsionForce, expectedForces, &expectedEnergy, log );
+   
+    State state                      = context.getState(State::Forces | State::Energy);
+    const std::vector<Vec3> forces   = state.getForces();
+
+    if( log ){
+        (void) fprintf( log, "computeAmoebaTorsionTorsionForces: 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 );
+    }
+
+    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 testOneTorsionTorsion( 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);
+
+    AmoebaTorsionTorsionForce* amoebaTorsionTorsionForce = new AmoebaTorsionTorsionForce();
+
+/*
+    double kTorsion = 6.85;
+    amoebaTorsionTorsionForce->addTorsionTorsion(0, 1, 2, 3, 4, 5, kTorsion );
+
+    system.addForce(amoebaTorsionTorsionForce);
+    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, *amoebaTorsionTorsionForce, TOL, "testOneTorsionTorsion", log );
+*/
+
+}
+
+int main( int numberOfArguments, char* argv[] ) {
+
+    try {
+        std::cout << "Running test..." << std::endl;
+        std::string openmmPluginDirectory = "/home/friedrim/src/openmm/trunk/OpenMM/bin";
+        Platform::loadPluginsFromDirectory( openmmPluginDirectory );
+        //FILE* log = stderr;
+
+        //FILE* log = fopen( "AmoebaTorsionTorsionForce1.log", "w" );;
+        //testOneTorsionTorsion( log );
+        //(void) fclose( log );
+    } catch(const std::exception& e) {
+        std::cout << "exception: " << e.what() << std::endl;
+        std::cout << "FAIL - ERROR.  Test failed." << std::endl;
+        return 1;
+    }
+    std::cout << "PASS - Test succeeded." << std::endl;
+    return 0;
+}