2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
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"
69 /*! \brief Return whether any atoms of two groups are below minimum distance.
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.
80 static bool groupsCloserThanCutoffWithPbc(t_pbc* pbc,
82 gmx::ArrayRef<const int> smallerGroupIndices,
83 gmx::ArrayRef<const int> largerGroupIndices,
86 const real minimumDistance2 = minimumDistance * minimumDistance;
87 for (int aIndex : largerGroupIndices)
89 for (int bIndex : smallerGroupIndices)
92 pbc_dx(pbc, x[aIndex], x[bIndex], dx);
93 if (norm2(dx) < minimumDistance2)
102 /*! \brief Calculate the solvent molecule atom indices from molecule number.
104 * \note the solvent group index has to be continuous
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
110 * \returns atom indices of the specified solvent molecule
112 static std::vector<int> solventMoleculeIndices(int solventMoleculeNumber,
113 int numberAtomsPerSolventMolecule,
114 gmx::ArrayRef<const int> solventGroupIndex)
116 std::vector<int> indices(numberAtomsPerSolventMolecule);
117 for (int solventAtomNumber = 0; solventAtomNumber < numberAtomsPerSolventMolecule; ++solventAtomNumber)
119 indices[solventAtomNumber] =
120 solventGroupIndex[numberAtomsPerSolventMolecule * solventMoleculeNumber + solventAtomNumber];
125 static void insert_ion(int nsa,
126 std::vector<int>* solventMoleculesForReplacement,
128 gmx::ArrayRef<const int> index,
136 std::vector<int>* notSolventGroup)
138 std::vector<int> solventMoleculeAtomsToBeReplaced =
139 solventMoleculeIndices(solventMoleculesForReplacement->back(), nsa, index);
143 // check for proximity to non-solvent
144 while (groupsCloserThanCutoffWithPbc(pbc, x, solventMoleculeAtomsToBeReplaced, *notSolventGroup, rmin)
145 && !solventMoleculesForReplacement->empty())
147 solventMoleculesForReplacement->pop_back();
148 solventMoleculeAtomsToBeReplaced =
149 solventMoleculeIndices(solventMoleculesForReplacement->back(), nsa, index);
153 if (solventMoleculesForReplacement->empty())
155 gmx_fatal(FARGS, "No more replaceable solvent!");
159 "Replacing solvent molecule %d (atom %d) with %s\n",
160 solventMoleculesForReplacement->back(),
161 solventMoleculeAtomsToBeReplaced[0],
164 /* Replace solvent molecule charges with ion charge */
165 notSolventGroup->push_back(solventMoleculeAtomsToBeReplaced[0]);
166 repl[solventMoleculesForReplacement->back()] = sign;
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)
175 atoms->atom[*replacedMoleculeAtom].q = 0;
177 solventMoleculesForReplacement->pop_back();
181 static char* aname(const char* mname)
186 str = gmx_strdup(mname);
187 i = std::strlen(str) - 1;
188 while (i > 1 && (std::isdigit(str[i]) || (str[i] == '+') || (str[i] == '-')))
197 static void sort_ions(int nsa,
200 gmx::ArrayRef<const int> index,
208 int i, j, k, r, np, nn, starta, startr, npi, nni;
213 /* Put all the solvent in front and count the added ions */
217 for (i = 0; i < nw; i++)
222 for (k = 0; k < nsa; k++)
224 copy_rvec(x[index[nsa * i + k]], xt[j++]);
239 /* Put the positive and negative ions at the end */
240 starta = index[nsa * (nw - np - nn)];
241 startr = atoms->atom[starta].resind;
245 for (i = 0; i < nw; i++)
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;
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;
269 for (i = index[nsa * nw - 1] + 1; i < atoms->nr; i++)
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]);
276 atoms->nr -= (nsa - 1) * (np + nn);
278 /* Copy the new positions back */
279 for (i = index[0]; i < atoms->nr; i++)
281 copy_rvec(xt[i], x[i]);
287 static void update_topol(const char* topinout, int p_num, int n_num, const char* p_name, const char* n_name, char* grpname)
290 char buf[STRLEN], buf2[STRLEN], *temp, **mol_line = nullptr;
291 int line, i, nmol_line, sol_line, nsol_last;
293 char temporary_filename[STRLEN];
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);
305 while (fgets(buf, STRLEN, fpin))
308 std::strcpy(buf2, buf);
309 if ((temp = std::strchr(buf2, '\n')) != nullptr)
317 if ((temp = std::strchr(buf2, '\n')) != nullptr)
322 if (buf2[std::strlen(buf2) - 1] == ']')
324 buf2[std::strlen(buf2) - 1] = '\0';
327 bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);
329 fprintf(fpout, "%s", buf);
331 else if (!bMolecules)
333 fprintf(fpout, "%s", buf);
337 /* Check if this is a line with solvent molecules */
338 sscanf(buf, "%s", buf2);
339 if (gmx_strcasecmp(buf2, grpname) == 0)
341 sol_line = nmol_line;
342 sscanf(buf, "%*s %d", &nsol_last);
344 /* Store this molecules section line */
345 srenew(mol_line, nmol_line + 1);
346 mol_line[nmol_line] = gmx_strdup(buf);
356 "No line with moleculetype '%s' found the [ molecules ] section of file '%s'",
360 if (nsol_last < p_num + n_num)
364 "The last entry for moleculetype '%s' in the [ molecules ] section of file '%s' "
365 "has less solvent molecules (%d) than were replaced (%d)",
372 /* Print all the molecule entries */
373 for (i = 0; i < nmol_line; i++)
377 fprintf(fpout, "%s", mol_line[i]);
381 printf("Replacing %d solute molecules in topology file (%s) "
382 " by %d %s and %d %s ions.\n",
389 nsol_last -= p_num + n_num;
392 fprintf(fpout, "%-10s %d\n", grpname, nsol_last);
396 fprintf(fpout, "%-15s %d\n", p_name, p_num);
400 fprintf(fpout, "%-15s %d\n", n_name, n_num);
405 make_backup(topinout);
406 gmx_file_rename(temporary_filename, topinout);
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
415 static std::vector<int> invertIndexGroup(int nrAtoms, std::vector<int> indexGroup)
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());
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)
428 const int firstToAddToInvertedGroup = *indexGroupIt + 1;
429 const int numIndicesToAdd = *(indexGroupIt + 1) - firstToAddToInvertedGroup;
430 if (numIndicesToAdd > 0)
432 invertedGroup.resize(invertedGroup.size() + numIndicesToAdd);
433 std::iota(std::end(invertedGroup) - numIndicesToAdd, std::end(invertedGroup), firstToAddToInvertedGroup);
437 return invertedGroup;
440 int gmx_genion(int argc, char* argv[])
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]."
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.",
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;
465 gmx_bool bNeutral = FALSE;
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)" },
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." },
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]. " }
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)
505 if (!parse_common_args(
506 &argc, argv, 0, NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs, &oenv))
510 output_env_done(oenv);
515 /* Check input for something sensible */
516 if ((p_num < 0) || (n_num < 0))
518 gmx_fatal(FARGS, "Negative number of ions to add?");
521 if (conc > 0 && (p_num > 0 || n_num > 0))
523 fprintf(stderr, "WARNING: -conc specified, overriding -nn and -np.\n");
526 /* Read atom positions and charges */
527 read_tps_conf(ftp2fn(efTPR, NFILE, fnm), &top, &pbcType, &x, nullptr, box, FALSE);
530 /* Compute total charge */
532 for (int i = 0; (i < atoms.nr); i++)
534 qtot += atoms.atom[i].q;
536 iqtot = gmx::roundToInt(qtot);
541 /* Compute number of ions to be added */
543 nsalt = gmx::roundToInt(conc * vol * AVOGADRO / 1e24);
544 p_num = abs(nsalt * n_q);
545 n_num = abs(nsalt * p_q);
549 int qdelta = p_num * p_q + n_num * n_q + iqtot;
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)
557 "Can't neutralize this system using -nq %d and"
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);
583 if ((p_num == 0) && (n_num == 0))
585 fprintf(stderr, "No ions to add, will just copy input configuration.\n");
589 char* grpname = nullptr;
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");
594 std::vector<int> solventGroup;
596 int* index = nullptr;
598 get_index(&atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nwa, &index, &grpname);
599 solventGroup.assign(index, index + nwa);
603 for (gmx::index i = 1; i < gmx::ssize(solventGroup); i++)
605 if (solventGroup[i] != solventGroup[i - 1] + 1)
608 "The solvent group %s is not continuous: "
609 "index[%d]=%d, index[%d]=%d",
612 solventGroup[i - 1] + 1,
614 solventGroup[i] + 1);
618 while ((nsa < gmx::ssize(solventGroup))
619 && (atoms.atom[solventGroup[nsa]].resind == atoms.atom[solventGroup[nsa - 1]].resind))
623 if (solventGroup.size() % nsa != 0)
626 "Your solvent group size (%td) is not a multiple of %d",
627 gmx::ssize(solventGroup),
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)
634 gmx_fatal(FARGS, "Not enough solvent for adding ions");
637 if (opt2bSet("-p", NFILE, fnm))
639 update_topol(opt2fn("-p", NFILE, fnm), p_num, n_num, p_name, n_name, grpname);
643 set_pbc(&pbc, pbcType, box);
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());
651 fprintf(stderr, "Using random seed %d.\n", seed);
654 std::vector<int> notSolventGroup = invertIndexGroup(atoms.nr, solventGroup);
656 std::vector<int> solventMoleculesForReplacement(nw);
657 std::iota(std::begin(solventMoleculesForReplacement), std::end(solventMoleculesForReplacement), 0);
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);
663 std::begin(solventMoleculesForReplacement), std::end(solventMoleculesForReplacement), rng);
665 /* Now loop over the ions that have to be placed */
669 nsa, &solventMoleculesForReplacement, repl, solventGroup, x, &pbc, 1, p_q, p_name, &atoms, rmin, ¬SolventGroup);
674 nsa, &solventMoleculesForReplacement, repl, solventGroup, x, &pbc, -1, n_q, n_name, &atoms, rmin, ¬SolventGroup);
676 fprintf(stderr, "\n");
680 sort_ions(nsa, nw, repl, solventGroup, &atoms, x, &pptr, &nptr, &paptr, &naptr);
687 sfree(atoms.pdbinfo);
688 atoms.pdbinfo = nullptr;
689 write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *top.name, &atoms, x, nullptr, pbcType, box);
697 output_env_done(oenv);