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