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, ePBC;
 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, &ePBC, &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 = gmx_greatest_common_divisor(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, ePBC, 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, ePBC, 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 }