src/gromacs/gmxana/anadih.c

   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, 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 "config.h"
  40
  41 #include <math.h>
  42 #include <stdio.h>
  43 #include <string.h>
  44 #include "gromacs/math/units.h"
  45 #include "gromacs/utility/smalloc.h"
  46 #include "gromacs/legacyheaders/macros.h"
  47 #include "gromacs/legacyheaders/txtdump.h"
  48 #include "gromacs/fileio/xvgr.h"
  49 #include "gromacs/legacyheaders/typedefs.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gstat.h"
  52 #include "gromacs/fileio/confio.h"
  53 #include "gromacs/fileio/trxio.h"
  54
  55 #include "gromacs/bonded/bonded.h"
  56 #include "gromacs/pbcutil/pbc.h"
  57 #include "gromacs/utility/fatalerror.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     gmx_ffclose(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, tt;
 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 = min(mind, dih[i][j]);
 249             maxd = 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         gmx_ffclose(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         gmx_ffclose(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             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 ( ((strstr(name, "PHE") != NULL) && (Dih == edChi2))  ||
 381                          ((strstr(name, "TYR") != NULL) && (Dih == edChi2))  ||
 382                          ((strstr(name, "PTR") != NULL) && (Dih == edChi2))  ||
 383                          ((strstr(name, "TRP") != NULL) && (Dih == edChi2))  ||
 384                          ((strstr(name, "HIS") != NULL) && (Dih == edChi2))  ||
 385                          ((strstr(name, "GLU") != NULL) && (Dih == edChi3))  ||
 386                          ((strstr(name, "ASP") != NULL) && (Dih == edChi2))  ||
 387                          ((strstr(name, "GLN") != NULL) && (Dih == edChi3))  ||
 388                          ((strstr(name, "ASN") != NULL) && (Dih == edChi2))  ||
 389                          ((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                 b        = 0;
 516                 bRotZero = TRUE;
 517                 bHaveChi = FALSE;
 518             }
 519             else
 520             {
 521                 b     = calc_bin(dih[index][j], multiplicity[index], core_frac);
 522                 accum = b - 1;
 523                 if (b == 0)
 524                 {
 525                     bRotZero = TRUE;
 526                 }
 527                 for (Xi = 1; Xi < maxchi; Xi++)
 528                 {
 529                     index = lookup[i][Xi]; /* chi_(Xi+1) of res i (-1 if off end) */
 530                     if (index >= 0)
 531                     {
 532                         n     = multiplicity[index];
 533                         b     = calc_bin(dih[index][j], n, core_frac);
 534                         accum = n * accum + b - 1;
 535                         if (b == 0)
 536                         {
 537                             bRotZero = TRUE;
 538                         }
 539                     }
 540                 }
 541                 accum++;
 542             }
 543             if (bRotZero)
 544             {
 545                 chi_prtrj[j] = 0.0;
 546             }
 547             else
 548             {
 549                 chi_prtrj[j] = accum;
 550                 if (accum+1 > nbin)
 551                 {
 552                     nbin = accum+1;
 553                 }
 554             }
 555         }
 556         if (bHaveChi)
 557         {
 558
 559             if (bAll)
 560             {
 561                 /* print cuml rotamer vs time */
 562                 print_one(oenv, "chiproduct", dlist[i].name, "chi product for",
 563                           "cumulative rotamer", nframes, time, chi_prtrj);
 564             }
 565
 566             /* make a histogram pf culm. rotamer occupancy too */
 567             snew(chi_prhist, nbin);
 568             make_histo(NULL, nframes, chi_prtrj, nbin, chi_prhist, 0, nbin);
 569             if (bAll)
 570             {
 571                 sprintf(hisfile, "histo-chiprod%s.xvg", dlist[i].name);
 572                 sprintf(histitle, "cumulative rotamer distribution for %s", dlist[i].name);
 573                 fprintf(stderr, "  and %s  ", hisfile);
 574                 fp = xvgropen(hisfile, histitle, "number", "", oenv);
 575                 if (output_env_get_print_xvgr_codes(oenv))
 576                 {
 577                     fprintf(fp, "@ xaxis tick on\n");
 578                     fprintf(fp, "@ xaxis tick major 1\n");
 579                     fprintf(fp, "@ type xy\n");
 580                 }
 581                 for (k = 0; (k < nbin); k++)
 582                 {
 583                     if (bNormalize)
 584                     {
 585                         fprintf(fp, "%5d  %10g\n", k, (1.0*chi_prhist[k])/nframes);
 586                     }
 587                     else
 588                     {
 589                         fprintf(fp, "%5d  %10d\n", k, chi_prhist[k]);
 590                     }
 591                 }
 592                 fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 593                 gmx_ffclose(fp);
 594             }
 595
 596             /* and finally print out occupancies to a single file */
 597             /* get the gmx from-1 res nr by setting a ptr to the number part
 598              * of dlist[i].name - potential bug for 4-letter res names... */
 599             namept = dlist[i].name + 3;
 600             fprintf(fpall, "%5s ", namept);
 601             for (k = 0; (k < nbin); k++)
 602             {
 603                 if (bNormalize)
 604                 {
 605                     fprintf(fpall, "  %10g", (1.0*chi_prhist[k])/nframes);
 606                 }
 607                 else
 608                 {
 609                     fprintf(fpall, "  %10d", chi_prhist[k]);
 610                 }
 611             }
 612             fprintf(fpall, "\n");
 613
 614             sfree(chi_prhist);
 615             /* histogram done */
 616         }
 617     }
 618
 619     sfree(chi_prtrj);
 620     gmx_ffclose(fpall);
 621     fprintf(stderr, "\n");
 622
 623 }
 624
 625 void calc_distribution_props(int nh, int histo[], real start,
 626                              int nkkk, t_karplus kkk[],
 627                              real *S2)
 628 {
 629     real d, dc, ds, c1, c2, tdc, tds;
 630     real fac, ang, invth, Jc;
 631     int  i, j, th;
 632
 633     if (nh == 0)
 634     {
 635         gmx_fatal(FARGS, "No points in histogram (%s, %d)", __FILE__, __LINE__);
 636     }
 637     fac = 2*M_PI/nh;
 638
 639     /* Compute normalisation factor */
 640     th = 0;
 641     for (j = 0; (j < nh); j++)
 642     {
 643         th += histo[j];
 644     }
 645     invth = 1.0/th;
 646
 647     for (i = 0; (i < nkkk); i++)
 648     {
 649         kkk[i].Jc    = 0;
 650         kkk[i].Jcsig = 0;
 651     }
 652     tdc = 0, tds = 0;
 653     for (j = 0; (j < nh); j++)
 654     {
 655         d    = invth*histo[j];
 656         ang  = j*fac-start;
 657         c1   = cos(ang);
 658         c2   = c1*c1;
 659         dc   = d*c1;
 660         ds   = d*sin(ang);
 661         tdc += dc;
 662         tds += ds;
 663         for (i = 0; (i < nkkk); i++)
 664         {
 665             c1            = cos(ang+kkk[i].offset);
 666             c2            = c1*c1;
 667             Jc            = (kkk[i].A*c2 + kkk[i].B*c1 + kkk[i].C);
 668             kkk[i].Jc    += histo[j]*Jc;
 669             kkk[i].Jcsig += histo[j]*sqr(Jc);
 670         }
 671     }
 672     for (i = 0; (i < nkkk); i++)
 673     {
 674         kkk[i].Jc    /= th;
 675         kkk[i].Jcsig  = sqrt(kkk[i].Jcsig/th-sqr(kkk[i].Jc));
 676     }
 677     *S2 = tdc*tdc+tds*tds;
 678 }
 679
 680 static void calc_angles(struct t_pbc *pbc,
 681                         int n3, atom_id index[], real ang[], rvec x_s[])
 682 {
 683     int  i, ix, t1, t2;
 684     rvec r_ij, r_kj;
 685     real costh;
 686
 687     for (i = ix = 0; (ix < n3); i++, ix += 3)
 688     {
 689         ang[i] = bond_angle(x_s[index[ix]], x_s[index[ix+1]], x_s[index[ix+2]],
 690                             pbc, r_ij, r_kj, &costh, &t1, &t2);
 691     }
 692     if (debug)
 693     {
 694         fprintf(debug, "Angle[0]=%g, costh=%g, index0 = %d, %d, %d\n",
 695                 ang[0], costh, index[0], index[1], index[2]);
 696         pr_rvec(debug, 0, "rij", r_ij, DIM, TRUE);
 697         pr_rvec(debug, 0, "rkj", r_kj, DIM, TRUE);
 698     }
 699 }
 700
 701 static real calc_fraction(real angles[], int nangles)
 702 {
 703     int  i;
 704     real trans = 0, gauche = 0;
 705     real angle;
 706
 707     for (i = 0; i < nangles; i++)
 708     {
 709         angle = angles[i] * RAD2DEG;
 710
 711         if (angle > 135 && angle < 225)
 712         {
 713             trans += 1.0;
 714         }
 715         else if (angle > 270 && angle < 330)
 716         {
 717             gauche += 1.0;
 718         }
 719         else if (angle < 90 && angle > 30)
 720         {
 721             gauche += 1.0;
 722         }
 723     }
 724     if (trans+gauche > 0)
 725     {
 726         return trans/(trans+gauche);
 727     }
 728     else
 729     {
 730         return 0;
 731     }
 732 }
 733
 734 static void calc_dihs(struct t_pbc *pbc,
 735                       int n4, atom_id index[], real ang[], rvec x_s[])
 736 {
 737     int  i, ix, t1, t2, t3;
 738     rvec r_ij, r_kj, r_kl, m, n;
 739     real sign, aaa;
 740
 741     for (i = ix = 0; (ix < n4); i++, ix += 4)
 742     {
 743         aaa = dih_angle(x_s[index[ix]], x_s[index[ix+1]], x_s[index[ix+2]],
 744                         x_s[index[ix+3]], pbc,
 745                         r_ij, r_kj, r_kl, m, n,
 746                         &sign, &t1, &t2, &t3);
 747
 748         ang[i] = aaa; /* not taking into account ryckaert bellemans yet */
 749     }
 750 }
 751
 752 void make_histo(FILE *log,
 753                 int ndata, real data[], int npoints, int histo[],
 754                 real minx, real maxx)
 755 {
 756     double dx;
 757     int    i, ind;
 758
 759     if (minx == maxx)
 760     {
 761         minx = maxx = data[0];
 762         for (i = 1; (i < ndata); i++)
 763         {
 764             minx = min(minx, data[i]);
 765             maxx = max(maxx, data[i]);
 766         }
 767         fprintf(log, "Min data: %10g  Max data: %10g\n", minx, maxx);
 768     }
 769     dx = (double)npoints/(maxx-minx);
 770     if (debug)
 771     {
 772         fprintf(debug,
 773                 "Histogramming: ndata=%d, nhisto=%d, minx=%g,maxx=%g,dx=%g\n",
 774                 ndata, npoints, minx, maxx, dx);
 775     }
 776     for (i = 0; (i < ndata); i++)
 777     {
 778         ind = (data[i]-minx)*dx;
 779         if ((ind >= 0) && (ind < npoints))
 780         {
 781             histo[ind]++;
 782         }
 783         else
 784         {
 785             fprintf(log, "index = %d, data[%d] = %g\n", ind, i, data[i]);
 786         }
 787     }
 788 }
 789
 790 void normalize_histo(int npoints, int histo[], real dx, real normhisto[])
 791 {
 792     int    i;
 793     double d, fac;
 794
 795     d = 0;
 796     for (i = 0; (i < npoints); i++)
 797     {
 798         d += dx*histo[i];
 799     }
 800     if (d == 0)
 801     {
 802         fprintf(stderr, "Empty histogram!\n");
 803         return;
 804     }
 805     fac = 1.0/d;
 806     for (i = 0; (i < npoints); i++)
 807     {
 808         normhisto[i] = fac*histo[i];
 809     }
 810 }
 811
 812 void read_ang_dih(const char *trj_fn,
 813                   gmx_bool bAngles, gmx_bool bSaveAll, gmx_bool bRb, gmx_bool bPBC,
 814                   int maxangstat, int angstat[],
 815                   int *nframes, real **time,
 816                   int isize, atom_id index[],
 817                   real **trans_frac,
 818                   real **aver_angle,
 819                   real *dih[],
 820                   const output_env_t oenv)
 821 {
 822     struct t_pbc *pbc;
 823     t_trxstatus  *status;
 824     int           i, angind, natoms, total, teller;
 825     int           nangles, n_alloc;
 826     real          t, fraction, pifac, aa, angle;
 827     real         *angles[2];
 828     matrix        box;
 829     rvec         *x;
 830     int           cur = 0;
 831 #define prev (1-cur)
 832
 833     snew(pbc, 1);
 834     natoms = read_first_x(oenv, &status, trj_fn, &t, &x, box);
 835
 836     if (bAngles)
 837     {
 838         nangles = isize/3;
 839         pifac   = M_PI;
 840     }
 841     else
 842     {
 843         nangles = isize/4;
 844         pifac   = 2.0*M_PI;
 845     }
 846     snew(angles[cur], nangles);
 847     snew(angles[prev], nangles);
 848
 849     /* Start the loop over frames */
 850     total       = 0;
 851     teller      = 0;
 852     n_alloc     = 0;
 853     *time       = NULL;
 854     *trans_frac = NULL;
 855     *aver_angle = NULL;
 856
 857     do
 858     {
 859         if (teller >= n_alloc)
 860         {
 861             n_alloc += 100;
 862             if (bSaveAll)
 863             {
 864                 for (i = 0; (i < nangles); i++)
 865                 {
 866                     srenew(dih[i], n_alloc);
 867                 }
 868             }
 869             srenew(*time, n_alloc);
 870             srenew(*trans_frac, n_alloc);
 871             srenew(*aver_angle, n_alloc);
 872         }
 873
 874         (*time)[teller] = t;
 875
 876         if (pbc)
 877         {
 878             set_pbc(pbc, -1, box);
 879         }
 880
 881         if (bAngles)
 882         {
 883             calc_angles(pbc, isize, index, angles[cur], x);
 884         }
 885         else
 886         {
 887             calc_dihs(pbc, isize, index, angles[cur], x);
 888
 889             /* Trans fraction */
 890             fraction              = calc_fraction(angles[cur], nangles);
 891             (*trans_frac)[teller] = fraction;
 892
 893             /* Change Ryckaert-Bellemans dihedrals to polymer convention
 894              * Modified 990913 by Erik:
 895              * We actually shouldn't change the convention, since it's
 896              * calculated from polymer above, but we change the intervall
 897              * from [-180,180] to [0,360].
 898              */
 899             if (bRb)
 900             {
 901                 for (i = 0; (i < nangles); i++)
 902                 {
 903                     if (angles[cur][i] <= 0.0)
 904                     {
 905                         angles[cur][i] += 2*M_PI;
 906                     }
 907                 }
 908             }
 909
 910             /* Periodicity in dihedral space... */
 911             if (bPBC)
 912             {
 913                 for (i = 0; (i < nangles); i++)
 914                 {
 915                     real dd = angles[cur][i];
 916                     angles[cur][i] = atan2(sin(dd), cos(dd));
 917                 }
 918             }
 919             else
 920             {
 921                 if (teller > 1)
 922                 {
 923                     for (i = 0; (i < nangles); i++)
 924                     {
 925                         while (angles[cur][i] <= angles[prev][i] - M_PI)
 926                         {
 927                             angles[cur][i] += 2*M_PI;
 928                         }
 929                         while (angles[cur][i] > angles[prev][i] + M_PI)
 930                         {
 931                             angles[cur][i] -= 2*M_PI;
 932                         }
 933                     }
 934                 }
 935             }
 936         }
 937
 938         /* Average angles */
 939         aa = 0;
 940         for (i = 0; (i < nangles); i++)
 941         {
 942             aa = aa+angles[cur][i];
 943
 944             /* angle in rad / 2Pi * max determines bin. bins go from 0 to maxangstat,
 945                even though scale goes from -pi to pi (dihedral) or -pi/2 to pi/2
 946                (angle) Basically: translate the x-axis by Pi. Translate it back by
 947                -Pi when plotting.
 948              */
 949
 950             angle = angles[cur][i];
 951             if (!bAngles)
 952             {
 953                 while (angle < -M_PI)
 954                 {
 955                     angle += 2*M_PI;
 956                 }
 957                 while (angle >= M_PI)
 958                 {
 959                     angle -= 2*M_PI;
 960                 }
 961
 962                 angle += M_PI;
 963             }
 964
 965             /* Update the distribution histogram */
 966             angind = (int) ((angle*maxangstat)/pifac + 0.5);
 967             if (angind == maxangstat)
 968             {
 969                 angind = 0;
 970             }
 971             if ( (angind < 0) || (angind >= maxangstat) )
 972             {
 973                 /* this will never happen */
 974                 gmx_fatal(FARGS, "angle (%f) index out of range (0..%d) : %d\n",
 975                           angle, maxangstat, angind);
 976             }
 977
 978             angstat[angind]++;
 979             if (angind == maxangstat)
 980             {
 981                 fprintf(stderr, "angle %d fr %d = %g\n", i, cur, angle);
 982             }
 983
 984             total++;
 985         }
 986
 987         /* average over all angles */
 988         (*aver_angle)[teller] = (aa/nangles);
 989
 990         /* this copies all current dih. angles to dih[i], teller is frame */
 991         if (bSaveAll)
 992         {
 993             for (i = 0; i < nangles; i++)
 994             {
 995                 dih[i][teller] = angles[cur][i];
 996             }
 997         }
 998
 999         /* Swap buffers */
1000         cur = prev;
1001
1002         /* Increment loop counter */
1003         teller++;
1004     }
1005     while (read_next_x(oenv, status, &t, x, box));
1006     close_trj(status);
1007
1008     sfree(x);
1009     sfree(angles[cur]);
1010     sfree(angles[prev]);
1011
1012     *nframes = teller;
1013 }