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, 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.
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.
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.
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.
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.
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.
44 #include "gromacs/utility/cstringutil.h"
45 #include "gromacs/utility/smalloc.h"
47 #include "gromacs/fileio/confio.h"
48 #include "gromacs/commandline/pargs.h"
51 #include "gmx_fatal.h"
52 #include "gromacs/fileio/futil.h"
53 #include "gromacs/math/utilities.h"
56 #include "gromacs/fileio/tpxio.h"
59 #include "gromacs/random/random.h"
61 #include "mtop_util.h"
64 static void insert_ion(int nsa, int *nwater,
65 gmx_bool bSet[], int repl[], atom_id index[],
67 int sign, int q, const char *ionname,
69 real rmin, gmx_rng_t rng)
83 ei = nw*gmx_rng_uniform_real(rng);
86 while (bSet[ei] && (maxrand > 0));
89 gmx_fatal(FARGS, "No more replaceable solvent!");
92 fprintf(stderr, "Replacing solvent molecule %d (atom %d) with %s\n",
93 ei, index[nsa*ei], ionname);
95 /* Replace solvent molecule charges with ion charge */
99 atoms->atom[index[nsa*ei]].q = q;
100 for (i = 1; i < nsa; i++)
102 atoms->atom[index[nsa*ei+i]].q = 0;
105 /* Mark all solvent molecules within rmin as unavailable for substitution */
109 for (i = 0; (i < nw); i++)
113 pbc_dx(pbc, x[index[nsa*ei]], x[index[nsa*i]], dx);
114 if (iprod(dx, dx) < rmin2)
124 static char *aname(const char *mname)
131 while (i > 1 && (isdigit(str[i]) || (str[i] == '+') || (str[i] == '-')))
140 void sort_ions(int nsa, int nw, int repl[], atom_id index[],
141 t_atoms *atoms, rvec x[],
142 const char *p_name, const char *n_name)
144 int i, j, k, r, np, nn, starta, startr, npi, nni;
146 char **pptr = NULL, **nptr = NULL, **paptr = NULL, **naptr = NULL;
150 /* Put all the solvent in front and count the added ions */
154 for (i = 0; i < nw; i++)
159 for (k = 0; k < nsa; k++)
161 copy_rvec(x[index[nsa*i+k]], xt[j++]);
176 /* Put the positive and negative ions at the end */
177 starta = index[nsa*(nw - np - nn)];
178 startr = atoms->atom[starta].resind;
183 pptr[0] = strdup(p_name);
185 paptr[0] = aname(p_name);
190 nptr[0] = strdup(n_name);
192 naptr[0] = aname(n_name);
196 for (i = 0; i < nw; i++)
203 copy_rvec(x[index[nsa*i]], xt[j]);
204 atoms->atomname[j] = paptr;
205 atoms->atom[j].resind = k;
206 atoms->resinfo[k].name = pptr;
213 copy_rvec(x[index[nsa*i]], xt[j]);
214 atoms->atomname[j] = naptr;
215 atoms->atom[j].resind = k;
216 atoms->resinfo[k].name = nptr;
220 for (i = index[nsa*nw-1]+1; i < atoms->nr; i++)
222 j = i-(nsa-1)*(np+nn);
223 atoms->atomname[j] = atoms->atomname[i];
224 atoms->atom[j] = atoms->atom[i];
225 copy_rvec(x[i], xt[j]);
227 atoms->nr -= (nsa-1)*(np+nn);
229 /* Copy the new positions back */
230 for (i = index[0]; i < atoms->nr; i++)
232 copy_rvec(xt[i], x[i]);
238 static void update_topol(const char *topinout, int p_num, int n_num,
239 const char *p_name, const char *n_name, char *grpname)
241 #define TEMP_FILENM "temp.top"
243 char buf[STRLEN], buf2[STRLEN], *temp, **mol_line = NULL;
244 int line, i, nsol, nmol_line, sol_line, nsol_last;
247 printf("\nProcessing topology\n");
248 fpin = gmx_ffopen(topinout, "r");
249 fpout = gmx_ffopen(TEMP_FILENM, "w");
256 while (fgets(buf, STRLEN, fpin))
260 if ((temp = strchr(buf2, '\n')) != NULL)
268 if ((temp = strchr(buf2, '\n')) != NULL)
273 if (buf2[strlen(buf2)-1] == ']')
275 buf2[strlen(buf2)-1] = '\0';
278 bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);
280 fprintf(fpout, "%s", buf);
282 else if (!bMolecules)
284 fprintf(fpout, "%s", buf);
288 /* Check if this is a line with solvent molecules */
289 sscanf(buf, "%s", buf2);
290 if (gmx_strcasecmp(buf2, grpname) == 0)
292 sol_line = nmol_line;
293 sscanf(buf, "%*s %d", &nsol_last);
295 /* Store this molecules section line */
296 srenew(mol_line, nmol_line+1);
297 mol_line[nmol_line] = strdup(buf);
306 gmx_fatal(FARGS, "No line with moleculetype '%s' found the [ molecules ] section of file '%s'", grpname, topinout);
308 if (nsol_last < p_num+n_num)
311 gmx_fatal(FARGS, "The last entry for moleculetype '%s' in the [ molecules ] section of file '%s' has less solvent molecules (%d) than were replaced (%d)", grpname, topinout, nsol_last, p_num+n_num);
314 /* Print all the molecule entries */
315 for (i = 0; i < nmol_line; i++)
319 fprintf(fpout, "%s", mol_line[i]);
323 printf("Replacing %d solute molecules in topology file (%s) "
324 " by %d %s and %d %s ions.\n",
325 p_num+n_num, topinout, p_num, p_name, n_num, n_name);
326 nsol_last -= p_num + n_num;
329 fprintf(fpout, "%-10s %d\n", grpname, nsol_last);
333 fprintf(fpout, "%-15s %d\n", p_name, p_num);
337 fprintf(fpout, "%-15s %d\n", n_name, n_num);
342 /* use gmx_ffopen to generate backup of topinout */
343 fpout = gmx_ffopen(topinout, "w");
345 rename(TEMP_FILENM, topinout);
349 int gmx_genion(int argc, char *argv[])
351 const char *desc[] = {
352 "[THISMODULE] randomly replaces solvent molecules with monoatomic ions.",
353 "The group of solvent molecules should be continuous and all molecules",
354 "should have the same number of atoms.",
355 "The user should add the ion molecules to the topology file or use",
356 "the [TT]-p[tt] option to automatically modify the topology.[PAR]",
357 "The ion molecule type, residue and atom names in all force fields",
358 "are the capitalized element names without sign. This molecule name",
359 "should be given with [TT]-pname[tt] or [TT]-nname[tt], and the",
360 "[TT][molecules][tt] section of your topology updated accordingly,",
361 "either by hand or with [TT]-p[tt]. Do not use an atom name instead!",
362 "[PAR]Ions which can have multiple charge states get the multiplicity",
363 "added, without sign, for the uncommon states only.[PAR]",
364 "For larger ions, e.g. sulfate we recommended using [gmx-insert-molecules]."
366 const char *bugs[] = {
367 "If you specify a salt concentration existing ions are not taken into "
368 "account. In effect you therefore specify the amount of salt to be added.",
370 static int p_num = 0, n_num = 0, p_q = 1, n_q = -1;
371 static const char *p_name = "NA", *n_name = "CL";
372 static real rmin = 0.6, conc = 0;
373 static int seed = 1993;
374 static gmx_bool bNeutral = FALSE;
375 static t_pargs pa[] = {
376 { "-np", FALSE, etINT, {&p_num}, "Number of positive ions" },
377 { "-pname", FALSE, etSTR, {&p_name}, "Name of the positive ion" },
378 { "-pq", FALSE, etINT, {&p_q}, "Charge of the positive ion" },
379 { "-nn", FALSE, etINT, {&n_num}, "Number of negative ions" },
380 { "-nname", FALSE, etSTR, {&n_name}, "Name of the negative ion" },
381 { "-nq", FALSE, etINT, {&n_q}, "Charge of the negative ion" },
382 { "-rmin", FALSE, etREAL, {&rmin}, "Minimum distance between ions" },
383 { "-seed", FALSE, etINT, {&seed}, "Seed for random number generator" },
384 { "-conc", FALSE, etREAL, {&conc},
385 "Specify salt concentration (mol/liter). This will add sufficient ions to reach up to the specified concentration as computed from the volume of the cell in the input [TT].tpr[tt] file. Overrides the [TT]-np[tt] and [TT]-nn[tt] options." },
386 { "-neutral", FALSE, etBOOL, {&bNeutral}, "This option will add enough ions to neutralize the system. These ions are added on top of those specified with [TT]-np[tt]/[TT]-nn[tt] or [TT]-conc[tt]. "}
396 char *grpname, title[STRLEN];
398 int i, nw, nwa, nsa, nsalt, iqtot;
402 { efTPX, NULL, NULL, ffREAD },
403 { efNDX, NULL, NULL, ffOPTRD },
404 { efSTO, "-o", NULL, ffWRITE },
405 { efTOP, "-p", "topol", ffOPTRW }
407 #define NFILE asize(fnm)
409 if (!parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa), pa,
410 asize(desc), desc, asize(bugs), bugs, &oenv))
415 /* Check input for something sensible */
416 if ((p_num < 0) || (n_num < 0))
418 gmx_fatal(FARGS, "Negative number of ions to add?");
421 if (conc > 0 && (p_num > 0 || n_num > 0))
423 fprintf(stderr, "WARNING: -conc specified, overriding -nn and -np.\n");
426 /* Read atom positions and charges */
427 read_tps_conf(ftp2fn(efTPX, NFILE, fnm), title, &top, &ePBC, &x, &v, box, FALSE);
430 /* Compute total charge */
432 for (i = 0; (i < atoms.nr); i++)
434 qtot += atoms.atom[i].q;
436 iqtot = gmx_nint(qtot);
441 /* Compute number of ions to be added */
443 nsalt = gmx_nint(conc*vol*AVOGADRO/1e24);
444 p_num = abs(nsalt*n_q);
445 n_num = abs(nsalt*p_q);
449 int qdelta = p_num*p_q + n_num*n_q + iqtot;
451 /* Check if the system is neutralizable
452 * is (qdelta == p_q*p_num + n_q*n_num) solvable for p_num and n_num? */
453 int gcd = gmx_greatest_common_divisor(n_q, p_q);
454 if ((qdelta % gcd) != 0)
456 gmx_fatal(FARGS, "Can't neutralize this system using -nq %d and"
457 " -pq %d.\n", n_q, p_q);
475 if ((p_num == 0) && (n_num == 0))
477 fprintf(stderr, "No ions to add.\n");
482 printf("Will try to add %d %s ions and %d %s ions.\n",
483 p_num, p_name, n_num, n_name);
484 printf("Select a continuous group of solvent molecules\n");
485 get_index(&atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nwa, &index, &grpname);
486 for (i = 1; i < nwa; i++)
488 if (index[i] != index[i-1]+1)
490 gmx_fatal(FARGS, "The solvent group %s is not continuous: "
491 "index[%d]=%d, index[%d]=%d",
492 grpname, i, index[i-1]+1, i+1, index[i]+1);
496 while ((nsa < nwa) &&
497 (atoms.atom[index[nsa]].resind ==
498 atoms.atom[index[nsa-1]].resind))
504 gmx_fatal(FARGS, "Your solvent group size (%d) is not a multiple of %d",
508 fprintf(stderr, "Number of (%d-atomic) solvent molecules: %d\n", nsa, nw);
509 if (p_num+n_num > nw)
511 gmx_fatal(FARGS, "Not enough solvent for adding ions");
515 if (opt2bSet("-p", NFILE, fnm))
517 update_topol(opt2fn("-p", NFILE, fnm), p_num, n_num, p_name, n_name, grpname);
524 snew(atoms.pdbinfo, atoms.nr);
526 set_pbc(&pbc, ePBC, box);
530 rng = gmx_rng_init(gmx_rng_make_seed());
534 rng = gmx_rng_init(seed);
536 /* Now loop over the ions that have to be placed */
539 insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
540 1, p_q, p_name, &atoms, rmin, rng);
544 insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
545 -1, n_q, n_name, &atoms, rmin, rng);
547 gmx_rng_destroy(rng);
548 fprintf(stderr, "\n");
552 sort_ions(nsa, nw, repl, index, &atoms, x, p_name, n_name);
555 sfree(atoms.pdbinfo);
556 atoms.pdbinfo = NULL;
557 write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *top.name, &atoms, x, NULL, ePBC,