src/tools/gmx_chi.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38 #include <stdio.h>
  39 #include <math.h>
  40
  41 #include "confio.h"
  42 #include "pdbio.h"
  43 #include "copyrite.h"
  44 #include "gmx_fatal.h"
  45 #include "futil.h"
  46 #include "gstat.h"
  47 #include "macros.h"
  48 #include "maths.h"
  49 #include "physics.h"
  50 #include "index.h"
  51 #include "smalloc.h"
  52 #include "statutil.h"
  53 #include "tpxio.h"
  54 #include "string.h"
  55 #include "sysstuff.h"
  56 #include "txtdump.h"
  57 #include "typedefs.h"
  58 #include "vec.h"
  59 #include "strdb.h"
  60 #include "xvgr.h"
  61 #include "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                   id[n++]=dl[i].atm.minC;
 154           else
 155                   id[n++]=dl[i].atm.H;
 156           id[n++]=dl[i].atm.N;
 157           id[n++]=dl[i].atm.Cn[1];
 158           id[n++]=dl[i].atm.C;
 159   }
 160   for(i=0; (i<nl); i++)
 161   {
 162           /* Psi, fake the last one */
 163           dl[i].j0[edPsi] = n/4;
 164           id[n++]=dl[i].atm.N;
 165           id[n++]=dl[i].atm.Cn[1];
 166           id[n++]=dl[i].atm.C;
 167           if ( i< (nl-1) )
 168                   id[n++]=dl[i+1].atm.N;
 169           else
 170                   id[n++]=dl[i].atm.O;
 171   }
 172   for(i=0; (i<nl); i++)
 173   {
 174           /* Omega */
 175           if (has_dihedral(edOmega,&(dl[i])))
 176           {
 177                   dl[i].j0[edOmega] = n/4;
 178                   id[n++]=dl[i].atm.minO;
 179                   id[n++]=dl[i].atm.minC;
 180                   id[n++]=dl[i].atm.N;
 181                   id[n++]=dl[i].atm.H;
 182           }
 183   }
 184   for(Xi=0; (Xi<MAXCHI); Xi++)
 185   {
 186           /* Chi# */
 187           for(i=0; (i<nl); i++)
 188           {
 189                   if (dl[i].atm.Cn[Xi+3] != -1)
 190                   {
 191                           dl[i].j0[edChi1+Xi] = n/4;
 192                           id[n++]=dl[i].atm.Cn[Xi];
 193                           id[n++]=dl[i].atm.Cn[Xi+1];
 194                           id[n++]=dl[i].atm.Cn[Xi+2];
 195                           id[n++]=dl[i].atm.Cn[Xi+3];
 196                   }
 197           }
 198   }
 199   *ndih=n/4;
 200
 201   return id;
 202 }
 203
 204 int bin(real chi,int mult)
 205 {
 206   mult=3;
 207
 208   return (int) (chi*mult/360.0);
 209 }
 210
 211
 212 static void do_dihcorr(const char *fn,int nf,int ndih,real **dih,real dt,
 213                        int nlist,t_dlist dlist[],real time[],int maxchi,
 214                        gmx_bool bPhi,gmx_bool bPsi,gmx_bool bChi,gmx_bool bOmega,
 215                        const output_env_t oenv)
 216 {
 217   char name1[256],name2[256];
 218   int  i,j,Xi;
 219
 220   do_autocorr(fn,oenv,"Dihedral Autocorrelation Function",
 221               nf,ndih,dih,dt,eacCos,FALSE);
 222   /* Dump em all */
 223   j=0;
 224   for(i=0; (i<nlist); i++) {
 225     if (bPhi)
 226       print_one(oenv,"corrphi",dlist[i].name,"Phi ACF for", "C(t)", nf/2,time,
 227                           dih[j]);
 228     j++;
 229   }
 230   for(i=0; (i<nlist); i++) {
 231     if (bPsi)
 232       print_one(oenv,"corrpsi",dlist[i].name,"Psi ACF for","C(t)",nf/2,time,
 233                 dih[j]);
 234     j++;
 235   }
 236   for(i=0; (i<nlist); i++) {
 237     if (has_dihedral(edOmega,&dlist[i])) {
 238       if (bOmega)
 239         print_one(oenv,"corromega",dlist[i].name,"Omega ACF for","C(t)",
 240                   nf/2,time,dih[j]);
 241       j++;
 242     }
 243   }
 244   for(Xi=0; (Xi<maxchi); Xi++) {
 245     sprintf(name1, "corrchi%d", Xi+1);
 246     sprintf(name2, "Chi%d ACF for", Xi+1);
 247     for(i=0; (i<nlist); i++) {
 248       if (dlist[i].atm.Cn[Xi+3] != -1) {
 249         if (bChi)
 250           print_one(oenv,name1,dlist[i].name,name2,"C(t)",nf/2,time,dih[j]);
 251         j++;
 252       }
 253     }
 254   }
 255   fprintf(stderr,"\n");
 256 }
 257
 258 static void copy_dih_data(real in[], real out[], int nf, gmx_bool bLEAVE)
 259 {
 260   /* if bLEAVE, do nothing to data in copying to out
 261    * otherwise multiply by 180/pi to convert rad to deg */
 262   int i ;
 263   real mult ;
 264   if (bLEAVE)
 265     mult = 1 ;
 266   else
 267     mult = (180.0/M_PI);
 268   for (i=0;(i<nf);i++){
 269     out[i]=in[i]*mult ;
 270   }
 271 }
 272
 273 static void dump_em_all(int nlist,t_dlist dlist[],int nf,real time[],
 274                         real **dih,int maxchi,
 275                         gmx_bool bPhi,gmx_bool bPsi,gmx_bool bChi,gmx_bool bOmega, gmx_bool bRAD,
 276                         const output_env_t oenv)
 277 {
 278   char name[256], titlestr[256], ystr[256];
 279   real *data ;
 280   int  i,j,Xi;
 281
 282   snew(data,nf);
 283   if (bRAD)
 284     strcpy(ystr,"Angle (rad)");
 285   else
 286     strcpy(ystr,"Angle (degrees)");
 287
 288   /* Dump em all */
 289   j = 0;
 290   for(i=0; (i<nlist); i++) {
 291       /* grs debug  printf("OK i %d j %d\n", i, j) ; */
 292     if (bPhi) {
 293       copy_dih_data(dih[j],data,nf,bRAD);
 294       print_one(oenv,"phi",dlist[i].name,"\\xf\\f{}",ystr, nf,time,data);
 295     }
 296     j++;
 297   }
 298   for(i=0; (i<nlist); i++) {
 299     if (bPsi) {
 300       copy_dih_data(dih[j],data,nf,bRAD);
 301       print_one(oenv,"psi",dlist[i].name,"\\xy\\f{}",ystr, nf,time,data);
 302     }
 303     j++;
 304   }
 305   for(i=0; (i<nlist); i++)
 306     if (has_dihedral(edOmega,&(dlist[i]))) {
 307       if (bOmega){
 308         copy_dih_data(dih[j],data,nf,bRAD);
 309         print_one(oenv,"omega",dlist[i].name,"\\xw\\f{}",ystr,nf,time,data);
 310       }
 311       j++;
 312     }
 313
 314   for(Xi=0; (Xi<maxchi); Xi++)
 315     for(i=0; (i<nlist); i++)
 316       if (dlist[i].atm.Cn[Xi+3] != -1) {
 317         if (bChi) {
 318           sprintf(name,"chi%d",Xi+1);
 319           sprintf(titlestr,"\\xc\\f{}\\s%d\\N",Xi+1);
 320           copy_dih_data(dih[j],data,nf,bRAD);
 321           print_one(oenv,name,dlist[i].name,titlestr,ystr, nf,time,data);
 322         }
 323         j++;
 324       }
 325   fprintf(stderr,"\n");
 326 }
 327
 328 static void reset_one(real dih[],int nf,real phase)
 329 {
 330   int j;
 331
 332   for(j=0; (j<nf); j++) {
 333     dih[j] += phase;
 334     while(dih[j] < -M_PI)
 335       dih[j] += 2*M_PI;
 336     while(dih[j] >= M_PI)
 337       dih[j] -= 2*M_PI;
 338   }
 339 }
 340
 341 static int reset_em_all(int nlist,t_dlist dlist[],int nf,
 342                          real **dih,int maxchi)
 343 {
 344   int  i,j,Xi;
 345
 346   /* Reset em all */
 347   j=0;
 348   /* Phi */
 349   for(i=0; (i<nlist); i++)
 350   {
 351           if (dlist[i].atm.minC == -1)
 352           {
 353                   reset_one(dih[j++],nf,M_PI);
 354           }
 355           else
 356       {
 357                   reset_one(dih[j++],nf,0);
 358           }
 359   }
 360   /* Psi */
 361   for(i=0; (i<nlist-1); i++)
 362   {
 363           reset_one(dih[j++],nf,0);
 364   }
 365   /* last Psi is faked from O */
 366   reset_one(dih[j++],nf,M_PI);
 367
 368   /* Omega */
 369   for(i=0; (i<nlist); i++)
 370           if (has_dihedral(edOmega,&dlist[i]))
 371                   reset_one(dih[j++],nf,0);
 372   /* Chi 1 thru maxchi */
 373   for(Xi=0; (Xi<maxchi); Xi++)
 374   {
 375           for(i=0; (i<nlist); i++)
 376           {
 377                   if (dlist[i].atm.Cn[Xi+3] != -1)
 378                   {
 379                           reset_one(dih[j],nf,0);
 380                           j++;
 381                   }
 382           }
 383   }
 384   fprintf(stderr,"j after resetting (nr. active dihedrals) = %d\n",j);
 385   return j ;
 386 }
 387
 388 static void histogramming(FILE *log,int nbin,gmx_residuetype_t rt,
 389                           int nf,int maxchi,real **dih,
 390                           int nlist,t_dlist dlist[],
 391                           atom_id index[],
 392                           gmx_bool bPhi,gmx_bool bPsi,gmx_bool bOmega,gmx_bool bChi,
 393                           gmx_bool bNormalize,gmx_bool bSSHisto,const char *ssdump,
 394                           real bfac_max,t_atoms *atoms,
 395                           gmx_bool bDo_jc, const char *fn,
 396                           const output_env_t oenv)
 397 {
 398   /* also gets 3J couplings and order parameters S2 */
 399   t_karplus kkkphi[] = {
 400     { "J_NHa1",    6.51, -1.76,  1.6, -M_PI/3,   0.0,  0.0 },
 401     { "J_NHa2",    6.51, -1.76,  1.6,  M_PI/3,   0.0,  0.0 },
 402     { "J_HaC'",    4.0,   1.1,   0.1,  0.0,      0.0,  0.0 },
 403     { "J_NHCb",    4.7,  -1.5,  -0.2,  M_PI/3,   0.0,  0.0 },
 404     { "J_Ci-1Hai", 4.5,  -1.3,  -1.2,  2*M_PI/3, 0.0,  0.0 }
 405   };
 406   t_karplus kkkpsi[] = {
 407     { "J_HaN",   -0.88, -0.61,-0.27,M_PI/3,  0.0,  0.0 }
 408   };
 409   t_karplus kkkchi1[] = {
 410     { "JHaHb2",       9.5, -1.6, 1.8, -M_PI/3, 0,  0.0 },
 411     { "JHaHb3",       9.5, -1.6, 1.8, 0, 0,  0.0 }
 412   };
 413 #define NKKKPHI asize(kkkphi)
 414 #define NKKKPSI asize(kkkpsi)
 415 #define NKKKCHI asize(kkkchi1)
 416 #define NJC (NKKKPHI+NKKKPSI+NKKKCHI)
 417
 418   FILE    *fp,*ssfp[3]={NULL,NULL,NULL};
 419   const char *sss[3] = { "sheet", "helix", "coil" };
 420   real    S2;
 421   real    *normhisto;
 422   real    **Jc,**Jcsig;
 423   int     ****his_aa_ss=NULL;
 424   int     ***his_aa,**his_aa1,*histmp;
 425   int     i,j,k,m,n,nn,Dih,nres,hindex,angle;
 426   gmx_bool    bBfac,bOccup;
 427   char    hisfile[256],hhisfile[256],sshisfile[256],title[256],*ss_str=NULL;
 428   char **leg;
 429   const char *residue_name;
 430   int     rt_size;
 431
 432   rt_size = gmx_residuetype_get_size(rt);
 433   if (bSSHisto) {
 434     fp = ffopen(ssdump,"r");
 435     if(1 != fscanf(fp,"%d",&nres))
 436     {
 437       gmx_fatal(FARGS,"Error reading from file %s",ssdump);
 438     }
 439
 440     snew(ss_str,nres+1);
 441     if( 1 != fscanf(fp,"%s",ss_str))
 442     {
 443       gmx_fatal(FARGS,"Error reading from file %s",ssdump);
 444     }
 445
 446     ffclose(fp);
 447     /* Four dimensional array... Very cool */
 448     snew(his_aa_ss,3);
 449     for(i=0; (i<3); i++) {
 450       snew(his_aa_ss[i],rt_size+1);
 451       for(j=0; (j<=rt_size); j++) {
 452         snew(his_aa_ss[i][j],edMax);
 453         for(Dih=0; (Dih<edMax); Dih++)
 454           snew(his_aa_ss[i][j][Dih],nbin+1);
 455       }
 456     }
 457   }
 458   snew(his_aa,edMax);
 459   for(Dih=0; (Dih<edMax); Dih++) {
 460     snew(his_aa[Dih],rt_size+1);
 461     for(i=0; (i<=rt_size); i++) {
 462       snew(his_aa[Dih][i],nbin+1);
 463     }
 464   }
 465   snew(histmp,nbin);
 466
 467   snew(Jc,nlist);
 468   snew(Jcsig,nlist);
 469   for(i=0; (i<nlist); i++) {
 470     snew(Jc[i],NJC);
 471     snew(Jcsig[i],NJC);
 472   }
 473
 474   j=0;
 475   n=0;
 476   for (Dih=0; (Dih<NONCHI+maxchi); Dih++) {
 477     for(i=0; (i<nlist); i++) {
 478       if (((Dih  < edOmega) ) ||
 479           ((Dih == edOmega) && (has_dihedral(edOmega,&(dlist[i])))) ||
 480           ((Dih  > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1))) {
 481         make_histo(log,nf,dih[j],nbin,histmp,-M_PI,M_PI);
 482
 483         if (bSSHisto) {
 484           /* Assume there is only one structure, the first.
 485            * Compute index in histogram.
 486            */
 487           /* Check the atoms to see whether their B-factors are low enough
 488            * Check atoms to see their occupancy is 1.
 489            */
 490           bBfac = bOccup = TRUE;
 491           for(nn=0; (nn<4); nn++,n++) {
 492             bBfac  = bBfac  && (atoms->pdbinfo[index[n]].bfac <= bfac_max);
 493             bOccup = bOccup && (atoms->pdbinfo[index[n]].occup == 1);
 494           }
 495           if (bOccup && ((bfac_max <= 0) || ((bfac_max > 0) && bBfac))) {
 496             hindex = ((dih[j][0]+M_PI)*nbin)/(2*M_PI);
 497             range_check(hindex,0,nbin);
 498
 499             /* Assign dihedral to either of the structure determined
 500              * histograms
 501              */
 502             switch(ss_str[dlist[i].resnr]) {
 503             case 'E':
 504               his_aa_ss[0][dlist[i].index][Dih][hindex]++;
 505               break;
 506             case 'H':
 507               his_aa_ss[1][dlist[i].index][Dih][hindex]++;
 508               break;
 509             default:
 510               his_aa_ss[2][dlist[i].index][Dih][hindex]++;
 511               break;
 512             }
 513           }
 514           else if (debug)
 515             fprintf(debug,"Res. %d has imcomplete occupancy or bfacs > %g\n",
 516                     dlist[i].resnr,bfac_max);
 517         }
 518         else
 519           n += 4;
 520
 521         switch (Dih) {
 522         case edPhi:
 523           calc_distribution_props(nbin,histmp,-M_PI,NKKKPHI,kkkphi,&S2);
 524
 525           for(m=0; (m<NKKKPHI); m++) {
 526             Jc[i][m]    = kkkphi[m].Jc;
 527             Jcsig[i][m] = kkkphi[m].Jcsig;
 528           }
 529           break;
 530         case edPsi:
 531           calc_distribution_props(nbin,histmp,-M_PI,NKKKPSI,kkkpsi,&S2);
 532
 533           for(m=0; (m<NKKKPSI); m++) {
 534             Jc[i][NKKKPHI+m]    = kkkpsi[m].Jc;
 535             Jcsig[i][NKKKPHI+m] = kkkpsi[m].Jcsig;
 536           }
 537           break;
 538         case edChi1:
 539           calc_distribution_props(nbin,histmp,-M_PI,NKKKCHI,kkkchi1,&S2);
 540           for(m=0; (m<NKKKCHI); m++) {
 541             Jc[i][NKKKPHI+NKKKPSI+m]    = kkkchi1[m].Jc;
 542             Jcsig[i][NKKKPHI+NKKKPSI+m] = kkkchi1[m].Jcsig;
 543           }
 544           break;
 545         default: /* covers edOmega and higher Chis than Chi1 */
 546           calc_distribution_props(nbin,histmp,-M_PI,0,NULL,&S2);
 547           break;
 548         }
 549         dlist[i].S2[Dih]        = S2;
 550
 551         /* Sum distribution per amino acid type as well */
 552         for(k=0; (k<nbin); k++) {
 553           his_aa[Dih][dlist[i].index][k] += histmp[k];
 554           histmp[k] = 0;
 555         }
 556         j++;
 557       } else { /* dihed not defined */
 558         dlist[i].S2[Dih] = 0.0 ;
 559       }
 560     }
 561   }
 562   sfree(histmp);
 563
 564   /* Print out Jcouplings */
 565   fprintf(log,"\n *** J-Couplings from simulation (plus std. dev.) ***\n\n");
 566   fprintf(log,"Residue   ");
 567   for(i=0; (i<NKKKPHI); i++)
 568     fprintf(log,"%7s   SD",kkkphi[i].name);
 569   for(i=0; (i<NKKKPSI); i++)
 570     fprintf(log,"%7s   SD",kkkpsi[i].name);
 571   for(i=0; (i<NKKKCHI); i++)
 572     fprintf(log,"%7s   SD",kkkchi1[i].name);
 573   fprintf(log,"\n");
 574   for(i=0; (i<NJC+1); i++)
 575     fprintf(log,"------------");
 576   fprintf(log,"\n");
 577   for(i=0; (i<nlist); i++) {
 578     fprintf(log,"%-10s",dlist[i].name);
 579     for(j=0; (j<NJC); j++)
 580       fprintf(log,"  %5.2f %4.2f",Jc[i][j],Jcsig[i][j]);
 581     fprintf(log,"\n");
 582   }
 583   fprintf(log,"\n");
 584
 585   /* and to -jc file... */
 586   if (bDo_jc) {
 587     fp=xvgropen(fn,"\\S3\\NJ-Couplings from Karplus Equation","Residue",
 588                 "Coupling",oenv);
 589     snew(leg,NJC);
 590     for(i=0; (i<NKKKPHI); i++){
 591                 leg[i] = strdup(kkkphi[i].name);
 592     }
 593     for(i=0; (i<NKKKPSI); i++){
 594                 leg[i+NKKKPHI]=strdup(kkkpsi[i].name);
 595     }
 596     for(i=0; (i<NKKKCHI); i++){
 597       leg[i+NKKKPHI+NKKKPSI]=strdup(kkkchi1[i].name);
 598     }
 599     xvgr_legend(fp,NJC,(const char**)leg,oenv);
 600     fprintf(fp,"%5s ","#Res.");
 601     for(i=0; (i<NJC); i++)
 602       fprintf(fp,"%10s ",leg[i]);
 603     fprintf(fp,"\n");
 604     for(i=0; (i<nlist); i++) {
 605       fprintf(fp,"%5d ",dlist[i].resnr);
 606       for(j=0; (j<NJC); j++)
 607         fprintf(fp,"  %8.3f",Jc[i][j]);
 608       fprintf(fp,"\n");
 609     }
 610     ffclose(fp);
 611     for(i=0; (i<NJC); i++)
 612       sfree(leg[i]);
 613   }
 614   /* finished -jc stuff */
 615
 616   snew(normhisto,nbin);
 617   for(i=0; (i<rt_size); i++) {
 618     for(Dih=0; (Dih<edMax); Dih++){
 619       /* First check whether something is in there */
 620       for(j=0; (j<nbin); j++)
 621         if (his_aa[Dih][i][j] != 0)
 622           break;
 623       if ((j < nbin) &&
 624           ((bPhi && (Dih==edPhi)) ||
 625            (bPsi && (Dih==edPsi)) ||
 626            (bOmega &&(Dih==edOmega)) ||
 627            (bChi && (Dih>=edChi1)))) {
 628         if (bNormalize)
 629           normalize_histo(nbin,his_aa[Dih][i],(360.0/nbin),normhisto);
 630
 631         residue_name = gmx_residuetype_get_name(rt,i);
 632         switch (Dih) {
 633         case edPhi:
 634           sprintf(hisfile,"histo-phi%s",residue_name);
 635           sprintf(title,"\\xf\\f{} Distribution for %s",residue_name);
 636           break;
 637         case edPsi:
 638           sprintf(hisfile,"histo-psi%s",residue_name);
 639           sprintf(title,"\\xy\\f{} Distribution for %s",residue_name);
 640           break;
 641         case edOmega:
 642           sprintf(hisfile,"histo-omega%s",residue_name);
 643           sprintf(title,"\\xw\\f{} Distribution for %s",residue_name);
 644           break;
 645         default:
 646           sprintf(hisfile,"histo-chi%d%s",Dih-NONCHI+1,residue_name);
 647           sprintf(title,"\\xc\\f{}\\s%d\\N Distribution for %s",
 648                   Dih-NONCHI+1,residue_name);
 649         }
 650         strcpy(hhisfile,hisfile);
 651         strcat(hhisfile,".xvg");
 652         fp=xvgropen(hhisfile,title,"Degrees","",oenv);
 653         fprintf(fp,"@ with g0\n");
 654         xvgr_world(fp,-180,0,180,0.1,oenv);
 655         fprintf(fp,"# this effort to set graph size fails unless you run with -autoscale none or -autoscale y flags\n");
 656         fprintf(fp,"@ xaxis tick on\n");
 657         fprintf(fp,"@ xaxis tick major 90\n");
 658         fprintf(fp,"@ xaxis tick minor 30\n");
 659         fprintf(fp,"@ xaxis ticklabel prec 0\n");
 660         fprintf(fp,"@ yaxis tick off\n");
 661         fprintf(fp,"@ yaxis ticklabel off\n");
 662         fprintf(fp,"@ type xy\n");
 663         if (bSSHisto) {
 664           for(k=0; (k<3); k++) {
 665             sprintf(sshisfile,"%s-%s.xvg",hisfile,sss[k]);
 666             ssfp[k] = ffopen(sshisfile,"w");
 667           }
 668         }
 669         for(j=0; (j<nbin); j++) {
 670           angle = -180 + (360/nbin)*j ;
 671           if (bNormalize)
 672             fprintf(fp,"%5d  %10g\n",angle,normhisto[j]);
 673           else
 674             fprintf(fp,"%5d  %10d\n",angle,his_aa[Dih][i][j]);
 675           if (bSSHisto)
 676             for(k=0; (k<3); k++)
 677               fprintf(ssfp[k],"%5d  %10d\n",angle,
 678                       his_aa_ss[k][i][Dih][j]);
 679         }
 680         fprintf(fp,"&\n");
 681         ffclose(fp);
 682         if (bSSHisto) {
 683           for(k=0; (k<3); k++) {
 684             fprintf(ssfp[k],"&\n");
 685             ffclose(ssfp[k]);
 686           }
 687         }
 688       }
 689     }
 690   }
 691   sfree(normhisto);
 692
 693   if (bSSHisto) {
 694     /* Four dimensional array... Very cool */
 695     for(i=0; (i<3); i++) {
 696       for(j=0; (j<=rt_size); j++) {
 697         for(Dih=0; (Dih<edMax); Dih++)
 698           sfree(his_aa_ss[i][j][Dih]);
 699         sfree(his_aa_ss[i][j]);
 700       }
 701       sfree(his_aa_ss[i]);
 702     }
 703     sfree(his_aa_ss);
 704     sfree(ss_str);
 705   }
 706 }
 707
 708 static FILE *rama_file(const char *fn,const char *title,const char *xaxis,
 709                        const char *yaxis,const output_env_t oenv)
 710 {
 711   FILE *fp;
 712
 713   fp = xvgropen(fn,title,xaxis,yaxis,oenv);
 714   fprintf(fp,"@ with g0\n");
 715   xvgr_world(fp,-180,-180,180,180,oenv);
 716   fprintf(fp,"@ xaxis tick on\n");
 717   fprintf(fp,"@ xaxis tick major 90\n");
 718   fprintf(fp,"@ xaxis tick minor 30\n");
 719   fprintf(fp,"@ xaxis ticklabel prec 0\n");
 720   fprintf(fp,"@ yaxis tick on\n");
 721   fprintf(fp,"@ yaxis tick major 90\n");
 722   fprintf(fp,"@ yaxis tick minor 30\n");
 723   fprintf(fp,"@ yaxis ticklabel prec 0\n");
 724   fprintf(fp,"@    s0 type xy\n");
 725   fprintf(fp,"@    s0 symbol 2\n");
 726   fprintf(fp,"@    s0 symbol size 0.410000\n");
 727   fprintf(fp,"@    s0 symbol fill 1\n");
 728   fprintf(fp,"@    s0 symbol color 1\n");
 729   fprintf(fp,"@    s0 symbol linewidth 1\n");
 730   fprintf(fp,"@    s0 symbol linestyle 1\n");
 731   fprintf(fp,"@    s0 symbol center false\n");
 732   fprintf(fp,"@    s0 symbol char 0\n");
 733   fprintf(fp,"@    s0 skip 0\n");
 734   fprintf(fp,"@    s0 linestyle 0\n");
 735   fprintf(fp,"@    s0 linewidth 1\n");
 736   fprintf(fp,"@ type xy\n");
 737
 738   return fp;
 739 }
 740
 741 static void do_rama(int nf,int nlist,t_dlist dlist[],real **dih,
 742                     gmx_bool bViol,gmx_bool bRamOmega,const output_env_t oenv)
 743 {
 744   FILE *fp,*gp=NULL;
 745   gmx_bool bOm;
 746   char fn[256];
 747   int  i,j,k,Xi1,Xi2,Phi,Psi,Om=0,nlevels;
 748 #define NMAT 120
 749   real **mat=NULL,phi,psi,omega,axis[NMAT],lo,hi;
 750   t_rgb rlo = { 1.0, 0.0, 0.0 };
 751   t_rgb rmid= { 1.0, 1.0, 1.0 };
 752   t_rgb rhi = { 0.0, 0.0, 1.0 };
 753
 754   for(i=0; (i<nlist); i++) {
 755     if ((has_dihedral(edPhi,&(dlist[i]))) &&
 756         (has_dihedral(edPsi,&(dlist[i])))) {
 757       sprintf(fn,"ramaPhiPsi%s.xvg",dlist[i].name);
 758       fp = rama_file(fn,"Ramachandran Plot",
 759                      "\\8f\\4 (deg)","\\8y\\4 (deg)",oenv);
 760       bOm = bRamOmega && has_dihedral(edOmega,&(dlist[i]));
 761       if (bOm) {
 762         Om = dlist[i].j0[edOmega];
 763         snew(mat,NMAT);
 764         for(j=0; (j<NMAT); j++) {
 765           snew(mat[j],NMAT);
 766           axis[j] = -180+(360*j)/NMAT;
 767         }
 768       }
 769       if (bViol) {
 770         sprintf(fn,"violPhiPsi%s.xvg",dlist[i].name);
 771         gp = ffopen(fn,"w");
 772       }
 773       Phi = dlist[i].j0[edPhi];
 774       Psi = dlist[i].j0[edPsi];
 775       for(j=0; (j<nf); j++) {
 776         phi = RAD2DEG*dih[Phi][j];
 777         psi = RAD2DEG*dih[Psi][j];
 778         fprintf(fp,"%10g  %10g\n",phi,psi);
 779         if (bViol)
 780           fprintf(gp,"%d\n",!bAllowed(dih[Phi][j],RAD2DEG*dih[Psi][j]));
 781         if (bOm) {
 782           omega = RAD2DEG*dih[Om][j];
 783           mat[(int)((phi*NMAT)/360)+NMAT/2][(int)((psi*NMAT)/360)+NMAT/2]
 784             += omega;
 785         }
 786       }
 787       if (bViol)
 788         ffclose(gp);
 789       ffclose(fp);
 790       if (bOm) {
 791         sprintf(fn,"ramomega%s.xpm",dlist[i].name);
 792         fp = ffopen(fn,"w");
 793         lo = hi = 0;
 794         for(j=0; (j<NMAT); j++)
 795           for(k=0; (k<NMAT); k++) {
 796             mat[j][k] /= nf;
 797             lo=min(mat[j][k],lo);
 798             hi=max(mat[j][k],hi);
 799           }
 800         /* Symmetrise */
 801         if (fabs(lo) > fabs(hi))
 802           hi = -lo;
 803         else
 804           lo = -hi;
 805         /* Add 180 */
 806         for(j=0; (j<NMAT); j++)
 807           for(k=0; (k<NMAT); k++)
 808             mat[j][k] += 180;
 809         lo += 180;
 810         hi += 180;
 811         nlevels = 20;
 812         write_xpm3(fp,0,"Omega/Ramachandran Plot","Deg","Phi","Psi",
 813                    NMAT,NMAT,axis,axis,mat,lo,180.0,hi,rlo,rmid,rhi,&nlevels);
 814         ffclose(fp);
 815         for(j=0; (j<NMAT); j++)
 816           sfree(mat[j]);
 817         sfree(mat);
 818       }
 819     }
 820     if ((has_dihedral(edChi1,&(dlist[i]))) &&
 821         (has_dihedral(edChi2,&(dlist[i])))) {
 822       sprintf(fn,"ramaX1X2%s.xvg",dlist[i].name);
 823       fp = rama_file(fn,"\\8c\\4\\s1\\N-\\8c\\4\\s2\\N Ramachandran Plot",
 824                      "\\8c\\4\\s1\\N (deg)","\\8c\\4\\s2\\N (deg)",oenv);
 825       Xi1 = dlist[i].j0[edChi1];
 826       Xi2 = dlist[i].j0[edChi2];
 827       for(j=0; (j<nf); j++)
 828         fprintf(fp,"%10g  %10g\n",RAD2DEG*dih[Xi1][j],RAD2DEG*dih[Xi2][j]);
 829       ffclose(fp);
 830     }
 831     else
 832       fprintf(stderr,"No chi1 & chi2 angle for %s\n",dlist[i].name);
 833   }
 834 }
 835
 836
 837 static void print_transitions(const char *fn,int maxchi,int nlist,
 838                               t_dlist dlist[], t_atoms *atoms,rvec x[],
 839                               matrix box, gmx_bool bPhi,gmx_bool bPsi,gmx_bool bChi,real dt,
 840                               const output_env_t oenv)
 841 {
 842   /* based on order_params below */
 843   FILE *fp;
 844   int  nh[edMax];
 845   int  i,Dih,Xi;
 846
 847   /*  must correspond with enum in pp2shift.h:38 */
 848   char *leg[edMax];
 849 #define NLEG asize(leg)
 850
 851   leg[0] = strdup("Phi");
 852   leg[1] = strdup("Psi");
 853   leg[2] = strdup("Omega");
 854   leg[3] = strdup("Chi1");
 855   leg[4] = strdup("Chi2");
 856   leg[5] = strdup("Chi3");
 857   leg[6] = strdup("Chi4");
 858   leg[7] = strdup("Chi5");
 859   leg[8] = strdup("Chi6");
 860
 861   /* Print order parameters */
 862   fp=xvgropen(fn,"Dihedral Rotamer Transitions","Residue","Transitions/ns",
 863                 oenv);
 864   xvgr_legend(fp,NONCHI+maxchi,(const char**)leg,oenv);
 865
 866   for (Dih=0; (Dih<edMax); Dih++)
 867     nh[Dih]=0;
 868
 869   fprintf(fp,"%5s ","#Res.");
 870   fprintf(fp,"%10s %10s %10s ",leg[edPhi],leg[edPsi],leg[edOmega]);
 871   for (Xi=0; Xi<maxchi; Xi++)
 872     fprintf(fp,"%10s ",leg[NONCHI+Xi]);
 873   fprintf(fp,"\n");
 874
 875   for(i=0; (i<nlist); i++) {
 876     fprintf(fp,"%5d ",dlist[i].resnr);
 877     for (Dih=0; (Dih<NONCHI+maxchi); Dih++)
 878       fprintf(fp,"%10.3f ",dlist[i].ntr[Dih]/dt);
 879     /* fprintf(fp,"%12s\n",dlist[i].name);  this confuses xmgrace */
 880     fprintf(fp,"\n");
 881   }
 882   ffclose(fp);
 883 }
 884
 885 static void order_params(FILE *log,
 886                          const char *fn,int maxchi,int nlist,t_dlist dlist[],
 887                          const char *pdbfn,real bfac_init,
 888                          t_atoms *atoms,rvec x[],int ePBC,matrix box,
 889                          gmx_bool bPhi,gmx_bool bPsi,gmx_bool bChi,const output_env_t oenv)
 890 {
 891   FILE *fp;
 892   int  nh[edMax];
 893   char buf[STRLEN];
 894   int  i,Dih,Xi;
 895   real S2Max, S2Min;
 896
 897   /* except for S2Min/Max, must correspond with enum in pp2shift.h:38 */
 898   const char *const_leg[2+edMax]= { "S2Min","S2Max","Phi","Psi","Omega",
 899                                     "Chi1", "Chi2", "Chi3", "Chi4", "Chi5",
 900                                     "Chi6" };
 901 #define NLEG asize(leg)
 902
 903   char *leg[2+edMax];
 904
 905   for(i=0;i<NLEG;i++)
 906     leg[i]=strdup(const_leg[i]);
 907
 908   /* Print order parameters */
 909   fp=xvgropen(fn,"Dihedral Order Parameters","Residue","S2",oenv);
 910   xvgr_legend(fp,2+NONCHI+maxchi,const_leg,oenv);
 911
 912   for (Dih=0; (Dih<edMax); Dih++)
 913     nh[Dih]=0;
 914
 915   fprintf(fp,"%5s ","#Res.");
 916   fprintf(fp,"%10s %10s ",leg[0],leg[1]);
 917   fprintf(fp,"%10s %10s %10s ",leg[2+edPhi],leg[2+edPsi],leg[2+edOmega]);
 918   for (Xi=0; Xi<maxchi; Xi++)
 919     fprintf(fp,"%10s ",leg[2+NONCHI+Xi]);
 920   fprintf(fp,"\n");
 921
 922   for(i=0; (i<nlist); i++) {
 923     S2Max=-10;
 924     S2Min=10;
 925     for (Dih=0; (Dih<NONCHI+maxchi); Dih++) {
 926       if (dlist[i].S2[Dih]!=0) {
 927         if (dlist[i].S2[Dih] > S2Max)
 928           S2Max=dlist[i].S2[Dih];
 929         if (dlist[i].S2[Dih] < S2Min)
 930           S2Min=dlist[i].S2[Dih];
 931       }
 932       if (dlist[i].S2[Dih] > 0.8)
 933         nh[Dih]++;
 934     }
 935     fprintf(fp,"%5d ",dlist[i].resnr);
 936     fprintf(fp,"%10.3f %10.3f ",S2Min,S2Max);
 937     for (Dih=0; (Dih<NONCHI+maxchi); Dih++)
 938       fprintf(fp,"%10.3f ",dlist[i].S2[Dih]);
 939     fprintf(fp,"\n");
 940     /* fprintf(fp,"%12s\n",dlist[i].name);  this confuses xmgrace */
 941   }
 942   ffclose(fp);
 943
 944   if (NULL != pdbfn) {
 945     real x0,y0,z0;
 946
 947     if (NULL == atoms->pdbinfo)
 948       snew(atoms->pdbinfo,atoms->nr);
 949     for(i=0; (i<atoms->nr); i++)
 950       atoms->pdbinfo[i].bfac=bfac_init;
 951
 952     for(i=0; (i<nlist); i++) {
 953       atoms->pdbinfo[dlist[i].atm.N].bfac=-dlist[i].S2[0];/* Phi */
 954       atoms->pdbinfo[dlist[i].atm.H].bfac=-dlist[i].S2[0];/* Phi */
 955       atoms->pdbinfo[dlist[i].atm.C].bfac=-dlist[i].S2[1];/* Psi */
 956       atoms->pdbinfo[dlist[i].atm.O].bfac=-dlist[i].S2[1];/* Psi */
 957       for (Xi=0; (Xi<maxchi); Xi++) {           /* Chi's */
 958         if (dlist[i].atm.Cn[Xi+3]!=-1) {
 959           atoms->pdbinfo[dlist[i].atm.Cn[Xi+1]].bfac=-dlist[i].S2[NONCHI+Xi];
 960         }
 961       }
 962     }
 963
 964     fp=ffopen(pdbfn,"w");
 965     fprintf(fp,"REMARK generated by g_chi\n");
 966     fprintf(fp,"REMARK "
 967             "B-factor field contains negative of dihedral order parameters\n");
 968     write_pdbfile(fp,NULL,atoms,x,ePBC,box,' ',0,NULL,TRUE);
 969     x0=y0=z0=1000.0;
 970     for (i=0; (i<atoms->nr); i++) {
 971       x0 = min(x0, x[i][XX]);
 972       y0 = min(y0, x[i][YY]);
 973       z0 = min(z0, x[i][ZZ]);
 974     }
 975     x0*=10.0;/* nm -> angstrom */
 976     y0*=10.0;/* nm -> angstrom */
 977     z0*=10.0;/* nm -> angstrom */
 978     sprintf(buf,"%s%%6.f%%6.2f\n",pdbformat);
 979     for (i=0; (i<10); i++) {
 980       fprintf(fp,buf,"ATOM  ", atoms->nr+1+i, "CA", "LEG",' ',
 981               atoms->nres+1, ' ',x0, y0, z0+(1.2*i), 0.0, -0.1*i);
 982     }
 983     ffclose(fp);
 984   }
 985
 986   fprintf(log,"Dihedrals with S2 > 0.8\n");
 987   fprintf(log,"Dihedral: ");
 988   if (bPhi) fprintf(log," Phi  ");
 989   if (bPsi) fprintf(log," Psi ");
 990   if (bChi)
 991     for(Xi=0; (Xi<maxchi); Xi++)
 992       fprintf(log," %s ", leg[2+NONCHI+Xi]);
 993   fprintf(log,"\nNumber:   ");
 994   if (bPhi) fprintf(log,"%4d  ",nh[0]);
 995   if (bPsi) fprintf(log,"%4d  ",nh[1]);
 996   if (bChi)
 997     for(Xi=0; (Xi<maxchi); Xi++)
 998       fprintf(log,"%4d  ",nh[NONCHI+Xi]);
 999   fprintf(log,"\n");
1000
1001   for(i=0;i<NLEG;i++)
1002     sfree(leg[i]);
1003
1004 }
1005
1006 int gmx_chi(int argc,char *argv[])
1007 {
1008   const char *desc[] = {
1009     "g_chi computes phi, psi, omega and chi dihedrals for all your ",
1010     "amino acid backbone and sidechains.",
1011     "It can compute dihedral angle as a function of time, and as",
1012     "histogram distributions.",
1013     "The distributions (histo-(dihedral)(RESIDUE).xvg) are cumulative over all residues of each type.[PAR]",
1014     "If option [TT]-corr[tt] is given, the program will",
1015     "calculate dihedral autocorrelation functions. The function used",
1016     "is C(t) = < cos(chi(tau)) cos(chi(tau+t)) >. The use of cosines",
1017     "rather than angles themselves, resolves the problem of periodicity.",
1018     "(Van der Spoel & Berendsen (1997), Biophys. J. 72, 2032-2041).",
1019     "Separate files for each dihedral of each residue",
1020     "(corr(dihedral)(RESIDUE)(nresnr).xvg) are output, as well as a",
1021     "file containing the information for all residues (argument of [TT]-corr[tt]).[PAR]",
1022     "With option [TT]-all[tt], the angles themselves as a function of time for",
1023     "each residue are printed to separate files (dihedral)(RESIDUE)(nresnr).xvg.",
1024     "These can be in radians or degrees.[PAR]",
1025     "A log file (argument [TT]-g[tt]) is also written. This contains [BR]",
1026     "(a) information about the number of residues of each type.[BR]",
1027     "(b) The NMR 3J coupling constants from the Karplus equation.[BR]",
1028     "(c) a table for each residue of the number of transitions between ",
1029     "rotamers per nanosecond,  and the order parameter S2 of each dihedral.[BR]",
1030     "(d) a table for each residue of the rotamer occupancy.[BR]",
1031     "All rotamers are taken as 3-fold, except for omegas and chi-dihedrals",
1032     "to planar groups (i.e. chi2 of aromatics Asp and Asn, chi3 of Glu",
1033     "and Gln, and chi4 of Arg), which are 2-fold. \"rotamer 0\" means ",
1034     "that the dihedral was not in the core region of each rotamer. ",
1035     "The width of the core region can be set with [TT]-core_rotamer[tt][PAR]",
1036
1037     "The S2 order parameters are also output to an xvg file",
1038     "(argument [TT]-o[tt] ) and optionally as a pdb file with",
1039     "the S2 values as B-factor (argument [TT]-p[tt]). ",
1040     "The total number of rotamer transitions per timestep",
1041     "(argument [TT]-ot[tt]), the number of transitions per rotamer",
1042     "(argument [TT]-rt[tt]), and the 3J couplings (argument [TT]-jc[tt]), ",
1043     "can also be written to .xvg files.[PAR]",
1044
1045     "If [TT]-chi_prod[tt] is set (and maxchi > 0), cumulative rotamers, e.g.",
1046     "1+9(chi1-1)+3(chi2-1)+(chi3-1) (if the residue has three 3-fold ",
1047     "dihedrals and maxchi >= 3)",
1048     "are calculated. As before, if any dihedral is not in the core region,",
1049     "the rotamer is taken to be 0. The occupancies of these cumulative ",
1050     "rotamers (starting with rotamer 0) are written to the file",
1051     "that is the argument of [TT]-cp[tt], and if the [TT]-all[tt] flag",
1052     "is given, the rotamers as functions of time",
1053     "are written to chiproduct(RESIDUE)(nresnr).xvg ",
1054     "and their occupancies to histo-chiproduct(RESIDUE)(nresnr).xvg.[PAR]",
1055
1056     "The option [TT]-r[tt] generates a contour plot of the average omega angle",
1057     "as a function of the phi and psi angles, that is, in a Ramachandran plot",
1058     "the average omega angle is plotted using color coding.",
1059
1060   };
1061
1062   const char *bugs[] = {
1063     "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.",
1064     "Phi and psi dihedrals are calculated in a non-standard way, using H-N-CA-C for phi instead of C(-)-N-CA-C, and N-CA-C-O for psi instead of N-CA-C-N(+). This causes (usually small) discrepancies with the output of other tools like g_rama.",
1065     "-r0 option does not work properly",
1066     "Rotamers with multiplicity 2 are printed in chi.log as if they had multiplicity 3, with the 3rd (g(+)) always having probability 0"
1067   };
1068
1069   /* defaults */
1070   static int  r0=1,ndeg=1,maxchi=2;
1071   static gmx_bool bAll=FALSE;
1072   static gmx_bool bPhi=FALSE,bPsi=FALSE,bOmega=FALSE;
1073   static real bfac_init=-1.0,bfac_max=0;
1074   static const char *maxchistr[] = { NULL, "0", "1", "2", "3",  "4", "5", "6", NULL };
1075   static gmx_bool bRama=FALSE,bShift=FALSE,bViol=FALSE,bRamOmega=FALSE;
1076   static gmx_bool bNormHisto=TRUE,bChiProduct=FALSE,bHChi=FALSE,bRAD=FALSE,bPBC=TRUE;
1077   static real core_frac=0.5 ;
1078   t_pargs pa[] = {
1079     { "-r0",  FALSE, etINT, {&r0},
1080       "starting residue" },
1081     { "-phi",  FALSE, etBOOL, {&bPhi},
1082       "Output for Phi dihedral angles" },
1083     { "-psi",  FALSE, etBOOL, {&bPsi},
1084       "Output for Psi dihedral angles" },
1085     { "-omega",FALSE, etBOOL, {&bOmega},
1086       "Output for Omega dihedrals (peptide bonds)" },
1087     { "-rama", FALSE, etBOOL, {&bRama},
1088       "Generate Phi/Psi and Chi1/Chi2 ramachandran plots" },
1089     { "-viol", FALSE, etBOOL, {&bViol},
1090       "Write a file that gives 0 or 1 for violated Ramachandran angles" },
1091     { "-periodic", FALSE, etBOOL, {&bPBC},
1092       "Print dihedral angles modulo 360 degrees" },
1093     { "-all",  FALSE, etBOOL, {&bAll},
1094       "Output separate files for every dihedral." },
1095     { "-rad",  FALSE, etBOOL, {&bRAD},
1096       "in angle vs time files, use radians rather than degrees."},
1097     { "-shift", FALSE, etBOOL, {&bShift},
1098         "Compute chemical shifts from Phi/Psi angles" },
1099     { "-binwidth", FALSE, etINT, {&ndeg},
1100       "bin width for histograms (degrees)" },
1101     { "-core_rotamer", FALSE, etREAL, {&core_frac},
1102       "only the central -core_rotamer*(360/multiplicity) belongs to each rotamer (the rest is assigned to rotamer 0)" },
1103     { "-maxchi", FALSE, etENUM, {maxchistr},
1104       "calculate first ndih Chi dihedrals" },
1105     { "-normhisto", FALSE, etBOOL, {&bNormHisto},
1106       "Normalize histograms" },
1107     { "-ramomega",FALSE,etBOOL, {&bRamOmega},
1108       "compute average omega as a function of phi/psi and plot it in an xpm plot" },
1109     { "-bfact", FALSE, etREAL, {&bfac_init},
1110       "B-factor value for pdb file for atoms with no calculated dihedral order parameter"},
1111     { "-chi_prod",FALSE,etBOOL, {&bChiProduct},
1112       "compute a single cumulative rotamer for each residue"},
1113     { "-HChi",FALSE,etBOOL, {&bHChi},
1114       "Include dihedrals to sidechain hydrogens"},
1115     { "-bmax",  FALSE, etREAL, {&bfac_max},
1116       "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. -bmax <= 0 means no limit." }
1117   };
1118
1119   FILE       *log;
1120   int        natoms,nlist,idum,nbin;
1121   t_atoms    atoms;
1122   rvec       *x;
1123   int        ePBC;
1124   matrix     box;
1125   char       title[256],grpname[256];
1126   t_dlist    *dlist;
1127   gmx_bool       bChi,bCorr,bSSHisto;
1128   gmx_bool       bDo_rt, bDo_oh, bDo_ot, bDo_jc ;
1129   real       dt=0, traj_t_ns;
1130   output_env_t oenv;
1131   gmx_residuetype_t rt;
1132
1133   atom_id    isize,*index;
1134   int        ndih,nactdih,nf;
1135   real       **dih,*trans_frac,*aver_angle,*time;
1136   int        i,j,**chi_lookup,*xity;
1137
1138   t_filenm  fnm[] = {
1139     { efSTX, "-s",  NULL,     ffREAD  },
1140     { efTRX, "-f",  NULL,     ffREAD  },
1141     { efXVG, "-o",  "order",  ffWRITE },
1142     { efPDB, "-p",  "order",  ffOPTWR },
1143     { efDAT, "-ss", "ssdump", ffOPTRD },
1144     { efXVG, "-jc", "Jcoupling", ffWRITE },
1145     { efXVG, "-corr",  "dihcorr",ffOPTWR },
1146     { efLOG, "-g",  "chi",    ffWRITE },
1147     /* add two more arguments copying from g_angle */
1148     { efXVG, "-ot", "dihtrans", ffOPTWR },
1149     { efXVG, "-oh", "trhisto",  ffOPTWR },
1150     { efXVG, "-rt", "restrans",  ffOPTWR },
1151     { efXVG, "-cp", "chiprodhisto",  ffOPTWR }
1152   };
1153 #define NFILE asize(fnm)
1154   int     npargs;
1155   t_pargs *ppa;
1156
1157   CopyRight(stderr,argv[0]);
1158   npargs = asize(pa);
1159   ppa    = add_acf_pargs(&npargs,pa);
1160   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
1161                     NFILE,fnm,npargs,ppa,asize(desc),desc,asize(bugs),bugs,
1162                     &oenv);
1163
1164   /* Handle result from enumerated type */
1165   sscanf(maxchistr[0],"%d",&maxchi);
1166   bChi = (maxchi > 0);
1167
1168   log=ffopen(ftp2fn(efLOG,NFILE,fnm),"w");
1169
1170   if (bRamOmega) {
1171     bOmega = TRUE;
1172     bPhi   = TRUE;
1173     bPsi   = TRUE;
1174   }
1175
1176   /* set some options */
1177   bDo_rt=(opt2bSet("-rt",NFILE,fnm));
1178   bDo_oh=(opt2bSet("-oh",NFILE,fnm));
1179   bDo_ot=(opt2bSet("-ot",NFILE,fnm));
1180   bDo_jc=(opt2bSet("-jc",NFILE,fnm));
1181   bCorr=(opt2bSet("-corr",NFILE,fnm));
1182   if (bCorr)
1183     fprintf(stderr,"Will calculate autocorrelation\n");
1184
1185   if (core_frac > 1.0 ) {
1186     fprintf(stderr, "core_rotamer fraction > 1.0 ; will use 1.0\n");
1187     core_frac=1.0 ;
1188   }
1189   if (core_frac < 0.0 ) {
1190     fprintf(stderr, "core_rotamer fraction < 0.0 ; will use 0.0\n");
1191     core_frac=0.0 ;
1192   }
1193
1194   if (maxchi > MAXCHI) {
1195     fprintf(stderr,
1196             "Will only calculate first %d Chi dihedrals in stead of %d.\n",
1197             MAXCHI, maxchi);
1198     maxchi=MAXCHI;
1199   }
1200   bSSHisto = ftp2bSet(efDAT,NFILE,fnm);
1201   nbin = 360/ndeg ;
1202
1203   /* Find the chi angles using atoms struct and a list of amino acids */
1204   get_stx_coordnum(ftp2fn(efSTX,NFILE,fnm),&natoms);
1205   init_t_atoms(&atoms,natoms,TRUE);
1206   snew(x,natoms);
1207   read_stx_conf(ftp2fn(efSTX,NFILE,fnm),title,&atoms,x,NULL,&ePBC,box);
1208   fprintf(log,"Title: %s\n",title);
1209
1210   gmx_residuetype_init(&rt);
1211   dlist=mk_dlist(log,&atoms,&nlist,bPhi,bPsi,bChi,bHChi,maxchi,r0,rt);
1212   fprintf(stderr,"%d residues with dihedrals found\n", nlist);
1213
1214   if (nlist == 0)
1215     gmx_fatal(FARGS,"No dihedrals in your structure!\n");
1216
1217   /* Make a linear index for reading all. */
1218   index=make_chi_ind(nlist,dlist,&ndih);
1219   isize=4*ndih;
1220   fprintf(stderr,"%d dihedrals found\n", ndih);
1221
1222   snew(dih,ndih);
1223
1224   /* COMPUTE ALL DIHEDRALS! */
1225   read_ang_dih(ftp2fn(efTRX,NFILE,fnm),FALSE,TRUE,FALSE,bPBC,1,&idum,
1226                &nf,&time,isize,index,&trans_frac,&aver_angle,dih,oenv);
1227
1228   dt=(time[nf-1]-time[0])/(nf-1); /* might want this for corr or n. transit*/
1229   if (bCorr)
1230   {
1231           if (nf < 2)
1232           {
1233                   gmx_fatal(FARGS,"Need at least 2 frames for correlation");
1234           }
1235   }
1236
1237   /* put angles in -M_PI to M_PI ! and correct phase factor for phi and psi
1238   * pass nactdih instead of ndih to low_ana_dih_trans and get_chi_product_traj
1239   * to prevent accessing off end of arrays when maxchi < 5 or 6. */
1240   nactdih = reset_em_all(nlist,dlist,nf,dih,maxchi);
1241
1242   if (bAll)
1243     dump_em_all(nlist,dlist,nf,time,dih,maxchi,bPhi,bPsi,bChi,bOmega,bRAD,oenv);
1244
1245   /* Histogramming & J coupling constants & calc of S2 order params */
1246   histogramming(log,nbin,rt,nf,maxchi,dih,nlist,dlist,index,
1247                 bPhi,bPsi,bOmega,bChi,
1248                 bNormHisto,bSSHisto,ftp2fn(efDAT,NFILE,fnm),bfac_max,&atoms,
1249                 bDo_jc,opt2fn("-jc",NFILE,fnm),oenv);
1250
1251   /* transitions
1252    *
1253    * added multiplicity */
1254
1255   snew(xity,ndih) ;
1256   mk_multiplicity_lookup(xity, maxchi, dih, nlist, dlist,ndih);
1257
1258   strcpy(grpname, "All residues, ");
1259   if(bPhi)
1260     strcat(grpname, "Phi ");
1261   if(bPsi)
1262     strcat(grpname, "Psi ");
1263   if(bOmega)
1264     strcat(grpname, "Omega ");
1265   if(bChi){
1266     strcat(grpname, "Chi 1-") ;
1267     sprintf(grpname + strlen(grpname), "%i", maxchi);
1268   }
1269
1270
1271   low_ana_dih_trans(bDo_ot, opt2fn("-ot",NFILE,fnm),
1272                     bDo_oh, opt2fn("-oh",NFILE,fnm),maxchi,
1273                     dih, nlist, dlist, nf, nactdih, grpname, xity,
1274                     *time,  dt, FALSE, core_frac,oenv) ;
1275
1276   /* Order parameters */
1277   order_params(log,opt2fn("-o",NFILE,fnm),maxchi,nlist,dlist,
1278                ftp2fn_null(efPDB,NFILE,fnm),bfac_init,
1279                &atoms,x,ePBC,box,bPhi,bPsi,bChi,oenv);
1280
1281   /* Print ramachandran maps! */
1282   if (bRama)
1283     do_rama(nf,nlist,dlist,dih,bViol,bRamOmega,oenv);
1284
1285   if (bShift)
1286     do_pp2shifts(log,nf,nlist,dlist,dih);
1287
1288   /* rprint S^2, transitions, and rotamer occupancies to log */
1289   traj_t_ns = 0.001 * (time[nf-1]-time[0]) ;
1290   pr_dlist(log,nlist,dlist,traj_t_ns,edPrintST,bPhi,bPsi,bChi,bOmega,maxchi);
1291   pr_dlist(log,nlist,dlist,traj_t_ns,edPrintRO,bPhi,bPsi,bChi,bOmega,maxchi);
1292   ffclose(log);
1293   /* transitions to xvg */
1294   if (bDo_rt)
1295     print_transitions(opt2fn("-rt",NFILE,fnm),maxchi,nlist,dlist,
1296                       &atoms,x,box,bPhi,bPsi,bChi,traj_t_ns,oenv);
1297
1298   /* chi_product trajectories (ie one "rotamer number" for each residue) */
1299   if (bChiProduct && bChi ) {
1300     snew(chi_lookup,nlist) ;
1301     for (i=0;i<nlist;i++)
1302       snew(chi_lookup[i], maxchi) ;
1303     mk_chi_lookup(chi_lookup, maxchi, dih, nlist, dlist);
1304
1305     get_chi_product_traj(dih,nf,nactdih,nlist,
1306                          maxchi,dlist,time,chi_lookup,xity,
1307                          FALSE,bNormHisto, core_frac,bAll,
1308                          opt2fn("-cp",NFILE,fnm),oenv);
1309
1310     for (i=0;i<nlist;i++)
1311       sfree(chi_lookup[i]);
1312   }
1313
1314   /* Correlation comes last because it fucks up the angles */
1315   if (bCorr)
1316     do_dihcorr(opt2fn("-corr",NFILE,fnm),nf,ndih,dih,dt,nlist,dlist,time,
1317                maxchi,bPhi,bPsi,bChi,bOmega,oenv);
1318
1319
1320   do_view(oenv,opt2fn("-o",NFILE,fnm),"-nxy");
1321   do_view(oenv,opt2fn("-jc",NFILE,fnm),"-nxy");
1322   if (bCorr)
1323     do_view(oenv,opt2fn("-corr",NFILE,fnm),"-nxy");
1324
1325   gmx_residuetype_destroy(rt);
1326
1327   thanx(stderr);
1328
1329   return 0;
1330 }