src/gromacs/gmxana/gmx_chi.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38 #include <stdio.h>
  39 #include <math.h>
  40
  41 #include "gromacs/fileio/confio.h"
  42 #include "gromacs/fileio/pdbio.h"
  43 #include "gmx_fatal.h"
  44 #include "gromacs/fileio/futil.h"
  45 #include "gstat.h"
  46 #include "macros.h"
  47 #include "maths.h"
  48 #include "physics.h"
  49 #include "index.h"
  50 #include "smalloc.h"
  51 #include "statutil.h"
  52 #include "gromacs/fileio/tpxio.h"
  53 #include <string.h>
  54 #include "sysstuff.h"
  55 #include "txtdump.h"
  56 #include "typedefs.h"
  57 #include "vec.h"
  58 #include "strdb.h"
  59 #include "xvgr.h"
  60 #include "gromacs/fileio/matio.h"
  61 #include "gmx_ana.h"
  62
  63 static gmx_bool bAllowed(real phi, real psi)
  64 {
  65     static const char *map[] = {
  66         "1100000000000000001111111000000000001111111111111111111111111",
  67         "1100000000000000001111110000000000011111111111111111111111111",
  68         "1100000000000000001111110000000000011111111111111111111111111",
  69         "1100000000000000001111100000000000111111111111111111111111111",
  70         "1100000000000000001111100000000000111111111111111111111111111",
  71         "1100000000000000001111100000000001111111111111111111111111111",
  72         "1100000000000000001111100000000001111111111111111111111111111",
  73         "1100000000000000001111100000000011111111111111111111111111111",
  74         "1110000000000000001111110000000111111111111111111111111111111",
  75         "1110000000000000001111110000001111111111111111111111111111111",
  76         "1110000000000000001111111000011111111111111111111111111111111",
  77         "1110000000000000001111111100111111111111111111111111111111111",
  78         "1110000000000000001111111111111111111111111111111111111111111",
  79         "1110000000000000001111111111111111111111111111111111111111111",
  80         "1110000000000000001111111111111111111111111111111111111111111",
  81         "1110000000000000001111111111111111111111111111111111111111111",
  82         "1110000000000000001111111111111110011111111111111111111111111",
  83         "1110000000000000001111111111111100000111111111111111111111111",
  84         "1110000000000000001111111111111000000000001111111111111111111",
  85         "1100000000000000001111111111110000000000000011111111111111111",
  86         "1100000000000000001111111111100000000000000011111111111111111",
  87         "1000000000000000001111111111000000000000000001111111111111110",
  88         "0000000000000000001111111110000000000000000000111111111111100",
  89         "0000000000000000000000000000000000000000000000000000000000000",
  90         "0000000000000000000000000000000000000000000000000000000000000",
  91         "0000000000000000000000000000000000000000000000000000000000000",
  92         "0000000000000000000000000000000000000000000000000000000000000",
  93         "0000000000000000000000000000000000000000000000000000000000000",
  94         "0000000000000000000000000000000000000000000000000000000000000",
  95         "0000000000000000000000000000000000000000000000000000000000000",
  96         "0000000000000000000000000000000000000000000000000000000000000",
  97         "0000000000000000000000000000000000000000000000000000000000000",
  98         "0000000000000000000000000000000000000000000000000000000000000",
  99         "0000000000000000000000000000000000000000000000000000000000000",
 100         "0000000000000000000000000000000000000000000000000000000000000",
 101         "0000000000000000000000000000000000000000000000000000000000000",
 102         "0000000000000000000000000000000000000000000000000000000000000",
 103         "0000000000000000000000000000000000000000000000000000000000000",
 104         "0000000000000000000000000000000000111111111111000000000000000",
 105         "1100000000000000000000000000000001111111111111100000000000111",
 106         "1100000000000000000000000000000001111111111111110000000000111",
 107         "0000000000000000000000000000000000000000000000000000000000000",
 108         "0000000000000000000000000000000000000000000000000000000000000",
 109         "0000000000000000000000000000000000000000000000000000000000000",
 110         "0000000000000000000000000000000000000000000000000000000000000",
 111         "0000000000000000000000000000000000000000000000000000000000000",
 112         "0000000000000000000000000000000000000000000000000000000000000",
 113         "0000000000000000000000000000000000000000000000000000000000000",
 114         "0000000000000000000000000000000000000000000000000000000000000",
 115         "0000000000000000000000000000000000000000000000000000000000000",
 116         "0000000000000000000000000000000000000000000000000000000000000",
 117         "0000000000000000000000000000000000000000000000000000000000000",
 118         "0000000000000000000000000000000000000000000000000000000000000",
 119         "0000000000000000000000000000000000000000000000000000000000000",
 120         "0000000000000000000000000000000000000000000000000000000000000",
 121         "0000000000000000000000000000000000000000000000000000000000000",
 122         "0000000000000000000000000000000000000000000000000000000000000",
 123         "0000000000000000000000000000000000000000000000000000000000000",
 124         "0000000000000000000000000000000000000000000000000000000000000",
 125         "0000000000000000000000000000000000000000000000000000000000000",
 126         "0000000000000000000000000000000000000000000000000000000000000"
 127     };
 128 #define NPP asize(map)
 129     int                x, y;
 130
 131 #define INDEX(ppp) ((((int) (360+ppp*RAD2DEG)) % 360)/6)
 132     x = INDEX(phi);
 133     y = INDEX(psi);
 134 #undef INDEX
 135     return (gmx_bool) map[x][y];
 136 }
 137
 138 atom_id *make_chi_ind(int nl, t_dlist dl[], int *ndih)
 139 {
 140     atom_id *id;
 141     int      i, Xi, n;
 142
 143     /* There are nl residues with max edMax dihedrals with 4 atoms each */
 144     snew(id, nl*edMax*4);
 145
 146     n = 0;
 147     for (i = 0; (i < nl); i++)
 148     {
 149         /* Phi, fake the first one */
 150         dl[i].j0[edPhi] = n/4;
 151         if (dl[i].atm.minC >= 0)
 152         {
 153             id[n++] = dl[i].atm.minC;
 154         }
 155         else
 156         {
 157             id[n++] = dl[i].atm.H;
 158         }
 159         id[n++] = dl[i].atm.N;
 160         id[n++] = dl[i].atm.Cn[1];
 161         id[n++] = dl[i].atm.C;
 162     }
 163     for (i = 0; (i < nl); i++)
 164     {
 165         /* Psi, fake the last one */
 166         dl[i].j0[edPsi] = n/4;
 167         id[n++]         = dl[i].atm.N;
 168         id[n++]         = dl[i].atm.Cn[1];
 169         id[n++]         = dl[i].atm.C;
 170         if (i < (nl-1) )
 171         {
 172             id[n++] = dl[i+1].atm.N;
 173         }
 174         else
 175         {
 176             id[n++] = dl[i].atm.O;
 177         }
 178     }
 179     for (i = 1; (i < nl); i++)
 180     {
 181         /* Omega */
 182         if (has_dihedral(edOmega, &(dl[i])))
 183         {
 184             dl[i].j0[edOmega] = n/4;
 185             id[n++]           = dl[i].atm.minCalpha;
 186             id[n++]           = dl[i].atm.minC;
 187             id[n++]           = dl[i].atm.N;
 188             id[n++]           = dl[i].atm.Cn[1];
 189         }
 190     }
 191     for (Xi = 0; (Xi < MAXCHI); Xi++)
 192     {
 193         /* Chi# */
 194         for (i = 0; (i < nl); i++)
 195         {
 196             if (dl[i].atm.Cn[Xi+3] != -1)
 197             {
 198                 dl[i].j0[edChi1+Xi] = n/4;
 199                 id[n++]             = dl[i].atm.Cn[Xi];
 200                 id[n++]             = dl[i].atm.Cn[Xi+1];
 201                 id[n++]             = dl[i].atm.Cn[Xi+2];
 202                 id[n++]             = dl[i].atm.Cn[Xi+3];
 203             }
 204         }
 205     }
 206     *ndih = n/4;
 207
 208     return id;
 209 }
 210
 211 int bin(real chi, int mult)
 212 {
 213     mult = 3;
 214
 215     return (int) (chi*mult/360.0);
 216 }
 217
 218
 219 static void do_dihcorr(const char *fn, int nf, int ndih, real **dih, real dt,
 220                        int nlist, t_dlist dlist[], real time[], int maxchi,
 221                        gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bOmega,
 222                        const output_env_t oenv)
 223 {
 224     char name1[256], name2[256];
 225     int  i, j, Xi;
 226
 227     do_autocorr(fn, oenv, "Dihedral Autocorrelation Function",
 228                 nf, ndih, dih, dt, eacCos, FALSE);
 229     /* Dump em all */
 230     j = 0;
 231     for (i = 0; (i < nlist); i++)
 232     {
 233         if (bPhi)
 234         {
 235             print_one(oenv, "corrphi", dlist[i].name, "Phi ACF for", "C(t)", nf/2, time,
 236                       dih[j]);
 237         }
 238         j++;
 239     }
 240     for (i = 0; (i < nlist); i++)
 241     {
 242         if (bPsi)
 243         {
 244             print_one(oenv, "corrpsi", dlist[i].name, "Psi ACF for", "C(t)", nf/2, time,
 245                       dih[j]);
 246         }
 247         j++;
 248     }
 249     for (i = 0; (i < nlist); i++)
 250     {
 251         if (has_dihedral(edOmega, &dlist[i]))
 252         {
 253             if (bOmega)
 254             {
 255                 print_one(oenv, "corromega", dlist[i].name, "Omega ACF for", "C(t)",
 256                           nf/2, time, dih[j]);
 257             }
 258             j++;
 259         }
 260     }
 261     for (Xi = 0; (Xi < maxchi); Xi++)
 262     {
 263         sprintf(name1, "corrchi%d", Xi+1);
 264         sprintf(name2, "Chi%d ACF for", Xi+1);
 265         for (i = 0; (i < nlist); i++)
 266         {
 267             if (dlist[i].atm.Cn[Xi+3] != -1)
 268             {
 269                 if (bChi)
 270                 {
 271                     print_one(oenv, name1, dlist[i].name, name2, "C(t)", nf/2, time, dih[j]);
 272                 }
 273                 j++;
 274             }
 275         }
 276     }
 277     fprintf(stderr, "\n");
 278 }
 279
 280 static void copy_dih_data(real in[], real out[], int nf, gmx_bool bLEAVE)
 281 {
 282     /* if bLEAVE, do nothing to data in copying to out
 283      * otherwise multiply by 180/pi to convert rad to deg */
 284     int  i;
 285     real mult;
 286     if (bLEAVE)
 287     {
 288         mult = 1;
 289     }
 290     else
 291     {
 292         mult = (180.0/M_PI);
 293     }
 294     for (i = 0; (i < nf); i++)
 295     {
 296         out[i] = in[i]*mult;
 297     }
 298 }
 299
 300 static void dump_em_all(int nlist, t_dlist dlist[], int nf, real time[],
 301                         real **dih, int maxchi,
 302                         gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bOmega, gmx_bool bRAD,
 303                         const output_env_t oenv)
 304 {
 305     char  name[256], titlestr[256], ystr[256];
 306     real *data;
 307     int   i, j, Xi;
 308
 309     snew(data, nf);
 310     if (bRAD)
 311     {
 312         strcpy(ystr, "Angle (rad)");
 313     }
 314     else
 315     {
 316         strcpy(ystr, "Angle (degrees)");
 317     }
 318
 319     /* Dump em all */
 320     j = 0;
 321     for (i = 0; (i < nlist); i++)
 322     {
 323         /* grs debug  printf("OK i %d j %d\n", i, j) ; */
 324         if (bPhi)
 325         {
 326             copy_dih_data(dih[j], data, nf, bRAD);
 327             print_one(oenv, "phi", dlist[i].name, "\\xf\\f{}", ystr, nf, time, data);
 328         }
 329         j++;
 330     }
 331     for (i = 0; (i < nlist); i++)
 332     {
 333         if (bPsi)
 334         {
 335             copy_dih_data(dih[j], data, nf, bRAD);
 336             print_one(oenv, "psi", dlist[i].name, "\\xy\\f{}", ystr, nf, time, data);
 337         }
 338         j++;
 339     }
 340     for (i = 0; (i < nlist); i++)
 341     {
 342         if (has_dihedral(edOmega, &(dlist[i])))
 343         {
 344             if (bOmega)
 345             {
 346                 copy_dih_data(dih[j], data, nf, bRAD);
 347                 print_one(oenv, "omega", dlist[i].name, "\\xw\\f{}", ystr, nf, time, data);
 348             }
 349             j++;
 350         }
 351     }
 352
 353     for (Xi = 0; (Xi < maxchi); Xi++)
 354     {
 355         for (i = 0; (i < nlist); i++)
 356         {
 357             if (dlist[i].atm.Cn[Xi+3] != -1)
 358             {
 359                 if (bChi)
 360                 {
 361                     sprintf(name, "chi%d", Xi+1);
 362                     sprintf(titlestr, "\\xc\\f{}\\s%d\\N", Xi+1);
 363                     copy_dih_data(dih[j], data, nf, bRAD);
 364                     print_one(oenv, name, dlist[i].name, titlestr, ystr, nf, time, data);
 365                 }
 366                 j++;
 367             }
 368         }
 369     }
 370     fprintf(stderr, "\n");
 371 }
 372
 373 static void reset_one(real dih[], int nf, real phase)
 374 {
 375     int j;
 376
 377     for (j = 0; (j < nf); j++)
 378     {
 379         dih[j] += phase;
 380         while (dih[j] < -M_PI)
 381         {
 382             dih[j] += 2*M_PI;
 383         }
 384         while (dih[j] >= M_PI)
 385         {
 386             dih[j] -= 2*M_PI;
 387         }
 388     }
 389 }
 390
 391 static int reset_em_all(int nlist, t_dlist dlist[], int nf,
 392                         real **dih, int maxchi)
 393 {
 394     int  i, j, Xi;
 395
 396     /* Reset em all */
 397     j = 0;
 398     /* Phi */
 399     for (i = 0; (i < nlist); i++)
 400     {
 401         if (dlist[i].atm.minC == -1)
 402         {
 403             reset_one(dih[j++], nf, M_PI);
 404         }
 405         else
 406         {
 407             reset_one(dih[j++], nf, 0);
 408         }
 409     }
 410     /* Psi */
 411     for (i = 0; (i < nlist-1); i++)
 412     {
 413         reset_one(dih[j++], nf, 0);
 414     }
 415     /* last Psi is faked from O */
 416     reset_one(dih[j++], nf, M_PI);
 417
 418     /* Omega */
 419     for (i = 0; (i < nlist); i++)
 420     {
 421         if (has_dihedral(edOmega, &dlist[i]))
 422         {
 423             reset_one(dih[j++], nf, 0);
 424         }
 425     }
 426     /* Chi 1 thru maxchi */
 427     for (Xi = 0; (Xi < maxchi); Xi++)
 428     {
 429         for (i = 0; (i < nlist); i++)
 430         {
 431             if (dlist[i].atm.Cn[Xi+3] != -1)
 432             {
 433                 reset_one(dih[j], nf, 0);
 434                 j++;
 435             }
 436         }
 437     }
 438     fprintf(stderr, "j after resetting (nr. active dihedrals) = %d\n", j);
 439     return j;
 440 }
 441
 442 static void histogramming(FILE *log, int nbin, gmx_residuetype_t rt,
 443                           int nf, int maxchi, real **dih,
 444                           int nlist, t_dlist dlist[],
 445                           atom_id index[],
 446                           gmx_bool bPhi, gmx_bool bPsi, gmx_bool bOmega, gmx_bool bChi,
 447                           gmx_bool bNormalize, gmx_bool bSSHisto, const char *ssdump,
 448                           real bfac_max, t_atoms *atoms,
 449                           gmx_bool bDo_jc, const char *fn,
 450                           const output_env_t oenv)
 451 {
 452     /* also gets 3J couplings and order parameters S2 */
 453     t_karplus kkkphi[] = {
 454         { "J_NHa1",    6.51, -1.76,  1.6, -M_PI/3,   0.0,  0.0 },
 455         { "J_NHa2",    6.51, -1.76,  1.6,  M_PI/3,   0.0,  0.0 },
 456         { "J_HaC'",    4.0,   1.1,   0.1,  0.0,      0.0,  0.0 },
 457         { "J_NHCb",    4.7,  -1.5,  -0.2,  M_PI/3,   0.0,  0.0 },
 458         { "J_Ci-1Hai", 4.5,  -1.3,  -1.2,  2*M_PI/3, 0.0,  0.0 }
 459     };
 460     t_karplus kkkpsi[] = {
 461         { "J_HaN",   -0.88, -0.61, -0.27, M_PI/3,  0.0,  0.0 }
 462     };
 463     t_karplus kkkchi1[] = {
 464         { "JHaHb2",       9.5, -1.6, 1.8, -M_PI/3, 0,  0.0 },
 465         { "JHaHb3",       9.5, -1.6, 1.8, 0, 0,  0.0 }
 466     };
 467 #define NKKKPHI asize(kkkphi)
 468 #define NKKKPSI asize(kkkpsi)
 469 #define NKKKCHI asize(kkkchi1)
 470 #define NJC (NKKKPHI+NKKKPSI+NKKKCHI)
 471
 472     FILE       *fp, *ssfp[3] = {NULL, NULL, NULL};
 473     const char *sss[3] = { "sheet", "helix", "coil" };
 474     real        S2;
 475     real       *normhisto;
 476     real      **Jc, **Jcsig;
 477     int     ****his_aa_ss = NULL;
 478     int      ***his_aa, **his_aa1, *histmp;
 479     int         i, j, k, m, n, nn, Dih, nres, hindex, angle;
 480     gmx_bool    bBfac, bOccup;
 481     char        hisfile[256], hhisfile[256], sshisfile[256], title[256], *ss_str = NULL;
 482     char      **leg;
 483     const char *residue_name;
 484     int         rt_size;
 485
 486     rt_size = gmx_residuetype_get_size(rt);
 487     if (bSSHisto)
 488     {
 489         fp = ffopen(ssdump, "r");
 490         if (1 != fscanf(fp, "%d", &nres))
 491         {
 492             gmx_fatal(FARGS, "Error reading from file %s", ssdump);
 493         }
 494
 495         snew(ss_str, nres+1);
 496         if (1 != fscanf(fp, "%s", ss_str))
 497         {
 498             gmx_fatal(FARGS, "Error reading from file %s", ssdump);
 499         }
 500
 501         ffclose(fp);
 502         /* Four dimensional array... Very cool */
 503         snew(his_aa_ss, 3);
 504         for (i = 0; (i < 3); i++)
 505         {
 506             snew(his_aa_ss[i], rt_size+1);
 507             for (j = 0; (j <= rt_size); j++)
 508             {
 509                 snew(his_aa_ss[i][j], edMax);
 510                 for (Dih = 0; (Dih < edMax); Dih++)
 511                 {
 512                     snew(his_aa_ss[i][j][Dih], nbin+1);
 513                 }
 514             }
 515         }
 516     }
 517     snew(his_aa, edMax);
 518     for (Dih = 0; (Dih < edMax); Dih++)
 519     {
 520         snew(his_aa[Dih], rt_size+1);
 521         for (i = 0; (i <= rt_size); i++)
 522         {
 523             snew(his_aa[Dih][i], nbin+1);
 524         }
 525     }
 526     snew(histmp, nbin);
 527
 528     snew(Jc, nlist);
 529     snew(Jcsig, nlist);
 530     for (i = 0; (i < nlist); i++)
 531     {
 532         snew(Jc[i], NJC);
 533         snew(Jcsig[i], NJC);
 534     }
 535
 536     j = 0;
 537     n = 0;
 538     for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
 539     {
 540         for (i = 0; (i < nlist); i++)
 541         {
 542             if (((Dih  < edOmega) ) ||
 543                 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
 544                 ((Dih  > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)))
 545             {
 546                 make_histo(log, nf, dih[j], nbin, histmp, -M_PI, M_PI);
 547
 548                 if (bSSHisto)
 549                 {
 550                     /* Assume there is only one structure, the first.
 551                      * Compute index in histogram.
 552                      */
 553                     /* Check the atoms to see whether their B-factors are low enough
 554                      * Check atoms to see their occupancy is 1.
 555                      */
 556                     bBfac = bOccup = TRUE;
 557                     for (nn = 0; (nn < 4); nn++, n++)
 558                     {
 559                         bBfac  = bBfac  && (atoms->pdbinfo[index[n]].bfac <= bfac_max);
 560                         bOccup = bOccup && (atoms->pdbinfo[index[n]].occup == 1);
 561                     }
 562                     if (bOccup && ((bfac_max <= 0) || ((bfac_max > 0) && bBfac)))
 563                     {
 564                         hindex = ((dih[j][0]+M_PI)*nbin)/(2*M_PI);
 565                         range_check(hindex, 0, nbin);
 566
 567                         /* Assign dihedral to either of the structure determined
 568                          * histograms
 569                          */
 570                         switch (ss_str[dlist[i].resnr])
 571                         {
 572                             case 'E':
 573                                 his_aa_ss[0][dlist[i].index][Dih][hindex]++;
 574                                 break;
 575                             case 'H':
 576                                 his_aa_ss[1][dlist[i].index][Dih][hindex]++;
 577                                 break;
 578                             default:
 579                                 his_aa_ss[2][dlist[i].index][Dih][hindex]++;
 580                                 break;
 581                         }
 582                     }
 583                     else if (debug)
 584                     {
 585                         fprintf(debug, "Res. %d has imcomplete occupancy or bfacs > %g\n",
 586                                 dlist[i].resnr, bfac_max);
 587                     }
 588                 }
 589                 else
 590                 {
 591                     n += 4;
 592                 }
 593
 594                 switch (Dih)
 595                 {
 596                     case edPhi:
 597                         calc_distribution_props(nbin, histmp, -M_PI, NKKKPHI, kkkphi, &S2);
 598
 599                         for (m = 0; (m < NKKKPHI); m++)
 600                         {
 601                             Jc[i][m]    = kkkphi[m].Jc;
 602                             Jcsig[i][m] = kkkphi[m].Jcsig;
 603                         }
 604                         break;
 605                     case edPsi:
 606                         calc_distribution_props(nbin, histmp, -M_PI, NKKKPSI, kkkpsi, &S2);
 607
 608                         for (m = 0; (m < NKKKPSI); m++)
 609                         {
 610                             Jc[i][NKKKPHI+m]    = kkkpsi[m].Jc;
 611                             Jcsig[i][NKKKPHI+m] = kkkpsi[m].Jcsig;
 612                         }
 613                         break;
 614                     case edChi1:
 615                         calc_distribution_props(nbin, histmp, -M_PI, NKKKCHI, kkkchi1, &S2);
 616                         for (m = 0; (m < NKKKCHI); m++)
 617                         {
 618                             Jc[i][NKKKPHI+NKKKPSI+m]    = kkkchi1[m].Jc;
 619                             Jcsig[i][NKKKPHI+NKKKPSI+m] = kkkchi1[m].Jcsig;
 620                         }
 621                         break;
 622                     default: /* covers edOmega and higher Chis than Chi1 */
 623                         calc_distribution_props(nbin, histmp, -M_PI, 0, NULL, &S2);
 624                         break;
 625                 }
 626                 dlist[i].S2[Dih]        = S2;
 627
 628                 /* Sum distribution per amino acid type as well */
 629                 for (k = 0; (k < nbin); k++)
 630                 {
 631                     his_aa[Dih][dlist[i].index][k] += histmp[k];
 632                     histmp[k] = 0;
 633                 }
 634                 j++;
 635             }
 636             else /* dihed not defined */
 637             {
 638                 dlist[i].S2[Dih] = 0.0;
 639             }
 640         }
 641     }
 642     sfree(histmp);
 643
 644     /* Print out Jcouplings */
 645     fprintf(log, "\n *** J-Couplings from simulation (plus std. dev.) ***\n\n");
 646     fprintf(log, "Residue   ");
 647     for (i = 0; (i < NKKKPHI); i++)
 648     {
 649         fprintf(log, "%7s   SD", kkkphi[i].name);
 650     }
 651     for (i = 0; (i < NKKKPSI); i++)
 652     {
 653         fprintf(log, "%7s   SD", kkkpsi[i].name);
 654     }
 655     for (i = 0; (i < NKKKCHI); i++)
 656     {
 657         fprintf(log, "%7s   SD", kkkchi1[i].name);
 658     }
 659     fprintf(log, "\n");
 660     for (i = 0; (i < NJC+1); i++)
 661     {
 662         fprintf(log, "------------");
 663     }
 664     fprintf(log, "\n");
 665     for (i = 0; (i < nlist); i++)
 666     {
 667         fprintf(log, "%-10s", dlist[i].name);
 668         for (j = 0; (j < NJC); j++)
 669         {
 670             fprintf(log, "  %5.2f %4.2f", Jc[i][j], Jcsig[i][j]);
 671         }
 672         fprintf(log, "\n");
 673     }
 674     fprintf(log, "\n");
 675
 676     /* and to -jc file... */
 677     if (bDo_jc)
 678     {
 679         fp = xvgropen(fn, "\\S3\\NJ-Couplings from Karplus Equation", "Residue",
 680                       "Coupling", oenv);
 681         snew(leg, NJC);
 682         for (i = 0; (i < NKKKPHI); i++)
 683         {
 684             leg[i] = strdup(kkkphi[i].name);
 685         }
 686         for (i = 0; (i < NKKKPSI); i++)
 687         {
 688             leg[i+NKKKPHI] = strdup(kkkpsi[i].name);
 689         }
 690         for (i = 0; (i < NKKKCHI); i++)
 691         {
 692             leg[i+NKKKPHI+NKKKPSI] = strdup(kkkchi1[i].name);
 693         }
 694         xvgr_legend(fp, NJC, (const char**)leg, oenv);
 695         fprintf(fp, "%5s ", "#Res.");
 696         for (i = 0; (i < NJC); i++)
 697         {
 698             fprintf(fp, "%10s ", leg[i]);
 699         }
 700         fprintf(fp, "\n");
 701         for (i = 0; (i < nlist); i++)
 702         {
 703             fprintf(fp, "%5d ", dlist[i].resnr);
 704             for (j = 0; (j < NJC); j++)
 705             {
 706                 fprintf(fp, "  %8.3f", Jc[i][j]);
 707             }
 708             fprintf(fp, "\n");
 709         }
 710         ffclose(fp);
 711         for (i = 0; (i < NJC); i++)
 712         {
 713             sfree(leg[i]);
 714         }
 715     }
 716     /* finished -jc stuff */
 717
 718     snew(normhisto, nbin);
 719     for (i = 0; (i < rt_size); i++)
 720     {
 721         for (Dih = 0; (Dih < edMax); Dih++)
 722         {
 723             /* First check whether something is in there */
 724             for (j = 0; (j < nbin); j++)
 725             {
 726                 if (his_aa[Dih][i][j] != 0)
 727                 {
 728                     break;
 729                 }
 730             }
 731             if ((j < nbin) &&
 732                 ((bPhi && (Dih == edPhi)) ||
 733                  (bPsi && (Dih == edPsi)) ||
 734                  (bOmega && (Dih == edOmega)) ||
 735                  (bChi && (Dih >= edChi1))))
 736             {
 737                 if (bNormalize)
 738                 {
 739                     normalize_histo(nbin, his_aa[Dih][i], (360.0/nbin), normhisto);
 740                 }
 741
 742                 residue_name = gmx_residuetype_get_name(rt, i);
 743                 switch (Dih)
 744                 {
 745                     case edPhi:
 746                         sprintf(hisfile, "histo-phi%s", residue_name);
 747                         sprintf(title, "\\xf\\f{} Distribution for %s", residue_name);
 748                         break;
 749                     case edPsi:
 750                         sprintf(hisfile, "histo-psi%s", residue_name);
 751                         sprintf(title, "\\xy\\f{} Distribution for %s", residue_name);
 752                         break;
 753                     case edOmega:
 754                         sprintf(hisfile, "histo-omega%s", residue_name);
 755                         sprintf(title, "\\xw\\f{} Distribution for %s", residue_name);
 756                         break;
 757                     default:
 758                         sprintf(hisfile, "histo-chi%d%s", Dih-NONCHI+1, residue_name);
 759                         sprintf(title, "\\xc\\f{}\\s%d\\N Distribution for %s",
 760                                 Dih-NONCHI+1, residue_name);
 761                 }
 762                 strcpy(hhisfile, hisfile);
 763                 strcat(hhisfile, ".xvg");
 764                 fp = xvgropen(hhisfile, title, "Degrees", "", oenv);
 765                 fprintf(fp, "@ with g0\n");
 766                 xvgr_world(fp, -180, 0, 180, 0.1, oenv);
 767                 fprintf(fp, "# this effort to set graph size fails unless you run with -autoscale none or -autoscale y flags\n");
 768                 fprintf(fp, "@ xaxis tick on\n");
 769                 fprintf(fp, "@ xaxis tick major 90\n");
 770                 fprintf(fp, "@ xaxis tick minor 30\n");
 771                 fprintf(fp, "@ xaxis ticklabel prec 0\n");
 772                 fprintf(fp, "@ yaxis tick off\n");
 773                 fprintf(fp, "@ yaxis ticklabel off\n");
 774                 fprintf(fp, "@ type xy\n");
 775                 if (bSSHisto)
 776                 {
 777                     for (k = 0; (k < 3); k++)
 778                     {
 779                         sprintf(sshisfile, "%s-%s.xvg", hisfile, sss[k]);
 780                         ssfp[k] = ffopen(sshisfile, "w");
 781                     }
 782                 }
 783                 for (j = 0; (j < nbin); j++)
 784                 {
 785                     angle = -180 + (360/nbin)*j;
 786                     if (bNormalize)
 787                     {
 788                         fprintf(fp, "%5d  %10g\n", angle, normhisto[j]);
 789                     }
 790                     else
 791                     {
 792                         fprintf(fp, "%5d  %10d\n", angle, his_aa[Dih][i][j]);
 793                     }
 794                     if (bSSHisto)
 795                     {
 796                         for (k = 0; (k < 3); k++)
 797                         {
 798                             fprintf(ssfp[k], "%5d  %10d\n", angle,
 799                                     his_aa_ss[k][i][Dih][j]);
 800                         }
 801                     }
 802                 }
 803                 fprintf(fp, "&\n");
 804                 ffclose(fp);
 805                 if (bSSHisto)
 806                 {
 807                     for (k = 0; (k < 3); k++)
 808                     {
 809                         fprintf(ssfp[k], "&\n");
 810                         ffclose(ssfp[k]);
 811                     }
 812                 }
 813             }
 814         }
 815     }
 816     sfree(normhisto);
 817
 818     if (bSSHisto)
 819     {
 820         /* Four dimensional array... Very cool */
 821         for (i = 0; (i < 3); i++)
 822         {
 823             for (j = 0; (j <= rt_size); j++)
 824             {
 825                 for (Dih = 0; (Dih < edMax); Dih++)
 826                 {
 827                     sfree(his_aa_ss[i][j][Dih]);
 828                 }
 829                 sfree(his_aa_ss[i][j]);
 830             }
 831             sfree(his_aa_ss[i]);
 832         }
 833         sfree(his_aa_ss);
 834         sfree(ss_str);
 835     }
 836 }
 837
 838 static FILE *rama_file(const char *fn, const char *title, const char *xaxis,
 839                        const char *yaxis, const output_env_t oenv)
 840 {
 841     FILE *fp;
 842
 843     fp = xvgropen(fn, title, xaxis, yaxis, oenv);
 844     fprintf(fp, "@ with g0\n");
 845     xvgr_world(fp, -180, -180, 180, 180, oenv);
 846     fprintf(fp, "@ xaxis tick on\n");
 847     fprintf(fp, "@ xaxis tick major 90\n");
 848     fprintf(fp, "@ xaxis tick minor 30\n");
 849     fprintf(fp, "@ xaxis ticklabel prec 0\n");
 850     fprintf(fp, "@ yaxis tick on\n");
 851     fprintf(fp, "@ yaxis tick major 90\n");
 852     fprintf(fp, "@ yaxis tick minor 30\n");
 853     fprintf(fp, "@ yaxis ticklabel prec 0\n");
 854     fprintf(fp, "@    s0 type xy\n");
 855     fprintf(fp, "@    s0 symbol 2\n");
 856     fprintf(fp, "@    s0 symbol size 0.410000\n");
 857     fprintf(fp, "@    s0 symbol fill 1\n");
 858     fprintf(fp, "@    s0 symbol color 1\n");
 859     fprintf(fp, "@    s0 symbol linewidth 1\n");
 860     fprintf(fp, "@    s0 symbol linestyle 1\n");
 861     fprintf(fp, "@    s0 symbol center false\n");
 862     fprintf(fp, "@    s0 symbol char 0\n");
 863     fprintf(fp, "@    s0 skip 0\n");
 864     fprintf(fp, "@    s0 linestyle 0\n");
 865     fprintf(fp, "@    s0 linewidth 1\n");
 866     fprintf(fp, "@ type xy\n");
 867
 868     return fp;
 869 }
 870
 871 static void do_rama(int nf, int nlist, t_dlist dlist[], real **dih,
 872                     gmx_bool bViol, gmx_bool bRamOmega, const output_env_t oenv)
 873 {
 874     FILE    *fp, *gp = NULL;
 875     gmx_bool bOm;
 876     char     fn[256];
 877     int      i, j, k, Xi1, Xi2, Phi, Psi, Om = 0, nlevels;
 878 #define NMAT 120
 879     real   **mat  = NULL, phi, psi, omega, axis[NMAT], lo, hi;
 880     t_rgb    rlo  = { 1.0, 0.0, 0.0 };
 881     t_rgb    rmid = { 1.0, 1.0, 1.0 };
 882     t_rgb    rhi  = { 0.0, 0.0, 1.0 };
 883
 884     for (i = 0; (i < nlist); i++)
 885     {
 886         if ((has_dihedral(edPhi, &(dlist[i]))) &&
 887             (has_dihedral(edPsi, &(dlist[i]))))
 888         {
 889             sprintf(fn, "ramaPhiPsi%s.xvg", dlist[i].name);
 890             fp = rama_file(fn, "Ramachandran Plot",
 891                            "\\8f\\4 (deg)", "\\8y\\4 (deg)", oenv);
 892             bOm = bRamOmega && has_dihedral(edOmega, &(dlist[i]));
 893             if (bOm)
 894             {
 895                 Om = dlist[i].j0[edOmega];
 896                 snew(mat, NMAT);
 897                 for (j = 0; (j < NMAT); j++)
 898                 {
 899                     snew(mat[j], NMAT);
 900                     axis[j] = -180+(360*j)/NMAT;
 901                 }
 902             }
 903             if (bViol)
 904             {
 905                 sprintf(fn, "violPhiPsi%s.xvg", dlist[i].name);
 906                 gp = ffopen(fn, "w");
 907             }
 908             Phi = dlist[i].j0[edPhi];
 909             Psi = dlist[i].j0[edPsi];
 910             for (j = 0; (j < nf); j++)
 911             {
 912                 phi = RAD2DEG*dih[Phi][j];
 913                 psi = RAD2DEG*dih[Psi][j];
 914                 fprintf(fp, "%10g  %10g\n", phi, psi);
 915                 if (bViol)
 916                 {
 917                     fprintf(gp, "%d\n", !bAllowed(dih[Phi][j], RAD2DEG*dih[Psi][j]));
 918                 }
 919                 if (bOm)
 920                 {
 921                     omega = RAD2DEG*dih[Om][j];
 922                     mat[(int)((phi*NMAT)/360)+NMAT/2][(int)((psi*NMAT)/360)+NMAT/2]
 923                         += omega;
 924                 }
 925             }
 926             if (bViol)
 927             {
 928                 ffclose(gp);
 929             }
 930             ffclose(fp);
 931             if (bOm)
 932             {
 933                 sprintf(fn, "ramomega%s.xpm", dlist[i].name);
 934                 fp = ffopen(fn, "w");
 935                 lo = hi = 0;
 936                 for (j = 0; (j < NMAT); j++)
 937                 {
 938                     for (k = 0; (k < NMAT); k++)
 939                     {
 940                         mat[j][k] /= nf;
 941                         lo         = min(mat[j][k], lo);
 942                         hi         = max(mat[j][k], hi);
 943                     }
 944                 }
 945                 /* Symmetrise */
 946                 if (fabs(lo) > fabs(hi))
 947                 {
 948                     hi = -lo;
 949                 }
 950                 else
 951                 {
 952                     lo = -hi;
 953                 }
 954                 /* Add 180 */
 955                 for (j = 0; (j < NMAT); j++)
 956                 {
 957                     for (k = 0; (k < NMAT); k++)
 958                     {
 959                         mat[j][k] += 180;
 960                     }
 961                 }
 962                 lo     += 180;
 963                 hi     += 180;
 964                 nlevels = 20;
 965                 write_xpm3(fp, 0, "Omega/Ramachandran Plot", "Deg", "Phi", "Psi",
 966                            NMAT, NMAT, axis, axis, mat, lo, 180.0, hi, rlo, rmid, rhi, &nlevels);
 967                 ffclose(fp);
 968                 for (j = 0; (j < NMAT); j++)
 969                 {
 970                     sfree(mat[j]);
 971                 }
 972                 sfree(mat);
 973             }
 974         }
 975         if ((has_dihedral(edChi1, &(dlist[i]))) &&
 976             (has_dihedral(edChi2, &(dlist[i]))))
 977         {
 978             sprintf(fn, "ramaX1X2%s.xvg", dlist[i].name);
 979             fp = rama_file(fn, "\\8c\\4\\s1\\N-\\8c\\4\\s2\\N Ramachandran Plot",
 980                            "\\8c\\4\\s1\\N (deg)", "\\8c\\4\\s2\\N (deg)", oenv);
 981             Xi1 = dlist[i].j0[edChi1];
 982             Xi2 = dlist[i].j0[edChi2];
 983             for (j = 0; (j < nf); j++)
 984             {
 985                 fprintf(fp, "%10g  %10g\n", RAD2DEG*dih[Xi1][j], RAD2DEG*dih[Xi2][j]);
 986             }
 987             ffclose(fp);
 988         }
 989         else
 990         {
 991             fprintf(stderr, "No chi1 & chi2 angle for %s\n", dlist[i].name);
 992         }
 993     }
 994 }
 995
 996
 997 static void print_transitions(const char *fn, int maxchi, int nlist,
 998                               t_dlist dlist[], real dt,
 999                               const output_env_t oenv)
