src/tools/gmx_chi.c

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