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