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 void insert_ion(int nsa, int *nwater,
64 gmx_bool bSet[], int repl[], atom_id index[],
65 real pot[], rvec x[], t_pbc *pbc,
66 int sign, int q, const char *ionname,
68 real rmin, gmx_bool bRandom, int *seed)
70 int i, ii, ei, owater, wlast, m, nw;
71 real extr_e, poti, rmin2;
73 gmx_bool bSub = FALSE;
74 gmx_large_int_t maxrand;
87 while (bSet[ei] && (maxrand > 0));
90 gmx_fatal(FARGS, "No more replaceable solvent!");
96 for (i = 0; (i < nw); i++)
104 if ((poti <= extr_e) || !bSub)
113 if ((poti >= extr_e) || !bSub)
124 gmx_fatal(FARGS, "No more replaceable solvent!");
127 fprintf(stderr, "Replacing solvent molecule %d (atom %d) with %s\n",
128 ei, index[nsa*ei], ionname);
130 /* Replace solvent molecule charges with ion charge */
133 mdatoms->chargeA[index[nsa*ei]] = q;
134 for (i = 1; i < nsa; i++)
136 mdatoms->chargeA[index[nsa*ei+i]] = 0;
139 /* Mark all solvent molecules within rmin as unavailable for substitution */
143 for (i = 0; (i < nw); i++)
147 pbc_dx(pbc, x[index[nsa*ei]], x[index[nsa*i]], dx);
148 if (iprod(dx, dx) < rmin2)
157 static char *aname(const char *mname)
164 while (i > 1 && (isdigit(str[i]) || (str[i] == '+') || (str[i] == '-')))
173 void sort_ions(int nsa, int nw, int repl[], atom_id index[],
174 t_atoms *atoms, rvec x[],
175 const char *p_name, const char *n_name)
177 int i, j, k, r, np, nn, starta, startr, npi, nni;
179 char **pptr = NULL, **nptr = NULL, **paptr = NULL, **naptr = NULL;
183 /* Put all the solvent in front and count the added ions */
187 for (i = 0; i < nw; i++)
192 for (k = 0; k < nsa; k++)
194 copy_rvec(x[index[nsa*i+k]], xt[j++]);
209 /* Put the positive and negative ions at the end */
210 starta = index[nsa*(nw - np - nn)];
211 startr = atoms->atom[starta].resind;
216 pptr[0] = strdup(p_name);
218 paptr[0] = aname(p_name);
223 nptr[0] = strdup(n_name);
225 naptr[0] = aname(n_name);
229 for (i = 0; i < nw; i++)
236 copy_rvec(x[index[nsa*i]], xt[j]);
237 atoms->atomname[j] = paptr;
238 atoms->atom[j].resind = k;
239 atoms->resinfo[k].name = pptr;
246 copy_rvec(x[index[nsa*i]], xt[j]);
247 atoms->atomname[j] = naptr;
248 atoms->atom[j].resind = k;
249 atoms->resinfo[k].name = nptr;
253 for (i = index[nsa*nw-1]+1; i < atoms->nr; i++)
255 j = i-(nsa-1)*(np+nn);
256 atoms->atomname[j] = atoms->atomname[i];
257 atoms->atom[j] = atoms->atom[i];
258 copy_rvec(x[i], xt[j]);
260 atoms->nr -= (nsa-1)*(np+nn);
262 /* Copy the new positions back */
263 for (i = index[0]; i < atoms->nr; i++)
265 copy_rvec(xt[i], x[i]);
271 static void update_topol(const char *topinout, int p_num, int n_num,
272 const char *p_name, const char *n_name, char *grpname)
274 #define TEMP_FILENM "temp.top"
276 char buf[STRLEN], buf2[STRLEN], *temp, **mol_line = NULL;
277 int line, i, nsol, nmol_line, sol_line, nsol_last;
280 printf("\nProcessing topology\n");
281 fpin = ffopen(topinout, "r");
282 fpout = ffopen(TEMP_FILENM, "w");
289 while (fgets(buf, STRLEN, fpin))
293 if ((temp = strchr(buf2, '\n')) != NULL)
301 if ((temp = strchr(buf2, '\n')) != NULL)
306 if (buf2[strlen(buf2)-1] == ']')
308 buf2[strlen(buf2)-1] = '\0';
311 bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);
313 fprintf(fpout, "%s", buf);
315 else if (!bMolecules)
317 fprintf(fpout, "%s", buf);
321 /* Check if this is a line with solvent molecules */
322 sscanf(buf, "%s", buf2);
323 if (gmx_strcasecmp(buf2, grpname) == 0)
325 sol_line = nmol_line;
326 sscanf(buf, "%*s %d", &nsol_last);
328 /* Store this molecules section line */
329 srenew(mol_line, nmol_line+1);
330 mol_line[nmol_line] = strdup(buf);
339 gmx_fatal(FARGS, "No line with moleculetype '%s' found the [ molecules ] section of file '%s'", grpname, topinout);
341 if (nsol_last < p_num+n_num)
344 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);
347 /* Print all the molecule entries */
348 for (i = 0; i < nmol_line; i++)
352 fprintf(fpout, "%s", mol_line[i]);
356 printf("Replacing %d solute molecules in topology file (%s) "
357 " by %d %s and %d %s ions.\n",
358 p_num+n_num, topinout, p_num, p_name, n_num, n_name);
359 nsol_last -= p_num + n_num;
362 fprintf(fpout, "%-10s %d\n", grpname, nsol_last);
366 fprintf(fpout, "%-15s %d\n", p_name, p_num);
370 fprintf(fpout, "%-15s %d\n", n_name, n_num);
375 /* use ffopen to generate backup of topinout */
376 fpout = ffopen(topinout, "w");
378 rename(TEMP_FILENM, topinout);
382 int gmx_genion(int argc, char *argv[])
384 const char *desc[] = {
385 "[TT]genion[tt] replaces solvent molecules by monoatomic ions at",
386 "the position of the first atoms with the most favorable electrostatic",
387 "potential or at random. The potential is calculated on all atoms, using",
388 "normal GROMACS particle-based methods (in contrast to other methods",
389 "based on solving the Poisson-Boltzmann equation).",
390 "The potential is recalculated after every ion insertion.",
391 "If specified in the run input file, a reaction field, shift function",
392 "or user function can be used. For the user function a table file",
393 "can be specified with the option [TT]-table[tt].",
394 "The group of solvent molecules should be continuous and all molecules",
395 "should have the same number of atoms.",
396 "The user should add the ion molecules to the topology file or use",
397 "the [TT]-p[tt] option to automatically modify the topology.[PAR]",
398 "The ion molecule type, residue and atom names in all force fields",
399 "are the capitalized element names without sign. This molecule name",
400 "should be given with [TT]-pname[tt] or [TT]-nname[tt], and the",
401 "[TT][molecules][tt] section of your topology updated accordingly,",
402 "either by hand or with [TT]-p[tt]. Do not use an atom name instead!",
403 "[PAR]Ions which can have multiple charge states get the multiplicity",
404 "added, without sign, for the uncommon states only.[PAR]",
405 "With the option [TT]-pot[tt] the potential can be written as B-factors",
406 "in a [TT].pdb[tt] file (for visualisation using e.g. Rasmol).",
407 "The unit of the potential is 1000 kJ/(mol e), the scaling be changed",
408 "with the [TT]-scale[tt] option.[PAR]",
409 "For larger ions, e.g. sulfate we recommended using [TT]genbox[tt]."
411 const char *bugs[] = {
412 "Calculation of the potential is not reliable, therefore the [TT]-random[tt] option is now turned on by default.",
413 "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."
415 static int p_num = 0, n_num = 0, p_q = 1, n_q = -1;
416 static const char *p_name = "NA", *n_name = "CL";
417 static real rmin = 0.6, scale = 0.001, conc = 0;
418 static int seed = 1993;
419 static gmx_bool bRandom = TRUE, bNeutral = FALSE;
420 static t_pargs pa[] = {
421 { "-np", FALSE, etINT, {&p_num}, "Number of positive ions" },
422 { "-pname", FALSE, etSTR, {&p_name}, "Name of the positive ion" },
423 { "-pq", FALSE, etINT, {&p_q}, "Charge of the positive ion" },
424 { "-nn", FALSE, etINT, {&n_num}, "Number of negative ions" },
425 { "-nname", FALSE, etSTR, {&n_name}, "Name of the negative ion" },
426 { "-nq", FALSE, etINT, {&n_q}, "Charge of the negative ion" },
427 { "-rmin", FALSE, etREAL, {&rmin}, "Minimum distance between ions" },
428 { "-random", FALSE, etBOOL, {&bRandom}, "Use random placement of ions instead of based on potential. The rmin option should still work" },
429 { "-seed", FALSE, etINT, {&seed}, "Seed for random number generator" },
430 { "-scale", FALSE, etREAL, {&scale}, "Scaling factor for the potential for [TT]-pot[tt]" },
431 { "-conc", FALSE, etREAL, {&conc},
432 "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." },
433 { "-neutral", FALSE, etBOOL, {&bNeutral},
434 "This option will add enough ions to neutralize the system. In combination with the concentration option a neutral system at a given salt concentration will be generated." }
441 gmx_enerdata_t enerd;
445 real *pot, vol, qtot;
452 gmx_bool *bSet, bPDB;
453 int i, nw, nwa, nsa, nsalt, iqtot;
457 { efTPX, NULL, NULL, ffREAD },
458 { efXVG, "-table", "table", ffOPTRD },
459 { efNDX, NULL, NULL, ffOPTRD },
460 { efSTO, "-o", NULL, ffWRITE },
461 { efLOG, "-g", "genion", ffWRITE },
462 { efPDB, "-pot", "pot", ffOPTWR },
463 { efTOP, "-p", "topol", ffOPTRW }
465 #define NFILE asize(fnm)
467 parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa), pa,
468 asize(desc), desc, asize(bugs), bugs, &oenv);
469 bPDB = ftp2bSet(efPDB, NFILE, fnm);
472 fprintf(stderr, "Not computing potential with random option!\n");
476 /* Check input for something sensible */
477 if ((p_num < 0) || (n_num < 0))
479 gmx_fatal(FARGS, "Negative number of ions to add?");
484 fplog = init_calcpot(ftp2fn(efLOG, NFILE, fnm), ftp2fn(efTPX, NFILE, fnm),
485 opt2fn("-table", NFILE, fnm), mtop, top, &inputrec, &cr,
486 &graph, &mdatoms, &fr, &enerd, &pot, box, &x, oenv);
488 atoms = gmx_mtop_global_atoms(mtop);
491 for (i = 0; (i < atoms.nr); i++)
493 qtot += atoms.atom[i].q;
495 iqtot = gmx_nint(qtot);
497 if ((conc > 0) || bNeutral)
499 /* Compute number of ions to be added */
503 nsalt = gmx_nint(conc*vol*AVOGADRO/1e24);
504 p_num = abs(nsalt*n_q);
505 n_num = abs(nsalt*p_q);
511 qdelta = (p_num*p_q + n_num*n_q + iqtot);
514 p_num += abs(qdelta/p_q);
515 qdelta = (p_num*p_q + n_num*n_q + iqtot);
519 n_num += abs(qdelta/n_q);
520 qdelta = (p_num*p_q + n_num*n_q + iqtot);
528 if ((p_num == 0) && (n_num == 0))
532 fprintf(stderr, "No ions to add and no potential to calculate.\n");
536 nsa = 0; /* to keep gcc happy */
540 printf("Will try to add %d %s ions and %d %s ions.\n",
541 p_num, p_name, n_num, n_name);
542 printf("Select a continuous group of solvent molecules\n");
543 get_index(&atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nwa, &index, &grpname);
544 for (i = 1; i < nwa; i++)
546 if (index[i] != index[i-1]+1)
548 gmx_fatal(FARGS, "The solvent group %s is not continuous: "
549 "index[%d]=%d, index[%d]=%d",
550 grpname, i, index[i-1]+1, i+1, index[i]+1);
554 while ((nsa < nwa) &&
555 (atoms.atom[index[nsa]].resind ==
556 atoms.atom[index[nsa-1]].resind))
562 gmx_fatal(FARGS, "Your solvent group size (%d) is not a multiple of %d",
566 fprintf(stderr, "Number of (%d-atomic) solvent molecules: %d\n", nsa, nw);
567 if (p_num+n_num > nw)
569 gmx_fatal(FARGS, "Not enough solvent for adding ions");
573 if (opt2bSet("-p", NFILE, fnm))
575 update_topol(opt2fn("-p", NFILE, fnm), p_num, n_num, p_name, n_name, grpname);
582 snew(atoms.pdbinfo, atoms.nr);
584 set_pbc(&pbc, inputrec.ePBC, box);
586 /* Now loop over the ions that have to be placed */
591 calc_pot(fplog, cr, mtop, &inputrec, top, x, fr, &enerd, mdatoms, pot, box, graph);
598 sprintf(buf, "%d_%s", p_num+n_num, ftp2fn(efPDB, NFILE, fnm));
602 strcpy(buf, ftp2fn(efPDB, NFILE, fnm));
604 for (i = 0; (i < atoms.nr); i++)
606 atoms.pdbinfo[i].bfac = pot[i]*scale;
608 write_sto_conf(buf, "Potential calculated by genion",
609 &atoms, x, v, inputrec.ePBC, box);
613 if ((p_num > 0) && (p_num >= n_num))
615 insert_ion(nsa, &nw, bSet, repl, index, pot, x, &pbc,
616 1, p_q, p_name, mdatoms, rmin, bRandom, &seed);
621 insert_ion(nsa, &nw, bSet, repl, index, pot, x, &pbc,
622 -1, n_q, n_name, mdatoms, rmin, bRandom, &seed);
626 while (p_num+n_num > 0);
627 fprintf(stderr, "\n");
631 sort_ions(nsa, nw, repl, index, &atoms, x, p_name, n_name);
634 sfree(atoms.pdbinfo);
635 atoms.pdbinfo = NULL;
636 write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *mtop->name, &atoms, x, NULL,
641 gmx_log_close(fplog);