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