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