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