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 <math.h>
  40 #include <stdio.h>
  41 #include <string.h>
  42
  43 #include "gromacs/bonded/bonded.h"
  44 #include "gromacs/fileio/confio.h"
  45 #include "gromacs/fileio/trxio.h"
  46 #include "gromacs/fileio/xvgr.h"
  47 #include "gromacs/gmxana/gstat.h"
  48 #include "gromacs/legacyheaders/macros.h"
  49 #include "gromacs/legacyheaders/txtdump.h"
  50 #include "gromacs/legacyheaders/typedefs.h"
  51 #include "gromacs/math/units.h"
  52 #include "gromacs/math/vec.h"
  53 #include "gromacs/pbcutil/pbc.h"
  54 #include "gromacs/utility/fatalerror.h"
  55 #include "gromacs/utility/smalloc.h"
  56
  57 void print_one(const output_env_t oenv, const char *base, const char *name,
  58                const char *title, const char *ylabel, int nf, real time[],
  59                real data[])
  60 {
  61     FILE *fp;
  62     char  buf[256], t2[256];
  63     int   k;
  64
  65     sprintf(buf, "%s%s.xvg", base, name);
  66     fprintf(stderr, "\rPrinting %s  ", buf);
  67     sprintf(t2, "%s %s", title, name);
  68     fp = xvgropen(buf, t2, "Time (ps)", ylabel, oenv);
  69     for (k = 0; (k < nf); k++)
  70     {
  71         fprintf(fp, "%10g  %10g\n", time[k], data[k]);
  72     }
  73     gmx_ffclose(fp);
  74 }
  75
  76 static int calc_RBbin(real phi, int gmx_unused multiplicity, real gmx_unused core_frac)
  77 {
  78     /* multiplicity and core_frac NOT used,
  79      * just given to enable use of pt-to-fn in caller low_ana_dih_trans*/
  80     static const real r30  = M_PI/6.0;
  81     static const real r90  = M_PI/2.0;
  82     static const real r150 = M_PI*5.0/6.0;
  83
  84     if ((phi < r30) && (phi > -r30))
  85     {
  86         return 1;
  87     }
  88     else if ((phi > -r150) && (phi < -r90))
  89     {
  90         return 2;
  91     }
  92     else if ((phi < r150) && (phi > r90))
  93     {
  94         return 3;
  95     }
  96     return 0;
  97 }
  98
  99 static int calc_Nbin(real phi, int multiplicity, real core_frac)
 100 {
 101     static const real r360 = 360*DEG2RAD;
 102     real              rot_width, core_width, core_offset, low, hi;
 103     int               bin;
 104     /* with multiplicity 3 and core_frac 0.5
 105      * 0<g(-)<120, 120<t<240, 240<g(+)<360
 106      * 0< bin0 < 30, 30<bin1<90, 90<bin0<150, 150<bin2<210, 210<bin0<270, 270<bin3<330, 330<bin0<360
 107      * so with multiplicity 3, bin1 is core g(-), bin2 is core t, bin3 is
 108        core g(+), bin0 is between rotamers */
 109     if (phi < 0)
 110     {
 111         phi += r360;
 112     }
 113
 114     rot_width   = 360/multiplicity;
 115     core_width  = core_frac * rot_width;
 116     core_offset = (rot_width - core_width)/2.0;
 117     for (bin = 1; bin <= multiplicity; bin++)
 118     {
 119         low  = ((bin - 1) * rot_width ) + core_offset;
 120         hi   = ((bin - 1) * rot_width ) + core_offset + core_width;
 121         low *= DEG2RAD;
 122         hi  *= DEG2RAD;
 123         if ((phi > low) && (phi < hi))
 124         {
 125             return bin;
 126         }
 127     }
 128     return 0;
 129 }
 130
 131 void ana_dih_trans(const char *fn_trans, const char *fn_histo,
 132                    real **dih, int nframes, int nangles,
 133                    const char *grpname, real *time, gmx_bool bRb,
 134                    const output_env_t oenv)
 135 {
 136     /* just a wrapper; declare extra args, then chuck away at end. */
 137     int      maxchi = 0;
 138     t_dlist *dlist;
 139     int     *multiplicity;
 140     int      nlist = nangles;
 141     int      k;
 142
 143     snew(dlist, nlist);
 144     snew(multiplicity, nangles);
 145     for (k = 0; (k < nangles); k++)
 146     {
 147         multiplicity[k] = 3;
 148     }
 149
 150     low_ana_dih_trans(TRUE, fn_trans, TRUE, fn_histo, maxchi,
 151                       dih, nlist, dlist, nframes,
 152                       nangles, grpname, multiplicity, time, bRb, 0.5, oenv);
 153     sfree(dlist);
 154     sfree(multiplicity);
 155
 156 }
 157
 158 void low_ana_dih_trans(gmx_bool bTrans, const char *fn_trans,
 159                        gmx_bool bHisto, const char *fn_histo, int maxchi,
 160                        real **dih, int nlist, t_dlist dlist[], int nframes,
 161                        int nangles, const char *grpname, int multiplicity[],
 162                        real *time, gmx_bool bRb, real core_frac,
 163                        const output_env_t oenv)
 164 {
 165     FILE *fp;
 166     int  *tr_f, *tr_h;
 167     char  title[256];
 168     int   i, j, k, Dih, ntrans;
 169     int   cur_bin, new_bin;
 170     real  ttime, tt;
 171     real *rot_occ[NROT];
 172     int   (*calc_bin)(real, int, real);
 173     real  dt;
 174
 175     if (1 <= nframes)
 176     {
 177         return;
 178     }
 179     /* Assumes the frames are equally spaced in time */
 180     dt = (time[nframes-1]-time[0])/(nframes-1);
 181
 182     /* Analysis of dihedral transitions */
 183     fprintf(stderr, "Now calculating transitions...\n");
 184
 185     if (bRb)
 186     {
 187         calc_bin = calc_RBbin;
 188     }
 189     else
 190     {
 191         calc_bin = calc_Nbin;
 192     }
 193
 194     for (k = 0; k < NROT; k++)
 195     {
 196         snew(rot_occ[k], nangles);
 197         for (i = 0; (i < nangles); i++)
 198         {
 199             rot_occ[k][i] = 0;
 200         }
 201     }
 202     snew(tr_h, nangles);
 203     snew(tr_f, nframes);
 204
 205     /* dih[i][j] is the dihedral angle i in frame j  */
 206     ntrans = 0;
 207     for (i = 0; (i < nangles); i++)
 208     {
 209
 210         /*#define OLDIE*/
 211 #ifdef OLDIE
 212         mind = maxd = prev = dih[i][0];
 213 #else
 214         cur_bin = calc_bin(dih[i][0], multiplicity[i], core_frac);
 215         rot_occ[cur_bin][i]++;
 216 #endif
 217         for (j = 1; (j < nframes); j++)
 218         {
 219             new_bin = calc_bin(dih[i][j], multiplicity[i], core_frac);
 220             rot_occ[new_bin][i]++;
 221 #ifndef OLDIE
 222             if (cur_bin == 0)
 223             {
 224                 cur_bin = new_bin;
 225             }
 226             else if ((new_bin != 0) && (cur_bin != new_bin))
 227             {
 228                 cur_bin = new_bin;
 229                 tr_f[j]++;
 230                 tr_h[i]++;
 231                 ntrans++;
 232             }
 233 #else
 234             /* why is all this md rubbish periodic? Remove 360 degree periodicity */
 235             if ( (dih[i][j] - prev) > M_PI)
 236             {
 237                 dih[i][j] -= 2*M_PI;
 238             }
 239             else if ( (dih[i][j] - prev) < -M_PI)
 240             {
 241                 dih[i][j] += 2*M_PI;
 242             }
 243
 244             prev = dih[i][j];
 245
 246             mind = min(mind, dih[i][j]);
 247             maxd = max(maxd, dih[i][j]);
 248             if ( (maxd - mind) > 2*M_PI/3) /* or 120 degrees, assuming       */
 249             {                              /* multiplicity 3. Not so general.*/
 250                 tr_f[j]++;
 251                 tr_h[i]++;
 252                 maxd = mind = dih[i][j]; /* get ready for next transition  */
 253                 ntrans++;
 254             }
 255 #endif
 256         } /* end j */
 257         for (k = 0; k < NROT; k++)
 258         {
 259             rot_occ[k][i] /= nframes;
 260         }
 261     } /* end i */
 262     fprintf(stderr, "Total number of transitions: %10d\n", ntrans);
 263     if (ntrans > 0)
 264     {
 265         ttime = (dt*nframes*nangles)/ntrans;
 266         fprintf(stderr, "Time between transitions:    %10.3f ps\n", ttime);
 267     }
 268
 269     /* new by grs - copy transitions from tr_h[] to dlist->ntr[]
 270      * and rotamer populations from rot_occ to dlist->rot_occ[]
 271      * based on fn histogramming in g_chi. diff roles for i and j here */
 272
 273     j = 0;
 274     for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
 275     {
 276         for (i = 0; (i < nlist); i++)
 277         {
 278             if (((Dih  < edOmega) ) ||
 279                 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
 280                 ((Dih  > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)))
 281             {
 282                 /* grs debug  printf("Not OK? i %d j %d Dih %d \n", i, j, Dih) ; */
 283                 dlist[i].ntr[Dih] = tr_h[j];
 284                 for (k = 0; k < NROT; k++)
 285                 {
 286                     dlist[i].rot_occ[Dih][k] = rot_occ[k][j];
 287                 }
 288                 j++;
 289             }
 290         }
 291     }
 292
 293     /* end addition by grs */
 294
 295     if (bTrans)
 296     {
 297         sprintf(title, "Number of transitions: %s", grpname);
 298         fp = xvgropen(fn_trans, title, "Time (ps)", "# transitions/timeframe", oenv);
 299         for (j = 0; (j < nframes); j++)
 300         {
 301             fprintf(fp, "%10.3f  %10d\n", time[j], tr_f[j]);
 302         }
 303         gmx_ffclose(fp);
 304     }
 305
 306     /* Compute histogram from # transitions per dihedral */
 307     /* Use old array */
 308     for (j = 0; (j < nframes); j++)
 309     {
 310         tr_f[j] = 0;
 311     }
 312     for (i = 0; (i < nangles); i++)
 313     {
 314         tr_f[tr_h[i]]++;
 315     }
 316     for (j = nframes; ((tr_f[j-1] == 0) && (j > 0)); j--)
 317     {
 318         ;
 319     }
 320
 321     ttime = dt*nframes;
 322     if (bHisto)
 323     {
 324         sprintf(title, "Transition time: %s", grpname);
 325         fp = xvgropen(fn_histo, title, "Time (ps)", "#", oenv);
 326         for (i = j-1; (i > 0); i--)
 327         {
 328             if (tr_f[i] != 0)
 329             {
 330                 fprintf(fp, "%10.3f  %10d\n", ttime/i, tr_f[i]);
 331             }
 332         }
 333         gmx_ffclose(fp);
 334     }
 335
 336     sfree(tr_f);
 337     sfree(tr_h);
 338     for (k = 0; k < NROT; k++)
 339     {
 340         sfree(rot_occ[k]);
 341     }
 342
 343 }
 344
 345 void mk_multiplicity_lookup (int *multiplicity, int maxchi,
 346                              int nlist, t_dlist dlist[], int nangles)
 347 {
 348     /* new by grs - for dihedral j (as in dih[j]) get multiplicity from dlist
 349      * and store in multiplicity[j]
 350      */
 351
 352     int  j, Dih, i;
 353     char name[4];
 354
 355     j = 0;
 356     for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
 357     {
 358         for (i = 0; (i < nlist); i++)
 359         {
 360             strncpy(name, dlist[i].name, 3);
 361             name[3] = '\0';
 362             if (((Dih  < edOmega) ) ||
 363                 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
 364                 ((Dih  > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)))
 365             {
 366                 /* default - we will correct the rest below */
 367                 multiplicity[j] = 3;
 368
 369                 /* make omegas 2fold, though doesn't make much more sense than 3 */
 370                 if (Dih == edOmega && (has_dihedral(edOmega, &(dlist[i]))))
 371                 {
 372                     multiplicity[j] = 2;
 373                 }
 374
 375                 /* dihedrals to aromatic rings, COO, CONH2 or guanidinium are 2fold*/
 376                 if (Dih > edOmega && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1))
 377                 {
 378                     if ( ((strstr(name, "PHE") != NULL) && (Dih == edChi2))  ||
 379                          ((strstr(name, "TYR") != NULL) && (Dih == edChi2))  ||
 380                          ((strstr(name, "PTR") != NULL) && (Dih == edChi2))  ||
 381                          ((strstr(name, "TRP") != NULL) && (Dih == edChi2))  ||
 382                          ((strstr(name, "HIS") != NULL) && (Dih == edChi2))  ||
 383                          ((strstr(name, "GLU") != NULL) && (Dih == edChi3))  ||
 384                          ((strstr(name, "ASP") != NULL) && (Dih == edChi2))  ||
 385                          ((strstr(name, "GLN") != NULL) && (Dih == edChi3))  ||
 386                          ((strstr(name, "ASN") != NULL) && (Dih == edChi2))  ||
 387                          ((strstr(name, "ARG") != NULL) && (Dih == edChi4))  )
 388                     {
 389                         multiplicity[j] = 2;
 390                     }
 391                 }
 392                 j++;
 393             }
 394         }
 395     }
 396     if (j < nangles)
 397     {
 398         fprintf(stderr, "WARNING: not all dihedrals found in topology (only %d out of %d)!\n",
 399                 j, nangles);
 400     }
 401     /* Check for remaining dihedrals */
 402     for (; (j < nangles); j++)
 403     {
 404         multiplicity[j] = 3;
 405     }
 406
 407 }
 408
 409 void mk_chi_lookup (int **lookup, int maxchi,
 410                     int nlist, t_dlist dlist[])
 411 {
 412
 413     /* by grs. should rewrite everything to use this. (but haven't,
 414      * and at mmt only used in get_chi_product_traj
 415      * returns the dihed number given the residue number (from-0)
 416      * and chi (from-0) nr. -1 for chi undefined for that res (eg gly, ala..)*/
 417
 418     int i, j, Dih, Chi;
 419
 420     j = 0;
 421     /* NONCHI points to chi1, therefore we have to start counting there. */
 422     for (Dih = NONCHI; (Dih < NONCHI+maxchi); Dih++)
 423     {
 424         for (i = 0; (i < nlist); i++)
 425         {
 426             Chi = Dih - NONCHI;
 427             if (((Dih  < edOmega) ) ||
 428                 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
 429                 ((Dih  > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)))
 430             {
 431                 /* grs debug  printf("Not OK? i %d j %d Dih %d \n", i, j, Dih) ; */
 432                 if (Dih > edOmega)
 433                 {
 434                     lookup[i][Chi] = j;
 435                 }
 436                 j++;
 437             }
 438             else
 439             {
 440                 lookup[i][Chi] = -1;
 441             }
 442         }
 443     }
 444
 445 }
 446
 447
 448 void get_chi_product_traj (real **dih, int nframes, int nlist,
 449                            int maxchi, t_dlist dlist[], real time[],
 450                            int **lookup, int *multiplicity, gmx_bool bRb, gmx_bool bNormalize,
 451                            real core_frac, gmx_bool bAll, const char *fnall,
 452                            const output_env_t oenv)
 453 {
 454
 455     gmx_bool bRotZero, bHaveChi = FALSE;
 456     int      accum = 0, index, i, j, k, Xi, n, b;
 457     real    *chi_prtrj;
 458     int     *chi_prhist;
 459     int      nbin;
 460     FILE    *fp, *fpall;
 461     char     hisfile[256], histitle[256], *namept;
 462
 463     int      (*calc_bin)(real, int, real);
 464
 465     /* Analysis of dihedral transitions */
 466     fprintf(stderr, "Now calculating Chi product trajectories...\n");
 467
 468     if (bRb)
 469     {
 470         calc_bin = calc_RBbin;
 471     }
 472     else
 473     {
 474         calc_bin = calc_Nbin;
 475     }
 476
 477     snew(chi_prtrj, nframes);
 478
 479     /* file for info on all residues */
 480     if (bNormalize)
 481     {
 482         fpall = xvgropen(fnall, "Cumulative Rotamers", "Residue", "Probability", oenv);
 483     }
 484     else
 485     {
 486         fpall = xvgropen(fnall, "Cumulative Rotamers", "Residue", "# Counts", oenv);
 487     }
 488
 489     for (i = 0; (i < nlist); i++)
 490     {
 491
 492         /* get nbin, the nr. of cumulative rotamers that need to be considered */
 493         nbin = 1;
 494         for (Xi = 0; Xi < maxchi; Xi++)
 495         {
 496             index = lookup[i][Xi]; /* chi_(Xi+1) of res i (-1 if off end) */
 497             if (index >= 0)
 498             {
 499                 n    = multiplicity[index];
 500                 nbin = n*nbin;
 501             }
 502         }
 503         nbin += 1; /* for the "zero rotamer", outside the core region */
 504
 505         for (j = 0; (j < nframes); j++)
 506         {
 507
 508             bRotZero = FALSE;
 509             bHaveChi = TRUE;
 510             index    = lookup[i][0]; /* index into dih of chi1 of res i */
 511             if (index == -1)
 512             {
 513                 b        = 0;
 514                 bRotZero = TRUE;
 515                 bHaveChi = FALSE;
 516             }
 517             else
 518             {
 519                 b     = calc_bin(dih[index][j], multiplicity[index], core_frac);
 520                 accum = b - 1;
 521                 if (b == 0)
 522                 {
 523                     bRotZero = TRUE;
 524                 }
 525                 for (Xi = 1; Xi < maxchi; Xi++)
 526                 {
 527                     index = lookup[i][Xi]; /* chi_(Xi+1) of res i (-1 if off end) */
 528                     if (index >= 0)
 529                     {
 530                         n     = multiplicity[index];
 531                         b     = calc_bin(dih[index][j], n, core_frac);
 532                         accum = n * accum + b - 1;
 533                         if (b == 0)
 534                         {
 535                             bRotZero = TRUE;
 536                         }
 537                     }
 538                 }
 539                 accum++;
 540             }
 541             if (bRotZero)
 542             {
 543                 chi_prtrj[j] = 0.0;
 544             }
 545             else
 546             {
 547                 chi_prtrj[j] = accum;
 548                 if (accum+1 > nbin)
 549                 {
 550                     nbin = accum+1;
 551                 }
 552             }
 553         }
 554         if (bHaveChi)
 555         {
 556
 557             if (bAll)
 558             {
 559                 /* print cuml rotamer vs time */
 560                 print_one(oenv, "chiproduct", dlist[i].name, "chi product for",
 561                           "cumulative rotamer", nframes, time, chi_prtrj);
 562             }
 563
 564             /* make a histogram pf culm. rotamer occupancy too */
 565             snew(chi_prhist, nbin);
 566             make_histo(NULL, nframes, chi_prtrj, nbin, chi_prhist, 0, nbin);
 567             if (bAll)
 568             {
 569                 sprintf(hisfile, "histo-chiprod%s.xvg", dlist[i].name);
 570                 sprintf(histitle, "cumulative rotamer distribution for %s", dlist[i].name);
 571                 fprintf(stderr, "  and %s  ", hisfile);
 572                 fp = xvgropen(hisfile, histitle, "number", "", oenv);
 573                 if (output_env_get_print_xvgr_codes(oenv))
 574                 {
 575                     fprintf(fp, "@ xaxis tick on\n");
 576                     fprintf(fp, "@ xaxis tick major 1\n");
 577                     fprintf(fp, "@ type xy\n");
 578                 }
 579                 for (k = 0; (k < nbin); k++)
 580                 {
 581                     if (bNormalize)
 582                     {
 583                         fprintf(fp, "%5d  %10g\n", k, (1.0*chi_prhist[k])/nframes);
 584                     }
 585                     else
 586                     {
 587                         fprintf(fp, "%5d  %10d\n", k, chi_prhist[k]);
 588                     }
 589                 }
 590                 fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 591                 gmx_ffclose(fp);
 592             }
 593
 594             /* and finally print out occupancies to a single file */
 595             /* get the gmx from-1 res nr by setting a ptr to the number part
 596              * of dlist[i].name - potential bug for 4-letter res names... */
 597             namept = dlist[i].name + 3;
 598             fprintf(fpall, "%5s ", namept);
 599             for (k = 0; (k < nbin); k++)
 600             {
 601                 if (bNormalize)
 602                 {
 603                     fprintf(fpall, "  %10g", (1.0*chi_prhist[k])/nframes);
 604                 }
 605                 else
 606                 {
 607                     fprintf(fpall, "  %10d", chi_prhist[k]);
 608                 }
 609             }
 610             fprintf(fpall, "\n");
 611
 612             sfree(chi_prhist);
 613             /* histogram done */
 614         }
 615     }
 616
 617     sfree(chi_prtrj);
 618     gmx_ffclose(fpall);
 619     fprintf(stderr, "\n");
 620
 621 }
 622
 623 void calc_distribution_props(int nh, int histo[], real start,
 624                              int nkkk, t_karplus kkk[],
 625                              real *S2)
 626 {
 627     real d, dc, ds, c1, c2, tdc, tds;
 628     real fac, ang, invth, Jc;
 629     int  i, j, th;
 630
 631     if (nh == 0)
 632     {
 633         gmx_fatal(FARGS, "No points in histogram (%s, %d)", __FILE__, __LINE__);
 634     }
 635     fac = 2*M_PI/nh;
 636
 637     /* Compute normalisation factor */
 638     th = 0;
 639     for (j = 0; (j < nh); j++)
 640     {
 641         th += histo[j];
 642     }
 643     invth = 1.0/th;
 644
 645     for (i = 0; (i < nkkk); i++)
 646     {
 647         kkk[i].Jc    = 0;
 648         kkk[i].Jcsig = 0;
 649     }
 650     tdc = 0, tds = 0;
 651     for (j = 0; (j < nh); j++)
 652     {
 653         d    = invth*histo[j];
 654         ang  = j*fac-start;
 655         c1   = cos(ang);
 656         c2   = c1*c1;
 657         dc   = d*c1;
 658         ds   = d*sin(ang);
 659         tdc += dc;
 660         tds += ds;
 661         for (i = 0; (i < nkkk); i++)
 662         {
 663             c1            = 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  = 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;
 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 = min(minx, data[i]);
 763             maxx = 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 = (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, natoms, 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     natoms = 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] = 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 = (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 }