3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
48 #include "gmx_fatal.h"
60 #include "mtop_util.h"
63 static int greatest_common_divisor(int p, int q)
75 static void insert_ion(int nsa, int *nwater,
76 gmx_bool bSet[], int repl[], atom_id index[],
77 real pot[], rvec x[], t_pbc *pbc,
78 int sign, int q, const char *ionname,
80 real rmin, gmx_bool bRandom, int *seed)
82 int i, ii, ei, owater, wlast, m, nw;
83 real extr_e, poti, rmin2;
85 gmx_bool bSub = FALSE;
86 gmx_large_int_t maxrand;
99 while (bSet[ei] && (maxrand > 0));
102 gmx_fatal(FARGS, "No more replaceable solvent!");
108 for (i = 0; (i < nw); i++)
116 if ((poti <= extr_e) || !bSub)
125 if ((poti >= extr_e) || !bSub)
136 gmx_fatal(FARGS, "No more replaceable solvent!");
139 fprintf(stderr, "Replacing solvent molecule %d (atom %d) with %s\n",
140 ei, index[nsa*ei], ionname);
142 /* Replace solvent molecule charges with ion charge */
145 mdatoms->chargeA[index[nsa*ei]] = q;
146 for (i = 1; i < nsa; i++)
148 mdatoms->chargeA[index[nsa*ei+i]] = 0;
151 /* Mark all solvent molecules within rmin as unavailable for substitution */
155 for (i = 0; (i < nw); i++)
159 pbc_dx(pbc, x[index[nsa*ei]], x[index[nsa*i]], dx);
160 if (iprod(dx, dx) < rmin2)
169 static char *aname(const char *mname)
176 while (i > 1 && (isdigit(str[i]) || (str[i] == '+') || (str[i] == '-')))
185 void sort_ions(int nsa, int nw, int repl[], atom_id index[],
186 t_atoms *atoms, rvec x[],
187 const char *p_name, const char *n_name)
189 int i, j, k, r, np, nn, starta, startr, npi, nni;
191 char **pptr = NULL, **nptr = NULL, **paptr = NULL, **naptr = NULL;
195 /* Put all the solvent in front and count the added ions */
199 for (i = 0; i < nw; i++)
204 for (k = 0; k < nsa; k++)
206 copy_rvec(x[index[nsa*i+k]], xt[j++]);
221 /* Put the positive and negative ions at the end */
222 starta = index[nsa*(nw - np - nn)];
223 startr = atoms->atom[starta].resind;
228 pptr[0] = strdup(p_name);
230 paptr[0] = aname(p_name);
235 nptr[0] = strdup(n_name);
237 naptr[0] = aname(n_name);
241 for (i = 0; i < nw; i++)
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;
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;
265 for (i = index[nsa*nw-1]+1; i < atoms->nr; i++)
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]);
272 atoms->nr -= (nsa-1)*(np+nn);
274 /* Copy the new positions back */
275 for (i = index[0]; i < atoms->nr; i++)
277 copy_rvec(xt[i], x[i]);
283 static void update_topol(const char *topinout, int p_num, int n_num,
284 const char *p_name, const char *n_name, char *grpname)
286 #define TEMP_FILENM "temp.top"
288 char buf[STRLEN], buf2[STRLEN], *temp, **mol_line = NULL;
289 int line, i, nsol, nmol_line, sol_line, nsol_last;
292 printf("\nProcessing topology\n");
293 fpin = ffopen(topinout, "r");
294 fpout = ffopen(TEMP_FILENM, "w");
301 while (fgets(buf, STRLEN, fpin))
305 if ((temp = strchr(buf2, '\n')) != NULL)
313 if ((temp = strchr(buf2, '\n')) != NULL)
318 if (buf2[strlen(buf2)-1] == ']')
320 buf2[strlen(buf2)-1] = '\0';
323 bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);
325 fprintf(fpout, "%s", buf);
327 else if (!bMolecules)
329 fprintf(fpout, "%s", buf);
333 /* Check if this is a line with solvent molecules */
334 sscanf(buf, "%s", buf2);
335 if (gmx_strcasecmp(buf2, grpname) == 0)
337 sol_line = nmol_line;
338 sscanf(buf, "%*s %d", &nsol_last);
340 /* Store this molecules section line */
341 srenew(mol_line, nmol_line+1);
342 mol_line[nmol_line] = strdup(buf);
351 gmx_fatal(FARGS, "No line with moleculetype '%s' found the [ molecules ] section of file '%s'", grpname, topinout);
353 if (nsol_last < p_num+n_num)
356 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);
359 /* Print all the molecule entries */
360 for (i = 0; i < nmol_line; i++)
364 fprintf(fpout, "%s", mol_line[i]);
368 printf("Replacing %d solute molecules in topology file (%s) "
369 " by %d %s and %d %s ions.\n",
370 p_num+n_num, topinout, p_num, p_name, n_num, n_name);
371 nsol_last -= p_num + n_num;
374 fprintf(fpout, "%-10s %d\n", grpname, nsol_last);
378 fprintf(fpout, "%-15s %d\n", p_name, p_num);
382 fprintf(fpout, "%-15s %d\n", n_name, n_num);
387 /* use ffopen to generate backup of topinout */
388 fpout = ffopen(topinout, "w");
390 rename(TEMP_FILENM, topinout);
394 int gmx_genion(int argc, char *argv[])
396 const char *desc[] = {
397 "[TT]genion[tt] replaces solvent molecules by monoatomic ions at",
398 "the position of the first atoms with the most favorable electrostatic",
399 "potential or at random. The potential is calculated on all atoms, using",
400 "normal GROMACS particle-based methods (in contrast to other methods",
401 "based on solving the Poisson-Boltzmann equation).",
402 "The potential is recalculated after every ion insertion.",
403 "If specified in the run input file, a reaction field, shift function",
404 "or user function can be used. For the user function a table file",
405 "can be specified with the option [TT]-table[tt].",
406 "The group of solvent molecules should be continuous and all molecules",
407 "should have the same number of atoms.",
408 "The user should add the ion molecules to the topology file or use",
409 "the [TT]-p[tt] option to automatically modify the topology.[PAR]",
410 "The ion molecule type, residue and atom names in all force fields",
411 "are the capitalized element names without sign. This molecule name",
412 "should be given with [TT]-pname[tt] or [TT]-nname[tt], and the",
413 "[TT][molecules][tt] section of your topology updated accordingly,",
414 "either by hand or with [TT]-p[tt]. Do not use an atom name instead!",
415 "[PAR]Ions which can have multiple charge states get the multiplicity",
416 "added, without sign, for the uncommon states only.[PAR]",
417 "With the option [TT]-pot[tt] the potential can be written as B-factors",
418 "in a [TT].pdb[tt] file (for visualisation using e.g. Rasmol).",
419 "The unit of the potential is 1000 kJ/(mol e), the scaling be changed",
420 "with the [TT]-scale[tt] option.[PAR]",
421 "For larger ions, e.g. sulfate we recommended using [TT]genbox[tt]."
423 const char *bugs[] = {
424 "Calculation of the potential is not reliable, therefore the [TT]-random[tt] option is now turned on by default.",
425 "If you specify a salt concentration existing ions are not taken into account. In effect you therefore specify the amount of salt to be added."
427 static int p_num = 0, n_num = 0, p_q = 1, n_q = -1;
428 static const char *p_name = "NA", *n_name = "CL";
429 static real rmin = 0.6, scale = 0.001, conc = 0;
430 static int seed = 1993;
431 static gmx_bool bRandom = TRUE, bNeutral = FALSE;
432 static t_pargs pa[] = {
433 { "-np", FALSE, etINT, {&p_num}, "Number of positive ions" },
434 { "-pname", FALSE, etSTR, {&p_name}, "Name of the positive ion" },
435 { "-pq", FALSE, etINT, {&p_q}, "Charge of the positive ion" },
436 { "-nn", FALSE, etINT, {&n_num}, "Number of negative ions" },
437 { "-nname", FALSE, etSTR, {&n_name}, "Name of the negative ion" },
438 { "-nq", FALSE, etINT, {&n_q}, "Charge of the negative ion" },
439 { "-rmin", FALSE, etREAL, {&rmin}, "Minimum distance between ions" },
440 { "-random", FALSE, etBOOL, {&bRandom}, "Use random placement of ions instead of based on potential. The rmin option should still work" },
441 { "-seed", FALSE, etINT, {&seed}, "Seed for random number generator" },
442 { "-scale", FALSE, etREAL, {&scale}, "Scaling factor for the potential for [TT]-pot[tt]" },
443 { "-conc", FALSE, etREAL, {&conc},
444 "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." },
445 { "-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]. "}
452 gmx_enerdata_t enerd;
456 real *pot, vol, qtot;
463 gmx_bool *bSet, bPDB;
464 int i, nw, nwa, nsa, nsalt, iqtot;
468 { efTPX, NULL, NULL, ffREAD },
469 { efXVG, "-table", "table", ffOPTRD },
470 { efNDX, NULL, NULL, ffOPTRD },
471 { efSTO, "-o", NULL, ffWRITE },
472 { efLOG, "-g", "genion", ffWRITE },
473 { efPDB, "-pot", "pot", ffOPTWR },
474 { efTOP, "-p", "topol", ffOPTRW }
476 #define NFILE asize(fnm)
478 parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa), pa,
479 asize(desc), desc, asize(bugs), bugs, &oenv);
480 bPDB = ftp2bSet(efPDB, NFILE, fnm);
483 fprintf(stderr, "Not computing potential with random option!\n");
487 /* Check input for something sensible */
488 if ((p_num < 0) || (n_num < 0))
490 gmx_fatal(FARGS, "Negative number of ions to add?");
495 fplog = init_calcpot(ftp2fn(efLOG, NFILE, fnm), ftp2fn(efTPX, NFILE, fnm),
496 opt2fn("-table", NFILE, fnm), mtop, top, &inputrec, &cr,
497 &graph, &mdatoms, &fr, &enerd, &pot, box, &x, oenv);
499 atoms = gmx_mtop_global_atoms(mtop);
502 for (i = 0; (i < atoms.nr); i++)
504 qtot += atoms.atom[i].q;
506 iqtot = gmx_nint(qtot);
511 /* Compute number of ions to be added */
513 nsalt = gmx_nint(conc*vol*AVOGADRO/1e24);
514 p_num = abs(nsalt*n_q);
515 n_num = abs(nsalt*p_q);
519 int qdelta = p_num*p_q + n_num*n_q + iqtot;
521 /* Check if the system is neutralizable
522 * is (qdelta == p_q*p_num + n_q*n_num) solvable for p_num and n_num? */
523 int gcd = greatest_common_divisor(n_q, p_q);
524 if ((qdelta % gcd) != 0)
526 gmx_fatal(FARGS, "Can't neutralize this system using -nq %d and"
527 " -pq %d.\n", n_q, p_q);
545 if ((p_num == 0) && (n_num == 0))
549 fprintf(stderr, "No ions to add and no potential to calculate.\n");
553 nsa = 0; /* to keep gcc happy */
557 printf("Will try to add %d %s ions and %d %s ions.\n",
558 p_num, p_name, n_num, n_name);
559 printf("Select a continuous group of solvent molecules\n");
560 get_index(&atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nwa, &index, &grpname);
561 for (i = 1; i < nwa; i++)
563 if (index[i] != index[i-1]+1)
565 gmx_fatal(FARGS, "The solvent group %s is not continuous: "
566 "index[%d]=%d, index[%d]=%d",
567 grpname, i, index[i-1]+1, i+1, index[i]+1);
571 while ((nsa < nwa) &&
572 (atoms.atom[index[nsa]].resind ==
573 atoms.atom[index[nsa-1]].resind))
579 gmx_fatal(FARGS, "Your solvent group size (%d) is not a multiple of %d",
583 fprintf(stderr, "Number of (%d-atomic) solvent molecules: %d\n", nsa, nw);
584 if (p_num+n_num > nw)
586 gmx_fatal(FARGS, "Not enough solvent for adding ions");
590 if (opt2bSet("-p", NFILE, fnm))
592 update_topol(opt2fn("-p", NFILE, fnm), p_num, n_num, p_name, n_name, grpname);
599 snew(atoms.pdbinfo, atoms.nr);
601 set_pbc(&pbc, inputrec.ePBC, box);
603 /* Now loop over the ions that have to be placed */
608 calc_pot(fplog, cr, mtop, &inputrec, top, x, fr, &enerd, mdatoms, pot, box, graph);
615 sprintf(buf, "%d_%s", p_num+n_num, ftp2fn(efPDB, NFILE, fnm));
619 strcpy(buf, ftp2fn(efPDB, NFILE, fnm));
621 for (i = 0; (i < atoms.nr); i++)
623 atoms.pdbinfo[i].bfac = pot[i]*scale;
625 write_sto_conf(buf, "Potential calculated by genion",
626 &atoms, x, v, inputrec.ePBC, box);
630 if ((p_num > 0) && (p_num >= n_num))
632 insert_ion(nsa, &nw, bSet, repl, index, pot, x, &pbc,
633 1, p_q, p_name, mdatoms, rmin, bRandom, &seed);
638 insert_ion(nsa, &nw, bSet, repl, index, pot, x, &pbc,
639 -1, n_q, n_name, mdatoms, rmin, bRandom, &seed);
643 while (p_num+n_num > 0);
644 fprintf(stderr, "\n");
648 sort_ions(nsa, nw, repl, index, &atoms, x, p_name, n_name);
651 sfree(atoms.pdbinfo);
652 atoms.pdbinfo = NULL;
653 write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *mtop->name, &atoms, x, NULL,
658 gmx_log_close(fplog);