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