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 The GROMACS development team.
7 * Copyright (c) 2018,2019,2020,2021, 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/commrec.h"
64 #include "gromacs/mdtypes/fcdata.h"
65 #include "gromacs/mdtypes/inputrec.h"
66 #include "gromacs/mdtypes/md_enums.h"
67 #include "gromacs/mdtypes/mdatom.h"
68 #include "gromacs/pbcutil/ishift.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,
151 const InteractionList& disres,
152 gmx::ArrayRef<const t_iparams> forceparams,
161 t_disresdata* disresdata)
163 int i, j, nat, n, type, nviol, ndr, label;
164 real rt, mviol, tviol, viol, lam, dvdl, drt;
166 static gmx_bool bFirst = TRUE;
178 gmx::ArrayRef<const int> forceatoms = disres.iatoms;
179 for (j = 0; (j < isize); j++)
183 nat = interaction_function[F_DISRES].nratoms + 1;
184 // Check internal consistency of disres.label
185 // The label for a distance restraint should be at most one larger
186 // than the previous label.
187 int label_old = forceparams[forceatoms[0]].disres.label;
188 for (i = 0; (i < disres.size()); i += nat)
190 type = forceatoms[i];
191 label = forceparams[type].disres.label;
192 if ((label == label_old) || (label == label_old + 1))
199 "Label mismatch in distance restrains. Label for restraint %d is %d, "
200 "expected it to be either %d or %d",
208 // Get offset for label index
209 label_old = forceparams[forceatoms[0]].disres.label;
210 for (i = 0; (i < disres.size());)
212 type = forceatoms[i];
214 label = forceparams[type].disres.label - label_old;
217 fprintf(debug, "DISRE: ndr = %d, label = %d i=%d, n =%d\n", ndr, label, i, n);
222 } while (((i + n) < disres.size())
223 && (forceparams[forceatoms[i + n]].disres.label == label + label_old));
225 calc_disres_R_6(nullptr, nullptr, n, &forceatoms[i], x, pbc, disresdata, nullptr);
227 if (disresdata->Rt_6[label] <= 0)
229 gmx_fatal(FARGS, "ndr = %d, rt_6 = %f", ndr, disresdata->Rt_6[label]);
232 rt = gmx::invsixthroot(disresdata->Rt_6[label]);
233 dr[clust_id].aver1[ndr] += rt;
234 dr[clust_id].aver2[ndr] += gmx::square(rt);
235 drt = 1.0 / gmx::power3(rt);
236 dr[clust_id].aver_3[ndr] += drt;
237 dr[clust_id].aver_6[ndr] += disresdata->Rt_6[label];
239 snew(fshift, SHIFTS);
240 ta_disres(n, &forceatoms[i], forceparams.data(), x, f, fshift, pbc, lam, &dvdl, nullptr, nullptr, disresdata, nullptr, nullptr);
242 viol = disresdata->sumviol;
249 add5(forceparams[type].disres.label, viol);
256 for (j = 0; (j < isize); j++)
258 if (index[j] == forceparams[type].disres.label)
260 vvindex[j] = gmx::invsixthroot(disresdata->Rt_6[label]);
267 dr[clust_id].nv = nviol;
268 dr[clust_id].maxv = mviol;
269 dr[clust_id].sumv = tviol;
270 dr[clust_id].averv = tviol / ndr;
271 dr[clust_id].nframes++;
275 fprintf(stderr, "\nThere are %d restraints and %d pairs\n", ndr, disres.size() / nat);
288 real up1, r, rT3, rT6, viol, violT3, violT6;
291 static void dump_dump(FILE* log, int ndr, t_dr_stats drs[])
293 static const char* core[] = { "All restraints", "Core restraints" };
294 static const char* tp[] = { "linear", "third power", "sixth power" };
295 real viol_tot, viol_max, viol = 0;
300 for (int iCore = 0; iCore < 2; iCore++)
302 bCore = (iCore == 1);
303 for (kkk = 0; (kkk < 3); kkk++)
309 for (i = 0; (i < ndr); i++)
311 if (!bCore || drs[i].bCore)
315 case 0: viol = drs[i].viol; break;
316 case 1: viol = drs[i].violT3; break;
317 case 2: viol = drs[i].violT6; break;
318 default: gmx_incons("Dumping violations");
320 viol_max = std::max(viol_max, viol);
329 if ((nrestr > 0) || (bCore && (nrestr < ndr)))
332 fprintf(log, "+++++++ %s ++++++++\n", core[bCore]);
333 fprintf(log, "+++++++ Using %s averaging: ++++++++\n", tp[kkk]);
334 fprintf(log, "Sum of violations: %8.3f nm\n", viol_tot);
337 fprintf(log, "Average violation: %8.3f nm\n", viol_tot / nrestr);
339 fprintf(log, "Largest violation: %8.3f nm\n", viol_max);
340 fprintf(log, "Number of violated restraints: %d/%d\n", nviol, nrestr);
346 static void dump_viol(FILE* log, int ndr, t_dr_stats* drs, gmx_bool bLinear)
350 fprintf(log, "Restr. Core Up1 <r> <rT3> <rT6> <viol><violT3><violT6>\n");
351 for (i = 0; (i < ndr); i++)
353 if (bLinear && (drs[i].viol == 0))
358 "%6d%5s%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n",
360 yesno_names[drs[i].bCore],
371 static gmx_bool is_core(int i, int isize, const int index[])
374 gmx_bool bIC = FALSE;
376 for (kk = 0; !bIC && (kk < isize); kk++)
378 bIC = (index[kk] == i);
384 static void dump_stats(FILE* log,
386 const t_disresdata& dd,
387 const InteractionList& disres,
388 gmx::ArrayRef<const t_iparams> ip,
397 fprintf(log, "++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
399 const int nra = interaction_function[F_DISRES].nratoms + 1;
400 for (int j = 0; j < disres.size(); j += nra)
402 // Note that the restraint i can be used by multiple pairs
403 const int i = disres.iatoms[j] - dd.type_min;
404 GMX_RELEASE_ASSERT(i >= 0 && i < dd.nres, "The restraint index should be in range");
406 drs[i].label = ip[disres.iatoms[j]].disres.label;
407 drs[i].bCore = is_core(drs[i].label, isize, index);
408 drs[i].up1 = ip[disres.iatoms[j]].disres.up1;
409 drs[i].r = dr->aver1[i] / nsteps;
410 drs[i].rT3 = gmx::invcbrt(dr->aver_3[i] / nsteps);
411 drs[i].rT6 = gmx::invsixthroot(dr->aver_6[i] / nsteps);
412 drs[i].viol = std::max(0.0, static_cast<double>(drs[i].r - drs[i].up1));
413 drs[i].violT3 = std::max(0.0, static_cast<double>(drs[i].rT3 - drs[i].up1));
414 drs[i].violT6 = std::max(0.0, static_cast<double>(drs[i].rT6 - drs[i].up1));
417 int j1 = disres.iatoms[j + 1];
418 int j2 = disres.iatoms[j + 2];
419 atoms->pdbinfo[j1].bfac += drs[i].violT3 * 5;
420 atoms->pdbinfo[j2].bfac += drs[i].violT3 * 5;
423 dump_viol(log, dd.nres, drs, FALSE);
425 fprintf(log, "+++ Sorted by linear averaged violations: +++\n");
426 std::sort(drs, drs + dd.nres, [](const t_dr_stats& a, const t_dr_stats& b) {
427 return a.viol > b.viol;
429 dump_viol(log, dd.nres, drs, TRUE);
431 dump_dump(log, dd.nres, drs);
436 static void dump_clust_stats(FILE* fp,
437 const t_disresdata& dd,
438 const InteractionList& disres,
439 gmx::ArrayRef<const t_iparams> ip,
447 double sumV, maxV, sumVT3, sumVT6, maxVT3, maxVT6;
451 fprintf(fp, "++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
452 fprintf(fp, "Cluster NFrames SumV MaxV SumVT MaxVT SumVS MaxVS\n");
456 for (k = 0; (k < clust->nr); k++)
458 if (dr[k].nframes == 0)
462 if (dr[k].nframes != (clust->index[k + 1] - clust->index[k]))
465 "Inconsistency in cluster %s.\n"
466 "Found %d frames in trajectory rather than the expected %d\n",
469 clust->index[k + 1] - clust->index[k]);
473 gmx_fatal(FARGS, "Inconsistency with cluster %d. Invalid name", k);
475 nra = interaction_function[F_DISRES].nratoms + 1;
476 sumV = sumVT3 = sumVT6 = maxV = maxVT3 = maxVT6 = 0;
478 // Use a map to process each restraint only once while looping over all pairs
479 std::unordered_map<int, bool> restraintHasBeenProcessed;
480 for (int j = 0; j < dd.nres; j += nra)
482 // Note that the restraint i can be used by multiple pairs
483 const int i = disres.iatoms[j] - dd.type_min;
485 if (restraintHasBeenProcessed[i])
490 drs[i].label = ip[disres.iatoms[j]].disres.label;
491 drs[i].bCore = is_core(drs[i].label, isize, index);
492 drs[i].up1 = ip[disres.iatoms[j]].disres.up1;
493 drs[i].r = dr[k].aver1[i] / dr[k].nframes;
494 if ((dr[k].aver_3[i] <= 0) || !std::isfinite(dr[k].aver_3[i]))
496 gmx_fatal(FARGS, "dr[%d].aver_3[%d] = %f", k, i, dr[k].aver_3[i]);
498 drs[i].rT3 = gmx::invcbrt(dr[k].aver_3[i] / dr[k].nframes);
499 drs[i].rT6 = gmx::invsixthroot(dr[k].aver_6[i] / dr[k].nframes);
500 drs[i].viol = std::max(0.0, static_cast<double>(drs[i].r - drs[i].up1));
501 drs[i].violT3 = std::max(0.0, static_cast<double>(drs[i].rT3 - drs[i].up1));
502 drs[i].violT6 = std::max(0.0, static_cast<double>(drs[i].rT6 - drs[i].up1));
504 sumVT3 += drs[i].violT3;
505 sumVT6 += drs[i].violT6;
506 maxV = std::max(maxV, static_cast<double>(drs[i].viol));
507 maxVT3 = std::max(maxVT3, static_cast<double>(drs[i].violT3));
508 maxVT6 = std::max(maxVT6, static_cast<double>(drs[i].violT6));
510 // We have processed restraint i, mark it as such
511 restraintHasBeenProcessed[i] = true;
513 if (std::strcmp(clust_name[k], "1000") == 0)
518 "%-10s%6d%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n",
532 static void init_dr_res(t_dr_result* dr, int ndr)
534 snew(dr->aver1, ndr + 1);
535 snew(dr->aver2, ndr + 1);
536 snew(dr->aver_3, ndr + 1);
537 snew(dr->aver_6, ndr + 1);
545 static void dump_disre_matrix(const char* fn,
549 const InteractionDefinitions& idef,
550 const gmx_mtop_t* mtop,
557 int n_res, a_offset, mol, a;
558 int i, j, nra, nratoms, tp, ri, rj, index, nlabel, label;
560 real **matrix, *t_res, hi, *w_dr, rav, rviol;
561 t_rgb rlo = { 1, 1, 1 };
562 t_rgb rhi = { 0, 0, 0 };
568 snew(resnr, mtop->natoms);
571 for (const gmx_molblock_t& molb : mtop->molblock)
573 const t_atoms& atoms = mtop->moltype[molb.type].atoms;
574 for (mol = 0; mol < molb.nmol; mol++)
576 for (a = 0; a < atoms.nr; a++)
578 resnr[a_offset + a] = n_res + atoms.atom[a].resind;
581 a_offset += atoms.nr;
586 for (i = 0; (i < n_res); i++)
591 for (i = 0; (i < n_res); i++)
593 snew(matrix[i], n_res);
595 nratoms = interaction_function[F_DISRES].nratoms;
596 nra = (idef.il[F_DISRES].size() / (nratoms + 1));
602 for (i = 0; (i < idef.il[F_DISRES].size()); i += nratoms + 1)
604 tp = idef.il[F_DISRES].iatoms[i];
605 label = idef.iparams[tp].disres.label;
609 /* Set index pointer */
613 gmx_fatal(FARGS, "nlabel is %d, label = %d", nlabel, label);
617 gmx_fatal(FARGS, "ndr = %d, index = %d", ndr, index);
619 /* Update the weight */
620 w_dr[index] = 1.0 / nlabel;
629 printf("nlabel = %d, index = %d, ndr = %d\n", nlabel, index, ndr);
631 for (i = 0; (i < ndr); i++)
633 for (j = ptr[i]; (j < ptr[i + 1]); j += nratoms + 1)
635 tp = idef.il[F_DISRES].iatoms[j];
636 ai = idef.il[F_DISRES].iatoms[j + 1];
637 aj = idef.il[F_DISRES].iatoms[j + 2];
643 rav = gmx::invcbrt(dr->aver_3[i] / nsteps);
647 rav = dr->aver1[i] / nsteps;
651 fprintf(debug, "DR %d, atoms %d, %d, distance %g\n", i, ai, aj, rav);
653 rviol = std::max(0.0_real, rav - idef.iparams[tp].disres.up1);
654 matrix[ri][rj] += w_dr[i] * rviol;
655 matrix[rj][ri] += w_dr[i] * rviol;
656 hi = std::max(hi, matrix[ri][rj]);
657 hi = std::max(hi, matrix[rj][ri]);
667 printf("Warning: the maxdr that you have specified (%g) is smaller than\nthe largest "
668 "value in your simulation (%g)\n",
674 printf("Highest level in the matrix will be %g\n", hi);
675 fp = gmx_ffopen(fn, "w");
678 "Distance Violations",
695 int gmx_disre(int argc, char* argv[])
697 const char* desc[] = {
698 "[THISMODULE] computes violations of distance restraints.",
699 "The program always",
700 "computes the instantaneous violations rather than time-averaged,",
701 "because this analysis is done from a trajectory file afterwards",
702 "it does not make sense to use time averaging. However,",
703 "the time averaged values per restraint are given in the log file.[PAR]",
704 "An index file may be used to select specific restraints by index group label for",
706 "When the optional [TT]-q[tt] flag is given a [REF].pdb[ref] file coloured by the",
707 "amount of average violations.[PAR]",
708 "When the [TT]-c[tt] option is given, an index file will be read",
709 "containing the frames in your trajectory corresponding to the clusters",
710 "(defined in another manner) that you want to analyze. For these clusters",
711 "the program will compute average violations using the third power",
712 "averaging algorithm and print them in the log file."
715 static int nlevels = 20;
716 static real max_dr = 0;
717 static gmx_bool bThird = TRUE;
723 "Number of large violations that are stored in the log file every step" },
728 "Maximum distance violation in matrix output. If less than or equal to 0 the "
729 "maximum will be determined by the data." },
730 { "-nlevels", FALSE, etINT, { &nlevels }, "Number of levels in the matrix output" },
735 "Use inverse third power averaging or linear for matrix output" }
738 FILE * out = nullptr, *aver = nullptr, *numv = nullptr, *maxxv = nullptr, *xvg = nullptr;
742 rvec * x, *xav = nullptr;
747 int * index = nullptr, *ind_fit = nullptr;
749 t_cluster_ndx* clust = nullptr;
750 t_dr_result dr, *dr_clust = nullptr;
752 real * vvindex = nullptr, *w_rls = nullptr;
753 t_pbc pbc, *pbc_null;
756 gmx_output_env_t* oenv;
757 gmx_rmpbc_t gpbc = nullptr;
759 t_filenm fnm[] = { { efTPR, nullptr, nullptr, ffREAD }, { efTRX, "-f", nullptr, ffREAD },
760 { efXVG, "-ds", "drsum", ffWRITE }, { efXVG, "-da", "draver", ffWRITE },
761 { efXVG, "-dn", "drnum", ffWRITE }, { efXVG, "-dm", "drmax", ffWRITE },
762 { efXVG, "-dr", "restr", ffWRITE }, { efLOG, "-l", "disres", ffWRITE },
763 { efNDX, nullptr, "viol", ffOPTRD }, { efPDB, "-q", "viol", ffOPTWR },
764 { efNDX, "-c", "clust", ffOPTRD }, { efXPM, "-x", "matrix", ffOPTWR } };
765 #define NFILE asize(fnm)
767 if (!parse_common_args(
768 &argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
773 fplog = ftp2FILE(efLOG, NFILE, fnm, "w");
780 t_inputrec irInstance;
781 t_inputrec* ir = &irInstance;
783 gmx::TopologyInformation topInfo;
784 topInfo.fillFromInputFile(ftp2fn(efTPR, NFILE, fnm));
785 int ntopatoms = topInfo.mtop()->natoms;
787 bPDB = opt2bSet("-q", NFILE, fnm);
790 snew(xav, ntopatoms);
791 snew(ind_fit, ntopatoms);
792 snew(w_rls, ntopatoms);
793 for (kkk = 0; (kkk < ntopatoms); kkk++)
799 atoms = topInfo.copyAtoms();
801 if (atoms->pdbinfo == nullptr)
803 snew(atoms->pdbinfo, atoms->nr);
805 atoms->havePdbInfo = TRUE;
808 gmx_localtop_t top(topInfo.mtop()->ffparams);
809 gmx_mtop_generate_local_top(*topInfo.mtop(), &top, ir->efep != efepNO);
810 const InteractionDefinitions& idef = top.idef;
813 if (ir->pbcType != PbcType::No)
818 if (ftp2bSet(efNDX, NFILE, fnm))
820 /* TODO: Nothing is written to this file if -c is provided, but it is
822 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
823 xvg = xvgropen(opt2fn("-dr", NFILE, fnm), "Individual Restraints", "Time (ps)", "nm", oenv);
824 snew(vvindex, isize);
826 for (i = 0; (i < isize); i++)
830 sprintf(leg[i], "index %d", index[i]);
832 xvgr_legend(xvg, isize, leg, oenv);
840 t_disresdata disresdata;
844 DisResRunMode::AnalysisTool,
853 int natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
856 init_dr_res(&dr, disresdata.nres);
857 if (opt2bSet("-c", NFILE, fnm))
859 clust = cluster_index(fplog, opt2fn("-c", NFILE, fnm));
860 snew(dr_clust, clust->clust->nr + 1);
861 for (i = 0; (i <= clust->clust->nr); i++)
863 init_dr_res(&dr_clust[i], disresdata.nres);
868 out = xvgropen(opt2fn("-ds", NFILE, fnm), "Sum of Violations", "Time (ps)", "nm", oenv);
869 aver = xvgropen(opt2fn("-da", NFILE, fnm), "Average Violation", "Time (ps)", "nm", oenv);
870 numv = xvgropen(opt2fn("-dn", NFILE, fnm), "# Violations", "Time (ps)", "#", oenv);
871 maxxv = xvgropen(opt2fn("-dm", NFILE, fnm), "Largest Violation", "Time (ps)", "nm", oenv);
874 auto mdAtoms = gmx::makeMDAtoms(fplog, *topInfo.mtop(), *ir, false);
875 atoms2md(topInfo.mtop(), ir, -1, {}, ntopatoms, mdAtoms.get());
876 update_mdatoms(mdAtoms->mdatoms(), ir->fepvals->init_lambda);
877 if (ir->pbcType != PbcType::No)
879 gpbc = gmx_rmpbc_init(idef, ir->pbcType, natoms);
885 if (ir->pbcType != PbcType::No)
887 if (ir->bPeriodicMols)
889 set_pbc(&pbc, ir->pbcType, box);
893 gmx_rmpbc(gpbc, natoms, box, x);
899 if (j > clust->maxframe)
902 "There are more frames in the trajectory than in the cluster index file. "
906 my_clust = clust->inv_clust[j];
907 range_check(my_clust, 0, clust->clust->nr);
909 fplog, idef.il[F_DISRES], idef.iparams, x, f, pbc_null, dr_clust, my_clust, isize, index, vvindex, &disresdata);
913 check_viol(fplog, idef.il[F_DISRES], idef.iparams, x, f, pbc_null, &dr, 0, isize, index, vvindex, &disresdata);
917 reset_x(atoms->nr, ind_fit, atoms->nr, nullptr, x, w_rls);
918 do_fit(atoms->nr, w_rls, x, x);
921 /* Store the first frame of the trajectory as 'characteristic'
922 * for colouring with violations.
924 for (kkk = 0; (kkk < atoms->nr); kkk++)
926 copy_rvec(x[kkk], xav[kkk]);
934 fprintf(xvg, "%10g", t);
935 for (i = 0; (i < isize); i++)
937 fprintf(xvg, " %10g", vvindex[i]);
941 fprintf(out, "%10g %10g\n", t, dr.sumv);
942 fprintf(aver, "%10g %10g\n", t, dr.averv);
943 fprintf(maxxv, "%10g %10g\n", t, dr.maxv);
944 fprintf(numv, "%10g %10d\n", t, dr.nv);
947 } while (read_next_x(oenv, status, &t, x, box));
949 if (ir->pbcType != PbcType::No)
951 gmx_rmpbc_done(gpbc);
957 fplog, disresdata, idef.il[F_DISRES], idef.iparams, clust->clust, dr_clust, clust->grpname, isize, index);
961 dump_stats(fplog, j, disresdata, idef.il[F_DISRES], idef.iparams, &dr, isize, index, bPDB ? atoms.get() : nullptr);
964 write_sto_conf(opt2fn("-q", NFILE, fnm),
965 "Coloured by average violation in Angstrom",
973 opt2fn_null("-x", NFILE, fnm), &dr, disresdata.nres, j, idef, topInfo.mtop(), max_dr, nlevels, bThird);
978 do_view(oenv, opt2fn("-dn", NFILE, fnm), "-nxy");
979 do_view(oenv, opt2fn("-da", NFILE, fnm), "-nxy");
980 do_view(oenv, opt2fn("-ds", NFILE, fnm), "-nxy");
981 do_view(oenv, opt2fn("-dm", NFILE, fnm), "-nxy");
988 do_view(oenv, opt2fn("-dr", NFILE, fnm), "-nxy");