diff --git a/src/ConfigSetup.cpp b/src/ConfigSetup.cpp index 86e1690c7..8eb1b2536 100644 --- a/src/ConfigSetup.cpp +++ b/src/ConfigSetup.cpp @@ -82,6 +82,8 @@ ConfigSetup::ConfigSetup(void) { sys.moves.rotate = DBL_MAX; sys.moves.intraSwap = DBL_MAX; sys.moves.multiParticleEnabled = false; + sys.moves.multiParticleLiquid = true; + sys.moves.multiParticleGas = false; sys.moves.multiParticle = DBL_MAX; sys.moves.multiParticleBrownian = DBL_MAX; sys.moves.regrowth = DBL_MAX; @@ -793,6 +795,14 @@ void ConfigSetup::Init(const char *fileName, MultiSim const *const &multisim) { } printf("%-40s %-4.4f \n", "Info: Multi-Particle move frequency", sys.moves.multiParticle); + } else if (CheckString(line[0], "MultiParticleLiquid")) { + sys.moves.multiParticleLiquid = checkBool(line[1]); + printf("%-40s %-s \n", "Info: Multi-Particle liquid box move", + sys.moves.multiParticleLiquid ? "Active" : "Inactive"); + } else if (CheckString(line[0], "MultiParticleGas")) { + sys.moves.multiParticleGas = checkBool(line[1]); + printf("%-40s %-s \n", "Info: Multi-Particle gas box move", + sys.moves.multiParticleGas ? "Active" : "Inactive"); } else if (CheckString(line[0], "MultiParticleBrownianFreq")) { sys.moves.multiParticleBrownian = stringtod(line[1]); if (sys.moves.multiParticleBrownian > 0.0) { @@ -1961,6 +1971,13 @@ void ConfigSetup::verifyInputs(void) { sys.moves.displace); } } + + if (sys.moves.multiParticleEnabled && !sys.moves.multiParticleLiquid && + !sys.moves.multiParticleGas) { + std::cout << "ERROR: Multi-Particle moves require either liquid or gas " + << "box active!" << std::endl; + exit(EXIT_FAILURE); + } #if ENSEMBLE == NPT if (sys.moves.volume == DBL_MAX) { if (!exptMode) { diff --git a/src/ConfigSetup.h b/src/ConfigSetup.h index c91b1ab56..a39b33085 100644 --- a/src/ConfigSetup.h +++ b/src/ConfigSetup.h @@ -5,8 +5,8 @@ A copy of the MIT License can be found in License.txt along with this program, also can be found at . ********************************************************************************/ -#ifndef CONFIGSETUP_H -#define CONFIGSETUP_H +#ifndef CONFIG_SETUP_H +#define CONFIG_SETUP_H #include //for cerr, cout; #include //for function handle storage. @@ -173,6 +173,7 @@ struct MovePercents { double displace, rotate, intraSwap, intraMemc, regrowth, crankShaft, multiParticle, multiParticleBrownian, intraTargetedSwap; bool multiParticleEnabled; // for both multiparticle and multiparticleBrownian + bool multiParticleLiquid, multiParticleGas; // GEMC: set boxes for MP moves #ifdef VARIABLE_VOLUME double volume; #endif @@ -918,4 +919,4 @@ class ConfigSetup { InputFileReader reader; }; -#endif +#endif /*CONFIG_SETUP_H*/ diff --git a/src/StaticVals.cpp b/src/StaticVals.cpp index 6b576765e..e1cdf8724 100644 --- a/src/StaticVals.cpp +++ b/src/StaticVals.cpp @@ -137,6 +137,8 @@ StaticVals::StaticVals(Setup &set) { multiParticleEnabled = set.config.sys.moves.multiParticleEnabled; + multiParticleLiquid = set.config.sys.moves.multiParticleLiquid; + multiParticleGas = set.config.sys.moves.multiParticleGas; isOrthogonal = true; if (set.config.in.restart.enable) { IsBoxOrthogonal(set.pdb.cryst.cellAngle); diff --git a/src/StaticVals.h b/src/StaticVals.h index f914f51b4..f8ddd191b 100644 --- a/src/StaticVals.h +++ b/src/StaticVals.h @@ -39,6 +39,7 @@ class StaticVals { #endif bool isOrthogonal; bool multiParticleEnabled; + bool multiParticleLiquid, multiParticleGas; Forcefield forcefield; SimEventFrequency simEventFreq; diff --git a/src/moves/MultiParticle.h b/src/moves/MultiParticle.h index 4d17b28f7..f92400c5c 100644 --- a/src/moves/MultiParticle.h +++ b/src/moves/MultiParticle.h @@ -52,6 +52,7 @@ class MultiParticle : public MoveBase { COM newCOMs; int moveType; bool allTranslate; + bool multiParticleLiquid, multiParticleGas; std::vector moleculeIndex; std::vector inForceRange; const MoleculeLookup &molLookup; @@ -63,6 +64,7 @@ class MultiParticle : public MoveBase { const Molecules &mols; double GetCoeff(); + uint ChooseBox(); void CalculateTrialDistRot(); void RotateForceBiased(uint molIndex); void TranslateForceBiased(uint molIndex); @@ -108,7 +110,8 @@ inline MultiParticle::MultiParticle(System &sys, StaticVals const &statV) // If we have only one atom in each kind, it means all molecules // in the system are monoatomic allTranslate = (numAtomsPerKind == molLookup.GetNumKind()); - + multiParticleLiquid = sys.statV.multiParticleLiquid; + multiParticleGas = sys.statV.multiParticleGas; #ifdef GOMC_CUDA cudaVars = sys.statV.forcefield.particles->getCUDAVars(); @@ -167,6 +170,11 @@ inline uint MultiParticle::Prep(const double subDraw, const double movPerc) { uint state = mv::fail_state::NO_FAIL; #if ENSEMBLE == GCMC bPick = mv::BOX0; +#elif ENSEMBLE == GEMC + if (multiParticleLiquid != multiParticleGas) // Only pick one of the two boxes + bPick = ChooseBox(); + else + prng.PickBox(bPick, subDraw, movPerc); #else prng.PickBox(bPick, subDraw, movPerc); #endif @@ -298,6 +306,31 @@ inline uint MultiParticle::PrepNEMTMC(const uint box, const uint midx, return state; } +inline uint MultiParticle::ChooseBox() { + double minDens = 1.0e20; + double maxDens = 0.0; + uint minB = 0, maxB = 0; + for (uint b = 0; b < BOX_TOTAL; ++b) { + double density = 0.0; + for (uint k = 0; k < molLookup.GetNumKind(); ++k) { + density += molLookup.NumKindInBox(k, b) * boxDimRef.volInv[b] * + molRef.kinds[k].molMass; + } + if (density > maxDens) { + maxDens = density; + maxB = b; + } + if (density < minDens) { + minDens = density; + minB = b; + } + } + if (multiParticleLiquid) + return maxB; + else + return minB; +} + inline uint MultiParticle::Transform() { GOMC_EVENT_START(1, GomcProfileEvent::TRANS_MULTIPARTICLE); // Based on the reference force decide whether to displace or rotate each diff --git a/src/moves/MultiParticleBrownianMotion.h b/src/moves/MultiParticleBrownianMotion.h index 3afe6540b..68c095cbb 100644 --- a/src/moves/MultiParticleBrownianMotion.h +++ b/src/moves/MultiParticleBrownianMotion.h @@ -58,12 +58,14 @@ class MultiParticleBrownian : public MoveBase { const MoleculeLookup &molLookup; Random123Wrapper &r123wrapper; bool allTranslate; + bool multiParticleLiquid, multiParticleGas; #ifdef GOMC_CUDA VariablesCUDA *cudaVars; bool isOrthogonal; #endif double GetCoeff(); + uint ChooseBox(); void CalculateTrialDistRot(); void RotateForceBiased(uint molIndex); void TranslateForceBiased(uint molIndex); @@ -100,7 +102,8 @@ inline MultiParticleBrownian::MultiParticleBrownian(System &sys, // If we have only one atom in each kind, it means all molecules // in the system are monoatomic allTranslate = (numAtomsPerKind == molLookup.GetNumKind()); - + multiParticleLiquid = sys.statV.multiParticleLiquid; + multiParticleGas = sys.statV.multiParticleGas; #ifdef GOMC_CUDA cudaVars = sys.statV.forcefield.particles->getCUDAVars(); isOrthogonal = statV.isOrthogonal; @@ -153,6 +156,11 @@ inline uint MultiParticleBrownian::Prep(const double subDraw, uint state = mv::fail_state::NO_FAIL; #if ENSEMBLE == GCMC bPick = mv::BOX0; +#elif ENSEMBLE == GEMC + if (multiParticleLiquid != multiParticleGas) // Only pick one of the two boxes + bPick = ChooseBox(); + else + prng.PickBox(bPick, subDraw, movPerc); #else prng.PickBox(bPick, subDraw, movPerc); #endif @@ -273,6 +281,31 @@ inline uint MultiParticleBrownian::PrepNEMTMC(const uint box, const uint midx, return state; } +inline uint MultiParticleBrownian::ChooseBox() { + double minDens = 1.0e20; + double maxDens = 0.0; + uint minB = 0, maxB = 0; + for (uint b = 0; b < BOX_TOTAL; ++b) { + double density = 0.0; + for (uint k = 0; k < molLookup.GetNumKind(); ++k) { + density += molLookup.NumKindInBox(k, b) * boxDimRef.volInv[b] * + molRef.kinds[k].molMass; + } + if (density > maxDens) { + maxDens = density; + maxB = b; + } + if (density < minDens) { + minDens = density; + minB = b; + } + } + if (multiParticleLiquid) + return maxB; + else + return minB; +} + inline uint MultiParticleBrownian::Transform() { GOMC_EVENT_START(1, GomcProfileEvent::TRANS_MULTIPARTICLE_BM); // Based on the reference force decided whether to displace or rotate each