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
  46 #include "gromacs/commandline/pargs.h"
  47 #include "gromacs/commandline/viewit.h"
  48 #include "gromacs/correlationfunctions/autocorr.h"
  49 #include "gromacs/fileio/confio.h"
  50 #include "gromacs/fileio/matio.h"
  51 #include "gromacs/fileio/pdbio.h"
  52 #include "gromacs/fileio/xvgr.h"
  53 #include "gromacs/gmxana/gmx_ana.h"
  54 #include "gromacs/gmxana/gstat.h"
  55 #include "gromacs/math/units.h"
  56 #include "gromacs/math/utilities.h"
  57 #include "gromacs/math/vec.h"
  58 #include "gromacs/topology/residuetypes.h"
  59 #include "gromacs/topology/topology.h"
  60 #include "gromacs/utility/arraysize.h"
  61 #include "gromacs/utility/cstringutil.h"
  62 #include "gromacs/utility/fatalerror.h"
  63 #include "gromacs/utility/futil.h"
  64 #include "gromacs/utility/smalloc.h"
  65 #include "gromacs/utility/stringutil.h"
  66
  67 static gmx_bool bAllowed(real phi, real psi)
  68 {
  69     static const char* map[] = { "1100000000000000001111111000000000001111111111111111111111111",
  70                                  "1100000000000000001111110000000000011111111111111111111111111",
  71                                  "1100000000000000001111110000000000011111111111111111111111111",
  72                                  "1100000000000000001111100000000000111111111111111111111111111",
  73                                  "1100000000000000001111100000000000111111111111111111111111111",
  74                                  "1100000000000000001111100000000001111111111111111111111111111",
  75                                  "1100000000000000001111100000000001111111111111111111111111111",
  76                                  "1100000000000000001111100000000011111111111111111111111111111",
  77                                  "1110000000000000001111110000000111111111111111111111111111111",
  78                                  "1110000000000000001111110000001111111111111111111111111111111",
  79                                  "1110000000000000001111111000011111111111111111111111111111111",
  80                                  "1110000000000000001111111100111111111111111111111111111111111",
  81                                  "1110000000000000001111111111111111111111111111111111111111111",
  82                                  "1110000000000000001111111111111111111111111111111111111111111",
  83                                  "1110000000000000001111111111111111111111111111111111111111111",
  84                                  "1110000000000000001111111111111111111111111111111111111111111",
  85                                  "1110000000000000001111111111111110011111111111111111111111111",
  86                                  "1110000000000000001111111111111100000111111111111111111111111",
  87                                  "1110000000000000001111111111111000000000001111111111111111111",
  88                                  "1100000000000000001111111111110000000000000011111111111111111",
  89                                  "1100000000000000001111111111100000000000000011111111111111111",
  90                                  "1000000000000000001111111111000000000000000001111111111111110",
  91                                  "0000000000000000001111111110000000000000000000111111111111100",
  92                                  "0000000000000000000000000000000000000000000000000000000000000",
  93                                  "0000000000000000000000000000000000000000000000000000000000000",
  94                                  "0000000000000000000000000000000000000000000000000000000000000",
  95                                  "0000000000000000000000000000000000000000000000000000000000000",
  96                                  "0000000000000000000000000000000000000000000000000000000000000",
  97                                  "0000000000000000000000000000000000000000000000000000000000000",
  98                                  "0000000000000000000000000000000000000000000000000000000000000",
  99                                  "0000000000000000000000000000000000000000000000000000000000000",
 100                                  "0000000000000000000000000000000000000000000000000000000000000",
 101                                  "0000000000000000000000000000000000000000000000000000000000000",
 102                                  "0000000000000000000000000000000000000000000000000000000000000",
 103                                  "0000000000000000000000000000000000000000000000000000000000000",
 104                                  "0000000000000000000000000000000000000000000000000000000000000",
 105                                  "0000000000000000000000000000000000000000000000000000000000000",
 106                                  "0000000000000000000000000000000000000000000000000000000000000",
 107                                  "0000000000000000000000000000000000111111111111000000000000000",
 108                                  "1100000000000000000000000000000001111111111111100000000000111",
 109                                  "1100000000000000000000000000000001111111111111110000000000111",
 110                                  "0000000000000000000000000000000000000000000000000000000000000",
 111                                  "0000000000000000000000000000000000000000000000000000000000000",
 112                                  "0000000000000000000000000000000000000000000000000000000000000",
 113                                  "0000000000000000000000000000000000000000000000000000000000000",
 114                                  "0000000000000000000000000000000000000000000000000000000000000",
 115                                  "0000000000000000000000000000000000000000000000000000000000000",
 116                                  "0000000000000000000000000000000000000000000000000000000000000",
 117                                  "0000000000000000000000000000000000000000000000000000000000000",
 118                                  "0000000000000000000000000000000000000000000000000000000000000",
 119                                  "0000000000000000000000000000000000000000000000000000000000000",
 120                                  "0000000000000000000000000000000000000000000000000000000000000",
 121                                  "0000000000000000000000000000000000000000000000000000000000000",
 122                                  "0000000000000000000000000000000000000000000000000000000000000",
 123                                  "0000000000000000000000000000000000000000000000000000000000000",
 124                                  "0000000000000000000000000000000000000000000000000000000000000",
 125                                  "0000000000000000000000000000000000000000000000000000000000000",
 126                                  "0000000000000000000000000000000000000000000000000000000000000",
 127                                  "0000000000000000000000000000000000000000000000000000000000000",
 128                                  "0000000000000000000000000000000000000000000000000000000000000",
 129                                  "0000000000000000000000000000000000000000000000000000000000000" };
 130     int                x, y;
 131
 132 #define INDEX(ppp) (((static_cast<int>(360 + (ppp)*gmx::c_rad2Deg)) % 360) / 6)
 133     x = INDEX(phi);
 134     y = INDEX(psi);
 135 #undef INDEX
 136
 137     return map[x][y] == '1';
 138 }
 139
 140 static int* make_chi_ind(int nl, t_dlist dl[], int* ndih)
 141 {
 142     int* id;
 143     int  i, Xi, n;
 144
 145     /* There are nl residues with max edMax dihedrals with 4 atoms each */
 146     snew(id, nl * edMax * 4);
 147
 148     n = 0;
 149     for (i = 0; (i < nl); i++)
 150     {
 151         /* Phi, fake the first one */
 152         dl[i].j0[edPhi] = n / 4;
 153         if (dl[i].atm.minC >= 0)
 154         {
 155             id[n++] = dl[i].atm.minC;
 156         }
 157         else
 158         {
 159             id[n++] = dl[i].atm.H;
 160         }
 161         id[n++] = dl[i].atm.N;
 162         id[n++] = dl[i].atm.Cn[1];
 163         id[n++] = dl[i].atm.C;
 164     }
 165     for (i = 0; (i < nl); i++)
 166     {
 167         /* Psi, fake the last one */
 168         dl[i].j0[edPsi] = n / 4;
 169         id[n++]         = dl[i].atm.N;
 170         id[n++]         = dl[i].atm.Cn[1];
 171         id[n++]         = dl[i].atm.C;
 172         if (i < (nl - 1))
 173         {
 174             id[n++] = dl[i + 1].atm.N;
 175         }
 176         else
 177         {
 178             id[n++] = dl[i].atm.O;
 179         }
 180     }
 181     for (i = 0; (i < nl); i++)
 182     {
 183         /* Omega */
 184         if (has_dihedral(edOmega, &(dl[i])))
 185         {
 186             dl[i].j0[edOmega] = n / 4;
 187             id[n++]           = dl[i].atm.minCalpha;
 188             id[n++]           = dl[i].atm.minC;
 189             id[n++]           = dl[i].atm.N;
 190             id[n++]           = dl[i].atm.Cn[1];
 191         }
 192     }
 193     for (Xi = 0; (Xi < MAXCHI); Xi++)
 194     {
 195         /* Chi# */
 196         for (i = 0; (i < nl); i++)
 197         {
 198             if (dl[i].atm.Cn[Xi + 3] != -1)
 199             {
 200                 dl[i].j0[edChi1 + Xi] = n / 4;
 201                 id[n++]               = dl[i].atm.Cn[Xi];
 202                 id[n++]               = dl[i].atm.Cn[Xi + 1];
 203                 id[n++]               = dl[i].atm.Cn[Xi + 2];
 204                 id[n++]               = dl[i].atm.Cn[Xi + 3];
 205             }
 206         }
 207     }
 208     *ndih = n / 4;
 209
 210     return id;
 211 }
 212
 213 static void do_dihcorr(const char*             fn,
 214                        int                     nf,
 215                        int                     ndih,
 216                        real**                  dih,
 217                        real                    dt,
 218                        int                     nlist,
 219                        t_dlist                 dlist[],
 220                        real                    time[],
 221                        int                     maxchi,
 222                        gmx_bool                bPhi,
 223                        gmx_bool                bPsi,
 224                        gmx_bool                bChi,
 225                        gmx_bool                bOmega,
 226                        const gmx_output_env_t* oenv)
 227 {
 228     char name1[256], name2[256];
 229     int  i, j, Xi;
 230
 231     do_autocorr(fn, oenv, "Dihedral Autocorrelation Function", nf, ndih, dih, dt, eacCos, FALSE);
 232     /* Dump em all */
 233     j = 0;
 234     for (i = 0; (i < nlist); i++)
 235     {
 236         if (bPhi)
 237         {
 238             print_one(oenv, "corrphi", dlist[i].name, "Phi ACF for", "C(t)", nf / 2, time, dih[j]);
 239         }
 240         j++;
 241     }
 242     for (i = 0; (i < nlist); i++)
 243     {
 244         if (bPsi)
 245         {
 246             print_one(oenv, "corrpsi", dlist[i].name, "Psi ACF for", "C(t)", nf / 2, time, dih[j]);
 247         }
 248         j++;
 249     }
 250     for (i = 0; (i < nlist); i++)
 251     {
 252         if (has_dihedral(edOmega, &dlist[i]))
 253         {
 254             if (bOmega)
 255             {
 256                 print_one(oenv, "corromega", dlist[i].name, "Omega ACF for", "C(t)", nf / 2, time, dih[j]);
 257             }
 258             j++;
 259         }
 260     }
 261     for (Xi = 0; (Xi < maxchi); Xi++)
 262     {
 263         sprintf(name1, "corrchi%d", Xi + 1);
 264         sprintf(name2, "Chi%d ACF for", Xi + 1);
 265         for (i = 0; (i < nlist); i++)
 266         {
 267             if (dlist[i].atm.Cn[Xi + 3] != -1)
 268             {
 269                 if (bChi)
 270                 {
 271                     print_one(oenv, name1, dlist[i].name, name2, "C(t)", nf / 2, time, dih[j]);
 272                 }
 273                 j++;
 274             }
 275         }
 276     }
 277     fprintf(stderr, "\n");
 278 }
 279
 280 static void copy_dih_data(const real in[], real out[], int nf, gmx_bool bLEAVE)
 281 {
 282     /* if bLEAVE, do nothing to data in copying to out
 283      * otherwise multiply by 180/pi to convert rad to deg */
 284     int  i;
 285     real mult;
 286     if (bLEAVE)
 287     {
 288         mult = 1;
 289     }
 290     else
 291     {
 292         mult = (180.0 / M_PI);
 293     }
 294     for (i = 0; (i < nf); i++)
 295     {
 296         out[i] = in[i] * mult;
 297     }
 298 }
 299
 300 static void dump_em_all(int                     nlist,
 301                         t_dlist                 dlist[],
 302                         int                     nf,
 303                         real                    time[],
 304                         real**                  dih,
 305                         int                     maxchi,
 306                         gmx_bool                bPhi,
 307                         gmx_bool                bPsi,
 308                         gmx_bool                bChi,
 309                         gmx_bool                bOmega,
 310                         gmx_bool                bRAD,
 311                         const gmx_output_env_t* oenv)
 312 {
 313     char  name[256], titlestr[256], ystr[256];
 314     real* data;
 315     int   i, j, Xi;
 316
 317     snew(data, nf);
 318     if (bRAD)
 319     {
 320         std::strcpy(ystr, "Angle (rad)");
 321     }
 322     else
 323     {
 324         std::strcpy(ystr, "Angle (degrees)");
 325     }
 326
 327     /* Dump em all */
 328     j = 0;
 329     for (i = 0; (i < nlist); i++)
 330     {
 331         /* grs debug  printf("OK i %d j %d\n", i, j) ; */
 332         if (bPhi)
 333         {
 334             copy_dih_data(dih[j], data, nf, bRAD);
 335             print_one(oenv, "phi", dlist[i].name, "\\xf\\f{}", ystr, nf, time, data);
 336         }
 337         j++;
 338     }
 339     for (i = 0; (i < nlist); i++)
 340     {
 341         if (bPsi)
 342         {
 343             copy_dih_data(dih[j], data, nf, bRAD);
 344             print_one(oenv, "psi", dlist[i].name, "\\xy\\f{}", ystr, nf, time, data);
 345         }
 346         j++;
 347     }
 348     for (i = 0; (i < nlist); i++)
 349     {
 350         if (has_dihedral(edOmega, &(dlist[i])))
 351         {
 352             if (bOmega)
 353             {
 354                 copy_dih_data(dih[j], data, nf, bRAD);
 355                 print_one(oenv, "omega", dlist[i].name, "\\xw\\f{}", ystr, nf, time, data);
 356             }
 357             j++;
 358         }
 359     }
 360
 361     for (Xi = 0; (Xi < maxchi); Xi++)
 362     {
 363         for (i = 0; (i < nlist); i++)
 364         {
 365             if (dlist[i].atm.Cn[Xi + 3] != -1)
 366             {
 367                 if (bChi)
 368                 {
 369                     sprintf(name, "chi%d", Xi + 1);
 370                     sprintf(titlestr, "\\xc\\f{}\\s%d\\N", Xi + 1);
 371                     copy_dih_data(dih[j], data, nf, bRAD);
 372                     print_one(oenv, name, dlist[i].name, titlestr, ystr, nf, time, data);
 373                 }
 374                 j++;
 375             }
 376         }
 377     }
 378     fprintf(stderr, "\n");
 379 }
 380
 381 static void reset_one(real dih[], int nf, real phase)
 382 {
 383     int j;
 384
 385     for (j = 0; (j < nf); j++)
 386     {
 387         dih[j] += phase;
 388         while (dih[j] < -M_PI)
 389         {
 390             dih[j] += 2 * M_PI;
 391         }
 392         while (dih[j] >= M_PI)
 393         {
 394             dih[j] -= 2 * M_PI;
 395         }
 396     }
 397 }
 398
 399 static int reset_em_all(int nlist, t_dlist dlist[], int nf, real** dih, int maxchi)
 400 {
 401     int i, j, Xi;
 402
 403     /* Reset em all */
 404     j = 0;
 405     /* Phi */
 406     for (i = 0; (i < nlist); i++)
 407     {
 408         if (dlist[i].atm.minC == -1)
 409         {
 410             reset_one(dih[j++], nf, M_PI);
 411         }
 412         else
 413         {
 414             reset_one(dih[j++], nf, 0);
 415         }
 416     }
 417     /* Psi */
 418     for (i = 0; (i < nlist - 1); i++)
 419     {
 420         reset_one(dih[j++], nf, 0);
 421     }
 422     /* last Psi is faked from O */
 423     reset_one(dih[j++], nf, M_PI);
 424
 425     /* Omega */
 426     for (i = 0; (i < nlist); i++)
 427     {
 428         if (has_dihedral(edOmega, &dlist[i]))
 429         {
 430             reset_one(dih[j++], nf, 0);
 431         }
 432     }
 433     /* Chi 1 thru maxchi */
 434     for (Xi = 0; (Xi < maxchi); Xi++)
 435     {
 436         for (i = 0; (i < nlist); i++)
 437         {
 438             if (dlist[i].atm.Cn[Xi + 3] != -1)
 439             {
 440                 reset_one(dih[j], nf, 0);
 441                 j++;
 442             }
 443         }
 444     }
 445     fprintf(stderr, "j after resetting (nr. active dihedrals) = %d\n", j);
 446     return j;
 447 }
 448
 449 static void histogramming(FILE*                   log,
 450                           int                     nbin,
 451                           ResidueType*            rt,
 452                           int                     nf,
 453                           int                     maxchi,
 454                           real**                  dih,
 455                           int                     nlist,
 456                           t_dlist                 dlist[],
 457                           const int               index[],
 458                           gmx_bool                bPhi,
 459                           gmx_bool                bPsi,
 460                           gmx_bool                bOmega,
 461                           gmx_bool                bChi,
 462                           gmx_bool                bNormalize,
 463                           gmx_bool                bSSHisto,
 464                           const char*             ssdump,
 465                           real                    bfac_max,
 466                           const t_atoms*          atoms,
 467                           gmx_bool                bDo_jc,
 468                           const char*             fn,
 469                           const gmx_output_env_t* oenv)
 470 {
 471     /* also gets 3J couplings and order parameters S2 */
 472     // Avoid warnings about narrowing conversions from double to real
 473 #ifdef _MSC_VER
 474 #    pragma warning(disable : 4838)
 475 #endif
 476     t_karplus kkkphi[]  = { { "J_NHa1", 6.51, -1.76, 1.6, -M_PI / 3, 0.0, 0.0 },
 477                            { "J_NHa2", 6.51, -1.76, 1.6, M_PI / 3, 0.0, 0.0 },
 478                            { "J_HaC'", 4.0, 1.1, 0.1, 0.0, 0.0, 0.0 },
 479                            { "J_NHCb", 4.7, -1.5, -0.2, M_PI / 3, 0.0, 0.0 },
 480                            { "J_Ci-1Hai", 4.5, -1.3, -1.2, 2 * M_PI / 3, 0.0, 0.0 } };
 481     t_karplus kkkpsi[]  = { { "J_HaN", -0.88, -0.61, -0.27, M_PI / 3, 0.0, 0.0 } };
 482     t_karplus kkkchi1[] = { { "JHaHb2", 9.5, -1.6, 1.8, -M_PI / 3, 0, 0.0 },
 483                             { "JHaHb3", 9.5, -1.6, 1.8, 0, 0, 0.0 } };
 484 #ifdef _MSC_VER
 485 #    pragma warning(default : 4838)
 486 #endif
 487 #define NKKKPHI asize(kkkphi)
 488 #define NKKKPSI asize(kkkpsi)
 489 #define NKKKCHI asize(kkkchi1)
 490 #define NJC (NKKKPHI + NKKKPSI + NKKKCHI)
 491
 492     FILE *      fp, *ssfp[3] = { nullptr, nullptr, nullptr };
 493     const char* sss[3] = { "sheet", "helix", "coil" };
 494     real        S2;
 495     real*       normhisto;
 496     real **     Jc, **Jcsig;
 497     int****     his_aa_ss = nullptr;
 498     int ***     his_aa, *histmp;
 499     int         i, j, k, m, n, nn, Dih, nres, hindex, angle;
 500     gmx_bool    bBfac, bOccup;
 501     char        hisfile[256], hhisfile[256], title[256], *ss_str = nullptr;
 502     char**      leg;
 503     const char* residue_name;
 504     int         rt_size;
 505
 506     rt_size = rt->numberOfEntries();
 507     if (bSSHisto)
 508     {
 509         fp = gmx_ffopen(ssdump, "r");
 510         if (1 != fscanf(fp, "%d", &nres))
 511         {
 512             gmx_fatal(FARGS, "Error reading from file %s", ssdump);
 513         }
 514
 515         snew(ss_str, nres + 1);
 516         if (1 != fscanf(fp, "%s", ss_str))
 517         {
 518             gmx_fatal(FARGS, "Error reading from file %s", ssdump);
 519         }
 520
 521         gmx_ffclose(fp);
 522         /* Four dimensional array... Very cool */
 523         snew(his_aa_ss, 3);
 524         for (i = 0; (i < 3); i++)
 525         {
 526             snew(his_aa_ss[i], rt_size + 1);
 527             for (j = 0; (j <= rt_size); j++)
 528             {
 529                 snew(his_aa_ss[i][j], edMax);
 530                 for (Dih = 0; (Dih < edMax); Dih++)
 531                 {
 532                     snew(his_aa_ss[i][j][Dih], nbin + 1);
 533                 }
 534             }
 535         }
 536     }
 537     snew(his_aa, edMax);
 538     for (Dih = 0; (Dih < edMax); Dih++)
 539     {
 540         snew(his_aa[Dih], rt_size + 1);
 541         for (i = 0; (i <= rt_size); i++)
 542         {
 543             snew(his_aa[Dih][i], nbin + 1);
 544         }
 545     }
 546     snew(histmp, nbin);
 547
 548     snew(Jc, nlist);
 549     snew(Jcsig, nlist);
 550     for (i = 0; (i < nlist); i++)
 551     {
 552         snew(Jc[i], NJC);
 553         snew(Jcsig[i], NJC);
 554     }
 555
 556     j = 0;
 557     n = 0;
 558     for (Dih = 0; (Dih < NONCHI + maxchi); Dih++)
 559     {
 560         for (i = 0; (i < nlist); i++)
 561         {
 562             if (((Dih < edOmega)) || ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i]))))
 563                 || ((Dih > edOmega) && (dlist[i].atm.Cn[Dih - NONCHI + 3] != -1)))
 564             {
 565                 make_histo(log, nf, dih[j], nbin, histmp, -M_PI, M_PI);
 566
 567                 if (bSSHisto)
 568                 {
 569                     /* Assume there is only one structure, the first.
 570                      * Compute index in histogram.
 571                      */
 572                     /* Check the atoms to see whether their B-factors are low enough
 573                      * Check atoms to see their occupancy is 1.
 574                      */
 575                     bBfac = bOccup = TRUE;
 576                     for (nn = 0; (nn < 4); nn++, n++)
 577                     {
 578                         bBfac  = bBfac && (atoms->pdbinfo[index[n]].bfac <= bfac_max);
 579                         bOccup = bOccup && (atoms->pdbinfo[index[n]].occup == 1);
 580                     }
 581                     if (bOccup && ((bfac_max <= 0) || bBfac))
 582                     {
 583                         hindex = static_cast<int>(((dih[j][0] + M_PI) * nbin) / (2 * M_PI));
 584                         range_check(hindex, 0, nbin);
 585
 586                         /* Assign dihedral to either of the structure determined
 587                          * histograms
 588                          */
 589                         switch (ss_str[dlist[i].resnr])
 590                         {
 591                             case 'E': his_aa_ss[0][dlist[i].index][Dih][hindex]++; break;
 592                             case 'H': his_aa_ss[1][dlist[i].index][Dih][hindex]++; break;
 593                             default: his_aa_ss[2][dlist[i].index][Dih][hindex]++; break;
 594                         }
 595                     }
 596                     else if (debug)
 597                     {
 598                         fprintf(debug, "Res. %d has imcomplete occupancy or bfacs > %g\n", dlist[i].resnr, bfac_max);
 599                     }
 600                 }
 601                 else
 602                 {
 603                     n += 4;
 604                 }
 605
 606                 switch (Dih)
 607                 {
 608                     case edPhi:
 609                         calc_distribution_props(nbin, histmp, -M_PI, NKKKPHI, kkkphi, &S2);
 610
 611                         for (m = 0; (m < NKKKPHI); m++)
 612                         {
 613                             Jc[i][m]    = kkkphi[m].Jc;
 614                             Jcsig[i][m] = kkkphi[m].Jcsig;
 615                         }
 616                         break;
 617                     case edPsi:
 618                         calc_distribution_props(nbin, histmp, -M_PI, NKKKPSI, kkkpsi, &S2);
 619
 620                         for (m = 0; (m < NKKKPSI); m++)
 621                         {
 622                             Jc[i][NKKKPHI + m]    = kkkpsi[m].Jc;
 623                             Jcsig[i][NKKKPHI + m] = kkkpsi[m].Jcsig;
 624                         }
 625                         break;
 626                     case edChi1:
 627                         calc_distribution_props(nbin, histmp, -M_PI, NKKKCHI, kkkchi1, &S2);
 628                         for (m = 0; (m < NKKKCHI); m++)
 629                         {
 630                             Jc[i][NKKKPHI + NKKKPSI + m]    = kkkchi1[m].Jc;
 631                             Jcsig[i][NKKKPHI + NKKKPSI + m] = kkkchi1[m].Jcsig;
 632                         }
 633                         break;
 634                     default: /* covers edOmega and higher Chis than Chi1 */
 635                         calc_distribution_props(nbin, histmp, -M_PI, 0, nullptr, &S2);
 636                         break;
 637                 }
 638                 dlist[i].S2[Dih] = S2;
 639
 640                 /* Sum distribution per amino acid type as well */
 641                 for (k = 0; (k < nbin); k++)
 642                 {
 643                     his_aa[Dih][dlist[i].index][k] += histmp[k];
 644                     histmp[k] = 0;
 645                 }
 646                 j++;
 647             }
 648             else /* dihed not defined */
 649             {
 650                 dlist[i].S2[Dih] = 0.0;
 651             }
 652         }
 653     }
 654     sfree(histmp);
 655
 656     /* Print out Jcouplings */
 657     fprintf(log, "\n *** J-Couplings from simulation (plus std. dev.) ***\n\n");
 658     fprintf(log, "Residue   ");
 659     for (i = 0; (i < NKKKPHI); i++)
 660     {
 661         fprintf(log, "%7s   SD", kkkphi[i].name);
 662     }
 663     for (i = 0; (i < NKKKPSI); i++)
 664     {
 665         fprintf(log, "%7s   SD", kkkpsi[i].name);
 666     }
 667     for (i = 0; (i < NKKKCHI); i++)
 668     {
 669         fprintf(log, "%7s   SD", kkkchi1[i].name);
 670     }
 671     fprintf(log, "\n");
 672     for (i = 0; (i < NJC + 1); i++)
 673     {
 674         fprintf(log, "------------");
 675     }
 676     fprintf(log, "\n");
 677     for (i = 0; (i < nlist); i++)
 678     {
 679         fprintf(log, "%-10s", dlist[i].name);
 680         for (j = 0; (j < NJC); j++)
 681         {
 682             fprintf(log, "  %5.2f %4.2f", Jc[i][j], Jcsig[i][j]);
 683         }
 684         fprintf(log, "\n");
 685     }
 686     fprintf(log, "\n");
 687
 688     /* and to -jc file... */
 689     if (bDo_jc)
 690     {
 691         fp = xvgropen(fn, "\\S3\\NJ-Couplings from Karplus Equation", "Residue", "Coupling", oenv);
 692         snew(leg, NJC);
 693         for (i = 0; (i < NKKKPHI); i++)
 694         {
 695             leg[i] = gmx_strdup(kkkphi[i].name);
 696         }
 697         for (i = 0; (i < NKKKPSI); i++)
 698         {
 699             leg[i + NKKKPHI] = gmx_strdup(kkkpsi[i].name);
 700         }
 701         for (i = 0; (i < NKKKCHI); i++)
 702         {
 703             leg[i + NKKKPHI + NKKKPSI] = gmx_strdup(kkkchi1[i].name);
 704         }
 705         xvgr_legend(fp, NJC, leg, oenv);
 706         fprintf(fp, "%5s ", "#Res.");
 707         for (i = 0; (i < NJC); i++)
 708         {
 709             fprintf(fp, "%10s ", leg[i]);
 710         }
 711         fprintf(fp, "\n");
 712         for (i = 0; (i < nlist); i++)
 713         {
 714             fprintf(fp, "%5d ", dlist[i].resnr);
 715             for (j = 0; (j < NJC); j++)
 716             {
 717                 fprintf(fp, "  %8.3f", Jc[i][j]);
 718             }
 719             fprintf(fp, "\n");
 720         }
 721         xvgrclose(fp);
 722         for (i = 0; (i < NJC); i++)
 723         {
 724             sfree(leg[i]);
 725         }
 726     }
 727     /* finished -jc stuff */
 728
 729     snew(normhisto, nbin);
 730     for (i = 0; (i < rt_size); i++)
 731     {
 732         for (Dih = 0; (Dih < edMax); Dih++)
 733         {
 734             /* First check whether something is in there */
 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(nbin, his_aa[Dih][i], (360.0 / nbin), normhisto);
 749                 }
 750
 751                 std::string residueName = rt->nameFromResidueIndex(i);
 752                 residue_name            = residueName.c_str();
 753                 switch (Dih)
 754                 {
 755                     case edPhi:
 756                         sprintf(hisfile, "histo-phi%s", residue_name);
 757                         sprintf(title, "\\xf\\f{} Distribution for %s", residue_name);
 758                         break;
 759                     case edPsi:
 760                         sprintf(hisfile, "histo-psi%s", residue_name);
 761                         sprintf(title, "\\xy\\f{} Distribution for %s", residue_name);
 762                         break;
 763                     case edOmega:
 764                         sprintf(hisfile, "histo-omega%s", residue_name);
 765                         sprintf(title, "\\xw\\f{} Distribution for %s", residue_name);
 766                         break;
 767                     default:
 768                         sprintf(hisfile, "histo-chi%d%s", Dih - NONCHI + 1, residue_name);
 769                         sprintf(title, "\\xc\\f{}\\s%d\\N Distribution for %s", Dih - NONCHI + 1, residue_name);
 770                 }
 771                 std::strcpy(hhisfile, hisfile);
 772                 std::strcat(hhisfile, ".xvg");
 773                 fp = xvgropen(hhisfile, title, "Degrees", "", oenv);
 774                 if (output_env_get_print_xvgr_codes(oenv))
 775                 {
 776                     fprintf(fp, "@ with g0\n");
 777                 }
 778                 xvgr_world(fp, -180, 0, 180, 0.1, oenv);
 779                 if (output_env_get_print_xvgr_codes(oenv))
 780                 {
 781                     fprintf(fp,
 782                             "# this effort to set graph size fails unless you run with -autoscale "
 783                             "none or -autoscale y flags\n");
 784                     fprintf(fp, "@ xaxis tick on\n");
 785                     fprintf(fp, "@ xaxis tick major 90\n");
 786                     fprintf(fp, "@ xaxis tick minor 30\n");
 787                     fprintf(fp, "@ xaxis ticklabel prec 0\n");
 788                     fprintf(fp, "@ yaxis tick off\n");
 789                     fprintf(fp, "@ yaxis ticklabel off\n");
 790                     fprintf(fp, "@ type xy\n");
 791                 }
 792                 if (bSSHisto)
 793                 {
 794                     for (k = 0; (k < 3); k++)
 795                     {
 796                         std::string sshisfile = gmx::formatString("%s-%s.xvg", hisfile, sss[k]);
 797                         ssfp[k]               = gmx_ffopen(sshisfile, "w");
 798                     }
 799                 }
 800                 for (j = 0; (j < nbin); j++)
 801                 {
 802                     angle = -180 + (360 / nbin) * j;
 803                     if (bNormalize)
 804                     {
 805                         fprintf(fp, "%5d  %10g\n", angle, normhisto[j]);
 806                     }
 807                     else
 808                     {
 809                         fprintf(fp, "%5d  %10d\n", angle, his_aa[Dih][i][j]);
 810                     }
 811                     if (bSSHisto)
 812                     {
 813                         for (k = 0; (k < 3); k++)
 814                         {
 815                             fprintf(ssfp[k], "%5d  %10d\n", angle, his_aa_ss[k][i][Dih][j]);
 816                         }
 817                     }
 818                 }
 819                 fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 820                 xvgrclose(fp);
 821                 if (bSSHisto)
 822                 {
 823                     for (k = 0; (k < 3); k++)
 824                     {
 825                         fprintf(ssfp[k], "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 826                         gmx_ffclose(ssfp[k]);
 827                     }
 828                 }
 829             }
 830         }
 831     }
 832     sfree(normhisto);
 833
 834     if (bSSHisto)
 835     {
 836         /* Four dimensional array... Very cool */
 837         for (i = 0; (i < 3); i++)
 838         {
 839             for (j = 0; (j <= rt_size); j++)
 840             {
 841                 for (Dih = 0; (Dih < edMax); Dih++)
 842                 {
 843                     sfree(his_aa_ss[i][j][Dih]);
 844                 }
 845                 sfree(his_aa_ss[i][j]);
 846             }
 847             sfree(his_aa_ss[i]);
 848         }
 849         sfree(his_aa_ss);
 850         sfree(ss_str);
 851     }
 852 }
 853
 854 static FILE* rama_file(const char* fn, const char* title, const char* xaxis, const char* yaxis, const gmx_output_env_t* oenv)
 855 {
 856     FILE* fp;
 857
 858     fp = xvgropen(fn, title, xaxis, yaxis, oenv);
 859     if (output_env_get_print_xvgr_codes(oenv))
 860     {
 861         fprintf(fp, "@ with g0\n");
 862     }
 863     xvgr_world(fp, -180, -180, 180, 180, oenv);
 864     if (output_env_get_print_xvgr_codes(oenv))
 865     {
 866         fprintf(fp, "@ xaxis tick on\n");
 867         fprintf(fp, "@ xaxis tick major 90\n");
 868         fprintf(fp, "@ xaxis tick minor 30\n");
 869         fprintf(fp, "@ xaxis ticklabel prec 0\n");
 870         fprintf(fp, "@ yaxis tick on\n");
 871         fprintf(fp, "@ yaxis tick major 90\n");
 872         fprintf(fp, "@ yaxis tick minor 30\n");
 873         fprintf(fp, "@ yaxis ticklabel prec 0\n");
 874         fprintf(fp, "@    s0 type xy\n");
 875         fprintf(fp, "@    s0 symbol 2\n");
 876         fprintf(fp, "@    s0 symbol size 0.410000\n");
 877         fprintf(fp, "@    s0 symbol fill 1\n");
 878         fprintf(fp, "@    s0 symbol color 1\n");
 879         fprintf(fp, "@    s0 symbol linewidth 1\n");
 880         fprintf(fp, "@    s0 symbol linestyle 1\n");
 881         fprintf(fp, "@    s0 symbol center false\n");
 882         fprintf(fp, "@    s0 symbol char 0\n");
 883         fprintf(fp, "@    s0 skip 0\n");
 884         fprintf(fp, "@    s0 linestyle 0\n");
 885         fprintf(fp, "@    s0 linewidth 1\n");
 886         fprintf(fp, "@ type xy\n");
 887     }
 888     return fp;
 889 }
 890
 891 static void do_rama(int                     nf,
 892                     int                     nlist,
 893                     t_dlist                 dlist[],
 894                     real**                  dih,
 895                     gmx_bool                bViol,
 896                     gmx_bool                bRamOmega,
 897                     const gmx_output_env_t* oenv)
 898 {
 899     FILE *        fp, *gp = nullptr;
 900     gmx_bool      bOm;
 901     char          fn[256];
 902     int           i, j, k, Xi1, Xi2, Phi, Psi, Om = 0, nlevels;
 903     constexpr int NMAT = 120;
 904     real **       mat  = nullptr, phi, psi, omega, axis[NMAT], lo, hi;
 905     t_rgb         rlo  = { 1.0, 0.0, 0.0 };
 906     t_rgb         rmid = { 1.0, 1.0, 1.0 };
 907     t_rgb         rhi  = { 0.0, 0.0, 1.0 };
 908
 909     for (i = 0; (i < nlist); i++)
 910     {
 911         if ((has_dihedral(edPhi, &(dlist[i]))) && (has_dihedral(edPsi, &(dlist[i]))))
 912         {
 913             sprintf(fn, "ramaPhiPsi%s.xvg", dlist[i].name);
 914             fp  = rama_file(fn, "Ramachandran Plot", "\\8f\\4 (deg)", "\\8y\\4 (deg)", oenv);
 915             bOm = bRamOmega && has_dihedral(edOmega, &(dlist[i]));
 916             if (bOm)
 917             {
 918                 Om = dlist[i].j0[edOmega];
 919                 snew(mat, NMAT);
 920                 for (j = 0; (j < NMAT); j++)
 921                 {
 922                     snew(mat[j], NMAT);
 923                     axis[j] = -180 + gmx::exactDiv(360 * j, NMAT);
 924                 }
 925             }
 926             if (bViol)
 927             {
 928                 sprintf(fn, "violPhiPsi%s.xvg", dlist[i].name);
 929                 gp = gmx_ffopen(fn, "w");
 930             }
 931             Phi = dlist[i].j0[edPhi];
 932             Psi = dlist[i].j0[edPsi];
 933             for (j = 0; (j < nf); j++)
 934             {
 935                 phi = gmx::c_rad2Deg * dih[Phi][j];
 936                 psi = gmx::c_rad2Deg * dih[Psi][j];
 937                 fprintf(fp, "%10g  %10g\n", phi, psi);
 938                 if (bViol)
 939                 {
 940                     fprintf(gp,
 941                             "%d\n",
 942                             static_cast<int>(!bAllowed(dih[Phi][j], gmx::c_rad2Deg * dih[Psi][j])));
 943                 }
 944                 if (bOm)
 945                 {
 946                     omega = gmx::c_rad2Deg * dih[Om][j];
 947                     mat[static_cast<int>(((phi * NMAT) / 360) + gmx::exactDiv(NMAT, 2))]
 948                        [static_cast<int>(((psi * NMAT) / 360) + gmx::exactDiv(NMAT, 2))] += omega;
 949                 }
 950             }
 951             if (bViol)
 952             {
 953                 gmx_ffclose(gp);
 954             }
 955             xvgrclose(fp);
 956             if (bOm)
 957             {
 958                 sprintf(fn, "ramomega%s.xpm", dlist[i].name);
 959                 fp = gmx_ffopen(fn, "w");
 960                 lo = hi = 0;
 961                 for (j = 0; (j < NMAT); j++)
 962                 {
 963                     for (k = 0; (k < NMAT); k++)
 964                     {
 965                         mat[j][k] /= nf;
 966                         lo = std::min(mat[j][k], lo);
 967                         hi = std::max(mat[j][k], hi);
 968                     }
 969                 }
 970                 /* Symmetrise */
 971                 if (std::abs(lo) > std::abs(hi))
 972                 {
 973                     hi = -lo;
 974                 }
 975                 else
 976                 {
 977                     lo = -hi;
 978                 }
 979                 /* Add 180 */
 980                 for (j = 0; (j < NMAT); j++)
 981                 {
 982                     for (k = 0; (k < NMAT); k++)
 983                     {
 984                         mat[j][k] += 180;
 985                     }
 986                 }
 987                 lo += 180;
 988                 hi += 180;
 989                 nlevels = 20;
 990                 write_xpm3(fp,
 991                            0,
 992                            "Omega/Ramachandran Plot",
 993                            "Deg",
 994                            "Phi",
 995                            "Psi",
 996                            NMAT,
 997                            NMAT,
 998                            axis,
 999                            axis,
1000                            mat,
1001                            lo,
1002                            180.0,
1003                            hi,
1004                            rlo,
1005                            rmid,
1006                            rhi,
1007                            &nlevels);
1008                 gmx_ffclose(fp);
1009                 for (j = 0; (j < NMAT); j++)
1010                 {
1011                     sfree(mat[j]);
1012                 }
1013                 sfree(mat);
1014             }
1015         }
1016         if ((has_dihedral(edChi1, &(dlist[i]))) && (has_dihedral(edChi2, &(dlist[i]))))
1017         {
1018             sprintf(fn, "ramaX1X2%s.xvg", dlist[i].name);
1019             fp  = rama_file(fn,
1020                            "\\8c\\4\\s1\\N-\\8c\\4\\s2\\N Ramachandran Plot",
1021                            "\\8c\\4\\s1\\N (deg)",
1022                            "\\8c\\4\\s2\\N (deg)",
1023                            oenv);
1024             Xi1 = dlist[i].j0[edChi1];
1025             Xi2 = dlist[i].j0[edChi2];
1026             for (j = 0; (j < nf); j++)
1027             {
1028                 fprintf(fp, "%10g  %10g\n", gmx::c_rad2Deg * dih[Xi1][j], gmx::c_rad2Deg * dih[Xi2][j]);
1029             }
1030             xvgrclose(fp);
1031         }
1032         else
1033         {
1034             fprintf(stderr, "No chi1 & chi2 angle for %s\n", dlist[i].name);
1035         }
1036     }
1037 }
1038
1039
1040 static void print_transitions(const char* fn, int maxchi, int nlist, t_dlist dlist[], real dt, const gmx_output_env_t* oenv)
1041 {
1042     /* based on order_params below */
1043     FILE* fp;
1044     int   i, Dih, Xi;
1045
1046     /*  must correspond with enum in pp2shift.h:38 */
1047     char* leg[edMax];
1048
1049     leg[0] = gmx_strdup("Phi");
1050     leg[1] = gmx_strdup("Psi");
1051     leg[2] = gmx_strdup("Omega");
1052     leg[3] = gmx_strdup("Chi1");
1053     leg[4] = gmx_strdup("Chi2");
1054     leg[5] = gmx_strdup("Chi3");
1055     leg[6] = gmx_strdup("Chi4");
1056     leg[7] = gmx_strdup("Chi5");
1057     leg[8] = gmx_strdup("Chi6");
1058
1059     /* Print order parameters */
1060     fp = xvgropen(fn, "Dihedral Rotamer Transitions", "Residue", "Transitions/ns", oenv);
1061     xvgr_legend(fp, NONCHI + maxchi, leg, oenv);
1062
1063     fprintf(fp, "%5s ", "#Res.");
1064     fprintf(fp, "%10s %10s %10s ", leg[edPhi], leg[edPsi], leg[edOmega]);
1065     for (Xi = 0; Xi < maxchi; Xi++)
1066     {
1067         fprintf(fp, "%10s ", leg[NONCHI + Xi]);
1068     }
1069     fprintf(fp, "\n");
1070
1071     for (i = 0; (i < nlist); i++)
1072     {
1073         fprintf(fp, "%5d ", dlist[i].resnr);
1074         for (Dih = 0; (Dih < NONCHI + maxchi); Dih++)
1075         {
1076             fprintf(fp, "%10.3f ", dlist[i].ntr[Dih] / dt);
1077         }
1078         /* fprintf(fp,"%12s\n",dlist[i].name);  this confuses xmgrace */
1079         fprintf(fp, "\n");
1080     }
1081     xvgrclose(fp);
1082 }
1083
1084 static void order_params(FILE*                   log,
1085                          const char*             fn,
1086                          int                     maxchi,
1087                          int                     nlist,
1088                          t_dlist                 dlist[],
1089                          const char*             pdbfn,
1090                          real                    bfac_init,
1091                          t_atoms*                atoms,
1092                          const rvec              x[],
1093                          PbcType                 pbcType,
1094                          matrix                  box,
1095                          gmx_bool                bPhi,
1096                          gmx_bool                bPsi,
1097                          gmx_bool                bChi,
1098                          const gmx_output_env_t* oenv)
1099 {
1100     FILE* fp;
1101     int   nh[edMax];
1102     int   i, Dih, Xi;
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 (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 (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 (Xi = 0; Xi < maxchi; Xi++)
1130     {
1131         fprintf(fp, "%10s ", leg[2 + NONCHI + Xi]);
1132     }
1133     fprintf(fp, "\n");
1134
1135     for (i = 0; (i < nlist); i++)
1136     {
1137         S2Max = -10;
1138         S2Min = 10;
1139         for (Dih = 0; (Dih < NONCHI + maxchi); Dih++)
1140         {
1141             if (dlist[i].S2[Dih] != 0)
1142             {
1143                 if (dlist[i].S2[Dih] > S2Max)
1144                 {
1145                     S2Max = dlist[i].S2[Dih];
1146                 }
1147                 if (dlist[i].S2[Dih] < S2Min)
1148                 {
1149                     S2Min = dlist[i].S2[Dih];
1150                 }
1151             }
1152             if (dlist[i].S2[Dih] > 0.8)
1153             {
1154                 nh[Dih]++;
1155             }
1156         }
1157         fprintf(fp, "%5d ", dlist[i].resnr);
1158         fprintf(fp, "%10.3f %10.3f ", S2Min, S2Max);
1159         for (Dih = 0; (Dih < NONCHI + maxchi); Dih++)
1160         {
1161             fprintf(fp, "%10.3f ", dlist[i].S2[Dih]);
1162         }
1163         fprintf(fp, "\n");
1164         /* fprintf(fp,"%12s\n",dlist[i].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 (i = 0; (i < atoms->nr); i++)
1179         {
1180             atoms->pdbinfo[i].bfac = bfac_init;
1181         }
1182
1183         for (i = 0; (i < nlist); i++)
1184         {
1185             atoms->pdbinfo[dlist[i].atm.N].bfac = -dlist[i].S2[0]; /* Phi */
1186             atoms->pdbinfo[dlist[i].atm.H].bfac = -dlist[i].S2[0]; /* Phi */
1187             atoms->pdbinfo[dlist[i].atm.C].bfac = -dlist[i].S2[1]; /* Psi */
1188             atoms->pdbinfo[dlist[i].atm.O].bfac = -dlist[i].S2[1]; /* Psi */
1189             for (Xi = 0; (Xi < maxchi); Xi++)                      /* Chi's */
1190             {
1191                 if (dlist[i].atm.Cn[Xi + 3] != -1)
1192                 {
1193                     atoms->pdbinfo[dlist[i].atm.Cn[Xi + 1]].bfac = -dlist[i].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 (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 (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 (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 (Xi = 0; (Xi < maxchi); Xi++)
1264         {
1265             fprintf(log, "%4d  ", nh[NONCHI + Xi]);
1266         }
1267     }
1268     fprintf(log, "\n");
1269
1270     for (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               nlist, idum, nbin;
1434     rvec*             x;
1435     PbcType           pbcType;
1436     matrix            box;
1437     char              grpname[256];
1438     t_dlist*          dlist;
1439     gmx_bool          bChi, bCorr, bSSHisto;
1440     gmx_bool          bDo_rt, bDo_oh, bDo_ot, bDo_jc;
1441     real              dt = 0, traj_t_ns;
1442     gmx_output_env_t* oenv;
1443
1444     int    isize, *index;
1445     int    ndih, nactdih, nf;
1446     real **dih, *trans_frac, *aver_angle, *time;
1447     int    i, **chi_lookup, *multiplicity;
1448
1449     t_filenm fnm[] = { { efSTX, "-s", nullptr, ffREAD },
1450                        { efTRX, "-f", nullptr, ffREAD },
1451                        { efXVG, "-o", "order", ffWRITE },
1452                        { efPDB, "-p", "order", ffOPTWR },
1453                        { efDAT, "-ss", "ssdump", ffOPTRD },
1454                        { efXVG, "-jc", "Jcoupling", ffWRITE },
1455                        { efXVG, "-corr", "dihcorr", ffOPTWR },
1456                        { efLOG, "-g", "chi", ffWRITE },
1457                        /* add two more arguments copying from g_angle */
1458                        { efXVG, "-ot", "dihtrans", ffOPTWR },
1459                        { efXVG, "-oh", "trhisto", ffOPTWR },
1460                        { efXVG, "-rt", "restrans", ffOPTWR },
1461                        { efXVG, "-cp", "chiprodhisto", ffOPTWR } };
1462 #define NFILE asize(fnm)
1463     int      npargs;
1464     t_pargs* ppa;
1465
1466     npargs = asize(pa);
1467     ppa    = add_acf_pargs(&npargs, pa);
1468     if (!parse_common_args(
1469                 &argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs, &oenv))
1470     {
1471         sfree(ppa);
1472         return 0;
1473     }
1474
1475     /* Handle result from enumerated type */
1476     sscanf(maxchistr[0], "%d", &maxchi);
1477     bChi = (maxchi > 0);
1478
1479     log = gmx_ffopen(ftp2fn(efLOG, NFILE, fnm), "w");
1480
1481     if (bRamOmega)
1482     {
1483         bOmega = TRUE;
1484         bPhi   = TRUE;
1485         bPsi   = TRUE;
1486     }
1487
1488     /* set some options */
1489     bDo_rt = (opt2bSet("-rt", NFILE, fnm));
1490     bDo_oh = (opt2bSet("-oh", NFILE, fnm));
1491     bDo_ot = (opt2bSet("-ot", NFILE, fnm));
1492     bDo_jc = (opt2bSet("-jc", NFILE, fnm));
1493     bCorr  = (opt2bSet("-corr", NFILE, fnm));
1494     if (bCorr)
1495     {
1496         fprintf(stderr, "Will calculate autocorrelation\n");
1497     }
1498
1499     if (core_frac > 1.0)
1500     {
1501         fprintf(stderr, "core_rotamer fraction > 1.0 ; will use 1.0\n");
1502         core_frac = 1.0;
1503     }
1504     if (core_frac < 0.0)
1505     {
1506         fprintf(stderr, "core_rotamer fraction < 0.0 ; will use 0.0\n");
1507         core_frac = 0.0;
1508     }
1509
1510     if (maxchi > MAXCHI)
1511     {
1512         fprintf(stderr, "Will only calculate first %d Chi dihedrals instead of %d.\n", MAXCHI, maxchi);
1513         maxchi = MAXCHI;
1514     }
1515     bSSHisto = ftp2bSet(efDAT, NFILE, fnm);
1516     nbin     = 360 / ndeg;
1517
1518     /* Find the chi angles using atoms struct and a list of amino acids */
1519     t_topology* top;
1520     snew(top, 1);
1521     read_tps_conf(ftp2fn(efSTX, NFILE, fnm), top, &pbcType, &x, nullptr, box, FALSE);
1522     t_atoms& atoms = top->atoms;
1523     if (atoms.pdbinfo == nullptr)
1524     {
1525         snew(atoms.pdbinfo, atoms.nr);
1526     }
1527     fprintf(log, "Title: %s\n", *top->name);
1528
1529     ResidueType rt;
1530     dlist = mk_dlist(log, &atoms, &nlist, bPhi, bPsi, bChi, bHChi, maxchi, r0, &rt);
1531     fprintf(stderr, "%d residues with dihedrals found\n", nlist);
1532
1533     if (nlist == 0)
1534     {
1535         gmx_fatal(FARGS, "No dihedrals in your structure!\n");
1536     }
1537
1538     /* Make a linear index for reading all. */
1539     index = make_chi_ind(nlist, dlist, &ndih);
1540     isize = 4 * ndih;
1541     fprintf(stderr, "%d dihedrals found\n", ndih);
1542
1543     snew(dih, ndih);
1544
1545     /* COMPUTE ALL DIHEDRALS! */
1546     read_ang_dih(
1547             ftp2fn(efTRX, NFILE, fnm), FALSE, TRUE, FALSE, bPBC, 1, &idum, &nf, &time, isize, index, &trans_frac, &aver_angle, dih, oenv);
1548
1549     dt = (time[nf - 1] - time[0]) / (nf - 1); /* might want this for corr or n. transit*/
1550     if (bCorr)
1551     {
1552         if (nf < 2)
1553         {
1554             gmx_fatal(FARGS, "Need at least 2 frames for correlation");
1555         }
1556     }
1557
1558     /* put angles in -M_PI to M_PI ! and correct phase factor for phi and psi
1559      * pass nactdih instead of ndih to low_ana_dih_trans
1560      * to prevent accessing off end of arrays when maxchi < 5 or 6. */
1561     nactdih = reset_em_all(nlist, dlist, nf, dih, maxchi);
1562
1563     if (bAll)
1564     {
1565         dump_em_all(nlist, dlist, nf, time, dih, maxchi, bPhi, bPsi, bChi, bOmega, bRAD, oenv);
1566     }
1567
1568     /* Histogramming & J coupling constants & calc of S2 order params */
1569     histogramming(log,
1570                   nbin,
1571                   &rt,
1572                   nf,
1573                   maxchi,
1574                   dih,
1575                   nlist,
1576                   dlist,
1577                   index,
1578                   bPhi,
1579                   bPsi,
1580                   bOmega,
1581                   bChi,
1582                   bNormHisto,
1583                   bSSHisto,
1584                   ftp2fn(efDAT, NFILE, fnm),
1585                   bfac_max,
1586                   &atoms,
1587                   bDo_jc,
1588                   opt2fn("-jc", NFILE, fnm),
1589                   oenv);
1590
1591     /* transitions
1592      *
1593      * added multiplicity */
1594
1595     snew(multiplicity, ndih);
1596     mk_multiplicity_lookup(multiplicity, maxchi, nlist, dlist, ndih);
1597
1598     std::strcpy(grpname, "All residues, ");
1599     if (bPhi)
1600     {
1601         std::strcat(grpname, "Phi ");
1602     }
1603     if (bPsi)
1604     {
1605         std::strcat(grpname, "Psi ");
1606     }
1607     if (bOmega)
1608     {
1609         std::strcat(grpname, "Omega ");
1610     }
1611     if (bChi)
1612     {
1613         std::strcat(grpname, "Chi 1-");
1614         sprintf(grpname + std::strlen(grpname), "%i", maxchi);
1615     }
1616
1617
1618     low_ana_dih_trans(bDo_ot,
1619                       opt2fn("-ot", NFILE, fnm),
1620                       bDo_oh,
1621                       opt2fn("-oh", NFILE, fnm),
1622                       maxchi,
1623                       dih,
1624                       nlist,
1625                       dlist,
1626                       nf,
1627                       nactdih,
1628                       grpname,
1629                       multiplicity,
1630                       time,
1631                       FALSE,
1632                       core_frac,
1633                       oenv);
1634
1635     /* Order parameters */
1636     order_params(log,
1637                  opt2fn("-o", NFILE, fnm),
1638                  maxchi,
1639                  nlist,
1640                  dlist,
1641                  ftp2fn_null(efPDB, NFILE, fnm),
1642                  bfac_init,
1643                  &atoms,
1644                  x,
1645                  pbcType,
1646                  box,
1647                  bPhi,
1648                  bPsi,
1649                  bChi,
1650                  oenv);
1651
1652     /* Print ramachandran maps! */
1653     if (bRama)
1654     {
1655         do_rama(nf, nlist, dlist, dih, bViol, bRamOmega, oenv);
1656     }
1657
1658     if (bShift)
1659     {
1660         do_pp2shifts(log, nf, nlist, dlist, dih);
1661     }
1662
1663     /* rprint S^2, transitions, and rotamer occupancies to log */
1664     traj_t_ns = 0.001 * (time[nf - 1] - time[0]);
1665     pr_dlist(log, nlist, dlist, traj_t_ns, edPrintST, bPhi, bPsi, bChi, bOmega, maxchi);
1666     pr_dlist(log, nlist, dlist, traj_t_ns, edPrintRO, bPhi, bPsi, bChi, bOmega, maxchi);
1667     gmx_ffclose(log);
1668     /* transitions to xvg */
1669     if (bDo_rt)
1670     {
1671         print_transitions(opt2fn("-rt", NFILE, fnm), maxchi, nlist, dlist, traj_t_ns, oenv);
1672     }
1673
1674     /* chi_product trajectories (ie one "rotamer number" for each residue) */
1675     if (bChiProduct && bChi)
1676     {
1677         snew(chi_lookup, nlist);
1678         for (i = 0; i < nlist; i++)
1679         {
1680             snew(chi_lookup[i], maxchi);
1681         }
1682         mk_chi_lookup(chi_lookup, maxchi, nlist, dlist);
1683
1684         get_chi_product_traj(dih,
1685                              nf,
1686                              nlist,
1687                              maxchi,
1688                              dlist,
1689                              time,
1690                              chi_lookup,
1691                              multiplicity,
1692                              FALSE,
1693                              bNormHisto,
1694                              core_frac,
1695                              bAll,
1696                              opt2fn("-cp", NFILE, fnm),
1697                              oenv);
1698
1699         for (i = 0; i < nlist; i++)
1700         {
1701             sfree(chi_lookup[i]);
1702         }
1703     }
1704
1705     /* Correlation comes last because it messes up the angles */
1706     if (bCorr)
1707     {
1708         do_dihcorr(opt2fn("-corr", NFILE, fnm), nf, ndih, dih, dt, nlist, dlist, time, maxchi, bPhi, bPsi, bChi, bOmega, oenv);
1709     }
1710
1711
1712     do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1713     do_view(oenv, opt2fn("-jc", NFILE, fnm), "-nxy");
1714     if (bCorr)
1715     {
1716         do_view(oenv, opt2fn("-corr", NFILE, fnm), "-nxy");
1717     }
1718
1719     return 0;
1720 }