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
47 #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 "gromacs/fileio/matio.h"
63 #include "mtop_util.h"
77 real sumv, averv, maxv;
78 real *aver1, *aver2, *aver_3, *aver_6;
81 static void init5(int n)
87 static void reset5(void)
91 for (i = 0; (i < ntop); i++)
98 int tpcomp(const void *a, const void *b)
106 return (1e7*(tpb->v-tpa->v));
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 qsort(top, ntop, sizeof(top[0]), tpcomp);
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, atom_id index[], real vvindex[],
154 int i, j, nat, n, type, nviol, ndr, label;
155 real ener, 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(n, &forceatoms[i], forceparams,
197 (const rvec*)x, pbc, fcd, NULL);
199 if (fcd->disres.Rt_6[0] <= 0)
201 gmx_fatal(FARGS, "ndr = %d, rt_6 = %f", ndr, fcd->disres.Rt_6[0]);
204 rt = pow(fcd->disres.Rt_6[0], -1.0/6.0);
205 dr[clust_id].aver1[ndr] += rt;
206 dr[clust_id].aver2[ndr] += sqr(rt);
208 dr[clust_id].aver_3[ndr] += drt;
209 dr[clust_id].aver_6[ndr] += fcd->disres.Rt_6[0];
211 snew(fshift, SHIFTS);
212 ener = interaction_function[F_DISRES].ifunc(n, &forceatoms[i], forceparams,
213 (const rvec*)x, f, fshift,
214 pbc, g, lam, &dvdl, NULL, fcd, NULL);
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] = pow(fcd->disres.Rt_6[0], -1.0/6.0);
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 int drs_comp(const void *a, const void *b)
269 da = (t_dr_stats *)a;
270 db = (t_dr_stats *)b;
272 if (da->viol > db->viol)
276 else if (da->viol < db->viol)
286 static void dump_dump(FILE *log, int ndr, t_dr_stats drs[])
288 static const char *core[] = { "All restraints", "Core restraints" };
289 static const char *tp[] = { "linear", "third power", "sixth power" };
290 real viol_tot, viol_max, viol = 0;
295 for (bCore = FALSE; (bCore <= TRUE); bCore++)
297 for (kkk = 0; (kkk < 3); kkk++)
303 for (i = 0; (i < ndr); i++)
305 if (!bCore || (bCore && drs[i].bCore))
313 viol = drs[i].violT3;
316 viol = drs[i].violT6;
319 gmx_incons("Dumping violations");
321 viol_max = max(viol_max, viol);
330 if ((nrestr > 0) || (bCore && (nrestr < ndr)))
333 fprintf(log, "+++++++ %s ++++++++\n", core[bCore]);
334 fprintf(log, "+++++++ Using %s averaging: ++++++++\n", tp[kkk]);
335 fprintf(log, "Sum of violations: %8.3f nm\n", viol_tot);
338 fprintf(log, "Average violation: %8.3f nm\n", viol_tot/nrestr);
340 fprintf(log, "Largest violation: %8.3f nm\n", viol_max);
341 fprintf(log, "Number of violated restraints: %d/%d\n", nviol, nrestr);
347 static void dump_viol(FILE *log, int ndr, t_dr_stats *drs, gmx_bool bLinear)
351 fprintf(log, "Restr. Core Up1 <r> <rT3> <rT6> <viol><violT3><violT6>\n");
352 for (i = 0; (i < ndr); i++)
354 if (bLinear && (drs[i].viol == 0))
358 fprintf(log, "%6d%5s%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n",
359 drs[i].index, yesno_names[drs[i].bCore],
360 drs[i].up1, drs[i].r, drs[i].rT3, drs[i].rT6,
361 drs[i].viol, drs[i].violT3, drs[i].violT6);
365 static gmx_bool is_core(int i, int isize, atom_id index[])
368 gmx_bool bIC = FALSE;
370 for (kk = 0; !bIC && (kk < isize); kk++)
372 bIC = (index[kk] == i);
378 static void dump_stats(FILE *log, int nsteps, int ndr, t_ilist *disres,
379 t_iparams ip[], t_dr_result *dr,
380 int isize, atom_id index[], t_atoms *atoms)
386 fprintf(log, "++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
389 nra = interaction_function[F_DISRES].nratoms+1;
390 for (i = 0; (i < ndr); i++)
392 /* Search for the appropriate restraint */
393 for (; (j < disres->nr); j += nra)
395 if (ip[disres->iatoms[j]].disres.label == i)
401 drs[i].bCore = is_core(i, isize, index);
402 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
403 drs[i].r = dr->aver1[i]/nsteps;
404 drs[i].rT3 = pow(dr->aver_3[i]/nsteps, -1.0/3.0);
405 drs[i].rT6 = pow(dr->aver_6[i]/nsteps, -1.0/6.0);
406 drs[i].viol = max(0, drs[i].r-drs[i].up1);
407 drs[i].violT3 = max(0, drs[i].rT3-drs[i].up1);
408 drs[i].violT6 = max(0, drs[i].rT6-drs[i].up1);
411 int j1 = disres->iatoms[j+1];
412 int j2 = disres->iatoms[j+2];
413 atoms->pdbinfo[j1].bfac += drs[i].violT3*5;
414 atoms->pdbinfo[j2].bfac += drs[i].violT3*5;
417 dump_viol(log, ndr, drs, FALSE);
419 fprintf(log, "+++ Sorted by linear averaged violations: +++\n");
420 qsort(drs, ndr, sizeof(drs[0]), drs_comp);
421 dump_viol(log, ndr, drs, TRUE);
423 dump_dump(log, ndr, drs);
428 static void dump_clust_stats(FILE *fp, int ndr, t_ilist *disres,
429 t_iparams ip[], t_blocka *clust, t_dr_result dr[],
430 char *clust_name[], int isize, atom_id index[])
432 int i, j, k, nra, mmm = 0;
433 double sumV, maxV, sumVT3, sumVT6, maxVT3, maxVT6;
437 fprintf(fp, "++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
438 fprintf(fp, "Cluster NFrames SumV MaxV SumVT MaxVT SumVS MaxVS\n");
442 for (k = 0; (k < clust->nr); k++)
444 if (dr[k].nframes == 0)
448 if (dr[k].nframes != (clust->index[k+1]-clust->index[k]))
450 gmx_fatal(FARGS, "Inconsistency in cluster %s.\n"
451 "Found %d frames in trajectory rather than the expected %d\n",
452 clust_name[k], dr[k].nframes,
453 clust->index[k+1]-clust->index[k]);
457 gmx_fatal(FARGS, "Inconsistency with cluster %d. Invalid name", k);
460 nra = interaction_function[F_DISRES].nratoms+1;
461 sumV = sumVT3 = sumVT6 = maxV = maxVT3 = maxVT6 = 0;
462 for (i = 0; (i < ndr); i++)
464 /* Search for the appropriate restraint */
465 for (; (j < disres->nr); j += nra)
467 if (ip[disres->iatoms[j]].disres.label == i)
473 drs[i].bCore = is_core(i, isize, index);
474 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
475 drs[i].r = dr[k].aver1[i]/dr[k].nframes;
476 if ((dr[k].aver_3[i] <= 0) || (dr[k].aver_3[i] != dr[k].aver_3[i]))
478 gmx_fatal(FARGS, "dr[%d].aver_3[%d] = %f", k, i, dr[k].aver_3[i]);
480 drs[i].rT3 = pow(dr[k].aver_3[i]/dr[k].nframes, -1.0/3.0);
481 drs[i].rT6 = pow(dr[k].aver_6[i]/dr[k].nframes, -1.0/6.0);
482 drs[i].viol = max(0, drs[i].r-drs[i].up1);
483 drs[i].violT3 = max(0, drs[i].rT3-drs[i].up1);
484 drs[i].violT6 = max(0, drs[i].rT6-drs[i].up1);
486 sumVT3 += drs[i].violT3;
487 sumVT6 += drs[i].violT6;
488 maxV = max(maxV, drs[i].viol);
489 maxVT3 = max(maxVT3, drs[i].violT3);
490 maxVT6 = max(maxVT6, drs[i].violT6);
492 if (strcmp(clust_name[k], "1000") == 0)
496 fprintf(fp, "%-10s%6d%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n",
498 dr[k].nframes, sumV, maxV, sumVT3, maxVT3, sumVT6, maxVT6);
505 static void init_dr_res(t_dr_result *dr, int ndr)
507 snew(dr->aver1, ndr+1);
508 snew(dr->aver2, ndr+1);
509 snew(dr->aver_3, ndr+1);
510 snew(dr->aver_6, ndr+1);
518 static void dump_disre_matrix(const char *fn, t_dr_result *dr, int ndr,
519 int nsteps, t_idef *idef, gmx_mtop_t *mtop,
520 real max_dr, int nlevels, gmx_bool bThird)
524 int n_res, a_offset, mb, mol, a;
526 int iii, i, j, nra, nratoms, tp, ri, rj, index, nlabel, label;
527 atom_id ai, aj, *ptr;
528 real **matrix, *t_res, hi, *w_dr, rav, rviol;
529 t_rgb rlo = { 1, 1, 1 };
530 t_rgb rhi = { 0, 0, 0 };
536 snew(resnr, mtop->natoms);
539 for (mb = 0; mb < mtop->nmolblock; mb++)
541 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
542 for (mol = 0; mol < mtop->molblock[mb].nmol; mol++)
544 for (a = 0; a < atoms->nr; a++)
546 resnr[a_offset+a] = n_res + atoms->atom[a].resind;
548 n_res += atoms->nres;
549 a_offset += atoms->nr;
554 for (i = 0; (i < n_res); i++)
559 for (i = 0; (i < n_res); i++)
561 snew(matrix[i], n_res);
563 nratoms = interaction_function[F_DISRES].nratoms;
564 nra = (idef->il[F_DISRES].nr/(nratoms+1));
570 for (i = 0; (i < idef->il[F_DISRES].nr); i += nratoms+1)
572 tp = idef->il[F_DISRES].iatoms[i];
573 label = idef->iparams[tp].disres.label;
577 /* Set index pointer */
581 gmx_fatal(FARGS, "nlabel is %d, label = %d", nlabel, label);
585 gmx_fatal(FARGS, "ndr = %d, index = %d", ndr, index);
587 /* Update the weight */
588 w_dr[index] = 1.0/nlabel;
597 printf("nlabel = %d, index = %d, ndr = %d\n", nlabel, index, ndr);
599 for (i = 0; (i < ndr); i++)
601 for (j = ptr[i]; (j < ptr[i+1]); j += nratoms+1)
603 tp = idef->il[F_DISRES].iatoms[j];
604 ai = idef->il[F_DISRES].iatoms[j+1];
605 aj = idef->il[F_DISRES].iatoms[j+2];
611 rav = pow(dr->aver_3[i]/nsteps, -1.0/3.0);
615 rav = dr->aver1[i]/nsteps;
619 fprintf(debug, "DR %d, atoms %d, %d, distance %g\n", i, ai, aj, rav);
621 rviol = max(0, rav-idef->iparams[tp].disres.up1);
622 matrix[ri][rj] += w_dr[i]*rviol;
623 matrix[rj][ri] += w_dr[i]*rviol;
624 hi = max(hi, matrix[ri][rj]);
625 hi = max(hi, matrix[rj][ri]);
635 printf("Warning: the maxdr that you have specified (%g) is smaller than\nthe largest value in your simulation (%g)\n", max_dr, hi);
639 printf("Highest level in the matrix will be %g\n", hi);
640 fp = ffopen(fn, "w");
641 write_xpm(fp, 0, "Distance Violations", "<V> (nm)", "Residue", "Residue",
642 n_res, n_res, t_res, t_res, matrix, 0, hi, rlo, rhi, &nlevels);
646 int gmx_disre(int argc, char *argv[])
648 const char *desc[] = {
649 "[TT]g_disre[tt] computes violations of distance restraints.",
650 "If necessary, all protons can be added to a protein molecule ",
651 "using the [TT]g_protonate[tt] program.[PAR]",
652 "The program always",
653 "computes the instantaneous violations rather than time-averaged,",
654 "because this analysis is done from a trajectory file afterwards",
655 "it does not make sense to use time averaging. However,",
656 "the time averaged values per restraint are given in the log file.[PAR]",
657 "An index file may be used to select specific restraints for",
659 "When the optional [TT]-q[tt] flag is given a [TT].pdb[tt] file coloured by the",
660 "amount of average violations.[PAR]",
661 "When the [TT]-c[tt] option is given, an index file will be read",
662 "containing the frames in your trajectory corresponding to the clusters",
663 "(defined in another manner) that you want to analyze. For these clusters",
664 "the program will compute average violations using the third power",
665 "averaging algorithm and print them in the log file."
668 static int nlevels = 20;
669 static real max_dr = 0;
670 static gmx_bool bThird = TRUE;
672 { "-ntop", FALSE, etINT, {&ntop},
673 "Number of large violations that are stored in the log file every step" },
674 { "-maxdr", FALSE, etREAL, {&max_dr},
675 "Maximum distance violation in matrix output. If less than or equal to 0 the maximum will be determined by the data." },
676 { "-nlevels", FALSE, etINT, {&nlevels},
677 "Number of levels in the matrix output" },
678 { "-third", FALSE, etBOOL, {&bThird},
679 "Use inverse third power averaging or linear for matrix output" }
682 FILE *out = NULL, *aver = NULL, *numv = NULL, *maxxv = NULL, *xvg = NULL;
688 t_atoms *atoms = NULL;
692 int ntopatoms, natoms, i, j, kkk;
695 rvec *x, *f, *xav = NULL;
699 atom_id *index = NULL, *ind_fit = NULL;
701 t_cluster_ndx *clust = NULL;
702 t_dr_result dr, *dr_clust = NULL;
704 real *vvindex = NULL, *w_rls = NULL;
706 t_pbc pbc, *pbc_null;
710 gmx_rmpbc_t gpbc = NULL;
713 { efTPX, NULL, NULL, ffREAD },
714 { efTRX, "-f", NULL, ffREAD },
715 { efXVG, "-ds", "drsum", ffWRITE },
716 { efXVG, "-da", "draver", ffWRITE },
717 { efXVG, "-dn", "drnum", ffWRITE },
718 { efXVG, "-dm", "drmax", ffWRITE },
719 { efXVG, "-dr", "restr", ffWRITE },
720 { efLOG, "-l", "disres", ffWRITE },
721 { efNDX, NULL, "viol", ffOPTRD },
722 { efPDB, "-q", "viol", ffOPTWR },
723 { efNDX, "-c", "clust", ffOPTRD },
724 { efXPM, "-x", "matrix", ffOPTWR }
726 #define NFILE asize(fnm)
728 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
729 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
734 fplog = ftp2FILE(efLOG, NFILE, fnm, "w");
741 read_tpxheader(ftp2fn(efTPX, NFILE, fnm), &header, FALSE, NULL, NULL);
742 snew(xtop, header.natoms);
743 read_tpx(ftp2fn(efTPX, NFILE, fnm), &ir, box, &ntopatoms, xtop, NULL, NULL, &mtop);
744 bPDB = opt2bSet("-q", NFILE, fnm);
747 snew(xav, ntopatoms);
748 snew(ind_fit, ntopatoms);
749 snew(w_rls, ntopatoms);
750 for (kkk = 0; (kkk < ntopatoms); kkk++)
757 *atoms = gmx_mtop_global_atoms(&mtop);
759 if (atoms->pdbinfo == NULL)
761 snew(atoms->pdbinfo, atoms->nr);
765 top = gmx_mtop_generate_local_top(&mtop, &ir);
769 if (ir.ePBC != epbcNONE)
771 if (ir.bPeriodicMols)
777 g = mk_graph(fplog, &top->idef, 0, mtop.natoms, FALSE, FALSE);
781 if (ftp2bSet(efNDX, NFILE, fnm))
783 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
784 xvg = xvgropen(opt2fn("-dr", NFILE, fnm), "Individual Restraints", "Time (ps)",
786 snew(vvindex, isize);
788 for (i = 0; (i < isize); i++)
792 sprintf(leg[i], "index %d", index[i]);
794 xvgr_legend(xvg, isize, (const char**)leg, oenv);
802 init_disres(fplog, &mtop, &ir, NULL, FALSE, &fcd, NULL, FALSE);
804 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
807 init_dr_res(&dr, fcd.disres.nres);
808 if (opt2bSet("-c", NFILE, fnm))
810 clust = cluster_index(fplog, opt2fn("-c", NFILE, fnm));
811 snew(dr_clust, clust->clust->nr+1);
812 for (i = 0; (i <= clust->clust->nr); i++)
814 init_dr_res(&dr_clust[i], fcd.disres.nres);
819 out = xvgropen(opt2fn("-ds", NFILE, fnm),
820 "Sum of Violations", "Time (ps)", "nm", oenv);
821 aver = xvgropen(opt2fn("-da", NFILE, fnm),
822 "Average Violation", "Time (ps)", "nm", oenv);
823 numv = xvgropen(opt2fn("-dn", NFILE, fnm),
824 "# Violations", "Time (ps)", "#", oenv);
825 maxxv = xvgropen(opt2fn("-dm", NFILE, fnm),
826 "Largest Violation", "Time (ps)", "nm", oenv);
829 mdatoms = init_mdatoms(fplog, &mtop, ir.efep != efepNO);
830 atoms2md(&mtop, &ir, 0, NULL, 0, mtop.natoms, mdatoms);
831 update_mdatoms(mdatoms, ir.fepvals->init_lambda);
833 if (ir.ePBC != epbcNONE)
835 gpbc = gmx_rmpbc_init(&top->idef, ir.ePBC, natoms);
841 if (ir.ePBC != epbcNONE)
843 if (ir.bPeriodicMols)
845 set_pbc(&pbc, ir.ePBC, box);
849 gmx_rmpbc(gpbc, natoms, box, x);
855 if (j > clust->maxframe)
857 gmx_fatal(FARGS, "There are more frames in the trajectory than in the cluster index file. t = %8f\n", t);
859 my_clust = clust->inv_clust[j];
860 range_check(my_clust, 0, clust->clust->nr);
861 check_viol(fplog, &(top->idef.il[F_DISRES]),
863 x, f, pbc_null, g, dr_clust, my_clust, isize, index, vvindex, &fcd);
867 check_viol(fplog, &(top->idef.il[F_DISRES]),
869 x, f, pbc_null, g, &dr, 0, isize, index, vvindex, &fcd);
873 reset_x(atoms->nr, ind_fit, atoms->nr, NULL, x, w_rls);
874 do_fit(atoms->nr, w_rls, x, x);
877 /* Store the first frame of the trajectory as 'characteristic'
878 * for colouring with violations.
880 for (kkk = 0; (kkk < atoms->nr); kkk++)
882 copy_rvec(x[kkk], xav[kkk]);
890 fprintf(xvg, "%10g", t);
891 for (i = 0; (i < isize); i++)
893 fprintf(xvg, " %10g", vvindex[i]);
897 fprintf(out, "%10g %10g\n", t, dr.sumv);
898 fprintf(aver, "%10g %10g\n", t, dr.averv);
899 fprintf(maxxv, "%10g %10g\n", t, dr.maxv);
900 fprintf(numv, "%10g %10d\n", t, dr.nv);
904 while (read_next_x(oenv, status, &t, x, box));
906 if (ir.ePBC != epbcNONE)
908 gmx_rmpbc_done(gpbc);
913 dump_clust_stats(fplog, fcd.disres.nres, &(top->idef.il[F_DISRES]),
914 top->idef.iparams, clust->clust, dr_clust,
915 clust->grpname, isize, index);
919 dump_stats(fplog, j, fcd.disres.nres, &(top->idef.il[F_DISRES]),
920 top->idef.iparams, &dr, isize, index,
921 bPDB ? atoms : NULL);
924 write_sto_conf(opt2fn("-q", NFILE, fnm),
925 "Coloured by average violation in Angstrom",
926 atoms, xav, NULL, ir.ePBC, box);
928 dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres,
929 j, &top->idef, &mtop, max_dr, nlevels, bThird);
937 do_view(oenv, opt2fn("-dr", NFILE, fnm), "-nxy");
939 do_view(oenv, opt2fn("-dn", NFILE, fnm), "-nxy");
940 do_view(oenv, opt2fn("-da", NFILE, fnm), "-nxy");
941 do_view(oenv, opt2fn("-ds", NFILE, fnm), "-nxy");
942 do_view(oenv, opt2fn("-dm", NFILE, fnm), "-nxy");
945 gmx_log_close(fplog);