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