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,2017, 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, const 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] = {nullptr, nullptr, nullptr};
 483     const char *sss[3] = { "sheet", "helix", "coil" };
 484     real        S2;
 485     real       *normhisto;
 486     real      **Jc, **Jcsig;
 487     int     ****his_aa_ss = nullptr;
 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 = nullptr;
 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, nullptr, &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 = nullptr;
 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  = nullptr, 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, const 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 (nullptr != pdbfn)
1144     {
1145         real x0, y0, z0;
1146
1147         atoms->havePdbInfo = TRUE;
1148
1149         if (nullptr == atoms->pdbinfo)
1150         {
1151             snew(atoms->pdbinfo, atoms->nr);
1152         }
1153         for (i = 0; (i < atoms->nr); i++)
1154         {
1155             atoms->pdbinfo[i].bfac = bfac_init;
1156         }
1157
1158         for (i = 0; (i < nlist); i++)
1159         {
1160             atoms->pdbinfo[dlist[i].atm.N].bfac = -dlist[i].S2[0]; /* Phi */
1161             atoms->pdbinfo[dlist[i].atm.H].bfac = -dlist[i].S2[0]; /* Phi */
1162             atoms->pdbinfo[dlist[i].atm.C].bfac = -dlist[i].S2[1]; /* Psi */
1163             atoms->pdbinfo[dlist[i].atm.O].bfac = -dlist[i].S2[1]; /* Psi */
1164             for (Xi = 0; (Xi < maxchi); Xi++)                      /* Chi's */
1165             {
1166                 if (dlist[i].atm.Cn[Xi+3] != -1)
1167                 {
1168                     atoms->pdbinfo[dlist[i].atm.Cn[Xi+1]].bfac = -dlist[i].S2[NONCHI+Xi];
1169                 }
1170             }
1171         }
1172
1173         fp = gmx_ffopen(pdbfn, "w");
1174         fprintf(fp, "REMARK generated by g_chi\n");
1175         fprintf(fp, "REMARK "
1176                 "B-factor field contains negative of dihedral order parameters\n");
1177         write_pdbfile(fp, nullptr, atoms, x, ePBC, box, ' ', 0, nullptr, TRUE);
1178         x0 = y0 = z0 = 1000.0;
1179         for (i = 0; (i < atoms->nr); i++)
1180         {
1181             x0 = std::min(x0, x[i][XX]);
1182             y0 = std::min(y0, x[i][YY]);
1183             z0 = std::min(z0, x[i][ZZ]);
1184         }
1185         x0 *= 10.0; /* nm -> angstrom */
1186         y0 *= 10.0; /* nm -> angstrom */
1187         z0 *= 10.0; /* nm -> angstrom */
1188         for (i = 0; (i < 10); i++)
1189         {
1190             gmx_fprintf_pdb_atomline(fp, epdbATOM, atoms->nr+1+i, "CA", ' ', "LEG", ' ', atoms->nres+1, ' ',
1191                                      x0, y0, z0+(1.2*i), 0.0, -0.1*i, "");
1192         }
1193         gmx_ffclose(fp);
1194     }
1195
1196     fprintf(log, "Dihedrals with S2 > 0.8\n");
1197     fprintf(log, "Dihedral: ");
1198     if (bPhi)
1199     {
1200         fprintf(log, " Phi  ");
1201     }
1202     if (bPsi)
1203     {
1204         fprintf(log, " Psi ");
1205     }
1206     if (bChi)
1207     {
1208         for (Xi = 0; (Xi < maxchi); Xi++)
1209         {
1210             fprintf(log, " %s ", leg[2+NONCHI+Xi]);
1211         }
1212     }
1213     fprintf(log, "\nNumber:   ");
1214     if (bPhi)
1215     {
1216         fprintf(log, "%4d  ", nh[0]);
1217     }
1218     if (bPsi)
1219     {
1220         fprintf(log, "%4d  ", nh[1]);
1221     }
1222     if (bChi)
1223     {
1224         for (Xi = 0; (Xi < maxchi); Xi++)
1225         {
1226             fprintf(log, "%4d  ", nh[NONCHI+Xi]);
1227         }
1228     }
1229     fprintf(log, "\n");
1230
1231     for (i = 0; i < NLEG; i++)
1232     {
1233         sfree(leg[i]);
1234     }
1235
1236 }
1237
1238 int gmx_chi(int argc, char *argv[])
1239 {
1240     const char *desc[] = {
1241         "[THISMODULE] computes [GRK]phi[grk], [GRK]psi[grk], [GRK]omega[grk],",
1242         "and [GRK]chi[grk] dihedrals for all your",
1243         "amino acid backbone and sidechains.",
1244         "It can compute dihedral angle as a function of time, and as",
1245         "histogram distributions.",
1246         "The distributions [TT](histo-(dihedral)(RESIDUE).xvg[tt]) are cumulative over all residues of each type.[PAR]",
1247         "If option [TT]-corr[tt] is given, the program will",
1248         "calculate dihedral autocorrelation functions. The function used",
1249         "is C(t) = [CHEVRON][COS][GRK]chi[grk]([GRK]tau[grk])[cos] [COS][GRK]chi[grk]([GRK]tau[grk]+t)[cos][chevron]. The use of cosines",
1250         "rather than angles themselves, resolves the problem of periodicity.",
1251         "(Van der Spoel & Berendsen (1997), Biophys. J. 72, 2032-2041).",
1252         "Separate files for each dihedral of each residue",
1253         "[TT](corr(dihedral)(RESIDUE)(nresnr).xvg[tt]) are output, as well as a",
1254         "file containing the information for all residues (argument of [TT]-corr[tt]).[PAR]",
1255         "With option [TT]-all[tt], the angles themselves as a function of time for",
1256         "each residue are printed to separate files [TT](dihedral)(RESIDUE)(nresnr).xvg[tt].",
1257         "These can be in radians or degrees.[PAR]",
1258         "A log file (argument [TT]-g[tt]) is also written. This contains",
1259         "",
1260         " * information about the number of residues of each type.",
1261         " * The NMR ^3J coupling constants from the Karplus equation.",
1262         " * a table for each residue of the number of transitions between ",
1263         "   rotamers per nanosecond,  and the order parameter S^2 of each dihedral.",
1264         " * a table for each residue of the rotamer occupancy.",
1265         "",
1266         "All rotamers are taken as 3-fold, except for [GRK]omega[grk] and [GRK]chi[grk] dihedrals",
1267         "to planar groups (i.e. [GRK]chi[grk][SUB]2[sub] of aromatics, Asp and Asn; [GRK]chi[grk][SUB]3[sub] of Glu",
1268         "and Gln; and [GRK]chi[grk][SUB]4[sub] of Arg), which are 2-fold. \"rotamer 0\" means ",
1269         "that the dihedral was not in the core region of each rotamer. ",
1270         "The width of the core region can be set with [TT]-core_rotamer[tt][PAR]",
1271
1272         "The S^2 order parameters are also output to an [REF].xvg[ref] file",
1273         "(argument [TT]-o[tt] ) and optionally as a [REF].pdb[ref] file with",
1274         "the S^2 values as B-factor (argument [TT]-p[tt]). ",
1275         "The total number of rotamer transitions per timestep",
1276         "(argument [TT]-ot[tt]), the number of transitions per rotamer",
1277         "(argument [TT]-rt[tt]), and the ^3J couplings (argument [TT]-jc[tt]), ",
1278         "can also be written to [REF].xvg[ref] files. Note that the analysis",
1279         "of rotamer transitions assumes that the supplied trajectory frames",
1280         "are equally spaced in time.[PAR]",
1281
1282         "If [TT]-chi_prod[tt] is set (and [TT]-maxchi[tt] > 0), cumulative rotamers, e.g.",
1283         "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 ",
1284         "dihedrals and [TT]-maxchi[tt] >= 3)",
1285         "are calculated. As before, if any dihedral is not in the core region,",
1286         "the rotamer is taken to be 0. The occupancies of these cumulative ",
1287         "rotamers (starting with rotamer 0) are written to the file",
1288         "that is the argument of [TT]-cp[tt], and if the [TT]-all[tt] flag",
1289         "is given, the rotamers as functions of time",
1290         "are written to [TT]chiproduct(RESIDUE)(nresnr).xvg[tt] ",
1291         "and their occupancies to [TT]histo-chiproduct(RESIDUE)(nresnr).xvg[tt].[PAR]",
1292
1293         "The option [TT]-r[tt] generates a contour plot of the average [GRK]omega[grk] angle",
1294         "as a function of the [GRK]phi[grk] and [GRK]psi[grk] angles, that is, in a Ramachandran plot",
1295         "the average [GRK]omega[grk] angle is plotted using color coding.",
1296
1297     };
1298
1299     const char *bugs[] = {
1300         "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.",
1301         "[GRK]phi[grk] and [GRK]psi[grk] dihedrals are calculated in a "
1302         "non-standard way, using H-N-CA-C for [GRK]phi[grk] instead of "
1303         "C(-)-N-CA-C, and N-CA-C-O for [GRK]psi[grk] instead of N-CA-C-N(+). "
1304         "This causes (usually small) discrepancies with the output of other "
1305         "tools like [gmx-rama].",
1306         "[TT]-r0[tt] option does not work properly",
1307         "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"
1308     };
1309
1310     /* defaults */
1311     static int         r0          = 1, ndeg = 1, maxchi = 2;
1312     static gmx_bool    bAll        = FALSE;
1313     static gmx_bool    bPhi        = FALSE, bPsi = FALSE, bOmega = FALSE;
1314     static real        bfac_init   = -1.0, bfac_max = 0;
1315     static const char *maxchistr[] = { nullptr, "0", "1", "2", "3",  "4", "5", "6", nullptr };
1316     static gmx_bool    bRama       = FALSE, bShift = FALSE, bViol = FALSE, bRamOmega = FALSE;
1317     static gmx_bool    bNormHisto  = TRUE, bChiProduct = FALSE, bHChi = FALSE, bRAD = FALSE, bPBC = TRUE;
1318     static real        core_frac   = 0.5;
1319     t_pargs            pa[]        = {
1320         { "-r0",  FALSE, etINT, {&r0},
1321           "starting residue" },
1322         { "-phi",  FALSE, etBOOL, {&bPhi},
1323           "Output for [GRK]phi[grk] dihedral angles" },
1324         { "-psi",  FALSE, etBOOL, {&bPsi},
1325           "Output for [GRK]psi[grk] dihedral angles" },
1326         { "-omega", FALSE, etBOOL, {&bOmega},
1327           "Output for [GRK]omega[grk] dihedrals (peptide bonds)" },
1328         { "-rama", FALSE, etBOOL, {&bRama},
1329           "Generate [GRK]phi[grk]/[GRK]psi[grk] and [GRK]chi[grk][SUB]1[sub]/[GRK]chi[grk][SUB]2[sub] Ramachandran plots" },
1330         { "-viol", FALSE, etBOOL, {&bViol},
1331           "Write a file that gives 0 or 1 for violated Ramachandran angles" },
1332         { "-periodic", FALSE, etBOOL, {&bPBC},
1333           "Print dihedral angles modulo 360 degrees" },
1334         { "-all",  FALSE, etBOOL, {&bAll},
1335           "Output separate files for every dihedral." },
1336         { "-rad",  FALSE, etBOOL, {&bRAD},
1337           "in angle vs time files, use radians rather than degrees."},
1338         { "-shift", FALSE, etBOOL, {&bShift},
1339           "Compute chemical shifts from [GRK]phi[grk]/[GRK]psi[grk] angles" },
1340         { "-binwidth", FALSE, etINT, {&ndeg},
1341           "bin width for histograms (degrees)" },
1342         { "-core_rotamer", FALSE, etREAL, {&core_frac},
1343           "only the central [TT]-core_rotamer[tt]\\*(360/multiplicity) belongs to each rotamer (the rest is assigned to rotamer 0)" },
1344         { "-maxchi", FALSE, etENUM, {maxchistr},
1345           "calculate first ndih [GRK]chi[grk] dihedrals" },
1346         { "-normhisto", FALSE, etBOOL, {&bNormHisto},
1347           "Normalize histograms" },
1348         { "-ramomega", FALSE, etBOOL, {&bRamOmega},
1349           "compute average omega as a function of [GRK]phi[grk]/[GRK]psi[grk] and plot it in an [REF].xpm[ref] plot" },
1350         { "-bfact", FALSE, etREAL, {&bfac_init},
1351           "B-factor value for [REF].pdb[ref] file for atoms with no calculated dihedral order parameter"},
1352         { "-chi_prod", FALSE, etBOOL, {&bChiProduct},
1353           "compute a single cumulative rotamer for each residue"},
1354         { "-HChi", FALSE, etBOOL, {&bHChi},
1355           "Include dihedrals to sidechain hydrogens"},
1356         { "-bmax",  FALSE, etREAL, {&bfac_max},
1357           "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." }
1358     };
1359
1360     FILE              *log;
1361     int                nlist, idum, nbin;
1362     rvec              *x;
1363     int                ePBC;
1364     matrix             box;
1365     char               grpname[256];
1366     t_dlist           *dlist;
1367     gmx_bool           bChi, bCorr, bSSHisto;
1368     gmx_bool           bDo_rt, bDo_oh, bDo_ot, bDo_jc;
1369     real               dt = 0, traj_t_ns;
1370     gmx_output_env_t  *oenv;
1371     gmx_residuetype_t *rt;
1372
1373     int                isize, *index;
1374     int                ndih, nactdih, nf;
1375     real             **dih, *trans_frac, *aver_angle, *time;
1376     int                i, **chi_lookup, *multiplicity;
1377
1378     t_filenm           fnm[] = {
1379         { efSTX, "-s",  nullptr,     ffREAD  },
1380         { efTRX, "-f",  nullptr,     ffREAD  },
1381         { efXVG, "-o",  "order",  ffWRITE },
1382         { efPDB, "-p",  "order",  ffOPTWR },
1383         { efDAT, "-ss", "ssdump", ffOPTRD },
1384         { efXVG, "-jc", "Jcoupling", ffWRITE },
1385         { efXVG, "-corr",  "dihcorr", ffOPTWR },
1386         { efLOG, "-g",  "chi",    ffWRITE },
1387         /* add two more arguments copying from g_angle */
1388         { efXVG, "-ot", "dihtrans", ffOPTWR },
1389         { efXVG, "-oh", "trhisto",  ffOPTWR },
1390         { efXVG, "-rt", "restrans",  ffOPTWR },
1391         { efXVG, "-cp", "chiprodhisto",  ffOPTWR }
1392     };
1393 #define NFILE asize(fnm)
1394     int                npargs;
1395     t_pargs           *ppa;
1396
1397     npargs = asize(pa);
1398     ppa    = add_acf_pargs(&npargs, pa);
1399     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
1400                            NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs,
1401                            &oenv))
1402     {
1403         sfree(ppa);
1404         return 0;
1405     }
1406
1407     /* Handle result from enumerated type */
1408     sscanf(maxchistr[0], "%d", &maxchi);
1409     bChi = (maxchi > 0);
1410
1411     log = gmx_ffopen(ftp2fn(efLOG, NFILE, fnm), "w");
1412
1413     if (bRamOmega)
1414     {
1415         bOmega = TRUE;
1416         bPhi   = TRUE;
1417         bPsi   = TRUE;
1418     }
1419
1420     /* set some options */
1421     bDo_rt = (opt2bSet("-rt", NFILE, fnm));
1422     bDo_oh = (opt2bSet("-oh", NFILE, fnm));
1423     bDo_ot = (opt2bSet("-ot", NFILE, fnm));
1424     bDo_jc = (opt2bSet("-jc", NFILE, fnm));
1425     bCorr  = (opt2bSet("-corr", NFILE, fnm));
1426     if (bCorr)
1427     {
1428         fprintf(stderr, "Will calculate autocorrelation\n");
1429     }
1430
1431     if (core_frac > 1.0)
1432     {
1433         fprintf(stderr, "core_rotamer fraction > 1.0 ; will use 1.0\n");
1434         core_frac = 1.0;
1435     }
1436     if (core_frac < 0.0)
1437     {
1438         fprintf(stderr, "core_rotamer fraction < 0.0 ; will use 0.0\n");
1439         core_frac = 0.0;
1440     }
1441
1442     if (maxchi > MAXCHI)
1443     {
1444         fprintf(stderr,
1445                 "Will only calculate first %d Chi dihedrals in stead of %d.\n",
1446                 MAXCHI, maxchi);
1447         maxchi = MAXCHI;
1448     }
1449     bSSHisto = ftp2bSet(efDAT, NFILE, fnm);
1450     nbin     = 360/ndeg;
1451
1452     /* Find the chi angles using atoms struct and a list of amino acids */
1453     t_topology *top;
1454     snew(top, 1);
1455     read_tps_conf(ftp2fn(efSTX, NFILE, fnm), top, &ePBC, &x, nullptr, box, FALSE);
1456     t_atoms    &atoms = top->atoms;
1457     if (atoms.pdbinfo == nullptr)
1458     {
1459         snew(atoms.pdbinfo, atoms.nr);
1460     }
1461     fprintf(log, "Title: %s\n", *top->name);
1462
1463     gmx_residuetype_init(&rt);
1464     dlist = mk_dlist(log, &atoms, &nlist, bPhi, bPsi, bChi, bHChi, maxchi, r0, rt);
1465     fprintf(stderr, "%d residues with dihedrals found\n", nlist);
1466
1467     if (nlist == 0)
1468     {
1469         gmx_fatal(FARGS, "No dihedrals in your structure!\n");
1470     }
1471
1472     /* Make a linear index for reading all. */
1473     index = make_chi_ind(nlist, dlist, &ndih);
1474     isize = 4*ndih;
1475     fprintf(stderr, "%d dihedrals found\n", ndih);
1476
1477     snew(dih, ndih);
1478
1479     /* COMPUTE ALL DIHEDRALS! */
1480     read_ang_dih(ftp2fn(efTRX, NFILE, fnm), FALSE, TRUE, FALSE, bPBC, 1, &idum,
1481                  &nf, &time, isize, index, &trans_frac, &aver_angle, dih, oenv);
1482
1483     dt = (time[nf-1]-time[0])/(nf-1); /* might want this for corr or n. transit*/
1484     if (bCorr)
1485     {
1486         if (nf < 2)
1487         {
1488             gmx_fatal(FARGS, "Need at least 2 frames for correlation");
1489         }
1490     }
1491
1492     /* put angles in -M_PI to M_PI ! and correct phase factor for phi and psi
1493      * pass nactdih instead of ndih to low_ana_dih_trans and get_chi_product_traj
1494      * to prevent accessing off end of arrays when maxchi < 5 or 6. */
1495     nactdih = reset_em_all(nlist, dlist, nf, dih, maxchi);
1496
1497     if (bAll)
1498     {
1499         dump_em_all(nlist, dlist, nf, time, dih, maxchi, bPhi, bPsi, bChi, bOmega, bRAD, oenv);
1500     }
1501
1502     /* Histogramming & J coupling constants & calc of S2 order params */
1503     histogramming(log, nbin, rt, nf, maxchi, dih, nlist, dlist, index,
1504                   bPhi, bPsi, bOmega, bChi,
1505                   bNormHisto, bSSHisto, ftp2fn(efDAT, NFILE, fnm), bfac_max, &atoms,
1506                   bDo_jc, opt2fn("-jc", NFILE, fnm), oenv);
1507
1508     /* transitions
1509      *
1510      * added multiplicity */
1511
1512     snew(multiplicity, ndih);
1513     mk_multiplicity_lookup(multiplicity, maxchi, nlist, dlist, ndih);
1514
1515     std::strcpy(grpname, "All residues, ");
1516     if (bPhi)
1517     {
1518         std::strcat(grpname, "Phi ");
1519     }
1520     if (bPsi)
1521     {
1522         std::strcat(grpname, "Psi ");
1523     }
1524     if (bOmega)
1525     {
1526         std::strcat(grpname, "Omega ");
1527     }
1528     if (bChi)
1529     {
1530         std::strcat(grpname, "Chi 1-");
1531         sprintf(grpname + std::strlen(grpname), "%i", maxchi);
1532     }
1533
1534
1535     low_ana_dih_trans(bDo_ot, opt2fn("-ot", NFILE, fnm),
1536                       bDo_oh, opt2fn("-oh", NFILE, fnm), maxchi,
1537                       dih, nlist, dlist, nf, nactdih, grpname, multiplicity,
1538                       time, FALSE, core_frac, oenv);
1539
1540     /* Order parameters */
1541     order_params(log, opt2fn("-o", NFILE, fnm), maxchi, nlist, dlist,
1542                  ftp2fn_null(efPDB, NFILE, fnm), bfac_init,
1543                  &atoms, x, ePBC, box, bPhi, bPsi, bChi, oenv);
1544
1545     /* Print ramachandran maps! */
1546     if (bRama)
1547     {
1548         do_rama(nf, nlist, dlist, dih, bViol, bRamOmega, oenv);
1549     }
1550
1551     if (bShift)
1552     {
1553         do_pp2shifts(log, nf, nlist, dlist, dih);
1554     }
1555
1556     /* rprint S^2, transitions, and rotamer occupancies to log */
1557     traj_t_ns = 0.001 * (time[nf-1]-time[0]);
1558     pr_dlist(log, nlist, dlist, traj_t_ns, edPrintST, bPhi, bPsi, bChi, bOmega, maxchi);
1559     pr_dlist(log, nlist, dlist, traj_t_ns, edPrintRO, bPhi, bPsi, bChi, bOmega, maxchi);
1560     gmx_ffclose(log);
1561     /* transitions to xvg */
1562     if (bDo_rt)
1563     {
1564         print_transitions(opt2fn("-rt", NFILE, fnm), maxchi, nlist, dlist, traj_t_ns, oenv);
1565     }
1566
1567     /* chi_product trajectories (ie one "rotamer number" for each residue) */
1568     if (bChiProduct && bChi)
1569     {
1570         snew(chi_lookup, nlist);
1571         for (i = 0; i < nlist; i++)
1572         {
1573             snew(chi_lookup[i], maxchi);
1574         }
1575         mk_chi_lookup(chi_lookup, maxchi, nlist, dlist);
1576
1577         get_chi_product_traj(dih, nf, nactdih,
1578                              maxchi, dlist, time, chi_lookup, multiplicity,
1579                              FALSE, bNormHisto, core_frac, bAll,
1580                              opt2fn("-cp", NFILE, fnm), oenv);
1581
1582         for (i = 0; i < nlist; i++)
1583         {
1584             sfree(chi_lookup[i]);
1585         }
1586     }
1587
1588     /* Correlation comes last because it messes up the angles */
1589     if (bCorr)
1590     {
1591         do_dihcorr(opt2fn("-corr", NFILE, fnm), nf, ndih, dih, dt, nlist, dlist, time,
1592                    maxchi, bPhi, bPsi, bChi, bOmega, oenv);
1593     }
1594
1595
1596     do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1597     do_view(oenv, opt2fn("-jc", NFILE, fnm), "-nxy");
1598     if (bCorr)
1599     {
1600         do_view(oenv, opt2fn("-corr", NFILE, fnm), "-nxy");
1601     }
1602
1603     gmx_residuetype_destroy(rt);
1604
1605     return 0;
1606 }