2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
45 #include "gromacs/legacyheaders/typedefs.h"
46 #include "gromacs/legacyheaders/macros.h"
47 #include "gromacs/pbcutil/mshift.h"
48 #include "gromacs/legacyheaders/viewit.h"
49 #include "gromacs/fileio/confio.h"
50 #include "gromacs/legacyheaders/nrnb.h"
51 #include "gromacs/legacyheaders/disre.h"
52 #include "gromacs/legacyheaders/force.h"
54 #include "gromacs/legacyheaders/main.h"
55 #include "gromacs/fileio/pdbio.h"
56 #include "gromacs/topology/index.h"
57 #include "gromacs/legacyheaders/mdatoms.h"
58 #include "gromacs/fileio/tpxio.h"
59 #include "gromacs/fileio/trxio.h"
60 #include "gromacs/legacyheaders/mdrun.h"
61 #include "gromacs/legacyheaders/names.h"
62 #include "gromacs/topology/mtop_util.h"
65 #include "gromacs/commandline/pargs.h"
66 #include "gromacs/fileio/matio.h"
67 #include "gromacs/fileio/xvgr.h"
68 #include "gromacs/math/do_fit.h"
69 #include "gromacs/math/vec.h"
70 #include "gromacs/pbcutil/ishift.h"
71 #include "gromacs/pbcutil/pbc.h"
72 #include "gromacs/pbcutil/rmpbc.h"
73 #include "gromacs/utility/fatalerror.h"
74 #include "gromacs/utility/smalloc.h"
86 real sumv, averv, maxv;
87 real *aver1, *aver2, *aver_3, *aver_6;
90 static void init5(int n)
96 static void reset5(void)
100 for (i = 0; (i < ntop); i++)
107 int tpcomp(const void *a, const void *b)
115 return (1e7*(tpb->v-tpa->v));
118 static void add5(int ndr, real viol)
123 for (i = 1; (i < ntop); i++)
125 if (top[i].v < top[mini].v)
130 if (viol > top[mini].v)
137 static void print5(FILE *fp)
141 qsort(top, ntop, sizeof(top[0]), tpcomp);
142 fprintf(fp, "Index:");
143 for (i = 0; (i < ntop); i++)
145 fprintf(fp, " %6d", top[i].n);
147 fprintf(fp, "\nViol: ");
148 for (i = 0; (i < ntop); i++)
150 fprintf(fp, " %6.3f", top[i].v);
155 static void check_viol(FILE *log,
156 t_ilist *disres, t_iparams forceparams[],
158 t_pbc *pbc, t_graph *g, t_dr_result dr[],
159 int clust_id, int isize, atom_id index[], real vvindex[],
163 int i, j, nat, n, type, nviol, ndr, label;
164 real ener, rt, mviol, tviol, viol, lam, dvdl, drt;
166 static gmx_bool bFirst = TRUE;
178 forceatoms = disres->iatoms;
179 for (j = 0; (j < isize); j++)
183 nat = interaction_function[F_DISRES].nratoms+1;
184 for (i = 0; (i < disres->nr); )
186 type = forceatoms[i];
188 label = forceparams[type].disres.label;
191 fprintf(debug, "DISRE: ndr = %d, label = %d i=%d, n =%d\n",
196 gmx_fatal(FARGS, "tpr inconsistency. ndr = %d, label = %d\n", ndr, label);
202 while (((i+n) < disres->nr) &&
203 (forceparams[forceatoms[i+n]].disres.label == label));
205 calc_disres_R_6(n, &forceatoms[i], forceparams,
206 (const rvec*)x, pbc, fcd, NULL);
208 if (fcd->disres.Rt_6[0] <= 0)
210 gmx_fatal(FARGS, "ndr = %d, rt_6 = %f", ndr, fcd->disres.Rt_6[0]);
213 rt = pow(fcd->disres.Rt_6[0], -1.0/6.0);
214 dr[clust_id].aver1[ndr] += rt;
215 dr[clust_id].aver2[ndr] += sqr(rt);
217 dr[clust_id].aver_3[ndr] += drt;
218 dr[clust_id].aver_6[ndr] += fcd->disres.Rt_6[0];
220 snew(fshift, SHIFTS);
221 ener = interaction_function[F_DISRES].ifunc(n, &forceatoms[i], forceparams,
222 (const rvec*)x, f, fshift,
223 pbc, g, lam, &dvdl, NULL, fcd, NULL);
225 viol = fcd->disres.sumviol;
232 add5(forceparams[type].disres.label, viol);
239 for (j = 0; (j < isize); j++)
241 if (index[j] == forceparams[type].disres.label)
243 vvindex[j] = pow(fcd->disres.Rt_6[0], -1.0/6.0);
250 dr[clust_id].nv = nviol;
251 dr[clust_id].maxv = mviol;
252 dr[clust_id].sumv = tviol;
253 dr[clust_id].averv = tviol/ndr;
254 dr[clust_id].nframes++;
258 fprintf(stderr, "\nThere are %d restraints and %d pairs\n",
259 ndr, disres->nr/nat);
271 real up1, r, rT3, rT6, viol, violT3, violT6;
274 static int drs_comp(const void *a, const void *b)
278 da = (t_dr_stats *)a;
279 db = (t_dr_stats *)b;
281 if (da->viol > db->viol)
285 else if (da->viol < db->viol)
295 static void dump_dump(FILE *log, int ndr, t_dr_stats drs[])
297 static const char *core[] = { "All restraints", "Core restraints" };
298 static const char *tp[] = { "linear", "third power", "sixth power" };
299 real viol_tot, viol_max, viol = 0;
304 for (bCore = FALSE; (bCore <= TRUE); bCore++)
306 for (kkk = 0; (kkk < 3); kkk++)
312 for (i = 0; (i < ndr); i++)
314 if (!bCore || (bCore && drs[i].bCore))
322 viol = drs[i].violT3;
325 viol = drs[i].violT6;
328 gmx_incons("Dumping violations");
330 viol_max = max(viol_max, viol);
339 if ((nrestr > 0) || (bCore && (nrestr < ndr)))
342 fprintf(log, "+++++++ %s ++++++++\n", core[bCore]);
343 fprintf(log, "+++++++ Using %s averaging: ++++++++\n", tp[kkk]);
344 fprintf(log, "Sum of violations: %8.3f nm\n", viol_tot);
347 fprintf(log, "Average violation: %8.3f nm\n", viol_tot/nrestr);
349 fprintf(log, "Largest violation: %8.3f nm\n", viol_max);
350 fprintf(log, "Number of violated restraints: %d/%d\n", nviol, nrestr);
356 static void dump_viol(FILE *log, int ndr, t_dr_stats *drs, gmx_bool bLinear)
360 fprintf(log, "Restr. Core Up1 <r> <rT3> <rT6> <viol><violT3><violT6>\n");
361 for (i = 0; (i < ndr); i++)
363 if (bLinear && (drs[i].viol == 0))
367 fprintf(log, "%6d%5s%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n",
368 drs[i].index, yesno_names[drs[i].bCore],
369 drs[i].up1, drs[i].r, drs[i].rT3, drs[i].rT6,
370 drs[i].viol, drs[i].violT3, drs[i].violT6);
374 static gmx_bool is_core(int i, int isize, atom_id index[])
377 gmx_bool bIC = FALSE;
379 for (kk = 0; !bIC && (kk < isize); kk++)
381 bIC = (index[kk] == i);
387 static void dump_stats(FILE *log, int nsteps, int ndr, t_ilist *disres,
388 t_iparams ip[], t_dr_result *dr,
389 int isize, atom_id index[], t_atoms *atoms)
395 fprintf(log, "++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
398 nra = interaction_function[F_DISRES].nratoms+1;
399 for (i = 0; (i < ndr); i++)
401 /* Search for the appropriate restraint */
402 for (; (j < disres->nr); j += nra)
404 if (ip[disres->iatoms[j]].disres.label == i)
410 drs[i].bCore = is_core(i, isize, index);
411 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
412 drs[i].r = dr->aver1[i]/nsteps;
413 drs[i].rT3 = pow(dr->aver_3[i]/nsteps, -1.0/3.0);
414 drs[i].rT6 = pow(dr->aver_6[i]/nsteps, -1.0/6.0);
415 drs[i].viol = max(0, drs[i].r-drs[i].up1);
416 drs[i].violT3 = max(0, drs[i].rT3-drs[i].up1);
417 drs[i].violT6 = max(0, drs[i].rT6-drs[i].up1);
420 int j1 = disres->iatoms[j+1];
421 int j2 = disres->iatoms[j+2];
422 atoms->pdbinfo[j1].bfac += drs[i].violT3*5;
423 atoms->pdbinfo[j2].bfac += drs[i].violT3*5;
426 dump_viol(log, ndr, drs, FALSE);
428 fprintf(log, "+++ Sorted by linear averaged violations: +++\n");
429 qsort(drs, ndr, sizeof(drs[0]), drs_comp);
430 dump_viol(log, ndr, drs, TRUE);
432 dump_dump(log, ndr, drs);
437 static void dump_clust_stats(FILE *fp, int ndr, t_ilist *disres,
438 t_iparams ip[], t_blocka *clust, t_dr_result dr[],
439 char *clust_name[], int isize, atom_id index[])
441 int i, j, k, nra, mmm = 0;
442 double sumV, maxV, sumVT3, sumVT6, maxVT3, maxVT6;
446 fprintf(fp, "++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
447 fprintf(fp, "Cluster NFrames SumV MaxV SumVT MaxVT SumVS MaxVS\n");
451 for (k = 0; (k < clust->nr); k++)
453 if (dr[k].nframes == 0)
457 if (dr[k].nframes != (clust->index[k+1]-clust->index[k]))
459 gmx_fatal(FARGS, "Inconsistency in cluster %s.\n"
460 "Found %d frames in trajectory rather than the expected %d\n",
461 clust_name[k], dr[k].nframes,
462 clust->index[k+1]-clust->index[k]);
466 gmx_fatal(FARGS, "Inconsistency with cluster %d. Invalid name", k);
469 nra = interaction_function[F_DISRES].nratoms+1;
470 sumV = sumVT3 = sumVT6 = maxV = maxVT3 = maxVT6 = 0;
471 for (i = 0; (i < ndr); i++)
473 /* Search for the appropriate restraint */
474 for (; (j < disres->nr); j += nra)
476 if (ip[disres->iatoms[j]].disres.label == i)
482 drs[i].bCore = is_core(i, isize, index);
483 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
484 drs[i].r = dr[k].aver1[i]/dr[k].nframes;
485 if ((dr[k].aver_3[i] <= 0) || (dr[k].aver_3[i] != dr[k].aver_3[i]))
487 gmx_fatal(FARGS, "dr[%d].aver_3[%d] = %f", k, i, dr[k].aver_3[i]);
489 drs[i].rT3 = pow(dr[k].aver_3[i]/dr[k].nframes, -1.0/3.0);
490 drs[i].rT6 = pow(dr[k].aver_6[i]/dr[k].nframes, -1.0/6.0);
491 drs[i].viol = max(0, drs[i].r-drs[i].up1);
492 drs[i].violT3 = max(0, drs[i].rT3-drs[i].up1);
493 drs[i].violT6 = max(0, drs[i].rT6-drs[i].up1);
495 sumVT3 += drs[i].violT3;
496 sumVT6 += drs[i].violT6;
497 maxV = max(maxV, drs[i].viol);
498 maxVT3 = max(maxVT3, drs[i].violT3);
499 maxVT6 = max(maxVT6, drs[i].violT6);
501 if (strcmp(clust_name[k], "1000") == 0)
505 fprintf(fp, "%-10s%6d%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n",
507 dr[k].nframes, sumV, maxV, sumVT3, maxVT3, sumVT6, maxVT6);
514 static void init_dr_res(t_dr_result *dr, int ndr)
516 snew(dr->aver1, ndr+1);
517 snew(dr->aver2, ndr+1);
518 snew(dr->aver_3, ndr+1);
519 snew(dr->aver_6, ndr+1);
527 static void dump_disre_matrix(const char *fn, t_dr_result *dr, int ndr,
528 int nsteps, t_idef *idef, gmx_mtop_t *mtop,
529 real max_dr, int nlevels, gmx_bool bThird)
533 int n_res, a_offset, mb, mol, a;
535 int iii, i, j, nra, nratoms, tp, ri, rj, index, nlabel, label;
536 atom_id ai, aj, *ptr;
537 real **matrix, *t_res, hi, *w_dr, rav, rviol;
538 t_rgb rlo = { 1, 1, 1 };
539 t_rgb rhi = { 0, 0, 0 };
545 snew(resnr, mtop->natoms);
548 for (mb = 0; mb < mtop->nmolblock; mb++)
550 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
551 for (mol = 0; mol < mtop->molblock[mb].nmol; mol++)
553 for (a = 0; a < atoms->nr; a++)
555 resnr[a_offset+a] = n_res + atoms->atom[a].resind;
557 n_res += atoms->nres;
558 a_offset += atoms->nr;
563 for (i = 0; (i < n_res); i++)
568 for (i = 0; (i < n_res); i++)
570 snew(matrix[i], n_res);
572 nratoms = interaction_function[F_DISRES].nratoms;
573 nra = (idef->il[F_DISRES].nr/(nratoms+1));
579 for (i = 0; (i < idef->il[F_DISRES].nr); i += nratoms+1)
581 tp = idef->il[F_DISRES].iatoms[i];
582 label = idef->iparams[tp].disres.label;
586 /* Set index pointer */
590 gmx_fatal(FARGS, "nlabel is %d, label = %d", nlabel, label);
594 gmx_fatal(FARGS, "ndr = %d, index = %d", ndr, index);
596 /* Update the weight */
597 w_dr[index] = 1.0/nlabel;
606 printf("nlabel = %d, index = %d, ndr = %d\n", nlabel, index, ndr);
608 for (i = 0; (i < ndr); i++)
610 for (j = ptr[i]; (j < ptr[i+1]); j += nratoms+1)
612 tp = idef->il[F_DISRES].iatoms[j];
613 ai = idef->il[F_DISRES].iatoms[j+1];
614 aj = idef->il[F_DISRES].iatoms[j+2];
620 rav = pow(dr->aver_3[i]/nsteps, -1.0/3.0);
624 rav = dr->aver1[i]/nsteps;
628 fprintf(debug, "DR %d, atoms %d, %d, distance %g\n", i, ai, aj, rav);
630 rviol = max(0, rav-idef->iparams[tp].disres.up1);
631 matrix[ri][rj] += w_dr[i]*rviol;
632 matrix[rj][ri] += w_dr[i]*rviol;
633 hi = max(hi, matrix[ri][rj]);
634 hi = max(hi, matrix[rj][ri]);
644 printf("Warning: the maxdr that you have specified (%g) is smaller than\nthe largest value in your simulation (%g)\n", max_dr, hi);
648 printf("Highest level in the matrix will be %g\n", hi);
649 fp = gmx_ffopen(fn, "w");
650 write_xpm(fp, 0, "Distance Violations", "<V> (nm)", "Residue", "Residue",
651 n_res, n_res, t_res, t_res, matrix, 0, hi, rlo, rhi, &nlevels);
655 int gmx_disre(int argc, char *argv[])
657 const char *desc[] = {
658 "[THISMODULE] computes violations of distance restraints.",
659 "If necessary, all protons can be added to a protein molecule ",
660 "using the [gmx-protonate] program.[PAR]",
661 "The program always",
662 "computes the instantaneous violations rather than time-averaged,",
663 "because this analysis is done from a trajectory file afterwards",
664 "it does not make sense to use time averaging. However,",
665 "the time averaged values per restraint are given in the log file.[PAR]",
666 "An index file may be used to select specific restraints for",
668 "When the optional [TT]-q[tt] flag is given a [TT].pdb[tt] file coloured by the",
669 "amount of average violations.[PAR]",
670 "When the [TT]-c[tt] option is given, an index file will be read",
671 "containing the frames in your trajectory corresponding to the clusters",
672 "(defined in another manner) that you want to analyze. For these clusters",
673 "the program will compute average violations using the third power",
674 "averaging algorithm and print them in the log file."
677 static int nlevels = 20;
678 static real max_dr = 0;
679 static gmx_bool bThird = TRUE;
681 { "-ntop", FALSE, etINT, {&ntop},
682 "Number of large violations that are stored in the log file every step" },
683 { "-maxdr", FALSE, etREAL, {&max_dr},
684 "Maximum distance violation in matrix output. If less than or equal to 0 the maximum will be determined by the data." },
685 { "-nlevels", FALSE, etINT, {&nlevels},
686 "Number of levels in the matrix output" },
687 { "-third", FALSE, etBOOL, {&bThird},
688 "Use inverse third power averaging or linear for matrix output" }
691 FILE *out = NULL, *aver = NULL, *numv = NULL, *maxxv = NULL, *xvg = NULL;
697 t_atoms *atoms = NULL;
701 int ntopatoms, natoms, i, j, kkk;
704 rvec *x, *f, *xav = NULL;
708 atom_id *index = NULL, *ind_fit = NULL;
710 t_cluster_ndx *clust = NULL;
711 t_dr_result dr, *dr_clust = NULL;
713 real *vvindex = NULL, *w_rls = NULL;
715 t_pbc pbc, *pbc_null;
719 gmx_rmpbc_t gpbc = NULL;
722 { efTPX, NULL, NULL, ffREAD },
723 { efTRX, "-f", NULL, ffREAD },
724 { efXVG, "-ds", "drsum", ffWRITE },
725 { efXVG, "-da", "draver", ffWRITE },
726 { efXVG, "-dn", "drnum", ffWRITE },
727 { efXVG, "-dm", "drmax", ffWRITE },
728 { efXVG, "-dr", "restr", ffWRITE },
729 { efLOG, "-l", "disres", ffWRITE },
730 { efNDX, NULL, "viol", ffOPTRD },
731 { efPDB, "-q", "viol", ffOPTWR },
732 { efNDX, "-c", "clust", ffOPTRD },
733 { efXPM, "-x", "matrix", ffOPTWR }
735 #define NFILE asize(fnm)
737 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
738 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
743 fplog = ftp2FILE(efLOG, NFILE, fnm, "w");
750 read_tpxheader(ftp2fn(efTPX, NFILE, fnm), &header, FALSE, NULL, NULL);
751 snew(xtop, header.natoms);
752 read_tpx(ftp2fn(efTPX, NFILE, fnm), &ir, box, &ntopatoms, xtop, NULL, NULL, &mtop);
753 bPDB = opt2bSet("-q", NFILE, fnm);
756 snew(xav, ntopatoms);
757 snew(ind_fit, ntopatoms);
758 snew(w_rls, ntopatoms);
759 for (kkk = 0; (kkk < ntopatoms); kkk++)
766 *atoms = gmx_mtop_global_atoms(&mtop);
768 if (atoms->pdbinfo == NULL)
770 snew(atoms->pdbinfo, atoms->nr);
774 top = gmx_mtop_generate_local_top(&mtop, &ir);
778 if (ir.ePBC != epbcNONE)
780 if (ir.bPeriodicMols)
786 g = mk_graph(fplog, &top->idef, 0, mtop.natoms, FALSE, FALSE);
790 if (ftp2bSet(efNDX, NFILE, fnm))
792 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
793 xvg = xvgropen(opt2fn("-dr", NFILE, fnm), "Individual Restraints", "Time (ps)",
795 snew(vvindex, isize);
797 for (i = 0; (i < isize); i++)
801 sprintf(leg[i], "index %d", index[i]);
803 xvgr_legend(xvg, isize, (const char**)leg, oenv);
811 init_disres(fplog, &mtop, &ir, NULL, &fcd, NULL, FALSE);
813 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
816 init_dr_res(&dr, fcd.disres.nres);
817 if (opt2bSet("-c", NFILE, fnm))
819 clust = cluster_index(fplog, opt2fn("-c", NFILE, fnm));
820 snew(dr_clust, clust->clust->nr+1);
821 for (i = 0; (i <= clust->clust->nr); i++)
823 init_dr_res(&dr_clust[i], fcd.disres.nres);
828 out = xvgropen(opt2fn("-ds", NFILE, fnm),
829 "Sum of Violations", "Time (ps)", "nm", oenv);
830 aver = xvgropen(opt2fn("-da", NFILE, fnm),
831 "Average Violation", "Time (ps)", "nm", oenv);
832 numv = xvgropen(opt2fn("-dn", NFILE, fnm),
833 "# Violations", "Time (ps)", "#", oenv);
834 maxxv = xvgropen(opt2fn("-dm", NFILE, fnm),
835 "Largest Violation", "Time (ps)", "nm", oenv);
838 mdatoms = init_mdatoms(fplog, &mtop, ir.efep != efepNO);
839 atoms2md(&mtop, &ir, 0, NULL, mtop.natoms, mdatoms);
840 update_mdatoms(mdatoms, ir.fepvals->init_lambda);
842 if (ir.ePBC != epbcNONE)
844 gpbc = gmx_rmpbc_init(&top->idef, ir.ePBC, natoms);
850 if (ir.ePBC != epbcNONE)
852 if (ir.bPeriodicMols)
854 set_pbc(&pbc, ir.ePBC, box);
858 gmx_rmpbc(gpbc, natoms, box, x);
864 if (j > clust->maxframe)
866 gmx_fatal(FARGS, "There are more frames in the trajectory than in the cluster index file. t = %8f\n", t);
868 my_clust = clust->inv_clust[j];
869 range_check(my_clust, 0, clust->clust->nr);
870 check_viol(fplog, &(top->idef.il[F_DISRES]),
872 x, f, pbc_null, g, dr_clust, my_clust, isize, index, vvindex, &fcd);
876 check_viol(fplog, &(top->idef.il[F_DISRES]),
878 x, f, pbc_null, g, &dr, 0, isize, index, vvindex, &fcd);
882 reset_x(atoms->nr, ind_fit, atoms->nr, NULL, x, w_rls);
883 do_fit(atoms->nr, w_rls, x, x);
886 /* Store the first frame of the trajectory as 'characteristic'
887 * for colouring with violations.
889 for (kkk = 0; (kkk < atoms->nr); kkk++)
891 copy_rvec(x[kkk], xav[kkk]);
899 fprintf(xvg, "%10g", t);
900 for (i = 0; (i < isize); i++)
902 fprintf(xvg, " %10g", vvindex[i]);
906 fprintf(out, "%10g %10g\n", t, dr.sumv);
907 fprintf(aver, "%10g %10g\n", t, dr.averv);
908 fprintf(maxxv, "%10g %10g\n", t, dr.maxv);
909 fprintf(numv, "%10g %10d\n", t, dr.nv);
913 while (read_next_x(oenv, status, &t, x, box));
915 if (ir.ePBC != epbcNONE)
917 gmx_rmpbc_done(gpbc);
922 dump_clust_stats(fplog, fcd.disres.nres, &(top->idef.il[F_DISRES]),
923 top->idef.iparams, clust->clust, dr_clust,
924 clust->grpname, isize, index);
928 dump_stats(fplog, j, fcd.disres.nres, &(top->idef.il[F_DISRES]),
929 top->idef.iparams, &dr, isize, index,
930 bPDB ? atoms : NULL);
933 write_sto_conf(opt2fn("-q", NFILE, fnm),
934 "Coloured by average violation in Angstrom",
935 atoms, xav, NULL, ir.ePBC, box);
937 dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres,
938 j, &top->idef, &mtop, max_dr, nlevels, bThird);
946 do_view(oenv, opt2fn("-dr", NFILE, fnm), "-nxy");
948 do_view(oenv, opt2fn("-dn", NFILE, fnm), "-nxy");
949 do_view(oenv, opt2fn("-da", NFILE, fnm), "-nxy");
950 do_view(oenv, opt2fn("-ds", NFILE, fnm), "-nxy");
951 do_view(oenv, opt2fn("-dm", NFILE, fnm), "-nxy");
954 gmx_log_close(fplog);