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