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 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
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.
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.
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.
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.
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.
45 #include <unordered_map>
47 #include "gromacs/commandline/pargs.h"
48 #include "gromacs/commandline/viewit.h"
49 #include "gromacs/fileio/confio.h"
50 #include "gromacs/fileio/matio.h"
51 #include "gromacs/fileio/pdbio.h"
52 #include "gromacs/fileio/tpxio.h"
53 #include "gromacs/fileio/trxio.h"
54 #include "gromacs/fileio/xvgr.h"
55 #include "gromacs/gmxana/gmx_ana.h"
56 #include "gromacs/gmxana/gstat.h"
57 #include "gromacs/listed_forces/disre.h"
58 #include "gromacs/math/do_fit.h"
59 #include "gromacs/math/functions.h"
60 #include "gromacs/math/vec.h"
61 #include "gromacs/mdlib/force.h"
62 #include "gromacs/mdlib/mdatoms.h"
63 #include "gromacs/mdtypes/fcdata.h"
64 #include "gromacs/mdtypes/inputrec.h"
65 #include "gromacs/mdtypes/md_enums.h"
66 #include "gromacs/mdtypes/mdatom.h"
67 #include "gromacs/pbcutil/ishift.h"
68 #include "gromacs/pbcutil/mshift.h"
69 #include "gromacs/pbcutil/pbc.h"
70 #include "gromacs/pbcutil/rmpbc.h"
71 #include "gromacs/topology/index.h"
72 #include "gromacs/topology/mtop_util.h"
73 #include "gromacs/topology/topology.h"
74 #include "gromacs/trajectoryanalysis/topologyinformation.h"
75 #include "gromacs/utility/arraysize.h"
76 #include "gromacs/utility/fatalerror.h"
77 #include "gromacs/utility/futil.h"
78 #include "gromacs/utility/smalloc.h"
86 static t_toppop* top = nullptr;
92 real sumv, averv, maxv;
93 real *aver1, *aver2, *aver_3, *aver_6;
96 static void init5(int n)
106 for (i = 0; (i < ntop); i++)
113 static void add5(int ndr, real viol)
118 for (i = 1; (i < ntop); i++)
120 if (top[i].v < top[mini].v)
125 if (viol > top[mini].v)
132 static void print5(FILE* fp)
136 std::sort(top, top + ntop, [](const t_toppop& a, const t_toppop& b) { return a.v > b.v; }); // reverse sort
137 fprintf(fp, "Index:");
138 for (i = 0; (i < ntop); i++)
140 fprintf(fp, " %6d", top[i].n);
142 fprintf(fp, "\nViol: ");
143 for (i = 0; (i < ntop); i++)
145 fprintf(fp, " %6.3f", top[i].v);
150 static void check_viol(FILE* log,
152 t_iparams forceparams[],
165 int i, j, nat, n, type, nviol, ndr, label;
166 real rt, mviol, tviol, viol, lam, dvdl, drt;
168 static gmx_bool bFirst = TRUE;
180 forceatoms = disres->iatoms;
181 for (j = 0; (j < isize); j++)
185 nat = interaction_function[F_DISRES].nratoms + 1;
186 for (i = 0; (i < disres->nr);)
188 type = forceatoms[i];
190 label = forceparams[type].disres.label;
193 fprintf(debug, "DISRE: ndr = %d, label = %d i=%d, n =%d\n", ndr, label, i, n);
198 } while (((i + n) < disres->nr) && (forceparams[forceatoms[i + n]].disres.label == label));
200 calc_disres_R_6(nullptr, nullptr, n, &forceatoms[i], x, pbc, fcd, nullptr);
202 if (fcd->disres.Rt_6[label] <= 0)
204 gmx_fatal(FARGS, "ndr = %d, rt_6 = %f", ndr, fcd->disres.Rt_6[label]);
207 rt = gmx::invsixthroot(fcd->disres.Rt_6[label]);
208 dr[clust_id].aver1[ndr] += rt;
209 dr[clust_id].aver2[ndr] += gmx::square(rt);
210 drt = 1.0 / gmx::power3(rt);
211 dr[clust_id].aver_3[ndr] += drt;
212 dr[clust_id].aver_6[ndr] += fcd->disres.Rt_6[label];
214 snew(fshift, SHIFTS);
215 ta_disres(n, &forceatoms[i], forceparams, x, f, fshift, pbc, g, lam, &dvdl, nullptr, fcd, nullptr);
217 viol = fcd->disres.sumviol;
224 add5(forceparams[type].disres.label, viol);
231 for (j = 0; (j < isize); j++)
233 if (index[j] == forceparams[type].disres.label)
235 vvindex[j] = gmx::invsixthroot(fcd->disres.Rt_6[label]);
242 dr[clust_id].nv = nviol;
243 dr[clust_id].maxv = mviol;
244 dr[clust_id].sumv = tviol;
245 dr[clust_id].averv = tviol / ndr;
246 dr[clust_id].nframes++;
250 fprintf(stderr, "\nThere are %d restraints and %d pairs\n", ndr, disres->nr / nat);
263 real up1, r, rT3, rT6, viol, violT3, violT6;
266 static void dump_dump(FILE* log, int ndr, t_dr_stats drs[])
268 static const char* core[] = { "All restraints", "Core restraints" };
269 static const char* tp[] = { "linear", "third power", "sixth power" };
270 real viol_tot, viol_max, viol = 0;
275 for (int iCore = 0; iCore < 2; iCore++)
277 bCore = (iCore == 1);
278 for (kkk = 0; (kkk < 3); kkk++)
284 for (i = 0; (i < ndr); i++)
286 if (!bCore || drs[i].bCore)
290 case 0: viol = drs[i].viol; break;
291 case 1: viol = drs[i].violT3; break;
292 case 2: viol = drs[i].violT6; break;
293 default: gmx_incons("Dumping violations");
295 viol_max = std::max(viol_max, viol);
304 if ((nrestr > 0) || (bCore && (nrestr < ndr)))
307 fprintf(log, "+++++++ %s ++++++++\n", core[bCore]);
308 fprintf(log, "+++++++ Using %s averaging: ++++++++\n", tp[kkk]);
309 fprintf(log, "Sum of violations: %8.3f nm\n", viol_tot);
312 fprintf(log, "Average violation: %8.3f nm\n", viol_tot / nrestr);
314 fprintf(log, "Largest violation: %8.3f nm\n", viol_max);
315 fprintf(log, "Number of violated restraints: %d/%d\n", nviol, nrestr);
321 static void dump_viol(FILE* log, int ndr, t_dr_stats* drs, gmx_bool bLinear)
325 fprintf(log, "Restr. Core Up1 <r> <rT3> <rT6> <viol><violT3><violT6>\n");
326 for (i = 0; (i < ndr); i++)
328 if (bLinear && (drs[i].viol == 0))
332 fprintf(log, "%6d%5s%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n", drs[i].label,
333 yesno_names[drs[i].bCore], drs[i].up1, drs[i].r, drs[i].rT3, drs[i].rT6,
334 drs[i].viol, drs[i].violT3, drs[i].violT6);
338 static gmx_bool is_core(int i, int isize, const int index[])
341 gmx_bool bIC = FALSE;
343 for (kk = 0; !bIC && (kk < isize); kk++)
345 bIC = (index[kk] == i);
351 static void dump_stats(FILE* log,
353 const t_disresdata& dd,
354 const t_ilist* disres,
364 fprintf(log, "++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
366 const int nra = interaction_function[F_DISRES].nratoms + 1;
367 for (int j = 0; j < disres->nr; j += nra)
369 // Note that the restraint i can be used by multiple pairs
370 const int i = disres->iatoms[j] - dd.type_min;
371 GMX_RELEASE_ASSERT(i >= 0 && i < dd.nres, "The restraint index should be in range");
373 drs[i].label = ip[disres->iatoms[j]].disres.label;
374 drs[i].bCore = is_core(drs[i].label, isize, index);
375 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
376 drs[i].r = dr->aver1[i] / nsteps;
377 drs[i].rT3 = gmx::invcbrt(dr->aver_3[i] / nsteps);
378 drs[i].rT6 = gmx::invsixthroot(dr->aver_6[i] / nsteps);
379 drs[i].viol = std::max(0.0, static_cast<double>(drs[i].r - drs[i].up1));
380 drs[i].violT3 = std::max(0.0, static_cast<double>(drs[i].rT3 - drs[i].up1));
381 drs[i].violT6 = std::max(0.0, static_cast<double>(drs[i].rT6 - drs[i].up1));
384 int j1 = disres->iatoms[j + 1];
385 int j2 = disres->iatoms[j + 2];
386 atoms->pdbinfo[j1].bfac += drs[i].violT3 * 5;
387 atoms->pdbinfo[j2].bfac += drs[i].violT3 * 5;
390 dump_viol(log, dd.nres, drs, FALSE);
392 fprintf(log, "+++ Sorted by linear averaged violations: +++\n");
393 std::sort(drs, drs + dd.nres,
394 [](const t_dr_stats& a, const t_dr_stats& b) { return a.viol > b.viol; }); // Reverse sort
395 dump_viol(log, dd.nres, drs, TRUE);
397 dump_dump(log, dd.nres, drs);
402 static void dump_clust_stats(FILE* fp,
403 const t_disresdata& dd,
404 const t_ilist* disres,
413 double sumV, maxV, sumVT3, sumVT6, maxVT3, maxVT6;
417 fprintf(fp, "++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
418 fprintf(fp, "Cluster NFrames SumV MaxV SumVT MaxVT SumVS MaxVS\n");
422 for (k = 0; (k < clust->nr); k++)
424 if (dr[k].nframes == 0)
428 if (dr[k].nframes != (clust->index[k + 1] - clust->index[k]))
431 "Inconsistency in cluster %s.\n"
432 "Found %d frames in trajectory rather than the expected %d\n",
433 clust_name[k], dr[k].nframes, clust->index[k + 1] - clust->index[k]);
437 gmx_fatal(FARGS, "Inconsistency with cluster %d. Invalid name", k);
439 nra = interaction_function[F_DISRES].nratoms + 1;
440 sumV = sumVT3 = sumVT6 = maxV = maxVT3 = maxVT6 = 0;
442 // Use a map to process each restraint only once while looping over all pairs
443 std::unordered_map<int, bool> restraintHasBeenProcessed;
444 for (int j = 0; j < dd.nres; j += nra)
446 // Note that the restraint i can be used by multiple pairs
447 const int i = disres->iatoms[j] - dd.type_min;
449 if (restraintHasBeenProcessed[i])
454 drs[i].label = ip[disres->iatoms[j]].disres.label;
455 drs[i].bCore = is_core(drs[i].label, isize, index);
456 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
457 drs[i].r = dr[k].aver1[i] / dr[k].nframes;
458 if ((dr[k].aver_3[i] <= 0) || !std::isfinite(dr[k].aver_3[i]))
460 gmx_fatal(FARGS, "dr[%d].aver_3[%d] = %f", k, i, dr[k].aver_3[i]);
462 drs[i].rT3 = gmx::invcbrt(dr[k].aver_3[i] / dr[k].nframes);
463 drs[i].rT6 = gmx::invsixthroot(dr[k].aver_6[i] / dr[k].nframes);
464 drs[i].viol = std::max(0.0, static_cast<double>(drs[i].r - drs[i].up1));
465 drs[i].violT3 = std::max(0.0, static_cast<double>(drs[i].rT3 - drs[i].up1));
466 drs[i].violT6 = std::max(0.0, static_cast<double>(drs[i].rT6 - drs[i].up1));
468 sumVT3 += drs[i].violT3;
469 sumVT6 += drs[i].violT6;
470 maxV = std::max(maxV, static_cast<double>(drs[i].viol));
471 maxVT3 = std::max(maxVT3, static_cast<double>(drs[i].violT3));
472 maxVT6 = std::max(maxVT6, static_cast<double>(drs[i].violT6));
474 // We have processed restraint i, mark it as such
475 restraintHasBeenProcessed[i] = true;
477 if (std::strcmp(clust_name[k], "1000") == 0)
481 fprintf(fp, "%-10s%6d%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n", clust_name[k],
482 dr[k].nframes, sumV, maxV, sumVT3, maxVT3, sumVT6, maxVT6);
488 static void init_dr_res(t_dr_result* dr, int ndr)
490 snew(dr->aver1, ndr + 1);
491 snew(dr->aver2, ndr + 1);
492 snew(dr->aver_3, ndr + 1);
493 snew(dr->aver_6, ndr + 1);
501 static void dump_disre_matrix(const char* fn,
506 const gmx_mtop_t* mtop,
513 int n_res, a_offset, mol, a;
514 int i, j, nra, nratoms, tp, ri, rj, index, nlabel, label;
516 real **matrix, *t_res, hi, *w_dr, rav, rviol;
517 t_rgb rlo = { 1, 1, 1 };
518 t_rgb rhi = { 0, 0, 0 };
524 snew(resnr, mtop->natoms);
527 for (const gmx_molblock_t& molb : mtop->molblock)
529 const t_atoms& atoms = mtop->moltype[molb.type].atoms;
530 for (mol = 0; mol < molb.nmol; mol++)
532 for (a = 0; a < atoms.nr; a++)
534 resnr[a_offset + a] = n_res + atoms.atom[a].resind;
537 a_offset += atoms.nr;
542 for (i = 0; (i < n_res); i++)
547 for (i = 0; (i < n_res); i++)
549 snew(matrix[i], n_res);
551 nratoms = interaction_function[F_DISRES].nratoms;
552 nra = (idef->il[F_DISRES].nr / (nratoms + 1));
558 for (i = 0; (i < idef->il[F_DISRES].nr); i += nratoms + 1)
560 tp = idef->il[F_DISRES].iatoms[i];
561 label = idef->iparams[tp].disres.label;
565 /* Set index pointer */
569 gmx_fatal(FARGS, "nlabel is %d, label = %d", nlabel, label);
573 gmx_fatal(FARGS, "ndr = %d, index = %d", ndr, index);
575 /* Update the weight */
576 w_dr[index] = 1.0 / nlabel;
585 printf("nlabel = %d, index = %d, ndr = %d\n", nlabel, index, ndr);
587 for (i = 0; (i < ndr); i++)
589 for (j = ptr[i]; (j < ptr[i + 1]); j += nratoms + 1)
591 tp = idef->il[F_DISRES].iatoms[j];
592 ai = idef->il[F_DISRES].iatoms[j + 1];
593 aj = idef->il[F_DISRES].iatoms[j + 2];
599 rav = gmx::invcbrt(dr->aver_3[i] / nsteps);
603 rav = dr->aver1[i] / nsteps;
607 fprintf(debug, "DR %d, atoms %d, %d, distance %g\n", i, ai, aj, rav);
609 rviol = std::max(0.0_real, rav - idef->iparams[tp].disres.up1);
610 matrix[ri][rj] += w_dr[i] * rviol;
611 matrix[rj][ri] += w_dr[i] * rviol;
612 hi = std::max(hi, matrix[ri][rj]);
613 hi = std::max(hi, matrix[rj][ri]);
623 printf("Warning: the maxdr that you have specified (%g) is smaller than\nthe largest "
624 "value in your simulation (%g)\n",
629 printf("Highest level in the matrix will be %g\n", hi);
630 fp = gmx_ffopen(fn, "w");
631 write_xpm(fp, 0, "Distance Violations", "<V> (nm)", "Residue", "Residue", n_res, n_res, t_res,
632 t_res, matrix, 0, hi, rlo, rhi, &nlevels);
636 int gmx_disre(int argc, char* argv[])
638 const char* desc[] = {
639 "[THISMODULE] computes violations of distance restraints.",
640 "The program always",
641 "computes the instantaneous violations rather than time-averaged,",
642 "because this analysis is done from a trajectory file afterwards",
643 "it does not make sense to use time averaging. However,",
644 "the time averaged values per restraint are given in the log file.[PAR]",
645 "An index file may be used to select specific restraints by index group label for",
647 "When the optional [TT]-q[tt] flag is given a [REF].pdb[ref] file coloured by the",
648 "amount of average violations.[PAR]",
649 "When the [TT]-c[tt] option is given, an index file will be read",
650 "containing the frames in your trajectory corresponding to the clusters",
651 "(defined in another manner) that you want to analyze. For these clusters",
652 "the program will compute average violations using the third power",
653 "averaging algorithm and print them in the log file."
656 static int nlevels = 20;
657 static real max_dr = 0;
658 static gmx_bool bThird = TRUE;
664 "Number of large violations that are stored in the log file every step" },
669 "Maximum distance violation in matrix output. If less than or equal to 0 the "
670 "maximum will be determined by the data." },
671 { "-nlevels", FALSE, etINT, { &nlevels }, "Number of levels in the matrix output" },
676 "Use inverse third power averaging or linear for matrix output" }
679 FILE *out = nullptr, *aver = nullptr, *numv = nullptr, *maxxv = nullptr, *xvg = nullptr;
686 rvec * x, *xav = nullptr;
691 int * index = nullptr, *ind_fit = nullptr;
693 t_cluster_ndx* clust = nullptr;
694 t_dr_result dr, *dr_clust = nullptr;
696 real * vvindex = nullptr, *w_rls = nullptr;
697 t_pbc pbc, *pbc_null;
700 gmx_output_env_t* oenv;
701 gmx_rmpbc_t gpbc = nullptr;
703 t_filenm fnm[] = { { efTPR, nullptr, nullptr, ffREAD }, { efTRX, "-f", nullptr, ffREAD },
704 { efXVG, "-ds", "drsum", ffWRITE }, { efXVG, "-da", "draver", ffWRITE },
705 { efXVG, "-dn", "drnum", ffWRITE }, { efXVG, "-dm", "drmax", ffWRITE },
706 { efXVG, "-dr", "restr", ffWRITE }, { efLOG, "-l", "disres", ffWRITE },
707 { efNDX, nullptr, "viol", ffOPTRD }, { efPDB, "-q", "viol", ffOPTWR },
708 { efNDX, "-c", "clust", ffOPTRD }, { efXPM, "-x", "matrix", ffOPTWR } };
709 #define NFILE asize(fnm)
711 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW, NFILE, fnm, asize(pa), pa,
712 asize(desc), desc, 0, nullptr, &oenv))
717 fplog = ftp2FILE(efLOG, NFILE, fnm, "w");
724 t_inputrec irInstance;
725 t_inputrec* ir = &irInstance;
727 gmx::TopologyInformation topInfo;
728 topInfo.fillFromInputFile(ftp2fn(efTPR, NFILE, fnm));
729 int ntopatoms = topInfo.mtop()->natoms;
731 bPDB = opt2bSet("-q", NFILE, fnm);
734 snew(xav, ntopatoms);
735 snew(ind_fit, ntopatoms);
736 snew(w_rls, ntopatoms);
737 for (kkk = 0; (kkk < ntopatoms); kkk++)
743 atoms = topInfo.copyAtoms();
745 if (atoms->pdbinfo == nullptr)
747 snew(atoms->pdbinfo, atoms->nr);
749 atoms->havePdbInfo = TRUE;
752 gmx_mtop_generate_local_top(*topInfo.mtop(), &top, ir->efep != efepNO);
756 if (ir->pbcType != PbcType::No)
758 if (ir->bPeriodicMols)
764 g = mk_graph(fplog, &top.idef, 0, ntopatoms, FALSE, FALSE);
768 if (ftp2bSet(efNDX, NFILE, fnm))
770 /* TODO: Nothing is written to this file if -c is provided, but it is
772 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
773 xvg = xvgropen(opt2fn("-dr", NFILE, fnm), "Individual Restraints", "Time (ps)", "nm", oenv);
774 snew(vvindex, isize);
776 for (i = 0; (i < isize); i++)
780 sprintf(leg[i], "index %d", index[i]);
782 xvgr_legend(xvg, isize, leg, oenv);
790 init_disres(fplog, topInfo.mtop(), ir, nullptr, nullptr, &fcd, nullptr, FALSE);
792 int natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
795 init_dr_res(&dr, fcd.disres.nres);
796 if (opt2bSet("-c", NFILE, fnm))
798 clust = cluster_index(fplog, opt2fn("-c", NFILE, fnm));
799 snew(dr_clust, clust->clust->nr + 1);
800 for (i = 0; (i <= clust->clust->nr); i++)
802 init_dr_res(&dr_clust[i], fcd.disres.nres);
807 out = xvgropen(opt2fn("-ds", NFILE, fnm), "Sum of Violations", "Time (ps)", "nm", oenv);
808 aver = xvgropen(opt2fn("-da", NFILE, fnm), "Average Violation", "Time (ps)", "nm", oenv);
809 numv = xvgropen(opt2fn("-dn", NFILE, fnm), "# Violations", "Time (ps)", "#", oenv);
810 maxxv = xvgropen(opt2fn("-dm", NFILE, fnm), "Largest Violation", "Time (ps)", "nm", oenv);
813 auto mdAtoms = gmx::makeMDAtoms(fplog, *topInfo.mtop(), *ir, false);
814 atoms2md(topInfo.mtop(), ir, -1, nullptr, ntopatoms, mdAtoms.get());
815 update_mdatoms(mdAtoms->mdatoms(), ir->fepvals->init_lambda);
816 if (ir->pbcType != PbcType::No)
818 gpbc = gmx_rmpbc_init(&top.idef, ir->pbcType, natoms);
824 if (ir->pbcType != PbcType::No)
826 if (ir->bPeriodicMols)
828 set_pbc(&pbc, ir->pbcType, box);
832 gmx_rmpbc(gpbc, natoms, box, x);
838 if (j > clust->maxframe)
841 "There are more frames in the trajectory than in the cluster index file. "
845 my_clust = clust->inv_clust[j];
846 range_check(my_clust, 0, clust->clust->nr);
847 check_viol(fplog, &(top.idef.il[F_DISRES]), top.idef.iparams, x, f, pbc_null, g,
848 dr_clust, my_clust, isize, index, vvindex, &fcd);
852 check_viol(fplog, &(top.idef.il[F_DISRES]), top.idef.iparams, x, f, pbc_null, g, &dr, 0,
853 isize, index, vvindex, &fcd);
857 reset_x(atoms->nr, ind_fit, atoms->nr, nullptr, x, w_rls);
858 do_fit(atoms->nr, w_rls, x, x);
861 /* Store the first frame of the trajectory as 'characteristic'
862 * for colouring with violations.
864 for (kkk = 0; (kkk < atoms->nr); kkk++)
866 copy_rvec(x[kkk], xav[kkk]);
874 fprintf(xvg, "%10g", t);
875 for (i = 0; (i < isize); i++)
877 fprintf(xvg, " %10g", vvindex[i]);
881 fprintf(out, "%10g %10g\n", t, dr.sumv);
882 fprintf(aver, "%10g %10g\n", t, dr.averv);
883 fprintf(maxxv, "%10g %10g\n", t, dr.maxv);
884 fprintf(numv, "%10g %10d\n", t, dr.nv);
887 } while (read_next_x(oenv, status, &t, x, box));
889 if (ir->pbcType != PbcType::No)
891 gmx_rmpbc_done(gpbc);
896 dump_clust_stats(fplog, fcd.disres, &(top.idef.il[F_DISRES]), top.idef.iparams,
897 clust->clust, dr_clust, clust->grpname, isize, index);
901 dump_stats(fplog, j, fcd.disres, &(top.idef.il[F_DISRES]), top.idef.iparams, &dr, isize,
902 index, bPDB ? atoms.get() : nullptr);
905 write_sto_conf(opt2fn("-q", NFILE, fnm), "Coloured by average violation in Angstrom",
906 atoms.get(), xav, nullptr, ir->pbcType, box);
908 dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres, j, &top.idef,
909 topInfo.mtop(), max_dr, nlevels, bThird);
914 do_view(oenv, opt2fn("-dn", NFILE, fnm), "-nxy");
915 do_view(oenv, opt2fn("-da", NFILE, fnm), "-nxy");
916 do_view(oenv, opt2fn("-ds", NFILE, fnm), "-nxy");
917 do_view(oenv, opt2fn("-dm", NFILE, fnm), "-nxy");
924 do_view(oenv, opt2fn("-dr", NFILE, fnm), "-nxy");