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