Particle locations reported through the API do not get moved into the first periodic box

platforms/cuda/src/CudaStreamImpl.h

@@ -150,12 +150,10 @@ void CudaStreamImpl<T>::loadFromArray(const void* array) {
         for (int j = 0; j < rowOffset; ++j)
             data[i*rowOffset+j] = paddingValues[j];
     stream->Upload();
-
-    // VisualStudio compiler did not like stream == gpu->psPosq4
-    //if( gpu && stream == gpu->psPosq4 ){
-
-    if( gpu && getName() == "particlePositions" ){
+    if (gpu && getName() == "particlePositions") {
         gpu->bRecalculateBornRadii = true;
+        for (int i = 0; i < gpu->posCellOffsets.size(); i++)
+            gpu->posCellOffsets[i] = make_int3(0, 0, 0);
     }
 }
 
@@ -175,6 +173,14 @@ void CudaStreamImpl<T>::saveToArray(void* array) {
         for (int i = 0; i < getSize(); ++i)
             for (int j = 0; j < width; ++j)
                 arrayData[order[i]*width+j] = data[i*rowOffset+j];
+        if (gpu && getName() == "particlePositions") {
+            for (int i = 0; i < getSize(); i++) {
+                int3 offset = gpu->posCellOffsets[i];
+                arrayData[order[i]*width] -= offset.x*gpu->sim.periodicBoxSizeX;
+                arrayData[order[i]*width+1] -= offset.y*gpu->sim.periodicBoxSizeY;
+                arrayData[order[i]*width+2] -= offset.z*gpu->sim.periodicBoxSizeZ;
+            }
+        }
     }
     else {
         int* arrayData = (int*) array;
@@ -209,6 +215,11 @@ void CudaStreamImpl<T>::fillWithValue(void* value) {
         for (int j = 0; j < rowOffset; ++j)
             data[i*rowOffset+j] = paddingValues[j];
     stream->Upload();
+    if (gpu && getName() == "particlePositions") {
+        gpu->bRecalculateBornRadii = true;
+        for (int i = 0; i < gpu->posCellOffsets.size(); i++)
+            gpu->posCellOffsets[i] = make_int3(0, 0, 0);
+    }
 }
 
 template <class T>

platforms/cuda/src/kernels/gpu.cpp

@@ -1333,6 +1333,7 @@ int gpuAllocateInitialBuffers(gpuContext gpu)
     for (int i = 0; i < (int) gpu->sim.paddedNumberOfAtoms; i++)
         (*gpu->psAtomIndex)[i] = i;
     gpu->psAtomIndex->Upload();
+    gpu->posCellOffsets.resize(gpu->natoms, make_int3(0, 0, 0));
     // Determine randoms
     gpu->seed                           = 1;
     gpu->sim.randomFrames               = 20;
@@ -2385,6 +2386,7 @@ void gpuReorderAtoms(gpuContext gpu)
     vector<int> originalIndex(numAtoms);
     vector<float4> newPosq(numAtoms);
     vector<float4> newVelm(numAtoms);
+    vector<int3> newCellOffsets(numAtoms);
     for (int group = 0; group < (int)gpu->moleculeGroups.size(); group++)
     {
         // Find the center of each molecule.
@@ -2415,9 +2417,12 @@ void gpuReorderAtoms(gpuContext gpu)
 
             for (int i = 0; i < numMolecules; i++)
             {
-                float dx = floor(molPos[i].x/gpu->sim.periodicBoxSizeX)*gpu->sim.periodicBoxSizeX;
-                float dy = floor(molPos[i].y/gpu->sim.periodicBoxSizeY)*gpu->sim.periodicBoxSizeY;
-                float dz = floor(molPos[i].z/gpu->sim.periodicBoxSizeZ)*gpu->sim.periodicBoxSizeZ;
+                int xcell = (int) floor(molPos[i].x/gpu->sim.periodicBoxSizeX);
+                int ycell = (int) floor(molPos[i].y/gpu->sim.periodicBoxSizeY);
+                int zcell = (int) floor(molPos[i].z/gpu->sim.periodicBoxSizeZ);
+                float dx = xcell*gpu->sim.periodicBoxSizeX;
+                float dy = ycell*gpu->sim.periodicBoxSizeY;
+                float dz = zcell*gpu->sim.periodicBoxSizeZ;
                 if (dx != 0.0f || dy != 0.0f || dz != 0.0f)
                 {
                     molPos[i].x -= dx;
@@ -2429,6 +2434,9 @@ void gpuReorderAtoms(gpuContext gpu)
                         posq[atom].x -= dx;
                         posq[atom].y -= dy;
                         posq[atom].z -= dz;
+                        gpu->posCellOffsets[atom].x -= xcell;
+                        gpu->posCellOffsets[atom].y -= ycell;
+                        gpu->posCellOffsets[atom].z -= zcell;
                     }
                 }
             }
@@ -2482,19 +2490,20 @@ void gpuReorderAtoms(gpuContext gpu)
                 originalIndex[newIndex] = (*gpu->psAtomIndex)[oldIndex];
                 newPosq[newIndex] = posq[oldIndex];
                 newVelm[newIndex] = velm[oldIndex];
+                newCellOffsets[newIndex] = gpu->posCellOffsets[oldIndex];
             }
         }
     }
 
     // Update the streams.
 
-    for (int i = 0; i < numAtoms; i++)
+    for (int i = 0; i < numAtoms; i++) {
         posq[i] = newPosq[i];
-    gpu->psPosq4->Upload();
-    for (int i = 0; i < numAtoms; i++)
         velm[i] = newVelm[i];
-    gpu->psVelm4->Upload();
-    for (int i = 0; i < numAtoms; i++)
         (*gpu->psAtomIndex)[i] = originalIndex[i];
+        gpu->posCellOffsets[i] = newCellOffsets[i];
+    }
+    gpu->psPosq4->Upload();
+    gpu->psVelm4->Upload();
     gpu->psAtomIndex->Upload();
 }

platforms/cuda/src/kernels/gputypes.h

@@ -74,6 +74,7 @@ struct _gpuContext {
     unsigned char* pAtomSymbol;
     std::vector<gpuMoleculeGroup> moleculeGroups;
     gpuTabulatedFunction tabulatedFunctions[MAX_TABULATED_FUNCTIONS];
+    std::vector<int3> posCellOffsets;
     float iterations;
     float epsfac;
     float solventDielectric;