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