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, 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.
47 #include "gromacs/pbcutil/mshift.h"
49 #include "gromacs/fileio/confio.h"
55 #include "gromacs/fileio/pdbio.h"
58 #include "gromacs/fileio/tpxio.h"
59 #include "gromacs/fileio/trxio.h"
62 #include "mtop_util.h"
65 #include "gromacs/commandline/pargs.h"
66 #include "gromacs/fileio/matio.h"
67 #include "gromacs/fileio/xvgr.h"
68 #include "gromacs/math/do_fit.h"
69 #include "gromacs/math/vec.h"
70 #include "gromacs/pbcutil/rmpbc.h"
71 #include "gromacs/utility/fatalerror.h"
72 #include "gromacs/utility/smalloc.h"
84 real sumv, averv, maxv;
85 real *aver1, *aver2, *aver_3, *aver_6;
88 static void init5(int n)
94 static void reset5(void)
98 for (i = 0; (i < ntop); i++)
105 int tpcomp(const void *a, const void *b)
113 return (1e7*(tpb->v-tpa->v));
116 static void add5(int ndr, real viol)
121 for (i = 1; (i < ntop); i++)
123 if (top[i].v < top[mini].v)
128 if (viol > top[mini].v)
135 static void print5(FILE *fp)
139 qsort(top, ntop, sizeof(top[0]), tpcomp);
140 fprintf(fp, "Index:");
141 for (i = 0; (i < ntop); i++)
143 fprintf(fp, " %6d", top[i].n);
145 fprintf(fp, "\nViol: ");
146 for (i = 0; (i < ntop); i++)
148 fprintf(fp, " %6.3f", top[i].v);
153 static void check_viol(FILE *log,
154 t_ilist *disres, t_iparams forceparams[],
156 t_pbc *pbc, t_graph *g, t_dr_result dr[],
157 int clust_id, int isize, atom_id index[], real vvindex[],
161 int i, j, nat, n, type, nviol, ndr, label;
162 real ener, rt, mviol, tviol, viol, lam, dvdl, drt;
164 static gmx_bool bFirst = TRUE;
176 forceatoms = disres->iatoms;
177 for (j = 0; (j < isize); j++)
181 nat = interaction_function[F_DISRES].nratoms+1;
182 for (i = 0; (i < disres->nr); )
184 type = forceatoms[i];
186 label = forceparams[type].disres.label;
189 fprintf(debug, "DISRE: ndr = %d, label = %d i=%d, n =%d\n",
194 gmx_fatal(FARGS, "tpr inconsistency. ndr = %d, label = %d\n", ndr, label);
200 while (((i+n) < disres->nr) &&
201 (forceparams[forceatoms[i+n]].disres.label == label));
203 calc_disres_R_6(n, &forceatoms[i], forceparams,
204 (const rvec*)x, pbc, fcd, NULL);
206 if (fcd->disres.Rt_6[0] <= 0)
208 gmx_fatal(FARGS, "ndr = %d, rt_6 = %f", ndr, fcd->disres.Rt_6[0]);
211 rt = pow(fcd->disres.Rt_6[0], -1.0/6.0);
212 dr[clust_id].aver1[ndr] += rt;
213 dr[clust_id].aver2[ndr] += sqr(rt);
215 dr[clust_id].aver_3[ndr] += drt;
216 dr[clust_id].aver_6[ndr] += fcd->disres.Rt_6[0];
218 snew(fshift, SHIFTS);
219 ener = interaction_function[F_DISRES].ifunc(n, &forceatoms[i], forceparams,
220 (const rvec*)x, f, fshift,
221 pbc, g, lam, &dvdl, NULL, fcd, NULL);
223 viol = fcd->disres.sumviol;
230 add5(forceparams[type].disres.label, viol);
237 for (j = 0; (j < isize); j++)
239 if (index[j] == forceparams[type].disres.label)
241 vvindex[j] = pow(fcd->disres.Rt_6[0], -1.0/6.0);
248 dr[clust_id].nv = nviol;
249 dr[clust_id].maxv = mviol;
250 dr[clust_id].sumv = tviol;
251 dr[clust_id].averv = tviol/ndr;
252 dr[clust_id].nframes++;
256 fprintf(stderr, "\nThere are %d restraints and %d pairs\n",
257 ndr, disres->nr/nat);
269 real up1, r, rT3, rT6, viol, violT3, violT6;
272 static int drs_comp(const void *a, const void *b)
276 da = (t_dr_stats *)a;
277 db = (t_dr_stats *)b;
279 if (da->viol > db->viol)
283 else if (da->viol < db->viol)
293 static void dump_dump(FILE *log, int ndr, t_dr_stats drs[])
295 static const char *core[] = { "All restraints", "Core restraints" };
296 static const char *tp[] = { "linear", "third power", "sixth power" };
297 real viol_tot, viol_max, viol = 0;
302 for (bCore = FALSE; (bCore <= TRUE); bCore++)
304 for (kkk = 0; (kkk < 3); kkk++)
310 for (i = 0; (i < ndr); i++)
312 if (!bCore || (bCore && drs[i].bCore))
320 viol = drs[i].violT3;
323 viol = drs[i].violT6;
326 gmx_incons("Dumping violations");
328 viol_max = max(viol_max, viol);
337 if ((nrestr > 0) || (bCore && (nrestr < ndr)))
340 fprintf(log, "+++++++ %s ++++++++\n", core[bCore]);
341 fprintf(log, "+++++++ Using %s averaging: ++++++++\n", tp[kkk]);
342 fprintf(log, "Sum of violations: %8.3f nm\n", viol_tot);
345 fprintf(log, "Average violation: %8.3f nm\n", viol_tot/nrestr);
347 fprintf(log, "Largest violation: %8.3f nm\n", viol_max);
348 fprintf(log, "Number of violated restraints: %d/%d\n", nviol, nrestr);
354 static void dump_viol(FILE *log, int ndr, t_dr_stats *drs, gmx_bool bLinear)
358 fprintf(log, "Restr. Core Up1 <r> <rT3> <rT6> <viol><violT3><violT6>\n");
359 for (i = 0; (i < ndr); i++)
361 if (bLinear && (drs[i].viol == 0))
365 fprintf(log, "%6d%5s%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n",
366 drs[i].index, yesno_names[drs[i].bCore],
367 drs[i].up1, drs[i].r, drs[i].rT3, drs[i].rT6,
368 drs[i].viol, drs[i].violT3, drs[i].violT6);
372 static gmx_bool is_core(int i, int isize, atom_id index[])
375 gmx_bool bIC = FALSE;
377 for (kk = 0; !bIC && (kk < isize); kk++)
379 bIC = (index[kk] == i);
385 static void dump_stats(FILE *log, int nsteps, int ndr, t_ilist *disres,
386 t_iparams ip[], t_dr_result *dr,
387 int isize, atom_id index[], t_atoms *atoms)
393 fprintf(log, "++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
396 nra = interaction_function[F_DISRES].nratoms+1;
397 for (i = 0; (i < ndr); i++)
399 /* Search for the appropriate restraint */
400 for (; (j < disres->nr); j += nra)
402 if (ip[disres->iatoms[j]].disres.label == i)
408 drs[i].bCore = is_core(i, isize, index);
409 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
410 drs[i].r = dr->aver1[i]/nsteps;
411 drs[i].rT3 = pow(dr->aver_3[i]/nsteps, -1.0/3.0);
412 drs[i].rT6 = pow(dr->aver_6[i]/nsteps, -1.0/6.0);
413 drs[i].viol = max(0, drs[i].r-drs[i].up1);
414 drs[i].violT3 = max(0, drs[i].rT3-drs[i].up1);
415 drs[i].violT6 = max(0, drs[i].rT6-drs[i].up1);
418 int j1 = disres->iatoms[j+1];
419 int j2 = disres->iatoms[j+2];
420 atoms->pdbinfo[j1].bfac += drs[i].violT3*5;
421 atoms->pdbinfo[j2].bfac += drs[i].violT3*5;
424 dump_viol(log, ndr, drs, FALSE);
426 fprintf(log, "+++ Sorted by linear averaged violations: +++\n");
427 qsort(drs, ndr, sizeof(drs[0]), drs_comp);
428 dump_viol(log, ndr, drs, TRUE);
430 dump_dump(log, ndr, drs);
435 static void dump_clust_stats(FILE *fp, int ndr, t_ilist *disres,
436 t_iparams ip[], t_blocka *clust, t_dr_result dr[],
437 char *clust_name[], int isize, atom_id index[])
439 int i, j, k, nra, mmm = 0;
440 double sumV, maxV, sumVT3, sumVT6, maxVT3, maxVT6;
444 fprintf(fp, "++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
445 fprintf(fp, "Cluster NFrames SumV MaxV SumVT MaxVT SumVS MaxVS\n");
449 for (k = 0; (k < clust->nr); k++)
451 if (dr[k].nframes == 0)
455 if (dr[k].nframes != (clust->index[k+1]-clust->index[k]))
457 gmx_fatal(FARGS, "Inconsistency in cluster %s.\n"
458 "Found %d frames in trajectory rather than the expected %d\n",
459 clust_name[k], dr[k].nframes,
460 clust->index[k+1]-clust->index[k]);
464 gmx_fatal(FARGS, "Inconsistency with cluster %d. Invalid name", k);
467 nra = interaction_function[F_DISRES].nratoms+1;
468 sumV = sumVT3 = sumVT6 = maxV = maxVT3 = maxVT6 = 0;
469 for (i = 0; (i < ndr); i++)
471 /* Search for the appropriate restraint */
472 for (; (j < disres->nr); j += nra)
474 if (ip[disres->iatoms[j]].disres.label == i)
480 drs[i].bCore = is_core(i, isize, index);
481 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
482 drs[i].r = dr[k].aver1[i]/dr[k].nframes;
483 if ((dr[k].aver_3[i] <= 0) || (dr[k].aver_3[i] != dr[k].aver_3[i]))
485 gmx_fatal(FARGS, "dr[%d].aver_3[%d] = %f", k, i, dr[k].aver_3[i]);
487 drs[i].rT3 = pow(dr[k].aver_3[i]/dr[k].nframes, -1.0/3.0);
488 drs[i].rT6 = pow(dr[k].aver_6[i]/dr[k].nframes, -1.0/6.0);
489 drs[i].viol = max(0, drs[i].r-drs[i].up1);
490 drs[i].violT3 = max(0, drs[i].rT3-drs[i].up1);
491 drs[i].violT6 = max(0, drs[i].rT6-drs[i].up1);
493 sumVT3 += drs[i].violT3;
494 sumVT6 += drs[i].violT6;
495 maxV = max(maxV, drs[i].viol);
496 maxVT3 = max(maxVT3, drs[i].violT3);
497 maxVT6 = max(maxVT6, drs[i].violT6);
499 if (strcmp(clust_name[k], "1000") == 0)
503 fprintf(fp, "%-10s%6d%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n",
505 dr[k].nframes, sumV, maxV, sumVT3, maxVT3, sumVT6, maxVT6);
512 static void init_dr_res(t_dr_result *dr, int ndr)
514 snew(dr->aver1, ndr+1);
515 snew(dr->aver2, ndr+1);
516 snew(dr->aver_3, ndr+1);
517 snew(dr->aver_6, ndr+1);
525 static void dump_disre_matrix(const char *fn, t_dr_result *dr, int ndr,
526 int nsteps, t_idef *idef, gmx_mtop_t *mtop,
527 real max_dr, int nlevels, gmx_bool bThird)
531 int n_res, a_offset, mb, mol, a;
533 int iii, i, j, nra, nratoms, tp, ri, rj, index, nlabel, label;
534 atom_id ai, aj, *ptr;
535 real **matrix, *t_res, hi, *w_dr, rav, rviol;
536 t_rgb rlo = { 1, 1, 1 };
537 t_rgb rhi = { 0, 0, 0 };
543 snew(resnr, mtop->natoms);
546 for (mb = 0; mb < mtop->nmolblock; mb++)
548 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
549 for (mol = 0; mol < mtop->molblock[mb].nmol; mol++)
551 for (a = 0; a < atoms->nr; a++)
553 resnr[a_offset+a] = n_res + atoms->atom[a].resind;
555 n_res += atoms->nres;
556 a_offset += atoms->nr;
561 for (i = 0; (i < n_res); i++)
566 for (i = 0; (i < n_res); i++)
568 snew(matrix[i], n_res);
570 nratoms = interaction_function[F_DISRES].nratoms;
571 nra = (idef->il[F_DISRES].nr/(nratoms+1));
577 for (i = 0; (i < idef->il[F_DISRES].nr); i += nratoms+1)
579 tp = idef->il[F_DISRES].iatoms[i];
580 label = idef->iparams[tp].disres.label;
584 /* Set index pointer */
588 gmx_fatal(FARGS, "nlabel is %d, label = %d", nlabel, label);
592 gmx_fatal(FARGS, "ndr = %d, index = %d", ndr, index);
594 /* Update the weight */
595 w_dr[index] = 1.0/nlabel;
604 printf("nlabel = %d, index = %d, ndr = %d\n", nlabel, index, ndr);
606 for (i = 0; (i < ndr); i++)
608 for (j = ptr[i]; (j < ptr[i+1]); j += nratoms+1)
610 tp = idef->il[F_DISRES].iatoms[j];
611 ai = idef->il[F_DISRES].iatoms[j+1];
612 aj = idef->il[F_DISRES].iatoms[j+2];
618 rav = pow(dr->aver_3[i]/nsteps, -1.0/3.0);
622 rav = dr->aver1[i]/nsteps;
626 fprintf(debug, "DR %d, atoms %d, %d, distance %g\n", i, ai, aj, rav);
628 rviol = max(0, rav-idef->iparams[tp].disres.up1);
629 matrix[ri][rj] += w_dr[i]*rviol;
630 matrix[rj][ri] += w_dr[i]*rviol;
631 hi = max(hi, matrix[ri][rj]);
632 hi = max(hi, matrix[rj][ri]);
642 printf("Warning: the maxdr that you have specified (%g) is smaller than\nthe largest value in your simulation (%g)\n", max_dr, hi);
646 printf("Highest level in the matrix will be %g\n", hi);
647 fp = gmx_ffopen(fn, "w");
648 write_xpm(fp, 0, "Distance Violations", "<V> (nm)", "Residue", "Residue",
649 n_res, n_res, t_res, t_res, matrix, 0, hi, rlo, rhi, &nlevels);
653 int gmx_disre(int argc, char *argv[])
655 const char *desc[] = {
656 "[THISMODULE] computes violations of distance restraints.",
657 "If necessary, all protons can be added to a protein molecule ",
658 "using the [gmx-protonate] program.[PAR]",
659 "The program always",
660 "computes the instantaneous violations rather than time-averaged,",
661 "because this analysis is done from a trajectory file afterwards",
662 "it does not make sense to use time averaging. However,",
663 "the time averaged values per restraint are given in the log file.[PAR]",
664 "An index file may be used to select specific restraints for",
666 "When the optional [TT]-q[tt] flag is given a [TT].pdb[tt] file coloured by the",
667 "amount of average violations.[PAR]",
668 "When the [TT]-c[tt] option is given, an index file will be read",
669 "containing the frames in your trajectory corresponding to the clusters",
670 "(defined in another manner) that you want to analyze. For these clusters",
671 "the program will compute average violations using the third power",
672 "averaging algorithm and print them in the log file."
675 static int nlevels = 20;
676 static real max_dr = 0;
677 static gmx_bool bThird = TRUE;
679 { "-ntop", FALSE, etINT, {&ntop},
680 "Number of large violations that are stored in the log file every step" },
681 { "-maxdr", FALSE, etREAL, {&max_dr},
682 "Maximum distance violation in matrix output. If less than or equal to 0 the maximum will be determined by the data." },
683 { "-nlevels", FALSE, etINT, {&nlevels},
684 "Number of levels in the matrix output" },
685 { "-third", FALSE, etBOOL, {&bThird},
686 "Use inverse third power averaging or linear for matrix output" }
689 FILE *out = NULL, *aver = NULL, *numv = NULL, *maxxv = NULL, *xvg = NULL;
695 t_atoms *atoms = NULL;
699 int ntopatoms, natoms, i, j, kkk;
702 rvec *x, *f, *xav = NULL;
706 atom_id *index = NULL, *ind_fit = NULL;
708 t_cluster_ndx *clust = NULL;
709 t_dr_result dr, *dr_clust = NULL;
711 real *vvindex = NULL, *w_rls = NULL;
713 t_pbc pbc, *pbc_null;
717 gmx_rmpbc_t gpbc = NULL;
720 { efTPX, NULL, NULL, ffREAD },
721 { efTRX, "-f", NULL, ffREAD },
722 { efXVG, "-ds", "drsum", ffWRITE },
723 { efXVG, "-da", "draver", ffWRITE },
724 { efXVG, "-dn", "drnum", ffWRITE },
725 { efXVG, "-dm", "drmax", ffWRITE },
726 { efXVG, "-dr", "restr", ffWRITE },
727 { efLOG, "-l", "disres", ffWRITE },
728 { efNDX, NULL, "viol", ffOPTRD },
729 { efPDB, "-q", "viol", ffOPTWR },
730 { efNDX, "-c", "clust", ffOPTRD },
731 { efXPM, "-x", "matrix", ffOPTWR }
733 #define NFILE asize(fnm)
735 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
736 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
741 fplog = ftp2FILE(efLOG, NFILE, fnm, "w");
748 read_tpxheader(ftp2fn(efTPX, NFILE, fnm), &header, FALSE, NULL, NULL);
749 snew(xtop, header.natoms);
750 read_tpx(ftp2fn(efTPX, NFILE, fnm), &ir, box, &ntopatoms, xtop, NULL, NULL, &mtop);
751 bPDB = opt2bSet("-q", NFILE, fnm);
754 snew(xav, ntopatoms);
755 snew(ind_fit, ntopatoms);
756 snew(w_rls, ntopatoms);
757 for (kkk = 0; (kkk < ntopatoms); kkk++)
764 *atoms = gmx_mtop_global_atoms(&mtop);
766 if (atoms->pdbinfo == NULL)
768 snew(atoms->pdbinfo, atoms->nr);
772 top = gmx_mtop_generate_local_top(&mtop, &ir);
776 if (ir.ePBC != epbcNONE)
778 if (ir.bPeriodicMols)
784 g = mk_graph(fplog, &top->idef, 0, mtop.natoms, FALSE, FALSE);
788 if (ftp2bSet(efNDX, NFILE, fnm))
790 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
791 xvg = xvgropen(opt2fn("-dr", NFILE, fnm), "Individual Restraints", "Time (ps)",
793 snew(vvindex, isize);
795 for (i = 0; (i < isize); i++)
799 sprintf(leg[i], "index %d", index[i]);
801 xvgr_legend(xvg, isize, (const char**)leg, oenv);
809 init_disres(fplog, &mtop, &ir, NULL, &fcd, NULL, FALSE);
811 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
814 init_dr_res(&dr, fcd.disres.nres);
815 if (opt2bSet("-c", NFILE, fnm))
817 clust = cluster_index(fplog, opt2fn("-c", NFILE, fnm));
818 snew(dr_clust, clust->clust->nr+1);
819 for (i = 0; (i <= clust->clust->nr); i++)
821 init_dr_res(&dr_clust[i], fcd.disres.nres);
826 out = xvgropen(opt2fn("-ds", NFILE, fnm),
827 "Sum of Violations", "Time (ps)", "nm", oenv);
828 aver = xvgropen(opt2fn("-da", NFILE, fnm),
829 "Average Violation", "Time (ps)", "nm", oenv);
830 numv = xvgropen(opt2fn("-dn", NFILE, fnm),
831 "# Violations", "Time (ps)", "#", oenv);
832 maxxv = xvgropen(opt2fn("-dm", NFILE, fnm),
833 "Largest Violation", "Time (ps)", "nm", oenv);
836 mdatoms = init_mdatoms(fplog, &mtop, ir.efep != efepNO);
837 atoms2md(&mtop, &ir, 0, NULL, mtop.natoms, mdatoms);
838 update_mdatoms(mdatoms, ir.fepvals->init_lambda);
840 if (ir.ePBC != epbcNONE)
842 gpbc = gmx_rmpbc_init(&top->idef, ir.ePBC, natoms);
848 if (ir.ePBC != epbcNONE)
850 if (ir.bPeriodicMols)
852 set_pbc(&pbc, ir.ePBC, box);
856 gmx_rmpbc(gpbc, natoms, box, x);
862 if (j > clust->maxframe)
864 gmx_fatal(FARGS, "There are more frames in the trajectory than in the cluster index file. t = %8f\n", t);
866 my_clust = clust->inv_clust[j];
867 range_check(my_clust, 0, clust->clust->nr);
868 check_viol(fplog, &(top->idef.il[F_DISRES]),
870 x, f, pbc_null, g, dr_clust, my_clust, isize, index, vvindex, &fcd);
874 check_viol(fplog, &(top->idef.il[F_DISRES]),
876 x, f, pbc_null, g, &dr, 0, isize, index, vvindex, &fcd);
880 reset_x(atoms->nr, ind_fit, atoms->nr, NULL, 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);
911 while (read_next_x(oenv, status, &t, x, box));
913 if (ir.ePBC != epbcNONE)
915 gmx_rmpbc_done(gpbc);
920 dump_clust_stats(fplog, fcd.disres.nres, &(top->idef.il[F_DISRES]),
921 top->idef.iparams, clust->clust, dr_clust,
922 clust->grpname, isize, index);
926 dump_stats(fplog, j, fcd.disres.nres, &(top->idef.il[F_DISRES]),
927 top->idef.iparams, &dr, isize, index,
928 bPDB ? atoms : NULL);
931 write_sto_conf(opt2fn("-q", NFILE, fnm),
932 "Coloured by average violation in Angstrom",
933 atoms, xav, NULL, ir.ePBC, box);
935 dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres,
936 j, &top->idef, &mtop, max_dr, nlevels, bThird);
944 do_view(oenv, opt2fn("-dr", NFILE, fnm), "-nxy");
946 do_view(oenv, opt2fn("-dn", NFILE, fnm), "-nxy");
947 do_view(oenv, opt2fn("-da", NFILE, fnm), "-nxy");
948 do_view(oenv, opt2fn("-ds", NFILE, fnm), "-nxy");
949 do_view(oenv, opt2fn("-dm", NFILE, fnm), "-nxy");
952 gmx_log_close(fplog);