3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
48 #include "gmx_fatal.h"
52 #include "gmx_statistics.h"
57 #define even(a) ( ( (a+1) / 2) == (a / 2) )
60 static void make_dist_leg(FILE *fp, int gnx, atom_id index[], t_atoms *atoms,
61 const output_env_t oenv)
67 for (i = 0; i < gnx; i += 2)
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);
78 xvgr_legend(fp, gnx/2, (const char**)leg, oenv);
79 for (i = 0; i < gnx/2; i++)
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)
93 FILE *out, *outd = NULL;
95 real b0 = 0, b1, db = 0;
97 gmx_stats_t b_one = NULL, *b_all = NULL;
98 /*real mean, mean2, sqrdev2, sigma2;
106 int bind, i, nframes, i0, i1;
109 real aver, sigma, error;
114 for (i = 0; (i < gnx/2); i++)
116 b_all[i] = gmx_stats_init();
121 b_one = gmx_stats_init();
122 snew(btab, MAXTAB+1);
125 natoms = read_first_x(oenv, &status, fn, &t, &x, box);
128 gmx_fatal(FARGS, "No atoms in trajectory!");
133 outd = xvgropen(fdist, bAverDist ? "Average distance" : "Distances",
134 "Time (ps)", "Distance (nm)", oenv);
137 make_dist_leg(outd, gnx, index, &(top->atoms), oenv);
144 set_pbc(&pbc, ePBC, box);
147 fprintf(outd, " %8.4f", t);
149 nframes++; /* count frames */
151 for (i = 0; (i < gnx); i += 2)
153 pbc_dx(&pbc, x[index[i]], x[index[i+1]], dx);
161 fprintf(outd, " %.3f", bond);
165 gmx_stats_add_point(b_one, t, bond, 0, 0);
178 db = (2.0*(b1-b0))/MAXTAB;
181 bind = (int)((bond-b0)/db+0.5);
182 if ((bind >= 0) && (bind <= MAXTAB))
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]);
196 gmx_stats_add_point(b_all[i/2], t, bond, 0, 0);
201 fprintf(outd, " %.5f", bav*2.0/gnx);
208 while (read_next_x(oenv, status, &t, natoms, x, box));
217 mean = mean / counter;
218 mean2 = mean2 / counter;
219 sqrdev2 = (mean2 - mean*mean);
220 sigma2 = sqrdev2*counter / (counter - 1);
222 /* For definitions see "Weet wat je meet" */
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);
235 out = xvgropen(fbond, "Bond Stretching Distribution",
236 "Bond Length (nm)", "", oenv);
238 for (i0 = 0; ((i0 < MAXTAB) && (btab[i0] == 0)); i0++)
243 for (i1 = MAXTAB; ((i1 > 0) && (btab[i1] == 0)); i1--)
247 i1 = min(MAXTAB, i1+1);
251 gmx_fatal(FARGS, "No distribution... (i0 = %d, i1 = %d)? ? ! ! ? !", i0, i1);
254 fac = 2.0/(nframes*gnx*db);
255 for (i = i0; (i <= i1); i++)
257 fprintf(out, "%8.5f %8.5f\n", b0+i*db, btab[i]*fac);
263 fprintf(log, "%5s %5s %8s %8s\n", "i", "j", "b_aver", "sigma");
264 for (i = 0; (i < gnx/2); i++)
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],
269 gmx_stats_done(b_all[i]);
276 int gmx_bond(int argc, char *argv[])
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."
291 const char *bugs[] = {
292 "It should be possible to get bond information from the topology."
294 static real blen = -1.0, tol = 0.1;
295 static gmx_bool bAver = TRUE, bAverDist = TRUE;
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])" }
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 }
326 #define NFILE asize(fnm)
328 parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
329 NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs,
334 fdist = opt2fn("-d", NFILE, fnm);
338 fdist = opt2fn_null("-d", NFILE, fnm);
341 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &x, NULL, box,
346 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &gnx, &index, &grpname);
349 fprintf(stderr, "WARNING: odd number of atoms (%d) in group!\n", gnx);
351 fprintf(stderr, "Will gather information on %d bonds\n", gnx/2);
355 fp = ftp2FILE(efLOG, NFILE, fnm, "w");
362 do_bonds(fp, ftp2fn(efTRX, NFILE, fnm), opt2fn("-o", NFILE, fnm), fdist, gnx, index,
363 blen, tol, bAver, &top, ePBC, bAverDist, oenv);
365 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
366 do_view(oenv, opt2fn_null("-d", NFILE, fnm), "-nxy");