src/tools/gmx_genion.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  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.
  14
  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.
  19  *
  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.
  26  *
  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.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <ctype.h>
  40 #include "copyrite.h"
  41 #include "string2.h"
  42 #include "smalloc.h"
  43 #include "sysstuff.h"
  44 #include "confio.h"
  45 #include "statutil.h"
  46 #include "pbc.h"
  47 #include "force.h"
  48 #include "gmx_fatal.h"
  49 #include "futil.h"
  50 #include "maths.h"
  51 #include "macros.h"
  52 #include "physics.h"
  53 #include "vec.h"
  54 #include "tpxio.h"
  55 #include "mdrun.h"
  56 #include "calcpot.h"
  57 #include "main.h"
  58 #include "random.h"
  59 #include "index.h"
  60 #include "mtop_util.h"
  61 #include "gmx_ana.h"
  62
  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,
  67                        t_mdatoms *mdatoms,
  68                        real rmin,gmx_bool bRandom,int *seed)
  69 {
  70   int  i,ii,ei,owater,wlast,m,nw;
  71   real extr_e,poti,rmin2;
  72   rvec xei,dx;
  73   gmx_bool bSub=FALSE;
  74   gmx_large_int_t maxrand;
  75
  76   ei=-1;
  77   nw = *nwater;
  78   maxrand  = nw;
  79   maxrand *= 1000;
  80   if (bRandom) {
  81     do {
  82       ei = nw*rando(seed);
  83       maxrand--;
  84     } while (bSet[ei] && (maxrand > 0));
  85     if (bSet[ei])
  86       gmx_fatal(FARGS,"No more replaceable solvent!");
  87   }
  88   else {
  89     extr_e = 0;
  90     for(i=0; (i<nw); i++) {
  91       if (!bSet[i]) {
  92         ii=index[nsa*i];
  93         poti=pot[ii];
  94         if (q > 0) {
  95           if ((poti <= extr_e) || !bSub) {
  96             extr_e = poti;
  97             ei     = i;
  98             bSub   = TRUE;
  99           }
 100         }
 101         else {
 102           if ((poti >= extr_e) || !bSub) {
 103             extr_e = poti;
 104             ei     = i;
 105             bSub   = TRUE;
 106           }
 107         }
 108       }
 109     }
 110     if (ei == -1)
 111       gmx_fatal(FARGS,"No more replaceable solvent!");
 112   }
 113   fprintf(stderr,"Replacing solvent molecule %d (atom %d) with %s\n",
 114           ei,index[nsa*ei],ionname);
 115
 116   /* Replace solvent molecule charges with ion charge */
 117   bSet[ei] = TRUE;
 118   repl[ei] = sign;
 119   mdatoms->chargeA[index[nsa*ei]] = q;
 120   for(i=1; i<nsa; i++)
 121     mdatoms->chargeA[index[nsa*ei+i]] = 0;
 122
 123   /* Mark all solvent molecules within rmin as unavailable for substitution */
 124   if (rmin > 0) {
 125     rmin2=rmin*rmin;
 126     for(i=0; (i<nw); i++) {
 127       if (!bSet[i]) {
 128         pbc_dx(pbc,x[index[nsa*ei]],x[index[nsa*i]],dx);
 129         if (iprod(dx,dx) < rmin2)
 130           bSet[i] = TRUE;
 131       }
 132     }
 133   }
 134 }
 135
 136 static char *aname(const char *mname)
 137 {
 138   char *str;
 139   int  i;
 140
 141   str = strdup(mname);
 142   i=strlen(str)-1;
 143   while (i>1 && (isdigit(str[i]) || (str[i]=='+') || (str[i]=='-'))) {
 144     str[i]='\0';
 145     i--;
 146   }
 147
 148   return str;
 149 }
 150
 151 void sort_ions(int nsa,int nw,int repl[],atom_id index[],
 152                t_atoms *atoms,rvec x[],
 153                const char *p_name,const char *n_name)
 154 {
 155   int i,j,k,r,np,nn,starta,startr,npi,nni;
 156   rvec *xt;
 157   char **pptr=NULL,**nptr=NULL,**paptr=NULL,**naptr=NULL;
 158
 159   snew(xt,atoms->nr);
 160
 161   /* Put all the solvent in front and count the added ions */
 162   np=0;
 163   nn=0;
 164   j=index[0];
 165   for(i=0; i<nw; i++) {
 166     r = repl[i];
 167     if (r == 0)
 168       for(k=0; k<nsa; k++)
 169         copy_rvec(x[index[nsa*i+k]],xt[j++]);
 170     else if (r>0)
 171       np++;
 172     else if (r<0)
 173       nn++;
 174   }
 175
 176   if (np+nn > 0) {
 177     /* Put the positive and negative ions at the end */
 178     starta = index[nsa*(nw - np - nn)];
 179     startr = atoms->atom[starta].resind;
 180
 181     if (np) {
 182       snew(pptr,1);
 183       pptr[0] = strdup(p_name);
 184       snew(paptr,1);
 185       paptr[0] = aname(p_name);
 186     }
 187     if (nn) {
 188       snew(nptr,1);
 189       nptr[0] = strdup(n_name);
 190       snew(naptr,1);
 191       naptr[0] = aname(n_name);
 192     }
 193     npi = 0;
 194     nni = 0;
 195     for(i=0; i<nw; i++) {
 196       r = repl[i];
 197       if (r > 0) {
 198         j = starta+npi;
 199         k = startr+npi;
 200         copy_rvec(x[index[nsa*i]],xt[j]);
 201         atoms->atomname[j] = paptr;
 202         atoms->atom[j].resind = k ;
 203         atoms->resinfo[k].name = pptr;
 204         npi++;
 205       } else if (r < 0) {
 206         j = starta+np+nni;
 207         k = startr+np+nni;
 208         copy_rvec(x[index[nsa*i]],xt[j]);
 209         atoms->atomname[j] = naptr;
 210         atoms->atom[j].resind = k;
 211         atoms->resinfo[k].name = nptr;
 212         nni++;
 213       }
 214     }
 215     for(i=index[nsa*nw-1]+1; i<atoms->nr; i++) {
 216       j = i-(nsa-1)*(np+nn);
 217       atoms->atomname[j] = atoms->atomname[i];
 218       atoms->atom[j] = atoms->atom[i];
 219       copy_rvec(x[i],xt[j]);
 220     }
 221     atoms->nr -= (nsa-1)*(np+nn);
 222
 223     /* Copy the new positions back */
 224     for(i=index[0]; i<atoms->nr; i++)
 225       copy_rvec(xt[i],x[i]);
 226     sfree(xt);
 227   }
 228 }
 229
 230 static void update_topol(const char *topinout,int p_num,int n_num,
 231                          const char *p_name,const char *n_name,char *grpname)
 232 {
 233 #define TEMP_FILENM "temp.top"
 234   FILE *fpin,*fpout;
 235   char  buf[STRLEN],buf2[STRLEN],*temp,**mol_line=NULL;
 236   int  line,i,nsol,nmol_line,sol_line,nsol_last;
 237   gmx_bool bMolecules;
 238
 239   printf("\nProcessing topology\n");
 240   fpin = ffopen(topinout,"r");
 241   fpout= ffopen(TEMP_FILENM,"w");
 242
 243   line=0;
 244   bMolecules = FALSE;
 245   nmol_line = 0;
 246   sol_line  = -1;
 247   nsol_last = -1;
 248   while (fgets(buf, STRLEN, fpin)) {
 249     line++;
 250     strcpy(buf2,buf);
 251     if ((temp=strchr(buf2,'\n')) != NULL)
 252       temp[0]='\0';
 253     ltrim(buf2);
 254     if (buf2[0]=='[') {
 255       buf2[0]=' ';
 256       if ((temp=strchr(buf2,'\n')) != NULL)
 257         temp[0]='\0';
 258       rtrim(buf2);
 259       if (buf2[strlen(buf2)-1]==']') {
 260         buf2[strlen(buf2)-1]='\0';
 261         ltrim(buf2);
 262         rtrim(buf2);
 263         bMolecules=(gmx_strcasecmp(buf2,"molecules")==0);
 264       }
 265       fprintf(fpout,"%s",buf);
 266     } else if (!bMolecules) {
 267       fprintf(fpout,"%s",buf);
 268     } else {
 269       /* Check if this is a line with solvent molecules */
 270       sscanf(buf,"%s",buf2);
 271       if (gmx_strcasecmp(buf2,grpname) == 0) {
 272         sol_line = nmol_line;
 273         sscanf(buf,"%*s %d",&nsol_last);
 274       }
 275       /* Store this molecules section line */
 276       srenew(mol_line,nmol_line+1);
 277       mol_line[nmol_line] = strdup(buf);
 278       nmol_line++;
 279     }
 280   }
 281   ffclose(fpin);
 282
 283   if (sol_line == -1) {
 284     ffclose(fpout);
 285     gmx_fatal(FARGS,"No line with moleculetype '%s' found the [ molecules ] section of file '%s'",grpname,topinout);
 286   }
 287   if (nsol_last < p_num+n_num) {
 288     ffclose(fpout);
 289     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);
 290   }
 291
 292   /* Print all the molecule entries */
 293   for(i=0; i<nmol_line; i++) {
 294     if (i != sol_line) {
 295       fprintf(fpout,"%s",mol_line[i]);
 296     } else {
 297       printf("Replacing %d solute molecules in topology file (%s) "
 298              " by %d %s and %d %s ions.\n",
 299              p_num+n_num,topinout,p_num,p_name,n_num,n_name);
 300       nsol_last -= p_num + n_num;
 301       if (nsol_last > 0) {
 302         fprintf(fpout,"%-10s  %d\n",grpname,nsol_last);
 303       }
 304       if (p_num > 0) {
 305         fprintf(fpout,"%-10s  %d\n",p_name,p_num);
 306       }
 307       if (n_num > 0) {
 308         fprintf(fpout,"%-10s  %d\n",n_name,n_num);
 309       }
 310     }
 311   }
 312   ffclose(fpout);
 313   /* use ffopen to generate backup of topinout */
 314   fpout = ffopen(topinout,"w");
 315   ffclose(fpout);
 316   rename(TEMP_FILENM,topinout);
 317 #undef TEMP_FILENM
 318 }
 319
 320 int gmx_genion(int argc, char *argv[])
 321 {
 322   const char *desc[] = {
 323     "[TT]genion[tt] replaces solvent molecules by monoatomic ions at",
 324     "the position of the first atoms with the most favorable electrostatic",
 325     "potential or at random. The potential is calculated on all atoms, using",
 326     "normal GROMACS particle-based methods (in contrast to other methods",
 327     "based on solving the Poisson-Boltzmann equation).",
 328     "The potential is recalculated after every ion insertion.",
 329     "If specified in the run input file, a reaction field, shift function",
 330     "or user function can be used. For the user function a table file",
 331     "can be specified with the option [TT]-table[tt].",
 332     "The group of solvent molecules should be continuous and all molecules",
 333     "should have the same number of atoms.",
 334     "The user should add the ion molecules to the topology file or use",
 335     "the [TT]-p[tt] option to automatically modify the topology.[PAR]",
 336     "The ion molecule type, residue and atom names in all force fields",
 337     "are the capitalized element names without sign. This molecule name",
 338     "should be given with [TT]-pname[tt] or [TT]-nname[tt], and the",
 339     "[TT][molecules][tt] section of your topology updated accordingly,",
 340     "either by hand or with [TT]-p[tt]. Do not use an atom name instead!",
 341     "[PAR]Ions which can have multiple charge states get the multiplicity",
 342     "added, without sign, for the uncommon states only.[PAR]",
 343     "With the option [TT]-pot[tt] the potential can be written as B-factors",
 344     "in a [TT].pdb[tt] file (for visualisation using e.g. rasmol).",
 345     "The unit of the potential is 1000 kJ/(mol e), the scaling be changed",
 346     "with the [TT]-scale[tt] option.[PAR]",
 347     "For larger ions, e.g. sulfate we recommended using [TT]genbox[tt]."
 348   };
 349   const char *bugs[] = {
 350     "Calculation of the potential is not reliable, therefore the [TT]-random[tt] option is now turned on by default.",
 351     "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."
 352   };
 353   static int  p_num=0,n_num=0,p_q=1,n_q=-1;
 354   static const char *p_name="NA",*n_name="CL";
 355   static real rmin=0.6,scale=0.001,conc=0;
 356   static int  seed=1993;
 357   static gmx_bool bRandom=TRUE,bNeutral=FALSE;
 358   static t_pargs pa[] = {
 359     { "-np",    FALSE, etINT,  {&p_num}, "Number of positive ions"       },
 360     { "-pname", FALSE, etSTR,  {&p_name},"Name of the positive ion"      },
 361     { "-pq",    FALSE, etINT,  {&p_q},   "Charge of the positive ion"    },
 362     { "-nn",    FALSE, etINT,  {&n_num}, "Number of negative ions"       },
 363     { "-nname", FALSE, etSTR,  {&n_name},"Name of the negative ion"      },
 364     { "-nq",    FALSE, etINT,  {&n_q},   "Charge of the negative ion"    },
 365     { "-rmin",  FALSE, etREAL, {&rmin},  "Minimum distance between ions" },
 366     { "-random",FALSE,etBOOL, {&bRandom},"Use random placement of ions instead of based on potential. The rmin option should still work" },
 367     { "-seed",  FALSE, etINT,  {&seed},  "Seed for random number generator" },
 368     { "-scale", FALSE, etREAL, {&scale}, "Scaling factor for the potential for [TT]-pot[tt]" },
 369     { "-conc",  FALSE, etREAL, {&conc},
 370       "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." },
 371     { "-neutral", FALSE, etBOOL, {&bNeutral},
 372       "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." }
 373   };
 374   gmx_mtop_t  *mtop;
 375   gmx_localtop_t *top;
 376   t_inputrec  inputrec;
 377   t_commrec   *cr;
 378   t_mdatoms   *mdatoms;
 379   gmx_enerdata_t enerd;
 380   t_graph     *graph;
 381   t_forcerec  *fr;
 382   rvec        *x,*v;
 383   real        *pot,vol,qtot;
 384   matrix      box;
 385   t_atoms     atoms;
 386   t_pbc       pbc;
 387   int         *repl;
 388   atom_id     *index;
 389   char        *grpname;
 390   gmx_bool        *bSet,bPDB;
 391   int         i,nw,nwa,nsa,nsalt,iqtot;
 392   FILE        *fplog;
 393   output_env_t oenv;
 394   t_filenm fnm[] = {
 395     { efTPX, NULL,  NULL,      ffREAD  },
 396     { efXVG, "-table","table", ffOPTRD },
 397     { efNDX, NULL,  NULL,      ffOPTRD },
 398     { efSTO, "-o",  NULL,      ffWRITE },
 399     { efLOG, "-g",  "genion",  ffWRITE },
 400     { efPDB, "-pot", "pot",    ffOPTWR },
 401     { efTOP, "-p",  "topol",   ffOPTRW }
 402   };
 403 #define NFILE asize(fnm)
 404
 405   CopyRight(stderr,argv[0]);
 406   parse_common_args(&argc,argv,PCA_BE_NICE,NFILE,fnm,asize(pa),pa,
 407                     asize(desc),desc, asize(bugs),bugs,&oenv);
 408   bPDB = ftp2bSet(efPDB,NFILE,fnm);
 409   if (bRandom && bPDB) {
 410     fprintf(stderr,"Not computing potential with random option!\n");
 411     bPDB = FALSE;
 412   }
 413
 414   /* Check input for something sensible */
 415   if ((p_num<0) || (n_num<0))
 416     gmx_fatal(FARGS,"Negative number of ions to add?");
 417
 418   snew(mtop,1);
 419   snew(top,1);
 420   fplog = init_calcpot(ftp2fn(efLOG,NFILE,fnm),ftp2fn(efTPX,NFILE,fnm),
 421                        opt2fn("-table",NFILE,fnm),mtop,top,&inputrec,&cr,
 422                        &graph,&mdatoms,&fr,&enerd,&pot,box,&x,oenv);
 423
 424   atoms = gmx_mtop_global_atoms(mtop);
 425
 426   qtot = 0;
 427   for(i=0; (i<atoms.nr); i++)
 428     qtot += atoms.atom[i].q;
 429   iqtot = gmx_nint(qtot);
 430
 431   if ((conc > 0) || bNeutral) {
 432     /* Compute number of ions to be added */
 433     vol = det(box);
 434     if (conc > 0) {
 435       nsalt = gmx_nint(conc*vol*AVOGADRO/1e24);
 436       p_num = abs(nsalt*n_q);
 437       n_num = abs(nsalt*p_q);
 438       if (bNeutral) {
 439         int qdelta = 0;
 440         do {
 441           qdelta = (p_num*p_q + n_num*n_q + iqtot);
 442           if (qdelta < 0) {
 443             p_num  += abs(qdelta/p_q);
 444             qdelta = (p_num*p_q + n_num*n_q + iqtot);
 445           }
 446           if (qdelta > 0) {
 447             n_num  += abs(qdelta/n_q);
 448             qdelta = (p_num*p_q + n_num*n_q + iqtot);
 449           }
 450         } while (qdelta != 0);
 451       }
 452     }
 453   }
 454
 455   if ((p_num == 0) && (n_num == 0)) {
 456     if (!bPDB) {
 457       fprintf(stderr,"No ions to add and no potential to calculate.\n");
 458       exit(0);
 459     }
 460     nw  = 0;
 461     nsa = 0; /* to keep gcc happy */
 462   } else {
 463     printf("Will try to add %d %s ions and %d %s ions.\n",
 464            p_num,p_name,n_num,n_name);
 465     printf("Select a continuous group of solvent molecules\n");
 466     get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&nwa,&index,&grpname);
 467     for(i=1; i<nwa; i++)
 468       if (index[i] != index[i-1]+1)
 469         gmx_fatal(FARGS,"The solvent group %s is not continuous: "
 470                   "index[%d]=%d, index[%d]=%d",
 471                   grpname,i,index[i-1]+1,i+1,index[i]+1);
 472     nsa = 1;
 473     while ((nsa<nwa) &&
 474            (atoms.atom[index[nsa]].resind ==
 475             atoms.atom[index[nsa-1]].resind))
 476       nsa++;
 477     if (nwa % nsa)
 478       gmx_fatal(FARGS,"Your solvent group size (%d) is not a multiple of %d",
 479                   nwa,nsa);
 480     nw = nwa/nsa;
 481     fprintf(stderr,"Number of (%d-atomic) solvent molecules: %d\n",nsa,nw);
 482         if (p_num+n_num > nw)
 483       gmx_fatal(FARGS,"Not enough solvent for adding ions");
 484   }
 485
 486   if (opt2bSet("-p",NFILE,fnm))
 487     update_topol(opt2fn("-p",NFILE,fnm),p_num,n_num,p_name,n_name,grpname);
 488
 489   snew(bSet,nw);
 490   snew(repl,nw);
 491
 492   snew(v,atoms.nr);
 493   snew(atoms.pdbinfo,atoms.nr);
 494
 495   set_pbc(&pbc,inputrec.ePBC,box);
 496
 497   /* Now loop over the ions that have to be placed */
 498   do {
 499     if (!bRandom) {
 500       calc_pot(fplog,cr,mtop,&inputrec,top,x,fr,&enerd,mdatoms,pot,box,graph);
 501       if (bPDB || debug) {
 502         char buf[STRLEN];
 503
 504         if (debug)
 505           sprintf(buf,"%d_%s",p_num+n_num,ftp2fn(efPDB,NFILE,fnm));
 506         else
 507           strcpy(buf,ftp2fn(efPDB,NFILE,fnm));
 508         for(i=0; (i<atoms.nr); i++)
 509             atoms.pdbinfo[i].bfac = pot[i]*scale;
 510         write_sto_conf(buf,"Potential calculated by genion",
 511                        &atoms,x,v,inputrec.ePBC,box);
 512         bPDB = FALSE;
 513       }
 514     }
 515     if ((p_num > 0) && (p_num >= n_num))  {
 516       insert_ion(nsa,&nw,bSet,repl,index,pot,x,&pbc,
 517                  1,p_q,p_name,mdatoms,rmin,bRandom,&seed);
 518       p_num--;
 519     }
 520     else if (n_num > 0) {
 521       insert_ion(nsa,&nw,bSet,repl,index,pot,x,&pbc,
 522                  -1,n_q,n_name,mdatoms,rmin,bRandom,&seed);
 523       n_num--;
 524     }
 525   } while (p_num+n_num > 0);
 526   fprintf(stderr,"\n");
 527
 528   if (nw)
 529     sort_ions(nsa,nw,repl,index,&atoms,x,p_name,n_name);
 530
 531   sfree(atoms.pdbinfo);
 532   atoms.pdbinfo = NULL;
 533   write_sto_conf(ftp2fn(efSTO,NFILE,fnm),*mtop->name,&atoms,x,NULL,
 534                  inputrec.ePBC,box);
 535
 536   thanx(stderr);
 537
 538   gmx_log_close(fplog);
 539
 540   return 0;
 541 }