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