3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
62 #include "mtop_util.h"
76 real sumv, averv, maxv;
77 real *aver1, *aver2, *aver_3, *aver_6;
80 static void init5(int n)
86 static void reset5(void)
90 for (i = 0; (i < ntop); i++)
97 int tpcomp(const void *a, const void *b)
105 return (1e7*(tpb->v-tpa->v));
108 static void add5(int ndr, real viol)
113 for (i = 1; (i < ntop); i++)
115 if (top[i].v < top[mini].v)
120 if (viol > top[mini].v)
127 static void print5(FILE *fp)
131 qsort(top, ntop, sizeof(top[0]), tpcomp);
132 fprintf(fp, "Index:");
133 for (i = 0; (i < ntop); i++)
135 fprintf(fp, " %6d", top[i].n);
137 fprintf(fp, "\nViol: ");
138 for (i = 0; (i < ntop); i++)
140 fprintf(fp, " %6.3f", top[i].v);
145 static void check_viol(FILE *log,
146 t_ilist *disres, t_iparams forceparams[],
148 t_pbc *pbc, t_graph *g, t_dr_result dr[],
149 int clust_id, int isize, atom_id index[], real vvindex[],
153 int i, j, nat, n, type, nviol, ndr, label;
154 real ener, rt, mviol, tviol, viol, lam, dvdl, drt;
156 static gmx_bool bFirst = TRUE;
168 forceatoms = disres->iatoms;
169 for (j = 0; (j < isize); j++)
173 nat = interaction_function[F_DISRES].nratoms+1;
174 for (i = 0; (i < disres->nr); )
176 type = forceatoms[i];
178 label = forceparams[type].disres.label;
181 fprintf(debug, "DISRE: ndr = %d, label = %d i=%d, n =%d\n",
186 gmx_fatal(FARGS, "tpr inconsistency. ndr = %d, label = %d\n", ndr, label);
192 while (((i+n) < disres->nr) &&
193 (forceparams[forceatoms[i+n]].disres.label == label));
195 calc_disres_R_6(n, &forceatoms[i], forceparams,
196 (const rvec*)x, pbc, fcd, NULL);
198 if (fcd->disres.Rt_6[0] <= 0)
200 gmx_fatal(FARGS, "ndr = %d, rt_6 = %f", ndr, fcd->disres.Rt_6[0]);
203 rt = pow(fcd->disres.Rt_6[0], -1.0/6.0);
204 dr[clust_id].aver1[ndr] += rt;
205 dr[clust_id].aver2[ndr] += sqr(rt);
207 dr[clust_id].aver_3[ndr] += drt;
208 dr[clust_id].aver_6[ndr] += fcd->disres.Rt_6[0];
210 snew(fshift, SHIFTS);
211 ener = interaction_function[F_DISRES].ifunc(n, &forceatoms[i], forceparams,
212 (const rvec*)x, f, fshift,
213 pbc, g, lam, &dvdl, NULL, fcd, NULL);
215 viol = fcd->disres.sumviol;
222 add5(forceparams[type].disres.label, viol);
229 for (j = 0; (j < isize); j++)
231 if (index[j] == forceparams[type].disres.label)
233 vvindex[j] = pow(fcd->disres.Rt_6[0], -1.0/6.0);
240 dr[clust_id].nv = nviol;
241 dr[clust_id].maxv = mviol;
242 dr[clust_id].sumv = tviol;
243 dr[clust_id].averv = tviol/ndr;
244 dr[clust_id].nframes++;
248 fprintf(stderr, "\nThere are %d restraints and %d pairs\n",
249 ndr, disres->nr/nat);
261 real up1, r, rT3, rT6, viol, violT3, violT6;
264 static int drs_comp(const void *a, const void *b)
268 da = (t_dr_stats *)a;
269 db = (t_dr_stats *)b;
271 if (da->viol > db->viol)
275 else if (da->viol < db->viol)
285 static void dump_dump(FILE *log, int ndr, t_dr_stats drs[])
287 static const char *core[] = { "All restraints", "Core restraints" };
288 static const char *tp[] = { "linear", "third power", "sixth power" };
289 real viol_tot, viol_max, viol = 0;
294 for (bCore = FALSE; (bCore <= TRUE); bCore++)
296 for (kkk = 0; (kkk < 3); kkk++)
302 for (i = 0; (i < ndr); i++)
304 if (!bCore || (bCore && drs[i].bCore))
312 viol = drs[i].violT3;
315 viol = drs[i].violT6;
318 gmx_incons("Dumping violations");
320 viol_max = max(viol_max, viol);
329 if ((nrestr > 0) || (bCore && (nrestr < ndr)))
332 fprintf(log, "+++++++ %s ++++++++\n", core[bCore]);
333 fprintf(log, "+++++++ Using %s averaging: ++++++++\n", tp[kkk]);
334 fprintf(log, "Sum of violations: %8.3f nm\n", viol_tot);
337 fprintf(log, "Average violation: %8.3f nm\n", viol_tot/nrestr);
339 fprintf(log, "Largest violation: %8.3f nm\n", viol_max);
340 fprintf(log, "Number of violated restraints: %d/%d\n", nviol, nrestr);
346 static void dump_viol(FILE *log, int ndr, t_dr_stats *drs, gmx_bool bLinear)
350 fprintf(log, "Restr. Core Up1 <r> <rT3> <rT6> <viol><violT3><violT6>\n");
351 for (i = 0; (i < ndr); i++)
353 if (bLinear && (drs[i].viol == 0))
357 fprintf(log, "%6d%5s%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n",
358 drs[i].index, yesno_names[drs[i].bCore],
359 drs[i].up1, drs[i].r, drs[i].rT3, drs[i].rT6,
360 drs[i].viol, drs[i].violT3, drs[i].violT6);
364 static gmx_bool is_core(int i, int isize, atom_id index[])
367 gmx_bool bIC = FALSE;
369 for (kk = 0; !bIC && (kk < isize); kk++)
371 bIC = (index[kk] == i);
377 static void dump_stats(FILE *log, int nsteps, int ndr, t_ilist *disres,
378 t_iparams ip[], t_dr_result *dr,
379 int isize, atom_id index[], t_atoms *atoms)
385 fprintf(log, "++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
388 nra = interaction_function[F_DISRES].nratoms+1;
389 for (i = 0; (i < ndr); i++)
391 /* Search for the appropriate restraint */
392 for (; (j < disres->nr); j += nra)
394 if (ip[disres->iatoms[j]].disres.label == i)
400 drs[i].bCore = is_core(i, isize, index);
401 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
402 drs[i].r = dr->aver1[i]/nsteps;
403 drs[i].rT3 = pow(dr->aver_3[i]/nsteps, -1.0/3.0);
404 drs[i].rT6 = pow(dr->aver_6[i]/nsteps, -1.0/6.0);
405 drs[i].viol = max(0, drs[i].r-drs[i].up1);
406 drs[i].violT3 = max(0, drs[i].rT3-drs[i].up1);
407 drs[i].violT6 = max(0, drs[i].rT6-drs[i].up1);
410 int j1 = disres->iatoms[j+1];
411 int j2 = disres->iatoms[j+2];
412 atoms->pdbinfo[j1].bfac += drs[i].violT3*5;
413 atoms->pdbinfo[j2].bfac += drs[i].violT3*5;
416 dump_viol(log, ndr, drs, FALSE);
418 fprintf(log, "+++ Sorted by linear averaged violations: +++\n");
419 qsort(drs, ndr, sizeof(drs[0]), drs_comp);
420 dump_viol(log, ndr, drs, TRUE);
422 dump_dump(log, ndr, drs);
427 static void dump_clust_stats(FILE *fp, int ndr, t_ilist *disres,
428 t_iparams ip[], t_blocka *clust, t_dr_result dr[],
429 char *clust_name[], int isize, atom_id index[])
431 int i, j, k, nra, mmm = 0;
432 double sumV, maxV, sumVT3, sumVT6, maxVT3, maxVT6;
436 fprintf(fp, "++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
437 fprintf(fp, "Cluster NFrames SumV MaxV SumVT MaxVT SumVS MaxVS\n");
441 for (k = 0; (k < clust->nr); k++)
443 if (dr[k].nframes == 0)
447 if (dr[k].nframes != (clust->index[k+1]-clust->index[k]))
449 gmx_fatal(FARGS, "Inconsistency in cluster %s.\n"
450 "Found %d frames in trajectory rather than the expected %d\n",
451 clust_name[k], dr[k].nframes,
452 clust->index[k+1]-clust->index[k]);
456 gmx_fatal(FARGS, "Inconsistency with cluster %d. Invalid name", k);
459 nra = interaction_function[F_DISRES].nratoms+1;
460 sumV = sumVT3 = sumVT6 = maxV = maxVT3 = maxVT6 = 0;
461 for (i = 0; (i < ndr); i++)
463 /* Search for the appropriate restraint */
464 for (; (j < disres->nr); j += nra)
466 if (ip[disres->iatoms[j]].disres.label == i)
472 drs[i].bCore = is_core(i, isize, index);
473 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
474 drs[i].r = dr[k].aver1[i]/dr[k].nframes;
475 if ((dr[k].aver_3[i] <= 0) || (dr[k].aver_3[i] != dr[k].aver_3[i]))
477 gmx_fatal(FARGS, "dr[%d].aver_3[%d] = %f", k, i, dr[k].aver_3[i]);
479 drs[i].rT3 = pow(dr[k].aver_3[i]/dr[k].nframes, -1.0/3.0);
480 drs[i].rT6 = pow(dr[k].aver_6[i]/dr[k].nframes, -1.0/6.0);
481 drs[i].viol = max(0, drs[i].r-drs[i].up1);
482 drs[i].violT3 = max(0, drs[i].rT3-drs[i].up1);
483 drs[i].violT6 = max(0, drs[i].rT6-drs[i].up1);
485 sumVT3 += drs[i].violT3;
486 sumVT6 += drs[i].violT6;
487 maxV = max(maxV, drs[i].viol);
488 maxVT3 = max(maxVT3, drs[i].violT3);
489 maxVT6 = max(maxVT6, drs[i].violT6);
491 if (strcmp(clust_name[k], "1000") == 0)
495 fprintf(fp, "%-10s%6d%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n",
497 dr[k].nframes, sumV, maxV, sumVT3, maxVT3, sumVT6, maxVT6);
504 static void init_dr_res(t_dr_result *dr, int ndr)
506 snew(dr->aver1, ndr+1);
507 snew(dr->aver2, ndr+1);
508 snew(dr->aver_3, ndr+1);
509 snew(dr->aver_6, ndr+1);
517 static void dump_disre_matrix(const char *fn, t_dr_result *dr, int ndr,
518 int nsteps, t_idef *idef, gmx_mtop_t *mtop,
519 real max_dr, int nlevels, gmx_bool bThird)
523 int n_res, a_offset, mb, mol, a;
525 int iii, i, j, nra, nratoms, tp, ri, rj, index, nlabel, label;
526 atom_id ai, aj, *ptr;
527 real **matrix, *t_res, hi, *w_dr, rav, rviol;
528 t_rgb rlo = { 1, 1, 1 };
529 t_rgb rhi = { 0, 0, 0 };
535 snew(resnr, mtop->natoms);
538 for (mb = 0; mb < mtop->nmolblock; mb++)
540 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
541 for (mol = 0; mol < mtop->molblock[mb].nmol; mol++)
543 for (a = 0; a < atoms->nr; a++)
545 resnr[a_offset+a] = n_res + atoms->atom[a].resind;
547 n_res += atoms->nres;
548 a_offset += atoms->nr;
553 for (i = 0; (i < n_res); i++)
558 for (i = 0; (i < n_res); i++)
560 snew(matrix[i], n_res);
562 nratoms = interaction_function[F_DISRES].nratoms;
563 nra = (idef->il[F_DISRES].nr/(nratoms+1));
569 for (i = 0; (i < idef->il[F_DISRES].nr); i += nratoms+1)
571 tp = idef->il[F_DISRES].iatoms[i];
572 label = idef->iparams[tp].disres.label;
576 /* Set index pointer */
580 gmx_fatal(FARGS, "nlabel is %d, label = %d", nlabel, label);
584 gmx_fatal(FARGS, "ndr = %d, index = %d", ndr, index);
586 /* Update the weight */
587 w_dr[index] = 1.0/nlabel;
596 printf("nlabel = %d, index = %d, ndr = %d\n", nlabel, index, ndr);
598 for (i = 0; (i < ndr); i++)
600 for (j = ptr[i]; (j < ptr[i+1]); j += nratoms+1)
602 tp = idef->il[F_DISRES].iatoms[j];
603 ai = idef->il[F_DISRES].iatoms[j+1];
604 aj = idef->il[F_DISRES].iatoms[j+2];
610 rav = pow(dr->aver_3[i]/nsteps, -1.0/3.0);
614 rav = dr->aver1[i]/nsteps;
618 fprintf(debug, "DR %d, atoms %d, %d, distance %g\n", i, ai, aj, rav);
620 rviol = max(0, rav-idef->iparams[tp].disres.up1);
621 matrix[ri][rj] += w_dr[i]*rviol;
622 matrix[rj][ri] += w_dr[i]*rviol;
623 hi = max(hi, matrix[ri][rj]);
624 hi = max(hi, matrix[rj][ri]);
634 printf("Warning: the maxdr that you have specified (%g) is smaller than\nthe largest value in your simulation (%g)\n", max_dr, hi);
638 printf("Highest level in the matrix will be %g\n", hi);
639 fp = ffopen(fn, "w");
640 write_xpm(fp, 0, "Distance Violations", "<V> (nm)", "Residue", "Residue",
641 n_res, n_res, t_res, t_res, matrix, 0, hi, rlo, rhi, &nlevels);
645 int gmx_disre(int argc, char *argv[])
647 const char *desc[] = {
648 "[TT]g_disre[tt] computes violations of distance restraints.",
649 "If necessary, all protons can be added to a protein molecule ",
650 "using the [TT]g_protonate[tt] program.[PAR]",
651 "The program always",
652 "computes the instantaneous violations rather than time-averaged,",
653 "because this analysis is done from a trajectory file afterwards",
654 "it does not make sense to use time averaging. However,",
655 "the time averaged values per restraint are given in the log file.[PAR]",
656 "An index file may be used to select specific restraints for",
658 "When the optional [TT]-q[tt] flag is given a [TT].pdb[tt] file coloured by the",
659 "amount of average violations.[PAR]",
660 "When the [TT]-c[tt] option is given, an index file will be read",
661 "containing the frames in your trajectory corresponding to the clusters",
662 "(defined in another manner) that you want to analyze. For these clusters",
663 "the program will compute average violations using the third power",
664 "averaging algorithm and print them in the log file."
667 static int nlevels = 20;
668 static real max_dr = 0;
669 static gmx_bool bThird = TRUE;
671 { "-ntop", FALSE, etINT, {&ntop},
672 "Number of large violations that are stored in the log file every step" },
673 { "-maxdr", FALSE, etREAL, {&max_dr},
674 "Maximum distance violation in matrix output. If less than or equal to 0 the maximum will be determined by the data." },
675 { "-nlevels", FALSE, etINT, {&nlevels},
676 "Number of levels in the matrix output" },
677 { "-third", FALSE, etBOOL, {&bThird},
678 "Use inverse third power averaging or linear for matrix output" }
681 FILE *out = NULL, *aver = NULL, *numv = NULL, *maxxv = NULL, *xvg = NULL;
687 t_atoms *atoms = NULL;
691 int ntopatoms, natoms, i, j, kkk;
694 rvec *x, *f, *xav = NULL;
698 atom_id *index = NULL, *ind_fit = NULL;
700 t_cluster_ndx *clust = NULL;
701 t_dr_result dr, *dr_clust = NULL;
703 real *vvindex = NULL, *w_rls = NULL;
705 t_pbc pbc, *pbc_null;
709 gmx_rmpbc_t gpbc = NULL;
712 { efTPX, NULL, NULL, ffREAD },
713 { efTRX, "-f", NULL, ffREAD },
714 { efXVG, "-ds", "drsum", ffWRITE },
715 { efXVG, "-da", "draver", ffWRITE },
716 { efXVG, "-dn", "drnum", ffWRITE },
717 { efXVG, "-dm", "drmax", ffWRITE },
718 { efXVG, "-dr", "restr", ffWRITE },
719 { efLOG, "-l", "disres", ffWRITE },
720 { efNDX, NULL, "viol", ffOPTRD },
721 { efPDB, "-q", "viol", ffOPTWR },
722 { efNDX, "-c", "clust", ffOPTRD },
723 { efXPM, "-x", "matrix", ffOPTWR }
725 #define NFILE asize(fnm)
727 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
728 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
733 fplog = ftp2FILE(efLOG, NFILE, fnm, "w");
740 read_tpxheader(ftp2fn(efTPX, NFILE, fnm), &header, FALSE, NULL, NULL);
741 snew(xtop, header.natoms);
742 read_tpx(ftp2fn(efTPX, NFILE, fnm), &ir, box, &ntopatoms, xtop, NULL, NULL, &mtop);
743 bPDB = opt2bSet("-q", NFILE, fnm);
746 snew(xav, ntopatoms);
747 snew(ind_fit, ntopatoms);
748 snew(w_rls, ntopatoms);
749 for (kkk = 0; (kkk < ntopatoms); kkk++)
756 *atoms = gmx_mtop_global_atoms(&mtop);
758 if (atoms->pdbinfo == NULL)
760 snew(atoms->pdbinfo, atoms->nr);
764 top = gmx_mtop_generate_local_top(&mtop, &ir);
768 if (ir.ePBC != epbcNONE)
770 if (ir.bPeriodicMols)
776 g = mk_graph(fplog, &top->idef, 0, mtop.natoms, FALSE, FALSE);
780 if (ftp2bSet(efNDX, NFILE, fnm))
782 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
783 xvg = xvgropen(opt2fn("-dr", NFILE, fnm), "Individual Restraints", "Time (ps)",
785 snew(vvindex, isize);
787 for (i = 0; (i < isize); i++)
791 sprintf(leg[i], "index %d", index[i]);
793 xvgr_legend(xvg, isize, (const char**)leg, oenv);
801 init_disres(fplog, &mtop, &ir, NULL, FALSE, &fcd, NULL, FALSE);
803 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
806 init_dr_res(&dr, fcd.disres.nres);
807 if (opt2bSet("-c", NFILE, fnm))
809 clust = cluster_index(fplog, opt2fn("-c", NFILE, fnm));
810 snew(dr_clust, clust->clust->nr+1);
811 for (i = 0; (i <= clust->clust->nr); i++)
813 init_dr_res(&dr_clust[i], fcd.disres.nres);
818 out = xvgropen(opt2fn("-ds", NFILE, fnm),
819 "Sum of Violations", "Time (ps)", "nm", oenv);
820 aver = xvgropen(opt2fn("-da", NFILE, fnm),
821 "Average Violation", "Time (ps)", "nm", oenv);
822 numv = xvgropen(opt2fn("-dn", NFILE, fnm),
823 "# Violations", "Time (ps)", "#", oenv);
824 maxxv = xvgropen(opt2fn("-dm", NFILE, fnm),
825 "Largest Violation", "Time (ps)", "nm", oenv);
828 mdatoms = init_mdatoms(fplog, &mtop, ir.efep != efepNO);
829 atoms2md(&mtop, &ir, 0, NULL, 0, mtop.natoms, mdatoms);
830 update_mdatoms(mdatoms, ir.fepvals->init_lambda);
832 if (ir.ePBC != epbcNONE)
834 gpbc = gmx_rmpbc_init(&top->idef, ir.ePBC, natoms);
840 if (ir.ePBC != epbcNONE)
842 if (ir.bPeriodicMols)
844 set_pbc(&pbc, ir.ePBC, box);
848 gmx_rmpbc(gpbc, natoms, box, x);
854 if (j > clust->maxframe)
856 gmx_fatal(FARGS, "There are more frames in the trajectory than in the cluster index file. t = %8f\n", t);
858 my_clust = clust->inv_clust[j];
859 range_check(my_clust, 0, clust->clust->nr);
860 check_viol(fplog, &(top->idef.il[F_DISRES]),
862 x, f, pbc_null, g, dr_clust, my_clust, isize, index, vvindex, &fcd);
866 check_viol(fplog, &(top->idef.il[F_DISRES]),
868 x, f, pbc_null, g, &dr, 0, isize, index, vvindex, &fcd);
872 reset_x(atoms->nr, ind_fit, atoms->nr, NULL, x, w_rls);
873 do_fit(atoms->nr, w_rls, x, x);
876 /* Store the first frame of the trajectory as 'characteristic'
877 * for colouring with violations.
879 for (kkk = 0; (kkk < atoms->nr); kkk++)
881 copy_rvec(x[kkk], xav[kkk]);
889 fprintf(xvg, "%10g", t);
890 for (i = 0; (i < isize); i++)
892 fprintf(xvg, " %10g", vvindex[i]);
896 fprintf(out, "%10g %10g\n", t, dr.sumv);
897 fprintf(aver, "%10g %10g\n", t, dr.averv);
898 fprintf(maxxv, "%10g %10g\n", t, dr.maxv);
899 fprintf(numv, "%10g %10d\n", t, dr.nv);
903 while (read_next_x(oenv, status, &t, x, box));
905 if (ir.ePBC != epbcNONE)
907 gmx_rmpbc_done(gpbc);
912 dump_clust_stats(fplog, fcd.disres.nres, &(top->idef.il[F_DISRES]),
913 top->idef.iparams, clust->clust, dr_clust,
914 clust->grpname, isize, index);
918 dump_stats(fplog, j, fcd.disres.nres, &(top->idef.il[F_DISRES]),
919 top->idef.iparams, &dr, isize, index,
920 bPDB ? atoms : NULL);
923 write_sto_conf(opt2fn("-q", NFILE, fnm),
924 "Coloured by average violation in Angstrom",
925 atoms, xav, NULL, ir.ePBC, box);
927 dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres,
928 j, &top->idef, &mtop, max_dr, nlevels, bThird);
936 do_view(oenv, opt2fn("-dr", NFILE, fnm), "-nxy");
938 do_view(oenv, opt2fn("-dn", NFILE, fnm), "-nxy");
939 do_view(oenv, opt2fn("-da", NFILE, fnm), "-nxy");
940 do_view(oenv, opt2fn("-ds", NFILE, fnm), "-nxy");
941 do_view(oenv, opt2fn("-dm", NFILE, fnm), "-nxy");
944 gmx_log_close(fplog);