src/gromacs/gmxana/gmx_chi.c

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