src/tools/g_bond.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  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_g_bond_c = "$Id$";
  30
  31 #include <math.h>
  32 #include <string.h>
  33 #include "sysstuff.h"
  34 #include "typedefs.h"
  35 #include "smalloc.h"
  36 #include "macros.h"
  37 #include "vec.h"
  38 #include "pbc.h"
  39 #include "xvgr.h"
  40 #include "copyrite.h"
  41 #include "fatal.h"
  42 #include "futil.h"
  43 #include "statutil.h"
  44 #include "rdgroup.h"
  45 #include "gstat.h"
  46
  47 void do_bonds(FILE *log,char *fn,char *outf,int gnx,atom_id index[],
  48               real blen,real tol,bool bAver)
  49 {
  50 #define MAXTAB 1000
  51   FILE   *out;
  52   int    *btab=NULL;
  53   real   b0=0,b1,db=0;
  54   real   bond;
  55   t_lsq  b_one,*b_all=NULL;
  56   /*real   mean, mean2, sqrdev2, sigma2;
  57     int    counter;*/
  58   rvec   *x;
  59   rvec   dx;
  60   int    status,natoms;
  61   matrix box;
  62   real   t,fac;
  63   int    bind,i,j,i0,i1;
  64
  65   if (!bAver) {
  66     snew(b_all,gnx/2);
  67     for(i=0; (i<gnx/2); i++)
  68       init_lsq(&(b_all[i]));
  69   }
  70   else {
  71     init_lsq(&b_one);
  72     snew(btab,MAXTAB+1);
  73   }
  74
  75   natoms=read_first_x(&status,fn,&t,&x,box);
  76   init_pbc(box,FALSE);
  77   if (natoms == 0)
  78     fatal_error(0,"No atoms in trajectory!");
  79
  80   j=0;
  81   do {
  82     j++; /* count frames */
  83     for(i=0; (i<gnx); i+=2) {
  84       pbc_dx(x[index[i]],x[index[i+1]],dx);
  85       bond   = norm(dx);
  86       if (bAver) {
  87         add_lsq(&b_one,t,bond);
  88         if (db == 0) {
  89           if (blen == -1) {
  90             b0 = 0;
  91             b1 = 0.2;
  92             db = (b1-b0)/MAXTAB;
  93           }
  94           else {
  95             b0   = (1.0-tol)*blen;
  96             b1   = (1.0+tol)*blen;
  97             db   = (2.0*(b1-b0))/MAXTAB;
  98           }
  99         }
 100         bind = ((bond-b0)/db);
 101         if ((bind >= 0) && (bind <= MAXTAB))
 102           btab[bind]++;
 103         else {
 104           /*
 105             printf("bond: %4d-%4d bond=%10.5e, dx=(%10.5e,%10.5e,%10.5e)\n",
 106             index[i],index[i+1],bond,dx[XX],dx[YY],dx[ZZ]);
 107           */
 108         }
 109       }
 110       else {
 111         add_lsq(&(b_all[i/2]),t,bond);
 112       }
 113     }
 114   } while (read_next_x(status,&t,natoms,x,box));
 115   close_trj(status);
 116
 117   /*
 118     mean = mean / counter;
 119     mean2 = mean2 / counter;
 120     sqrdev2 = (mean2 - mean*mean);
 121     sigma2 = sqrdev2*counter / (counter - 1);
 122   */
 123   /* For definitions see "Weet wat je meet" */
 124   if (bAver) {
 125     fprintf(stderr,"\n");
 126     fprintf(stderr,"Total number of samples               : %d\n",
 127             npoints_lsq(&b_one));
 128     fprintf(stderr,"Mean                                  : %g\n",
 129             aver_lsq(&b_one));
 130     fprintf(stderr,"Standard deviation of the distribution: %g\n",
 131             sigma_lsq(&b_one));
 132     fprintf(stderr,"Standard deviation of the mean        : %g\n",
 133             error_lsq(&b_one));
 134
 135     out=xvgropen(outf,"Bond Stretching Distribution",
 136                  "Bond Length (nm)","Time (ps)");
 137
 138     for(i0=0;      ((i0 < MAXTAB) && (btab[i0]==0)); i0++)
 139       ;
 140     i0=max(0,i0-1);
 141     for(i1=MAXTAB; ((i1 > 0)      && (btab[i1]==0)); i1--)
 142       ;
 143     i1=min(MAXTAB,i1+1);
 144
 145     if (i0 >= i1)
 146       fatal_error(0,"No distribution... (i0 = %d, i1 = %d)? ? ! ! ? !",i0,i1);
 147
 148     fac=((2.0/j)/gnx);
 149     for(i=i0; (i<=i1); i++)
 150       fprintf(out,"%16.10e  %16.10e  %10d\n",b0+i*db,btab[i]*fac,btab[i]);
 151     fclose(out);
 152   }
 153   else {
 154     fprintf(log,"%5s  %5s  %8s  %8s\n","i","j","b_aver","sigma");
 155     for(i=0; (i<gnx/2); i++) {
 156       fprintf(log,"%5u  %5u  %8.5f  %8.5f\n",1+index[2*i],1+index[2*i+1],
 157               aver_lsq(&b_all[i]),sigma_lsq(&b_all[i]));
 158     }
 159   }
 160 }
 161
 162 int main(int argc,char *argv[])
 163 {
 164   static char *desc[] = {
 165     "g_bond makes a distribution of bond lengths. If all is well a",
 166     "gaussian distribution should be made when using a harmonic potential.",
 167     "bonds are read from a single group in the index file in order i1-j1",
 168     "i2-j2 thru in-jn.[PAR]",
 169     "[TT]-tol[tt] gives the half-width of the distribution as a fraction",
 170     "of the bondlength ([TT]-blen[tt]). That means, for a bond of 0.2",
 171     "a tol of 0.1 gives a distribution from 0.18 to 0.22"
 172   };
 173   static char *bugs[] = {
 174     "It should be possible to get bond information from the topology."
 175   };
 176   static real blen=-1.0,tol=0.1;
 177   static bool bAver=TRUE;
 178   t_pargs pa[] = {
 179     { "-blen", FALSE, etREAL, {&blen},
 180       "Bond length. By default length of first bond" },
 181     { "-tol",  FALSE, etREAL, {&tol},
 182       "Half width of distribution as fraction of blen" },
 183     { "-aver", FALSE, etBOOL, {&bAver},
 184       "Sum up distributions" }
 185   };
 186   FILE      *fp;
 187   char      *grpname;
 188   int       gnx;
 189   atom_id   *index;
 190   t_filenm fnm[] = {
 191     { efTRX, "-f", NULL, ffREAD  },
 192     { efNDX, NULL, NULL, ffREAD  },
 193     { efXVG, NULL, "bonds", ffWRITE },
 194     { efLOG, NULL, "bonds", ffOPTWR }
 195   };
 196 #define NFILE asize(fnm)
 197
 198   CopyRight(stderr,argv[0]);
 199   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME,TRUE,
 200                     NFILE,fnm,asize(pa),pa,asize(desc),desc,asize(bugs),bugs);
 201
 202   rd_index(ftp2fn(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);
 203   if ( !even(gnx) )
 204     fprintf(stderr,"WARNING: odd number of atoms (%d) in group!\n",gnx);
 205   fprintf(stderr,"Will gather information on %d bonds\n",gnx/2);
 206
 207   if (!bAver)
 208     fp = ftp2FILE(efLOG,NFILE,fnm,"w");
 209   else
 210     fp = NULL;
 211
 212   do_bonds(fp,ftp2fn(efTRX,NFILE,fnm),ftp2fn(efXVG,NFILE,fnm),gnx,index,
 213            blen,tol,bAver);
 214
 215   xvgr_file(ftp2fn(efXVG,NFILE,fnm),NULL);
 216
 217   thanx(stdout);
 218
 219   return 0;
 220 }