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"
71 #include "gromacs/pbcutil/rmpbc.h"
72 #include "gromacs/utility/fatalerror.h"
73 #include "gromacs/utility/smalloc.h"
85 real sumv, averv, maxv;
86 real *aver1, *aver2, *aver_3, *aver_6;
89 static void init5(int n)
95 static void reset5(void)
99 for (i = 0; (i < ntop); i++)
106 int tpcomp(const void *a, const void *b)
114 return (1e7*(tpb->v-tpa->v));
117 static void add5(int ndr, real viol)
122 for (i = 1; (i < ntop); i++)
124 if (top[i].v < top[mini].v)
129 if (viol > top[mini].v)
136 static void print5(FILE *fp)
140 qsort(top, ntop, sizeof(top[0]), tpcomp);
141 fprintf(fp, "Index:");
142 for (i = 0; (i < ntop); i++)
144 fprintf(fp, " %6d", top[i].n);
146 fprintf(fp, "\nViol: ");
147 for (i = 0; (i < ntop); i++)
149 fprintf(fp, " %6.3f", top[i].v);
154 static void check_viol(FILE *log,
155 t_ilist *disres, t_iparams forceparams[],
157 t_pbc *pbc, t_graph *g, t_dr_result dr[],
158 int clust_id, int isize, atom_id index[], real vvindex[],
162 int i, j, nat, n, type, nviol, ndr, label;
163 real ener, rt, mviol, tviol, viol, lam, dvdl, drt;
165 static gmx_bool bFirst = TRUE;
177 forceatoms = disres->iatoms;
178 for (j = 0; (j < isize); j++)
182 nat = interaction_function[F_DISRES].nratoms+1;
183 for (i = 0; (i < disres->nr); )
185 type = forceatoms[i];
187 label = forceparams[type].disres.label;
190 fprintf(debug, "DISRE: ndr = %d, label = %d i=%d, n =%d\n",
195 gmx_fatal(FARGS, "tpr inconsistency. ndr = %d, label = %d\n", ndr, label);
201 while (((i+n) < disres->nr) &&
202 (forceparams[forceatoms[i+n]].disres.label == label));
204 calc_disres_R_6(n, &forceatoms[i], forceparams,
205 (const rvec*)x, pbc, fcd, NULL);
207 if (fcd->disres.Rt_6[0] <= 0)
209 gmx_fatal(FARGS, "ndr = %d, rt_6 = %f", ndr, fcd->disres.Rt_6[0]);
212 rt = pow(fcd->disres.Rt_6[0], -1.0/6.0);
213 dr[clust_id].aver1[ndr] += rt;
214 dr[clust_id].aver2[ndr] += sqr(rt);
216 dr[clust_id].aver_3[ndr] += drt;
217 dr[clust_id].aver_6[ndr] += fcd->disres.Rt_6[0];
219 snew(fshift, SHIFTS);
220 ener = interaction_function[F_DISRES].ifunc(n, &forceatoms[i], forceparams,
221 (const rvec*)x, f, fshift,
222 pbc, g, lam, &dvdl, NULL, fcd, NULL);
224 viol = fcd->disres.sumviol;
231 add5(forceparams[type].disres.label, viol);
238 for (j = 0; (j < isize); j++)
240 if (index[j] == forceparams[type].disres.label)
242 vvindex[j] = pow(fcd->disres.Rt_6[0], -1.0/6.0);
249 dr[clust_id].nv = nviol;
250 dr[clust_id].maxv = mviol;
251 dr[clust_id].sumv = tviol;
252 dr[clust_id].averv = tviol/ndr;
253 dr[clust_id].nframes++;
257 fprintf(stderr, "\nThere are %d restraints and %d pairs\n",
258 ndr, disres->nr/nat);
270 real up1, r, rT3, rT6, viol, violT3, violT6;
273 static int drs_comp(const void *a, const void *b)
277 da = (t_dr_stats *)a;
278 db = (t_dr_stats *)b;
280 if (da->viol > db->viol)
284 else if (da->viol < db->viol)
294 static void dump_dump(FILE *log, int ndr, t_dr_stats drs[])
296 static const char *core[] = { "All restraints", "Core restraints" };
297 static const char *tp[] = { "linear", "third power", "sixth power" };
298 real viol_tot, viol_max, viol = 0;
303 for (bCore = FALSE; (bCore <= TRUE); bCore++)
305 for (kkk = 0; (kkk < 3); kkk++)
311 for (i = 0; (i < ndr); i++)
313 if (!bCore || (bCore && drs[i].bCore))
321 viol = drs[i].violT3;
324 viol = drs[i].violT6;
327 gmx_incons("Dumping violations");
329 viol_max = max(viol_max, viol);
338 if ((nrestr > 0) || (bCore && (nrestr < ndr)))
341 fprintf(log, "+++++++ %s ++++++++\n", core[bCore]);
342 fprintf(log, "+++++++ Using %s averaging: ++++++++\n", tp[kkk]);
343 fprintf(log, "Sum of violations: %8.3f nm\n", viol_tot);
346 fprintf(log, "Average violation: %8.3f nm\n", viol_tot/nrestr);
348 fprintf(log, "Largest violation: %8.3f nm\n", viol_max);
349 fprintf(log, "Number of violated restraints: %d/%d\n", nviol, nrestr);
355 static void dump_viol(FILE *log, int ndr, t_dr_stats *drs, gmx_bool bLinear)
359 fprintf(log, "Restr. Core Up1 <r> <rT3> <rT6> <viol><violT3><violT6>\n");
360 for (i = 0; (i < ndr); i++)
362 if (bLinear && (drs[i].viol == 0))
366 fprintf(log, "%6d%5s%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n",
367 drs[i].index, yesno_names[drs[i].bCore],
368 drs[i].up1, drs[i].r, drs[i].rT3, drs[i].rT6,
369 drs[i].viol, drs[i].violT3, drs[i].violT6);
373 static gmx_bool is_core(int i, int isize, atom_id index[])
376 gmx_bool bIC = FALSE;
378 for (kk = 0; !bIC && (kk < isize); kk++)
380 bIC = (index[kk] == i);
386 static void dump_stats(FILE *log, int nsteps, int ndr, t_ilist *disres,
387 t_iparams ip[], t_dr_result *dr,
388 int isize, atom_id index[], t_atoms *atoms)
394 fprintf(log, "++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
397 nra = interaction_function[F_DISRES].nratoms+1;
398 for (i = 0; (i < ndr); i++)
400 /* Search for the appropriate restraint */
401 for (; (j < disres->nr); j += nra)
403 if (ip[disres->iatoms[j]].disres.label == i)
409 drs[i].bCore = is_core(i, isize, index);
410 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
411 drs[i].r = dr->aver1[i]/nsteps;
412 drs[i].rT3 = pow(dr->aver_3[i]/nsteps, -1.0/3.0);
413 drs[i].rT6 = pow(dr->aver_6[i]/nsteps, -1.0/6.0);
414 drs[i].viol = max(0, drs[i].r-drs[i].up1);
415 drs[i].violT3 = max(0, drs[i].rT3-drs[i].up1);
416 drs[i].violT6 = max(0, drs[i].rT6-drs[i].up1);
419 int j1 = disres->iatoms[j+1];
420 int j2 = disres->iatoms[j+2];
421 atoms->pdbinfo[j1].bfac += drs[i].violT3*5;
422 atoms->pdbinfo[j2].bfac += drs[i].violT3*5;
425 dump_viol(log, ndr, drs, FALSE);
427 fprintf(log, "+++ Sorted by linear averaged violations: +++\n");
428 qsort(drs, ndr, sizeof(drs[0]), drs_comp);
429 dump_viol(log, ndr, drs, TRUE);
431 dump_dump(log, ndr, drs);
436 static void dump_clust_stats(FILE *fp, int ndr, t_ilist *disres,
437 t_iparams ip[], t_blocka *clust, t_dr_result dr[],
438 char *clust_name[], int isize, atom_id index[])
440 int i, j, k, nra, mmm = 0;
441 double sumV, maxV, sumVT3, sumVT6, maxVT3, maxVT6;
445 fprintf(fp, "++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
446 fprintf(fp, "Cluster NFrames SumV MaxV SumVT MaxVT SumVS MaxVS\n");
450 for (k = 0; (k < clust->nr); k++)
452 if (dr[k].nframes == 0)
456 if (dr[k].nframes != (clust->index[k+1]-clust->index[k]))
458 gmx_fatal(FARGS, "Inconsistency in cluster %s.\n"
459 "Found %d frames in trajectory rather than the expected %d\n",
460 clust_name[k], dr[k].nframes,
461 clust->index[k+1]-clust->index[k]);
465 gmx_fatal(FARGS, "Inconsistency with cluster %d. Invalid name", k);
468 nra = interaction_function[F_DISRES].nratoms+1;
469 sumV = sumVT3 = sumVT6 = maxV = maxVT3 = maxVT6 = 0;
470 for (i = 0; (i < ndr); i++)
472 /* Search for the appropriate restraint */
473 for (; (j < disres->nr); j += nra)
475 if (ip[disres->iatoms[j]].disres.label == i)
481 drs[i].bCore = is_core(i, isize, index);
482 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
483 drs[i].r = dr[k].aver1[i]/dr[k].nframes;
484 if ((dr[k].aver_3[i] <= 0) || (dr[k].aver_3[i] != dr[k].aver_3[i]))
486 gmx_fatal(FARGS, "dr[%d].aver_3[%d] = %f", k, i, dr[k].aver_3[i]);
488 drs[i].rT3 = pow(dr[k].aver_3[i]/dr[k].nframes, -1.0/3.0);
489 drs[i].rT6 = pow(dr[k].aver_6[i]/dr[k].nframes, -1.0/6.0);
490 drs[i].viol = max(0, drs[i].r-drs[i].up1);
491 drs[i].violT3 = max(0, drs[i].rT3-drs[i].up1);
492 drs[i].violT6 = max(0, drs[i].rT6-drs[i].up1);
494 sumVT3 += drs[i].violT3;
495 sumVT6 += drs[i].violT6;
496 maxV = max(maxV, drs[i].viol);
497 maxVT3 = max(maxVT3, drs[i].violT3);
498 maxVT6 = max(maxVT6, drs[i].violT6);
500 if (strcmp(clust_name[k], "1000") == 0)
504 fprintf(fp, "%-10s%6d%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n",
506 dr[k].nframes, sumV, maxV, sumVT3, maxVT3, sumVT6, maxVT6);
513 static void init_dr_res(t_dr_result *dr, int ndr)
515 snew(dr->aver1, ndr+1);
516 snew(dr->aver2, ndr+1);
517 snew(dr->aver_3, ndr+1);
518 snew(dr->aver_6, ndr+1);
526 static void dump_disre_matrix(const char *fn, t_dr_result *dr, int ndr,
527 int nsteps, t_idef *idef, gmx_mtop_t *mtop,
528 real max_dr, int nlevels, gmx_bool bThird)
532 int n_res, a_offset, mb, mol, a;
534 int iii, i, j, nra, nratoms, tp, ri, rj, index, nlabel, label;
535 atom_id ai, aj, *ptr;
536 real **matrix, *t_res, hi, *w_dr, rav, rviol;
537 t_rgb rlo = { 1, 1, 1 };
538 t_rgb rhi = { 0, 0, 0 };
544 snew(resnr, mtop->natoms);
547 for (mb = 0; mb < mtop->nmolblock; mb++)
549 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
550 for (mol = 0; mol < mtop->molblock[mb].nmol; mol++)
552 for (a = 0; a < atoms->nr; a++)
554 resnr[a_offset+a] = n_res + atoms->atom[a].resind;
556 n_res += atoms->nres;
557 a_offset += atoms->nr;
562 for (i = 0; (i < n_res); i++)
567 for (i = 0; (i < n_res); i++)
569 snew(matrix[i], n_res);
571 nratoms = interaction_function[F_DISRES].nratoms;
572 nra = (idef->il[F_DISRES].nr/(nratoms+1));
578 for (i = 0; (i < idef->il[F_DISRES].nr); i += nratoms+1)
580 tp = idef->il[F_DISRES].iatoms[i];
581 label = idef->iparams[tp].disres.label;
585 /* Set index pointer */
589 gmx_fatal(FARGS, "nlabel is %d, label = %d", nlabel, label);
593 gmx_fatal(FARGS, "ndr = %d, index = %d", ndr, index);
595 /* Update the weight */
596 w_dr[index] = 1.0/nlabel;
605 printf("nlabel = %d, index = %d, ndr = %d\n", nlabel, index, ndr);
607 for (i = 0; (i < ndr); i++)
609 for (j = ptr[i]; (j < ptr[i+1]); j += nratoms+1)
611 tp = idef->il[F_DISRES].iatoms[j];
612 ai = idef->il[F_DISRES].iatoms[j+1];
613 aj = idef->il[F_DISRES].iatoms[j+2];
619 rav = pow(dr->aver_3[i]/nsteps, -1.0/3.0);
623 rav = dr->aver1[i]/nsteps;
627 fprintf(debug, "DR %d, atoms %d, %d, distance %g\n", i, ai, aj, rav);
629 rviol = max(0, rav-idef->iparams[tp].disres.up1);
630 matrix[ri][rj] += w_dr[i]*rviol;
631 matrix[rj][ri] += w_dr[i]*rviol;
632 hi = max(hi, matrix[ri][rj]);
633 hi = max(hi, matrix[rj][ri]);
643 printf("Warning: the maxdr that you have specified (%g) is smaller than\nthe largest value in your simulation (%g)\n", max_dr, hi);
647 printf("Highest level in the matrix will be %g\n", hi);
648 fp = gmx_ffopen(fn, "w");
649 write_xpm(fp, 0, "Distance Violations", "<V> (nm)", "Residue", "Residue",
650 n_res, n_res, t_res, t_res, matrix, 0, hi, rlo, rhi, &nlevels);
654 int gmx_disre(int argc, char *argv[])
656 const char *desc[] = {
657 "[THISMODULE] computes violations of distance restraints.",
658 "If necessary, all protons can be added to a protein molecule ",
659 "using the [gmx-protonate] program.[PAR]",
660 "The program always",
661 "computes the instantaneous violations rather than time-averaged,",
662 "because this analysis is done from a trajectory file afterwards",
663 "it does not make sense to use time averaging. However,",
664 "the time averaged values per restraint are given in the log file.[PAR]",
665 "An index file may be used to select specific restraints for",
667 "When the optional [TT]-q[tt] flag is given a [TT].pdb[tt] file coloured by the",
668 "amount of average violations.[PAR]",
669 "When the [TT]-c[tt] option is given, an index file will be read",
670 "containing the frames in your trajectory corresponding to the clusters",
671 "(defined in another manner) that you want to analyze. For these clusters",
672 "the program will compute average violations using the third power",
673 "averaging algorithm and print them in the log file."
676 static int nlevels = 20;
677 static real max_dr = 0;
678 static gmx_bool bThird = TRUE;
680 { "-ntop", FALSE, etINT, {&ntop},
681 "Number of large violations that are stored in the log file every step" },
682 { "-maxdr", FALSE, etREAL, {&max_dr},
683 "Maximum distance violation in matrix output. If less than or equal to 0 the maximum will be determined by the data." },
684 { "-nlevels", FALSE, etINT, {&nlevels},
685 "Number of levels in the matrix output" },
686 { "-third", FALSE, etBOOL, {&bThird},
687 "Use inverse third power averaging or linear for matrix output" }
690 FILE *out = NULL, *aver = NULL, *numv = NULL, *maxxv = NULL, *xvg = NULL;
696 t_atoms *atoms = NULL;
700 int ntopatoms, natoms, i, j, kkk;
703 rvec *x, *f, *xav = NULL;
707 atom_id *index = NULL, *ind_fit = NULL;
709 t_cluster_ndx *clust = NULL;
710 t_dr_result dr, *dr_clust = NULL;
712 real *vvindex = NULL, *w_rls = NULL;
714 t_pbc pbc, *pbc_null;
718 gmx_rmpbc_t gpbc = NULL;
721 { efTPX, NULL, NULL, ffREAD },
722 { efTRX, "-f", NULL, ffREAD },
723 { efXVG, "-ds", "drsum", ffWRITE },
724 { efXVG, "-da", "draver", ffWRITE },
725 { efXVG, "-dn", "drnum", ffWRITE },
726 { efXVG, "-dm", "drmax", ffWRITE },
727 { efXVG, "-dr", "restr", ffWRITE },
728 { efLOG, "-l", "disres", ffWRITE },
729 { efNDX, NULL, "viol", ffOPTRD },
730 { efPDB, "-q", "viol", ffOPTWR },
731 { efNDX, "-c", "clust", ffOPTRD },
732 { efXPM, "-x", "matrix", ffOPTWR }
734 #define NFILE asize(fnm)
736 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
737 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
742 fplog = ftp2FILE(efLOG, NFILE, fnm, "w");
749 read_tpxheader(ftp2fn(efTPX, NFILE, fnm), &header, FALSE, NULL, NULL);
750 snew(xtop, header.natoms);
751 read_tpx(ftp2fn(efTPX, NFILE, fnm), &ir, box, &ntopatoms, xtop, NULL, NULL, &mtop);
752 bPDB = opt2bSet("-q", NFILE, fnm);
755 snew(xav, ntopatoms);
756 snew(ind_fit, ntopatoms);
757 snew(w_rls, ntopatoms);
758 for (kkk = 0; (kkk < ntopatoms); kkk++)
765 *atoms = gmx_mtop_global_atoms(&mtop);
767 if (atoms->pdbinfo == NULL)
769 snew(atoms->pdbinfo, atoms->nr);
773 top = gmx_mtop_generate_local_top(&mtop, &ir);
777 if (ir.ePBC != epbcNONE)
779 if (ir.bPeriodicMols)
785 g = mk_graph(fplog, &top->idef, 0, mtop.natoms, FALSE, FALSE);
789 if (ftp2bSet(efNDX, NFILE, fnm))
791 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
792 xvg = xvgropen(opt2fn("-dr", NFILE, fnm), "Individual Restraints", "Time (ps)",
794 snew(vvindex, isize);
796 for (i = 0; (i < isize); i++)
800 sprintf(leg[i], "index %d", index[i]);
802 xvgr_legend(xvg, isize, (const char**)leg, oenv);
810 init_disres(fplog, &mtop, &ir, NULL, &fcd, NULL, FALSE);
812 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
815 init_dr_res(&dr, fcd.disres.nres);
816 if (opt2bSet("-c", NFILE, fnm))
818 clust = cluster_index(fplog, opt2fn("-c", NFILE, fnm));
819 snew(dr_clust, clust->clust->nr+1);
820 for (i = 0; (i <= clust->clust->nr); i++)
822 init_dr_res(&dr_clust[i], fcd.disres.nres);
827 out = xvgropen(opt2fn("-ds", NFILE, fnm),
828 "Sum of Violations", "Time (ps)", "nm", oenv);
829 aver = xvgropen(opt2fn("-da", NFILE, fnm),
830 "Average Violation", "Time (ps)", "nm", oenv);
831 numv = xvgropen(opt2fn("-dn", NFILE, fnm),
832 "# Violations", "Time (ps)", "#", oenv);
833 maxxv = xvgropen(opt2fn("-dm", NFILE, fnm),
834 "Largest Violation", "Time (ps)", "nm", oenv);
837 mdatoms = init_mdatoms(fplog, &mtop, ir.efep != efepNO);
838 atoms2md(&mtop, &ir, 0, NULL, mtop.natoms, mdatoms);
839 update_mdatoms(mdatoms, ir.fepvals->init_lambda);
841 if (ir.ePBC != epbcNONE)
843 gpbc = gmx_rmpbc_init(&top->idef, ir.ePBC, natoms);
849 if (ir.ePBC != epbcNONE)
851 if (ir.bPeriodicMols)
853 set_pbc(&pbc, ir.ePBC, box);
857 gmx_rmpbc(gpbc, natoms, box, x);
863 if (j > clust->maxframe)
865 gmx_fatal(FARGS, "There are more frames in the trajectory than in the cluster index file. t = %8f\n", t);
867 my_clust = clust->inv_clust[j];
868 range_check(my_clust, 0, clust->clust->nr);
869 check_viol(fplog, &(top->idef.il[F_DISRES]),
871 x, f, pbc_null, g, dr_clust, my_clust, isize, index, vvindex, &fcd);
875 check_viol(fplog, &(top->idef.il[F_DISRES]),
877 x, f, pbc_null, g, &dr, 0, isize, index, vvindex, &fcd);
881 reset_x(atoms->nr, ind_fit, atoms->nr, NULL, x, w_rls);
882 do_fit(atoms->nr, w_rls, x, x);
885 /* Store the first frame of the trajectory as 'characteristic'
886 * for colouring with violations.
888 for (kkk = 0; (kkk < atoms->nr); kkk++)
890 copy_rvec(x[kkk], xav[kkk]);
898 fprintf(xvg, "%10g", t);
899 for (i = 0; (i < isize); i++)
901 fprintf(xvg, " %10g", vvindex[i]);
905 fprintf(out, "%10g %10g\n", t, dr.sumv);
906 fprintf(aver, "%10g %10g\n", t, dr.averv);
907 fprintf(maxxv, "%10g %10g\n", t, dr.maxv);
908 fprintf(numv, "%10g %10d\n", t, dr.nv);
912 while (read_next_x(oenv, status, &t, x, box));
914 if (ir.ePBC != epbcNONE)
916 gmx_rmpbc_done(gpbc);
921 dump_clust_stats(fplog, fcd.disres.nres, &(top->idef.il[F_DISRES]),
922 top->idef.iparams, clust->clust, dr_clust,
923 clust->grpname, isize, index);
927 dump_stats(fplog, j, fcd.disres.nres, &(top->idef.il[F_DISRES]),
928 top->idef.iparams, &dr, isize, index,
929 bPDB ? atoms : NULL);
932 write_sto_conf(opt2fn("-q", NFILE, fnm),
933 "Coloured by average violation in Angstrom",
934 atoms, xav, NULL, ir.ePBC, box);
936 dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres,
937 j, &top->idef, &mtop, max_dr, nlevels, bThird);
945 do_view(oenv, opt2fn("-dr", NFILE, fnm), "-nxy");
947 do_view(oenv, opt2fn("-dn", NFILE, fnm), "-nxy");
948 do_view(oenv, opt2fn("-da", NFILE, fnm), "-nxy");
949 do_view(oenv, opt2fn("-ds", NFILE, fnm), "-nxy");
950 do_view(oenv, opt2fn("-dm", NFILE, fnm), "-nxy");
953 gmx_log_close(fplog);