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