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