src/tools/gmx_genion.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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,2013, 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.
  11  *
  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.
  16  *
  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.
  21  *
  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.
  26  *
  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.
  34  *
  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.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <ctype.h>
  43 #include "copyrite.h"
  44 #include "string2.h"
  45 #include "smalloc.h"
  46 #include "sysstuff.h"
  47 #include "confio.h"
  48 #include "statutil.h"
  49 #include "pbc.h"
  50 #include "force.h"
  51 #include "gmx_fatal.h"
  52 #include "futil.h"
  53 #include "maths.h"
  54 #include "macros.h"
  55 #include "physics.h"
  56 #include "vec.h"
  57 #include "tpxio.h"
  58 #include "mdrun.h"
  59 #include "calcpot.h"
  60 #include "main.h"
  61 #include "random.h"
  62 #include "index.h"
  63 #include "mtop_util.h"
  64 #include "gmx_ana.h"
  65
  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,
  70                        t_mdatoms *mdatoms,
  71                        real rmin,gmx_bool bRandom,int *seed)
  72 {
  73   int  i,ii,ei,owater,wlast,m,nw;
  74   real extr_e,poti,rmin2;
  75   rvec xei,dx;
  76   gmx_bool bSub=FALSE;
  77   gmx_large_int_t maxrand;
  78
  79   ei=-1;
  80   nw = *nwater;
  81   maxrand  = nw;
  82   maxrand *= 1000;
  83   if (bRandom) {
  84     do {
  85       ei = nw*rando(seed);
  86       maxrand--;
  87     } while (bSet[ei] && (maxrand > 0));
  88     if (bSet[ei])
  89       gmx_fatal(FARGS,"No more replaceable solvent!");
  90   }
  91   else {
  92     extr_e = 0;
  93     for(i=0; (i<nw); i++) {
  94       if (!bSet[i]) {
  95         ii=index[nsa*i];
  96         poti=pot[ii];
  97         if (q > 0) {
  98           if ((poti <= extr_e) || !bSub) {
  99             extr_e = poti;
 100             ei     = i;
 101             bSub   = TRUE;
 102           }
 103         }
 104         else {
 105           if ((poti >= extr_e) || !bSub) {
 106             extr_e = poti;
 107             ei     = i;
 108             bSub   = TRUE;
 109           }
 110         }
 111       }
 112     }
 113     if (ei == -1)
 114       gmx_fatal(FARGS,"No more replaceable solvent!");
 115   }
 116   fprintf(stderr,"Replacing solvent molecule %d (atom %d) with %s\n",
 117           ei,index[nsa*ei],ionname);
 118
 119   /* Replace solvent molecule charges with ion charge */
 120   bSet[ei] = TRUE;
 121   repl[ei] = sign;
 122   mdatoms->chargeA[index[nsa*ei]] = q;
 123   for(i=1; i<nsa; i++)
 124     mdatoms->chargeA[index[nsa*ei+i]] = 0;
 125
 126   /* Mark all solvent molecules within rmin as unavailable for substitution */
 127   if (rmin > 0) {
 128     rmin2=rmin*rmin;
 129     for(i=0; (i<nw); i++) {
 130       if (!bSet[i]) {
 131         pbc_dx(pbc,x[index[nsa*ei]],x[index[nsa*i]],dx);
 132         if (iprod(dx,dx) < rmin2)
 133           bSet[i] = TRUE;
 134       }
 135     }
 136   }
 137 }
 138
 139 static char *aname(const char *mname)
 140 {
 141   char *str;
 142   int  i;
 143
 144   str = strdup(mname);
 145   i=strlen(str)-1;
 146   while (i>1 && (isdigit(str[i]) || (str[i]=='+') || (str[i]=='-'))) {
 147     str[i]='\0';
 148     i--;
 149   }
 150
 151   return str;
 152 }
 153
 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)
 157 {
 158   int i,j,k,r,np,nn,starta,startr,npi,nni;
 159   rvec *xt;
 160   char **pptr=NULL,**nptr=NULL,**paptr=NULL,**naptr=NULL;
 161
 162   snew(xt,atoms->nr);
 163
 164   /* Put all the solvent in front and count the added ions */
 165   np=0;
 166   nn=0;
 167   j=index[0];
 168   for(i=0; i<nw; i++) {
 169     r = repl[i];
 170     if (r == 0)
 171       for(k=0; k<nsa; k++)
 172         copy_rvec(x[index[nsa*i+k]],xt[j++]);
 173     else if (r>0)
 174       np++;
 175     else if (r<0)
 176       nn++;
 177   }
 178
 179   if (np+nn > 0) {
 180     /* Put the positive and negative ions at the end */
 181     starta = index[nsa*(nw - np - nn)];
 182     startr = atoms->atom[starta].resind;
 183
 184     if (np) {
 185       snew(pptr,1);
 186       pptr[0] = strdup(p_name);
 187       snew(paptr,1);
 188       paptr[0] = aname(p_name);
 189     }
 190     if (nn) {
 191       snew(nptr,1);
 192       nptr[0] = strdup(n_name);
 193       snew(naptr,1);
 194       naptr[0] = aname(n_name);
 195     }
 196     npi = 0;
 197     nni = 0;
 198     for(i=0; i<nw; i++) {
 199       r = repl[i];
 200       if (r > 0) {
 201         j = starta+npi;
 202         k = startr+npi;
 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;
 207         npi++;
 208       } else if (r < 0) {
 209         j = starta+np+nni;
 210         k = startr+np+nni;
 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;
 215         nni++;
 216       }
 217     }
 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]);
 223     }
 224     atoms->nr -= (nsa-1)*(np+nn);
 225
 226     /* Copy the new positions back */
 227     for(i=index[0]; i<atoms->nr; i++)
 228       copy_rvec(xt[i],x[i]);
 229     sfree(xt);
 230   }
 231 }
 232
 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)
 235 {
 236 #define TEMP_FILENM "temp.top"
 237   FILE *fpin,*fpout;
 238   char  buf[STRLEN],buf2[STRLEN],*temp,**mol_line=NULL;
 239   int  line,i,nsol,nmol_line,sol_line,nsol_last;
 240   gmx_bool bMolecules;
 241
 242   printf("\nProcessing topology\n");
 243   fpin = ffopen(topinout,"r");
 244   fpout= ffopen(TEMP_FILENM,"w");
 245
 246   line=0;
 247   bMolecules = FALSE;
 248   nmol_line = 0;
 249   sol_line  = -1;
 250   nsol_last = -1;
 251   while (fgets(buf, STRLEN, fpin)) {
 252     line++;
 253     strcpy(buf2,buf);
 254     if ((temp=strchr(buf2,'\n')) != NULL)
 255       temp[0]='\0';
 256     ltrim(buf2);
 257     if (buf2[0]=='[') {
 258       buf2[0]=' ';
 259       if ((temp=strchr(buf2,'\n')) != NULL)
 260         temp[0]='\0';
 261       rtrim(buf2);
 262       if (buf2[strlen(buf2)-1]==']') {
 263         buf2[strlen(buf2)-1]='\0';
 264         ltrim(buf2);
 265         rtrim(buf2);
 266         bMolecules=(gmx_strcasecmp(buf2,"molecules")==0);
 267       }
 268       fprintf(fpout,"%s",buf);
 269     } else if (!bMolecules) {
 270       fprintf(fpout,"%s",buf);
 271     } else {
 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);
 277       }
 278       /* Store this molecules section line */
 279       srenew(mol_line,nmol_line+1);
 280       mol_line[nmol_line] = strdup(buf);
 281       nmol_line++;
 282     }
 283   }
 284   ffclose(fpin);
 285
 286   if (sol_line == -1) {
 287     ffclose(fpout);
 288     gmx_fatal(FARGS,"No line with moleculetype '%s' found the [ molecules ] section of file '%s'",grpname,topinout);
 289   }
 290   if (nsol_last < p_num+n_num) {
 291     ffclose(fpout);
 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);
 293   }
 294
 295   /* Print all the molecule entries */
 296   for(i=0; i<nmol_line; i++) {
 297     if (i != sol_line) {
 298       fprintf(fpout,"%s",mol_line[i]);
 299     } else {
 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;
 304       if (nsol_last > 0) {
 305         fprintf(fpout,"%-10s  %d\n",grpname,nsol_last);
 306       }
 307       if (p_num > 0) {
 308         fprintf(fpout,"%-15s  %d\n",p_name,p_num);
 309       }
 310       if (n_num > 0) {
 311         fprintf(fpout,"%-15s  %d\n",n_name,n_num);
 312       }
 313     }
 314   }
 315   ffclose(fpout);
 316   /* use ffopen to generate backup of topinout */
 317   fpout = ffopen(topinout,"w");
 318   ffclose(fpout);
 319   rename(TEMP_FILENM,topinout);
 320 #undef TEMP_FILENM
 321 }
 322
 323 int gmx_genion(int argc, char *argv[])
 324 {
 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]."
 351   };
 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."
 355   };
 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." }
 376   };
 377   gmx_mtop_t  *mtop;
 378   gmx_localtop_t *top;
 379   t_inputrec  inputrec;
 380   t_commrec   *cr;
 381   t_mdatoms   *mdatoms;
 382   gmx_enerdata_t enerd;
 383   t_graph     *graph;
 384   t_forcerec  *fr;
 385   rvec        *x,*v;
 386   real        *pot,vol,qtot;
 387   matrix      box;
 388   t_atoms     atoms;
 389   t_pbc       pbc;
 390   int         *repl;
 391   atom_id     *index;
 392   char        *grpname;
 393   gmx_bool        *bSet,bPDB;
 394   int         i,nw,nwa,nsa,nsalt,iqtot;
 395   FILE        *fplog;
 396   output_env_t oenv;
 397   t_filenm fnm[] = {
 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 }
 405   };
 406 #define NFILE asize(fnm)
 407
 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");
 414     bPDB = FALSE;
 415   }
 416
 417   /* Check input for something sensible */
 418   if ((p_num<0) || (n_num<0))
 419     gmx_fatal(FARGS,"Negative number of ions to add?");
 420
 421   snew(mtop,1);
 422   snew(top,1);
 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);
 426
 427   atoms = gmx_mtop_global_atoms(mtop);
 428
 429   qtot = 0;
 430   for(i=0; (i<atoms.nr); i++)
 431     qtot += atoms.atom[i].q;
 432   iqtot = gmx_nint(qtot);
 433
 434   if ((conc > 0) || bNeutral) {
 435     /* Compute number of ions to be added */
 436     vol = det(box);
 437     if (conc > 0) {
 438       nsalt = gmx_nint(conc*vol*AVOGADRO/1e24);
 439       p_num = abs(nsalt*n_q);
 440       n_num = abs(nsalt*p_q);
 441       if (bNeutral) {
 442         int qdelta = 0;
 443         do {
 444           qdelta = (p_num*p_q + n_num*n_q + iqtot);
 445           if (qdelta < 0) {
 446             p_num  += abs(qdelta/p_q);
 447             qdelta = (p_num*p_q + n_num*n_q + iqtot);
 448           }
 449           if (qdelta > 0) {
 450             n_num  += abs(qdelta/n_q);
 451             qdelta = (p_num*p_q + n_num*n_q + iqtot);
 452           }
 453         } while (qdelta != 0);
 454       }
 455     }
 456   }
 457
 458   if ((p_num == 0) && (n_num == 0)) {
 459     if (!bPDB) {
 460       fprintf(stderr,"No ions to add and no potential to calculate.\n");
 461       exit(0);
 462     }
 463     nw  = 0;
 464     nsa = 0; /* to keep gcc happy */
 465   } else {
 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);
 470     for(i=1; i<nwa; i++)
 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);
 475     nsa = 1;
 476     while ((nsa<nwa) &&
 477            (atoms.atom[index[nsa]].resind ==
 478             atoms.atom[index[nsa-1]].resind))
 479       nsa++;
 480     if (nwa % nsa)
 481       gmx_fatal(FARGS,"Your solvent group size (%d) is not a multiple of %d",
 482                   nwa,nsa);
 483     nw = nwa/nsa;
 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");
 487   }
 488
 489   if (opt2bSet("-p",NFILE,fnm))
 490     update_topol(opt2fn("-p",NFILE,fnm),p_num,n_num,p_name,n_name,grpname);
 491
 492   snew(bSet,nw);
 493   snew(repl,nw);
 494
 495   snew(v,atoms.nr);
 496   snew(atoms.pdbinfo,atoms.nr);
 497
 498   set_pbc(&pbc,inputrec.ePBC,box);
 499
 500   /* Now loop over the ions that have to be placed */
 501   do {
 502     if (!bRandom) {
 503       calc_pot(fplog,cr,mtop,&inputrec,top,x,fr,&enerd,mdatoms,pot,box,graph);
 504       if (bPDB || debug) {
 505         char buf[STRLEN];
 506
 507         if (debug)
 508           sprintf(buf,"%d_%s",p_num+n_num,ftp2fn(efPDB,NFILE,fnm));
 509         else
 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);
 515         bPDB = FALSE;
 516       }
 517     }
 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);
 521       p_num--;
 522     }
 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);
 526       n_num--;
 527     }
 528   } while (p_num+n_num > 0);
 529   fprintf(stderr,"\n");
 530
 531   if (nw)
 532     sort_ions(nsa,nw,repl,index,&atoms,x,p_name,n_name);
 533
 534   sfree(atoms.pdbinfo);
 535   atoms.pdbinfo = NULL;
 536   write_sto_conf(ftp2fn(efSTO,NFILE,fnm),*mtop->name,&atoms,x,NULL,
 537                  inputrec.ePBC,box);
 538
 539   thanx(stderr);
 540
 541   gmx_log_close(fplog);
 542
 543   return 0;
 544 }