src/gromacs/gmxana/gmx_chi.cpp

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