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, 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,
151 const InteractionList& disres,
152 gmx::ArrayRef<const t_iparams> forceparams,
164 int i, j, nat, n, type, nviol, ndr, label;
165 real rt, mviol, tviol, viol, lam, dvdl, drt;
167 static gmx_bool bFirst = TRUE;
179 gmx::ArrayRef<const int> forceatoms = disres.iatoms;
180 for (j = 0; (j < isize); j++)
184 nat = interaction_function[F_DISRES].nratoms + 1;
185 // Check internal consistency of disres.label
186 // The label for a distance restraint should be at most one larger
187 // than the previous label.
188 int label_old = forceparams[forceatoms[0]].disres.label;
189 for (i = 0; (i < disres.size()); i += nat)
191 type = forceatoms[i];
192 label = forceparams[type].disres.label;
193 if ((label == label_old) || (label == label_old + 1))
200 "Label mismatch in distance restrains. Label for restraint %d is %d, "
201 "expected it to be either %d or %d",
202 i / nat, label, label_old, label_old + 1);
205 // Get offset for label index
206 label_old = forceparams[forceatoms[0]].disres.label;
207 for (i = 0; (i < disres.size());)
209 type = forceatoms[i];
211 label = forceparams[type].disres.label - label_old;
214 fprintf(debug, "DISRE: ndr = %d, label = %d i=%d, n =%d\n", ndr, label, i, n);
219 } while (((i + n) < disres.size())
220 && (forceparams[forceatoms[i + n]].disres.label == label + label_old));
222 calc_disres_R_6(nullptr, nullptr, n, &forceatoms[i], x, pbc, fcd, nullptr);
224 if (fcd->disres.Rt_6[label] <= 0)
226 gmx_fatal(FARGS, "ndr = %d, rt_6 = %f", ndr, fcd->disres.Rt_6[label]);
229 rt = gmx::invsixthroot(fcd->disres.Rt_6[label]);
230 dr[clust_id].aver1[ndr] += rt;
231 dr[clust_id].aver2[ndr] += gmx::square(rt);
232 drt = 1.0 / gmx::power3(rt);
233 dr[clust_id].aver_3[ndr] += drt;
234 dr[clust_id].aver_6[ndr] += fcd->disres.Rt_6[label];
236 snew(fshift, SHIFTS);
237 ta_disres(n, &forceatoms[i], forceparams.data(), x, f, fshift, pbc, g, lam, &dvdl, nullptr,
240 viol = fcd->disres.sumviol;
247 add5(forceparams[type].disres.label, viol);
254 for (j = 0; (j < isize); j++)
256 if (index[j] == forceparams[type].disres.label)
258 vvindex[j] = gmx::invsixthroot(fcd->disres.Rt_6[label]);
265 dr[clust_id].nv = nviol;
266 dr[clust_id].maxv = mviol;
267 dr[clust_id].sumv = tviol;
268 dr[clust_id].averv = tviol / ndr;
269 dr[clust_id].nframes++;
273 fprintf(stderr, "\nThere are %d restraints and %d pairs\n", ndr, disres.size() / nat);
286 real up1, r, rT3, rT6, viol, violT3, violT6;
289 static void dump_dump(FILE* log, int ndr, t_dr_stats drs[])
291 static const char* core[] = { "All restraints", "Core restraints" };
292 static const char* tp[] = { "linear", "third power", "sixth power" };
293 real viol_tot, viol_max, viol = 0;
298 for (int iCore = 0; iCore < 2; iCore++)
300 bCore = (iCore == 1);
301 for (kkk = 0; (kkk < 3); kkk++)
307 for (i = 0; (i < ndr); i++)
309 if (!bCore || drs[i].bCore)
313 case 0: viol = drs[i].viol; break;
314 case 1: viol = drs[i].violT3; break;
315 case 2: viol = drs[i].violT6; break;
316 default: gmx_incons("Dumping violations");
318 viol_max = std::max(viol_max, viol);
327 if ((nrestr > 0) || (bCore && (nrestr < ndr)))
330 fprintf(log, "+++++++ %s ++++++++\n", core[bCore]);
331 fprintf(log, "+++++++ Using %s averaging: ++++++++\n", tp[kkk]);
332 fprintf(log, "Sum of violations: %8.3f nm\n", viol_tot);
335 fprintf(log, "Average violation: %8.3f nm\n", viol_tot / nrestr);
337 fprintf(log, "Largest violation: %8.3f nm\n", viol_max);
338 fprintf(log, "Number of violated restraints: %d/%d\n", nviol, nrestr);
344 static void dump_viol(FILE* log, int ndr, t_dr_stats* drs, gmx_bool bLinear)
348 fprintf(log, "Restr. Core Up1 <r> <rT3> <rT6> <viol><violT3><violT6>\n");
349 for (i = 0; (i < ndr); i++)
351 if (bLinear && (drs[i].viol == 0))
355 fprintf(log, "%6d%5s%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n", drs[i].label,
356 yesno_names[drs[i].bCore], drs[i].up1, drs[i].r, drs[i].rT3, drs[i].rT6,
357 drs[i].viol, drs[i].violT3, drs[i].violT6);
361 static gmx_bool is_core(int i, int isize, const int index[])
364 gmx_bool bIC = FALSE;
366 for (kk = 0; !bIC && (kk < isize); kk++)
368 bIC = (index[kk] == i);
374 static void dump_stats(FILE* log,
376 const t_disresdata& dd,
377 const InteractionList& disres,
378 gmx::ArrayRef<const t_iparams> ip,
387 fprintf(log, "++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
389 const int nra = interaction_function[F_DISRES].nratoms + 1;
390 for (int j = 0; j < disres.size(); j += nra)
392 // Note that the restraint i can be used by multiple pairs
393 const int i = disres.iatoms[j] - dd.type_min;
394 GMX_RELEASE_ASSERT(i >= 0 && i < dd.nres, "The restraint index should be in range");
396 drs[i].label = ip[disres.iatoms[j]].disres.label;
397 drs[i].bCore = is_core(drs[i].label, isize, index);
398 drs[i].up1 = ip[disres.iatoms[j]].disres.up1;
399 drs[i].r = dr->aver1[i] / nsteps;
400 drs[i].rT3 = gmx::invcbrt(dr->aver_3[i] / nsteps);
401 drs[i].rT6 = gmx::invsixthroot(dr->aver_6[i] / nsteps);
402 drs[i].viol = std::max(0.0, static_cast<double>(drs[i].r - drs[i].up1));
403 drs[i].violT3 = std::max(0.0, static_cast<double>(drs[i].rT3 - drs[i].up1));
404 drs[i].violT6 = std::max(0.0, static_cast<double>(drs[i].rT6 - drs[i].up1));
407 int j1 = disres.iatoms[j + 1];
408 int j2 = disres.iatoms[j + 2];
409 atoms->pdbinfo[j1].bfac += drs[i].violT3 * 5;
410 atoms->pdbinfo[j2].bfac += drs[i].violT3 * 5;
413 dump_viol(log, dd.nres, drs, FALSE);
415 fprintf(log, "+++ Sorted by linear averaged violations: +++\n");
416 std::sort(drs, drs + dd.nres,
417 [](const t_dr_stats& a, const t_dr_stats& b) { return a.viol > b.viol; }); // Reverse sort
418 dump_viol(log, dd.nres, drs, TRUE);
420 dump_dump(log, dd.nres, drs);
425 static void dump_clust_stats(FILE* fp,
426 const t_disresdata& dd,
427 const InteractionList& disres,
428 gmx::ArrayRef<const t_iparams> ip,
436 double sumV, maxV, sumVT3, sumVT6, maxVT3, maxVT6;
440 fprintf(fp, "++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
441 fprintf(fp, "Cluster NFrames SumV MaxV SumVT MaxVT SumVS MaxVS\n");
445 for (k = 0; (k < clust->nr); k++)
447 if (dr[k].nframes == 0)
451 if (dr[k].nframes != (clust->index[k + 1] - clust->index[k]))
454 "Inconsistency in cluster %s.\n"
455 "Found %d frames in trajectory rather than the expected %d\n",
456 clust_name[k], dr[k].nframes, clust->index[k + 1] - clust->index[k]);
460 gmx_fatal(FARGS, "Inconsistency with cluster %d. Invalid name", k);
462 nra = interaction_function[F_DISRES].nratoms + 1;
463 sumV = sumVT3 = sumVT6 = maxV = maxVT3 = maxVT6 = 0;
465 // Use a map to process each restraint only once while looping over all pairs
466 std::unordered_map<int, bool> restraintHasBeenProcessed;
467 for (int j = 0; j < dd.nres; j += nra)
469 // Note that the restraint i can be used by multiple pairs
470 const int i = disres.iatoms[j] - dd.type_min;
472 if (restraintHasBeenProcessed[i])
477 drs[i].label = ip[disres.iatoms[j]].disres.label;
478 drs[i].bCore = is_core(drs[i].label, isize, index);
479 drs[i].up1 = ip[disres.iatoms[j]].disres.up1;
480 drs[i].r = dr[k].aver1[i] / dr[k].nframes;
481 if ((dr[k].aver_3[i] <= 0) || !std::isfinite(dr[k].aver_3[i]))
483 gmx_fatal(FARGS, "dr[%d].aver_3[%d] = %f", k, i, dr[k].aver_3[i]);
485 drs[i].rT3 = gmx::invcbrt(dr[k].aver_3[i] / dr[k].nframes);
486 drs[i].rT6 = gmx::invsixthroot(dr[k].aver_6[i] / dr[k].nframes);
487 drs[i].viol = std::max(0.0, static_cast<double>(drs[i].r - drs[i].up1));
488 drs[i].violT3 = std::max(0.0, static_cast<double>(drs[i].rT3 - drs[i].up1));
489 drs[i].violT6 = std::max(0.0, static_cast<double>(drs[i].rT6 - drs[i].up1));
491 sumVT3 += drs[i].violT3;
492 sumVT6 += drs[i].violT6;
493 maxV = std::max(maxV, static_cast<double>(drs[i].viol));
494 maxVT3 = std::max(maxVT3, static_cast<double>(drs[i].violT3));
495 maxVT6 = std::max(maxVT6, static_cast<double>(drs[i].violT6));
497 // We have processed restraint i, mark it as such
498 restraintHasBeenProcessed[i] = true;
500 if (std::strcmp(clust_name[k], "1000") == 0)
504 fprintf(fp, "%-10s%6d%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n", clust_name[k],
505 dr[k].nframes, sumV, maxV, sumVT3, maxVT3, sumVT6, maxVT6);
511 static void init_dr_res(t_dr_result* dr, int ndr)
513 snew(dr->aver1, ndr + 1);
514 snew(dr->aver2, ndr + 1);
515 snew(dr->aver_3, ndr + 1);
516 snew(dr->aver_6, ndr + 1);
524 static void dump_disre_matrix(const char* fn,
528 const InteractionDefinitions& idef,
529 const gmx_mtop_t* mtop,
536 int n_res, a_offset, mol, a;
537 int i, j, nra, nratoms, tp, ri, rj, index, nlabel, label;
539 real **matrix, *t_res, hi, *w_dr, rav, rviol;
540 t_rgb rlo = { 1, 1, 1 };
541 t_rgb rhi = { 0, 0, 0 };
547 snew(resnr, mtop->natoms);
550 for (const gmx_molblock_t& molb : mtop->molblock)
552 const t_atoms& atoms = mtop->moltype[molb.type].atoms;
553 for (mol = 0; mol < molb.nmol; mol++)
555 for (a = 0; a < atoms.nr; a++)
557 resnr[a_offset + a] = n_res + atoms.atom[a].resind;
560 a_offset += atoms.nr;
565 for (i = 0; (i < n_res); i++)
570 for (i = 0; (i < n_res); i++)
572 snew(matrix[i], n_res);
574 nratoms = interaction_function[F_DISRES].nratoms;
575 nra = (idef.il[F_DISRES].size() / (nratoms + 1));
581 for (i = 0; (i < idef.il[F_DISRES].size()); i += nratoms + 1)
583 tp = idef.il[F_DISRES].iatoms[i];
584 label = idef.iparams[tp].disres.label;
588 /* Set index pointer */
592 gmx_fatal(FARGS, "nlabel is %d, label = %d", nlabel, label);
596 gmx_fatal(FARGS, "ndr = %d, index = %d", ndr, index);
598 /* Update the weight */
599 w_dr[index] = 1.0 / nlabel;
608 printf("nlabel = %d, index = %d, ndr = %d\n", nlabel, index, ndr);
610 for (i = 0; (i < ndr); i++)
612 for (j = ptr[i]; (j < ptr[i + 1]); j += nratoms + 1)
614 tp = idef.il[F_DISRES].iatoms[j];
615 ai = idef.il[F_DISRES].iatoms[j + 1];
616 aj = idef.il[F_DISRES].iatoms[j + 2];
622 rav = gmx::invcbrt(dr->aver_3[i] / nsteps);
626 rav = dr->aver1[i] / nsteps;
630 fprintf(debug, "DR %d, atoms %d, %d, distance %g\n", i, ai, aj, rav);
632 rviol = std::max(0.0_real, rav - idef.iparams[tp].disres.up1);
633 matrix[ri][rj] += w_dr[i] * rviol;
634 matrix[rj][ri] += w_dr[i] * rviol;
635 hi = std::max(hi, matrix[ri][rj]);
636 hi = std::max(hi, matrix[rj][ri]);
646 printf("Warning: the maxdr that you have specified (%g) is smaller than\nthe largest "
647 "value in your simulation (%g)\n",
652 printf("Highest level in the matrix will be %g\n", hi);
653 fp = gmx_ffopen(fn, "w");
654 write_xpm(fp, 0, "Distance Violations", "<V> (nm)", "Residue", "Residue", n_res, n_res, t_res,
655 t_res, matrix, 0, hi, rlo, rhi, &nlevels);
659 int gmx_disre(int argc, char* argv[])
661 const char* desc[] = {
662 "[THISMODULE] computes violations of distance restraints.",
663 "The program always",
664 "computes the instantaneous violations rather than time-averaged,",
665 "because this analysis is done from a trajectory file afterwards",
666 "it does not make sense to use time averaging. However,",
667 "the time averaged values per restraint are given in the log file.[PAR]",
668 "An index file may be used to select specific restraints by index group label for",
670 "When the optional [TT]-q[tt] flag is given a [REF].pdb[ref] file coloured by the",
671 "amount of average violations.[PAR]",
672 "When the [TT]-c[tt] option is given, an index file will be read",
673 "containing the frames in your trajectory corresponding to the clusters",
674 "(defined in another manner) that you want to analyze. For these clusters",
675 "the program will compute average violations using the third power",
676 "averaging algorithm and print them in the log file."
679 static int nlevels = 20;
680 static real max_dr = 0;
681 static gmx_bool bThird = TRUE;
687 "Number of large violations that are stored in the log file every step" },
692 "Maximum distance violation in matrix output. If less than or equal to 0 the "
693 "maximum will be determined by the data." },
694 { "-nlevels", FALSE, etINT, { &nlevels }, "Number of levels in the matrix output" },
699 "Use inverse third power averaging or linear for matrix output" }
702 FILE * out = nullptr, *aver = nullptr, *numv = nullptr, *maxxv = nullptr, *xvg = nullptr;
708 rvec * x, *xav = nullptr;
713 int * index = nullptr, *ind_fit = nullptr;
715 t_cluster_ndx* clust = nullptr;
716 t_dr_result dr, *dr_clust = nullptr;
718 real * vvindex = nullptr, *w_rls = nullptr;
719 t_pbc pbc, *pbc_null;
722 gmx_output_env_t* oenv;
723 gmx_rmpbc_t gpbc = nullptr;
725 t_filenm fnm[] = { { efTPR, nullptr, nullptr, ffREAD }, { efTRX, "-f", nullptr, ffREAD },
726 { efXVG, "-ds", "drsum", ffWRITE }, { efXVG, "-da", "draver", ffWRITE },
727 { efXVG, "-dn", "drnum", ffWRITE }, { efXVG, "-dm", "drmax", ffWRITE },
728 { efXVG, "-dr", "restr", ffWRITE }, { efLOG, "-l", "disres", ffWRITE },
729 { efNDX, nullptr, "viol", ffOPTRD }, { efPDB, "-q", "viol", ffOPTWR },
730 { efNDX, "-c", "clust", ffOPTRD }, { efXPM, "-x", "matrix", ffOPTWR } };
731 #define NFILE asize(fnm)
733 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW, NFILE, fnm, asize(pa), pa,
734 asize(desc), desc, 0, nullptr, &oenv))
739 fplog = ftp2FILE(efLOG, NFILE, fnm, "w");
746 t_inputrec irInstance;
747 t_inputrec* ir = &irInstance;
749 gmx::TopologyInformation topInfo;
750 topInfo.fillFromInputFile(ftp2fn(efTPR, NFILE, fnm));
751 int ntopatoms = topInfo.mtop()->natoms;
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++)
765 atoms = topInfo.copyAtoms();
767 if (atoms->pdbinfo == nullptr)
769 snew(atoms->pdbinfo, atoms->nr);
771 atoms->havePdbInfo = TRUE;
774 gmx_localtop_t top(topInfo.mtop()->ffparams);
775 gmx_mtop_generate_local_top(*topInfo.mtop(), &top, ir->efep != efepNO);
779 if (ir->pbcType != PbcType::No)
781 if (ir->bPeriodicMols)
787 g = mk_graph(fplog, top.idef, 0, ntopatoms, FALSE, FALSE);
791 if (ftp2bSet(efNDX, NFILE, fnm))
793 /* TODO: Nothing is written to this file if -c is provided, but it is
795 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
796 xvg = xvgropen(opt2fn("-dr", NFILE, fnm), "Individual Restraints", "Time (ps)", "nm", oenv);
797 snew(vvindex, isize);
799 for (i = 0; (i < isize); i++)
803 sprintf(leg[i], "index %d", index[i]);
805 xvgr_legend(xvg, isize, leg, oenv);
813 init_disres(fplog, topInfo.mtop(), ir, nullptr, nullptr, &fcd, nullptr, FALSE);
815 int natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
818 init_dr_res(&dr, fcd.disres.nres);
819 if (opt2bSet("-c", NFILE, fnm))
821 clust = cluster_index(fplog, opt2fn("-c", NFILE, fnm));
822 snew(dr_clust, clust->clust->nr + 1);
823 for (i = 0; (i <= clust->clust->nr); i++)
825 init_dr_res(&dr_clust[i], fcd.disres.nres);
830 out = xvgropen(opt2fn("-ds", NFILE, fnm), "Sum of Violations", "Time (ps)", "nm", oenv);
831 aver = xvgropen(opt2fn("-da", NFILE, fnm), "Average Violation", "Time (ps)", "nm", oenv);
832 numv = xvgropen(opt2fn("-dn", NFILE, fnm), "# Violations", "Time (ps)", "#", oenv);
833 maxxv = xvgropen(opt2fn("-dm", NFILE, fnm), "Largest Violation", "Time (ps)", "nm", oenv);
836 auto mdAtoms = gmx::makeMDAtoms(fplog, *topInfo.mtop(), *ir, false);
837 atoms2md(topInfo.mtop(), ir, -1, nullptr, ntopatoms, mdAtoms.get());
838 update_mdatoms(mdAtoms->mdatoms(), ir->fepvals->init_lambda);
839 if (ir->pbcType != PbcType::No)
841 gpbc = gmx_rmpbc_init(top.idef, ir->pbcType, natoms);
847 if (ir->pbcType != PbcType::No)
849 if (ir->bPeriodicMols)
851 set_pbc(&pbc, ir->pbcType, box);
855 gmx_rmpbc(gpbc, natoms, box, x);
861 if (j > clust->maxframe)
864 "There are more frames in the trajectory than in the cluster index file. "
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], top.idef.iparams, x, f, pbc_null, g, dr_clust,
871 my_clust, isize, index, vvindex, &fcd);
875 check_viol(fplog, top.idef.il[F_DISRES], top.idef.iparams, x, f, pbc_null, g, &dr, 0,
876 isize, index, vvindex, &fcd);
880 reset_x(atoms->nr, ind_fit, atoms->nr, nullptr, x, w_rls);
881 do_fit(atoms->nr, w_rls, x, x);
884 /* Store the first frame of the trajectory as 'characteristic'
885 * for colouring with violations.
887 for (kkk = 0; (kkk < atoms->nr); kkk++)
889 copy_rvec(x[kkk], xav[kkk]);
897 fprintf(xvg, "%10g", t);
898 for (i = 0; (i < isize); i++)
900 fprintf(xvg, " %10g", vvindex[i]);
904 fprintf(out, "%10g %10g\n", t, dr.sumv);
905 fprintf(aver, "%10g %10g\n", t, dr.averv);
906 fprintf(maxxv, "%10g %10g\n", t, dr.maxv);
907 fprintf(numv, "%10g %10d\n", t, dr.nv);
910 } while (read_next_x(oenv, status, &t, x, box));
912 if (ir->pbcType != PbcType::No)
914 gmx_rmpbc_done(gpbc);
919 dump_clust_stats(fplog, fcd.disres, top.idef.il[F_DISRES], top.idef.iparams, clust->clust,
920 dr_clust, clust->grpname, isize, index);
924 dump_stats(fplog, j, fcd.disres, top.idef.il[F_DISRES], top.idef.iparams, &dr, isize, index,
925 bPDB ? atoms.get() : nullptr);
928 write_sto_conf(opt2fn("-q", NFILE, fnm), "Coloured by average violation in Angstrom",
929 atoms.get(), xav, nullptr, ir->pbcType, box);
931 dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres, j, top.idef,
932 topInfo.mtop(), max_dr, nlevels, bThird);
937 do_view(oenv, opt2fn("-dn", NFILE, fnm), "-nxy");
938 do_view(oenv, opt2fn("-da", NFILE, fnm), "-nxy");
939 do_view(oenv, opt2fn("-ds", NFILE, fnm), "-nxy");
940 do_view(oenv, opt2fn("-dm", NFILE, fnm), "-nxy");
947 do_view(oenv, opt2fn("-dr", NFILE, fnm), "-nxy");