src/tools/gmx_bond.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 #include "sysstuff.h"
  41 #include "typedefs.h"
  42 #include "smalloc.h"
  43 #include "macros.h"
  44 #include "vec.h"
  45 #include "pbc.h"
  46 #include "xvgr.h"
  47 #include "copyrite.h"
  48 #include "gmx_fatal.h"
  49 #include "futil.h"
  50 #include "statutil.h"
  51 #include "index.h"
  52 #include "gmx_statistics.h"
  53 #include "tpxio.h"
  54
  55 static void make_dist_leg(FILE *fp,int gnx,atom_id index[],t_atoms *atoms,
  56                           const output_env_t oenv)
  57 {
  58   char **leg;
  59   int  i;
  60
  61   snew(leg,gnx/2);
  62   for(i=0; i<gnx; i+=2) {
  63     snew(leg[i/2],256);
  64     sprintf(leg[i/2],"%s %s%d - %s %s%d",
  65             *(atoms->atomname[index[i]]),
  66             *(atoms->resinfo[atoms->atom[index[i]].resind].name),
  67             atoms->resinfo[atoms->atom[index[i]].resind].nr,
  68             *(atoms->atomname[index[i+1]]),
  69             *(atoms->resinfo[atoms->atom[index[i+1]].resind].name),
  70             atoms->resinfo[atoms->atom[index[i+1]].resind].nr);
  71   }
  72   xvgr_legend(fp,gnx/2,leg,oenv);
  73   for(i=0; i<gnx/2; i++)
  74     sfree(leg[i]);
  75   sfree(leg);
  76 }
  77
  78 static void do_bonds(FILE *log,const char *fn,const char *fbond,
  79                      const char *fdist, int gnx,atom_id index[],
  80                      real blen,real tol,bool bAver,
  81                      t_topology *top,int ePBC,bool bAverDist,
  82                      const output_env_t oenv)
  83 {
  84 #define MAXTAB 1000
  85   FILE   *out,*outd=NULL;
  86   int    *btab=NULL;
  87   real   b0=0,b1,db=0;
  88   real   bond,bav;
  89   gmx_stats_t b_one=NULL,*b_all=NULL;
  90   /*real   mean, mean2, sqrdev2, sigma2;
  91     int    counter;*/
  92   rvec   *x;
  93   rvec   dx;
  94   int    status,natoms;
  95   matrix box;
  96   real   t,fac;
  97   int    bind,i,nframes,i0,i1;
  98   t_pbc  pbc;
  99   int    N;
 100   real   aver,sigma,error;
 101
 102   if (!bAver) {
 103     snew(b_all,gnx/2);
 104     for(i=0; (i<gnx/2); i++)
 105       b_all[i] = gmx_stats_init();
 106   }
 107   else {
 108     b_one = gmx_stats_init();
 109     snew(btab,MAXTAB+1);
 110   }
 111
 112   natoms=read_first_x(oenv,&status,fn,&t,&x,box);
 113   if (natoms == 0)
 114     gmx_fatal(FARGS,"No atoms in trajectory!");
 115
 116   if (fdist) {
 117     outd = xvgropen(fdist,bAverDist ? "Average distance" : "Distances",
 118                     "Time (ps)","Distance (nm)",oenv);
 119     if (!bAverDist)
 120       make_dist_leg(outd,gnx,index,&(top->atoms),oenv);
 121   }
 122
 123   nframes=0;
 124   do {
 125     set_pbc(&pbc,ePBC,box);
 126     if (fdist)
 127       fprintf(outd," %8.4f",t);
 128     nframes++; /* count frames */
 129     bav = 0.0;
 130     for(i=0; (i<gnx); i+=2) {
 131       pbc_dx(&pbc,x[index[i]],x[index[i+1]],dx);
 132       bond   = norm(dx);
 133       if (bAverDist)
 134         bav += bond;
 135       else if (fdist)
 136         fprintf(outd," %.3f",bond);
 137       if (bAver) {
 138         gmx_stats_add_point(b_one,t,bond,0,0);
 139         if (db == 0) {
 140           if (blen == -1) {
 141             b0 = 0;
 142             b1 = 0.2;
 143             db = (b1-b0)/MAXTAB;
 144           }
 145           else {
 146             b0   = (1.0-tol)*blen;
 147             b1   = (1.0+tol)*blen;
 148             db   = (2.0*(b1-b0))/MAXTAB;
 149           }
 150         }
 151         bind = (int)((bond-b0)/db+0.5);
 152         if ((bind >= 0) && (bind <= MAXTAB))
 153           btab[bind]++;
 154         else {
 155           /*
 156             printf("bond: %4d-%4d bond=%10.5e, dx=(%10.5e,%10.5e,%10.5e)\n",
 157             index[i],index[i+1],bond,dx[XX],dx[YY],dx[ZZ]);
 158           */
 159         }
 160       }
 161       else {
 162         gmx_stats_add_point(b_all[i/2],t,bond,0,0);
 163       }
 164     }
 165     if (bAverDist)
 166       fprintf(outd," %.5f",bav*2.0/gnx);
 167     if (fdist)
 168       fprintf(outd,"\n");
 169   } while (read_next_x(oenv,status,&t,natoms,x,box));
 170   close_trj(status);
 171
 172   if (fdist)
 173     ffclose(outd);
 174
 175   /*
 176     mean = mean / counter;
 177     mean2 = mean2 / counter;
 178     sqrdev2 = (mean2 - mean*mean);
 179     sigma2 = sqrdev2*counter / (counter - 1);
 180   */
 181   /* For definitions see "Weet wat je meet" */
 182   if (bAver) {
 183     printf("\n");
 184     gmx_stats_get_npoints(b_one,&N);
 185     printf("Total number of samples               : %d\n",N);
 186     gmx_stats_get_ase(b_one,&aver,&sigma,&error);
 187     printf("Mean                                  : %g\n",aver);
 188     printf("Standard deviation of the distribution: %g\n",sigma);
 189     printf("Standard deviation of the mean        : %g\n",error);
 190     gmx_stats_done(b_one);
 191     sfree(b_one);
 192
 193     out=xvgropen(fbond,"Bond Stretching Distribution",
 194                  "Bond Length (nm)","",oenv);
 195
 196     for(i0=0;      ((i0 < MAXTAB) && (btab[i0]==0)); i0++)
 197       ;
 198     i0=max(0,i0-1);
 199     for(i1=MAXTAB; ((i1 > 0)      && (btab[i1]==0)); i1--)
 200       ;
 201     i1=min(MAXTAB,i1+1);
 202
 203     if (i0 >= i1)
 204       gmx_fatal(FARGS,"No distribution... (i0 = %d, i1 = %d)? ? ! ! ? !",i0,i1);
 205
 206     fac=2.0/(nframes*gnx*db);
 207     for(i=i0; (i<=i1); i++)
 208       fprintf(out,"%8.5f  %8.5f\n",b0+i*db,btab[i]*fac);
 209     fclose(out);
 210   }
 211   else {
 212     fprintf(log,"%5s  %5s  %8s  %8s\n","i","j","b_aver","sigma");
 213     for(i=0; (i<gnx/2); i++) {
 214       gmx_stats_get_ase(b_all[i],&aver,&sigma,NULL);
 215       fprintf(log,"%5u  %5u  %8.5f  %8.5f\n",1+index[2*i],1+index[2*i+1],
 216               aver,sigma);
 217       gmx_stats_done(b_all[i]);
 218       sfree(b_all[i]);
 219     }
 220     sfree(b_all);
 221   }
 222 }
 223
 224 int gmx_bond(int argc,char *argv[])
 225 {
 226   const char *desc[] = {
 227     "g_bond makes a distribution of bond lengths. If all is well a",
 228     "gaussian distribution should be made when using a harmonic potential.",
 229     "bonds are read from a single group in the index file in order i1-j1",
 230     "i2-j2 thru in-jn.[PAR]",
 231     "[TT]-tol[tt] gives the half-width of the distribution as a fraction",
 232     "of the bondlength ([TT]-blen[tt]). That means, for a bond of 0.2",
 233     "a tol of 0.1 gives a distribution from 0.18 to 0.22.[PAR]",
 234     "Option [TT]-d[tt] plots all the distances as a function of time.",
 235     "This requires a structure file for the atom and residue names in",
 236     "the output. If however the option [TT]-averdist[tt] is given (as well",
 237     "or separately) the average bond length is plotted instead."
 238   };
 239   const char *bugs[] = {
 240     "It should be possible to get bond information from the topology."
 241   };
 242   static real blen=-1.0,tol=0.1;
 243   static bool bAver=TRUE,bAverDist=TRUE;
 244   t_pargs pa[] = {
 245     { "-blen", FALSE, etREAL, {&blen},
 246       "Bond length. By default length of first bond" },
 247     { "-tol",  FALSE, etREAL, {&tol},
 248       "Half width of distribution as fraction of blen" },
 249     { "-aver", FALSE, etBOOL, {&bAver},
 250       "Average bond length distributions" },
 251     { "-averdist", FALSE, etBOOL, {&bAverDist},
 252       "Average distances (turns on -d)" }
 253   };
 254   FILE      *fp;
 255   char      *grpname;
 256   const char *fdist;
 257   int       gnx;
 258   atom_id   *index;
 259   char      title[STRLEN];
 260   t_topology top;
 261   int       ePBC=-1;
 262   rvec      *x;
 263   matrix    box;
 264   output_env_t oenv;
 265
 266   t_filenm fnm[] = {
 267     { efTRX, "-f", NULL, ffREAD  },
 268     { efNDX, NULL, NULL, ffREAD  },
 269     { efTPS, NULL, NULL, ffOPTRD },
 270     { efXVG, "-o", "bonds", ffWRITE },
 271     { efLOG, NULL, "bonds", ffOPTWR },
 272     { efXVG, "-d", "distance", ffOPTWR }
 273   };
 274 #define NFILE asize(fnm)
 275
 276   CopyRight(stderr,argv[0]);
 277   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE ,
 278                     NFILE,fnm,asize(pa),pa,asize(desc),desc,asize(bugs),bugs,
 279                     &oenv);
 280
 281   if (bAverDist)
 282     fdist = opt2fn("-d",NFILE,fnm);
 283   else {
 284     fdist = opt2fn_null("-d",NFILE,fnm);
 285     if (fdist)
 286       read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&x,NULL,box,
 287                     FALSE);
 288   }
 289
 290   rd_index(ftp2fn(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);
 291   if ( !even(gnx) )
 292     fprintf(stderr,"WARNING: odd number of atoms (%d) in group!\n",gnx);
 293   fprintf(stderr,"Will gather information on %d bonds\n",gnx/2);
 294
 295   if (!bAver)
 296     fp = ftp2FILE(efLOG,NFILE,fnm,"w");
 297   else
 298     fp = NULL;
 299
 300   do_bonds(fp,ftp2fn(efTRX,NFILE,fnm),opt2fn("-o",NFILE,fnm),fdist,gnx,index,
 301            blen,tol,bAver,&top,ePBC,bAverDist,oenv);
 302
 303   do_view(oenv,opt2fn("-o",NFILE,fnm),"-nxy");
 304   do_view(oenv,opt2fn_null("-d",NFILE,fnm),"-nxy");
 305
 306   thanx(stderr);
 307
 308   return 0;
 309 }
 310