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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * 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.
51 #include "gmx_fatal.h"
63 #include "mtop_util.h"
66 static void insert_ion(int nsa,int *nwater,
67 gmx_bool bSet[],int repl[],atom_id index[],
68 real pot[],rvec x[],t_pbc *pbc,
69 int sign,int q,const char *ionname,
71 real rmin,gmx_bool bRandom,int *seed)
73 int i,ii,ei,owater,wlast,m,nw;
74 real extr_e,poti,rmin2;
77 gmx_large_int_t maxrand;
87 } while (bSet[ei] && (maxrand > 0));
89 gmx_fatal(FARGS,"No more replaceable solvent!");
93 for(i=0; (i<nw); i++) {
98 if ((poti <= extr_e) || !bSub) {
105 if ((poti >= extr_e) || !bSub) {
114 gmx_fatal(FARGS,"No more replaceable solvent!");
116 fprintf(stderr,"Replacing solvent molecule %d (atom %d) with %s\n",
117 ei,index[nsa*ei],ionname);
119 /* Replace solvent molecule charges with ion charge */
122 mdatoms->chargeA[index[nsa*ei]] = q;
124 mdatoms->chargeA[index[nsa*ei+i]] = 0;
126 /* Mark all solvent molecules within rmin as unavailable for substitution */
129 for(i=0; (i<nw); i++) {
131 pbc_dx(pbc,x[index[nsa*ei]],x[index[nsa*i]],dx);
132 if (iprod(dx,dx) < rmin2)
139 static char *aname(const char *mname)
146 while (i>1 && (isdigit(str[i]) || (str[i]=='+') || (str[i]=='-'))) {
154 void sort_ions(int nsa,int nw,int repl[],atom_id index[],
155 t_atoms *atoms,rvec x[],
156 const char *p_name,const char *n_name)
158 int i,j,k,r,np,nn,starta,startr,npi,nni;
160 char **pptr=NULL,**nptr=NULL,**paptr=NULL,**naptr=NULL;
164 /* Put all the solvent in front and count the added ions */
168 for(i=0; i<nw; i++) {
172 copy_rvec(x[index[nsa*i+k]],xt[j++]);
180 /* Put the positive and negative ions at the end */
181 starta = index[nsa*(nw - np - nn)];
182 startr = atoms->atom[starta].resind;
186 pptr[0] = strdup(p_name);
188 paptr[0] = aname(p_name);
192 nptr[0] = strdup(n_name);
194 naptr[0] = aname(n_name);
198 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;
211 copy_rvec(x[index[nsa*i]],xt[j]);
212 atoms->atomname[j] = naptr;
213 atoms->atom[j].resind = k;
214 atoms->resinfo[k].name = nptr;
218 for(i=index[nsa*nw-1]+1; i<atoms->nr; i++) {
219 j = i-(nsa-1)*(np+nn);
220 atoms->atomname[j] = atoms->atomname[i];
221 atoms->atom[j] = atoms->atom[i];
222 copy_rvec(x[i],xt[j]);
224 atoms->nr -= (nsa-1)*(np+nn);
226 /* Copy the new positions back */
227 for(i=index[0]; i<atoms->nr; i++)
228 copy_rvec(xt[i],x[i]);
233 static void update_topol(const char *topinout,int p_num,int n_num,
234 const char *p_name,const char *n_name,char *grpname)
236 #define TEMP_FILENM "temp.top"
238 char buf[STRLEN],buf2[STRLEN],*temp,**mol_line=NULL;
239 int line,i,nsol,nmol_line,sol_line,nsol_last;
242 printf("\nProcessing topology\n");
243 fpin = ffopen(topinout,"r");
244 fpout= ffopen(TEMP_FILENM,"w");
251 while (fgets(buf, STRLEN, fpin)) {
254 if ((temp=strchr(buf2,'\n')) != NULL)
259 if ((temp=strchr(buf2,'\n')) != NULL)
262 if (buf2[strlen(buf2)-1]==']') {
263 buf2[strlen(buf2)-1]='\0';
266 bMolecules=(gmx_strcasecmp(buf2,"molecules")==0);
268 fprintf(fpout,"%s",buf);
269 } else if (!bMolecules) {
270 fprintf(fpout,"%s",buf);
272 /* Check if this is a line with solvent molecules */
273 sscanf(buf,"%s",buf2);
274 if (gmx_strcasecmp(buf2,grpname) == 0) {
275 sol_line = nmol_line;
276 sscanf(buf,"%*s %d",&nsol_last);
278 /* Store this molecules section line */
279 srenew(mol_line,nmol_line+1);
280 mol_line[nmol_line] = strdup(buf);
286 if (sol_line == -1) {
288 gmx_fatal(FARGS,"No line with moleculetype '%s' found the [ molecules ] section of file '%s'",grpname,topinout);
290 if (nsol_last < p_num+n_num) {
292 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);
295 /* Print all the molecule entries */
296 for(i=0; i<nmol_line; i++) {
298 fprintf(fpout,"%s",mol_line[i]);
300 printf("Replacing %d solute molecules in topology file (%s) "
301 " by %d %s and %d %s ions.\n",
302 p_num+n_num,topinout,p_num,p_name,n_num,n_name);
303 nsol_last -= p_num + n_num;
305 fprintf(fpout,"%-10s %d\n",grpname,nsol_last);
308 fprintf(fpout,"%-15s %d\n",p_name,p_num);
311 fprintf(fpout,"%-15s %d\n",n_name,n_num);
316 /* use ffopen to generate backup of topinout */
317 fpout = ffopen(topinout,"w");
319 rename(TEMP_FILENM,topinout);
323 int gmx_genion(int argc, char *argv[])
325 const char *desc[] = {
326 "[TT]genion[tt] replaces solvent molecules by monoatomic ions at",
327 "the position of the first atoms with the most favorable electrostatic",
328 "potential or at random. The potential is calculated on all atoms, using",
329 "normal GROMACS particle-based methods (in contrast to other methods",
330 "based on solving the Poisson-Boltzmann equation).",
331 "The potential is recalculated after every ion insertion.",
332 "If specified in the run input file, a reaction field, shift function",
333 "or user function can be used. For the user function a table file",
334 "can be specified with the option [TT]-table[tt].",
335 "The group of solvent molecules should be continuous and all molecules",
336 "should have the same number of atoms.",
337 "The user should add the ion molecules to the topology file or use",
338 "the [TT]-p[tt] option to automatically modify the topology.[PAR]",
339 "The ion molecule type, residue and atom names in all force fields",
340 "are the capitalized element names without sign. This molecule name",
341 "should be given with [TT]-pname[tt] or [TT]-nname[tt], and the",
342 "[TT][molecules][tt] section of your topology updated accordingly,",
343 "either by hand or with [TT]-p[tt]. Do not use an atom name instead!",
344 "[PAR]Ions which can have multiple charge states get the multiplicity",
345 "added, without sign, for the uncommon states only.[PAR]",
346 "With the option [TT]-pot[tt] the potential can be written as B-factors",
347 "in a [TT].pdb[tt] file (for visualisation using e.g. Rasmol).",
348 "The unit of the potential is 1000 kJ/(mol e), the scaling be changed",
349 "with the [TT]-scale[tt] option.[PAR]",
350 "For larger ions, e.g. sulfate we recommended using [TT]genbox[tt]."
352 const char *bugs[] = {
353 "Calculation of the potential is not reliable, therefore the [TT]-random[tt] option is now turned on by default.",
354 "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."
356 static int p_num=0,n_num=0,p_q=1,n_q=-1;
357 static const char *p_name="NA",*n_name="CL";
358 static real rmin=0.6,scale=0.001,conc=0;
359 static int seed=1993;
360 static gmx_bool bRandom=TRUE,bNeutral=FALSE;
361 static t_pargs pa[] = {
362 { "-np", FALSE, etINT, {&p_num}, "Number of positive ions" },
363 { "-pname", FALSE, etSTR, {&p_name},"Name of the positive ion" },
364 { "-pq", FALSE, etINT, {&p_q}, "Charge of the positive ion" },
365 { "-nn", FALSE, etINT, {&n_num}, "Number of negative ions" },
366 { "-nname", FALSE, etSTR, {&n_name},"Name of the negative ion" },
367 { "-nq", FALSE, etINT, {&n_q}, "Charge of the negative ion" },
368 { "-rmin", FALSE, etREAL, {&rmin}, "Minimum distance between ions" },
369 { "-random",FALSE,etBOOL, {&bRandom},"Use random placement of ions instead of based on potential. The rmin option should still work" },
370 { "-seed", FALSE, etINT, {&seed}, "Seed for random number generator" },
371 { "-scale", FALSE, etREAL, {&scale}, "Scaling factor for the potential for [TT]-pot[tt]" },
372 { "-conc", FALSE, etREAL, {&conc},
373 "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." },
374 { "-neutral", FALSE, etBOOL, {&bNeutral},
375 "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." }
382 gmx_enerdata_t enerd;
394 int i,nw,nwa,nsa,nsalt,iqtot;
398 { efTPX, NULL, NULL, ffREAD },
399 { efXVG, "-table","table", ffOPTRD },
400 { efNDX, NULL, NULL, ffOPTRD },
401 { efSTO, "-o", NULL, ffWRITE },
402 { efLOG, "-g", "genion", ffWRITE },
403 { efPDB, "-pot", "pot", ffOPTWR },
404 { efTOP, "-p", "topol", ffOPTRW }
406 #define NFILE asize(fnm)
408 CopyRight(stderr,argv[0]);
409 parse_common_args(&argc,argv,PCA_BE_NICE,NFILE,fnm,asize(pa),pa,
410 asize(desc),desc, asize(bugs),bugs,&oenv);
411 bPDB = ftp2bSet(efPDB,NFILE,fnm);
412 if (bRandom && bPDB) {
413 fprintf(stderr,"Not computing potential with random option!\n");
417 /* Check input for something sensible */
418 if ((p_num<0) || (n_num<0))
419 gmx_fatal(FARGS,"Negative number of ions to add?");
423 fplog = init_calcpot(ftp2fn(efLOG,NFILE,fnm),ftp2fn(efTPX,NFILE,fnm),
424 opt2fn("-table",NFILE,fnm),mtop,top,&inputrec,&cr,
425 &graph,&mdatoms,&fr,&enerd,&pot,box,&x,oenv);
427 atoms = gmx_mtop_global_atoms(mtop);
430 for(i=0; (i<atoms.nr); i++)
431 qtot += atoms.atom[i].q;
432 iqtot = gmx_nint(qtot);
434 if ((conc > 0) || bNeutral) {
435 /* Compute number of ions to be added */
438 nsalt = gmx_nint(conc*vol*AVOGADRO/1e24);
439 p_num = abs(nsalt*n_q);
440 n_num = abs(nsalt*p_q);
444 qdelta = (p_num*p_q + n_num*n_q + iqtot);
446 p_num += abs(qdelta/p_q);
447 qdelta = (p_num*p_q + n_num*n_q + iqtot);
450 n_num += abs(qdelta/n_q);
451 qdelta = (p_num*p_q + n_num*n_q + iqtot);
453 } while (qdelta != 0);
458 if ((p_num == 0) && (n_num == 0)) {
460 fprintf(stderr,"No ions to add and no potential to calculate.\n");
464 nsa = 0; /* to keep gcc happy */
466 printf("Will try to add %d %s ions and %d %s ions.\n",
467 p_num,p_name,n_num,n_name);
468 printf("Select a continuous group of solvent molecules\n");
469 get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&nwa,&index,&grpname);
471 if (index[i] != index[i-1]+1)
472 gmx_fatal(FARGS,"The solvent group %s is not continuous: "
473 "index[%d]=%d, index[%d]=%d",
474 grpname,i,index[i-1]+1,i+1,index[i]+1);
477 (atoms.atom[index[nsa]].resind ==
478 atoms.atom[index[nsa-1]].resind))
481 gmx_fatal(FARGS,"Your solvent group size (%d) is not a multiple of %d",
484 fprintf(stderr,"Number of (%d-atomic) solvent molecules: %d\n",nsa,nw);
485 if (p_num+n_num > nw)
486 gmx_fatal(FARGS,"Not enough solvent for adding ions");
489 if (opt2bSet("-p",NFILE,fnm))
490 update_topol(opt2fn("-p",NFILE,fnm),p_num,n_num,p_name,n_name,grpname);
496 snew(atoms.pdbinfo,atoms.nr);
498 set_pbc(&pbc,inputrec.ePBC,box);
500 /* Now loop over the ions that have to be placed */
503 calc_pot(fplog,cr,mtop,&inputrec,top,x,fr,&enerd,mdatoms,pot,box,graph);
508 sprintf(buf,"%d_%s",p_num+n_num,ftp2fn(efPDB,NFILE,fnm));
510 strcpy(buf,ftp2fn(efPDB,NFILE,fnm));
511 for(i=0; (i<atoms.nr); i++)
512 atoms.pdbinfo[i].bfac = pot[i]*scale;
513 write_sto_conf(buf,"Potential calculated by genion",
514 &atoms,x,v,inputrec.ePBC,box);
518 if ((p_num > 0) && (p_num >= n_num)) {
519 insert_ion(nsa,&nw,bSet,repl,index,pot,x,&pbc,
520 1,p_q,p_name,mdatoms,rmin,bRandom,&seed);
523 else if (n_num > 0) {
524 insert_ion(nsa,&nw,bSet,repl,index,pot,x,&pbc,
525 -1,n_q,n_name,mdatoms,rmin,bRandom,&seed);
528 } while (p_num+n_num > 0);
529 fprintf(stderr,"\n");
532 sort_ions(nsa,nw,repl,index,&atoms,x,p_name,n_name);
534 sfree(atoms.pdbinfo);
535 atoms.pdbinfo = NULL;
536 write_sto_conf(ftp2fn(efSTO,NFILE,fnm),*mtop->name,&atoms,x,NULL,
541 gmx_log_close(fplog);