Merge pull request #334 from msoroush/brownian

Fix issues with MultiParticle Brownian.

src/CellList.cpp

@@ -123,7 +123,14 @@ void CellList::AddMol(const int molIndex, const int box, const XYZArray& pos)
   // so list should point to that
   // if this is the first particle in a particular cell.
   while(p != end) {
-    int cell = PositionToCell(pos[p], box);
+    int cell = PositionToCell(pos[p], box); 
+#ifndef NDEBUG
+    if(cell >= static_cast<int>(head[box].size())) {
+      std::cout << "CellList.cpp:129: box " << box << ", pos out of cell: " << pos[p]
+                << std::endl;
+      std::cout << "AxisDimensions: " << dimensions->GetAxis(box) << std::endl;
+    }
+#endif
     //Make the current head index the index the new head points at.
     list[p] = head[box][cell];
     //Assign the new head as our particle index

src/GPU/TransformParticlesCUDAKernel.cu

@@ -50,9 +50,9 @@ __device__ inline void ApplyRotation(double &x, double &y, double &z,
   double axisy = roty * (1.0 / rotLen);
   double axisz = rotz * (1.0 / rotLen);
   double matrix[3][3], cross[3][3], tensor[3][3];
-  double halfAxx = axisx * 0.5;
-  double halfAxy = axisy * 0.5;
-  double halfAxz = axisz * 0.5;
+  double halfAxx = axx * 0.5;
+  double halfAxy = axy * 0.5;
+  double halfAxz = axz * 0.5;
 
   // build cross
   cross[0][0] = 0.0;
@@ -476,14 +476,15 @@ void BrownianMotionRotateParticlesGPU(
   unsigned int step,
   unsigned int seed,
   const int box,
-  bool isOrthogonal)
+  bool isOrthogonal,
+  int *kill)
 {
   int atomCount = newMolPos.Count();
   int molCount = newCOMs.Count();
   int molCountInBox = moleculeInvolved.size();
   int *gpu_moleculeInvolved;
   // Each block would handle one molecule
-  int threadsPerBlock = 64;
+  int threadsPerBlock = 32;
   int blocksPerGrid = molCountInBox;
 
   CUMALLOC((void **) &gpu_moleculeInvolved, molCountInBox * sizeof(int));
@@ -530,7 +531,8 @@ void BrownianMotionRotateParticlesGPU(
       r_max,
       step,
       seed,
-      BETA);
+      BETA,
+      kill);
   } else {
     BrownianMotionRotateKernel<false><<< blocksPerGrid, threadsPerBlock>>>(
       vars->gpu_startAtomIdx,
@@ -559,7 +561,8 @@ void BrownianMotionRotateParticlesGPU(
       r_max,
       step,
       seed,
-      BETA);
+      BETA,
+      kill);
   }
 
   cudaDeviceSynchronize();
