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