Use random water removal in gmx solvate

Given a number of solvent molecules to remove, gmx solvate
previously used a striding selection. This patch addresses a
TODO to use random selection

Change-Id: Icb065d61349a55b816337fa49470050eac194aa9

diff --git a/src/gromacs/gmxpreprocess/solvate.cpp b/src/gromacs/gmxpreprocess/solvate.cpp
index bf44bd395e803a013bc8a71a4191f255e7942e83..ec1b3360bcde5218ef8cc140511610ff367a2425 100644 (file)
--- a/src/gromacs/gmxpreprocess/solvate.cpp
+++ b/src/gromacs/gmxpreprocess/solvate.cpp
@@ -41,6 +41,7 @@
 #include <cstring>
 
 #include <algorithm>
+#include <random>
 #include <vector>
 
 #include "gromacs/commandline/pargs.h"
@@ -612,22 +613,18 @@ static void removeExtraSolventMolecules(t_atoms *atoms, std::vector<RVec> *x,
                                         std::vector<RVec> *v,
                                         int numberToRemove)
 {
-    gmx::AtomsRemover remover(*atoms);
-    // TODO: It might be nicer to remove a random set of residues, but
-    // in practice this should give a roughly uniform spatial distribution.
-    const int stride = atoms->nr / numberToRemove;
-    for (int i = 0; i < numberToRemove; ++i)
+    gmx::AtomsRemover               remover(*atoms);
+    std::random_device              rd;
+    std::mt19937                    randomNumberGenerator(rd());
+    std::uniform_int_distribution<> randomDistribution(0, atoms->nr - 1);
+    while (numberToRemove > 0)
     {
-        int atomIndex = (i+1)*stride - 1;
-        while (remover.isMarked(atomIndex))
+        int atomIndex = randomDistribution(randomNumberGenerator);
+        if (!remover.isMarked(atomIndex))
         {
-            ++atomIndex;
-            if (atomIndex == atoms->nr)
-            {
-                atomIndex = 0;
-            }
+            remover.markResidue(*atoms, atomIndex, true);
+            numberToRemove--;
         }
-        remover.markResidue(*atoms, atomIndex, true);
     }
     remover.removeMarkedElements(x);
     if (!v->empty())