@@ -603,7 +606,8 @@ __global__ void BrownianMotionRotateKernel(
   double r_max,
   unsigned int step,
   unsigned int seed,
-  double BETA)
+  double BETA,
+  int *kill)
 {
   //Each grid take cares of one molecule
   int molIndex = moleculeInvolved[blockIdx.x];
@@ -630,6 +634,10 @@ __global__ void BrownianMotionRotateKernel(
     gpu_r_k_x[molIndex] = rot_x;
     gpu_r_k_y[molIndex] = rot_y;
     gpu_r_k_z[molIndex] = rot_z;
+    //check for bad configuration
+    if(!isfinite(rot_x + rot_y + rot_z)) {
+      atomicAdd(kill, 1);
+    }
     // build rotation matrix
     double cross[3][3], tensor[3][3];
     double rotLen = sqrt(rot_x * rot_x + rot_y * rot_y + rot_z * rot_z);
@@ -669,7 +677,7 @@ __global__ void BrownianMotionRotateKernel(
   }
 
   __syncthreads();
-  // use strid of blockDim.x, which is 64
+  // use strid of blockDim.x, which is 32
   // each thread handles one atom rotation
   for(atomIdx = startIdx + threadIdx.x; atomIdx < endIdx; atomIdx += blockDim.x) {
     double3 coor = make_double3(gpu_x[atomIdx], gpu_y[atomIdx], gpu_z[atomIdx]);
@@ -727,14 +735,15 @@ void BrownianMotionTranslateParticlesGPU(
   unsigned int step,
   unsigned int seed,
   const int box,
-  bool isOrthogonal)
+  bool isOrthogonal,
+  int *kill)
 {
   int atomCount = newMolPos.Count();
   int molCount = newCOMs.Count();
   int molCountInBox = moleculeInvolved.size();
   int *gpu_moleculeInvolved;
   // Each block would handle one molecule
-  int threadsPerBlock = 64;
+  int threadsPerBlock = 32;
   int blocksPerGrid = molCountInBox;
 
   CUMALLOC((void **) &gpu_moleculeInvolved, molCountInBox * sizeof(int));
@@ -787,7 +796,8 @@ void BrownianMotionTranslateParticlesGPU(
       t_max,
       step,
       seed,
-      BETA);
+      BETA,
+      kill);
   } else {
     BrownianMotionTranslateKernel<false><<< blocksPerGrid, threadsPerBlock>>>(
       vars->gpu_startAtomIdx,
@@ -819,7 +829,8 @@ void BrownianMotionTranslateParticlesGPU(
       t_max,
       step,
       seed,
-      BETA);
+      BETA,
+      kill);
   }
 
   cudaDeviceSynchronize();
@@ -870,7 +881,8 @@ __global__ void BrownianMotionTranslateKernel(
   double t_max,
   unsigned int step,
   unsigned int seed,
-  double BETA)
+  double BETA,
+  int *kill)
 {
   //Each grid take cares of one molecule
   int molIndex = moleculeInvolved[blockIdx.x];
@@ -913,10 +925,16 @@ __global__ void BrownianMotionTranslateKernel(
     gpu_comx[molIndex] = com.x;
     gpu_comy[molIndex] = com.y;
     gpu_comz[molIndex] = com.z;
+    //check for bad configuration
+    if(!isfinite(shift.x + shift.y + shift.z)) {
+      atomicAdd(kill, 1);
+    } else if (shift.x > halfAx.x || shift.y > halfAx.y || shift.z > halfAx.z) {
+      atomicAdd(kill, 1);
+    }
   }
 
   __syncthreads();
-  // use strid of blockDim.x, which is 64
+  // use strid of blockDim.x, which is 32
   // each thread handles one atom translation
   for(atomIdx = startIdx + threadIdx.x; atomIdx < endIdx; atomIdx += blockDim.x) {
     double3 coor = make_double3(gpu_x[atomIdx], gpu_y[atomIdx], gpu_z[atomIdx]);

src/GPU/TransformParticlesCUDAKernel.cuh

@@ -14,6 +14,7 @@ typedef r123::Philox4x32 RNG;
 #include <cuda_runtime.h>
 #include "VariablesCUDA.cuh"
 #include "XYZArray.h"
+#include "math.h"
 
 void CallTranslateParticlesGPU(VariablesCUDA *vars,
                                const std::vector<int8_t> &isMoleculeInvolved,
@@ -119,7 +120,8 @@ void BrownianMotionRotateParticlesGPU(
   unsigned int step,
   unsigned int seed,
   const int box,
-  bool isOrthogonal);
+  bool isOrthogonal,
+  int *kill);
 
 
 void BrownianMotionTranslateParticlesGPU(
@@ -136,7 +138,8 @@ void BrownianMotionTranslateParticlesGPU(
   unsigned int step,
   unsigned int seed,
   const int box,
-  bool isOrthogonal);
+  bool isOrthogonal,
+  int *kill);
 
 
 template<bool isOrthogonal>
@@ -167,7 +170,8 @@ __global__ void BrownianMotionRotateKernel(
   double r_max,
   unsigned int step,
   unsigned int seed,
-  double BETA);
+  double BETA,
+  int *kill);
 
 
 template<bool isOrthogonal>
@@ -201,6 +205,7 @@ __global__ void BrownianMotionTranslateKernel(
   double t_max,
   unsigned int step,
   unsigned int seed,
-  double BETA);
+  double BETA,
+  int *kill);
 
 #endif

src/System.cpp

