From 10229b7a20452f36714747a85a6922c8cdd16a2d Mon Sep 17 00:00:00 2001 From: Mohammad Soroush Date: Mon, 15 Mar 2021 12:02:36 -0400 Subject: [PATCH 1/4] Fix to the issue #330 --- src/GPU/TransformParticlesCUDAKernel.cu | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/GPU/TransformParticlesCUDAKernel.cu b/src/GPU/TransformParticlesCUDAKernel.cu index a2ccf3068..9a9bd2a41 100644 --- a/src/GPU/TransformParticlesCUDAKernel.cu +++ b/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; @@ -483,7 +483,7 @@ void BrownianMotionRotateParticlesGPU( 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)); @@ -669,7 +669,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]); @@ -734,7 +734,7 @@ void BrownianMotionTranslateParticlesGPU( 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)); @@ -916,7 +916,7 @@ __global__ void BrownianMotionTranslateKernel( } __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]); From b13096e983ecb0c7f57829eb8fd4456b63d26718 Mon Sep 17 00:00:00 2001 From: Mohammad Soroush Date: Mon, 15 Mar 2021 13:17:58 -0400 Subject: [PATCH 2/4] Fix to the issue #328. Fix the recip energy term for brownian Motion GPU --- src/moves/MultiParticle.h | 11 +++++++---- src/moves/MultiParticleBrownianMotion.h | 15 +++++++++------ 2 files changed, 16 insertions(+), 10 deletions(-) diff --git a/src/moves/MultiParticle.h b/src/moves/MultiParticle.h index f546a5ead..b9b3a4029 100644 --- a/src/moves/MultiParticle.h +++ b/src/moves/MultiParticle.h @@ -216,13 +216,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. diff --git a/src/moves/MultiParticleBrownianMotion.h b/src/moves/MultiParticleBrownianMotion.h index f766e7063..d7a2fa5e6 100644 --- a/src/moves/MultiParticleBrownianMotion.h +++ b/src/moves/MultiParticleBrownianMotion.h @@ -198,13 +198,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. @@ -295,8 +298,8 @@ inline void MultiParticleBrownian::CalcEn() //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 From 477c5c072befbcdc212d3269dcb6d27f40369d7e Mon Sep 17 00:00:00 2001 From: Mohammad Soroush Date: Mon, 15 Mar 2021 17:26:24 -0400 Subject: [PATCH 3/4] Fix to the issue #331 --- src/CellList.cpp | 9 +++++- src/GPU/TransformParticlesCUDAKernel.cu | 34 +++++++++++++++----- src/GPU/TransformParticlesCUDAKernel.cuh | 13 +++++--- src/System.cpp | 1 - src/moves/MultiParticle.h | 1 + src/moves/MultiParticleBrownianMotion.h | 41 ++++++++++++++++++------ 6 files changed, 75 insertions(+), 24 deletions(-) diff --git a/src/CellList.cpp b/src/CellList.cpp index 19d35d556..85b90d47e 100644 --- a/src/CellList.cpp +++ b/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(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 diff --git a/src/GPU/TransformParticlesCUDAKernel.cu b/src/GPU/TransformParticlesCUDAKernel.cu index 9a9bd2a41..6bfa7ef4c 100644 --- a/src/GPU/TransformParticlesCUDAKernel.cu +++ b/src/GPU/TransformParticlesCUDAKernel.cu @@ -476,7 +476,8 @@ 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(); @@ -530,7 +531,8 @@ void BrownianMotionRotateParticlesGPU( r_max, step, seed, - BETA); + BETA, + kill); } else { BrownianMotionRotateKernel<<< 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); @@ -727,7 +735,8 @@ 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(); @@ -787,7 +796,8 @@ void BrownianMotionTranslateParticlesGPU( t_max, step, seed, - BETA); + BETA, + kill); } else { BrownianMotionTranslateKernel<<< 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,6 +925,12 @@ __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(); diff --git a/src/GPU/TransformParticlesCUDAKernel.cuh b/src/GPU/TransformParticlesCUDAKernel.cuh index ac2e9d67d..177dd79d3 100644 --- a/src/GPU/TransformParticlesCUDAKernel.cuh +++ b/src/GPU/TransformParticlesCUDAKernel.cuh @@ -14,6 +14,7 @@ typedef r123::Philox4x32 RNG; #include #include "VariablesCUDA.cuh" #include "XYZArray.h" +#include "math.h" void CallTranslateParticlesGPU(VariablesCUDA *vars, const std::vector &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 @@ -167,7 +170,8 @@ __global__ void BrownianMotionRotateKernel( double r_max, unsigned int step, unsigned int seed, - double BETA); + double BETA, + int *kill); template @@ -201,6 +205,7 @@ __global__ void BrownianMotionTranslateKernel( double t_max, unsigned int step, unsigned int seed, - double BETA); + double BETA, + int *kill); #endif diff --git a/src/System.cpp b/src/System.cpp index e9b1dbc6c..1df46e5c9 100644 --- a/src/System.cpp +++ b/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 diff --git a/src/moves/MultiParticle.h b/src/moves/MultiParticle.h index b9b3a4029..47e87b61b 100644 --- a/src/moves/MultiParticle.h +++ b/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(); diff --git a/src/moves/MultiParticleBrownianMotion.h b/src/moves/MultiParticleBrownianMotion.h index d7a2fa5e6..6a4fa0873 100644 --- a/src/moves/MultiParticleBrownianMotion.h +++ b/src/moves/MultiParticleBrownianMotion.h @@ -13,6 +13,9 @@ along with this program, also can be found at . #include #include "Random123Wrapper.h" #ifdef GOMC_CUDA +#include +#include +#include "CUDAMemoryManager.cuh" #include "TransformParticlesCUDAKernel.cuh" #include "VariablesCUDA.cuh" #endif @@ -21,6 +24,12 @@ class MultiParticleBrownian : public MoveBase { public: MultiParticleBrownian(System &sys, StaticVals const& statV); + ~MultiParticleBrownian() { +#ifdef GOMC_CUDA + cudaVars = NULL; + cudaFreeHost(kill); +#endif + } virtual uint Prep(const double subDraw, const double movPerc); virtual void CalcEn(); @@ -50,6 +59,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 +103,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 } @@ -258,8 +270,8 @@ inline uint MultiParticleBrownian::Transform() r123wrapper.GetStep(), r123wrapper.GetSeedValue(), bPick, - isOrthogonal); - + isOrthogonal, + kill); } else { double t_max = moveSetRef.GetTMAX(bPick); BrownianMotionTranslateParticlesGPU( @@ -276,9 +288,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 @@ -293,7 +316,8 @@ 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(); @@ -417,11 +441,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); @@ -501,6 +521,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"; From b9f38cf3a0fabc0684e74f5fdf24187ad42ae73d Mon Sep 17 00:00:00 2001 From: Mohammad Soroush Date: Tue, 16 Mar 2021 11:35:05 -0400 Subject: [PATCH 4/4] Set the kill value at every step --- src/moves/MultiParticleBrownianMotion.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/moves/MultiParticleBrownianMotion.h b/src/moves/MultiParticleBrownianMotion.h index 6a4fa0873..0f0b5c4b0 100644 --- a/src/moves/MultiParticleBrownianMotion.h +++ b/src/moves/MultiParticleBrownianMotion.h @@ -28,6 +28,7 @@ class MultiParticleBrownian : public MoveBase #ifdef GOMC_CUDA cudaVars = NULL; cudaFreeHost(kill); + kill = NULL; #endif } @@ -254,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);