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