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.
53 #include "gromacs/fileio/confio.h"
54 #include "gromacs/fileio/trxio.h"
56 void print_one(const output_env_t oenv, const char *base, const char *name,
57 const char *title, const char *ylabel, int nf, real time[],
61 char buf[256], t2[256];
64 sprintf(buf, "%s%s.xvg", base, name);
65 fprintf(stderr, "\rPrinting %s ", buf);
66 sprintf(t2, "%s %s", title, name);
67 fp = xvgropen(buf, t2, "Time (ps)", ylabel, oenv);
68 for (k = 0; (k < nf); k++)
70 fprintf(fp, "%10g %10g\n", time[k], data[k]);
75 static int calc_RBbin(real phi, int gmx_unused multiplicity, real gmx_unused core_frac)
77 /* multiplicity and core_frac NOT used,
78 * just given to enable use of pt-to-fn in caller low_ana_dih_trans*/
79 static const real r30 = M_PI/6.0;
80 static const real r90 = M_PI/2.0;
81 static const real r150 = M_PI*5.0/6.0;
83 if ((phi < r30) && (phi > -r30))
87 else if ((phi > -r150) && (phi < -r90))
91 else if ((phi < r150) && (phi > r90))
98 static int calc_Nbin(real phi, int multiplicity, real core_frac)
100 static const real r360 = 360*DEG2RAD;
101 real rot_width, core_width, core_offset, low, hi;
103 /* with multiplicity 3 and core_frac 0.5
104 * 0<g(-)<120, 120<t<240, 240<g(+)<360
105 * 0< bin0 < 30, 30<bin1<90, 90<bin0<150, 150<bin2<210, 210<bin0<270, 270<bin3<330, 330<bin0<360
106 * so with multiplicity 3, bin1 is core g(-), bin2 is core t, bin3 is
107 core g(+), bin0 is between rotamers */
113 rot_width = 360/multiplicity;
114 core_width = core_frac * rot_width;
115 core_offset = (rot_width - core_width)/2.0;
116 for (bin = 1; bin <= multiplicity; bin++)
118 low = ((bin - 1) * rot_width ) + core_offset;
119 hi = ((bin - 1) * rot_width ) + core_offset + core_width;
122 if ((phi > low) && (phi < hi))
130 void ana_dih_trans(const char *fn_trans, const char *fn_histo,
131 real **dih, int nframes, int nangles,
132 const char *grpname, real *time, gmx_bool bRb,
133 const output_env_t oenv)
135 /* just a wrapper; declare extra args, then chuck away at end. */
143 snew(multiplicity, nangles);
144 for (k = 0; (k < nangles); k++)
149 low_ana_dih_trans(TRUE, fn_trans, TRUE, fn_histo, maxchi,
150 dih, nlist, dlist, nframes,
151 nangles, grpname, multiplicity, time, bRb, 0.5, oenv);
157 void low_ana_dih_trans(gmx_bool bTrans, const char *fn_trans,
158 gmx_bool bHisto, const char *fn_histo, int maxchi,
159 real **dih, int nlist, t_dlist dlist[], int nframes,
160 int nangles, const char *grpname, int multiplicity[],
161 real *time, gmx_bool bRb, real core_frac,
162 const output_env_t oenv)
167 int i, j, k, Dih, ntrans;
168 int cur_bin, new_bin;
171 int (*calc_bin)(real, int, real);
178 /* Assumes the frames are equally spaced in time */
179 dt = (time[nframes-1]-time[0])/(nframes-1);
181 /* Analysis of dihedral transitions */
182 fprintf(stderr, "Now calculating transitions...\n");
186 calc_bin = calc_RBbin;
190 calc_bin = calc_Nbin;
193 for (k = 0; k < NROT; k++)
195 snew(rot_occ[k], nangles);
196 for (i = 0; (i < nangles); i++)
204 /* dih[i][j] is the dihedral angle i in frame j */
206 for (i = 0; (i < nangles); i++)
211 mind = maxd = prev = dih[i][0];
213 cur_bin = calc_bin(dih[i][0], multiplicity[i], core_frac);
214 rot_occ[cur_bin][i]++;
216 for (j = 1; (j < nframes); j++)
218 new_bin = calc_bin(dih[i][j], multiplicity[i], core_frac);
219 rot_occ[new_bin][i]++;
225 else if ((new_bin != 0) && (cur_bin != new_bin))
233 /* why is all this md rubbish periodic? Remove 360 degree periodicity */
234 if ( (dih[i][j] - prev) > M_PI)
238 else if ( (dih[i][j] - prev) < -M_PI)
245 mind = min(mind, dih[i][j]);
246 maxd = max(maxd, dih[i][j]);
247 if ( (maxd - mind) > 2*M_PI/3) /* or 120 degrees, assuming */
248 { /* multiplicity 3. Not so general.*/
251 maxd = mind = dih[i][j]; /* get ready for next transition */
256 for (k = 0; k < NROT; k++)
258 rot_occ[k][i] /= nframes;
261 fprintf(stderr, "Total number of transitions: %10d\n", ntrans);
264 ttime = (dt*nframes*nangles)/ntrans;
265 fprintf(stderr, "Time between transitions: %10.3f ps\n", ttime);
268 /* new by grs - copy transitions from tr_h[] to dlist->ntr[]
269 * and rotamer populations from rot_occ to dlist->rot_occ[]
270 * based on fn histogramming in g_chi. diff roles for i and j here */
273 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
275 for (i = 0; (i < nlist); i++)
277 if (((Dih < edOmega) ) ||
278 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
279 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)))
281 /* grs debug printf("Not OK? i %d j %d Dih %d \n", i, j, Dih) ; */
282 dlist[i].ntr[Dih] = tr_h[j];
283 for (k = 0; k < NROT; k++)
285 dlist[i].rot_occ[Dih][k] = rot_occ[k][j];
292 /* end addition by grs */
296 sprintf(title, "Number of transitions: %s", grpname);
297 fp = xvgropen(fn_trans, title, "Time (ps)", "# transitions/timeframe", oenv);
298 for (j = 0; (j < nframes); j++)
300 fprintf(fp, "%10.3f %10d\n", time[j], tr_f[j]);
305 /* Compute histogram from # transitions per dihedral */
307 for (j = 0; (j < nframes); j++)
311 for (i = 0; (i < nangles); i++)
315 for (j = nframes; ((tr_f[j-1] == 0) && (j > 0)); j--)
323 sprintf(title, "Transition time: %s", grpname);
324 fp = xvgropen(fn_histo, title, "Time (ps)", "#", oenv);
325 for (i = j-1; (i > 0); i--)
329 fprintf(fp, "%10.3f %10d\n", ttime/i, tr_f[i]);
337 for (k = 0; k < NROT; k++)
344 void mk_multiplicity_lookup (int *multiplicity, int maxchi,
345 int nlist, t_dlist dlist[], int nangles)
347 /* new by grs - for dihedral j (as in dih[j]) get multiplicity from dlist
348 * and store in multiplicity[j]
355 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
357 for (i = 0; (i < nlist); i++)
359 strncpy(name, dlist[i].name, 3);
361 if (((Dih < edOmega) ) ||
362 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
363 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)))
365 /* default - we will correct the rest below */
368 /* make omegas 2fold, though doesn't make much more sense than 3 */
369 if (Dih == edOmega && (has_dihedral(edOmega, &(dlist[i]))))
374 /* dihedrals to aromatic rings, COO, CONH2 or guanidinium are 2fold*/
375 if (Dih > edOmega && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1))
377 if ( ((strstr(name, "PHE") != NULL) && (Dih == edChi2)) ||
378 ((strstr(name, "TYR") != NULL) && (Dih == edChi2)) ||
379 ((strstr(name, "PTR") != NULL) && (Dih == edChi2)) ||
380 ((strstr(name, "TRP") != NULL) && (Dih == edChi2)) ||
381 ((strstr(name, "HIS") != NULL) && (Dih == edChi2)) ||
382 ((strstr(name, "GLU") != NULL) && (Dih == edChi3)) ||
383 ((strstr(name, "ASP") != NULL) && (Dih == edChi2)) ||
384 ((strstr(name, "GLN") != NULL) && (Dih == edChi3)) ||
385 ((strstr(name, "ASN") != NULL) && (Dih == edChi2)) ||
386 ((strstr(name, "ARG") != NULL) && (Dih == edChi4)) )
397 fprintf(stderr, "WARNING: not all dihedrals found in topology (only %d out of %d)!\n",
400 /* Check for remaining dihedrals */
401 for (; (j < nangles); j++)
408 void mk_chi_lookup (int **lookup, int maxchi,
409 int nlist, t_dlist dlist[])
412 /* by grs. should rewrite everything to use this. (but haven't,
413 * and at mmt only used in get_chi_product_traj
414 * returns the dihed number given the residue number (from-0)
415 * and chi (from-0) nr. -1 for chi undefined for that res (eg gly, ala..)*/
420 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
422 for (i = 0; (i < nlist); i++)
425 if (((Dih < edOmega) ) ||
426 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
427 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)))
429 /* grs debug printf("Not OK? i %d j %d Dih %d \n", i, j, Dih) ; */
446 void get_chi_product_traj (real **dih, int nframes, int nlist,
447 int maxchi, t_dlist dlist[], real time[],
448 int **lookup, int *multiplicity, gmx_bool bRb, gmx_bool bNormalize,
449 real core_frac, gmx_bool bAll, const char *fnall,
450 const output_env_t oenv)
453 gmx_bool bRotZero, bHaveChi = FALSE;
454 int accum = 0, index, i, j, k, Xi, n, b;
459 char hisfile[256], histitle[256], *namept;
461 int (*calc_bin)(real, int, real);
463 /* Analysis of dihedral transitions */
464 fprintf(stderr, "Now calculating Chi product trajectories...\n");
468 calc_bin = calc_RBbin;
472 calc_bin = calc_Nbin;
475 snew(chi_prtrj, nframes);
477 /* file for info on all residues */
480 fpall = xvgropen(fnall, "Cumulative Rotamers", "Residue", "Probability", oenv);
484 fpall = xvgropen(fnall, "Cumulative Rotamers", "Residue", "# Counts", oenv);
487 for (i = 0; (i < nlist); i++)
490 /* get nbin, the nr. of cumulative rotamers that need to be considered */
492 for (Xi = 0; Xi < maxchi; Xi++)
494 index = lookup[i][Xi]; /* chi_(Xi+1) of res i (-1 if off end) */
497 n = multiplicity[index];
501 nbin += 1; /* for the "zero rotamer", outside the core region */
503 for (j = 0; (j < nframes); j++)
508 index = lookup[i][0]; /* index into dih of chi1 of res i */
517 b = calc_bin(dih[index][j], multiplicity[index], core_frac);
523 for (Xi = 1; Xi < maxchi; Xi++)
525 index = lookup[i][Xi]; /* chi_(Xi+1) of res i (-1 if off end) */
528 n = multiplicity[index];
529 b = calc_bin(dih[index][j], n, core_frac);
530 accum = n * accum + b - 1;
545 chi_prtrj[j] = accum;
557 /* print cuml rotamer vs time */
558 print_one(oenv, "chiproduct", dlist[i].name, "chi product for",
559 "cumulative rotamer", nframes, time, chi_prtrj);
562 /* make a histogram pf culm. rotamer occupancy too */
563 snew(chi_prhist, nbin);
564 make_histo(NULL, nframes, chi_prtrj, nbin, chi_prhist, 0, nbin);
567 sprintf(hisfile, "histo-chiprod%s.xvg", dlist[i].name);
568 sprintf(histitle, "cumulative rotamer distribution for %s", dlist[i].name);
569 fprintf(stderr, " and %s ", hisfile);
570 fp = xvgropen(hisfile, histitle, "number", "", oenv);
571 fprintf(fp, "@ xaxis tick on\n");
572 fprintf(fp, "@ xaxis tick major 1\n");
573 fprintf(fp, "@ type xy\n");
574 for (k = 0; (k < nbin); k++)
578 fprintf(fp, "%5d %10g\n", k, (1.0*chi_prhist[k])/nframes);
582 fprintf(fp, "%5d %10d\n", k, chi_prhist[k]);
589 /* and finally print out occupancies to a single file */
590 /* get the gmx from-1 res nr by setting a ptr to the number part
591 * of dlist[i].name - potential bug for 4-letter res names... */
592 namept = dlist[i].name + 3;
593 fprintf(fpall, "%5s ", namept);
594 for (k = 0; (k < nbin); k++)
598 fprintf(fpall, " %10g", (1.0*chi_prhist[k])/nframes);
602 fprintf(fpall, " %10d", chi_prhist[k]);
605 fprintf(fpall, "\n");
614 fprintf(stderr, "\n");
618 void calc_distribution_props(int nh, int histo[], real start,
619 int nkkk, t_karplus kkk[],
622 real d, dc, ds, c1, c2, tdc, tds;
623 real fac, ang, invth, Jc;
628 gmx_fatal(FARGS, "No points in histogram (%s, %d)", __FILE__, __LINE__);
632 /* Compute normalisation factor */
634 for (j = 0; (j < nh); j++)
640 for (i = 0; (i < nkkk); i++)
646 for (j = 0; (j < nh); j++)
656 for (i = 0; (i < nkkk); i++)
658 c1 = cos(ang+kkk[i].offset);
660 Jc = (kkk[i].A*c2 + kkk[i].B*c1 + kkk[i].C);
661 kkk[i].Jc += histo[j]*Jc;
662 kkk[i].Jcsig += histo[j]*sqr(Jc);
665 for (i = 0; (i < nkkk); i++)
668 kkk[i].Jcsig = sqrt(kkk[i].Jcsig/th-sqr(kkk[i].Jc));
670 *S2 = tdc*tdc+tds*tds;
673 static void calc_angles(t_pbc *pbc,
674 int n3, atom_id index[], real ang[], rvec x_s[])
680 for (i = ix = 0; (ix < n3); i++, ix += 3)
682 ang[i] = bond_angle(x_s[index[ix]], x_s[index[ix+1]], x_s[index[ix+2]],
683 pbc, r_ij, r_kj, &costh, &t1, &t2);
687 fprintf(debug, "Angle[0]=%g, costh=%g, index0 = %d, %d, %d\n",
688 ang[0], costh, index[0], index[1], index[2]);
689 pr_rvec(debug, 0, "rij", r_ij, DIM, TRUE);
690 pr_rvec(debug, 0, "rkj", r_kj, DIM, TRUE);
694 static real calc_fraction(real angles[], int nangles)
697 real trans = 0, gauche = 0;
700 for (i = 0; i < nangles; i++)
702 angle = angles[i] * RAD2DEG;
704 if (angle > 135 && angle < 225)
708 else if (angle > 270 && angle < 330)
712 else if (angle < 90 && angle > 30)
717 if (trans+gauche > 0)
719 return trans/(trans+gauche);
727 static void calc_dihs(t_pbc *pbc,
728 int n4, atom_id index[], real ang[], rvec x_s[])
730 int i, ix, t1, t2, t3;
731 rvec r_ij, r_kj, r_kl, m, n;
734 for (i = ix = 0; (ix < n4); i++, ix += 4)
736 aaa = dih_angle(x_s[index[ix]], x_s[index[ix+1]], x_s[index[ix+2]],
737 x_s[index[ix+3]], pbc,
738 r_ij, r_kj, r_kl, m, n,
739 &sign, &t1, &t2, &t3);
741 ang[i] = aaa; /* not taking into account ryckaert bellemans yet */
745 void make_histo(FILE *log,
746 int ndata, real data[], int npoints, int histo[],
747 real minx, real maxx)
754 minx = maxx = data[0];
755 for (i = 1; (i < ndata); i++)
757 minx = min(minx, data[i]);
758 maxx = max(maxx, data[i]);
760 fprintf(log, "Min data: %10g Max data: %10g\n", minx, maxx);
762 dx = (double)npoints/(maxx-minx);
766 "Histogramming: ndata=%d, nhisto=%d, minx=%g,maxx=%g,dx=%g\n",
767 ndata, npoints, minx, maxx, dx);
769 for (i = 0; (i < ndata); i++)
771 ind = (data[i]-minx)*dx;
772 if ((ind >= 0) && (ind < npoints))
778 fprintf(log, "index = %d, data[%d] = %g\n", ind, i, data[i]);
783 void normalize_histo(int npoints, int histo[], real dx, real normhisto[])
789 for (i = 0; (i < npoints); i++)
795 fprintf(stderr, "Empty histogram!\n");
799 for (i = 0; (i < npoints); i++)
801 normhisto[i] = fac*histo[i];
805 void read_ang_dih(const char *trj_fn,
806 gmx_bool bAngles, gmx_bool bSaveAll, gmx_bool bRb, gmx_bool bPBC,
807 int maxangstat, int angstat[],
808 int *nframes, real **time,
809 int isize, atom_id index[],
813 const output_env_t oenv)
817 int i, angind, natoms, total, teller;
818 int nangles, n_alloc;
819 real t, fraction, pifac, aa, angle;
827 natoms = read_first_x(oenv, &status, trj_fn, &t, &x, box);
839 snew(angles[cur], nangles);
840 snew(angles[prev], nangles);
842 /* Start the loop over frames */
852 if (teller >= n_alloc)
857 for (i = 0; (i < nangles); i++)
859 srenew(dih[i], n_alloc);
862 srenew(*time, n_alloc);
863 srenew(*trans_frac, n_alloc);
864 srenew(*aver_angle, n_alloc);
871 set_pbc(pbc, -1, box);
876 calc_angles(pbc, isize, index, angles[cur], x);
880 calc_dihs(pbc, isize, index, angles[cur], x);
883 fraction = calc_fraction(angles[cur], nangles);
884 (*trans_frac)[teller] = fraction;
886 /* Change Ryckaert-Bellemans dihedrals to polymer convention
887 * Modified 990913 by Erik:
888 * We actually shouldn't change the convention, since it's
889 * calculated from polymer above, but we change the intervall
890 * from [-180,180] to [0,360].
894 for (i = 0; (i < nangles); i++)
896 if (angles[cur][i] <= 0.0)
898 angles[cur][i] += 2*M_PI;
903 /* Periodicity in dihedral space... */
906 for (i = 0; (i < nangles); i++)
908 real dd = angles[cur][i];
909 angles[cur][i] = atan2(sin(dd), cos(dd));
916 for (i = 0; (i < nangles); i++)
918 while (angles[cur][i] <= angles[prev][i] - M_PI)
920 angles[cur][i] += 2*M_PI;
922 while (angles[cur][i] > angles[prev][i] + M_PI)
924 angles[cur][i] -= 2*M_PI;
933 for (i = 0; (i < nangles); i++)
935 aa = aa+angles[cur][i];
937 /* angle in rad / 2Pi * max determines bin. bins go from 0 to maxangstat,
938 even though scale goes from -pi to pi (dihedral) or -pi/2 to pi/2
939 (angle) Basically: translate the x-axis by Pi. Translate it back by
943 angle = angles[cur][i];
946 while (angle < -M_PI)
950 while (angle >= M_PI)
958 /* Update the distribution histogram */
959 angind = (int) ((angle*maxangstat)/pifac + 0.5);
960 if (angind == maxangstat)
964 if ( (angind < 0) || (angind >= maxangstat) )
966 /* this will never happen */
967 gmx_fatal(FARGS, "angle (%f) index out of range (0..%d) : %d\n",
968 angle, maxangstat, angind);
972 if (angind == maxangstat)
974 fprintf(stderr, "angle %d fr %d = %g\n", i, cur, angle);
980 /* average over all angles */
981 (*aver_angle)[teller] = (aa/nangles);
983 /* this copies all current dih. angles to dih[i], teller is frame */
986 for (i = 0; i < nangles; i++)
988 dih[i][teller] = angles[cur][i];
995 /* Increment loop counter */
998 while (read_next_x(oenv, status, &t, x, box));
1003 sfree(angles[prev]);