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, int naa,char **aa,
 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
 430   if (bSSHisto) {
 431     fp = ffopen(ssdump,"r");
 432     if(1 != fscanf(fp,"%d",&nres))
 433     {
 434       gmx_fatal(FARGS,"Error reading from file %s",ssdump);
 435     }
 436
 437     snew(ss_str,nres+1);
 438     if( 1 != fscanf(fp,"%s",ss_str))
 439     {
 440       gmx_fatal(FARGS,"Error reading from file %s",ssdump);
 441     }
 442
 443     ffclose(fp);
 444     /* Four dimensional array... Very cool */
 445     snew(his_aa_ss,3);
 446     for(i=0; (i<3); i++) {
 447       snew(his_aa_ss[i],naa+1);
 448       for(j=0; (j<=naa); j++) {
 449         snew(his_aa_ss[i][j],edMax);
 450         for(Dih=0; (Dih<edMax); Dih++)
 451           snew(his_aa_ss[i][j][Dih],nbin+1);
 452       }
 453     }
 454   }
 455   snew(his_aa,edMax);
 456   for(Dih=0; (Dih<edMax); Dih++) {
 457     snew(his_aa[Dih],naa+1);
 458     for(i=0; (i<=naa); i++) {
 459       snew(his_aa[Dih][i],nbin+1);
 460     }
 461   }
 462   snew(histmp,nbin);
 463
 464   snew(Jc,nlist);
 465   snew(Jcsig,nlist);
 466   for(i=0; (i<nlist); i++) {
 467     snew(Jc[i],NJC);
 468     snew(Jcsig[i],NJC);
 469   }
 470
 471   j=0;
 472   n=0;
 473   for (Dih=0; (Dih<NONCHI+maxchi); Dih++) {
 474     for(i=0; (i<nlist); i++) {
 475       if (((Dih  < edOmega) ) ||
 476           ((Dih == edOmega) && (has_dihedral(edOmega,&(dlist[i])))) ||
 477           ((Dih  > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1))) {
 478         make_histo(log,nf,dih[j],nbin,histmp,-M_PI,M_PI);
 479
 480         if (bSSHisto) {
 481           /* Assume there is only one structure, the first.
 482            * Compute index in histogram.
 483            */
 484           /* Check the atoms to see whether their B-factors are low enough
 485            * Check atoms to see their occupancy is 1.
 486            */
 487           bBfac = bOccup = TRUE;
 488           for(nn=0; (nn<4); nn++,n++) {
 489             bBfac  = bBfac  && (atoms->pdbinfo[index[n]].bfac <= bfac_max);
 490             bOccup = bOccup && (atoms->pdbinfo[index[n]].occup == 1);
 491           }
 492           if (bOccup && ((bfac_max <= 0) || ((bfac_max > 0) && bBfac))) {
 493             hindex = ((dih[j][0]+M_PI)*nbin)/(2*M_PI);
 494             range_check(hindex,0,nbin);
 495
 496             /* Assign dihedral to either of the structure determined
 497              * histograms
 498              */
 499             switch(ss_str[dlist[i].resnr]) {
 500             case 'E':
 501               his_aa_ss[0][dlist[i].index][Dih][hindex]++;
 502               break;
 503             case 'H':
 504               his_aa_ss[1][dlist[i].index][Dih][hindex]++;
 505               break;
 506             default:
 507               his_aa_ss[2][dlist[i].index][Dih][hindex]++;
 508               break;
 509             }
 510           }
 511           else if (debug)
 512             fprintf(debug,"Res. %d has imcomplete occupancy or bfacs > %g\n",
 513                     dlist[i].resnr,bfac_max);
 514         }
 515         else
 516           n += 4;
 517
 518         switch (Dih) {
 519         case edPhi:
 520           calc_distribution_props(nbin,histmp,-M_PI,NKKKPHI,kkkphi,&S2);
 521
 522           for(m=0; (m<NKKKPHI); m++) {
 523             Jc[i][m]    = kkkphi[m].Jc;
 524             Jcsig[i][m] = kkkphi[m].Jcsig;
 525           }
 526           break;
 527         case edPsi:
 528           calc_distribution_props(nbin,histmp,-M_PI,NKKKPSI,kkkpsi,&S2);
 529
 530           for(m=0; (m<NKKKPSI); m++) {
 531             Jc[i][NKKKPHI+m]    = kkkpsi[m].Jc;
 532             Jcsig[i][NKKKPHI+m] = kkkpsi[m].Jcsig;
 533           }
 534           break;
 535         case edChi1:
 536           calc_distribution_props(nbin,histmp,-M_PI,NKKKCHI,kkkchi1,&S2);
 537           for(m=0; (m<NKKKCHI); m++) {
 538             Jc[i][NKKKPHI+NKKKPSI+m]    = kkkchi1[m].Jc;
 539             Jcsig[i][NKKKPHI+NKKKPSI+m] = kkkchi1[m].Jcsig;
 540           }
 541           break;
 542         default: /* covers edOmega and higher Chis than Chi1 */
 543           calc_distribution_props(nbin,histmp,-M_PI,0,NULL,&S2);
 544           break;
 545         }
 546         dlist[i].S2[Dih]        = S2;
 547
 548         /* Sum distribution per amino acid type as well */
 549         for(k=0; (k<nbin); k++) {
 550           his_aa[Dih][dlist[i].index][k] += histmp[k];
 551           histmp[k] = 0;
 552         }
 553         j++;
 554       } else { /* dihed not defined */
 555         dlist[i].S2[Dih] = 0.0 ;
 556       }
 557     }
 558   }
 559   sfree(histmp);
 560
 561   /* Print out Jcouplings */
 562   fprintf(log,"\n *** J-Couplings from simulation (plus std. dev.) ***\n\n");
 563   fprintf(log,"Residue   ");
 564   for(i=0; (i<NKKKPHI); i++)
 565     fprintf(log,"%7s   SD",kkkphi[i].name);
 566   for(i=0; (i<NKKKPSI); i++)
 567     fprintf(log,"%7s   SD",kkkpsi[i].name);
 568   for(i=0; (i<NKKKCHI); i++)
 569     fprintf(log,"%7s   SD",kkkchi1[i].name);
 570   fprintf(log,"\n");
 571   for(i=0; (i<NJC+1); i++)
 572     fprintf(log,"------------");
 573   fprintf(log,"\n");
 574   for(i=0; (i<nlist); i++) {
 575     fprintf(log,"%-10s",dlist[i].name);
 576     for(j=0; (j<NJC); j++)
 577       fprintf(log,"  %5.2f %4.2f",Jc[i][j],Jcsig[i][j]);
 578     fprintf(log,"\n");
 579   }
 580   fprintf(log,"\n");
 581
 582   /* and to -jc file... */
 583   if (bDo_jc) {
 584     fp=xvgropen(fn,"\\S3\\NJ-Couplings from Karplus Equation","Residue",
 585                 "Coupling",oenv);
 586     snew(leg,NJC);
 587     for(i=0; (i<NKKKPHI); i++){
 588                 leg[i] = strdup(kkkphi[i].name);
 589     }
 590     for(i=0; (i<NKKKPSI); i++){
 591                 leg[i+NKKKPHI]=strdup(kkkpsi[i].name);
 592     }
 593     for(i=0; (i<NKKKCHI); i++){
 594       leg[i+NKKKPHI+NKKKPSI]=strdup(kkkchi1[i].name);
 595     }
 596     xvgr_legend(fp,NJC,(const char**)leg,oenv);
 597     fprintf(fp,"%5s ","#Res.");
 598     for(i=0; (i<NJC); i++)
 599       fprintf(fp,"%10s ",leg[i]);
 600     fprintf(fp,"\n");
 601     for(i=0; (i<nlist); i++) {
 602       fprintf(fp,"%5d ",dlist[i].resnr);
 603       for(j=0; (j<NJC); j++)
 604         fprintf(fp,"  %8.3f",Jc[i][j]);
 605       fprintf(fp,"\n");
 606     }
 607     ffclose(fp);
 608     for(i=0; (i<NJC); i++)
 609       sfree(leg[i]);
 610   }
 611   /* finished -jc stuff */
 612
 613   snew(normhisto,nbin);
 614   for(i=0; (i<naa); i++) {
 615     for(Dih=0; (Dih<edMax); Dih++){
 616       /* First check whether something is in there */
 617       for(j=0; (j<nbin); j++)
 618         if (his_aa[Dih][i][j] != 0)
 619           break;
 620       if ((j < nbin) &&
 621           ((bPhi && (Dih==edPhi)) ||
 622            (bPsi && (Dih==edPsi)) ||
 623            (bOmega &&(Dih==edOmega)) ||
 624            (bChi && (Dih>=edChi1)))) {
 625         if (bNormalize)
 626           normalize_histo(nbin,his_aa[Dih][i],(360.0/nbin),normhisto);
 627
 628         switch (Dih) {
 629         case edPhi:
 630           sprintf(hisfile,"histo-phi%s",aa[i]);
 631           sprintf(title,"\\xf\\f{} Distribution for %s",aa[i]);
 632           break;
 633         case edPsi:
 634           sprintf(hisfile,"histo-psi%s",aa[i]);
 635           sprintf(title,"\\xy\\f{} Distribution for %s",aa[i]);
 636           break;
 637         case edOmega:
 638           sprintf(hisfile,"histo-omega%s",aa[i]);
 639           sprintf(title,"\\xw\\f{} Distribution for %s",aa[i]);
 640           break;
 641         default:
 642           sprintf(hisfile,"histo-chi%d%s",Dih-NONCHI+1,aa[i]);
 643           sprintf(title,"\\xc\\f{}\\s%d\\N Distribution for %s",
 644                   Dih-NONCHI+1,aa[i]);
 645         }
 646         strcpy(hhisfile,hisfile);
 647         strcat(hhisfile,".xvg");
 648         fp=xvgropen(hhisfile,title,"Degrees","",oenv);
 649         fprintf(fp,"@ with g0\n");
 650         xvgr_world(fp,-180,0,180,0.1,oenv);
 651         fprintf(fp,"# this effort to set graph size fails unless you run with -autoscale none or -autoscale y flags\n");
 652         fprintf(fp,"@ xaxis tick on\n");
 653         fprintf(fp,"@ xaxis tick major 90\n");
 654         fprintf(fp,"@ xaxis tick minor 30\n");
 655         fprintf(fp,"@ xaxis ticklabel prec 0\n");
 656         fprintf(fp,"@ yaxis tick off\n");
 657         fprintf(fp,"@ yaxis ticklabel off\n");
 658         fprintf(fp,"@ type xy\n");
 659         if (bSSHisto) {
 660           for(k=0; (k<3); k++) {
 661             sprintf(sshisfile,"%s-%s.xvg",hisfile,sss[k]);
 662             ssfp[k] = ffopen(sshisfile,"w");
 663           }
 664         }
 665         for(j=0; (j<nbin); j++) {
 666           angle = -180 + (360/nbin)*j ;
 667           if (bNormalize)
 668             fprintf(fp,"%5d  %10g\n",angle,normhisto[j]);
 669           else
 670             fprintf(fp,"%5d  %10d\n",angle,his_aa[Dih][i][j]);
 671           if (bSSHisto)
 672             for(k=0; (k<3); k++)
 673               fprintf(ssfp[k],"%5d  %10d\n",angle,
 674                       his_aa_ss[k][i][Dih][j]);
 675         }
 676         fprintf(fp,"&\n");
 677         ffclose(fp);
 678         if (bSSHisto) {
 679           for(k=0; (k<3); k++) {
 680             fprintf(ssfp[k],"&\n");
 681             ffclose(ssfp[k]);
 682           }
 683         }
 684       }
 685     }
 686   }
 687   sfree(normhisto);
 688
 689   if (bSSHisto) {
 690     /* Four dimensional array... Very cool */
 691     for(i=0; (i<3); i++) {
 692       for(j=0; (j<=naa); j++) {
 693         for(Dih=0; (Dih<edMax); Dih++)
 694           sfree(his_aa_ss[i][j][Dih]);
 695         sfree(his_aa_ss[i][j]);
 696       }
 697       sfree(his_aa_ss[i]);
 698     }
 699     sfree(his_aa_ss);
 700     sfree(ss_str);
 701   }
 702 }
 703
 704 static FILE *rama_file(const char *fn,const char *title,const char *xaxis,
 705                        const char *yaxis,const output_env_t oenv)
 706 {
 707   FILE *fp;
 708
 709   fp = xvgropen(fn,title,xaxis,yaxis,oenv);
 710   fprintf(fp,"@ with g0\n");
 711   xvgr_world(fp,-180,-180,180,180,oenv);
 712   fprintf(fp,"@ xaxis tick on\n");
 713   fprintf(fp,"@ xaxis tick major 90\n");
 714   fprintf(fp,"@ xaxis tick minor 30\n");
 715   fprintf(fp,"@ xaxis ticklabel prec 0\n");
 716   fprintf(fp,"@ yaxis tick on\n");
 717   fprintf(fp,"@ yaxis tick major 90\n");
 718   fprintf(fp,"@ yaxis tick minor 30\n");
 719   fprintf(fp,"@ yaxis ticklabel prec 0\n");
 720   fprintf(fp,"@    s0 type xy\n");
 721   fprintf(fp,"@    s0 symbol 2\n");
 722   fprintf(fp,"@    s0 symbol size 0.410000\n");
 723   fprintf(fp,"@    s0 symbol fill 1\n");
 724   fprintf(fp,"@    s0 symbol color 1\n");
 725   fprintf(fp,"@    s0 symbol linewidth 1\n");
 726   fprintf(fp,"@    s0 symbol linestyle 1\n");
 727   fprintf(fp,"@    s0 symbol center false\n");
 728   fprintf(fp,"@    s0 symbol char 0\n");
 729   fprintf(fp,"@    s0 skip 0\n");
 730   fprintf(fp,"@    s0 linestyle 0\n");
 731   fprintf(fp,"@    s0 linewidth 1\n");
 732   fprintf(fp,"@ type xy\n");
 733
 734   return fp;
 735 }
 736
 737 static void do_rama(int nf,int nlist,t_dlist dlist[],real **dih,
 738                     gmx_bool bViol,gmx_bool bRamOmega,const output_env_t oenv)
 739 {
 740   FILE *fp,*gp=NULL;
 741   gmx_bool bOm;
 742   char fn[256];
 743   int  i,j,k,Xi1,Xi2,Phi,Psi,Om=0,nlevels;
 744 #define NMAT 120
 745   real **mat=NULL,phi,psi,omega,axis[NMAT],lo,hi;
 746   t_rgb rlo = { 1.0, 0.0, 0.0 };
 747   t_rgb rmid= { 1.0, 1.0, 1.0 };
 748   t_rgb rhi = { 0.0, 0.0, 1.0 };
 749
 750   for(i=0; (i<nlist); i++) {
 751     if ((has_dihedral(edPhi,&(dlist[i]))) &&
 752         (has_dihedral(edPsi,&(dlist[i])))) {
 753       sprintf(fn,"ramaPhiPsi%s.xvg",dlist[i].name);
 754       fp = rama_file(fn,"Ramachandran Plot",
 755                      "\\8f\\4 (deg)","\\8y\\4 (deg)",oenv);
 756       bOm = bRamOmega && has_dihedral(edOmega,&(dlist[i]));
 757       if (bOm) {
 758         Om = dlist[i].j0[edOmega];
 759         snew(mat,NMAT);
 760         for(j=0; (j<NMAT); j++) {
 761           snew(mat[j],NMAT);
 762           axis[j] = -180+(360*j)/NMAT;
 763         }
 764       }
 765       if (bViol) {
 766         sprintf(fn,"violPhiPsi%s.xvg",dlist[i].name);
 767         gp = ffopen(fn,"w");
 768       }
 769       Phi = dlist[i].j0[edPhi];
 770       Psi = dlist[i].j0[edPsi];
 771       for(j=0; (j<nf); j++) {
 772         phi = RAD2DEG*dih[Phi][j];
 773         psi = RAD2DEG*dih[Psi][j];
 774         fprintf(fp,"%10g  %10g\n",phi,psi);
 775         if (bViol)
 776           fprintf(gp,"%d\n",!bAllowed(dih[Phi][j],RAD2DEG*dih[Psi][j]));
 777         if (bOm) {
 778           omega = RAD2DEG*dih[Om][j];
 779           mat[(int)((phi*NMAT)/360)+NMAT/2][(int)((psi*NMAT)/360)+NMAT/2]
 780             += omega;
 781         }
 782       }
 783       if (bViol)
 784         ffclose(gp);
 785       ffclose(fp);
 786       if (bOm) {
 787         sprintf(fn,"ramomega%s.xpm",dlist[i].name);
 788         fp = ffopen(fn,"w");
 789         lo = hi = 0;
 790         for(j=0; (j<NMAT); j++)
 791           for(k=0; (k<NMAT); k++) {
 792             mat[j][k] /= nf;
 793             lo=min(mat[j][k],lo);
 794             hi=max(mat[j][k],hi);
 795           }
 796         /* Symmetrise */
 797         if (fabs(lo) > fabs(hi))
 798           hi = -lo;
 799         else
 800           lo = -hi;
 801         /* Add 180 */
 802         for(j=0; (j<NMAT); j++)
 803           for(k=0; (k<NMAT); k++)
 804             mat[j][k] += 180;
 805         lo += 180;
 806         hi += 180;
 807         nlevels = 20;
 808         write_xpm3(fp,0,"Omega/Ramachandran Plot","Deg","Phi","Psi",
 809                    NMAT,NMAT,axis,axis,mat,lo,180.0,hi,rlo,rmid,rhi,&nlevels);
 810         ffclose(fp);
 811         for(j=0; (j<NMAT); j++)
 812           sfree(mat[j]);
 813         sfree(mat);
 814       }
 815     }
 816     if ((has_dihedral(edChi1,&(dlist[i]))) &&
 817         (has_dihedral(edChi2,&(dlist[i])))) {
 818       sprintf(fn,"ramaX1X2%s.xvg",dlist[i].name);
 819       fp = rama_file(fn,"\\8c\\4\\s1\\N-\\8c\\4\\s2\\N Ramachandran Plot",
 820                      "\\8c\\4\\s1\\N (deg)","\\8c\\4\\s2\\N (deg)",oenv);
 821       Xi1 = dlist[i].j0[edChi1];
 822       Xi2 = dlist[i].j0[edChi2];
 823       for(j=0; (j<nf); j++)
 824         fprintf(fp,"%10g  %10g\n",RAD2DEG*dih[Xi1][j],RAD2DEG*dih[Xi2][j]);
 825       ffclose(fp);
 826     }
 827     else
 828       fprintf(stderr,"No chi1 & chi2 angle for %s\n",dlist[i].name);
 829   }
 830 }
 831
 832
 833 static void print_transitions(const char *fn,int maxchi,int nlist,
 834                               t_dlist dlist[], t_atoms *atoms,rvec x[],
 835                               matrix box, gmx_bool bPhi,gmx_bool bPsi,gmx_bool bChi,real dt,
 836                               const output_env_t oenv)
 837 {
 838   /* based on order_params below */
 839   FILE *fp;
 840   int  nh[edMax];
 841   int  i,Dih,Xi;
 842
 843   /*  must correspond with enum in pp2shift.h:38 */
 844   char *leg[edMax];
 845 #define NLEG asize(leg)
 846
 847   leg[0] = strdup("Phi");
 848   leg[1] = strdup("Psi");
 849   leg[2] = strdup("Omega");
 850   leg[3] = strdup("Chi1");
 851   leg[4] = strdup("Chi2");
 852   leg[5] = strdup("Chi3");
 853   leg[6] = strdup("Chi4");
 854   leg[7] = strdup("Chi5");
 855   leg[8] = strdup("Chi6");
 856
 857   /* Print order parameters */
 858   fp=xvgropen(fn,"Dihedral Rotamer Transitions","Residue","Transitions/ns",
 859                 oenv);
 860   xvgr_legend(fp,NONCHI+maxchi,(const char**)leg,oenv);
 861
 862   for (Dih=0; (Dih<edMax); Dih++)
 863     nh[Dih]=0;
 864
 865   fprintf(fp,"%5s ","#Res.");
 866   fprintf(fp,"%10s %10s %10s ",leg[edPhi],leg[edPsi],leg[edOmega]);
 867   for (Xi=0; Xi<maxchi; Xi++)
 868     fprintf(fp,"%10s ",leg[NONCHI+Xi]);
 869   fprintf(fp,"\n");
 870
 871   for(i=0; (i<nlist); i++) {
 872     fprintf(fp,"%5d ",dlist[i].resnr);
 873     for (Dih=0; (Dih<NONCHI+maxchi); Dih++)
 874       fprintf(fp,"%10.3f ",dlist[i].ntr[Dih]/dt);
 875     /* fprintf(fp,"%12s\n",dlist[i].name);  this confuses xmgrace */
 876     fprintf(fp,"\n");
 877   }
 878   ffclose(fp);
 879 }
 880
 881 static void order_params(FILE *log,
 882                          const char *fn,int maxchi,int nlist,t_dlist dlist[],
 883                          const char *pdbfn,real bfac_init,
 884                          t_atoms *atoms,rvec x[],int ePBC,matrix box,
 885                          gmx_bool bPhi,gmx_bool bPsi,gmx_bool bChi,const output_env_t oenv)
 886 {
 887   FILE *fp;
 888   int  nh[edMax];
 889   char buf[STRLEN];
 890   int  i,Dih,Xi;
 891   real S2Max, S2Min;
 892
 893   /* except for S2Min/Max, must correspond with enum in pp2shift.h:38 */
 894   const char *const_leg[2+edMax]= { "S2Min","S2Max","Phi","Psi","Omega",
 895                                     "Chi1", "Chi2", "Chi3", "Chi4", "Chi5",
 896                                     "Chi6" };
 897 #define NLEG asize(leg)
 898
 899   char *leg[2+edMax];
 900
 901   for(i=0;i<NLEG;i++)
 902     leg[i]=strdup(const_leg[i]);
 903
 904   /* Print order parameters */
 905   fp=xvgropen(fn,"Dihedral Order Parameters","Residue","S2",oenv);
 906   xvgr_legend(fp,2+NONCHI+maxchi,const_leg,oenv);
 907
 908   for (Dih=0; (Dih<edMax); Dih++)
 909     nh[Dih]=0;
 910
 911   fprintf(fp,"%5s ","#Res.");
 912   fprintf(fp,"%10s %10s ",leg[0],leg[1]);
 913   fprintf(fp,"%10s %10s %10s ",leg[2+edPhi],leg[2+edPsi],leg[2+edOmega]);
 914   for (Xi=0; Xi<maxchi; Xi++)
 915     fprintf(fp,"%10s ",leg[2+NONCHI+Xi]);
 916   fprintf(fp,"\n");
 917
 918   for(i=0; (i<nlist); i++) {
 919     S2Max=-10;
 920     S2Min=10;
 921     for (Dih=0; (Dih<NONCHI+maxchi); Dih++) {
 922       if (dlist[i].S2[Dih]!=0) {
 923         if (dlist[i].S2[Dih] > S2Max)
 924           S2Max=dlist[i].S2[Dih];
 925         if (dlist[i].S2[Dih] < S2Min)
 926           S2Min=dlist[i].S2[Dih];
 927       }
 928       if (dlist[i].S2[Dih] > 0.8)
 929         nh[Dih]++;
 930     }
 931     fprintf(fp,"%5d ",dlist[i].resnr);
 932     fprintf(fp,"%10.3f %10.3f ",S2Min,S2Max);
 933     for (Dih=0; (Dih<NONCHI+maxchi); Dih++)
 934       fprintf(fp,"%10.3f ",dlist[i].S2[Dih]);
 935     fprintf(fp,"\n");
 936     /* fprintf(fp,"%12s\n",dlist[i].name);  this confuses xmgrace */
 937   }
 938   ffclose(fp);
 939
 940   if (NULL != pdbfn) {
 941     real x0,y0,z0;
 942
 943     if (NULL == atoms->pdbinfo)
 944       snew(atoms->pdbinfo,atoms->nr);
 945     for(i=0; (i<atoms->nr); i++)
 946       atoms->pdbinfo[i].bfac=bfac_init;
 947
 948     for(i=0; (i<nlist); i++) {
 949       atoms->pdbinfo[dlist[i].atm.N].bfac=-dlist[i].S2[0];/* Phi */
 950       atoms->pdbinfo[dlist[i].atm.H].bfac=-dlist[i].S2[0];/* Phi */
 951       atoms->pdbinfo[dlist[i].atm.C].bfac=-dlist[i].S2[1];/* Psi */
 952       atoms->pdbinfo[dlist[i].atm.O].bfac=-dlist[i].S2[1];/* Psi */
 953       for (Xi=0; (Xi<maxchi); Xi++) {           /* Chi's */
 954         if (dlist[i].atm.Cn[Xi+3]!=-1) {
 955           atoms->pdbinfo[dlist[i].atm.Cn[Xi+1]].bfac=-dlist[i].S2[NONCHI+Xi];
 956         }
 957       }
 958     }
 959
 960     fp=ffopen(pdbfn,"w");
 961     fprintf(fp,"REMARK generated by g_chi\n");
 962     fprintf(fp,"REMARK "
 963             "B-factor field contains negative of dihedral order parameters\n");
 964     write_pdbfile(fp,NULL,atoms,x,ePBC,box,' ',0,NULL,TRUE);
 965     x0=y0=z0=1000.0;
 966     for (i=0; (i<atoms->nr); i++) {
 967       x0 = min(x0, x[i][XX]);
 968       y0 = min(y0, x[i][YY]);
 969       z0 = min(z0, x[i][ZZ]);
 970     }
 971     x0*=10.0;/* nm -> angstrom */
 972     y0*=10.0;/* nm -> angstrom */
 973     z0*=10.0;/* nm -> angstrom */
 974     sprintf(buf,"%s%%6.f%%6.2f\n",pdbformat);
 975     for (i=0; (i<10); i++) {
 976       fprintf(fp,buf,"ATOM  ", atoms->nr+1+i, "CA", "LEG",' ',
 977               atoms->nres+1, ' ',x0, y0, z0+(1.2*i), 0.0, -0.1*i);
 978     }
 979     ffclose(fp);
 980   }
 981
 982   fprintf(log,"Dihedrals with S2 > 0.8\n");
 983   fprintf(log,"Dihedral: ");
 984   if (bPhi) fprintf(log," Phi  ");
 985   if (bPsi) fprintf(log," Psi ");
 986   if (bChi)
 987     for(Xi=0; (Xi<maxchi); Xi++)
 988       fprintf(log," %s ", leg[2+NONCHI+Xi]);
 989   fprintf(log,"\nNumber:   ");
 990   if (bPhi) fprintf(log,"%4d  ",nh[0]);
 991   if (bPsi) fprintf(log,"%4d  ",nh[1]);
 992   if (bChi)
 993     for(Xi=0; (Xi<maxchi); Xi++)
 994       fprintf(log,"%4d  ",nh[NONCHI+Xi]);
 995   fprintf(log,"\n");
 996
 997   for(i=0;i<NLEG;i++)
 998     sfree(leg[i]);
 999
1000 }
1001
1002 int gmx_chi(int argc,char *argv[])
1003 {
1004   const char *desc[] = {
1005     "g_chi computes phi, psi, omega and chi dihedrals for all your ",
1006     "amino acid backbone and sidechains.",
1007     "It can compute dihedral angle as a function of time, and as",
1008     "histogram distributions.",
1009     "The distributions (histo-(dihedral)(RESIDUE).xvg) are cumulative over all residues of each type.[PAR]",
1010     "If option [TT]-corr[tt] is given, the program will",
1011     "calculate dihedral autocorrelation functions. The function used",
1012     "is C(t) = < cos(chi(tau)) cos(chi(tau+t)) >. The use of cosines",
1013     "rather than angles themselves, resolves the problem of periodicity.",
1014     "(Van der Spoel & Berendsen (1997), [BB]Biophys. J. 72[bb], 2032-2041).",
1015     "Separate files for each dihedral of each residue",
1016     "(corr(dihedral)(RESIDUE)(nresnr).xvg) are output, as well as a",
1017     "file containing the information for all residues (argument of [TT]-corr[tt]).[PAR]",
1018     "With option [TT]-all[tt], the angles themselves as a function of time for",
1019     "each residue are printed to separate files (dihedral)(RESIDUE)(nresnr).xvg.",
1020     "These can be in radians or degrees.[PAR]",
1021     "A log file (argument [TT]-g[tt]) is also written. This contains [BR]",
1022     "(a) information about the number of residues of each type.[BR]",
1023     "(b) The NMR 3J coupling constants from the Karplus equation.[BR]",
1024     "(c) a table for each residue of the number of transitions between ",
1025     "rotamers per nanosecond,  and the order parameter S2 of each dihedral.[BR]",
1026     "(d) a table for each residue of the rotamer occupancy.[BR]",
1027     "All rotamers are taken as 3-fold, except for omegas and chi-dihedrals",
1028     "to planar groups (i.e. chi2 of aromatics asp and asn, chi3 of glu",
1029     "and gln, and chi4 of arg), which are 2-fold. \"rotamer 0\" means ",
1030     "that the dihedral was not in the core region of each rotamer. ",
1031     "The width of the core region can be set with [TT]-core_rotamer[tt][PAR]",
1032
1033     "The S2 order parameters are also output to an xvg file",
1034     "(argument [TT]-o[tt] ) and optionally as a pdb file with",
1035     "the S2 values as B-factor (argument [TT]-p[tt]). ",
1036     "The total number of rotamer transitions per timestep",
1037     "(argument [TT]-ot[tt]), the number of transitions per rotamer",
1038     "(argument [TT]-rt[tt]), and the 3J couplings (argument [TT]-jc[tt]), ",
1039     "can also be written to .xvg files.[PAR]",
1040
1041     "If [TT]-chi_prod[tt] is set (and maxchi > 0), cumulative rotamers, e.g.",
1042     "1+9(chi1-1)+3(chi2-1)+(chi3-1) (if the residue has three 3-fold ",
1043     "dihedrals and maxchi >= 3)",
1044     "are calculated. As before, if any dihedral is not in the core region,",
1045     "the rotamer is taken to be 0. The occupancies of these cumulative ",
1046     "rotamers (starting with rotamer 0) are written to the file",
1047     "that is the argument of [TT]-cp[tt], and if the [TT]-all[tt] flag",
1048     "is given, the rotamers as functions of time",
1049     "are written to chiproduct(RESIDUE)(nresnr).xvg ",
1050     "and their occupancies to histo-chiproduct(RESIDUE)(nresnr).xvg.[PAR]",
1051
1052     "The option [TT]-r[tt] generates a contour plot of the average omega angle",
1053     "as a function of the phi and psi angles, that is, in a Ramachandran plot",
1054     "the average omega angle is plotted using color coding.",
1055
1056   };
1057
1058   const char *bugs[] = {
1059     "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.",
1060     "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.",
1061     "-r0 option does not work properly",
1062     "Rotamers with multiplicity 2 are printed in chi.log as if they had multiplicity 3, with the 3rd (g(+)) always having probability 0"
1063   };
1064
1065   /* defaults */
1066   static int  r0=1,ndeg=1,maxchi=2;
1067   static gmx_bool bAll=FALSE;
1068   static gmx_bool bPhi=FALSE,bPsi=FALSE,bOmega=FALSE;
1069   static real bfac_init=-1.0,bfac_max=0;
1070   static const char *maxchistr[] = { NULL, "0", "1", "2", "3",  "4", "5", "6", NULL };
1071   static gmx_bool bRama=FALSE,bShift=FALSE,bViol=FALSE,bRamOmega=FALSE;
1072   static gmx_bool bNormHisto=TRUE,bChiProduct=FALSE,bHChi=FALSE,bRAD=FALSE,bPBC=TRUE;
1073   static real core_frac=0.5 ;
1074   t_pargs pa[] = {
1075     { "-r0",  FALSE, etINT, {&r0},
1076       "starting residue" },
1077     { "-phi",  FALSE, etBOOL, {&bPhi},
1078       "Output for Phi dihedral angles" },
1079     { "-psi",  FALSE, etBOOL, {&bPsi},
1080       "Output for Psi dihedral angles" },
1081     { "-omega",FALSE, etBOOL, {&bOmega},
1082       "Output for Omega dihedrals (peptide bonds)" },
1083     { "-rama", FALSE, etBOOL, {&bRama},
1084       "Generate Phi/Psi and Chi1/Chi2 ramachandran plots" },
1085     { "-viol", FALSE, etBOOL, {&bViol},
1086       "Write a file that gives 0 or 1 for violated Ramachandran angles" },
1087     { "-periodic", FALSE, etBOOL, {&bPBC},
1088       "Print dihedral angles modulo 360 degrees" },
1089     { "-all",  FALSE, etBOOL, {&bAll},
1090       "Output separate files for every dihedral." },
1091     { "-rad",  FALSE, etBOOL, {&bRAD},
1092       "in angle vs time files, use radians rather than degrees."},
1093     { "-shift", FALSE, etBOOL, {&bShift},
1094         "Compute chemical shifts from Phi/Psi angles" },
1095     { "-binwidth", FALSE, etINT, {&ndeg},
1096       "bin width for histograms (degrees)" },
1097     { "-core_rotamer", FALSE, etREAL, {&core_frac},
1098       "only the central -core_rotamer*(360/multiplicity) belongs to each rotamer (the rest is assigned to rotamer 0)" },
1099     { "-maxchi", FALSE, etENUM, {maxchistr},
1100       "calculate first ndih Chi dihedrals" },
1101     { "-normhisto", FALSE, etBOOL, {&bNormHisto},
1102       "Normalize histograms" },
1103     { "-ramomega",FALSE,etBOOL, {&bRamOmega},
1104       "compute average omega as a function of phi/psi and plot it in an xpm plot" },
1105     { "-bfact", FALSE, etREAL, {&bfac_init},
1106       "B-factor value for pdb file for atoms with no calculated dihedral order parameter"},
1107     { "-chi_prod",FALSE,etBOOL, {&bChiProduct},
1108       "compute a single cumulative rotamer for each residue"},
1109     { "-HChi",FALSE,etBOOL, {&bHChi},
1110       "Include dihedrals to sidechain hydrogens"},
1111     { "-bmax",  FALSE, etREAL, {&bfac_max},
1112       "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." }
1113   };
1114
1115   FILE       *log;
1116   int        natoms,nlist,naa,idum,nbin;
1117   t_atoms    atoms;
1118   rvec       *x;
1119   int        ePBC;
1120   matrix     box;
1121   char       title[256],grpname[256];
1122   t_dlist    *dlist;
1123   char       **aa;
1124   gmx_bool       bChi,bCorr,bSSHisto;
1125   gmx_bool       bDo_rt, bDo_oh, bDo_ot, bDo_jc ;
1126   real       dt=0, traj_t_ns;
1127   output_env_t oenv;
1128
1129   atom_id    isize,*index;
1130   int        ndih,nactdih,nf;
1131   real       **dih,*trans_frac,*aver_angle,*time;
1132   int        i,j,**chi_lookup,*xity;
1133
1134   t_filenm  fnm[] = {
1135     { efSTX, "-s",  NULL,     ffREAD  },
1136     { efTRX, "-f",  NULL,     ffREAD  },
1137     { efXVG, "-o",  "order",  ffWRITE },
1138     { efPDB, "-p",  "order",  ffOPTWR },
1139     { efDAT, "-ss", "ssdump", ffOPTRD },
1140     { efXVG, "-jc", "Jcoupling", ffWRITE },
1141     { efXVG, "-corr",  "dihcorr",ffOPTWR },
1142     { efLOG, "-g",  "chi",    ffWRITE },
1143     /* add two more arguments copying from g_angle */
1144     { efXVG, "-ot", "dihtrans", ffOPTWR },
1145     { efXVG, "-oh", "trhisto",  ffOPTWR },
1146     { efXVG, "-rt", "restrans",  ffOPTWR },
1147     { efXVG, "-cp", "chiprodhisto",  ffOPTWR }
1148   };
1149 #define NFILE asize(fnm)
1150   int     npargs;
1151   t_pargs *ppa;
1152
1153   CopyRight(stderr,argv[0]);
1154   npargs = asize(pa);
1155   ppa    = add_acf_pargs(&npargs,pa);
1156   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
1157                     NFILE,fnm,npargs,ppa,asize(desc),desc,asize(bugs),bugs,
1158                     &oenv);
1159
1160   /* Handle result from enumerated type */
1161   sscanf(maxchistr[0],"%d",&maxchi);
1162   bChi = (maxchi > 0);
1163
1164   log=ffopen(ftp2fn(efLOG,NFILE,fnm),"w");
1165
1166   if (bRamOmega) {
1167     bOmega = TRUE;
1168     bPhi   = TRUE;
1169     bPsi   = TRUE;
1170   }
1171
1172   /* set some options */
1173   bDo_rt=(opt2bSet("-rt",NFILE,fnm));
1174   bDo_oh=(opt2bSet("-oh",NFILE,fnm));
1175   bDo_ot=(opt2bSet("-ot",NFILE,fnm));
1176   bDo_jc=(opt2bSet("-jc",NFILE,fnm));
1177   bCorr=(opt2bSet("-corr",NFILE,fnm));
1178   if (bCorr)
1179     fprintf(stderr,"Will calculate autocorrelation\n");
1180
1181   if (core_frac > 1.0 ) {
1182     fprintf(stderr, "core_rotamer fraction > 1.0 ; will use 1.0\n");
1183     core_frac=1.0 ;
1184   }
1185   if (core_frac < 0.0 ) {
1186     fprintf(stderr, "core_rotamer fraction < 0.0 ; will use 0.0\n");
1187     core_frac=0.0 ;
1188   }
1189
1190   if (maxchi > MAXCHI) {
1191     fprintf(stderr,
1192             "Will only calculate first %d Chi dihedrals in stead of %d.\n",
1193             MAXCHI, maxchi);
1194     maxchi=MAXCHI;
1195   }
1196   bSSHisto = ftp2bSet(efDAT,NFILE,fnm);
1197   nbin = 360/ndeg ;
1198
1199   /* Find the chi angles using atoms struct and a list of amino acids */
1200   get_stx_coordnum(ftp2fn(efSTX,NFILE,fnm),&natoms);
1201   init_t_atoms(&atoms,natoms,TRUE);
1202   snew(x,natoms);
1203   read_stx_conf(ftp2fn(efSTX,NFILE,fnm),title,&atoms,x,NULL,&ePBC,box);
1204   fprintf(log,"Title: %s\n",title);
1205
1206   naa=get_strings("aminoacids.dat",&aa);
1207   dlist=mk_dlist(log,&atoms,&nlist,bPhi,bPsi,bChi,bHChi,maxchi,r0,naa,aa);
1208   fprintf(stderr,"%d residues with dihedrals found\n", nlist);
1209
1210   if (nlist == 0)
1211     gmx_fatal(FARGS,"No dihedrals in your structure!\n");
1212
1213   /* Make a linear index for reading all. */
1214   index=make_chi_ind(nlist,dlist,&ndih);
1215   isize=4*ndih;
1216   fprintf(stderr,"%d dihedrals found\n", ndih);
1217
1218   snew(dih,ndih);
1219
1220   /* COMPUTE ALL DIHEDRALS! */
1221   read_ang_dih(ftp2fn(efTRX,NFILE,fnm),FALSE,TRUE,FALSE,bPBC,1,&idum,
1222                &nf,&time,isize,index,&trans_frac,&aver_angle,dih,oenv);
1223
1224   dt=(time[nf-1]-time[0])/(nf-1); /* might want this for corr or n. transit*/
1225   if (bCorr)
1226   {
1227           if (nf < 2)
1228           {
1229                   gmx_fatal(FARGS,"Need at least 2 frames for correlation");
1230           }
1231   }
1232
1233   /* put angles in -M_PI to M_PI ! and correct phase factor for phi and psi
1234   * pass nactdih instead of ndih to low_ana_dih_trans and get_chi_product_traj
1235   * to prevent accessing off end of arrays when maxchi < 5 or 6. */
1236   nactdih = reset_em_all(nlist,dlist,nf,dih,maxchi);
1237
1238   if (bAll)
1239     dump_em_all(nlist,dlist,nf,time,dih,maxchi,bPhi,bPsi,bChi,bOmega,bRAD,oenv);
1240
1241   /* Histogramming & J coupling constants & calc of S2 order params */
1242   histogramming(log,nbin,naa,aa,nf,maxchi,dih,nlist,dlist,index,
1243                 bPhi,bPsi,bOmega,bChi,
1244                 bNormHisto,bSSHisto,ftp2fn(efDAT,NFILE,fnm),bfac_max,&atoms,
1245                 bDo_jc,opt2fn("-jc",NFILE,fnm),oenv);
1246
1247   /* transitions
1248    *
1249    * added multiplicity */
1250
1251   snew(xity,ndih) ;
1252   mk_multiplicity_lookup(xity, maxchi, dih, nlist, dlist,ndih);
1253
1254   strcpy(grpname, "All residues, ");
1255   if(bPhi)
1256     strcat(grpname, "Phi ");
1257   if(bPsi)
1258     strcat(grpname, "Psi ");
1259   if(bOmega)
1260     strcat(grpname, "Omega ");
1261   if(bChi){
1262     strcat(grpname, "Chi 1-") ;
1263     sprintf(grpname + strlen(grpname), "%i", maxchi);
1264   }
1265
1266
1267   low_ana_dih_trans(bDo_ot, opt2fn("-ot",NFILE,fnm),
1268                     bDo_oh, opt2fn("-oh",NFILE,fnm),maxchi,
1269                     dih, nlist, dlist, nf, nactdih, grpname, xity,
1270                     *time,  dt, FALSE, core_frac,oenv) ;
1271
1272   /* Order parameters */
1273   order_params(log,opt2fn("-o",NFILE,fnm),maxchi,nlist,dlist,
1274                ftp2fn_null(efPDB,NFILE,fnm),bfac_init,
1275                &atoms,x,ePBC,box,bPhi,bPsi,bChi,oenv);
1276
1277   /* Print ramachandran maps! */
1278   if (bRama)
1279     do_rama(nf,nlist,dlist,dih,bViol,bRamOmega,oenv);
1280
1281   if (bShift)
1282     do_pp2shifts(log,nf,nlist,dlist,dih);
1283
1284   /* rprint S^2, transitions, and rotamer occupancies to log */
1285   traj_t_ns = 0.001 * (time[nf-1]-time[0]) ;
1286   pr_dlist(log,nlist,dlist,traj_t_ns,edPrintST,bPhi,bPsi,bChi,bOmega,maxchi);
1287   pr_dlist(log,nlist,dlist,traj_t_ns,edPrintRO,bPhi,bPsi,bChi,bOmega,maxchi);
1288   ffclose(log);
1289   /* transitions to xvg */
1290   if (bDo_rt)
1291     print_transitions(opt2fn("-rt",NFILE,fnm),maxchi,nlist,dlist,
1292                       &atoms,x,box,bPhi,bPsi,bChi,traj_t_ns,oenv);
1293
1294   /* chi_product trajectories (ie one "rotamer number" for each residue) */
1295   if (bChiProduct && bChi ) {
1296     snew(chi_lookup,nlist) ;
1297     for (i=0;i<nlist;i++)
1298       snew(chi_lookup[i], maxchi) ;
1299     mk_chi_lookup(chi_lookup, maxchi, dih, nlist, dlist);
1300
1301     get_chi_product_traj(dih,nf,nactdih,nlist,
1302                          maxchi,dlist,time,chi_lookup,xity,
1303                          FALSE,bNormHisto, core_frac,bAll,
1304                          opt2fn("-cp",NFILE,fnm),oenv);
1305
1306     for (i=0;i<nlist;i++)
1307       sfree(chi_lookup[i]);
1308   }
1309
1310   /* Correlation comes last because it fucks up the angles */
1311   if (bCorr)
1312     do_dihcorr(opt2fn("-corr",NFILE,fnm),nf,ndih,dih,dt,nlist,dlist,time,
1313                maxchi,bPhi,bPsi,bChi,bOmega,oenv);
1314
1315
1316   do_view(oenv,opt2fn("-o",NFILE,fnm),"-nxy");
1317   do_view(oenv,opt2fn("-jc",NFILE,fnm),"-nxy");
1318   if (bCorr)
1319     do_view(oenv,opt2fn("-corr",NFILE,fnm),"-nxy");
1320
1321   thanx(stderr);
1322
1323   return 0;
1324 }