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, 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/mdrunutility/mdmodules.h"
64 #include "gromacs/mdtypes/fcdata.h"
65 #include "gromacs/mdtypes/inputrec.h"
66 #include "gromacs/mdtypes/md_enums.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/utility/arraysize.h"
75 #include "gromacs/utility/fatalerror.h"
76 #include "gromacs/utility/futil.h"
77 #include "gromacs/utility/smalloc.h"
89 real sumv, averv, maxv;
90 real *aver1, *aver2, *aver_3, *aver_6;
93 static void init5(int n)
99 static void reset5(void)
103 for (i = 0; (i < ntop); i++)
110 int tpcomp(const void *a, const void *b)
118 return static_cast<int>(1e7*(tpb->v-tpa->v));
121 static void add5(int ndr, real viol)
126 for (i = 1; (i < ntop); i++)
128 if (top[i].v < top[mini].v)
133 if (viol > top[mini].v)
140 static void print5(FILE *fp)
144 qsort(top, ntop, sizeof(top[0]), tpcomp);
145 fprintf(fp, "Index:");
146 for (i = 0; (i < ntop); i++)
148 fprintf(fp, " %6d", top[i].n);
150 fprintf(fp, "\nViol: ");
151 for (i = 0; (i < ntop); i++)
153 fprintf(fp, " %6.3f", top[i].v);
158 static void check_viol(FILE *log,
159 t_ilist *disres, t_iparams forceparams[],
161 t_pbc *pbc, t_graph *g, t_dr_result dr[],
162 int clust_id, int isize, int index[], real vvindex[],
166 int i, j, nat, n, type, nviol, ndr, label;
167 real rt, mviol, tviol, viol, lam, dvdl, drt;
169 static gmx_bool bFirst = TRUE;
181 forceatoms = disres->iatoms;
182 for (j = 0; (j < isize); j++)
186 nat = interaction_function[F_DISRES].nratoms+1;
187 for (i = 0; (i < disres->nr); )
189 type = forceatoms[i];
191 label = forceparams[type].disres.label;
194 fprintf(debug, "DISRE: ndr = %d, label = %d i=%d, n =%d\n",
199 gmx_fatal(FARGS, "tpr inconsistency. ndr = %d, label = %d\n", ndr, label);
205 while (((i+n) < disres->nr) &&
206 (forceparams[forceatoms[i+n]].disres.label == label));
208 calc_disres_R_6(NULL, n, &forceatoms[i],
209 (const rvec*)x, pbc, fcd, NULL);
211 if (fcd->disres.Rt_6[0] <= 0)
213 gmx_fatal(FARGS, "ndr = %d, rt_6 = %f", ndr, fcd->disres.Rt_6[0]);
216 rt = gmx::invsixthroot(fcd->disres.Rt_6[0]);
217 dr[clust_id].aver1[ndr] += rt;
218 dr[clust_id].aver2[ndr] += gmx::square(rt);
219 drt = 1.0/gmx::power3(rt);
220 dr[clust_id].aver_3[ndr] += drt;
221 dr[clust_id].aver_6[ndr] += fcd->disres.Rt_6[0];
223 snew(fshift, SHIFTS);
224 interaction_function[F_DISRES].ifunc(n, &forceatoms[i], forceparams,
225 (const rvec*)x, f, fshift,
226 pbc, g, lam, &dvdl, NULL, fcd, NULL);
228 viol = fcd->disres.sumviol;
235 add5(forceparams[type].disres.label, viol);
242 for (j = 0; (j < isize); j++)
244 if (index[j] == forceparams[type].disres.label)
246 vvindex[j] = gmx::invsixthroot(fcd->disres.Rt_6[0]);
253 dr[clust_id].nv = nviol;
254 dr[clust_id].maxv = mviol;
255 dr[clust_id].sumv = tviol;
256 dr[clust_id].averv = tviol/ndr;
257 dr[clust_id].nframes++;
261 fprintf(stderr, "\nThere are %d restraints and %d pairs\n",
262 ndr, disres->nr/nat);
274 real up1, r, rT3, rT6, viol, violT3, violT6;
277 static int drs_comp(const void *a, const void *b)
281 da = (t_dr_stats *)a;
282 db = (t_dr_stats *)b;
284 if (da->viol > db->viol)
288 else if (da->viol < db->viol)
298 static void dump_dump(FILE *log, int ndr, t_dr_stats drs[])
300 static const char *core[] = { "All restraints", "Core restraints" };
301 static const char *tp[] = { "linear", "third power", "sixth power" };
302 real viol_tot, viol_max, viol = 0;
307 for (bCore = FALSE; (bCore <= TRUE); bCore++)
309 for (kkk = 0; (kkk < 3); kkk++)
315 for (i = 0; (i < ndr); i++)
317 if (!bCore || drs[i].bCore)
325 viol = drs[i].violT3;
328 viol = drs[i].violT6;
331 gmx_incons("Dumping violations");
333 viol_max = std::max(viol_max, viol);
342 if ((nrestr > 0) || (bCore && (nrestr < ndr)))
345 fprintf(log, "+++++++ %s ++++++++\n", core[bCore]);
346 fprintf(log, "+++++++ Using %s averaging: ++++++++\n", tp[kkk]);
347 fprintf(log, "Sum of violations: %8.3f nm\n", viol_tot);
350 fprintf(log, "Average violation: %8.3f nm\n", viol_tot/nrestr);
352 fprintf(log, "Largest violation: %8.3f nm\n", viol_max);
353 fprintf(log, "Number of violated restraints: %d/%d\n", nviol, nrestr);
359 static void dump_viol(FILE *log, int ndr, t_dr_stats *drs, gmx_bool bLinear)
363 fprintf(log, "Restr. Core Up1 <r> <rT3> <rT6> <viol><violT3><violT6>\n");
364 for (i = 0; (i < ndr); i++)
366 if (bLinear && (drs[i].viol == 0))
370 fprintf(log, "%6d%5s%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n",
371 drs[i].index, yesno_names[drs[i].bCore],
372 drs[i].up1, drs[i].r, drs[i].rT3, drs[i].rT6,
373 drs[i].viol, drs[i].violT3, drs[i].violT6);
377 static gmx_bool is_core(int i, int isize, int index[])
380 gmx_bool bIC = FALSE;
382 for (kk = 0; !bIC && (kk < isize); kk++)
384 bIC = (index[kk] == i);
390 static void dump_stats(FILE *log, int nsteps, int ndr, t_ilist *disres,
391 t_iparams ip[], t_dr_result *dr,
392 int isize, int index[], t_atoms *atoms)
398 fprintf(log, "++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
401 nra = interaction_function[F_DISRES].nratoms+1;
402 for (i = 0; (i < ndr); i++)
404 /* Search for the appropriate restraint */
405 for (; (j < disres->nr); j += nra)
407 if (ip[disres->iatoms[j]].disres.label == i)
413 drs[i].bCore = is_core(i, isize, index);
414 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
415 drs[i].r = dr->aver1[i]/nsteps;
416 drs[i].rT3 = gmx::invcbrt(dr->aver_3[i]/nsteps);
417 drs[i].rT6 = gmx::invsixthroot(dr->aver_6[i]/nsteps);
418 drs[i].viol = std::max(0.0, static_cast<double>(drs[i].r-drs[i].up1));
419 drs[i].violT3 = std::max(0.0, static_cast<double>(drs[i].rT3-drs[i].up1));
420 drs[i].violT6 = std::max(0.0, static_cast<double>(drs[i].rT6-drs[i].up1));
423 int j1 = disres->iatoms[j+1];
424 int j2 = disres->iatoms[j+2];
425 atoms->pdbinfo[j1].bfac += drs[i].violT3*5;
426 atoms->pdbinfo[j2].bfac += drs[i].violT3*5;
429 dump_viol(log, ndr, drs, FALSE);
431 fprintf(log, "+++ Sorted by linear averaged violations: +++\n");
432 qsort(drs, ndr, sizeof(drs[0]), drs_comp);
433 dump_viol(log, ndr, drs, TRUE);
435 dump_dump(log, ndr, drs);
440 static void dump_clust_stats(FILE *fp, int ndr, t_ilist *disres,
441 t_iparams ip[], t_blocka *clust, t_dr_result dr[],
442 char *clust_name[], int isize, int index[])
444 int i, j, k, nra, mmm = 0;
445 double sumV, maxV, sumVT3, sumVT6, maxVT3, maxVT6;
449 fprintf(fp, "++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
450 fprintf(fp, "Cluster NFrames SumV MaxV SumVT MaxVT SumVS MaxVS\n");
454 for (k = 0; (k < clust->nr); k++)
456 if (dr[k].nframes == 0)
460 if (dr[k].nframes != (clust->index[k+1]-clust->index[k]))
462 gmx_fatal(FARGS, "Inconsistency in cluster %s.\n"
463 "Found %d frames in trajectory rather than the expected %d\n",
464 clust_name[k], dr[k].nframes,
465 clust->index[k+1]-clust->index[k]);
469 gmx_fatal(FARGS, "Inconsistency with cluster %d. Invalid name", k);
472 nra = interaction_function[F_DISRES].nratoms+1;
473 sumV = sumVT3 = sumVT6 = maxV = maxVT3 = maxVT6 = 0;
474 for (i = 0; (i < ndr); i++)
476 /* Search for the appropriate restraint */
477 for (; (j < disres->nr); j += nra)
479 if (ip[disres->iatoms[j]].disres.label == i)
485 drs[i].bCore = is_core(i, isize, index);
486 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
487 drs[i].r = dr[k].aver1[i]/dr[k].nframes;
488 if ((dr[k].aver_3[i] <= 0) || (dr[k].aver_3[i] != dr[k].aver_3[i]))
490 gmx_fatal(FARGS, "dr[%d].aver_3[%d] = %f", k, i, dr[k].aver_3[i]);
492 drs[i].rT3 = gmx::invcbrt(dr[k].aver_3[i]/dr[k].nframes);
493 drs[i].rT6 = gmx::invsixthroot(dr[k].aver_6[i]/dr[k].nframes);
494 drs[i].viol = std::max(0.0, static_cast<double>(drs[i].r-drs[i].up1));
495 drs[i].violT3 = std::max(0.0, static_cast<double>(drs[i].rT3-drs[i].up1));
496 drs[i].violT6 = std::max(0.0, static_cast<double>(drs[i].rT6-drs[i].up1));
498 sumVT3 += drs[i].violT3;
499 sumVT6 += drs[i].violT6;
500 maxV = std::max(maxV, static_cast<double>(drs[i].viol));
501 maxVT3 = std::max(maxVT3, static_cast<double>(drs[i].violT3));
502 maxVT6 = std::max(maxVT6, static_cast<double>(drs[i].violT6));
504 if (std::strcmp(clust_name[k], "1000") == 0)
508 fprintf(fp, "%-10s%6d%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n",
510 dr[k].nframes, sumV, maxV, sumVT3, maxVT3, sumVT6, maxVT6);
517 static void init_dr_res(t_dr_result *dr, int ndr)
519 snew(dr->aver1, ndr+1);
520 snew(dr->aver2, ndr+1);
521 snew(dr->aver_3, ndr+1);
522 snew(dr->aver_6, ndr+1);
530 static void dump_disre_matrix(const char *fn, t_dr_result *dr, int ndr,
531 int nsteps, t_idef *idef, const gmx_mtop_t *mtop,
532 real max_dr, int nlevels, gmx_bool bThird)
536 int n_res, a_offset, mb, mol, a;
538 int i, j, nra, nratoms, tp, ri, rj, index, nlabel, label;
540 real **matrix, *t_res, hi, *w_dr, rav, rviol;
541 t_rgb rlo = { 1, 1, 1 };
542 t_rgb rhi = { 0, 0, 0 };
548 snew(resnr, mtop->natoms);
551 for (mb = 0; mb < mtop->nmolblock; mb++)
553 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
554 for (mol = 0; mol < mtop->molblock[mb].nmol; mol++)
556 for (a = 0; a < atoms->nr; a++)
558 resnr[a_offset+a] = n_res + atoms->atom[a].resind;
560 n_res += atoms->nres;
561 a_offset += atoms->nr;
566 for (i = 0; (i < n_res); i++)
571 for (i = 0; (i < n_res); i++)
573 snew(matrix[i], n_res);
575 nratoms = interaction_function[F_DISRES].nratoms;
576 nra = (idef->il[F_DISRES].nr/(nratoms+1));
582 for (i = 0; (i < idef->il[F_DISRES].nr); i += nratoms+1)
584 tp = idef->il[F_DISRES].iatoms[i];
585 label = idef->iparams[tp].disres.label;
589 /* Set index pointer */
593 gmx_fatal(FARGS, "nlabel is %d, label = %d", nlabel, label);
597 gmx_fatal(FARGS, "ndr = %d, index = %d", ndr, index);
599 /* Update the weight */
600 w_dr[index] = 1.0/nlabel;
609 printf("nlabel = %d, index = %d, ndr = %d\n", nlabel, index, ndr);
611 for (i = 0; (i < ndr); i++)
613 for (j = ptr[i]; (j < ptr[i+1]); j += nratoms+1)
615 tp = idef->il[F_DISRES].iatoms[j];
616 ai = idef->il[F_DISRES].iatoms[j+1];
617 aj = idef->il[F_DISRES].iatoms[j+2];
623 rav = gmx::invcbrt(dr->aver_3[i]/nsteps);
627 rav = dr->aver1[i]/nsteps;
631 fprintf(debug, "DR %d, atoms %d, %d, distance %g\n", i, ai, aj, rav);
633 rviol = std::max(static_cast<real>(0.0), rav-idef->iparams[tp].disres.up1);
634 matrix[ri][rj] += w_dr[i]*rviol;
635 matrix[rj][ri] += w_dr[i]*rviol;
636 hi = std::max(hi, matrix[ri][rj]);
637 hi = std::max(hi, matrix[rj][ri]);
647 printf("Warning: the maxdr that you have specified (%g) is smaller than\nthe largest value in your simulation (%g)\n", max_dr, hi);
651 printf("Highest level in the matrix will be %g\n", hi);
652 fp = gmx_ffopen(fn, "w");
653 write_xpm(fp, 0, "Distance Violations", "<V> (nm)", "Residue", "Residue",
654 n_res, n_res, t_res, t_res, matrix, 0, hi, rlo, rhi, &nlevels);
658 int gmx_disre(int argc, char *argv[])
660 const char *desc[] = {
661 "[THISMODULE] computes violations of distance restraints.",
662 "The program always",
663 "computes the instantaneous violations rather than time-averaged,",
664 "because this analysis is done from a trajectory file afterwards",
665 "it does not make sense to use time averaging. However,",
666 "the time averaged values per restraint are given in the log file.[PAR]",
667 "An index file may be used to select specific restraints for",
669 "When the optional [TT]-q[tt] flag is given a [REF].pdb[ref] file coloured by the",
670 "amount of average violations.[PAR]",
671 "When the [TT]-c[tt] option is given, an index file will be read",
672 "containing the frames in your trajectory corresponding to the clusters",
673 "(defined in another manner) that you want to analyze. For these clusters",
674 "the program will compute average violations using the third power",
675 "averaging algorithm and print them in the log file."
678 static int nlevels = 20;
679 static real max_dr = 0;
680 static gmx_bool bThird = TRUE;
682 { "-ntop", FALSE, etINT, {&ntop},
683 "Number of large violations that are stored in the log file every step" },
684 { "-maxdr", FALSE, etREAL, {&max_dr},
685 "Maximum distance violation in matrix output. If less than or equal to 0 the maximum will be determined by the data." },
686 { "-nlevels", FALSE, etINT, {&nlevels},
687 "Number of levels in the matrix output" },
688 { "-third", FALSE, etBOOL, {&bThird},
689 "Use inverse third power averaging or linear for matrix output" }
692 FILE *out = NULL, *aver = NULL, *numv = NULL, *maxxv = NULL, *xvg = NULL;
697 t_atoms *atoms = NULL;
701 int ntopatoms, natoms, i, j, kkk;
704 rvec *x, *xav = NULL;
709 int *index = NULL, *ind_fit = NULL;
711 t_cluster_ndx *clust = NULL;
712 t_dr_result dr, *dr_clust = NULL;
714 real *vvindex = NULL, *w_rls = NULL;
716 t_pbc pbc, *pbc_null;
719 gmx_output_env_t *oenv;
720 gmx_rmpbc_t gpbc = NULL;
723 { efTPR, NULL, NULL, ffREAD },
724 { efTRX, "-f", NULL, ffREAD },
725 { efXVG, "-ds", "drsum", ffWRITE },
726 { efXVG, "-da", "draver", ffWRITE },
727 { efXVG, "-dn", "drnum", ffWRITE },
728 { efXVG, "-dm", "drmax", ffWRITE },
729 { efXVG, "-dr", "restr", ffWRITE },
730 { efLOG, "-l", "disres", ffWRITE },
731 { efNDX, NULL, "viol", ffOPTRD },
732 { efPDB, "-q", "viol", ffOPTWR },
733 { efNDX, "-c", "clust", ffOPTRD },
734 { efXPM, "-x", "matrix", ffOPTWR }
736 #define NFILE asize(fnm)
738 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
739 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
744 fplog = ftp2FILE(efLOG, NFILE, fnm, "w");
751 gmx::MDModules mdModules;
752 t_inputrec *ir = mdModules.inputrec();
754 read_tpxheader(ftp2fn(efTPR, NFILE, fnm), &header, FALSE);
755 snew(xtop, header.natoms);
756 read_tpx(ftp2fn(efTPR, NFILE, fnm), ir, box, &ntopatoms, xtop, NULL, &mtop);
757 bPDB = opt2bSet("-q", NFILE, fnm);
760 snew(xav, ntopatoms);
761 snew(ind_fit, ntopatoms);
762 snew(w_rls, ntopatoms);
763 for (kkk = 0; (kkk < ntopatoms); kkk++)
770 *atoms = gmx_mtop_global_atoms(&mtop);
772 if (atoms->pdbinfo == NULL)
774 snew(atoms->pdbinfo, atoms->nr);
776 atoms->havePdbInfo = TRUE;
779 top = gmx_mtop_generate_local_top(&mtop, ir->efep != efepNO);
783 if (ir->ePBC != epbcNONE)
785 if (ir->bPeriodicMols)
791 g = mk_graph(fplog, &top->idef, 0, mtop.natoms, FALSE, FALSE);
795 if (ftp2bSet(efNDX, NFILE, fnm))
797 /* TODO: Nothing is written to this file if -c is provided, but it is
799 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
800 xvg = xvgropen(opt2fn("-dr", NFILE, fnm), "Individual Restraints", "Time (ps)",
802 snew(vvindex, isize);
804 for (i = 0; (i < isize); i++)
808 sprintf(leg[i], "index %d", index[i]);
810 xvgr_legend(xvg, isize, (const char**)leg, oenv);
818 init_disres(fplog, &mtop, ir, NULL, &fcd, NULL, FALSE);
820 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
823 init_dr_res(&dr, fcd.disres.nres);
824 if (opt2bSet("-c", NFILE, fnm))
826 clust = cluster_index(fplog, opt2fn("-c", NFILE, fnm));
827 snew(dr_clust, clust->clust->nr+1);
828 for (i = 0; (i <= clust->clust->nr); i++)
830 init_dr_res(&dr_clust[i], fcd.disres.nres);
835 out = xvgropen(opt2fn("-ds", NFILE, fnm),
836 "Sum of Violations", "Time (ps)", "nm", oenv);
837 aver = xvgropen(opt2fn("-da", NFILE, fnm),
838 "Average Violation", "Time (ps)", "nm", oenv);
839 numv = xvgropen(opt2fn("-dn", NFILE, fnm),
840 "# Violations", "Time (ps)", "#", oenv);
841 maxxv = xvgropen(opt2fn("-dm", NFILE, fnm),
842 "Largest Violation", "Time (ps)", "nm", oenv);
845 mdatoms = init_mdatoms(fplog, &mtop, ir->efep != efepNO);
846 atoms2md(&mtop, ir, -1, NULL, mtop.natoms, mdatoms);
847 update_mdatoms(mdatoms, ir->fepvals->init_lambda);
849 if (ir->ePBC != epbcNONE)
851 gpbc = gmx_rmpbc_init(&top->idef, ir->ePBC, natoms);
857 if (ir->ePBC != epbcNONE)
859 if (ir->bPeriodicMols)
861 set_pbc(&pbc, ir->ePBC, box);
865 gmx_rmpbc(gpbc, natoms, box, x);
871 if (j > clust->maxframe)
873 gmx_fatal(FARGS, "There are more frames in the trajectory than in the cluster index file. t = %8f\n", t);
875 my_clust = clust->inv_clust[j];
876 range_check(my_clust, 0, clust->clust->nr);
877 check_viol(fplog, &(top->idef.il[F_DISRES]),
879 x, f, pbc_null, g, dr_clust, my_clust, isize, index, vvindex, &fcd);
883 check_viol(fplog, &(top->idef.il[F_DISRES]),
885 x, f, pbc_null, g, &dr, 0, isize, index, vvindex, &fcd);
889 reset_x(atoms->nr, ind_fit, atoms->nr, NULL, x, w_rls);
890 do_fit(atoms->nr, w_rls, x, x);
893 /* Store the first frame of the trajectory as 'characteristic'
894 * for colouring with violations.
896 for (kkk = 0; (kkk < atoms->nr); kkk++)
898 copy_rvec(x[kkk], xav[kkk]);
906 fprintf(xvg, "%10g", t);
907 for (i = 0; (i < isize); i++)
909 fprintf(xvg, " %10g", vvindex[i]);
913 fprintf(out, "%10g %10g\n", t, dr.sumv);
914 fprintf(aver, "%10g %10g\n", t, dr.averv);
915 fprintf(maxxv, "%10g %10g\n", t, dr.maxv);
916 fprintf(numv, "%10g %10d\n", t, dr.nv);
920 while (read_next_x(oenv, status, &t, x, box));
922 if (ir->ePBC != epbcNONE)
924 gmx_rmpbc_done(gpbc);
929 dump_clust_stats(fplog, fcd.disres.nres, &(top->idef.il[F_DISRES]),
930 top->idef.iparams, clust->clust, dr_clust,
931 clust->grpname, isize, index);
935 dump_stats(fplog, j, fcd.disres.nres, &(top->idef.il[F_DISRES]),
936 top->idef.iparams, &dr, isize, index,
937 bPDB ? atoms : NULL);
940 write_sto_conf(opt2fn("-q", NFILE, fnm),
941 "Coloured by average violation in Angstrom",
942 atoms, xav, NULL, ir->ePBC, box);
944 dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres,
945 j, &top->idef, &mtop, max_dr, nlevels, bThird);
950 do_view(oenv, opt2fn("-dn", NFILE, fnm), "-nxy");
951 do_view(oenv, opt2fn("-da", NFILE, fnm), "-nxy");
952 do_view(oenv, opt2fn("-ds", NFILE, fnm), "-nxy");
953 do_view(oenv, opt2fn("-dm", NFILE, fnm), "-nxy");
960 do_view(oenv, opt2fn("-dr", NFILE, fnm), "-nxy");