Moved timing source to C++
[alexxy/gromacs.git] / src / gromacs / gmxana / anadih.cpp
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5  * Copyright (c) 2001-2004, The GROMACS development team.
6  * Copyright (c) 2013,2014,2015, 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.
10  *
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.
15  *
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.
20  *
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.
25  *
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.
33  *
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.
36  */
37 #include "gmxpre.h"
38
39 #include <cmath>
40 #include <cstdio>
41 #include <cstring>
42
43 #include <algorithm>
44
45 #include "gromacs/fileio/confio.h"
46 #include "gromacs/fileio/trxio.h"
47 #include "gromacs/fileio/xvgr.h"
48 #include "gromacs/gmxana/gstat.h"
49 #include "gromacs/legacyheaders/macros.h"
50 #include "gromacs/legacyheaders/txtdump.h"
51 #include "gromacs/legacyheaders/typedefs.h"
52 #include "gromacs/listed-forces/bonded.h"
53 #include "gromacs/math/units.h"
54 #include "gromacs/math/vec.h"
55 #include "gromacs/pbcutil/pbc.h"
56 #include "gromacs/utility/fatalerror.h"
57 #include "gromacs/utility/smalloc.h"
58
59 void print_one(const output_env_t oenv, const char *base, const char *name,
60                const char *title, const char *ylabel, int nf, real time[],
61                real data[])
62 {
63     FILE *fp;
64     char  buf[256], t2[256];
65     int   k;
66
67     sprintf(buf, "%s%s.xvg", base, name);
68     fprintf(stderr, "\rPrinting %s  ", buf);
69     sprintf(t2, "%s %s", title, name);
70     fp = xvgropen(buf, t2, "Time (ps)", ylabel, oenv);
71     for (k = 0; (k < nf); k++)
72     {
73         fprintf(fp, "%10g  %10g\n", time[k], data[k]);
74     }
75     xvgrclose(fp);
76 }
77
78 static int calc_RBbin(real phi, int gmx_unused multiplicity, real gmx_unused core_frac)
79 {
80     /* multiplicity and core_frac NOT used,
81      * just given to enable use of pt-to-fn in caller low_ana_dih_trans*/
82     static const real r30  = M_PI/6.0;
83     static const real r90  = M_PI/2.0;
84     static const real r150 = M_PI*5.0/6.0;
85
86     if ((phi < r30) && (phi > -r30))
87     {
88         return 1;
89     }
90     else if ((phi > -r150) && (phi < -r90))
91     {
92         return 2;
93     }
94     else if ((phi < r150) && (phi > r90))
95     {
96         return 3;
97     }
98     return 0;
99 }
100
101 static int calc_Nbin(real phi, int multiplicity, real core_frac)
102 {
103     static const real r360 = 360*DEG2RAD;
104     real              rot_width, core_width, core_offset, low, hi;
105     int               bin;
106     /* with multiplicity 3 and core_frac 0.5
107      * 0<g(-)<120, 120<t<240, 240<g(+)<360
108      * 0< bin0 < 30, 30<bin1<90, 90<bin0<150, 150<bin2<210, 210<bin0<270, 270<bin3<330, 330<bin0<360
109      * so with multiplicity 3, bin1 is core g(-), bin2 is core t, bin3 is
110        core g(+), bin0 is between rotamers */
111     if (phi < 0)
112     {
113         phi += r360;
114     }
115
116     rot_width   = 360/multiplicity;
117     core_width  = core_frac * rot_width;
118     core_offset = (rot_width - core_width)/2.0;
119     for (bin = 1; bin <= multiplicity; bin++)
120     {
121         low  = ((bin - 1) * rot_width ) + core_offset;
122         hi   = ((bin - 1) * rot_width ) + core_offset + core_width;
123         low *= DEG2RAD;
124         hi  *= DEG2RAD;
125         if ((phi > low) && (phi < hi))
126         {
127             return bin;
128         }
129     }
130     return 0;
131 }
132
133 void ana_dih_trans(const char *fn_trans, const char *fn_histo,
134                    real **dih, int nframes, int nangles,
135                    const char *grpname, real *time, gmx_bool bRb,
136                    const output_env_t oenv)
137 {
138     /* just a wrapper; declare extra args, then chuck away at end. */
139     int      maxchi = 0;
140     t_dlist *dlist;
141     int     *multiplicity;
142     int      nlist = nangles;
143     int      k;
144
145     snew(dlist, nlist);
146     snew(multiplicity, nangles);
147     for (k = 0; (k < nangles); k++)
148     {
149         multiplicity[k] = 3;
150     }
151
152     low_ana_dih_trans(TRUE, fn_trans, TRUE, fn_histo, maxchi,
153                       dih, nlist, dlist, nframes,
154                       nangles, grpname, multiplicity, time, bRb, 0.5, oenv);
155     sfree(dlist);
156     sfree(multiplicity);
157
158 }
159
160 void low_ana_dih_trans(gmx_bool bTrans, const char *fn_trans,
161                        gmx_bool bHisto, const char *fn_histo, int maxchi,
162                        real **dih, int nlist, t_dlist dlist[], int nframes,
163                        int nangles, const char *grpname, int multiplicity[],
164                        real *time, gmx_bool bRb, real core_frac,
165                        const output_env_t oenv)
166 {
167     FILE *fp;
168     int  *tr_f, *tr_h;
169     char  title[256];
170     int   i, j, k, Dih, ntrans;
171     int   cur_bin, new_bin;
172     real  ttime;
173     real *rot_occ[NROT];
174     int   (*calc_bin)(real, int, real);
175     real  dt;
176
177     if (1 <= nframes)
178     {
179         return;
180     }
181     /* Assumes the frames are equally spaced in time */
182     dt = (time[nframes-1]-time[0])/(nframes-1);
183
184     /* Analysis of dihedral transitions */
185     fprintf(stderr, "Now calculating transitions...\n");
186
187     if (bRb)
188     {
189         calc_bin = calc_RBbin;
190     }
191     else
192     {
193         calc_bin = calc_Nbin;
194     }
195
196     for (k = 0; k < NROT; k++)
197     {
198         snew(rot_occ[k], nangles);
199         for (i = 0; (i < nangles); i++)
200         {
201             rot_occ[k][i] = 0;
202         }
203     }
204     snew(tr_h, nangles);
205     snew(tr_f, nframes);
206
207     /* dih[i][j] is the dihedral angle i in frame j  */
208     ntrans = 0;
209     for (i = 0; (i < nangles); i++)
210     {
211
212         /*#define OLDIE*/
213 #ifdef OLDIE
214         mind = maxd = prev = dih[i][0];
215 #else
216         cur_bin = calc_bin(dih[i][0], multiplicity[i], core_frac);
217         rot_occ[cur_bin][i]++;
218 #endif
219         for (j = 1; (j < nframes); j++)
220         {
221             new_bin = calc_bin(dih[i][j], multiplicity[i], core_frac);
222             rot_occ[new_bin][i]++;
223 #ifndef OLDIE
224             if (cur_bin == 0)
225             {
226                 cur_bin = new_bin;
227             }
228             else if ((new_bin != 0) && (cur_bin != new_bin))
229             {
230                 cur_bin = new_bin;
231                 tr_f[j]++;
232                 tr_h[i]++;
233                 ntrans++;
234             }
235 #else
236             /* why is all this md rubbish periodic? Remove 360 degree periodicity */
237             if ( (dih[i][j] - prev) > M_PI)
238             {
239                 dih[i][j] -= 2*M_PI;
240             }
241             else if ( (dih[i][j] - prev) < -M_PI)
242             {
243                 dih[i][j] += 2*M_PI;
244             }
245
246             prev = dih[i][j];
247
248             mind = std::min(mind, dih[i][j]);
249             maxd = std::max(maxd, dih[i][j]);
250             if ( (maxd - mind) > 2*M_PI/3) /* or 120 degrees, assuming       */
251             {                              /* multiplicity 3. Not so general.*/
252                 tr_f[j]++;
253                 tr_h[i]++;
254                 maxd = mind = dih[i][j]; /* get ready for next transition  */
255                 ntrans++;
256             }
257 #endif
258         } /* end j */
259         for (k = 0; k < NROT; k++)
260         {
261             rot_occ[k][i] /= nframes;
262         }
263     } /* end i */
264     fprintf(stderr, "Total number of transitions: %10d\n", ntrans);
265     if (ntrans > 0)
266     {
267         ttime = (dt*nframes*nangles)/ntrans;
268         fprintf(stderr, "Time between transitions:    %10.3f ps\n", ttime);
269     }
270
271     /* new by grs - copy transitions from tr_h[] to dlist->ntr[]
272      * and rotamer populations from rot_occ to dlist->rot_occ[]
273      * based on fn histogramming in g_chi. diff roles for i and j here */
274
275     j = 0;
276     for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
277     {
278         for (i = 0; (i < nlist); i++)
279         {
280             if (((Dih  < edOmega) ) ||
281                 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
282                 ((Dih  > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)))
283             {
284                 /* grs debug  printf("Not OK? i %d j %d Dih %d \n", i, j, Dih) ; */
285                 dlist[i].ntr[Dih] = tr_h[j];
286                 for (k = 0; k < NROT; k++)
287                 {
288                     dlist[i].rot_occ[Dih][k] = rot_occ[k][j];
289                 }
290                 j++;
291             }
292         }
293     }
294
295     /* end addition by grs */
296
297     if (bTrans)
298     {
299         sprintf(title, "Number of transitions: %s", grpname);
300         fp = xvgropen(fn_trans, title, "Time (ps)", "# transitions/timeframe", oenv);
301         for (j = 0; (j < nframes); j++)
302         {
303             fprintf(fp, "%10.3f  %10d\n", time[j], tr_f[j]);
304         }
305         xvgrclose(fp);
306     }
307
308     /* Compute histogram from # transitions per dihedral */
309     /* Use old array */
310     for (j = 0; (j < nframes); j++)
311     {
312         tr_f[j] = 0;
313     }
314     for (i = 0; (i < nangles); i++)
315     {
316         tr_f[tr_h[i]]++;
317     }
318     for (j = nframes; ((tr_f[j-1] == 0) && (j > 0)); j--)
319     {
320         ;
321     }
322
323     ttime = dt*nframes;
324     if (bHisto)
325     {
326         sprintf(title, "Transition time: %s", grpname);
327         fp = xvgropen(fn_histo, title, "Time (ps)", "#", oenv);
328         for (i = j-1; (i > 0); i--)
329         {
330             if (tr_f[i] != 0)
331             {
332                 fprintf(fp, "%10.3f  %10d\n", ttime/i, tr_f[i]);
333             }
334         }
335         xvgrclose(fp);
336     }
337
338     sfree(tr_f);
339     sfree(tr_h);
340     for (k = 0; k < NROT; k++)
341     {
342         sfree(rot_occ[k]);
343     }
344
345 }
346
347 void mk_multiplicity_lookup (int *multiplicity, int maxchi,
348                              int nlist, t_dlist dlist[], int nangles)
349 {
350     /* new by grs - for dihedral j (as in dih[j]) get multiplicity from dlist
351      * and store in multiplicity[j]
352      */
353
354     int  j, Dih, i;
355     char name[4];
356
357     j = 0;
358     for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
359     {
360         for (i = 0; (i < nlist); i++)
361         {
362             std::strncpy(name, dlist[i].name, 3);
363             name[3] = '\0';
364             if (((Dih  < edOmega) ) ||
365                 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
366                 ((Dih  > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)))
367             {
368                 /* default - we will correct the rest below */
369                 multiplicity[j] = 3;
370
371                 /* make omegas 2fold, though doesn't make much more sense than 3 */
372                 if (Dih == edOmega && (has_dihedral(edOmega, &(dlist[i]))))
373                 {
374                     multiplicity[j] = 2;
375                 }
376
377                 /* dihedrals to aromatic rings, COO, CONH2 or guanidinium are 2fold*/
378                 if (Dih > edOmega && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1))
379                 {
380                     if ( ((std::strstr(name, "PHE") != NULL) && (Dih == edChi2))  ||
381                          ((std::strstr(name, "TYR") != NULL) && (Dih == edChi2))  ||
382                          ((std::strstr(name, "PTR") != NULL) && (Dih == edChi2))  ||
383                          ((std::strstr(name, "TRP") != NULL) && (Dih == edChi2))  ||
384                          ((std::strstr(name, "HIS") != NULL) && (Dih == edChi2))  ||
385                          ((std::strstr(name, "GLU") != NULL) && (Dih == edChi3))  ||
386                          ((std::strstr(name, "ASP") != NULL) && (Dih == edChi2))  ||
387                          ((std::strstr(name, "GLN") != NULL) && (Dih == edChi3))  ||
388                          ((std::strstr(name, "ASN") != NULL) && (Dih == edChi2))  ||
389                          ((std::strstr(name, "ARG") != NULL) && (Dih == edChi4))  )
390                     {
391                         multiplicity[j] = 2;
392                     }
393                 }
394                 j++;
395             }
396         }
397     }
398     if (j < nangles)
399     {
400         fprintf(stderr, "WARNING: not all dihedrals found in topology (only %d out of %d)!\n",
401                 j, nangles);
402     }
403     /* Check for remaining dihedrals */
404     for (; (j < nangles); j++)
405     {
406         multiplicity[j] = 3;
407     }
408
409 }
410
411 void mk_chi_lookup (int **lookup, int maxchi,
412                     int nlist, t_dlist dlist[])
413 {
414
415     /* by grs. should rewrite everything to use this. (but haven't,
416      * and at mmt only used in get_chi_product_traj
417      * returns the dihed number given the residue number (from-0)
418      * and chi (from-0) nr. -1 for chi undefined for that res (eg gly, ala..)*/
419
420     int i, j, Dih, Chi;
421
422     j = 0;
423     /* NONCHI points to chi1, therefore we have to start counting there. */
424     for (Dih = NONCHI; (Dih < NONCHI+maxchi); Dih++)
425     {
426         for (i = 0; (i < nlist); i++)
427         {
428             Chi = Dih - NONCHI;
429             if (((Dih  < edOmega) ) ||
430                 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
431                 ((Dih  > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)))
432             {
433                 /* grs debug  printf("Not OK? i %d j %d Dih %d \n", i, j, Dih) ; */
434                 if (Dih > edOmega)
435                 {
436                     lookup[i][Chi] = j;
437                 }
438                 j++;
439             }
440             else
441             {
442                 lookup[i][Chi] = -1;
443             }
444         }
445     }
446
447 }
448
449
450 void get_chi_product_traj (real **dih, int nframes, int nlist,
451                            int maxchi, t_dlist dlist[], real time[],
452                            int **lookup, int *multiplicity, gmx_bool bRb, gmx_bool bNormalize,
453                            real core_frac, gmx_bool bAll, const char *fnall,
454                            const output_env_t oenv)
455 {
456
457     gmx_bool bRotZero, bHaveChi = FALSE;
458     int      accum = 0, index, i, j, k, Xi, n, b;
459     real    *chi_prtrj;
460     int     *chi_prhist;
461     int      nbin;
462     FILE    *fp, *fpall;
463     char     hisfile[256], histitle[256], *namept;
464
465     int      (*calc_bin)(real, int, real);
466
467     /* Analysis of dihedral transitions */
468     fprintf(stderr, "Now calculating Chi product trajectories...\n");
469
470     if (bRb)
471     {
472         calc_bin = calc_RBbin;
473     }
474     else
475     {
476         calc_bin = calc_Nbin;
477     }
478
479     snew(chi_prtrj, nframes);
480
481     /* file for info on all residues */
482     if (bNormalize)
483     {
484         fpall = xvgropen(fnall, "Cumulative Rotamers", "Residue", "Probability", oenv);
485     }
486     else
487     {
488         fpall = xvgropen(fnall, "Cumulative Rotamers", "Residue", "# Counts", oenv);
489     }
490
491     for (i = 0; (i < nlist); i++)
492     {
493
494         /* get nbin, the nr. of cumulative rotamers that need to be considered */
495         nbin = 1;
496         for (Xi = 0; Xi < maxchi; Xi++)
497         {
498             index = lookup[i][Xi]; /* chi_(Xi+1) of res i (-1 if off end) */
499             if (index >= 0)
500             {
501                 n    = multiplicity[index];
502                 nbin = n*nbin;
503             }
504         }
505         nbin += 1; /* for the "zero rotamer", outside the core region */
506
507         for (j = 0; (j < nframes); j++)
508         {
509
510             bRotZero = FALSE;
511             bHaveChi = TRUE;
512             index    = lookup[i][0]; /* index into dih of chi1 of res i */
513             if (index == -1)
514             {
515                 bRotZero = TRUE;
516                 bHaveChi = FALSE;
517             }
518             else
519             {
520                 b     = calc_bin(dih[index][j], multiplicity[index], core_frac);
521                 accum = b - 1;
522                 if (b == 0)
523                 {
524                     bRotZero = TRUE;
525                 }
526                 for (Xi = 1; Xi < maxchi; Xi++)
527                 {
528                     index = lookup[i][Xi]; /* chi_(Xi+1) of res i (-1 if off end) */
529                     if (index >= 0)
530                     {
531                         n     = multiplicity[index];
532                         b     = calc_bin(dih[index][j], n, core_frac);
533                         accum = n * accum + b - 1;
534                         if (b == 0)
535                         {
536                             bRotZero = TRUE;
537                         }
538                     }
539                 }
540                 accum++;
541             }
542             if (bRotZero)
543             {
544                 chi_prtrj[j] = 0.0;
545             }
546             else
547             {
548                 chi_prtrj[j] = accum;
549                 if (accum+1 > nbin)
550                 {
551                     nbin = accum+1;
552                 }
553             }
554         }
555         if (bHaveChi)
556         {
557
558             if (bAll)
559             {
560                 /* print cuml rotamer vs time */
561                 print_one(oenv, "chiproduct", dlist[i].name, "chi product for",
562                           "cumulative rotamer", nframes, time, chi_prtrj);
563             }
564
565             /* make a histogram pf culm. rotamer occupancy too */
566             snew(chi_prhist, nbin);
567             make_histo(NULL, nframes, chi_prtrj, nbin, chi_prhist, 0, nbin);
568             if (bAll)
569             {
570                 sprintf(hisfile, "histo-chiprod%s.xvg", dlist[i].name);
571                 sprintf(histitle, "cumulative rotamer distribution for %s", dlist[i].name);
572                 fprintf(stderr, "  and %s  ", hisfile);
573                 fp = xvgropen(hisfile, histitle, "number", "", oenv);
574                 if (output_env_get_print_xvgr_codes(oenv))
575                 {
576                     fprintf(fp, "@ xaxis tick on\n");
577                     fprintf(fp, "@ xaxis tick major 1\n");
578                     fprintf(fp, "@ type xy\n");
579                 }
580                 for (k = 0; (k < nbin); k++)
581                 {
582                     if (bNormalize)
583                     {
584                         fprintf(fp, "%5d  %10g\n", k, (1.0*chi_prhist[k])/nframes);
585                     }
586                     else
587                     {
588                         fprintf(fp, "%5d  %10d\n", k, chi_prhist[k]);
589                     }
590                 }
591                 fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
592                 xvgrclose(fp);
593             }
594
595             /* and finally print out occupancies to a single file */
596             /* get the gmx from-1 res nr by setting a ptr to the number part
597              * of dlist[i].name - potential bug for 4-letter res names... */
598             namept = dlist[i].name + 3;
599             fprintf(fpall, "%5s ", namept);
600             for (k = 0; (k < nbin); k++)
601             {
602                 if (bNormalize)
603                 {
604                     fprintf(fpall, "  %10g", (1.0*chi_prhist[k])/nframes);
605                 }
606                 else
607                 {
608                     fprintf(fpall, "  %10d", chi_prhist[k]);
609                 }
610             }
611             fprintf(fpall, "\n");
612
613             sfree(chi_prhist);
614             /* histogram done */
615         }
616     }
617
618     sfree(chi_prtrj);
619     xvgrclose(fpall);
620     fprintf(stderr, "\n");
621
622 }
623
624 void calc_distribution_props(int nh, int histo[], real start,
625                              int nkkk, t_karplus kkk[],
626                              real *S2)
627 {
628     real d, dc, ds, c1, c2, tdc, tds;
629     real fac, ang, invth, Jc;
630     int  i, j, th;
631
632     if (nh == 0)
633     {
634         gmx_fatal(FARGS, "No points in histogram (%s, %d)", __FILE__, __LINE__);
635     }
636     fac = 2*M_PI/nh;
637
638     /* Compute normalisation factor */
639     th = 0;
640     for (j = 0; (j < nh); j++)
641     {
642         th += histo[j];
643     }
644     invth = 1.0/th;
645
646     for (i = 0; (i < nkkk); i++)
647     {
648         kkk[i].Jc    = 0;
649         kkk[i].Jcsig = 0;
650     }
651     tdc = 0, tds = 0;
652     for (j = 0; (j < nh); j++)
653     {
654         d    = invth*histo[j];
655         ang  = j*fac-start;
656         c1   = std::cos(ang);
657         dc   = d*c1;
658         ds   = d*std::sin(ang);
659         tdc += dc;
660         tds += ds;
661         for (i = 0; (i < nkkk); i++)
662         {
663             c1            = std::cos(ang+kkk[i].offset);
664             c2            = c1*c1;
665             Jc            = (kkk[i].A*c2 + kkk[i].B*c1 + kkk[i].C);
666             kkk[i].Jc    += histo[j]*Jc;
667             kkk[i].Jcsig += histo[j]*sqr(Jc);
668         }
669     }
670     for (i = 0; (i < nkkk); i++)
671     {
672         kkk[i].Jc    /= th;
673         kkk[i].Jcsig  = std::sqrt(kkk[i].Jcsig/th-sqr(kkk[i].Jc));
674     }
675     *S2 = tdc*tdc+tds*tds;
676 }
677
678 static void calc_angles(struct t_pbc *pbc,
679                         int n3, atom_id index[], real ang[], rvec x_s[])
680 {
681     int  i, ix, t1, t2;
682     rvec r_ij, r_kj;
683     real costh = 0.0;
684
685     for (i = ix = 0; (ix < n3); i++, ix += 3)
686     {
687         ang[i] = bond_angle(x_s[index[ix]], x_s[index[ix+1]], x_s[index[ix+2]],
688                             pbc, r_ij, r_kj, &costh, &t1, &t2);
689     }
690     if (debug)
691     {
692         fprintf(debug, "Angle[0]=%g, costh=%g, index0 = %d, %d, %d\n",
693                 ang[0], costh, index[0], index[1], index[2]);
694         pr_rvec(debug, 0, "rij", r_ij, DIM, TRUE);
695         pr_rvec(debug, 0, "rkj", r_kj, DIM, TRUE);
696     }
697 }
698
699 static real calc_fraction(real angles[], int nangles)
700 {
701     int  i;
702     real trans = 0, gauche = 0;
703     real angle;
704
705     for (i = 0; i < nangles; i++)
706     {
707         angle = angles[i] * RAD2DEG;
708
709         if (angle > 135 && angle < 225)
710         {
711             trans += 1.0;
712         }
713         else if (angle > 270 && angle < 330)
714         {
715             gauche += 1.0;
716         }
717         else if (angle < 90 && angle > 30)
718         {
719             gauche += 1.0;
720         }
721     }
722     if (trans+gauche > 0)
723     {
724         return trans/(trans+gauche);
725     }
726     else
727     {
728         return 0;
729     }
730 }
731
732 static void calc_dihs(struct t_pbc *pbc,
733                       int n4, atom_id index[], real ang[], rvec x_s[])
734 {
735     int  i, ix, t1, t2, t3;
736     rvec r_ij, r_kj, r_kl, m, n;
737     real sign, aaa;
738
739     for (i = ix = 0; (ix < n4); i++, ix += 4)
740     {
741         aaa = dih_angle(x_s[index[ix]], x_s[index[ix+1]], x_s[index[ix+2]],
742                         x_s[index[ix+3]], pbc,
743                         r_ij, r_kj, r_kl, m, n,
744                         &sign, &t1, &t2, &t3);
745
746         ang[i] = aaa; /* not taking into account ryckaert bellemans yet */
747     }
748 }
749
750 void make_histo(FILE *log,
751                 int ndata, real data[], int npoints, int histo[],
752                 real minx, real maxx)
753 {
754     double dx;
755     int    i, ind;
756
757     if (minx == maxx)
758     {
759         minx = maxx = data[0];
760         for (i = 1; (i < ndata); i++)
761         {
762             minx = std::min(minx, data[i]);
763             maxx = std::max(maxx, data[i]);
764         }
765         fprintf(log, "Min data: %10g  Max data: %10g\n", minx, maxx);
766     }
767     dx = (double)npoints/(maxx-minx);
768     if (debug)
769     {
770         fprintf(debug,
771                 "Histogramming: ndata=%d, nhisto=%d, minx=%g,maxx=%g,dx=%g\n",
772                 ndata, npoints, minx, maxx, dx);
773     }
774     for (i = 0; (i < ndata); i++)
775     {
776         ind = static_cast<int>((data[i]-minx)*dx);
777         if ((ind >= 0) && (ind < npoints))
778         {
779             histo[ind]++;
780         }
781         else
782         {
783             fprintf(log, "index = %d, data[%d] = %g\n", ind, i, data[i]);
784         }
785     }
786 }
787
788 void normalize_histo(int npoints, int histo[], real dx, real normhisto[])
789 {
790     int    i;
791     double d, fac;
792
793     d = 0;
794     for (i = 0; (i < npoints); i++)
795     {
796         d += dx*histo[i];
797     }
798     if (d == 0)
799     {
800         fprintf(stderr, "Empty histogram!\n");
801         return;
802     }
803     fac = 1.0/d;
804     for (i = 0; (i < npoints); i++)
805     {
806         normhisto[i] = fac*histo[i];
807     }
808 }
809
810 void read_ang_dih(const char *trj_fn,
811                   gmx_bool bAngles, gmx_bool bSaveAll, gmx_bool bRb, gmx_bool bPBC,
812                   int maxangstat, int angstat[],
813                   int *nframes, real **time,
814                   int isize, atom_id index[],
815                   real **trans_frac,
816                   real **aver_angle,
817                   real *dih[],
818                   const output_env_t oenv)
819 {
820     struct t_pbc *pbc;
821     t_trxstatus  *status;
822     int           i, angind, total, teller;
823     int           nangles, n_alloc;
824     real          t, fraction, pifac, aa, angle;
825     real         *angles[2];
826     matrix        box;
827     rvec         *x;
828     int           cur = 0;
829 #define prev (1-cur)
830
831     snew(pbc, 1);
832     read_first_x(oenv, &status, trj_fn, &t, &x, box);
833
834     if (bAngles)
835     {
836         nangles = isize/3;
837         pifac   = M_PI;
838     }
839     else
840     {
841         nangles = isize/4;
842         pifac   = 2.0*M_PI;
843     }
844     snew(angles[cur], nangles);
845     snew(angles[prev], nangles);
846
847     /* Start the loop over frames */
848     total       = 0;
849     teller      = 0;
850     n_alloc     = 0;
851     *time       = NULL;
852     *trans_frac = NULL;
853     *aver_angle = NULL;
854
855     do
856     {
857         if (teller >= n_alloc)
858         {
859             n_alloc += 100;
860             if (bSaveAll)
861             {
862                 for (i = 0; (i < nangles); i++)
863                 {
864                     srenew(dih[i], n_alloc);
865                 }
866             }
867             srenew(*time, n_alloc);
868             srenew(*trans_frac, n_alloc);
869             srenew(*aver_angle, n_alloc);
870         }
871
872         (*time)[teller] = t;
873
874         if (pbc)
875         {
876             set_pbc(pbc, -1, box);
877         }
878
879         if (bAngles)
880         {
881             calc_angles(pbc, isize, index, angles[cur], x);
882         }
883         else
884         {
885             calc_dihs(pbc, isize, index, angles[cur], x);
886
887             /* Trans fraction */
888             fraction              = calc_fraction(angles[cur], nangles);
889             (*trans_frac)[teller] = fraction;
890
891             /* Change Ryckaert-Bellemans dihedrals to polymer convention
892              * Modified 990913 by Erik:
893              * We actually shouldn't change the convention, since it's
894              * calculated from polymer above, but we change the intervall
895              * from [-180,180] to [0,360].
896              */
897             if (bRb)
898             {
899                 for (i = 0; (i < nangles); i++)
900                 {
901                     if (angles[cur][i] <= 0.0)
902                     {
903                         angles[cur][i] += 2*M_PI;
904                     }
905                 }
906             }
907
908             /* Periodicity in dihedral space... */
909             if (bPBC)
910             {
911                 for (i = 0; (i < nangles); i++)
912                 {
913                     real dd = angles[cur][i];
914                     angles[cur][i] = std::atan2(sin(dd), cos(dd));
915                 }
916             }
917             else
918             {
919                 if (teller > 1)
920                 {
921                     for (i = 0; (i < nangles); i++)
922                     {
923                         while (angles[cur][i] <= angles[prev][i] - M_PI)
924                         {
925                             angles[cur][i] += 2*M_PI;
926                         }
927                         while (angles[cur][i] > angles[prev][i] + M_PI)
928                         {
929                             angles[cur][i] -= 2*M_PI;
930                         }
931                     }
932                 }
933             }
934         }
935
936         /* Average angles */
937         aa = 0;
938         for (i = 0; (i < nangles); i++)
939         {
940             aa = aa+angles[cur][i];
941
942             /* angle in rad / 2Pi * max determines bin. bins go from 0 to maxangstat,
943                even though scale goes from -pi to pi (dihedral) or -pi/2 to pi/2
944                (angle) Basically: translate the x-axis by Pi. Translate it back by
945                -Pi when plotting.
946              */
947
948             angle = angles[cur][i];
949             if (!bAngles)
950             {
951                 while (angle < -M_PI)
952                 {
953                     angle += 2*M_PI;
954                 }
955                 while (angle >= M_PI)
956                 {
957                     angle -= 2*M_PI;
958                 }
959
960                 angle += M_PI;
961             }
962
963             /* Update the distribution histogram */
964             angind = static_cast<int>((angle*maxangstat)/pifac + 0.5);
965             if (angind == maxangstat)
966             {
967                 angind = 0;
968             }
969             if ( (angind < 0) || (angind >= maxangstat) )
970             {
971                 /* this will never happen */
972                 gmx_fatal(FARGS, "angle (%f) index out of range (0..%d) : %d\n",
973                           angle, maxangstat, angind);
974             }
975
976             angstat[angind]++;
977             if (angind == maxangstat)
978             {
979                 fprintf(stderr, "angle %d fr %d = %g\n", i, cur, angle);
980             }
981
982             total++;
983         }
984
985         /* average over all angles */
986         (*aver_angle)[teller] = (aa/nangles);
987
988         /* this copies all current dih. angles to dih[i], teller is frame */
989         if (bSaveAll)
990         {
991             for (i = 0; i < nangles; i++)
992             {
993                 dih[i][teller] = angles[cur][i];
994             }
995         }
996
997         /* Swap buffers */
998         cur = prev;
999
1000         /* Increment loop counter */
1001         teller++;
1002     }
1003     while (read_next_x(oenv, status, &t, x, box));
1004     close_trj(status);
1005
1006     sfree(x);
1007     sfree(angles[cur]);
1008     sfree(angles[prev]);
1009
1010     *nframes = teller;
1011 }