@@ -75,7 +75,6 @@ System::~System()
   delete moves[mv::REGROWTH];
   delete moves[mv::INTRA_MEMC];
   delete moves[mv::CRANKSHAFT];
-  delete moves[mv::MULTIPARTICLE];
 #if ENSEMBLE == GEMC || ENSEMBLE == NPT
   delete moves[mv::VOL_TRANSFER];
 #endif

src/moves/MultiParticle.h

@@ -25,6 +25,7 @@ class MultiParticle : public MoveBase
 {
 public:
   MultiParticle(System &sys, StaticVals const& statV);
+  ~MultiParticle() {}
 
   virtual uint Prep(const double subDraw, const double movPerc);
   virtual void CalcEn();
@@ -216,13 +217,16 @@ inline uint MultiParticle::Prep(const double subDraw, const double movPerc)
 inline void MultiParticle::PrepCFCMC(const uint box)
 {
   bPick = box;
-  moveType = prng.randIntExc(mp::MPTOTALTYPES);
-  SetMolInBox(bPick);
-  uint length = molRef.GetKind(moleculeIndex[0]).NumAtoms();
-  if(length == 1) {
+  // In each step, we perform either:
+  // 1- All displacement move.
+  // 2- All rotation move.
+  if(allTranslate) {
     moveType = mp::MPDISPLACE;
+  } else {
+    moveType = prng.randIntExc(mp::MPTOTALTYPES);
   }
 
+  SetMolInBox(bPick);
   //We don't use forces for non-MP moves, so we need to calculate them for the
   //current system if any other moves, besides other MP moves, have been accepted.
   //Or, if this is the first MP move, which is handled with the same flag.

src/moves/MultiParticleBrownianMotion.h

@@ -13,6 +13,9 @@ along with this program, also can be found at <http://www.gnu.org/licenses/>.
 #include <cmath>
 #include "Random123Wrapper.h"
 #ifdef GOMC_CUDA
+#include <cuda.h>
+#include <cuda_runtime.h>
+#include "CUDAMemoryManager.cuh"
 #include "TransformParticlesCUDAKernel.cuh"
 #include "VariablesCUDA.cuh"
 #endif
@@ -21,6 +24,13 @@ class MultiParticleBrownian : public MoveBase
 {
 public:
   MultiParticleBrownian(System &sys, StaticVals const& statV);
+  ~MultiParticleBrownian() {
+#ifdef GOMC_CUDA
+    cudaVars = NULL;
+    cudaFreeHost(kill);
+    kill = NULL;
+#endif
+  }
 
   virtual uint Prep(const double subDraw, const double movPerc);
   virtual void CalcEn();
@@ -50,6 +60,7 @@ class MultiParticleBrownian : public MoveBase
 #ifdef GOMC_CUDA
   VariablesCUDA *cudaVars;
   bool isOrthogonal;
+  int *kill; // kill the simulation if we started with bad configuration
 #endif
 
   double GetCoeff();
@@ -93,6 +104,8 @@ inline MultiParticleBrownian::MultiParticleBrownian(System &sys, StaticVals cons
 #ifdef GOMC_CUDA
   cudaVars = sys.statV.forcefield.particles->getCUDAVars();
   isOrthogonal = statV.isOrthogonal;
+  cudaMallocHost((void**) &kill, sizeof(int));
+  checkLastErrorCUDA(__FILE__, __LINE__);
 #endif
 }
 
@@ -198,13 +211,16 @@ inline uint MultiParticleBrownian::Prep(const double subDraw, const double movPe
 inline void MultiParticleBrownian::PrepCFCMC(const uint box)
 {
   bPick = box;
-  moveType = prng.randIntExc(mp::MPTOTALTYPES);
-  SetMolInBox(bPick);
-  uint length = molRef.GetKind(moleculeIndex[0]).NumAtoms();
-  if(length == 1) {
+  // In each step, we perform either:
+  // 1- All displacement move.
+  // 2- All rotation move.
+  if(allTranslate) {
     moveType = mp::MPDISPLACE;
+  } else {
+    moveType = prng.randIntExc(mp::MPTOTALTYPES);
   }
 
+  SetMolInBox(bPick);
   //We don't use forces for non-MP moves, so we need to calculate them for the
   //current system if any other moves, besides other MP moves, have been accepted.
   //Or, if this is the first MP move, which is handled with the same flag.
@@ -239,6 +255,7 @@ inline uint MultiParticleBrownian::Transform()
   uint state = mv::fail_state::NO_FAIL;
 
 #ifdef GOMC_CUDA
+  kill[0] = 0;
   // This kernel will calculate translation/rotation amount + shifting/rotating
   if(moveType == mp::MPROTATE) {
     double r_max = moveSetRef.GetRMAX(bPick);
@@ -255,8 +272,8 @@ inline uint MultiParticleBrownian::Transform()
       r123wrapper.GetStep(), 
       r123wrapper.GetSeedValue(),
       bPick,
-      isOrthogonal);
-
+      isOrthogonal,
+      kill);
   } else {
     double t_max = moveSetRef.GetTMAX(bPick);
     BrownianMotionTranslateParticlesGPU(
@@ -273,9 +290,20 @@ inline uint MultiParticleBrownian::Transform()
       r123wrapper.GetStep(), 
       r123wrapper.GetSeedValue(),
       bPick,
-      isOrthogonal);
-
+      isOrthogonal,
+      kill);
   }
+  // kill the simulation if we had bad initial configuration
+  if(kill[0]) {
+    std::cout << "Error: Transformation of " << kill[0] << "in Multiparticle Brownian Motion move failed!\n"
+              << "       Either trial rotate/translate is not finite number or their translation \n" 
+              << "       exceeded half of the box length! \n\n";
+    std::cout << "This might be due to bad initial configuration, where atom of the molecules \n" 
+              << "are too close to each other or overlaps. Please equilibrate your system using \n"
+              << "rigid body translation or rotation MC moves, before using the Multiparticle \n"
+              << "Brownian Motion move. \n\n";
+    exit(EXIT_FAILURE);
+} 
 #else
   // Calculate trial translate and rotate
   // move particles according to force and torque and store them in the new pos
@@ -290,13 +318,14 @@ inline void MultiParticleBrownian::CalcEn()
   GOMC_EVENT_START(1, GomcProfileEvent::CALC_EN_MULTIPARTICLE_BM);
   // Calculate the new force and energy and we will compare that to the
   // reference values in Accept() function
-  cellList.GridAll(boxDimRef, newMolsPos, molLookup);
+  //cellList.GridAll(boxDimRef, newMolsPos, molLookup);
+  cellList.GridBox(boxDimRef, newMolsPos, molLookup, bPick);
 
   //back up cached fourier term
   calcEwald->backupMolCache();
 
-  //setup reciprocate vectors for new positions
-  calcEwald->BoxReciprocalSetup(bPick, newMolsPos);
+  //setup reciprocal vectors for new positions
+  calcEwald->BoxReciprocalSums(bPick, newMolsPos);
 
   sysPotNew = sysPotRef;
   //calculate short range energy and force
@@ -414,11 +443,7 @@ inline XYZ MultiParticleBrownian::CalcRandomTransform(XYZ const &lb, double cons
   num.z = lbmax.z + r123wrapper.GetGaussianNumber(molIndex * 3 + 2, 0.0, stdDev);
 
 
-  if(num.Length() >= boxDimRef.axis.Min(bPick)) {
-    std::cout << "Trial Displacement exceeds half of the box length in Brownian Motion MultiParticle move.\n";
-    std::cout << "Trial transform: " << num;
-    exit(EXIT_FAILURE);
-  } else if (!std::isfinite(num.Length())) {
+  if (!std::isfinite(num.Length())) {
     std::cout << "Error: Trial transform is not a finite number in Brownian Motion Multiparticle move.\n";
     std::cout << "       Trial transform: " << num;
     exit(EXIT_FAILURE);
@@ -498,6 +523,7 @@ inline void MultiParticleBrownian::RotateForceBiased(uint molIndex)
 inline void MultiParticleBrownian::TranslateForceBiased(uint molIndex)
 {
   XYZ shift = t_k.Get(molIndex);
+  // check for PBC error and bad initial configuration
   if(shift > boxDimRef.GetHalfAxis(bPick)) {
     std::cout << "Error: Trial Displacement exceed half of the box length in Multiparticle \n" 
               << "       Brownian Motion move!\n";