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", dlist[i].resnr, bfac_max);
 600                     }
 601                 }
 602                 else
 603                 {
 604                     n += 4;
 605                 }
 606
 607                 switch (Dih)
 608                 {
 609                     case edPhi:
 610                         calc_distribution_props(nbin, histmp, -M_PI, NKKKPHI, kkkphi, &S2);
 611
 612                         for (m = 0; (m < NKKKPHI); m++)
 613                         {
 614                             Jc[i][m]    = kkkphi[m].Jc;
 615                             Jcsig[i][m] = kkkphi[m].Jcsig;
 616                         }
 617                         break;
 618                     case edPsi:
 619                         calc_distribution_props(nbin, histmp, -M_PI, NKKKPSI, kkkpsi, &S2);
 620
 621                         for (m = 0; (m < NKKKPSI); m++)
 622                         {
 623                             Jc[i][NKKKPHI + m]    = kkkpsi[m].Jc;
 624                             Jcsig[i][NKKKPHI + m] = kkkpsi[m].Jcsig;
 625                         }
 626                         break;
 627                     case edChi1:
 628                         calc_distribution_props(nbin, histmp, -M_PI, NKKKCHI, kkkchi1, &S2);
 629                         for (m = 0; (m < NKKKCHI); m++)
 630                         {
 631                             Jc[i][NKKKPHI + NKKKPSI + m]    = kkkchi1[m].Jc;
 632                             Jcsig[i][NKKKPHI + NKKKPSI + m] = kkkchi1[m].Jcsig;
 633                         }
 634                         break;
 635                     default: /* covers edOmega and higher Chis than Chi1 */
 636                         calc_distribution_props(nbin, histmp, -M_PI, 0, nullptr, &S2);
 637                         break;
 638                 }
 639                 dlist[i].S2[Dih] = S2;
 640
 641                 /* Sum distribution per amino acid type as well */
 642                 for (k = 0; (k < nbin); k++)
 643                 {
 644                     his_aa[Dih][dlist[i].index][k] += histmp[k];
 645                     histmp[k] = 0;
 646                 }
 647                 j++;
 648             }
 649             else /* dihed not defined */
 650             {
 651                 dlist[i].S2[Dih] = 0.0;
 652             }
 653         }
 654     }
 655     sfree(histmp);
 656
 657     /* Print out Jcouplings */
 658     fprintf(log, "\n *** J-Couplings from simulation (plus std. dev.) ***\n\n");
 659     fprintf(log, "Residue   ");
 660     for (i = 0; (i < NKKKPHI); i++)
 661     {
 662         fprintf(log, "%7s   SD", kkkphi[i].name);
 663     }
 664     for (i = 0; (i < NKKKPSI); i++)
 665     {
 666         fprintf(log, "%7s   SD", kkkpsi[i].name);
 667     }
 668     for (i = 0; (i < NKKKCHI); i++)
 669     {
 670         fprintf(log, "%7s   SD", kkkchi1[i].name);
 671     }
 672     fprintf(log, "\n");
 673     for (i = 0; (i < NJC + 1); i++)
 674     {
 675         fprintf(log, "------------");
 676     }
 677     fprintf(log, "\n");
 678     for (i = 0; (i < nlist); i++)
 679     {
 680         fprintf(log, "%-10s", dlist[i].name);
 681         for (j = 0; (j < NJC); j++)
 682         {
 683             fprintf(log, "  %5.2f %4.2f", Jc[i][j], Jcsig[i][j]);
 684         }
 685         fprintf(log, "\n");
 686     }
 687     fprintf(log, "\n");
 688
 689     /* and to -jc file... */
 690     if (bDo_jc)
 691     {
 692         fp = xvgropen(fn, "\\S3\\NJ-Couplings from Karplus Equation", "Residue", "Coupling", oenv);
 693         snew(leg, NJC);
 694         for (i = 0; (i < NKKKPHI); i++)
 695         {
 696             leg[i] = gmx_strdup(kkkphi[i].name);
 697         }
 698         for (i = 0; (i < NKKKPSI); i++)
 699         {
 700             leg[i + NKKKPHI] = gmx_strdup(kkkpsi[i].name);
 701         }
 702         for (i = 0; (i < NKKKCHI); i++)
 703         {
 704             leg[i + NKKKPHI + NKKKPSI] = gmx_strdup(kkkchi1[i].name);
 705         }
 706         xvgr_legend(fp, NJC, leg, oenv);
 707         fprintf(fp, "%5s ", "#Res.");
 708         for (i = 0; (i < NJC); i++)
 709         {
 710             fprintf(fp, "%10s ", leg[i]);
 711         }
 712         fprintf(fp, "\n");
 713         for (i = 0; (i < nlist); i++)
 714         {
 715             fprintf(fp, "%5d ", dlist[i].resnr);
 716             for (j = 0; (j < NJC); j++)
 717             {
 718                 fprintf(fp, "  %8.3f", Jc[i][j]);
 719             }
 720             fprintf(fp, "\n");
 721         }
 722         xvgrclose(fp);
 723         for (i = 0; (i < NJC); i++)
 724         {
 725             sfree(leg[i]);
 726         }
 727     }
 728     /* finished -jc stuff */
 729
 730     snew(normhisto, nbin);
 731     for (i = 0; (i < rt_size); i++)
 732     {
 733         for (Dih = 0; (Dih < edMax); Dih++)
 734         {
 735             /* First check whether something is in there */
 736             for (j = 0; (j < nbin); j++)
 737             {
 738                 if (his_aa[Dih][i][j] != 0)
 739                 {
 740                     break;
 741                 }
 742             }
 743             if ((j < nbin)
 744                 && ((bPhi && (Dih == edPhi)) || (bPsi && (Dih == edPsi))
 745                     || (bOmega && (Dih == edOmega)) || (bChi && (Dih >= edChi1))))
 746             {
 747                 if (bNormalize)
 748                 {
 749                     normalize_histo(nbin, his_aa[Dih][i], (360.0 / nbin), normhisto);
 750                 }
 751
 752                 residue_name = rt->nameFromResidueIndex(i).c_str();
 753                 switch (Dih)
 754                 {
 755                     case edPhi:
 756                         sprintf(hisfile, "histo-phi%s", residue_name);
 757                         sprintf(title, "\\xf\\f{} Distribution for %s", residue_name);
 758                         break;
 759                     case edPsi:
 760                         sprintf(hisfile, "histo-psi%s", residue_name);
 761                         sprintf(title, "\\xy\\f{} Distribution for %s", residue_name);
 762                         break;
 763                     case edOmega:
 764                         sprintf(hisfile, "histo-omega%s", residue_name);
 765                         sprintf(title, "\\xw\\f{} Distribution for %s", residue_name);
 766                         break;
 767                     default:
 768                         sprintf(hisfile, "histo-chi%d%s", Dih - NONCHI + 1, residue_name);
 769                         sprintf(title, "\\xc\\f{}\\s%d\\N Distribution for %s", Dih - NONCHI + 1, residue_name);
 770                 }
 771                 std::strcpy(hhisfile, hisfile);
 772                 std::strcat(hhisfile, ".xvg");
 773                 fp = xvgropen(hhisfile, title, "Degrees", "", oenv);
 774                 if (output_env_get_print_xvgr_codes(oenv))
 775                 {
 776                     fprintf(fp, "@ with g0\n");
 777                 }
 778                 xvgr_world(fp, -180, 0, 180, 0.1, oenv);
 779                 if (output_env_get_print_xvgr_codes(oenv))
 780                 {
 781                     fprintf(fp,
 782                             "# this effort to set graph size fails unless you run with -autoscale "
 783                             "none or -autoscale y flags\n");
 784                     fprintf(fp, "@ xaxis tick on\n");
 785                     fprintf(fp, "@ xaxis tick major 90\n");
 786                     fprintf(fp, "@ xaxis tick minor 30\n");
 787                     fprintf(fp, "@ xaxis ticklabel prec 0\n");
 788                     fprintf(fp, "@ yaxis tick off\n");
 789                     fprintf(fp, "@ yaxis ticklabel off\n");
 790                     fprintf(fp, "@ type xy\n");
 791                 }
 792                 if (bSSHisto)
 793                 {
 794                     for (k = 0; (k < 3); k++)
 795                     {
 796                         std::string sshisfile = gmx::formatString("%s-%s.xvg", hisfile, sss[k]);
 797                         ssfp[k]               = gmx_ffopen(sshisfile, "w");
 798                     }
 799                 }
 800                 for (j = 0; (j < nbin); j++)
 801                 {
 802                     angle = -180 + (360 / nbin) * j;
 803                     if (bNormalize)
 804                     {
 805                         fprintf(fp, "%5d  %10g\n", angle, normhisto[j]);
 806                     }
 807                     else
 808                     {
 809                         fprintf(fp, "%5d  %10d\n", angle, his_aa[Dih][i][j]);
 810                     }
 811                     if (bSSHisto)
 812                     {
 813                         for (k = 0; (k < 3); k++)
 814                         {
 815                             fprintf(ssfp[k], "%5d  %10d\n", angle, his_aa_ss[k][i][Dih][j]);
 816                         }
 817                     }
 818                 }
 819                 fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 820                 xvgrclose(fp);
 821                 if (bSSHisto)
 822                 {
 823                     for (k = 0; (k < 3); k++)
 824                     {
 825                         fprintf(ssfp[k], "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 826                         gmx_ffclose(ssfp[k]);
 827                     }
 828                 }
 829             }
 830         }
 831     }
 832     sfree(normhisto);
 833
 834     if (bSSHisto)
 835     {
 836         /* Four dimensional array... Very cool */
 837         for (i = 0; (i < 3); i++)
 838         {
 839             for (j = 0; (j <= rt_size); j++)
 840             {
 841                 for (Dih = 0; (Dih < edMax); Dih++)
 842                 {
 843                     sfree(his_aa_ss[i][j][Dih]);
 844                 }
 845                 sfree(his_aa_ss[i][j]);
 846             }
 847             sfree(his_aa_ss[i]);
 848         }
 849         sfree(his_aa_ss);
 850         sfree(ss_str);
 851     }
 852 }
 853
 854 static FILE* rama_file(const char* fn, const char* title, const char* xaxis, const char* yaxis, const gmx_output_env_t* oenv)
 855 {
 856     FILE* fp;
 857
 858     fp = xvgropen(fn, title, xaxis, yaxis, oenv);
 859     if (output_env_get_print_xvgr_codes(oenv))
 860     {
 861         fprintf(fp, "@ with g0\n");
 862     }
 863     xvgr_world(fp, -180, -180, 180, 180, oenv);
 864     if (output_env_get_print_xvgr_codes(oenv))
 865     {
 866         fprintf(fp, "@ xaxis tick on\n");
 867         fprintf(fp, "@ xaxis tick major 90\n");
 868         fprintf(fp, "@ xaxis tick minor 30\n");
 869         fprintf(fp, "@ xaxis ticklabel prec 0\n");
 870         fprintf(fp, "@ yaxis tick on\n");
 871         fprintf(fp, "@ yaxis tick major 90\n");
 872         fprintf(fp, "@ yaxis tick minor 30\n");
 873         fprintf(fp, "@ yaxis ticklabel prec 0\n");
 874         fprintf(fp, "@    s0 type xy\n");
 875         fprintf(fp, "@    s0 symbol 2\n");
 876         fprintf(fp, "@    s0 symbol size 0.410000\n");
 877         fprintf(fp, "@    s0 symbol fill 1\n");
 878         fprintf(fp, "@    s0 symbol color 1\n");
 879         fprintf(fp, "@    s0 symbol linewidth 1\n");
 880         fprintf(fp, "@    s0 symbol linestyle 1\n");
 881         fprintf(fp, "@    s0 symbol center false\n");
 882         fprintf(fp, "@    s0 symbol char 0\n");
 883         fprintf(fp, "@    s0 skip 0\n");
 884         fprintf(fp, "@    s0 linestyle 0\n");
 885         fprintf(fp, "@    s0 linewidth 1\n");
 886         fprintf(fp, "@ type xy\n");
 887     }
 888     return fp;
 889 }
 890
 891 static void do_rama(int                     nf,
 892                     int                     nlist,
 893                     t_dlist                 dlist[],
 894                     real**                  dih,
 895                     gmx_bool                bViol,
 896                     gmx_bool                bRamOmega,
 897                     const gmx_output_env_t* oenv)
 898 {
 899     FILE *        fp, *gp = nullptr;
 900     gmx_bool      bOm;
 901     char          fn[256];
 902     int           i, j, k, Xi1, Xi2, Phi, Psi, Om = 0, nlevels;
 903     constexpr int NMAT = 120;
 904     real **       mat  = nullptr, phi, psi, omega, axis[NMAT], lo, hi;
 905     t_rgb         rlo  = { 1.0, 0.0, 0.0 };
 906     t_rgb         rmid = { 1.0, 1.0, 1.0 };
 907     t_rgb         rhi  = { 0.0, 0.0, 1.0 };
 908
 909     for (i = 0; (i < nlist); i++)
 910     {
 911         if ((has_dihedral(edPhi, &(dlist[i]))) && (has_dihedral(edPsi, &(dlist[i]))))
 912         {
 913             sprintf(fn, "ramaPhiPsi%s.xvg", dlist[i].name);
 914             fp  = rama_file(fn, "Ramachandran Plot", "\\8f\\4 (deg)", "\\8y\\4 (deg)", oenv);
 915             bOm = bRamOmega && has_dihedral(edOmega, &(dlist[i]));
 916             if (bOm)
 917             {
 918                 Om = dlist[i].j0[edOmega];
 919                 snew(mat, NMAT);
 920                 for (j = 0; (j < NMAT); j++)
 921                 {
 922                     snew(mat[j], NMAT);
 923                     axis[j] = -180 + gmx::exactDiv(360 * j, NMAT);
 924                 }
 925             }
 926             if (bViol)
 927             {
 928                 sprintf(fn, "violPhiPsi%s.xvg", dlist[i].name);
 929                 gp = gmx_ffopen(fn, "w");
 930             }
 931             Phi = dlist[i].j0[edPhi];
 932             Psi = dlist[i].j0[edPsi];
 933             for (j = 0; (j < nf); j++)
 934             {
 935                 phi = RAD2DEG * dih[Phi][j];
 936                 psi = RAD2DEG * dih[Psi][j];
 937                 fprintf(fp, "%10g  %10g\n", phi, psi);
 938                 if (bViol)
 939                 {
 940                     fprintf(gp, "%d\n", static_cast<int>(!bAllowed(dih[Phi][j], RAD2DEG * dih[Psi][j])));
 941                 }
 942                 if (bOm)
 943                 {
 944                     omega = RAD2DEG * dih[Om][j];
 945                     mat[static_cast<int>(((phi * NMAT) / 360) + gmx::exactDiv(NMAT, 2))]
 946                        [static_cast<int>(((psi * NMAT) / 360) + gmx::exactDiv(NMAT, 2))] += omega;
 947                 }
 948             }
 949             if (bViol)
 950             {
 951                 gmx_ffclose(gp);
 952             }
 953             xvgrclose(fp);
 954             if (bOm)
 955             {
 956                 sprintf(fn, "ramomega%s.xpm", dlist[i].name);
 957                 fp = gmx_ffopen(fn, "w");
 958                 lo = hi = 0;
 959                 for (j = 0; (j < NMAT); j++)
 960                 {
 961                     for (k = 0; (k < NMAT); k++)
 962                     {
 963                         mat[j][k] /= nf;
 964                         lo = std::min(mat[j][k], lo);
 965                         hi = std::max(mat[j][k], hi);
 966                     }
 967                 }
 968                 /* Symmetrise */
 969                 if (std::abs(lo) > std::abs(hi))
 970                 {
 971                     hi = -lo;
 972                 }
 973                 else
 974                 {
 975                     lo = -hi;
 976                 }
 977                 /* Add 180 */
 978                 for (j = 0; (j < NMAT); j++)
 979                 {
 980                     for (k = 0; (k < NMAT); k++)
 981                     {
 982                         mat[j][k] += 180;
 983                     }
 984                 }
 985                 lo += 180;
 986                 hi += 180;
 987                 nlevels = 20;
 988                 write_xpm3(fp,
 989                            0,
 990                            "Omega/Ramachandran Plot",
 991                            "Deg",
 992                            "Phi",
 993                            "Psi",
 994                            NMAT,
 995                            NMAT,
 996                            axis,
 997                            axis,
 998                            mat,
 999                            lo,
1000                            180.0,
1001                            hi,
1002                            rlo,
1003                            rmid,
1004                            rhi,
1005                            &nlevels);
1006                 gmx_ffclose(fp);
1007                 for (j = 0; (j < NMAT); j++)
1008                 {
1009                     sfree(mat[j]);
1010                 }
1011                 sfree(mat);
1012             }
1013         }
1014         if ((has_dihedral(edChi1, &(dlist[i]))) && (has_dihedral(edChi2, &(dlist[i]))))
1015         {
1016             sprintf(fn, "ramaX1X2%s.xvg", dlist[i].name);
1017             fp  = rama_file(fn,
1018                            "\\8c\\4\\s1\\N-\\8c\\4\\s2\\N Ramachandran Plot",
1019                            "\\8c\\4\\s1\\N (deg)",
1020                            "\\8c\\4\\s2\\N (deg)",
1021                            oenv);
1022             Xi1 = dlist[i].j0[edChi1];
1023             Xi2 = dlist[i].j0[edChi2];
1024             for (j = 0; (j < nf); j++)
1025             {
1026                 fprintf(fp, "%10g  %10g\n", RAD2DEG * dih[Xi1][j], RAD2DEG * dih[Xi2][j]);
1027             }
1028             xvgrclose(fp);
1029         }
1030         else
1031         {
1032             fprintf(stderr, "No chi1 & chi2 angle for %s\n", dlist[i].name);
1033         }
1034     }
1035 }
1036
1037
1038 static void print_transitions(const char* fn, int maxchi, int nlist, t_dlist dlist[], real dt, const gmx_output_env_t* oenv)
1039 {
1040     /* based on order_params below */
1041     FILE* fp;
1042     int   i, Dih, Xi;
1043
1044     /*  must correspond with enum in pp2shift.h:38 */
1045     char* leg[edMax];
1046 #define NLEG asize(leg)
1047
1048     leg[0] = gmx_strdup("Phi");
1049     leg[1] = gmx_strdup("Psi");
1050     leg[2] = gmx_strdup("Omega");
1051     leg[3] = gmx_strdup("Chi1");
1052     leg[4] = gmx_strdup("Chi2");
1053     leg[5] = gmx_strdup("Chi3");
1054     leg[6] = gmx_strdup("Chi4");
1055     leg[7] = gmx_strdup("Chi5");
1056     leg[8] = gmx_strdup("Chi6");
1057
1058     /* Print order parameters */
1059     fp = xvgropen(fn, "Dihedral Rotamer Transitions", "Residue", "Transitions/ns", oenv);
1060     xvgr_legend(fp, NONCHI + maxchi, leg, oenv);
1061
1062     fprintf(fp, "%5s ", "#Res.");
1063     fprintf(fp, "%10s %10s %10s ", leg[edPhi], leg[edPsi], leg[edOmega]);
1064     for (Xi = 0; Xi < maxchi; Xi++)
1065     {
1066         fprintf(fp, "%10s ", leg[NONCHI + Xi]);
1067     }
1068     fprintf(fp, "\n");
1069
1070     for (i = 0; (i < nlist); i++)
1071     {
1072         fprintf(fp, "%5d ", dlist[i].resnr);
1073         for (Dih = 0; (Dih < NONCHI + maxchi); Dih++)
1074         {
1075             fprintf(fp, "%10.3f ", dlist[i].ntr[Dih] / dt);
1076         }
1077         /* fprintf(fp,"%12s\n",dlist[i].name);  this confuses xmgrace */
1078         fprintf(fp, "\n");
1079     }
1080     xvgrclose(fp);
1081 }
1082
1083 static void order_params(FILE*                   log,
1084                          const char*             fn,
1085                          int                     maxchi,
1086                          int                     nlist,
1087                          t_dlist                 dlist[],
1088                          const char*             pdbfn,
1089                          real                    bfac_init,
1090                          t_atoms*                atoms,
1091                          const rvec              x[],
1092                          PbcType                 pbcType,
1093                          matrix                  box,
1094                          gmx_bool                bPhi,
1095                          gmx_bool                bPsi,
1096                          gmx_bool                bChi,
1097                          const gmx_output_env_t* oenv)
1098 {
1099     FILE* fp;
1100     int   nh[edMax];
1101     int   i, Dih, Xi;
1102     real  S2Max, S2Min;
1103
1104     /* except for S2Min/Max, must correspond with enum in pp2shift.h:38 */
1105     const char* const_leg[2 + edMax] = { "S2Min", "S2Max", "Phi",  "Psi",  "Omega", "Chi1",
1106                                          "Chi2",  "Chi3",  "Chi4", "Chi5", "Chi6" };
1107 #define NLEG asize(leg)
1108
1109     char* leg[2 + edMax];
1110
1111     for (i = 0; i < NLEG; i++)
1112     {
1113         leg[i] = gmx_strdup(const_leg[i]);
1114     }
1115
1116     /* Print order parameters */
1117     fp = xvgropen(fn, "Dihedral Order Parameters", "Residue", "S2", oenv);
1118     xvgr_legend(fp, 2 + NONCHI + maxchi, const_leg, oenv);
1119
1120     for (Dih = 0; (Dih < edMax); Dih++)
1121     {
1122         nh[Dih] = 0;
1123     }
1124
1125     fprintf(fp, "%5s ", "#Res.");
1126     fprintf(fp, "%10s %10s ", leg[0], leg[1]);
1127     fprintf(fp, "%10s %10s %10s ", leg[2 + edPhi], leg[2 + edPsi], leg[2 + edOmega]);
1128     for (Xi = 0; Xi < maxchi; Xi++)
1129     {
1130         fprintf(fp, "%10s ", leg[2 + NONCHI + Xi]);
1131     }
1132     fprintf(fp, "\n");
1133
1134     for (i = 0; (i < nlist); i++)
1135     {
1136         S2Max = -10;
1137         S2Min = 10;
1138         for (Dih = 0; (Dih < NONCHI + maxchi); Dih++)
1139         {
1140             if (dlist[i].S2[Dih] != 0)
1141             {
1142                 if (dlist[i].S2[Dih] > S2Max)
1143                 {
1144                     S2Max = dlist[i].S2[Dih];
1145                 }
1146                 if (dlist[i].S2[Dih] < S2Min)
1147                 {
1148                     S2Min = dlist[i].S2[Dih];
1149                 }
1150             }
1151             if (dlist[i].S2[Dih] > 0.8)
1152             {
1153                 nh[Dih]++;
1154             }
1155         }
1156         fprintf(fp, "%5d ", dlist[i].resnr);
1157         fprintf(fp, "%10.3f %10.3f ", S2Min, S2Max);
1158         for (Dih = 0; (Dih < NONCHI + maxchi); Dih++)
1159         {
1160             fprintf(fp, "%10.3f ", dlist[i].S2[Dih]);
1161         }
1162         fprintf(fp, "\n");
1163         /* fprintf(fp,"%12s\n",dlist[i].name);  this confuses xmgrace */
1164     }
1165     xvgrclose(fp);
1166
1167     if (nullptr != pdbfn)
1168     {
1169         real x0, y0, z0;
1170
1171         atoms->havePdbInfo = TRUE;
1172
1173         if (nullptr == atoms->pdbinfo)
1174         {
1175             snew(atoms->pdbinfo, atoms->nr);
1176         }
1177         for (i = 0; (i < atoms->nr); i++)
1178         {
1179             atoms->pdbinfo[i].bfac = bfac_init;
1180         }
1181
1182         for (i = 0; (i < nlist); i++)
1183         {
1184             atoms->pdbinfo[dlist[i].atm.N].bfac = -dlist[i].S2[0]; /* Phi */
1185             atoms->pdbinfo[dlist[i].atm.H].bfac = -dlist[i].S2[0]; /* Phi */
1186             atoms->pdbinfo[dlist[i].atm.C].bfac = -dlist[i].S2[1]; /* Psi */
1187             atoms->pdbinfo[dlist[i].atm.O].bfac = -dlist[i].S2[1]; /* Psi */
1188             for (Xi = 0; (Xi < maxchi); Xi++)                      /* Chi's */
1189             {
1190                 if (dlist[i].atm.Cn[Xi + 3] != -1)
1191                 {
1192                     atoms->pdbinfo[dlist[i].atm.Cn[Xi + 1]].bfac = -dlist[i].S2[NONCHI + Xi];
1193                 }
1194             }
1195         }
1196
1197         fp = gmx_ffopen(pdbfn, "w");
1198         fprintf(fp, "REMARK generated by g_chi\n");
1199         fprintf(fp,
1200                 "REMARK "
1201                 "B-factor field contains negative of dihedral order parameters\n");
1202         write_pdbfile(fp, nullptr, atoms, x, pbcType, box, ' ', 0, nullptr);
1203         x0 = y0 = z0 = 1000.0;
1204         for (i = 0; (i < atoms->nr); i++)
1205         {
1206             x0 = std::min(x0, x[i][XX]);
1207             y0 = std::min(y0, x[i][YY]);
1208             z0 = std::min(z0, x[i][ZZ]);
1209         }
1210         x0 *= 10.0; /* nm -> angstrom */
1211         y0 *= 10.0; /* nm -> angstrom */
1212         z0 *= 10.0; /* nm -> angstrom */
1213         for (i = 0; (i < 10); i++)
1214         {
1215             gmx_fprintf_pdb_atomline(fp,
1216                                      epdbATOM,
1217                                      atoms->nr + 1 + i,
1218                                      "CA",
1219                                      ' ',
1220                                      "LEG",
1221                                      ' ',
1222                                      atoms->nres + 1,
1223                                      ' ',
1224                                      x0,
1225                                      y0,
1226                                      z0 + (1.2 * i),
1227                                      0.0,
1228                                      -0.1 * i,
1229                                      "");
1230         }
1231         gmx_ffclose(fp);
1232     }
1233
1234     fprintf(log, "Dihedrals with S2 > 0.8\n");
1235     fprintf(log, "Dihedral: ");
1236     if (bPhi)
1237     {
1238         fprintf(log, " Phi  ");
1239     }
1240     if (bPsi)
1241     {
1242         fprintf(log, " Psi ");
1243     }
1244     if (bChi)
1245     {
1246         for (Xi = 0; (Xi < maxchi); Xi++)
1247         {
1248             fprintf(log, " %s ", leg[2 + NONCHI + Xi]);
1249         }
1250     }
1251     fprintf(log, "\nNumber:   ");
1252     if (bPhi)
1253     {
1254         fprintf(log, "%4d  ", nh[0]);
1255     }
1256     if (bPsi)
1257     {
1258         fprintf(log, "%4d  ", nh[1]);
1259     }
1260     if (bChi)
1261     {
1262         for (Xi = 0; (Xi < maxchi); Xi++)
1263         {
1264             fprintf(log, "%4d  ", nh[NONCHI + Xi]);
1265         }
1266     }
1267     fprintf(log, "\n");
1268
1269     for (i = 0; i < NLEG; i++)
1270     {
1271         sfree(leg[i]);
1272     }
1273 }
1274
1275 int gmx_chi(int argc, char* argv[])
1276 {
1277     const char* desc[] = {
1278         "[THISMODULE] computes [GRK]phi[grk], [GRK]psi[grk], [GRK]omega[grk],",
1279         "and [GRK]chi[grk] dihedrals for all your",
1280         "amino acid backbone and sidechains.",
1281         "It can compute dihedral angle as a function of time, and as",
1282         "histogram distributions.",
1283         "The distributions [TT](histo-(dihedral)(RESIDUE).xvg[tt]) are cumulative over all ",
1284         "residues of each type.[PAR]",
1285         "If option [TT]-corr[tt] is given, the program will",
1286         "calculate dihedral autocorrelation functions. The function used",
1287         "is C(t) = [CHEVRON][COS][GRK]chi[grk]([GRK]tau[grk])[cos] ",
1288         "[COS][GRK]chi[grk]([GRK]tau[grk]+t)[cos][chevron]. The use of cosines",
1289         "rather than angles themselves, resolves the problem of periodicity.",
1290         "(Van der Spoel & Berendsen (1997), Biophys. J. 72, 2032-2041).",
1291         "Separate files for each dihedral of each residue",
1292         "[TT](corr(dihedral)(RESIDUE)(nresnr).xvg[tt]) are output, as well as a",
1293         "file containing the information for all residues (argument of [TT]-corr[tt]).[PAR]",
1294         "With option [TT]-all[tt], the angles themselves as a function of time for",
1295         "each residue are printed to separate files [TT](dihedral)(RESIDUE)(nresnr).xvg[tt].",
1296         "These can be in radians or degrees.[PAR]",
1297         "A log file (argument [TT]-g[tt]) is also written. This contains",
1298         "",
1299         " * information about the number of residues of each type.",
1300         " * The NMR ^3J coupling constants from the Karplus equation.",
1301         " * a table for each residue of the number of transitions between ",
1302         "   rotamers per nanosecond,  and the order parameter S^2 of each dihedral.",
1303         " * a table for each residue of the rotamer occupancy.",
1304         "",
1305         "All rotamers are taken as 3-fold, except for [GRK]omega[grk] and [GRK]chi[grk] dihedrals",
1306         "to planar groups (i.e. [GRK]chi[grk][SUB]2[sub] of aromatics, Asp and Asn; ",
1307         "[GRK]chi[grk][SUB]3[sub] of Glu",
1308         "and Gln; and [GRK]chi[grk][SUB]4[sub] of Arg), which are 2-fold. \"rotamer 0\" means ",
1309         "that the dihedral was not in the core region of each rotamer. ",
1310         "The width of the core region can be set with [TT]-core_rotamer[tt][PAR]",
1311
1312         "The S^2 order parameters are also output to an [REF].xvg[ref] file",
1313         "(argument [TT]-o[tt] ) and optionally as a [REF].pdb[ref] file with",
1314         "the S^2 values as B-factor (argument [TT]-p[tt]). ",
1315         "The total number of rotamer transitions per timestep",
1316         "(argument [TT]-ot[tt]), the number of transitions per rotamer",
1317         "(argument [TT]-rt[tt]), and the ^3J couplings (argument [TT]-jc[tt]), ",
1318         "can also be written to [REF].xvg[ref] files. Note that the analysis",
1319         "of rotamer transitions assumes that the supplied trajectory frames",
1320         "are equally spaced in time.[PAR]",
1321
1322         "If [TT]-chi_prod[tt] is set (and [TT]-maxchi[tt] > 0), cumulative rotamers, e.g.",
1323         "1+9([GRK]chi[grk][SUB]1[sub]-1)+3([GRK]chi[grk][SUB]2[sub]-1)+",
1324         "([GRK]chi[grk][SUB]3[sub]-1) (if the residue has three 3-fold ",
1325         "dihedrals and [TT]-maxchi[tt] >= 3)",
1326         "are calculated. As before, if any dihedral is not in the core region,",
1327         "the rotamer is taken to be 0. The occupancies of these cumulative ",
1328         "rotamers (starting with rotamer 0) are written to the file",
1329         "that is the argument of [TT]-cp[tt], and if the [TT]-all[tt] flag",
1330         "is given, the rotamers as functions of time",
1331         "are written to [TT]chiproduct(RESIDUE)(nresnr).xvg[tt] ",
1332         "and their occupancies to [TT]histo-chiproduct(RESIDUE)(nresnr).xvg[tt].[PAR]",
1333
1334         "The option [TT]-r[tt] generates a contour plot of the average [GRK]omega[grk] angle",
1335         "as a function of the [GRK]phi[grk] and [GRK]psi[grk] angles, that is, in a Ramachandran ",
1336         "plot the average [GRK]omega[grk] angle is plotted using color coding.",
1337
1338     };
1339
1340     const char* bugs[] = {
1341         "Produces MANY output files (up to about 4 times the number of residues in the "
1342         "protein, twice that if autocorrelation functions are calculated). Typically "
1343         "several hundred files are output.",
1344         "[GRK]phi[grk] and [GRK]psi[grk] dihedrals are calculated in a "
1345         "non-standard way, using H-N-CA-C for [GRK]phi[grk] instead of "
1346         "C(-)-N-CA-C, and N-CA-C-O for [GRK]psi[grk] instead of N-CA-C-N(+). "
1347         "This causes (usually small) discrepancies with the output of other "
1348         "tools like [gmx-rama].",
1349         "[TT]-r0[tt] option does not work properly",
1350         "Rotamers with multiplicity 2 are printed in [TT]chi.log[tt] as if they had ",
1351         "multiplicity 3, with the 3rd (g(+)) always having probability 0"
1352     };
1353
1354     /* defaults */
1355     static int         r0 = 1, ndeg = 1, maxchi = 2;
1356     static gmx_bool    bAll = FALSE;
1357     static gmx_bool    bPhi = FALSE, bPsi = FALSE, bOmega = FALSE;
1358     static real        bfac_init = -1.0, bfac_max = 0;
1359     static const char* maxchistr[] = { nullptr, "0", "1", "2", "3", "4", "5", "6", nullptr };
1360     static gmx_bool    bRama = FALSE, bShift = FALSE, bViol = FALSE, bRamOmega = FALSE;
1361     static gmx_bool bNormHisto = TRUE, bChiProduct = FALSE, bHChi = FALSE, bRAD = FALSE, bPBC = TRUE;
1362     static real     core_frac = 0.5;
1363     t_pargs         pa[]      = {
1364         { "-r0", FALSE, etINT, { &r0 }, "starting residue" },
1365         { "-phi", FALSE, etBOOL, { &bPhi }, "Output for [GRK]phi[grk] dihedral angles" },
1366         { "-psi", FALSE, etBOOL, { &bPsi }, "Output for [GRK]psi[grk] dihedral angles" },
1367         { "-omega",
1368           FALSE,
1369           etBOOL,
1370           { &bOmega },
1371           "Output for [GRK]omega[grk] dihedrals (peptide bonds)" },
1372         { "-rama",
1373           FALSE,
1374           etBOOL,
1375           { &bRama },
1376           "Generate [GRK]phi[grk]/[GRK]psi[grk] and "
1377           "[GRK]chi[grk][SUB]1[sub]/[GRK]chi[grk][SUB]2[sub] Ramachandran plots" },
1378         { "-viol",
1379           FALSE,
1380           etBOOL,
1381           { &bViol },
1382           "Write a file that gives 0 or 1 for violated Ramachandran angles" },
1383         { "-periodic", FALSE, etBOOL, { &bPBC }, "Print dihedral angles modulo 360 degrees" },
1384         { "-all", FALSE, etBOOL, { &bAll }, "Output separate files for every dihedral." },
1385         { "-rad",
1386           FALSE,
1387           etBOOL,
1388           { &bRAD },
1389           "in angle vs time files, use radians rather than degrees." },
1390         { "-shift",
1391           FALSE,
1392           etBOOL,
1393           { &bShift },
1394           "Compute chemical shifts from [GRK]phi[grk]/[GRK]psi[grk] angles" },
1395         { "-binwidth", FALSE, etINT, { &ndeg }, "bin width for histograms (degrees)" },
1396         { "-core_rotamer",
1397           FALSE,
1398           etREAL,
1399           { &core_frac },
1400           "only the central [TT]-core_rotamer[tt]\\*(360/multiplicity) belongs to each rotamer "
1401           "(the rest is assigned to rotamer 0)" },
1402         { "-maxchi", FALSE, etENUM, { maxchistr }, "calculate first ndih [GRK]chi[grk] dihedrals" },
1403         { "-normhisto", FALSE, etBOOL, { &bNormHisto }, "Normalize histograms" },
1404         { "-ramomega",
1405           FALSE,
1406           etBOOL,
1407           { &bRamOmega },
1408           "compute average omega as a function of [GRK]phi[grk]/[GRK]psi[grk] and plot it in an "
1409           "[REF].xpm[ref] plot" },
1410         { "-bfact",
1411           FALSE,
1412           etREAL,
1413           { &bfac_init },
1414           "B-factor value for [REF].pdb[ref] file for atoms with no calculated dihedral order "
1415           "parameter" },
1416         { "-chi_prod",
1417           FALSE,
1418           etBOOL,
1419           { &bChiProduct },
1420           "compute a single cumulative rotamer for each residue" },
1421         { "-HChi", FALSE, etBOOL, { &bHChi }, "Include dihedrals to sidechain hydrogens" },
1422         { "-bmax",
1423           FALSE,
1424           etREAL,
1425           { &bfac_max },
1426           "Maximum B-factor on any of the atoms that make up a dihedral, for the dihedral angle to "
1427           "be considere in the statistics. Applies to database work where a number of X-Ray "
1428           "structures is analyzed. [TT]-bmax[tt] <= 0 means no limit." }
1429     };
1430
1431     FILE*             log;
1432     int               nlist, idum, nbin;
1433     rvec*             x;
1434     PbcType           pbcType;
1435     matrix            box;
1436     char              grpname[256];
1437     t_dlist*          dlist;
1438     gmx_bool          bChi, bCorr, bSSHisto;
1439     gmx_bool          bDo_rt, bDo_oh, bDo_ot, bDo_jc;
1440     real              dt = 0, traj_t_ns;
1441     gmx_output_env_t* oenv;
1442
1443     int    isize, *index;
1444     int    ndih, nactdih, nf;
1445     real **dih, *trans_frac, *aver_angle, *time;
1446     int    i, **chi_lookup, *multiplicity;
1447
1448     t_filenm fnm[] = { { efSTX, "-s", nullptr, ffREAD },
1449                        { efTRX, "-f", nullptr, ffREAD },
1450                        { efXVG, "-o", "order", ffWRITE },
1451                        { efPDB, "-p", "order", ffOPTWR },
1452                        { efDAT, "-ss", "ssdump", ffOPTRD },
1453                        { efXVG, "-jc", "Jcoupling", ffWRITE },
1454                        { efXVG, "-corr", "dihcorr", ffOPTWR },
1455                        { efLOG, "-g", "chi", ffWRITE },
1456                        /* add two more arguments copying from g_angle */
1457                        { efXVG, "-ot", "dihtrans", ffOPTWR },
1458                        { efXVG, "-oh", "trhisto", ffOPTWR },
1459                        { efXVG, "-rt", "restrans", ffOPTWR },
1460                        { efXVG, "-cp", "chiprodhisto", ffOPTWR } };
1461 #define NFILE asize(fnm)
1462     int      npargs;
1463     t_pargs* ppa;
1464
1465     npargs = asize(pa);
1466     ppa    = add_acf_pargs(&npargs, pa);
1467     if (!parse_common_args(
1468                 &argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs, &oenv))
1469     {
1470         sfree(ppa);
1471         return 0;
1472     }
1473
1474     /* Handle result from enumerated type */
1475     sscanf(maxchistr[0], "%d", &maxchi);
1476     bChi = (maxchi > 0);
1477
1478     log = gmx_ffopen(ftp2fn(efLOG, NFILE, fnm), "w");
1479
1480     if (bRamOmega)
1481     {
1482         bOmega = TRUE;
1483         bPhi   = TRUE;
1484         bPsi   = TRUE;
1485     }
1486
1487     /* set some options */
1488     bDo_rt = (opt2bSet("-rt", NFILE, fnm));
1489     bDo_oh = (opt2bSet("-oh", NFILE, fnm));
1490     bDo_ot = (opt2bSet("-ot", NFILE, fnm));
1491     bDo_jc = (opt2bSet("-jc", NFILE, fnm));
1492     bCorr  = (opt2bSet("-corr", NFILE, fnm));
1493     if (bCorr)
1494     {
1495         fprintf(stderr, "Will calculate autocorrelation\n");
1496     }
1497
1498     if (core_frac > 1.0)
1499     {
1500         fprintf(stderr, "core_rotamer fraction > 1.0 ; will use 1.0\n");
1501         core_frac = 1.0;
1502     }
1503     if (core_frac < 0.0)
1504     {
1505         fprintf(stderr, "core_rotamer fraction < 0.0 ; will use 0.0\n");
1506         core_frac = 0.0;
1507     }
1508
1509     if (maxchi > MAXCHI)
1510     {
1511         fprintf(stderr, "Will only calculate first %d Chi dihedrals instead of %d.\n", MAXCHI, maxchi);
1512         maxchi = MAXCHI;
1513     }
1514     bSSHisto = ftp2bSet(efDAT, NFILE, fnm);
1515     nbin     = 360 / ndeg;
1516
1517     /* Find the chi angles using atoms struct and a list of amino acids */
1518     t_topology* top;
1519     snew(top, 1);
1520     read_tps_conf(ftp2fn(efSTX, NFILE, fnm), top, &pbcType, &x, nullptr, box, FALSE);
1521     t_atoms& atoms = top->atoms;
1522     if (atoms.pdbinfo == nullptr)
1523     {
1524         snew(atoms.pdbinfo, atoms.nr);
1525     }
1526     fprintf(log, "Title: %s\n", *top->name);
1527
1528     ResidueType rt;
1529     dlist = mk_dlist(log, &atoms, &nlist, bPhi, bPsi, bChi, bHChi, maxchi, r0, &rt);
1530     fprintf(stderr, "%d residues with dihedrals found\n", nlist);
1531
1532     if (nlist == 0)
1533     {
1534         gmx_fatal(FARGS, "No dihedrals in your structure!\n");
1535     }
1536
1537     /* Make a linear index for reading all. */
1538     index = make_chi_ind(nlist, dlist, &ndih);
1539     isize = 4 * ndih;
1540     fprintf(stderr, "%d dihedrals found\n", ndih);
1541
1542     snew(dih, ndih);
1543
1544     /* COMPUTE ALL DIHEDRALS! */
1545     read_ang_dih(
1546             ftp2fn(efTRX, NFILE, fnm), FALSE, TRUE, FALSE, bPBC, 1, &idum, &nf, &time, isize, index, &trans_frac, &aver_angle, dih, oenv);
1547
1548     dt = (time[nf - 1] - time[0]) / (nf - 1); /* might want this for corr or n. transit*/
1549     if (bCorr)
1550     {
1551         if (nf < 2)
1552         {
1553             gmx_fatal(FARGS, "Need at least 2 frames for correlation");
1554         }
1555     }
1556
1557     /* put angles in -M_PI to M_PI ! and correct phase factor for phi and psi
1558      * pass nactdih instead of ndih to low_ana_dih_trans and get_chi_product_traj
1559      * to prevent accessing off end of arrays when maxchi < 5 or 6. */
1560     nactdih = reset_em_all(nlist, dlist, nf, dih, maxchi);
1561
1562     if (bAll)
1563     {
1564         dump_em_all(nlist, dlist, nf, time, dih, maxchi, bPhi, bPsi, bChi, bOmega, bRAD, oenv);
1565     }
1566
1567     /* Histogramming & J coupling constants & calc of S2 order params */
1568     histogramming(log,
1569                   nbin,
1570                   &rt,
1571                   nf,
1572                   maxchi,
1573                   dih,
1574                   nlist,
1575                   dlist,
1576                   index,
1577                   bPhi,
1578                   bPsi,
1579                   bOmega,
1580                   bChi,
1581                   bNormHisto,
1582                   bSSHisto,
1583                   ftp2fn(efDAT, NFILE, fnm),
1584                   bfac_max,
1585                   &atoms,
1586                   bDo_jc,
1587                   opt2fn("-jc", NFILE, fnm),
1588                   oenv);
1589
1590     /* transitions
1591      *
1592      * added multiplicity */
1593
1594     snew(multiplicity, ndih);
1595     mk_multiplicity_lookup(multiplicity, maxchi, nlist, dlist, ndih);
1596
1597     std::strcpy(grpname, "All residues, ");
1598     if (bPhi)
1599     {
1600         std::strcat(grpname, "Phi ");
1601     }
1602     if (bPsi)
1603     {
1604         std::strcat(grpname, "Psi ");
1605     }
1606     if (bOmega)
1607     {
1608         std::strcat(grpname, "Omega ");
1609     }
1610     if (bChi)
1611     {
1612         std::strcat(grpname, "Chi 1-");
1613         sprintf(grpname + std::strlen(grpname), "%i", maxchi);
1614     }
1615
1616
1617     low_ana_dih_trans(bDo_ot,
1618                       opt2fn("-ot", NFILE, fnm),
1619                       bDo_oh,
1620                       opt2fn("-oh", NFILE, fnm),
1621                       maxchi,
1622                       dih,
1623                       nlist,
1624                       dlist,
1625                       nf,
1626                       nactdih,
1627                       grpname,
1628                       multiplicity,
1629                       time,
1630                       FALSE,
1631                       core_frac,
1632                       oenv);
1633
1634     /* Order parameters */
1635     order_params(log,
1636                  opt2fn("-o", NFILE, fnm),
1637                  maxchi,
1638                  nlist,
1639                  dlist,
1640                  ftp2fn_null(efPDB, NFILE, fnm),
1641                  bfac_init,
1642                  &atoms,
1643                  x,
1644                  pbcType,
1645                  box,
1646                  bPhi,
1647                  bPsi,
1648                  bChi,
1649                  oenv);
1650
1651     /* Print ramachandran maps! */
1652     if (bRama)
1653     {
1654         do_rama(nf, nlist, dlist, dih, bViol, bRamOmega, oenv);
1655     }
1656
1657     if (bShift)
1658     {
1659         do_pp2shifts(log, nf, nlist, dlist, dih);
1660     }
1661
1662     /* rprint S^2, transitions, and rotamer occupancies to log */
1663     traj_t_ns = 0.001 * (time[nf - 1] - time[0]);
1664     pr_dlist(log, nlist, dlist, traj_t_ns, edPrintST, bPhi, bPsi, bChi, bOmega, maxchi);
1665     pr_dlist(log, nlist, dlist, traj_t_ns, edPrintRO, bPhi, bPsi, bChi, bOmega, maxchi);
1666     gmx_ffclose(log);
1667     /* transitions to xvg */
1668     if (bDo_rt)
1669     {
1670         print_transitions(opt2fn("-rt", NFILE, fnm), maxchi, nlist, dlist, traj_t_ns, oenv);
1671     }
1672
1673     /* chi_product trajectories (ie one "rotamer number" for each residue) */
1674     if (bChiProduct && bChi)
1675     {
1676         snew(chi_lookup, nlist);
1677         for (i = 0; i < nlist; i++)
1678         {
1679             snew(chi_lookup[i], maxchi);
1680         }
1681         mk_chi_lookup(chi_lookup, maxchi, nlist, dlist);
1682
1683         get_chi_product_traj(dih,
1684                              nf,
1685                              nactdih,
1686                              maxchi,
1687                              dlist,
1688                              time,
1689                              chi_lookup,
1690                              multiplicity,
1691                              FALSE,
1692                              bNormHisto,
1693                              core_frac,
1694                              bAll,
1695                              opt2fn("-cp", NFILE, fnm),
1696                              oenv);
1697
1698         for (i = 0; i < nlist; i++)
1699         {
1700             sfree(chi_lookup[i]);
1701         }
1702     }
1703
1704     /* Correlation comes last because it messes up the angles */
1705     if (bCorr)
1706     {
1707         do_dihcorr(opt2fn("-corr", NFILE, fnm), nf, ndih, dih, dt, nlist, dlist, time, maxchi, bPhi, bPsi, bChi, bOmega, oenv);
1708     }
1709
1710
1711     do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1712     do_view(oenv, opt2fn("-jc", NFILE, fnm), "-nxy");
1713     if (bCorr)
1714     {
1715         do_view(oenv, opt2fn("-corr", NFILE, fnm), "-nxy");
1716     }
1717
1718     return 0;
1719 }