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