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