src/tools/genion.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_genion_c = "$Id$";
  30
  31 #include "copyrite.h"
  32 #include "string2.h"
  33 #include "smalloc.h"
  34 #include "sysstuff.h"
  35 #include "confio.h"
  36 #include "assert.h"
  37 #include "statutil.h"
  38 #include "pbc.h"
  39 #include "force.h"
  40 #include "fatal.h"
  41 #include "futil.h"
  42 #include "maths.h"
  43 #include "macros.h"
  44 #include "physics.h"
  45 #include "vec.h"
  46 #include "tpxio.h"
  47 #include "mdrun.h"
  48 #include "calcpot.h"
  49 #include "main.h"
  50 #include "random.h"
  51 #include "rdgroup.h"
  52
  53 static void insert_ion(int nsa,int *nwater,
  54                        bool bSet[],int repl[],atom_id index[],
  55                        real pot[],rvec x[],
  56                        int sign,real q,char *ionname,
  57                        t_mdatoms *mdatoms,
  58                        real rmin,bool bRandom,int *seed)
  59 {
  60   int  i,ii,ei,owater,wlast,m,nw;
  61   real extr_e,poti,rmin2;
  62   rvec xei,dx;
  63   bool bSub=FALSE;
  64   int  maxrand;
  65
  66   ei=-1;
  67   nw = *nwater;
  68   maxrand = 100*nw;
  69   if (bRandom) {
  70     do {
  71       ei = nw*rando(seed);
  72       maxrand--;
  73     } while (bSet[ei] && (maxrand > 0));
  74   }
  75   else {
  76     extr_e = 0;
  77     for(i=0; (i<nw); i++) {
  78       if (!bSet[i]) {
  79         ii=index[nsa*i];
  80         poti=pot[nsa*ii];
  81         if (q > 0) {
  82           if ((poti <= extr_e) || !bSub) {
  83             extr_e = poti;
  84             ei     = i;
  85             bSub   = TRUE;
  86           }
  87         }
  88         else {
  89           if ((poti >= extr_e) || !bSub) {
  90             extr_e = poti;
  91             ei     = i;
  92             bSub   = TRUE;
  93           }
  94         }
  95       }
  96     }
  97   }
  98   if (ei == -1)
  99     fatal_error(0,"No more replaceable solvent!");
 100   fprintf(stderr,"Replacing solvent molecule %d (atom %d) with %s\n",
 101           ei,index[nsa*ei],ionname);
 102
 103   /* Replace solvent molecule charges with ion charge */
 104   bSet[ei] = TRUE;
 105   repl[ei] = sign;
 106   mdatoms->chargeA[index[nsa*ei]] = q;
 107   for(i=1; i<nsa; i++)
 108     mdatoms->chargeA[index[nsa*ei+i]] = 0;
 109
 110   /* Mark all solvent molecules within rmin as unavailable for substitution */
 111   if (rmin > 0) {
 112     rmin2=rmin*rmin;
 113     for(i=0; (i<nw); i++) {
 114       if (!bSet[i]) {
 115         pbc_dx(x[index[nsa*ei]],x[index[nsa*i]],dx);
 116         if (iprod(dx,dx) < rmin2)
 117           bSet[i] = TRUE;
 118       }
 119     }
 120   }
 121 }
 122
 123 static char *aname(char *mname)
 124 {
 125   char *str;
 126   int  i;
 127
 128   str = strdup(mname);
 129   i=strlen(str)-1;
 130   while (i>1 && (isdigit(str[i]) || (str[i]=='+') || (str[i]=='-'))) {
 131     str[i]='\0';
 132     i--;
 133   }
 134
 135   return str;
 136 }
 137
 138 void sort_ions(int nsa,int nw,int repl[],atom_id index[],
 139                t_atoms *atoms,rvec x[],
 140                char *p_name,char *n_name)
 141 {
 142   int i,j,k,r,np,nn,starta,startr,npi,nni;
 143   rvec *xt;
 144   char **pptr=NULL,**nptr=NULL,**paptr=NULL,**naptr=NULL;
 145
 146   snew(xt,atoms->nr);
 147
 148   /* Put all the solvent in front and count the added ions */
 149   np=0;
 150   nn=0;
 151   j=index[0];
 152   for(i=0; i<nw; i++) {
 153     r = repl[i];
 154     if (r == 0)
 155       for(k=0; k<nsa; k++)
 156         copy_rvec(x[index[nsa*i+k]],xt[j++]);
 157     else if (r>0)
 158       np++;
 159     else if (r<0)
 160       nn++;
 161   }
 162
 163   if (np+nn > 0) {
 164     /* Put the positive and negative ions at the end */
 165     starta = index[nsa*(nw - np - nn)];
 166     startr = atoms->atom[starta].resnr;
 167
 168     if (np) {
 169       snew(pptr,1);
 170       pptr[0] = p_name;
 171       snew(paptr,1);
 172       paptr[0] = aname(p_name);
 173     }
 174     if (nn) {
 175       snew(nptr,1);
 176       nptr[0] = n_name;
 177       snew(naptr,1);
 178       naptr[0] = aname(n_name);
 179     }
 180     npi = 0;
 181     nni = 0;
 182     for(i=0; i<nw; i++) {
 183       r = repl[i];
 184       if (r > 0) {
 185         j = starta+npi;
 186         k = startr+npi;
 187         copy_rvec(x[index[nsa*i]],xt[j]);
 188         atoms->atomname[j] = paptr;
 189         atoms->atom[j].resnr = k ;
 190         atoms->resname[k] = pptr;
 191         npi++;
 192       } else if (r < 0) {
 193         j = starta+np+nni;
 194         k = startr+np+nni;
 195         copy_rvec(x[index[nsa*i]],xt[j]);
 196         atoms->atomname[j] = naptr;
 197         atoms->atom[j].resnr = k;
 198         atoms->resname[k] = nptr;
 199         nni++;
 200       }
 201     }
 202     for(i=index[nsa*nw-1]+1; i<atoms->nr; i++) {
 203       j = i-(nsa-1)*(np+nn);
 204       atoms->atomname[j] = atoms->atomname[i];
 205       atoms->atom[j] = atoms->atom[i];
 206       copy_rvec(x[i],xt[j]);
 207     }
 208     atoms->nr -= (nsa-1)*(np+nn);
 209
 210     /* Copy the new positions back */
 211     for(i=index[0]; i<atoms->nr; i++)
 212       copy_rvec(xt[i],x[i]);
 213     sfree(xt);
 214   }
 215 }
 216
 217 int main(int argc, char *argv[])
 218 {
 219   static char *desc[] = {
 220     "genion replaces solvent molecules by monoatomic ions at",
 221     "the position of the first atoms with the most favorable electrostatic",
 222     "potential or at random. The potential is calculated on all atoms, using",
 223     "normal GROMACS particle based methods (in contrast to other methods",
 224     "based on solving the Poisson-Boltzmann equation).",
 225     "The potential is recalculated after every ion insertion.",
 226     "If specified in the run input file, a reaction field or shift function",
 227     "can be used.",
 228     "The group of solvent molecules should be continuous and all molecules",
 229     "should have the same number of atoms.",
 230     "The user should add the ion molecules to the topology file and include",
 231     "the file [TT]ions.itp[tt].",
 232     "Ion names for Gromos96 should include the charge.[PAR]",
 233     "The potential can be written as B-factors",
 234     "in a pdb file (for visualisation using e.g. rasmol)[PAR]"
 235     "For larger ions, e.g. sulfate we recommended to use genbox."
 236   };
 237   static int  p_num=0,n_num=0;
 238   static char *p_name="Na",*n_name="Cl";
 239   static real p_q=1,n_q=-1,rmin=0.6;
 240   static int  seed=1993;
 241   static bool bRandom=FALSE;
 242   static t_pargs pa[] = {
 243     { "-np",    FALSE, etINT,  {&p_num}, "Number of positive ions"       },
 244     { "-pname", FALSE, etSTR,  {&p_name},"Name of the positive ion"      },
 245     { "-pq",    FALSE, etREAL, {&p_q},   "Charge of the positive ion"    },
 246     { "-nn",    FALSE, etINT,  {&n_num}, "Number of negative ions"       },
 247     { "-nname", FALSE, etSTR,  {&n_name},"Name of the negative ion"      },
 248     { "-nq",    FALSE, etREAL, {&n_q},   "Charge of the negative ion"    },
 249     { "-rmin",  FALSE, etREAL, {&rmin},  "Minimum distance between ions" },
 250     { "-random",FALSE,etBOOL, {&bRandom},"Use random placement of ions instead of based on potential. The rmin option should still work" },
 251     { "-seed",  FALSE, etINT,  {&seed},  "Seed for random number generator" }
 252   };
 253   t_topology  *top;
 254   t_parm      parm;
 255   t_commrec   cr;
 256   t_mdatoms   *mdatoms;
 257   t_nsborder  nsb;
 258   t_groups    grps;
 259   t_graph     *graph;
 260   t_forcerec  *fr;
 261   rvec        *x,*v;
 262   real        *pot;
 263   matrix      box;
 264   int         *repl;
 265   atom_id     *index;
 266   char        *grpname;
 267   bool        *bSet,bPDB;
 268   int         i,nw,nwa,nsa;
 269   t_filenm fnm[] = {
 270     { efTPX, NULL,  NULL,      ffREAD  },
 271     { efNDX, NULL,  NULL,      ffOPTRD },
 272     { efSTO, "-o",  NULL,      ffWRITE },
 273     { efLOG, "-g",  "genion",  ffWRITE },
 274     { efPDB, "-pot", "pot",    ffOPTWR }
 275   };
 276 #define NFILE asize(fnm)
 277
 278   CopyRight(stdout,argv[0]);
 279   parse_common_args(&argc,argv,0,TRUE,NFILE,fnm,asize(pa),pa,asize(desc),desc,
 280                     0,NULL);
 281   bPDB = ftp2bSet(efPDB,NFILE,fnm);
 282   if (bRandom && bPDB) {
 283     fprintf(stderr,"Not computing potential with random option!\n");
 284     bPDB = FALSE;
 285   }
 286
 287   /* Check input for something sensible */
 288   if ((p_num<0) || (n_num<0))
 289     fatal_error(0,"Negative number of ions to add?");
 290
 291   snew(top,1);
 292   init_calcpot(NFILE,fnm,top,&x,&parm,&cr,
 293                &graph,&mdatoms,&nsb,&grps,&fr,&pot,box);
 294
 295   if ((p_num == 0) && (n_num == 0)) {
 296     if (!bPDB) {
 297       fprintf(stderr,"No ions to add and no potential to calculate.\n");
 298       exit(0);
 299     }
 300     nw = 0;
 301   } else {
 302     printf("Select a continuous group of solvent molecules\n");
 303     get_index(&top->atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&nwa,&index,&grpname);
 304     for(i=1; i<nwa; i++)
 305       if (index[i] != index[i-1]+1)
 306         fatal_error(0,"The solvent group is not continuous: index[%d]=%d, "
 307                     "index[%d]=%d",i,index[i-1]+1,i+1,index[i]+1);
 308     nsa = 1;
 309     while ((nsa<nwa) &&
 310            (top->atoms.atom[index[nsa]].resnr ==
 311             top->atoms.atom[index[nsa-1]].resnr))
 312       nsa++;
 313     if (nwa % nsa)
 314       fatal_error(0,"Your solvent group size (%d) is not a multiple of %d",
 315                   nwa,nsa);
 316     nw = nwa/nsa;
 317     fprintf(stderr,"Number of (%d-atomic) solvent molecules: %d\n",nsa,nw);
 318         if (p_num+n_num > nw)
 319       fatal_error(0,"Not enough solvent for adding ions");
 320   }
 321   snew(bSet,nw);
 322   snew(repl,nw);
 323
 324   snew(v,top->atoms.nr);
 325   snew(top->atoms.pdbinfo,top->atoms.nr);
 326
 327   /* Now loop over the ions that have to be placed */
 328   do {
 329     if (!bRandom) {
 330       calc_pot(stdlog,&nsb,&cr,&grps,&parm,top,x,fr,mdatoms,pot);
 331       if (bPDB || debug) {
 332         char buf[STRLEN];
 333
 334         if (debug)
 335           sprintf(buf,"%d_%s",p_num+n_num,ftp2fn(efPDB,NFILE,fnm));
 336         else
 337           strcpy(buf,ftp2fn(efPDB,NFILE,fnm));
 338         for(i=0; (i<top->atoms.nr); i++)
 339             top->atoms.pdbinfo[i].bfac = pot[i]*0.001;
 340         write_sto_conf(buf,"Potential calculated by genion",
 341                        &top->atoms,x,v,box);
 342         bPDB = FALSE;
 343       }
 344     }
 345     if ((p_num > 0) && (p_num >= n_num))  {
 346       insert_ion(nsa,&nw,bSet,repl,index,pot,x,
 347                  1,p_q,p_name,mdatoms,rmin,bRandom,&seed);
 348       p_num--;
 349     }
 350     else if (n_num > 0) {
 351       insert_ion(nsa,&nw,bSet,repl,index,pot,x,
 352                  -1,n_q,n_name,mdatoms,rmin,bRandom,&seed);
 353       n_num--;
 354     }
 355   } while (p_num+n_num > 0);
 356   fprintf(stderr,"\n");
 357
 358   if (nw)
 359     sort_ions(nsa,nw,repl,index,&top->atoms,x,p_name,n_name);
 360
 361   sfree(top->atoms.pdbinfo);
 362   top->atoms.pdbinfo = NULL;
 363   write_sto_conf(ftp2fn(efSTO,NFILE,fnm),*top->name,&top->atoms,x,NULL,box);
 364
 365   thanx(stdout);
 366
 367   return 0;
 368 }