src/tools/gmx_saltbr.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, 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 #include <math.h>
  42 #include <string.h>
  43
  44 #include "macros.h"
  45 #include "vec.h"
  46 #include "sysstuff.h"
  47 #include "typedefs.h"
  48 #include "filenm.h"
  49 #include "statutil.h"
  50 #include "copyrite.h"
  51 #include "futil.h"
  52 #include "gmx_fatal.h"
  53 #include "smalloc.h"
  54 #include "pbc.h"
  55 #include "xvgr.h"
  56 #include "gmx_ana.h"
  57
  58
  59 typedef struct {
  60   char *label;
  61   int  cg;
  62   real q;
  63 } t_charge;
  64
  65 static t_charge *mk_charge(t_atoms *atoms,t_block *cgs,int *nncg)
  66 {
  67   t_charge *cg=NULL;
  68   char     buf[32];
  69   int      i,j,ncg,resnr,anr;
  70   real     qq;
  71
  72   /* Find the charged groups */
  73   ncg=0;
  74   for(i=0; (i<cgs->nr); i++) {
  75     qq=0.0;
  76     for(j=cgs->index[i]; (j<cgs->index[i+1]); j++) {
  77       qq+=atoms->atom[j].q;
  78     }
  79     if (fabs(qq) > 1.0e-5) {
  80       srenew(cg,ncg+1);
  81       cg[ncg].q=qq;
  82       cg[ncg].cg=i;
  83       anr=cgs->index[i];
  84       resnr=atoms->atom[anr].resind;
  85       sprintf(buf,"%s%d-%d",
  86               *(atoms->resinfo[resnr].name),
  87               atoms->resinfo[resnr].nr,
  88               anr+1);
  89       cg[ncg].label=strdup(buf);
  90       ncg++;
  91     }
  92   }
  93   *nncg=ncg;
  94
  95   for(i=0; (i<ncg); i++) {
  96     printf("CG: %10s Q: %6g  Atoms:",
  97            cg[i].label,cg[i].q);
  98     for(j=cgs->index[cg[i].cg]; (j<cgs->index[cg[i].cg+1]); j++)
  99       printf(" %4u",j);
 100     printf("\n");
 101   }
 102
 103   return cg;
 104 }
 105
 106 static real calc_dist(t_pbc *pbc,rvec x[],t_block *cgs,int icg,int jcg)
 107 {
 108   int  i,j;
 109   rvec dx;
 110   real d2,mindist2=1000;
 111
 112   for(i=cgs->index[icg]; (i<cgs->index[icg+1]); i++) {
 113     for(j=cgs->index[jcg]; (j<cgs->index[jcg+1]); j++) {
 114       pbc_dx(pbc,x[i],x[j],dx);
 115       d2 = norm2(dx);
 116       if (d2 < mindist2)
 117         mindist2 = d2;
 118     }
 119   }
 120   return sqrt(mindist2);
 121 }
 122
 123 int gmx_saltbr(int argc,char *argv[])
 124 {
 125   const char *desc[] = {
 126     "[TT]g_saltbr[tt] plots the distance between all combination of charged groups",
 127     "as a function of time. The groups are combined in different ways.",
 128     "A minimum distance can be given (i.e. a cut-off), such that groups",
 129     "that are never closer than that distance will not be plotted.[PAR]",
 130     "Output will be in a number of fixed filenames, [TT]min-min.xvg[tt], [TT]plus-min.xvg[tt]",
 131     "and [TT]plus-plus.xvg[tt], or files for every individual ion pair if the [TT]-sep[tt]",
 132     "option is selected. In this case, files are named as [TT]sb-(Resname)(Resnr)-(Atomnr)[tt].",
 133     "There may be [BB]many[bb] such files."
 134   };
 135   static gmx_bool bSep=FALSE;
 136   static real truncate=1000.0;
 137   t_pargs pa[] = {
 138     { "-t",   FALSE, etREAL, {&truncate},
 139       "Groups that are never closer than this distance are not plotted" },
 140     { "-sep", FALSE, etBOOL, {&bSep},
 141       "Use separate files for each interaction (may be MANY)" }
 142   };
 143   t_filenm   fnm[] = {
 144     { efTRX, "-f",  NULL, ffREAD },
 145     { efTPX, NULL,  NULL, ffREAD },
 146   };
 147 #define NFILE asize(fnm)
 148
 149   FILE       *out[3],*fp;
 150   static const char *title[3] = {
 151     "Distance between positively charged groups",
 152     "Distance between negatively charged groups",
 153     "Distance between oppositely charged groups"
 154   };
 155   static const char *fn[3] = {
 156     "plus-plus.xvg",
 157     "min-min.xvg",
 158     "plus-min.xvg"
 159   };
 160   int        nset[3]={0,0,0};
 161
 162   t_topology *top;
 163   int        ePBC;
 164   char       *buf;
 165   t_trxstatus *status;
 166   int        i,j,k,m,nnn,teller,ncg,n1,n2,n3,natoms;
 167   real       t,*time,qi,qj;
 168   t_charge   *cg;
 169   real       ***cgdist;
 170   int        **nWithin;
 171
 172   double     t0,dt;
 173   char       label[234];
 174   t_pbc      pbc;
 175   rvec       *x;
 176   matrix     box;
 177   output_env_t oenv;
 178
 179   CopyRight(stderr,argv[0]);
 180
 181   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,
 182                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 183
 184   top=read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);
 185   cg=mk_charge(&top->atoms,&(top->cgs),&ncg);
 186   snew(cgdist,ncg);
 187   snew(nWithin,ncg);
 188   for(i=0; (i<ncg); i++) {
 189     snew(cgdist[i],ncg);
 190     snew(nWithin[i],ncg);
 191   }
 192
 193   natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 194
 195   teller=0;
 196   time=NULL;
 197   do {
 198     srenew(time,teller+1);
 199     time[teller]=t;
 200
 201     set_pbc(&pbc,ePBC,box);
 202
 203     for(i=0; (i<ncg); i++) {
 204       for(j=i+1; (j<ncg); j++) {
 205         srenew(cgdist[i][j],teller+1);
 206         cgdist[i][j][teller]=
 207           calc_dist(&pbc,x,&(top->cgs),cg[i].cg,cg[j].cg);
 208         if (cgdist[i][j][teller] < truncate)
 209           nWithin[i][j]=1;
 210       }
 211     }
 212
 213     teller++;
 214   } while (read_next_x(oenv,status,&t,natoms,x,box));
 215   fprintf(stderr,"\n");
 216   close_trj(status);
 217
 218   if (bSep) {
 219     snew(buf,256);
 220     for(i=0; (i<ncg); i++)
 221       for(j=i+1; (j<ncg); j++) {
 222         if (nWithin[i][j]) {
 223           sprintf(buf,"sb-%s:%s.xvg",cg[i].label,cg[j].label);
 224           fp=xvgropen(buf,buf,"Time (ps)","Distance (nm)",oenv);
 225           for(k=0; (k<teller); k++)
 226             fprintf(fp,"%10g  %10g\n",time[k],cgdist[i][j][k]);
 227           ffclose(fp);
 228         }
 229       }
 230     sfree(buf);
 231   }
 232   else {
 233
 234     for(m=0; (m<3); m++)
 235       out[m]=xvgropen(fn[m],title[m],"Time (ps)","Distance (nm)",oenv);
 236
 237     snew(buf,256);
 238     for(i=0; (i<ncg); i++) {
 239       qi=cg[i].q;
 240       for(j=i+1; (j<ncg); j++) {
 241         qj=cg[j].q;
 242         if (nWithin[i][j]) {
 243           sprintf(buf,"%s:%s",cg[i].label,cg[j].label);
 244           if (qi*qj < 0)
 245             nnn=2;
 246           else if (qi+qj > 0)
 247             nnn=0;
 248           else
 249             nnn=1;
 250
 251           if (nset[nnn] == 0)
 252             xvgr_legend(out[nnn],1,(const char**)&buf,oenv);
 253           else {
 254             if (output_env_get_xvg_format(oenv) == exvgXMGR) {
 255               fprintf(out[nnn],"@ legend string %d \"%s\"\n",nset[nnn],buf);
 256             } else if (output_env_get_xvg_format(oenv) == exvgXMGRACE) {
 257               fprintf(out[nnn],"@ s%d legend \"%s\"\n",nset[nnn],buf);
 258             }
 259           }
 260           nset[nnn]++;
 261           nWithin[i][j]=nnn+1;
 262         }
 263       }
 264     }
 265     for(k=0; (k<teller); k++) {
 266       for(m=0; (m<3); m++)
 267         fprintf(out[m],"%10g",time[k]);
 268
 269       for(i=0; (i<ncg); i++) {
 270         for(j=i+1; (j<ncg); j++) {
 271           nnn=nWithin[i][j];
 272           if (nnn >0)
 273             fprintf(out[nnn-1],"  %10g",cgdist[i][j][k]);
 274         }
 275       }
 276       for(m=0; (m<3); m++)
 277         fprintf(out[m],"\n");
 278     }
 279     for(m=0; (m<3); m++) {
 280       ffclose(out[m]);
 281       if (nset[m] == 0)
 282         remove(fn[m]);
 283     }
 284   }
 285   thanx(stderr);
 286
 287   return 0;
 288 }