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