1000 {
1001     /* based on order_params below */
1002     FILE *fp;
1003     int   nh[edMax];
1004     int   i, Dih, Xi;
1005
1006     /*  must correspond with enum in pp2shift.h:38 */
1007     char *leg[edMax];
1008 #define NLEG asize(leg)
1009
1010     leg[0] = strdup("Phi");
1011     leg[1] = strdup("Psi");
1012     leg[2] = strdup("Omega");
1013     leg[3] = strdup("Chi1");
1014     leg[4] = strdup("Chi2");
1015     leg[5] = strdup("Chi3");
1016     leg[6] = strdup("Chi4");
1017     leg[7] = strdup("Chi5");
1018     leg[8] = strdup("Chi6");
1019
1020     /* Print order parameters */
1021     fp = xvgropen(fn, "Dihedral Rotamer Transitions", "Residue", "Transitions/ns",
1022                   oenv);
1023     xvgr_legend(fp, NONCHI+maxchi, (const char**)leg, oenv);
1024
1025     for (Dih = 0; (Dih < edMax); Dih++)
1026     {
1027         nh[Dih] = 0;
1028     }
1029
1030     fprintf(fp, "%5s ", "#Res.");
1031     fprintf(fp, "%10s %10s %10s ", leg[edPhi], leg[edPsi], leg[edOmega]);
1032     for (Xi = 0; Xi < maxchi; Xi++)
1033     {
1034         fprintf(fp, "%10s ", leg[NONCHI+Xi]);
1035     }
1036     fprintf(fp, "\n");
1037
1038     for (i = 0; (i < nlist); i++)
1039     {
1040         fprintf(fp, "%5d ", dlist[i].resnr);
1041         for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1042         {
1043             fprintf(fp, "%10.3f ", dlist[i].ntr[Dih]/dt);
1044         }
1045         /* fprintf(fp,"%12s\n",dlist[i].name);  this confuses xmgrace */
1046         fprintf(fp, "\n");
1047     }
1048     ffclose(fp);
1049 }
1050
1051 static void order_params(FILE *log,
1052                          const char *fn, int maxchi, int nlist, t_dlist dlist[],
1053                          const char *pdbfn, real bfac_init,
1054                          t_atoms *atoms, rvec x[], int ePBC, matrix box,
1055                          gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, const output_env_t oenv)
1056 {
1057     FILE *fp;
1058     int   nh[edMax];
1059     char  buf[STRLEN];
1060     int   i, Dih, Xi;
1061     real  S2Max, S2Min;
1062
1063     /* except for S2Min/Max, must correspond with enum in pp2shift.h:38 */
1064     const char *const_leg[2+edMax] = {
1065         "S2Min", "S2Max", "Phi", "Psi", "Omega",
1066         "Chi1", "Chi2", "Chi3", "Chi4", "Chi5",
1067         "Chi6"
1068     };
1069 #define NLEG asize(leg)
1070
1071     char *leg[2+edMax];
1072
1073     for (i = 0; i < NLEG; i++)
1074     {
1075         leg[i] = strdup(const_leg[i]);
1076     }
1077
1078     /* Print order parameters */
1079     fp = xvgropen(fn, "Dihedral Order Parameters", "Residue", "S2", oenv);
1080     xvgr_legend(fp, 2+NONCHI+maxchi, const_leg, oenv);
1081
1082     for (Dih = 0; (Dih < edMax); Dih++)
1083     {
1084         nh[Dih] = 0;
1085     }
1086
1087     fprintf(fp, "%5s ", "#Res.");
1088     fprintf(fp, "%10s %10s ", leg[0], leg[1]);
1089     fprintf(fp, "%10s %10s %10s ", leg[2+edPhi], leg[2+edPsi], leg[2+edOmega]);
1090     for (Xi = 0; Xi < maxchi; Xi++)
1091     {
1092         fprintf(fp, "%10s ", leg[2+NONCHI+Xi]);
1093     }
1094     fprintf(fp, "\n");
1095
1096     for (i = 0; (i < nlist); i++)
1097     {
1098         S2Max = -10;
1099         S2Min = 10;
1100         for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1101         {
1102             if (dlist[i].S2[Dih] != 0)
1103             {
1104                 if (dlist[i].S2[Dih] > S2Max)
1105                 {
1106                     S2Max = dlist[i].S2[Dih];
1107                 }
1108                 if (dlist[i].S2[Dih] < S2Min)
1109                 {
1110                     S2Min = dlist[i].S2[Dih];
1111                 }
1112             }
1113             if (dlist[i].S2[Dih] > 0.8)
1114             {
1115                 nh[Dih]++;
1116             }
1117         }
1118         fprintf(fp, "%5d ", dlist[i].resnr);
1119         fprintf(fp, "%10.3f %10.3f ", S2Min, S2Max);
1120         for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1121         {
1122             fprintf(fp, "%10.3f ", dlist[i].S2[Dih]);
1123         }
1124         fprintf(fp, "\n");
1125         /* fprintf(fp,"%12s\n",dlist[i].name);  this confuses xmgrace */
1126     }
1127     ffclose(fp);
1128
1129     if (NULL != pdbfn)
1130     {
1131         real x0, y0, z0;
1132
1133         if (NULL == atoms->pdbinfo)
1134         {
1135             snew(atoms->pdbinfo, atoms->nr);
1136         }
1137         for (i = 0; (i < atoms->nr); i++)
1138         {
1139             atoms->pdbinfo[i].bfac = bfac_init;
1140         }
1141
1142         for (i = 0; (i < nlist); i++)
1143         {
1144             atoms->pdbinfo[dlist[i].atm.N].bfac = -dlist[i].S2[0]; /* Phi */
1145             atoms->pdbinfo[dlist[i].atm.H].bfac = -dlist[i].S2[0]; /* Phi */
1146             atoms->pdbinfo[dlist[i].atm.C].bfac = -dlist[i].S2[1]; /* Psi */
1147             atoms->pdbinfo[dlist[i].atm.O].bfac = -dlist[i].S2[1]; /* Psi */
1148             for (Xi = 0; (Xi < maxchi); Xi++)                      /* Chi's */
1149             {
1150                 if (dlist[i].atm.Cn[Xi+3] != -1)
1151                 {
1152                     atoms->pdbinfo[dlist[i].atm.Cn[Xi+1]].bfac = -dlist[i].S2[NONCHI+Xi];
1153                 }
1154             }
1155         }
1156
1157         fp = ffopen(pdbfn, "w");
1158         fprintf(fp, "REMARK generated by g_chi\n");
1159         fprintf(fp, "REMARK "
1160                 "B-factor field contains negative of dihedral order parameters\n");
1161         write_pdbfile(fp, NULL, atoms, x, ePBC, box, ' ', 0, NULL, TRUE);
1162         x0 = y0 = z0 = 1000.0;
1163         for (i = 0; (i < atoms->nr); i++)
1164         {
1165             x0 = min(x0, x[i][XX]);
1166             y0 = min(y0, x[i][YY]);
1167             z0 = min(z0, x[i][ZZ]);
1168         }
1169         x0 *= 10.0; /* nm -> angstrom */
1170         y0 *= 10.0; /* nm -> angstrom */
1171         z0 *= 10.0; /* nm -> angstrom */
1172         sprintf(buf, "%s%%6.f%%6.2f\n", get_pdbformat());
1173         for (i = 0; (i < 10); i++)
1174         {
1175             fprintf(fp, buf, "ATOM  ", atoms->nr+1+i, "CA", "LEG", ' ',
1176                     atoms->nres+1, ' ', x0, y0, z0+(1.2*i), 0.0, -0.1*i);
1177         }
1178         ffclose(fp);
1179     }
1180
1181     fprintf(log, "Dihedrals with S2 > 0.8\n");
1182     fprintf(log, "Dihedral: ");
1183     if (bPhi)
1184     {
1185         fprintf(log, " Phi  ");
1186     }
1187     if (bPsi)
1188     {
1189         fprintf(log, " Psi ");
1190     }
1191     if (bChi)
1192     {
1193         for (Xi = 0; (Xi < maxchi); Xi++)
1194         {
1195             fprintf(log, " %s ", leg[2+NONCHI+Xi]);
1196         }
1197     }
1198     fprintf(log, "\nNumber:   ");
1199     if (bPhi)
1200     {
1201         fprintf(log, "%4d  ", nh[0]);
1202     }
1203     if (bPsi)
1204     {
1205         fprintf(log, "%4d  ", nh[1]);
1206     }
1207     if (bChi)
1208     {
1209         for (Xi = 0; (Xi < maxchi); Xi++)
1210         {
1211             fprintf(log, "%4d  ", nh[NONCHI+Xi]);
1212         }
1213     }
1214     fprintf(log, "\n");
1215
1216     for (i = 0; i < NLEG; i++)
1217     {
1218         sfree(leg[i]);
1219     }
1220
1221 }
1222
1223 int gmx_chi(int argc, char *argv[])
1224 {
1225     const char *desc[] = {
1226         "[TT]g_chi[tt] computes [GRK]phi[grk], [GRK]psi[grk], [GRK]omega[grk], and [GRK]chi[grk] dihedrals for all your ",
1227         "amino acid backbone and sidechains.",
1228         "It can compute dihedral angle as a function of time, and as",
1229         "histogram distributions.",
1230         "The distributions [TT](histo-(dihedral)(RESIDUE).xvg[tt]) are cumulative over all residues of each type.[PAR]",
1231         "If option [TT]-corr[tt] is given, the program will",
1232         "calculate dihedral autocorrelation functions. The function used",
1233         "is C(t) = [CHEVRON][COS][GRK]chi[grk]([GRK]tau[grk])[cos] [COS][GRK]chi[grk]([GRK]tau[grk]+t)[cos][chevron]. The use of cosines",
1234         "rather than angles themselves, resolves the problem of periodicity.",
1235         "(Van der Spoel & Berendsen (1997), Biophys. J. 72, 2032-2041).",
1236         "Separate files for each dihedral of each residue",
1237         "[TT](corr(dihedral)(RESIDUE)(nresnr).xvg[tt]) are output, as well as a",
1238         "file containing the information for all residues (argument of [TT]-corr[tt]).[PAR]",
1239         "With option [TT]-all[tt], the angles themselves as a function of time for",
1240         "each residue are printed to separate files [TT](dihedral)(RESIDUE)(nresnr).xvg[tt].",
1241         "These can be in radians or degrees.[PAR]",
1242         "A log file (argument [TT]-g[tt]) is also written. This contains [BR]",
1243         "(a) information about the number of residues of each type.[BR]",
1244         "(b) The NMR ^3J coupling constants from the Karplus equation.[BR]",
1245         "(c) a table for each residue of the number of transitions between ",
1246         "rotamers per nanosecond,  and the order parameter S^2 of each dihedral.[BR]",
1247         "(d) a table for each residue of the rotamer occupancy.[PAR]",
1248         "All rotamers are taken as 3-fold, except for [GRK]omega[grk] and [GRK]chi[grk] dihedrals",
1249         "to planar groups (i.e. [GRK]chi[grk][SUB]2[sub] of aromatics, Asp and Asn; [GRK]chi[grk][SUB]3[sub] of Glu",
1250         "and Gln; and [GRK]chi[grk][SUB]4[sub] of Arg), which are 2-fold. \"rotamer 0\" means ",
1251         "that the dihedral was not in the core region of each rotamer. ",
1252         "The width of the core region can be set with [TT]-core_rotamer[tt][PAR]",
1253
1254         "The S^2 order parameters are also output to an [TT].xvg[tt] file",
1255         "(argument [TT]-o[tt] ) and optionally as a [TT].pdb[tt] file with",
1256         "the S^2 values as B-factor (argument [TT]-p[tt]). ",
1257         "The total number of rotamer transitions per timestep",
1258         "(argument [TT]-ot[tt]), the number of transitions per rotamer",
1259         "(argument [TT]-rt[tt]), and the ^3J couplings (argument [TT]-jc[tt]), ",
1260         "can also be written to [TT].xvg[tt] files. Note that the analysis",
1261         "of rotamer transitions assumes that the supplied trajectory frames",
1262         "are equally spaced in time.[PAR]",
1263
1264         "If [TT]-chi_prod[tt] is set (and [TT]-maxchi[tt] > 0), cumulative rotamers, e.g.",
1265         "1+9([GRK]chi[grk][SUB]1[sub]-1)+3([GRK]chi[grk][SUB]2[sub]-1)+([GRK]chi[grk][SUB]3[sub]-1) (if the residue has three 3-fold ",
1266         "dihedrals and [TT]-maxchi[tt] >= 3)",
1267         "are calculated. As before, if any dihedral is not in the core region,",
1268         "the rotamer is taken to be 0. The occupancies of these cumulative ",
1269         "rotamers (starting with rotamer 0) are written to the file",
1270         "that is the argument of [TT]-cp[tt], and if the [TT]-all[tt] flag",
1271         "is given, the rotamers as functions of time",
1272         "are written to [TT]chiproduct(RESIDUE)(nresnr).xvg[tt] ",
1273         "and their occupancies to [TT]histo-chiproduct(RESIDUE)(nresnr).xvg[tt].[PAR]",
1274
1275         "The option [TT]-r[tt] generates a contour plot of the average [GRK]omega[grk] angle",
1276         "as a function of the [GRK]phi[grk] and [GRK]psi[grk] angles, that is, in a Ramachandran plot",
1277         "the average [GRK]omega[grk] angle is plotted using color coding.",
1278
1279     };
1280
1281     const char *bugs[] = {
1282         "Produces MANY output files (up to about 4 times the number of residues in the protein, twice that if autocorrelation functions are calculated). Typically several hundred files are output.",
1283         "[GRK]phi[grk] and [GRK]psi[grk] dihedrals are calculated in a non-standard way, using H-N-CA-C for [GRK]phi[grk] instead of C(-)-N-CA-C, and N-CA-C-O for [GRK]psi[grk] instead of N-CA-C-N(+). This causes (usually small) discrepancies with the output of other tools like [TT]g_rama[tt].",
1284         "[TT]-r0[tt] option does not work properly",
1285         "Rotamers with multiplicity 2 are printed in [TT]chi.log[tt] as if they had multiplicity 3, with the 3rd (g(+)) always having probability 0"
1286     };
1287
1288     /* defaults */
1289     static int         r0          = 1, ndeg = 1, maxchi = 2;
1290     static gmx_bool    bAll        = FALSE;
1291     static gmx_bool    bPhi        = FALSE, bPsi = FALSE, bOmega = FALSE;
1292     static real        bfac_init   = -1.0, bfac_max = 0;
1293     static const char *maxchistr[] = { NULL, "0", "1", "2", "3",  "4", "5", "6", NULL };
1294     static gmx_bool    bRama       = FALSE, bShift = FALSE, bViol = FALSE, bRamOmega = FALSE;
1295     static gmx_bool    bNormHisto  = TRUE, bChiProduct = FALSE, bHChi = FALSE, bRAD = FALSE, bPBC = TRUE;
1296     static real        core_frac   = 0.5;
1297     t_pargs            pa[]        = {
1298         { "-r0",  FALSE, etINT, {&r0},
1299           "starting residue" },
1300         { "-phi",  FALSE, etBOOL, {&bPhi},
1301           "Output for [GRK]phi[grk] dihedral angles" },
1302         { "-psi",  FALSE, etBOOL, {&bPsi},
1303           "Output for [GRK]psi[grk] dihedral angles" },
1304         { "-omega", FALSE, etBOOL, {&bOmega},
1305           "Output for [GRK]omega[grk] dihedrals (peptide bonds)" },
1306         { "-rama", FALSE, etBOOL, {&bRama},
1307           "Generate [GRK]phi[grk]/[GRK]psi[grk] and [GRK]chi[grk][SUB]1[sub]/[GRK]chi[grk][SUB]2[sub] Ramachandran plots" },
1308         { "-viol", FALSE, etBOOL, {&bViol},
1309           "Write a file that gives 0 or 1 for violated Ramachandran angles" },
1310         { "-periodic", FALSE, etBOOL, {&bPBC},
1311           "Print dihedral angles modulo 360 degrees" },
1312         { "-all",  FALSE, etBOOL, {&bAll},
1313           "Output separate files for every dihedral." },
1314         { "-rad",  FALSE, etBOOL, {&bRAD},
1315           "in angle vs time files, use radians rather than degrees."},
1316         { "-shift", FALSE, etBOOL, {&bShift},
1317           "Compute chemical shifts from [GRK]phi[grk]/[GRK]psi[grk] angles" },
1318         { "-binwidth", FALSE, etINT, {&ndeg},
1319           "bin width for histograms (degrees)" },
1320         { "-core_rotamer", FALSE, etREAL, {&core_frac},
1321           "only the central [TT]-core_rotamer[tt]*(360/multiplicity) belongs to each rotamer (the rest is assigned to rotamer 0)" },
1322         { "-maxchi", FALSE, etENUM, {maxchistr},
1323           "calculate first ndih [GRK]chi[grk] dihedrals" },
1324         { "-normhisto", FALSE, etBOOL, {&bNormHisto},
1325           "Normalize histograms" },
1326         { "-ramomega", FALSE, etBOOL, {&bRamOmega},
1327           "compute average omega as a function of [GRK]phi[grk]/[GRK]psi[grk] and plot it in an [TT].xpm[tt] plot" },
1328         { "-bfact", FALSE, etREAL, {&bfac_init},
1329           "B-factor value for [TT].pdb[tt] file for atoms with no calculated dihedral order parameter"},
1330         { "-chi_prod", FALSE, etBOOL, {&bChiProduct},
1331           "compute a single cumulative rotamer for each residue"},
1332         { "-HChi", FALSE, etBOOL, {&bHChi},
1333           "Include dihedrals to sidechain hydrogens"},
1334         { "-bmax",  FALSE, etREAL, {&bfac_max},
1335           "Maximum B-factor on any of the atoms that make up a dihedral, for the dihedral angle to be considere in the statistics. Applies to database work where a number of X-Ray structures is analyzed. [TT]-bmax[tt] <= 0 means no limit." }
1336     };
1337
1338     FILE              *log;
1339     int                natoms, nlist, idum, nbin;
1340     t_atoms            atoms;
1341     rvec              *x;
1342     int                ePBC;
1343     matrix             box;
1344     char               title[256], grpname[256];
1345     t_dlist           *dlist;
1346     gmx_bool           bChi, bCorr, bSSHisto;
1347     gmx_bool           bDo_rt, bDo_oh, bDo_ot, bDo_jc;
1348     real               dt = 0, traj_t_ns;
1349     output_env_t       oenv;
1350     gmx_residuetype_t  rt;
1351
1352     atom_id            isize, *index;
1353     int                ndih, nactdih, nf;
1354     real             **dih, *trans_frac, *aver_angle, *time;
1355     int                i, j, **chi_lookup, *multiplicity;
1356
1357     t_filenm           fnm[] = {
1358         { efSTX, "-s",  NULL,     ffREAD  },
1359         { efTRX, "-f",  NULL,     ffREAD  },
1360         { efXVG, "-o",  "order",  ffWRITE },
1361         { efPDB, "-p",  "order",  ffOPTWR },
1362         { efDAT, "-ss", "ssdump", ffOPTRD },
1363         { efXVG, "-jc", "Jcoupling", ffWRITE },
1364         { efXVG, "-corr",  "dihcorr", ffOPTWR },
1365         { efLOG, "-g",  "chi",    ffWRITE },
1366         /* add two more arguments copying from g_angle */
1367         { efXVG, "-ot", "dihtrans", ffOPTWR },
1368         { efXVG, "-oh", "trhisto",  ffOPTWR },
1369         { efXVG, "-rt", "restrans",  ffOPTWR },
1370         { efXVG, "-cp", "chiprodhisto",  ffOPTWR }
1371     };
1372 #define NFILE asize(fnm)
1373     int                npargs;
1374     t_pargs           *ppa;
1375
1376     npargs = asize(pa);
1377     ppa    = add_acf_pargs(&npargs, pa);
1378     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
1379                            NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs,
1380                            &oenv))
1381     {
1382         return 0;
1383     }
1384
1385     /* Handle result from enumerated type */
1386     sscanf(maxchistr[0], "%d", &maxchi);
1387     bChi = (maxchi > 0);
1388
1389     log = ffopen(ftp2fn(efLOG, NFILE, fnm), "w");
1390
1391     if (bRamOmega)
1392     {
1393         bOmega = TRUE;
1394         bPhi   = TRUE;
1395         bPsi   = TRUE;
1396     }
1397
1398     /* set some options */
1399     bDo_rt = (opt2bSet("-rt", NFILE, fnm));
1400     bDo_oh = (opt2bSet("-oh", NFILE, fnm));
1401     bDo_ot = (opt2bSet("-ot", NFILE, fnm));
1402     bDo_jc = (opt2bSet("-jc", NFILE, fnm));
1403     bCorr  = (opt2bSet("-corr", NFILE, fnm));
1404     if (bCorr)
1405     {
1406         fprintf(stderr, "Will calculate autocorrelation\n");
1407     }
1408
1409     if (core_frac > 1.0)
1410     {
1411         fprintf(stderr, "core_rotamer fraction > 1.0 ; will use 1.0\n");
1412         core_frac = 1.0;
1413     }
1414     if (core_frac < 0.0)
1415     {
1416         fprintf(stderr, "core_rotamer fraction < 0.0 ; will use 0.0\n");
1417         core_frac = 0.0;
1418     }
1419
1420     if (maxchi > MAXCHI)
1421     {
1422         fprintf(stderr,
1423                 "Will only calculate first %d Chi dihedrals in stead of %d.\n",
1424                 MAXCHI, maxchi);
1425         maxchi = MAXCHI;
1426     }
1427     bSSHisto = ftp2bSet(efDAT, NFILE, fnm);
1428     nbin     = 360/ndeg;
1429
1430     /* Find the chi angles using atoms struct and a list of amino acids */
1431     get_stx_coordnum(ftp2fn(efSTX, NFILE, fnm), &natoms);
1432     init_t_atoms(&atoms, natoms, TRUE);
1433     snew(x, natoms);
1434     read_stx_conf(ftp2fn(efSTX, NFILE, fnm), title, &atoms, x, NULL, &ePBC, box);
1435     fprintf(log, "Title: %s\n", title);
1436
1437     gmx_residuetype_init(&rt);
1438     dlist = mk_dlist(log, &atoms, &nlist, bPhi, bPsi, bChi, bHChi, maxchi, r0, rt);
1439     fprintf(stderr, "%d residues with dihedrals found\n", nlist);
1440
1441     if (nlist == 0)
1442     {
1443         gmx_fatal(FARGS, "No dihedrals in your structure!\n");
1444     }
1445
1446     /* Make a linear index for reading all. */
1447     index = make_chi_ind(nlist, dlist, &ndih);
1448     isize = 4*ndih;
1449     fprintf(stderr, "%d dihedrals found\n", ndih);
1450
1451     snew(dih, ndih);
1452
1453     /* COMPUTE ALL DIHEDRALS! */
1454     read_ang_dih(ftp2fn(efTRX, NFILE, fnm), FALSE, TRUE, FALSE, bPBC, 1, &idum,
1455                  &nf, &time, isize, index, &trans_frac, &aver_angle, dih, oenv);
1456
1457     dt = (time[nf-1]-time[0])/(nf-1); /* might want this for corr or n. transit*/
1458     if (bCorr)
1459     {
1460         if (nf < 2)
1461         {
1462             gmx_fatal(FARGS, "Need at least 2 frames for correlation");
1463         }
1464     }
1465
1466     /* put angles in -M_PI to M_PI ! and correct phase factor for phi and psi
1467      * pass nactdih instead of ndih to low_ana_dih_trans and get_chi_product_traj
1468      * to prevent accessing off end of arrays when maxchi < 5 or 6. */
1469     nactdih = reset_em_all(nlist, dlist, nf, dih, maxchi);
1470
1471     if (bAll)
1472     {
1473         dump_em_all(nlist, dlist, nf, time, dih, maxchi, bPhi, bPsi, bChi, bOmega, bRAD, oenv);
1474     }
1475
1476     /* Histogramming & J coupling constants & calc of S2 order params */
1477     histogramming(log, nbin, rt, nf, maxchi, dih, nlist, dlist, index,
1478                   bPhi, bPsi, bOmega, bChi,
1479                   bNormHisto, bSSHisto, ftp2fn(efDAT, NFILE, fnm), bfac_max, &atoms,
1480                   bDo_jc, opt2fn("-jc", NFILE, fnm), oenv);
1481
1482     /* transitions
1483      *
1484      * added multiplicity */
1485
1486     snew(multiplicity, ndih);
1487     mk_multiplicity_lookup(multiplicity, maxchi, nlist, dlist, ndih);
1488
1489     strcpy(grpname, "All residues, ");
1490     if (bPhi)
1491     {
1492         strcat(grpname, "Phi ");
1493     }
1494     if (bPsi)
1495     {
1496         strcat(grpname, "Psi ");
1497     }
1498     if (bOmega)
1499     {
1500         strcat(grpname, "Omega ");
1501     }
1502     if (bChi)
1503     {
1504         strcat(grpname, "Chi 1-");
1505         sprintf(grpname + strlen(grpname), "%i", maxchi);
1506     }
1507
1508
1509     low_ana_dih_trans(bDo_ot, opt2fn("-ot", NFILE, fnm),
1510                       bDo_oh, opt2fn("-oh", NFILE, fnm), maxchi,
1511                       dih, nlist, dlist, nf, nactdih, grpname, multiplicity,
1512                       time, FALSE, core_frac, oenv);
1513
1514     /* Order parameters */
1515     order_params(log, opt2fn("-o", NFILE, fnm), maxchi, nlist, dlist,
1516                  ftp2fn_null(efPDB, NFILE, fnm), bfac_init,
1517                  &atoms, x, ePBC, box, bPhi, bPsi, bChi, oenv);
1518
1519     /* Print ramachandran maps! */
1520     if (bRama)
1521     {
1522         do_rama(nf, nlist, dlist, dih, bViol, bRamOmega, oenv);
1523     }
1524
1525     if (bShift)
1526     {
1527         do_pp2shifts(log, nf, nlist, dlist, dih);
1528     }
1529
1530     /* rprint S^2, transitions, and rotamer occupancies to log */
1531     traj_t_ns = 0.001 * (time[nf-1]-time[0]);
1532     pr_dlist(log, nlist, dlist, traj_t_ns, edPrintST, bPhi, bPsi, bChi, bOmega, maxchi);
1533     pr_dlist(log, nlist, dlist, traj_t_ns, edPrintRO, bPhi, bPsi, bChi, bOmega, maxchi);
1534     ffclose(log);
1535     /* transitions to xvg */
1536     if (bDo_rt)
1537     {
1538         print_transitions(opt2fn("-rt", NFILE, fnm), maxchi, nlist, dlist, traj_t_ns, oenv);
1539     }
1540
1541     /* chi_product trajectories (ie one "rotamer number" for each residue) */
1542     if (bChiProduct && bChi)
1543     {
1544         snew(chi_lookup, nlist);
1545         for (i = 0; i < nlist; i++)
1546         {
1547             snew(chi_lookup[i], maxchi);
1548         }
1549         mk_chi_lookup(chi_lookup, maxchi, nlist, dlist);
1550
1551         get_chi_product_traj(dih, nf, nactdih,
1552                              maxchi, dlist, time, chi_lookup, multiplicity,
1553                              FALSE, bNormHisto, core_frac, bAll,
1554                              opt2fn("-cp", NFILE, fnm), oenv);
1555
1556         for (i = 0; i < nlist; i++)
1557         {
1558             sfree(chi_lookup[i]);
1559         }
1560     }
1561
1562     /* Correlation comes last because it fucks up the angles */
1563     if (bCorr)
1564     {
1565         do_dihcorr(opt2fn("-corr", NFILE, fnm), nf, ndih, dih, dt, nlist, dlist, time,
1566                    maxchi, bPhi, bPsi, bChi, bOmega, oenv);
1567     }
1568
1569
1570     do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1571     do_view(oenv, opt2fn("-jc", NFILE, fnm), "-nxy");
1572     if (bCorr)
1573     {
1574         do_view(oenv, opt2fn("-corr", NFILE, fnm), "-nxy");
1575     }
1576
1577     gmx_residuetype_destroy(rt);
1578
1579     return 0;
1580 }