Allows GEMC simulations to control which boxes run MP moves

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 %-4.4f \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 %-4.4f \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) {

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
 <https://opensource.org/licenses/MIT>.
 ********************************************************************************/
-#ifndef CONFIGSETUP_H
-#define CONFIGSETUP_H
+#ifndef CONFIG_SETUP_H
+#define CONFIG_SETUP_H
 
 #include <iostream> //for cerr, cout;
 #include <map>      //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*/

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);

src/StaticVals.h

@@ -39,6 +39,7 @@ class StaticVals {
 #endif
   bool isOrthogonal;
   bool multiParticleEnabled;
+  bool multiParticleLiquid, multiParticleGas;
 
   Forcefield forcefield;
   SimEventFrequency simEventFreq;

src/moves/MultiParticle.h

@@ -52,6 +52,7 @@ class MultiParticle : public MoveBase {
   COM newCOMs;
   int moveType;
   bool allTranslate;
+  bool multiParticleLiquid, multiParticleGas;
   std::vector<uint> moleculeIndex;
   std::vector<int> 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

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