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,2018,2019, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "genion.h"
  40
  41 #include <cctype>
  42 #include <cmath>
  43 #include <cstdlib>
  44 #include <cstring>
  45
  46 #include <numeric>
  47 #include <vector>
  48
  49 #include "gromacs/commandline/pargs.h"
  50 #include "gromacs/fileio/confio.h"
  51 #include "gromacs/math/units.h"
  52 #include "gromacs/math/utilities.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, ePBC;
 484     int               nw, nsa, nsalt, iqtot;
 485     gmx_output_env_t* oenv  = nullptr;
 486     t_filenm          fnm[] = { { efTPR, nullptr, nullptr, ffREAD },
 487                        { efNDX, nullptr, nullptr, ffOPTRD },
 488                        { efSTO, "-o", nullptr, ffWRITE },
 489                        { efTOP, "-p", "topol", ffOPTRW } };
 490 #define NFILE asize(fnm)
 491
 492     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa, asize(desc), desc,
 493                            asize(bugs), bugs, &oenv))
 494     {
 495         if (oenv != nullptr)
 496         {
 497             output_env_done(oenv);
 498         }
 499         return 0;
 500     }
 501
 502     /* Check input for something sensible */
 503     if ((p_num < 0) || (n_num < 0))
 504     {
 505         gmx_fatal(FARGS, "Negative number of ions to add?");
 506     }
 507
 508     if (conc > 0 && (p_num > 0 || n_num > 0))
 509     {
 510         fprintf(stderr, "WARNING: -conc specified, overriding -nn and -np.\n");
 511     }
 512
 513     /* Read atom positions and charges */
 514     read_tps_conf(ftp2fn(efTPR, NFILE, fnm), &top, &ePBC, &x, nullptr, box, FALSE);
 515     atoms = top.atoms;
 516
 517     /* Compute total charge */
 518     double qtot = 0;
 519     for (int i = 0; (i < atoms.nr); i++)
 520     {
 521         qtot += atoms.atom[i].q;
 522     }
 523     iqtot = gmx::roundToInt(qtot);
 524
 525
 526     if (conc > 0)
 527     {
 528         /* Compute number of ions to be added */
 529         vol   = det(box);
 530         nsalt = gmx::roundToInt(conc * vol * AVOGADRO / 1e24);
 531         p_num = abs(nsalt * n_q);
 532         n_num = abs(nsalt * p_q);
 533     }
 534     if (bNeutral)
 535     {
 536         int qdelta = p_num * p_q + n_num * n_q + iqtot;
 537
 538         /* Check if the system is neutralizable
 539          * is (qdelta == p_q*p_num + n_q*n_num) solvable for p_num and n_num? */
 540         int gcd = gmx_greatest_common_divisor(n_q, p_q);
 541         if ((qdelta % gcd) != 0)
 542         {
 543             gmx_fatal(FARGS,
 544                       "Can't neutralize this system using -nq %d and"
 545                       " -pq %d.\n",
 546                       n_q, p_q);
 547         }
 548
 549         while (qdelta != 0)
 550         {
 551             while (qdelta < 0)
 552             {
 553                 p_num++;
 554                 qdelta += p_q;
 555             }
 556             while (qdelta > 0)
 557             {
 558                 n_num++;
 559                 qdelta += n_q;
 560             }
 561         }
 562     }
 563
 564     char* pptr  = gmx_strdup(p_name);
 565     char* paptr = aname(p_name);
 566     char* nptr  = gmx_strdup(n_name);
 567     char* naptr = aname(n_name);
 568
 569     if ((p_num == 0) && (n_num == 0))
 570     {
 571         fprintf(stderr, "No ions to add, will just copy input configuration.\n");
 572     }
 573     else
 574     {
 575         char* grpname = nullptr;
 576
 577         printf("Will try to add %d %s ions and %d %s ions.\n", p_num, p_name, n_num, n_name);
 578         printf("Select a continuous group of solvent molecules\n");
 579
 580         std::vector<int> solventGroup;
 581         {
 582             int* index = nullptr;
 583             int  nwa;
 584             get_index(&atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nwa, &index, &grpname);
 585             solventGroup.assign(index, index + nwa);
 586             sfree(index);
 587         }
 588
 589         for (gmx::index i = 1; i < gmx::ssize(solventGroup); i++)
 590         {
 591             if (solventGroup[i] != solventGroup[i - 1] + 1)
 592             {
 593                 gmx_fatal(FARGS,
 594                           "The solvent group %s is not continuous: "
 595                           "index[%d]=%d, index[%d]=%d",
 596                           grpname, int(i), solventGroup[i - 1] + 1, int(i + 1), solventGroup[i] + 1);
 597             }
 598         }
 599         nsa = 1;
 600         while ((nsa < gmx::ssize(solventGroup))
 601                && (atoms.atom[solventGroup[nsa]].resind == atoms.atom[solventGroup[nsa - 1]].resind))
 602         {
 603             nsa++;
 604         }
 605         if (solventGroup.size() % nsa != 0)
 606         {
 607             gmx_fatal(FARGS, "Your solvent group size (%td) is not a multiple of %d",
 608                       gmx::ssize(solventGroup), nsa);
 609         }
 610         nw = solventGroup.size() / nsa;
 611         fprintf(stderr, "Number of (%d-atomic) solvent molecules: %d\n", nsa, nw);
 612         if (p_num + n_num > nw)
 613         {
 614             gmx_fatal(FARGS, "Not enough solvent for adding ions");
 615         }
 616
 617         if (opt2bSet("-p", NFILE, fnm))
 618         {
 619             update_topol(opt2fn("-p", NFILE, fnm), p_num, n_num, p_name, n_name, grpname);
 620         }
 621
 622         snew(repl, nw);
 623         set_pbc(&pbc, ePBC, box);
 624
 625
 626         if (seed == 0)
 627         {
 628             // For now we make do with 32 bits to avoid changing the user input to 64 bit hex
 629             seed = static_cast<int>(gmx::makeRandomSeed());
 630         }
 631         fprintf(stderr, "Using random seed %d.\n", seed);
 632
 633
 634         std::vector<int> notSolventGroup = invertIndexGroup(atoms.nr, solventGroup);
 635
 636         std::vector<int> solventMoleculesForReplacement(nw);
 637         std::iota(std::begin(solventMoleculesForReplacement), std::end(solventMoleculesForReplacement), 0);
 638
 639         // Randomly shuffle the solvent molecules that shall be replaced by ions
 640         // then pick molecules from the back of the list as replacement candidates
 641         gmx::DefaultRandomEngine rng(seed);
 642         std::shuffle(std::begin(solventMoleculesForReplacement),
 643                      std::end(solventMoleculesForReplacement), rng);
 644
 645         /* Now loop over the ions that have to be placed */
 646         while (p_num-- > 0)
 647         {
 648             insert_ion(nsa, &solventMoleculesForReplacement, repl, solventGroup, x, &pbc, 1, p_q,
 649                        p_name, &atoms, rmin, &notSolventGroup);
 650         }
 651         while (n_num-- > 0)
 652         {
 653             insert_ion(nsa, &solventMoleculesForReplacement, repl, solventGroup, x, &pbc, -1, n_q,
 654                        n_name, &atoms, rmin, &notSolventGroup);
 655         }
 656         fprintf(stderr, "\n");
 657
 658         if (nw)
 659         {
 660             sort_ions(nsa, nw, repl, solventGroup, &atoms, x, &pptr, &nptr, &paptr, &naptr);
 661         }
 662
 663         sfree(repl);
 664         sfree(grpname);
 665     }
 666
 667     sfree(atoms.pdbinfo);
 668     atoms.pdbinfo = nullptr;
 669     write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *top.name, &atoms, x, nullptr, ePBC, box);
 670
 671     sfree(pptr);
 672     sfree(paptr);
 673     sfree(nptr);
 674     sfree(naptr);
 675
 676     sfree(x);
 677     output_env_done(oenv);
 678     done_top(&top);
 679
 680     return 0;
 681 }