Fixed bug in RPMD contractions

plugins/rpmd/platforms/cuda/src/kernels/rpmdContraction.cu

@@ -114,8 +114,8 @@ extern "C" __global__ void contractForces(long long* force, long long* contracte
         
         // Store results.
         
-        force[forceIndex] = (long long) (FORCE_SCALE*freal[indexInBlock].x);
-        force[forceIndex+PADDED_NUM_ATOMS] = (long long) (FORCE_SCALE*freal[indexInBlock].y);
-        force[forceIndex+PADDED_NUM_ATOMS*2] = (long long) (FORCE_SCALE*freal[indexInBlock].z);
+        force[forceIndex] += (long long) (FORCE_SCALE*freal[indexInBlock].x);
+        force[forceIndex+PADDED_NUM_ATOMS] += (long long) (FORCE_SCALE*freal[indexInBlock].y);
+        force[forceIndex+PADDED_NUM_ATOMS*2] += (long long) (FORCE_SCALE*freal[indexInBlock].z);
     }
 }

plugins/rpmd/platforms/opencl/src/kernels/rpmdContraction.cl

@@ -113,6 +113,6 @@ __kernel void contractForces(__global real4* force, __global real4* contracted)
         
         // Store results.
 
-        force[index] = convert_real4(FORCE_SCALE*freal[indexInBlock]);
+        force[index] += convert_real4(FORCE_SCALE*freal[indexInBlock]);
     }
 }

plugins/rpmd/platforms/reference/src/ReferenceRpmdKernels.cpp

@@ -336,7 +336,7 @@ void ReferenceIntegrateRPMDStepKernel::computeForces(ContextImpl& context, const
                     q[k] = t_complex(0, 0);
                 fftpack_exec_1d(fft, FFTPACK_BACKWARD, &q[0], &q[0]);
                 for (int k = 0; k < totalCopies; k++)
-                    forces[k][particle][component] = scale2*q[k].re;
+                    forces[k][particle][component] += scale2*q[k].re;
             }
         }
     }