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,2019, 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/mdtypes/fcdata.h"
63 #include "gromacs/mdtypes/inputrec.h"
64 #include "gromacs/mdtypes/md_enums.h"
65 #include "gromacs/pbcutil/ishift.h"
66 #include "gromacs/pbcutil/mshift.h"
67 #include "gromacs/pbcutil/pbc.h"
68 #include "gromacs/pbcutil/rmpbc.h"
69 #include "gromacs/topology/index.h"
70 #include "gromacs/topology/mtop_util.h"
71 #include "gromacs/topology/topology.h"
72 #include "gromacs/trajectoryanalysis/topologyinformation.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 ta_disres(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;
668 rvec *x, *xav = nullptr;
673 int *index = nullptr, *ind_fit = nullptr;
675 t_cluster_ndx *clust = nullptr;
676 t_dr_result dr, *dr_clust = nullptr;
678 real *vvindex = nullptr, *w_rls = nullptr;
679 t_pbc pbc, *pbc_null;
682 gmx_output_env_t *oenv;
683 gmx_rmpbc_t gpbc = nullptr;
686 { efTPR, nullptr, nullptr, ffREAD },
687 { efTRX, "-f", nullptr, ffREAD },
688 { efXVG, "-ds", "drsum", ffWRITE },
689 { efXVG, "-da", "draver", ffWRITE },
690 { efXVG, "-dn", "drnum", ffWRITE },
691 { efXVG, "-dm", "drmax", ffWRITE },
692 { efXVG, "-dr", "restr", ffWRITE },
693 { efLOG, "-l", "disres", ffWRITE },
694 { efNDX, nullptr, "viol", ffOPTRD },
695 { efPDB, "-q", "viol", ffOPTWR },
696 { efNDX, "-c", "clust", ffOPTRD },
697 { efXPM, "-x", "matrix", ffOPTWR }
699 #define NFILE asize(fnm)
701 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
702 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
707 fplog = ftp2FILE(efLOG, NFILE, fnm, "w");
714 t_inputrec irInstance;
715 t_inputrec *ir = &irInstance;
717 gmx::TopologyInformation topInfo;
718 topInfo.fillFromInputFile(ftp2fn(efTPR, NFILE, fnm));
719 int ntopatoms = topInfo.mtop()->natoms;
721 bPDB = opt2bSet("-q", NFILE, fnm);
724 snew(xav, ntopatoms);
725 snew(ind_fit, ntopatoms);
726 snew(w_rls, ntopatoms);
727 for (kkk = 0; (kkk < ntopatoms); kkk++)
733 atoms = topInfo.copyAtoms();
735 if (atoms->pdbinfo == nullptr)
737 snew(atoms->pdbinfo, atoms->nr);
739 atoms->havePdbInfo = TRUE;
742 gmx_mtop_generate_local_top(*topInfo.mtop(), &top, ir->efep != efepNO);
746 if (ir->ePBC != epbcNONE)
748 if (ir->bPeriodicMols)
754 g = mk_graph(fplog, &top.idef, 0, ntopatoms, FALSE, FALSE);
758 if (ftp2bSet(efNDX, NFILE, fnm))
760 /* TODO: Nothing is written to this file if -c is provided, but it is
762 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
763 xvg = xvgropen(opt2fn("-dr", NFILE, fnm), "Individual Restraints", "Time (ps)",
765 snew(vvindex, isize);
767 for (i = 0; (i < isize); i++)
771 sprintf(leg[i], "index %d", index[i]);
773 xvgr_legend(xvg, isize, leg, oenv);
781 init_disres(fplog, topInfo.mtop(), ir, nullptr, nullptr, &fcd, nullptr, FALSE);
783 int natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
786 init_dr_res(&dr, fcd.disres.nres);
787 if (opt2bSet("-c", NFILE, fnm))
789 clust = cluster_index(fplog, opt2fn("-c", NFILE, fnm));
790 snew(dr_clust, clust->clust->nr+1);
791 for (i = 0; (i <= clust->clust->nr); i++)
793 init_dr_res(&dr_clust[i], fcd.disres.nres);
798 out = xvgropen(opt2fn("-ds", NFILE, fnm),
799 "Sum of Violations", "Time (ps)", "nm", oenv);
800 aver = xvgropen(opt2fn("-da", NFILE, fnm),
801 "Average Violation", "Time (ps)", "nm", oenv);
802 numv = xvgropen(opt2fn("-dn", NFILE, fnm),
803 "# Violations", "Time (ps)", "#", oenv);
804 maxxv = xvgropen(opt2fn("-dm", NFILE, fnm),
805 "Largest Violation", "Time (ps)", "nm", oenv);
808 auto mdAtoms = gmx::makeMDAtoms(fplog, *topInfo.mtop(), *ir, false);
809 atoms2md(topInfo.mtop(), ir, -1, nullptr, ntopatoms, mdAtoms.get());
810 update_mdatoms(mdAtoms->mdatoms(), ir->fepvals->init_lambda);
812 if (ir->ePBC != epbcNONE)
814 gpbc = gmx_rmpbc_init(&top.idef, ir->ePBC, natoms);
820 if (ir->ePBC != epbcNONE)
822 if (ir->bPeriodicMols)
824 set_pbc(&pbc, ir->ePBC, box);
828 gmx_rmpbc(gpbc, natoms, box, x);
834 if (j > clust->maxframe)
836 gmx_fatal(FARGS, "There are more frames in the trajectory than in the cluster index file. t = %8f\n", t);
838 my_clust = clust->inv_clust[j];
839 range_check(my_clust, 0, clust->clust->nr);
840 check_viol(fplog, &(top.idef.il[F_DISRES]),
842 x, f, pbc_null, g, dr_clust, my_clust, isize, index, vvindex, &fcd);
846 check_viol(fplog, &(top.idef.il[F_DISRES]),
848 x, f, pbc_null, g, &dr, 0, isize, index, vvindex, &fcd);
852 reset_x(atoms->nr, ind_fit, atoms->nr, nullptr, x, w_rls);
853 do_fit(atoms->nr, w_rls, x, x);
856 /* Store the first frame of the trajectory as 'characteristic'
857 * for colouring with violations.
859 for (kkk = 0; (kkk < atoms->nr); kkk++)
861 copy_rvec(x[kkk], xav[kkk]);
869 fprintf(xvg, "%10g", t);
870 for (i = 0; (i < isize); i++)
872 fprintf(xvg, " %10g", vvindex[i]);
876 fprintf(out, "%10g %10g\n", t, dr.sumv);
877 fprintf(aver, "%10g %10g\n", t, dr.averv);
878 fprintf(maxxv, "%10g %10g\n", t, dr.maxv);
879 fprintf(numv, "%10g %10d\n", t, dr.nv);
883 while (read_next_x(oenv, status, &t, x, box));
885 if (ir->ePBC != epbcNONE)
887 gmx_rmpbc_done(gpbc);
892 dump_clust_stats(fplog, fcd.disres.nres, &(top.idef.il[F_DISRES]),
893 top.idef.iparams, clust->clust, dr_clust,
894 clust->grpname, isize, index);
898 dump_stats(fplog, j, fcd.disres.nres, &(top.idef.il[F_DISRES]),
899 top.idef.iparams, &dr, isize, index,
900 bPDB ? atoms.get() : nullptr);
903 write_sto_conf(opt2fn("-q", NFILE, fnm),
904 "Coloured by average violation in Angstrom",
905 atoms.get(), xav, nullptr, ir->ePBC, box);
907 dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres,
908 j, &top.idef, topInfo.mtop(), max_dr, nlevels, bThird);
913 do_view(oenv, opt2fn("-dn", NFILE, fnm), "-nxy");
914 do_view(oenv, opt2fn("-da", NFILE, fnm), "-nxy");
915 do_view(oenv, opt2fn("-ds", NFILE, fnm), "-nxy");
916 do_view(oenv, opt2fn("-dm", NFILE, fnm), "-nxy");
923 do_view(oenv, opt2fn("-dr", NFILE, fnm), "-nxy");