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