src/gromacs/gmxpreprocess/genion.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "genion.h"
  41
  42 #include <cctype>
  43 #include <cmath>
  44 #include <cstdlib>
  45 #include <cstring>
  46
  47 #include <numeric>
  48 #include <vector>
  49
  50 #include "gromacs/commandline/pargs.h"
  51 #include "gromacs/fileio/confio.h"
  52 #include "gromacs/math/units.h"
  53 #include "gromacs/math/utilities.h"
  54 #include "gromacs/math/vec.h"
  55 #include "gromacs/mdlib/force.h"
  56 #include "gromacs/pbcutil/pbc.h"
  57 #include "gromacs/random/threefry.h"
  58 #include "gromacs/random/uniformintdistribution.h"
  59 #include "gromacs/topology/index.h"
  60 #include "gromacs/topology/topology.h"
  61 #include "gromacs/utility/arrayref.h"
  62 #include "gromacs/utility/arraysize.h"
  63 #include "gromacs/utility/basedefinitions.h"
  64 #include "gromacs/utility/cstringutil.h"
  65 #include "gromacs/utility/fatalerror.h"
  66 #include "gromacs/utility/futil.h"
  67 #include "gromacs/utility/smalloc.h"
  68
  69
  70 /*! \brief Return whether any atoms of two groups are below minimum distance.
  71  *
  72  * \param[in] pbc the periodic boundary conditions
  73  * \param[in] x the coordinates
  74  * \param[in] smallerGroupIndices the atom indices of the first group to check
  75  * \param[in] largerGroupIndices the atom indices of the second group to check
  76  * \param[in] minimumDistance the minimum required distance betwenn any atom in
  77  *                            the first and second group
  78  * \returns true if any distance between an atom from group A and group B is
  79  *               smaller than a minimum distance.
  80  */
  81 static bool groupsCloserThanCutoffWithPbc(t_pbc*                   pbc,
  82                                           rvec                     x[],
  83                                           gmx::ArrayRef<const int> smallerGroupIndices,
  84                                           gmx::ArrayRef<const int> largerGroupIndices,
  85                                           real                     minimumDistance)
  86 {
  87     const real minimumDistance2 = minimumDistance * minimumDistance;
  88     for (int aIndex : largerGroupIndices)
  89     {
  90         for (int bIndex : smallerGroupIndices)
  91         {
  92             rvec dx;
  93             pbc_dx(pbc, x[aIndex], x[bIndex], dx);
  94             if (norm2(dx) < minimumDistance2)
  95             {
  96                 return true;
  97             }
  98         }
  99     }
 100     return false;
 101 }
 102
 103 /*! \brief Calculate the solvent molecule atom indices from molecule number.
 104  *
 105  * \note the solvent group index has to be continuous
 106  *
 107  * \param[in] solventMoleculeNumber the number of the solvent molecule
 108  * \param[in] numberAtomsPerSolventMolecule how many atoms each solvent molecule contains
 109  * \param[in] solventGroupIndex continuous index of solvent atoms
 110  *
 111  * \returns atom indices of the specified solvent molecule
 112  */
 113 static std::vector<int> solventMoleculeIndices(int solventMoleculeNumber,
 114                                                int numberAtomsPerSolventMolecule,
 115                                                gmx::ArrayRef<const int> solventGroupIndex)
 116 {
 117     std::vector<int> indices(numberAtomsPerSolventMolecule);
 118     for (int solventAtomNumber = 0; solventAtomNumber < numberAtomsPerSolventMolecule; ++solventAtomNumber)
 119     {
 120         indices[solventAtomNumber] =
 121                 solventGroupIndex[numberAtomsPerSolventMolecule * solventMoleculeNumber + solventAtomNumber];
 122     }
 123     return indices;
 124 }
 125
 126 static void insert_ion(int                      nsa,
 127                        std::vector<int>*        solventMoleculesForReplacement,
 128                        int                      repl[],
 129                        gmx::ArrayRef<const int> index,
 130                        rvec                     x[],
 131                        t_pbc*                   pbc,
 132                        int                      sign,
 133                        int                      q,
 134                        const char*              ionname,
 135                        t_atoms*                 atoms,
 136                        real                     rmin,
 137                        std::vector<int>*        notSolventGroup)
 138 {
 139     std::vector<int> solventMoleculeAtomsToBeReplaced =
 140             solventMoleculeIndices(solventMoleculesForReplacement->back(), nsa, index);
 141
 142     if (rmin > 0.0)
 143     {
 144         // check for proximity to non-solvent
 145         while (groupsCloserThanCutoffWithPbc(pbc, x, solventMoleculeAtomsToBeReplaced, *notSolventGroup, rmin)
 146                && !solventMoleculesForReplacement->empty())
 147         {
 148             solventMoleculesForReplacement->pop_back();
 149             solventMoleculeAtomsToBeReplaced =
 150                     solventMoleculeIndices(solventMoleculesForReplacement->back(), nsa, index);
 151         }
 152     }
 153
 154     if (solventMoleculesForReplacement->empty())
 155     {
 156         gmx_fatal(FARGS, "No more replaceable solvent!");
 157     }
 158
 159     fprintf(stderr, "Replacing solvent molecule %d (atom %d) with %s\n",
 160             solventMoleculesForReplacement->back(), solventMoleculeAtomsToBeReplaced[0], ionname);
 161
 162     /* Replace solvent molecule charges with ion charge */
 163     notSolventGroup->push_back(solventMoleculeAtomsToBeReplaced[0]);
 164     repl[solventMoleculesForReplacement->back()] = sign;
 165
 166     // The first solvent molecule atom is replaced with an ion and the respective
 167     // charge while the rest of the solvent molecule atoms is set to 0 charge.
 168     atoms->atom[solventMoleculeAtomsToBeReplaced.front()].q = q;
 169     for (auto replacedMoleculeAtom = solventMoleculeAtomsToBeReplaced.begin() + 1;
 170          replacedMoleculeAtom != solventMoleculeAtomsToBeReplaced.end(); ++replacedMoleculeAtom)
 171     {
 172         atoms->atom[*replacedMoleculeAtom].q = 0;
 173     }
 174     solventMoleculesForReplacement->pop_back();
 175 }
 176
 177
 178 static char* aname(const char* mname)
 179 {
 180     char* str;
 181     int   i;
 182
 183     str = gmx_strdup(mname);
 184     i   = std::strlen(str) - 1;
 185     while (i > 1 && (std::isdigit(str[i]) || (str[i] == '+') || (str[i] == '-')))
 186     {
 187         str[i] = '\0';
 188         i--;
 189     }
 190
 191     return str;
 192 }
 193
 194 static void sort_ions(int                      nsa,
 195                       int                      nw,
 196                       const int                repl[],
 197                       gmx::ArrayRef<const int> index,
 198                       t_atoms*                 atoms,
 199                       rvec                     x[],
 200                       char**                   pptr,
 201                       char**                   nptr,
 202                       char**                   paptr,
 203                       char**                   naptr)
 204 {
 205     int   i, j, k, r, np, nn, starta, startr, npi, nni;
 206     rvec* xt;
 207
 208     snew(xt, atoms->nr);
 209
 210     /* Put all the solvent in front and count the added ions */
 211     np = 0;
 212     nn = 0;
 213     j  = index[0];
 214     for (i = 0; i < nw; i++)
 215     {
 216         r = repl[i];
 217         if (r == 0)
 218         {
 219             for (k = 0; k < nsa; k++)
 220             {
 221                 copy_rvec(x[index[nsa * i + k]], xt[j++]);
 222             }
 223         }
 224         else if (r > 0)
 225         {
 226             np++;
 227         }
 228         else if (r < 0)
 229         {
 230             nn++;
 231         }
 232     }
 233
 234     if (np + nn > 0)
 235     {
 236         /* Put the positive and negative ions at the end */
 237         starta = index[nsa * (nw - np - nn)];
 238         startr = atoms->atom[starta].resind;
 239
 240         npi = 0;
 241         nni = 0;
 242         for (i = 0; i < nw; i++)
 243         {
 244             r = repl[i];
 245             if (r > 0)
 246             {
 247                 j = starta + npi;
 248                 k = startr + npi;
 249                 copy_rvec(x[index[nsa * i]], xt[j]);
 250                 atoms->atomname[j]     = paptr;
 251                 atoms->atom[j].resind  = k;
 252                 atoms->resinfo[k].name = pptr;
 253                 npi++;
 254             }
 255             else if (r < 0)
 256             {
 257                 j = starta + np + nni;
 258                 k = startr + np + nni;
 259                 copy_rvec(x[index[nsa * i]], xt[j]);
 260                 atoms->atomname[j]     = naptr;
 261                 atoms->atom[j].resind  = k;
 262                 atoms->resinfo[k].name = nptr;
 263                 nni++;
 264             }
 265         }
 266         for (i = index[nsa * nw - 1] + 1; i < atoms->nr; i++)
 267         {
 268             j                  = i - (nsa - 1) * (np + nn);
 269             atoms->atomname[j] = atoms->atomname[i];
 270             atoms->atom[j]     = atoms->atom[i];
 271             copy_rvec(x[i], xt[j]);
 272         }
 273         atoms->nr -= (nsa - 1) * (np + nn);
 274
 275         /* Copy the new positions back */
 276         for (i = index[0]; i < atoms->nr; i++)
 277         {
 278             copy_rvec(xt[i], x[i]);
 279         }
 280         sfree(xt);
 281     }
 282 }
 283
 284 static void update_topol(const char* topinout, int p_num, int n_num, const char* p_name, const char* n_name, char* grpname)
 285 {
 286     FILE *   fpin, *fpout;
 287     char     buf[STRLEN], buf2[STRLEN], *temp, **mol_line = nullptr;
 288     int      line, i, nmol_line, sol_line, nsol_last;
 289     gmx_bool bMolecules;
 290     char     temporary_filename[STRLEN];
 291
 292     printf("\nProcessing topology\n");
 293     fpin = gmx_ffopen(topinout, "r");
 294     std::strncpy(temporary_filename, "temp.topXXXXXX", STRLEN);
 295     fpout = gmx_fopen_temporary(temporary_filename);
 296
 297     line       = 0;
 298     bMolecules = FALSE;
 299     nmol_line  = 0;
 300     sol_line   = -1;
 301     nsol_last  = -1;
 302     while (fgets(buf, STRLEN, fpin))
 303     {
 304         line++;
 305         std::strcpy(buf2, buf);
 306         if ((temp = std::strchr(buf2, '\n')) != nullptr)
 307         {
 308             temp[0] = '\0';
 309         }
 310         ltrim(buf2);
 311         if (buf2[0] == '[')
 312         {
 313             buf2[0] = ' ';
 314             if ((temp = std::strchr(buf2, '\n')) != nullptr)
 315             {
 316                 temp[0] = '\0';
 317             }
 318             rtrim(buf2);
 319             if (buf2[std::strlen(buf2) - 1] == ']')
 320             {
 321                 buf2[std::strlen(buf2) - 1] = '\0';
 322                 ltrim(buf2);
 323                 rtrim(buf2);
 324                 bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);
 325             }
 326             fprintf(fpout, "%s", buf);
 327         }
 328         else if (!bMolecules)
 329         {
 330             fprintf(fpout, "%s", buf);
 331         }
 332         else
 333         {
 334             /* Check if this is a line with solvent molecules */
 335             sscanf(buf, "%s", buf2);
 336             if (gmx_strcasecmp(buf2, grpname) == 0)
 337             {
 338                 sol_line = nmol_line;
 339                 sscanf(buf, "%*s %d", &nsol_last);
 340             }
 341             /* Store this molecules section line */
 342             srenew(mol_line, nmol_line + 1);
 343             mol_line[nmol_line] = gmx_strdup(buf);
 344             nmol_line++;
 345         }
 346     }
 347     gmx_ffclose(fpin);
 348
 349     if (sol_line == -1)
 350     {
 351         gmx_ffclose(fpout);
 352         gmx_fatal(FARGS,
 353                   "No line with moleculetype '%s' found the [ molecules ] section of file '%s'",
 354                   grpname, topinout);
 355     }
 356     if (nsol_last < p_num + n_num)
 357     {
 358         gmx_ffclose(fpout);
 359         gmx_fatal(FARGS,
 360                   "The last entry for moleculetype '%s' in the [ molecules ] section of file '%s' "
 361                   "has less solvent molecules (%d) than were replaced (%d)",
 362                   grpname, topinout, nsol_last, p_num + n_num);
 363     }
 364
 365     /* Print all the molecule entries */
 366     for (i = 0; i < nmol_line; i++)
 367     {
 368         if (i != sol_line)
 369         {
 370             fprintf(fpout, "%s", mol_line[i]);
 371         }
 372         else
 373         {
 374             printf("Replacing %d solute molecules in topology file (%s) "
 375                    " by %d %s and %d %s ions.\n",
 376                    p_num + n_num, topinout, p_num, p_name, n_num, n_name);
 377             nsol_last -= p_num + n_num;
 378             if (nsol_last > 0)
 379             {
 380                 fprintf(fpout, "%-10s  %d\n", grpname, nsol_last);
 381             }
 382             if (p_num > 0)
 383             {
 384                 fprintf(fpout, "%-15s  %d\n", p_name, p_num);
 385             }
 386             if (n_num > 0)
 387             {
 388                 fprintf(fpout, "%-15s  %d\n", n_name, n_num);
 389             }
 390         }
 391     }
 392     gmx_ffclose(fpout);
 393     make_backup(topinout);
 394     gmx_file_rename(temporary_filename, topinout);
 395 }
 396
 397 /*! \brief Return all atom indices that do not belong to an index group.
 398  * \param[in] nrAtoms the total number of atoms
 399  * \param[in] indexGroup the index group to be inverted. Note that a copy is
 400  *                       made in order to sort the group indices
 401  * \returns the inverted group indices
 402  */
 403 static std::vector<int> invertIndexGroup(int nrAtoms, std::vector<int> indexGroup)
 404 {
 405     // Add the indices -1 and nrAtoms, so that all intervals 0 .. firstofIndexGroup
 406     // as well as lastOfIndexGroup .. nrAtoms are covered for the inverted indexgroup
 407     indexGroup.push_back(-1);
 408     indexGroup.push_back(nrAtoms);
 409     std::sort(indexGroup.begin(), indexGroup.end());
 410
 411     // construct the inverted index group by adding all indicies between two
 412     // indices of indexGroup
 413     std::vector<int> invertedGroup;
 414     for (auto indexGroupIt = std::begin(indexGroup); indexGroupIt != std::end(indexGroup) - 1; ++indexGroupIt)
 415     {
 416         const int firstToAddToInvertedGroup = *indexGroupIt + 1;
 417         const int numIndicesToAdd           = *(indexGroupIt + 1) - firstToAddToInvertedGroup;
 418         if (numIndicesToAdd > 0)
 419         {
 420             invertedGroup.resize(invertedGroup.size() + numIndicesToAdd);
 421             std::iota(std::end(invertedGroup) - numIndicesToAdd, std::end(invertedGroup),
 422                       firstToAddToInvertedGroup);
 423         }
 424     }
 425
 426     return invertedGroup;
 427 }
 428
 429 int gmx_genion(int argc, char* argv[])
 430 {
 431     const char* desc[] = {
 432         "[THISMODULE] randomly replaces solvent molecules with monoatomic ions.",
 433         "The group of solvent molecules should be continuous and all molecules",
 434         "should have the same number of atoms.",
 435         "The user should add the ion molecules to the topology file or use",
 436         "the [TT]-p[tt] option to automatically modify the topology.[PAR]",
 437         "The ion molecule type, residue and atom names in all force fields",
 438         "are the capitalized element names without sign. This molecule name",
 439         "should be given with [TT]-pname[tt] or [TT]-nname[tt], and the",
 440         "[TT][molecules][tt] section of your topology updated accordingly,",
 441         "either by hand or with [TT]-p[tt]. Do not use an atom name instead!",
 442         "[PAR]Ions which can have multiple charge states get the multiplicity",
 443         "added, without sign, for the uncommon states only.[PAR]",
 444         "For larger ions, e.g. sulfate we recommended using [gmx-insert-molecules]."
 445     };
 446     const char* bugs[] = {
 447         "If you specify a salt concentration existing ions are not taken into "
 448         "account. In effect you therefore specify the amount of salt to be added.",
 449     };
 450     int         p_num = 0, n_num = 0, p_q = 1, n_q = -1;
 451     const char *p_name = "NA", *n_name = "CL";
 452     real        rmin = 0.6, conc = 0;
 453     int         seed     = 0;
 454     gmx_bool    bNeutral = FALSE;
 455     t_pargs     pa[]     = {
 456         { "-np", FALSE, etINT, { &p_num }, "Number of positive ions" },
 457         { "-pname", FALSE, etSTR, { &p_name }, "Name of the positive ion" },
 458         { "-pq", FALSE, etINT, { &p_q }, "Charge of the positive ion" },
 459         { "-nn", FALSE, etINT, { &n_num }, "Number of negative ions" },
 460         { "-nname", FALSE, etSTR, { &n_name }, "Name of the negative ion" },
 461         { "-nq", FALSE, etINT, { &n_q }, "Charge of the negative ion" },
 462         { "-rmin", FALSE, etREAL, { &rmin }, "Minimum distance between ions and non-solvent" },
 463         { "-seed", FALSE, etINT, { &seed }, "Seed for random number generator (0 means generate)" },
 464         { "-conc",
 465           FALSE,
 466           etREAL,
 467           { &conc },
 468           "Specify salt concentration (mol/liter). This will add sufficient ions to reach up to "
 469           "the specified concentration as computed from the volume of the cell in the input "
 470           "[REF].tpr[ref] file. Overrides the [TT]-np[tt] and [TT]-nn[tt] options." },
 471         { "-neutral",
 472           FALSE,
 473           etBOOL,
 474           { &bNeutral },
 475           "This option will add enough ions to neutralize the system. These ions are added on top "
 476           "of those specified with [TT]-np[tt]/[TT]-nn[tt] or [TT]-conc[tt]. " }
 477     };
 478     t_topology        top;
 479     rvec*             x;
 480     real              vol;
 481     matrix            box;
 482     t_atoms           atoms;
 483     t_pbc             pbc;
 484     int*              repl;
 485     PbcType           pbcType;
 486     int               nw, nsa, nsalt, iqtot;
 487     gmx_output_env_t* oenv  = nullptr;
 488     t_filenm          fnm[] = { { efTPR, nullptr, nullptr, ffREAD },
 489                        { efNDX, nullptr, nullptr, ffOPTRD },
 490                        { efSTO, "-o", nullptr, ffWRITE },
 491                        { efTOP, "-p", "topol", ffOPTRW } };
 492 #define NFILE asize(fnm)
 493
 494     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa, asize(desc), desc,
 495                            asize(bugs), bugs, &oenv))
 496     {
 497         if (oenv != nullptr)
 498         {
 499             output_env_done(oenv);
 500         }
 501         return 0;
 502     }
 503
 504     /* Check input for something sensible */
 505     if ((p_num < 0) || (n_num < 0))
 506     {
 507         gmx_fatal(FARGS, "Negative number of ions to add?");
 508     }
 509
 510     if (conc > 0 && (p_num > 0 || n_num > 0))
 511     {
 512         fprintf(stderr, "WARNING: -conc specified, overriding -nn and -np.\n");
 513     }
 514
 515     /* Read atom positions and charges */
 516     read_tps_conf(ftp2fn(efTPR, NFILE, fnm), &top, &pbcType, &x, nullptr, box, FALSE);
 517     atoms = top.atoms;
 518
 519     /* Compute total charge */
 520     double qtot = 0;
 521     for (int i = 0; (i < atoms.nr); i++)
 522     {
 523         qtot += atoms.atom[i].q;
 524     }
 525     iqtot = gmx::roundToInt(qtot);
 526
 527
 528     if (conc > 0)
 529     {
 530         /* Compute number of ions to be added */
 531         vol   = det(box);
 532         nsalt = gmx::roundToInt(conc * vol * AVOGADRO / 1e24);
 533         p_num = abs(nsalt * n_q);
 534         n_num = abs(nsalt * p_q);
 535     }
 536     if (bNeutral)
 537     {
 538         int qdelta = p_num * p_q + n_num * n_q + iqtot;
 539
 540         /* Check if the system is neutralizable
 541          * is (qdelta == p_q*p_num + n_q*n_num) solvable for p_num and n_num? */
 542         int gcd = gmx_greatest_common_divisor(n_q, p_q);
 543         if ((qdelta % gcd) != 0)
 544         {
 545             gmx_fatal(FARGS,
 546                       "Can't neutralize this system using -nq %d and"
 547                       " -pq %d.\n",
 548                       n_q, p_q);
 549         }
 550
 551         while (qdelta != 0)
 552         {
 553             while (qdelta < 0)
 554             {
 555                 p_num++;
 556                 qdelta += p_q;
 557             }
 558             while (qdelta > 0)
 559             {
 560                 n_num++;
 561                 qdelta += n_q;
 562             }
 563         }
 564     }
 565
 566     char* pptr  = gmx_strdup(p_name);
 567     char* paptr = aname(p_name);
 568     char* nptr  = gmx_strdup(n_name);
 569     char* naptr = aname(n_name);
 570
 571     if ((p_num == 0) && (n_num == 0))
 572     {
 573         fprintf(stderr, "No ions to add, will just copy input configuration.\n");
 574     }
 575     else
 576     {
 577         char* grpname = nullptr;
 578
 579         printf("Will try to add %d %s ions and %d %s ions.\n", p_num, p_name, n_num, n_name);
 580         printf("Select a continuous group of solvent molecules\n");
 581
 582         std::vector<int> solventGroup;
 583         {
 584             int* index = nullptr;
 585             int  nwa;
 586             get_index(&atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nwa, &index, &grpname);
 587             solventGroup.assign(index, index + nwa);
 588             sfree(index);
 589         }
 590
 591         for (gmx::index i = 1; i < gmx::ssize(solventGroup); i++)
 592         {
 593             if (solventGroup[i] != solventGroup[i - 1] + 1)
 594             {
 595                 gmx_fatal(FARGS,
 596                           "The solvent group %s is not continuous: "
 597                           "index[%d]=%d, index[%d]=%d",
 598                           grpname, int(i), solventGroup[i - 1] + 1, int(i + 1), solventGroup[i] + 1);
 599             }
 600         }
 601         nsa = 1;
 602         while ((nsa < gmx::ssize(solventGroup))
 603                && (atoms.atom[solventGroup[nsa]].resind == atoms.atom[solventGroup[nsa - 1]].resind))
 604         {
 605             nsa++;
 606         }
 607         if (solventGroup.size() % nsa != 0)
 608         {
 609             gmx_fatal(FARGS, "Your solvent group size (%td) is not a multiple of %d",
 610                       gmx::ssize(solventGroup), nsa);
 611         }
 612         nw = solventGroup.size() / nsa;
 613         fprintf(stderr, "Number of (%d-atomic) solvent molecules: %d\n", nsa, nw);
 614         if (p_num + n_num > nw)
 615         {
 616             gmx_fatal(FARGS, "Not enough solvent for adding ions");
 617         }
 618
 619         if (opt2bSet("-p", NFILE, fnm))
 620         {
 621             update_topol(opt2fn("-p", NFILE, fnm), p_num, n_num, p_name, n_name, grpname);
 622         }
 623
 624         snew(repl, nw);
 625         set_pbc(&pbc, pbcType, box);
 626
 627
 628         if (seed == 0)
 629         {
 630             // For now we make do with 32 bits to avoid changing the user input to 64 bit hex
 631             seed = static_cast<int>(gmx::makeRandomSeed());
 632         }
 633         fprintf(stderr, "Using random seed %d.\n", seed);
 634
 635
 636         std::vector<int> notSolventGroup = invertIndexGroup(atoms.nr, solventGroup);
 637
 638         std::vector<int> solventMoleculesForReplacement(nw);
 639         std::iota(std::begin(solventMoleculesForReplacement), std::end(solventMoleculesForReplacement), 0);
 640
 641         // Randomly shuffle the solvent molecules that shall be replaced by ions
 642         // then pick molecules from the back of the list as replacement candidates
 643         gmx::DefaultRandomEngine rng(seed);
 644         std::shuffle(std::begin(solventMoleculesForReplacement),
 645                      std::end(solventMoleculesForReplacement), rng);
 646
 647         /* Now loop over the ions that have to be placed */
 648         while (p_num-- > 0)
 649         {
 650             insert_ion(nsa, &solventMoleculesForReplacement, repl, solventGroup, x, &pbc, 1, p_q,
 651                        p_name, &atoms, rmin, &notSolventGroup);
 652         }
 653         while (n_num-- > 0)
 654         {
 655             insert_ion(nsa, &solventMoleculesForReplacement, repl, solventGroup, x, &pbc, -1, n_q,
 656                        n_name, &atoms, rmin, &notSolventGroup);
 657         }
 658         fprintf(stderr, "\n");
 659
 660         if (nw)
 661         {
 662             sort_ions(nsa, nw, repl, solventGroup, &atoms, x, &pptr, &nptr, &paptr, &naptr);
 663         }
 664
 665         sfree(repl);
 666         sfree(grpname);
 667     }
 668
 669     sfree(atoms.pdbinfo);
 670     atoms.pdbinfo = nullptr;
 671     write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *top.name, &atoms, x, nullptr, pbcType, box);
 672
 673     sfree(pptr);
 674     sfree(paptr);
 675     sfree(nptr);
 676     sfree(naptr);
 677
 678     sfree(x);
 679     output_env_done(oenv);
 680     done_top(&top);
 681
 682     return 0;
 683 }