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,2015,2016,2017,2018, 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/commandline/pargs.h"
46 #include "gromacs/commandline/viewit.h"
47 #include "gromacs/fileio/confio.h"
48 #include "gromacs/fileio/matio.h"
49 #include "gromacs/fileio/pdbio.h"
50 #include "gromacs/fileio/tpxio.h"
51 #include "gromacs/fileio/trxio.h"
52 #include "gromacs/fileio/xvgr.h"
53 #include "gromacs/gmxana/gmx_ana.h"
54 #include "gromacs/gmxana/gstat.h"
55 #include "gromacs/gmxlib/nrnb.h"
56 #include "gromacs/listed-forces/disre.h"
57 #include "gromacs/math/do_fit.h"
58 #include "gromacs/math/functions.h"
59 #include "gromacs/math/vec.h"
60 #include "gromacs/mdlib/force.h"
61 #include "gromacs/mdlib/mdatoms.h"
62 #include "gromacs/mdlib/mdrun.h"
63 #include "gromacs/mdtypes/fcdata.h"
64 #include "gromacs/mdtypes/inputrec.h"
65 #include "gromacs/mdtypes/md_enums.h"
66 #include "gromacs/pbcutil/ishift.h"
67 #include "gromacs/pbcutil/mshift.h"
68 #include "gromacs/pbcutil/pbc.h"
69 #include "gromacs/pbcutil/rmpbc.h"
70 #include "gromacs/topology/index.h"
71 #include "gromacs/topology/mtop_util.h"
72 #include "gromacs/topology/topology.h"
73 #include "gromacs/utility/arraysize.h"
74 #include "gromacs/utility/fatalerror.h"
75 #include "gromacs/utility/futil.h"
76 #include "gromacs/utility/smalloc.h"
83 static t_toppop *top = nullptr;
88 real sumv, averv, maxv;
89 real *aver1, *aver2, *aver_3, *aver_6;
92 static void init5(int n)
102 for (i = 0; (i < ntop); i++)
109 static void add5(int ndr, real viol)
114 for (i = 1; (i < ntop); i++)
116 if (top[i].v < top[mini].v)
121 if (viol > top[mini].v)
128 static void print5(FILE *fp)
132 std::sort(top, top+ntop, [](const t_toppop &a, const t_toppop &b) {return a.v > b.v; }); //reverse sort
133 fprintf(fp, "Index:");
134 for (i = 0; (i < ntop); i++)
136 fprintf(fp, " %6d", top[i].n);
138 fprintf(fp, "\nViol: ");
139 for (i = 0; (i < ntop); i++)
141 fprintf(fp, " %6.3f", top[i].v);
146 static void check_viol(FILE *log,
147 t_ilist *disres, t_iparams forceparams[],
149 t_pbc *pbc, t_graph *g, t_dr_result dr[],
150 int clust_id, int isize, const int index[], real vvindex[],
154 int i, j, nat, n, type, nviol, ndr, label;
155 real rt, mviol, tviol, viol, lam, dvdl, drt;
157 static gmx_bool bFirst = TRUE;
169 forceatoms = disres->iatoms;
170 for (j = 0; (j < isize); j++)
174 nat = interaction_function[F_DISRES].nratoms+1;
175 for (i = 0; (i < disres->nr); )
177 type = forceatoms[i];
179 label = forceparams[type].disres.label;
182 fprintf(debug, "DISRE: ndr = %d, label = %d i=%d, n =%d\n",
187 gmx_fatal(FARGS, "tpr inconsistency. ndr = %d, label = %d\n", ndr, label);
193 while (((i+n) < disres->nr) &&
194 (forceparams[forceatoms[i+n]].disres.label == label));
196 calc_disres_R_6(nullptr, nullptr, n, &forceatoms[i],
197 x, pbc, fcd, nullptr);
199 if (fcd->disres.Rt_6[label] <= 0)
201 gmx_fatal(FARGS, "ndr = %d, rt_6 = %f", ndr, fcd->disres.Rt_6[label]);
204 rt = gmx::invsixthroot(fcd->disres.Rt_6[label]);
205 dr[clust_id].aver1[ndr] += rt;
206 dr[clust_id].aver2[ndr] += gmx::square(rt);
207 drt = 1.0/gmx::power3(rt);
208 dr[clust_id].aver_3[ndr] += drt;
209 dr[clust_id].aver_6[ndr] += fcd->disres.Rt_6[label];
211 snew(fshift, SHIFTS);
212 interaction_function[F_DISRES].ifunc(n, &forceatoms[i], forceparams,
214 pbc, g, lam, &dvdl, nullptr, fcd, nullptr);
216 viol = fcd->disres.sumviol;
223 add5(forceparams[type].disres.label, viol);
230 for (j = 0; (j < isize); j++)
232 if (index[j] == forceparams[type].disres.label)
234 vvindex[j] = gmx::invsixthroot(fcd->disres.Rt_6[label]);
241 dr[clust_id].nv = nviol;
242 dr[clust_id].maxv = mviol;
243 dr[clust_id].sumv = tviol;
244 dr[clust_id].averv = tviol/ndr;
245 dr[clust_id].nframes++;
249 fprintf(stderr, "\nThere are %d restraints and %d pairs\n",
250 ndr, disres->nr/nat);
262 real up1, r, rT3, rT6, viol, violT3, violT6;
265 static void dump_dump(FILE *log, int ndr, t_dr_stats drs[])
267 static const char *core[] = { "All restraints", "Core restraints" };
268 static const char *tp[] = { "linear", "third power", "sixth power" };
269 real viol_tot, viol_max, viol = 0;
274 for (int iCore = 0; iCore < 2; iCore++)
276 bCore = (iCore == 1);
277 for (kkk = 0; (kkk < 3); kkk++)
283 for (i = 0; (i < ndr); i++)
285 if (!bCore || drs[i].bCore)
293 viol = drs[i].violT3;
296 viol = drs[i].violT6;
299 gmx_incons("Dumping violations");
301 viol_max = std::max(viol_max, viol);
310 if ((nrestr > 0) || (bCore && (nrestr < ndr)))
313 fprintf(log, "+++++++ %s ++++++++\n", core[bCore]);
314 fprintf(log, "+++++++ Using %s averaging: ++++++++\n", tp[kkk]);
315 fprintf(log, "Sum of violations: %8.3f nm\n", viol_tot);
318 fprintf(log, "Average violation: %8.3f nm\n", viol_tot/nrestr);
320 fprintf(log, "Largest violation: %8.3f nm\n", viol_max);
321 fprintf(log, "Number of violated restraints: %d/%d\n", nviol, nrestr);
327 static void dump_viol(FILE *log, int ndr, t_dr_stats *drs, gmx_bool bLinear)
331 fprintf(log, "Restr. Core Up1 <r> <rT3> <rT6> <viol><violT3><violT6>\n");
332 for (i = 0; (i < ndr); i++)
334 if (bLinear && (drs[i].viol == 0))
338 fprintf(log, "%6d%5s%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n",
339 drs[i].index, yesno_names[drs[i].bCore],
340 drs[i].up1, drs[i].r, drs[i].rT3, drs[i].rT6,
341 drs[i].viol, drs[i].violT3, drs[i].violT6);
345 static gmx_bool is_core(int i, int isize, const int index[])
348 gmx_bool bIC = FALSE;
350 for (kk = 0; !bIC && (kk < isize); kk++)
352 bIC = (index[kk] == i);
358 static void dump_stats(FILE *log, int nsteps, int ndr, t_ilist *disres,
359 t_iparams ip[], t_dr_result *dr,
360 int isize, int index[], t_atoms *atoms)
366 fprintf(log, "++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
369 nra = interaction_function[F_DISRES].nratoms+1;
370 for (i = 0; (i < ndr); i++)
372 /* Search for the appropriate restraint */
373 for (; (j < disres->nr); j += nra)
375 if (ip[disres->iatoms[j]].disres.label == i)
381 drs[i].bCore = is_core(i, isize, index);
382 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
383 drs[i].r = dr->aver1[i]/nsteps;
384 drs[i].rT3 = gmx::invcbrt(dr->aver_3[i]/nsteps);
385 drs[i].rT6 = gmx::invsixthroot(dr->aver_6[i]/nsteps);
386 drs[i].viol = std::max(0.0, static_cast<double>(drs[i].r-drs[i].up1));
387 drs[i].violT3 = std::max(0.0, static_cast<double>(drs[i].rT3-drs[i].up1));
388 drs[i].violT6 = std::max(0.0, static_cast<double>(drs[i].rT6-drs[i].up1));
391 int j1 = disres->iatoms[j+1];
392 int j2 = disres->iatoms[j+2];
393 atoms->pdbinfo[j1].bfac += drs[i].violT3*5;
394 atoms->pdbinfo[j2].bfac += drs[i].violT3*5;
397 dump_viol(log, ndr, drs, FALSE);
399 fprintf(log, "+++ Sorted by linear averaged violations: +++\n");
400 std::sort(drs, drs+ndr, [](const t_dr_stats &a, const t_dr_stats &b)
401 {return a.viol > b.viol; }); //Reverse sort
402 dump_viol(log, ndr, drs, TRUE);
404 dump_dump(log, ndr, drs);
409 static void dump_clust_stats(FILE *fp, int ndr, t_ilist *disres,
410 t_iparams ip[], t_blocka *clust, t_dr_result dr[],
411 char *clust_name[], int isize, int index[])
413 int i, j, k, nra, mmm = 0;
414 double sumV, maxV, sumVT3, sumVT6, maxVT3, maxVT6;
418 fprintf(fp, "++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
419 fprintf(fp, "Cluster NFrames SumV MaxV SumVT MaxVT SumVS MaxVS\n");
423 for (k = 0; (k < clust->nr); k++)
425 if (dr[k].nframes == 0)
429 if (dr[k].nframes != (clust->index[k+1]-clust->index[k]))
431 gmx_fatal(FARGS, "Inconsistency in cluster %s.\n"
432 "Found %d frames in trajectory rather than the expected %d\n",
433 clust_name[k], dr[k].nframes,
434 clust->index[k+1]-clust->index[k]);
438 gmx_fatal(FARGS, "Inconsistency with cluster %d. Invalid name", k);
441 nra = interaction_function[F_DISRES].nratoms+1;
442 sumV = sumVT3 = sumVT6 = maxV = maxVT3 = maxVT6 = 0;
443 for (i = 0; (i < ndr); i++)
445 /* Search for the appropriate restraint */
446 for (; (j < disres->nr); j += nra)
448 if (ip[disres->iatoms[j]].disres.label == i)
454 drs[i].bCore = is_core(i, isize, index);
455 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
456 drs[i].r = dr[k].aver1[i]/dr[k].nframes;
457 if ((dr[k].aver_3[i] <= 0) || !std::isfinite(dr[k].aver_3[i]))
459 gmx_fatal(FARGS, "dr[%d].aver_3[%d] = %f", k, i, dr[k].aver_3[i]);
461 drs[i].rT3 = gmx::invcbrt(dr[k].aver_3[i]/dr[k].nframes);
462 drs[i].rT6 = gmx::invsixthroot(dr[k].aver_6[i]/dr[k].nframes);
463 drs[i].viol = std::max(0.0, static_cast<double>(drs[i].r-drs[i].up1));
464 drs[i].violT3 = std::max(0.0, static_cast<double>(drs[i].rT3-drs[i].up1));
465 drs[i].violT6 = std::max(0.0, static_cast<double>(drs[i].rT6-drs[i].up1));
467 sumVT3 += drs[i].violT3;
468 sumVT6 += drs[i].violT6;
469 maxV = std::max(maxV, static_cast<double>(drs[i].viol));
470 maxVT3 = std::max(maxVT3, static_cast<double>(drs[i].violT3));
471 maxVT6 = std::max(maxVT6, static_cast<double>(drs[i].violT6));
473 if (std::strcmp(clust_name[k], "1000") == 0)
477 fprintf(fp, "%-10s%6d%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n",
479 dr[k].nframes, sumV, maxV, sumVT3, maxVT3, sumVT6, maxVT6);
486 static void init_dr_res(t_dr_result *dr, int ndr)
488 snew(dr->aver1, ndr+1);
489 snew(dr->aver2, ndr+1);
490 snew(dr->aver_3, ndr+1);
491 snew(dr->aver_6, ndr+1);
499 static void dump_disre_matrix(const char *fn, t_dr_result *dr, int ndr,
500 int nsteps, t_idef *idef, const gmx_mtop_t *mtop,
501 real max_dr, int nlevels, gmx_bool bThird)
505 int n_res, a_offset, mol, a;
506 int i, j, nra, nratoms, tp, ri, rj, index, nlabel, label;
508 real **matrix, *t_res, hi, *w_dr, rav, rviol;
509 t_rgb rlo = { 1, 1, 1 };
510 t_rgb rhi = { 0, 0, 0 };
516 snew(resnr, mtop->natoms);
519 for (const gmx_molblock_t &molb : mtop->molblock)
521 const t_atoms &atoms = mtop->moltype[molb.type].atoms;
522 for (mol = 0; mol < molb.nmol; mol++)
524 for (a = 0; a < atoms.nr; a++)
526 resnr[a_offset + a] = n_res + atoms.atom[a].resind;
529 a_offset += atoms.nr;
534 for (i = 0; (i < n_res); i++)
539 for (i = 0; (i < n_res); i++)
541 snew(matrix[i], n_res);
543 nratoms = interaction_function[F_DISRES].nratoms;
544 nra = (idef->il[F_DISRES].nr/(nratoms+1));
550 for (i = 0; (i < idef->il[F_DISRES].nr); i += nratoms+1)
552 tp = idef->il[F_DISRES].iatoms[i];
553 label = idef->iparams[tp].disres.label;
557 /* Set index pointer */
561 gmx_fatal(FARGS, "nlabel is %d, label = %d", nlabel, label);
565 gmx_fatal(FARGS, "ndr = %d, index = %d", ndr, index);
567 /* Update the weight */
568 w_dr[index] = 1.0/nlabel;
577 printf("nlabel = %d, index = %d, ndr = %d\n", nlabel, index, ndr);
579 for (i = 0; (i < ndr); i++)
581 for (j = ptr[i]; (j < ptr[i+1]); j += nratoms+1)
583 tp = idef->il[F_DISRES].iatoms[j];
584 ai = idef->il[F_DISRES].iatoms[j+1];
585 aj = idef->il[F_DISRES].iatoms[j+2];
591 rav = gmx::invcbrt(dr->aver_3[i]/nsteps);
595 rav = dr->aver1[i]/nsteps;
599 fprintf(debug, "DR %d, atoms %d, %d, distance %g\n", i, ai, aj, rav);
601 rviol = std::max(0.0_real, rav-idef->iparams[tp].disres.up1);
602 matrix[ri][rj] += w_dr[i]*rviol;
603 matrix[rj][ri] += w_dr[i]*rviol;
604 hi = std::max(hi, matrix[ri][rj]);
605 hi = std::max(hi, matrix[rj][ri]);
615 printf("Warning: the maxdr that you have specified (%g) is smaller than\nthe largest value in your simulation (%g)\n", max_dr, hi);
619 printf("Highest level in the matrix will be %g\n", hi);
620 fp = gmx_ffopen(fn, "w");
621 write_xpm(fp, 0, "Distance Violations", "<V> (nm)", "Residue", "Residue",
622 n_res, n_res, t_res, t_res, matrix, 0, hi, rlo, rhi, &nlevels);
626 int gmx_disre(int argc, char *argv[])
628 const char *desc[] = {
629 "[THISMODULE] computes violations of distance restraints.",
630 "The program always",
631 "computes the instantaneous violations rather than time-averaged,",
632 "because this analysis is done from a trajectory file afterwards",
633 "it does not make sense to use time averaging. However,",
634 "the time averaged values per restraint are given in the log file.[PAR]",
635 "An index file may be used to select specific restraints for",
637 "When the optional [TT]-q[tt] flag is given a [REF].pdb[ref] file coloured by the",
638 "amount of average violations.[PAR]",
639 "When the [TT]-c[tt] option is given, an index file will be read",
640 "containing the frames in your trajectory corresponding to the clusters",
641 "(defined in another manner) that you want to analyze. For these clusters",
642 "the program will compute average violations using the third power",
643 "averaging algorithm and print them in the log file."
646 static int nlevels = 20;
647 static real max_dr = 0;
648 static gmx_bool bThird = TRUE;
650 { "-ntop", FALSE, etINT, {&ntop},
651 "Number of large violations that are stored in the log file every step" },
652 { "-maxdr", FALSE, etREAL, {&max_dr},
653 "Maximum distance violation in matrix output. If less than or equal to 0 the maximum will be determined by the data." },
654 { "-nlevels", FALSE, etINT, {&nlevels},
655 "Number of levels in the matrix output" },
656 { "-third", FALSE, etBOOL, {&bThird},
657 "Use inverse third power averaging or linear for matrix output" }
660 FILE *out = nullptr, *aver = nullptr, *numv = nullptr, *maxxv = nullptr, *xvg = nullptr;
665 t_atoms *atoms = nullptr;
669 int ntopatoms, natoms, i, j, kkk;
672 rvec *x, *xav = nullptr;
677 int *index = nullptr, *ind_fit = nullptr;
679 t_cluster_ndx *clust = nullptr;
680 t_dr_result dr, *dr_clust = nullptr;
682 real *vvindex = nullptr, *w_rls = nullptr;
683 t_pbc pbc, *pbc_null;
686 gmx_output_env_t *oenv;
687 gmx_rmpbc_t gpbc = nullptr;
690 { efTPR, nullptr, nullptr, ffREAD },
691 { efTRX, "-f", nullptr, ffREAD },
692 { efXVG, "-ds", "drsum", ffWRITE },
693 { efXVG, "-da", "draver", ffWRITE },
694 { efXVG, "-dn", "drnum", ffWRITE },
695 { efXVG, "-dm", "drmax", ffWRITE },
696 { efXVG, "-dr", "restr", ffWRITE },
697 { efLOG, "-l", "disres", ffWRITE },
698 { efNDX, nullptr, "viol", ffOPTRD },
699 { efPDB, "-q", "viol", ffOPTWR },
700 { efNDX, "-c", "clust", ffOPTRD },
701 { efXPM, "-x", "matrix", ffOPTWR }
703 #define NFILE asize(fnm)
705 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
706 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
711 fplog = ftp2FILE(efLOG, NFILE, fnm, "w");
718 t_inputrec irInstance;
719 t_inputrec *ir = &irInstance;
721 read_tpxheader(ftp2fn(efTPR, NFILE, fnm), &header, FALSE);
722 snew(xtop, header.natoms);
723 read_tpx(ftp2fn(efTPR, NFILE, fnm), ir, box, &ntopatoms, xtop, nullptr, &mtop);
724 bPDB = opt2bSet("-q", NFILE, fnm);
727 snew(xav, ntopatoms);
728 snew(ind_fit, ntopatoms);
729 snew(w_rls, ntopatoms);
730 for (kkk = 0; (kkk < ntopatoms); kkk++)
737 *atoms = gmx_mtop_global_atoms(&mtop);
739 if (atoms->pdbinfo == nullptr)
741 snew(atoms->pdbinfo, atoms->nr);
743 atoms->havePdbInfo = TRUE;
746 top = gmx_mtop_generate_local_top(&mtop, ir->efep != efepNO);
750 if (ir->ePBC != epbcNONE)
752 if (ir->bPeriodicMols)
758 g = mk_graph(fplog, &top->idef, 0, mtop.natoms, FALSE, FALSE);
762 if (ftp2bSet(efNDX, NFILE, fnm))
764 /* TODO: Nothing is written to this file if -c is provided, but it is
766 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
767 xvg = xvgropen(opt2fn("-dr", NFILE, fnm), "Individual Restraints", "Time (ps)",
769 snew(vvindex, isize);
771 for (i = 0; (i < isize); i++)
775 sprintf(leg[i], "index %d", index[i]);
777 xvgr_legend(xvg, isize, leg, oenv);
785 init_disres(fplog, &mtop, ir, nullptr, nullptr, &fcd, nullptr, FALSE);
787 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
790 init_dr_res(&dr, fcd.disres.nres);
791 if (opt2bSet("-c", NFILE, fnm))
793 clust = cluster_index(fplog, opt2fn("-c", NFILE, fnm));
794 snew(dr_clust, clust->clust->nr+1);
795 for (i = 0; (i <= clust->clust->nr); i++)
797 init_dr_res(&dr_clust[i], fcd.disres.nres);
802 out = xvgropen(opt2fn("-ds", NFILE, fnm),
803 "Sum of Violations", "Time (ps)", "nm", oenv);
804 aver = xvgropen(opt2fn("-da", NFILE, fnm),
805 "Average Violation", "Time (ps)", "nm", oenv);
806 numv = xvgropen(opt2fn("-dn", NFILE, fnm),
807 "# Violations", "Time (ps)", "#", oenv);
808 maxxv = xvgropen(opt2fn("-dm", NFILE, fnm),
809 "Largest Violation", "Time (ps)", "nm", oenv);
812 auto mdAtoms = gmx::makeMDAtoms(fplog, mtop, *ir, false);
813 atoms2md(&mtop, ir, -1, nullptr, mtop.natoms, mdAtoms.get());
814 update_mdatoms(mdAtoms->mdatoms(), ir->fepvals->init_lambda);
816 if (ir->ePBC != epbcNONE)
818 gpbc = gmx_rmpbc_init(&top->idef, ir->ePBC, natoms);
824 if (ir->ePBC != epbcNONE)
826 if (ir->bPeriodicMols)
828 set_pbc(&pbc, ir->ePBC, box);
832 gmx_rmpbc(gpbc, natoms, box, x);
838 if (j > clust->maxframe)
840 gmx_fatal(FARGS, "There are more frames in the trajectory than in the cluster index file. t = %8f\n", t);
842 my_clust = clust->inv_clust[j];
843 range_check(my_clust, 0, clust->clust->nr);
844 check_viol(fplog, &(top->idef.il[F_DISRES]),
846 x, f, pbc_null, g, dr_clust, my_clust, isize, index, vvindex, &fcd);
850 check_viol(fplog, &(top->idef.il[F_DISRES]),
852 x, f, pbc_null, g, &dr, 0, isize, index, vvindex, &fcd);
856 reset_x(atoms->nr, ind_fit, atoms->nr, nullptr, x, w_rls);
857 do_fit(atoms->nr, w_rls, x, x);
860 /* Store the first frame of the trajectory as 'characteristic'
861 * for colouring with violations.
863 for (kkk = 0; (kkk < atoms->nr); kkk++)
865 copy_rvec(x[kkk], xav[kkk]);
873 fprintf(xvg, "%10g", t);
874 for (i = 0; (i < isize); i++)
876 fprintf(xvg, " %10g", vvindex[i]);
880 fprintf(out, "%10g %10g\n", t, dr.sumv);
881 fprintf(aver, "%10g %10g\n", t, dr.averv);
882 fprintf(maxxv, "%10g %10g\n", t, dr.maxv);
883 fprintf(numv, "%10g %10d\n", t, dr.nv);
887 while (read_next_x(oenv, status, &t, x, box));
889 if (ir->ePBC != epbcNONE)
891 gmx_rmpbc_done(gpbc);
896 dump_clust_stats(fplog, fcd.disres.nres, &(top->idef.il[F_DISRES]),
897 top->idef.iparams, clust->clust, dr_clust,
898 clust->grpname, isize, index);
902 dump_stats(fplog, j, fcd.disres.nres, &(top->idef.il[F_DISRES]),
903 top->idef.iparams, &dr, isize, index,
904 bPDB ? atoms : nullptr);
907 write_sto_conf(opt2fn("-q", NFILE, fnm),
908 "Coloured by average violation in Angstrom",
909 atoms, xav, nullptr, ir->ePBC, box);
911 dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres,
912 j, &top->idef, &mtop, max_dr, nlevels, bThird);
917 do_view(oenv, opt2fn("-dn", NFILE, fnm), "-nxy");
918 do_view(oenv, opt2fn("-da", NFILE, fnm), "-nxy");
919 do_view(oenv, opt2fn("-ds", NFILE, fnm), "-nxy");
920 do_view(oenv, opt2fn("-dm", NFILE, fnm), "-nxy");
927 do_view(oenv, opt2fn("-dr", NFILE, fnm), "-nxy");