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     "[TT]g_chi[tt] computes [GRK]phi[grk], [GRK]psi[grk], [GRK]omega[grk], and [GRK]chi[grk] 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 [TT](histo-(dihedral)(RESIDUE).xvg[tt]) 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) = [CHEVRON][COS][GRK]chi[grk]([GRK]tau[grk])[cos] [COS][GRK]chi[grk]([GRK]tau[grk]+t)[cos][chevron]. 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     "[TT](corr(dihedral)(RESIDUE)(nresnr).xvg[tt]) 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 [TT](dihedral)(RESIDUE)(nresnr).xvg[tt].",
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 S^2 of each dihedral.[BR]",
1030     "(d) a table for each residue of the rotamer occupancy.[PAR]",
1031     "All rotamers are taken as 3-fold, except for [GRK]omega[grk] and [GRK]chi[grk] dihedrals",
1032     "to planar groups (i.e. [GRK]chi[grk][SUB]2[sub] of aromatics, Asp and Asn; [GRK]chi[grk][SUB]3[sub] of Glu",
1033     "and Gln; and [GRK]chi[grk][SUB]4[sub] 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 S^2 order parameters are also output to an [TT].xvg[tt] file",
1038     "(argument [TT]-o[tt] ) and optionally as a [TT].pdb[tt] file with",
1039     "the S^2 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 [TT].xvg[tt] files. Note that the analysis",
1044     "of rotamer transitions assumes that the supplied trajectory frames",
1045     "are equally spaced in time.[PAR]",
1046
1047     "If [TT]-chi_prod[tt] is set (and [TT]-maxchi[tt] > 0), cumulative rotamers, e.g.",
1048     "1+9([GRK]chi[grk][SUB]1[sub]-1)+3([GRK]chi[grk][SUB]2[sub]-1)+([GRK]chi[grk][SUB]3[sub]-1) (if the residue has three 3-fold ",
1049     "dihedrals and [TT]-maxchi[tt] >= 3)",
1050     "are calculated. As before, if any dihedral is not in the core region,",
1051     "the rotamer is taken to be 0. The occupancies of these cumulative ",
1052     "rotamers (starting with rotamer 0) are written to the file",
1053     "that is the argument of [TT]-cp[tt], and if the [TT]-all[tt] flag",
1054     "is given, the rotamers as functions of time",
1055     "are written to [TT]chiproduct(RESIDUE)(nresnr).xvg[tt] ",
1056     "and their occupancies to [TT]histo-chiproduct(RESIDUE)(nresnr).xvg[tt].[PAR]",
1057
1058     "The option [TT]-r[tt] generates a contour plot of the average [GRK]omega[grk] angle",
1059     "as a function of the [GRK]phi[grk] and [GRK]psi[grk] angles, that is, in a Ramachandran plot",
1060     "the average [GRK]omega[grk] angle is plotted using color coding.",
1061
1062   };
1063
1064   const char *bugs[] = {
1065     "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.",
1066     "[GRK]phi[grk] and [GRK]psi[grk] dihedrals are calculated in a non-standard way, using H-N-CA-C for [GRK]phi[grk] instead of C(-)-N-CA-C, and N-CA-C-O for [GRK]psi[grk] instead of N-CA-C-N(+). This causes (usually small) discrepancies with the output of other tools like [TT]g_rama[tt].",
1067     "[TT]-r0[tt] option does not work properly",
1068     "Rotamers with multiplicity 2 are printed in [TT]chi.log[tt] as if they had multiplicity 3, with the 3rd (g(+)) always having probability 0"
1069   };
1070
1071   /* defaults */
1072   static int  r0=1,ndeg=1,maxchi=2;
1073   static gmx_bool bAll=FALSE;
1074   static gmx_bool bPhi=FALSE,bPsi=FALSE,bOmega=FALSE;
1075   static real bfac_init=-1.0,bfac_max=0;
1076   static const char *maxchistr[] = { NULL, "0", "1", "2", "3",  "4", "5", "6", NULL };
1077   static gmx_bool bRama=FALSE,bShift=FALSE,bViol=FALSE,bRamOmega=FALSE;
1078   static gmx_bool bNormHisto=TRUE,bChiProduct=FALSE,bHChi=FALSE,bRAD=FALSE,bPBC=TRUE;
1079   static real core_frac=0.5 ;
1080   t_pargs pa[] = {
1081     { "-r0",  FALSE, etINT, {&r0},
1082       "starting residue" },
1083     { "-phi",  FALSE, etBOOL, {&bPhi},
1084       "Output for [GRK]phi[grk] dihedral angles" },
1085     { "-psi",  FALSE, etBOOL, {&bPsi},
1086       "Output for [GRK]psi[grk] dihedral angles" },
1087     { "-omega",FALSE, etBOOL, {&bOmega},
1088       "Output for [GRK]omega[grk] dihedrals (peptide bonds)" },
1089     { "-rama", FALSE, etBOOL, {&bRama},
1090       "Generate [GRK]phi[grk]/[GRK]psi[grk] and [GRK]chi[grk][SUB]1[sub]/[GRK]chi[grk][SUB]2[sub] Ramachandran plots" },
1091     { "-viol", FALSE, etBOOL, {&bViol},
1092       "Write a file that gives 0 or 1 for violated Ramachandran angles" },
1093     { "-periodic", FALSE, etBOOL, {&bPBC},
1094       "Print dihedral angles modulo 360 degrees" },
1095     { "-all",  FALSE, etBOOL, {&bAll},
1096       "Output separate files for every dihedral." },
1097     { "-rad",  FALSE, etBOOL, {&bRAD},
1098       "in angle vs time files, use radians rather than degrees."},
1099     { "-shift", FALSE, etBOOL, {&bShift},
1100         "Compute chemical shifts from [GRK]phi[grk]/[GRK]psi[grk] angles" },
1101     { "-binwidth", FALSE, etINT, {&ndeg},
1102       "bin width for histograms (degrees)" },
1103     { "-core_rotamer", FALSE, etREAL, {&core_frac},
1104       "only the central [TT]-core_rotamer[tt]*(360/multiplicity) belongs to each rotamer (the rest is assigned to rotamer 0)" },
1105     { "-maxchi", FALSE, etENUM, {maxchistr},
1106       "calculate first ndih [GRK]chi[grk] dihedrals" },
1107     { "-normhisto", FALSE, etBOOL, {&bNormHisto},
1108       "Normalize histograms" },
1109     { "-ramomega",FALSE,etBOOL, {&bRamOmega},
1110       "compute average omega as a function of [GRK]phi[grk]/[GRK]psi[grk] and plot it in an [TT].xpm[tt] plot" },
1111     { "-bfact", FALSE, etREAL, {&bfac_init},
1112       "B-factor value for [TT].pdb[tt] file for atoms with no calculated dihedral order parameter"},
1113     { "-chi_prod",FALSE,etBOOL, {&bChiProduct},
1114       "compute a single cumulative rotamer for each residue"},
1115     { "-HChi",FALSE,etBOOL, {&bHChi},
1116       "Include dihedrals to sidechain hydrogens"},
1117     { "-bmax",  FALSE, etREAL, {&bfac_max},
1118       "Maximum B-factor on any of the atoms that make up a dihedral, for the dihedral angle to be considere in the statistics. Applies to database work where a number of X-Ray structures is analyzed. [TT]-bmax[tt] <= 0 means no limit." }
1119   };
1120
1121   FILE       *log;
1122   int        natoms,nlist,idum,nbin;
1123   t_atoms    atoms;
1124   rvec       *x;
1125   int        ePBC;
1126   matrix     box;
1127   char       title[256],grpname[256];
1128   t_dlist    *dlist;
1129   gmx_bool       bChi,bCorr,bSSHisto;
1130   gmx_bool       bDo_rt, bDo_oh, bDo_ot, bDo_jc ;
1131   real       dt=0, traj_t_ns;
1132   output_env_t oenv;
1133   gmx_residuetype_t rt;
1134
1135   atom_id    isize,*index;
1136   int        ndih,nactdih,nf;
1137   real       **dih,*trans_frac,*aver_angle,*time;
1138   int        i,j,**chi_lookup,*multiplicity;
1139
1140   t_filenm  fnm[] = {
1141     { efSTX, "-s",  NULL,     ffREAD  },
1142     { efTRX, "-f",  NULL,     ffREAD  },
1143     { efXVG, "-o",  "order",  ffWRITE },
1144     { efPDB, "-p",  "order",  ffOPTWR },
1145     { efDAT, "-ss", "ssdump", ffOPTRD },
1146     { efXVG, "-jc", "Jcoupling", ffWRITE },
1147     { efXVG, "-corr",  "dihcorr",ffOPTWR },
1148     { efLOG, "-g",  "chi",    ffWRITE },
1149     /* add two more arguments copying from g_angle */
1150     { efXVG, "-ot", "dihtrans", ffOPTWR },
1151     { efXVG, "-oh", "trhisto",  ffOPTWR },
1152     { efXVG, "-rt", "restrans",  ffOPTWR },
1153     { efXVG, "-cp", "chiprodhisto",  ffOPTWR }
1154   };
1155 #define NFILE asize(fnm)
1156   int     npargs;
1157   t_pargs *ppa;
1158
1159   CopyRight(stderr,argv[0]);
1160   npargs = asize(pa);
1161   ppa    = add_acf_pargs(&npargs,pa);
1162   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
1163                     NFILE,fnm,npargs,ppa,asize(desc),desc,asize(bugs),bugs,
1164                     &oenv);
1165
1166   /* Handle result from enumerated type */
1167   sscanf(maxchistr[0],"%d",&maxchi);
1168   bChi = (maxchi > 0);
1169
1170   log=ffopen(ftp2fn(efLOG,NFILE,fnm),"w");
1171
1172   if (bRamOmega) {
1173     bOmega = TRUE;
1174     bPhi   = TRUE;
1175     bPsi   = TRUE;
1176   }
1177
1178   /* set some options */
1179   bDo_rt=(opt2bSet("-rt",NFILE,fnm));
1180   bDo_oh=(opt2bSet("-oh",NFILE,fnm));
1181   bDo_ot=(opt2bSet("-ot",NFILE,fnm));
1182   bDo_jc=(opt2bSet("-jc",NFILE,fnm));
1183   bCorr=(opt2bSet("-corr",NFILE,fnm));
1184   if (bCorr)
1185     fprintf(stderr,"Will calculate autocorrelation\n");
1186
1187   if (core_frac > 1.0 ) {
1188     fprintf(stderr, "core_rotamer fraction > 1.0 ; will use 1.0\n");
1189     core_frac=1.0 ;
1190   }
1191   if (core_frac < 0.0 ) {
1192     fprintf(stderr, "core_rotamer fraction < 0.0 ; will use 0.0\n");
1193     core_frac=0.0 ;
1194   }
1195
1196   if (maxchi > MAXCHI) {
1197     fprintf(stderr,
1198             "Will only calculate first %d Chi dihedrals in stead of %d.\n",
1199             MAXCHI, maxchi);
1200     maxchi=MAXCHI;
1201   }
1202   bSSHisto = ftp2bSet(efDAT,NFILE,fnm);
1203   nbin = 360/ndeg ;
1204
1205   /* Find the chi angles using atoms struct and a list of amino acids */
1206   get_stx_coordnum(ftp2fn(efSTX,NFILE,fnm),&natoms);
1207   init_t_atoms(&atoms,natoms,TRUE);
1208   snew(x,natoms);
1209   read_stx_conf(ftp2fn(efSTX,NFILE,fnm),title,&atoms,x,NULL,&ePBC,box);
1210   fprintf(log,"Title: %s\n",title);
1211
1212   gmx_residuetype_init(&rt);
1213   dlist=mk_dlist(log,&atoms,&nlist,bPhi,bPsi,bChi,bHChi,maxchi,r0,rt);
1214   fprintf(stderr,"%d residues with dihedrals found\n", nlist);
1215
1216   if (nlist == 0)
1217     gmx_fatal(FARGS,"No dihedrals in your structure!\n");
1218
1219   /* Make a linear index for reading all. */
1220   index=make_chi_ind(nlist,dlist,&ndih);
1221   isize=4*ndih;
1222   fprintf(stderr,"%d dihedrals found\n", ndih);
1223
1224   snew(dih,ndih);
1225
1226   /* COMPUTE ALL DIHEDRALS! */
1227   read_ang_dih(ftp2fn(efTRX,NFILE,fnm),FALSE,TRUE,FALSE,bPBC,1,&idum,
1228                &nf,&time,isize,index,&trans_frac,&aver_angle,dih,oenv);
1229
1230   dt=(time[nf-1]-time[0])/(nf-1); /* might want this for corr or n. transit*/
1231   if (bCorr)
1232   {
1233           if (nf < 2)
1234           {
1235                   gmx_fatal(FARGS,"Need at least 2 frames for correlation");
1236           }
1237   }
1238
1239   /* put angles in -M_PI to M_PI ! and correct phase factor for phi and psi
1240   * pass nactdih instead of ndih to low_ana_dih_trans and get_chi_product_traj
1241   * to prevent accessing off end of arrays when maxchi < 5 or 6. */
1242   nactdih = reset_em_all(nlist,dlist,nf,dih,maxchi);
1243
1244   if (bAll)
1245     dump_em_all(nlist,dlist,nf,time,dih,maxchi,bPhi,bPsi,bChi,bOmega,bRAD,oenv);
1246
1247   /* Histogramming & J coupling constants & calc of S2 order params */
1248   histogramming(log,nbin,rt,nf,maxchi,dih,nlist,dlist,index,
1249                 bPhi,bPsi,bOmega,bChi,
1250                 bNormHisto,bSSHisto,ftp2fn(efDAT,NFILE,fnm),bfac_max,&atoms,
1251                 bDo_jc,opt2fn("-jc",NFILE,fnm),oenv);
1252
1253   /* transitions
1254    *
1255    * added multiplicity */
1256
1257   snew(multiplicity,ndih) ;
1258   mk_multiplicity_lookup(multiplicity, maxchi, dih, nlist, dlist,ndih);
1259
1260   strcpy(grpname, "All residues, ");
1261   if(bPhi)
1262     strcat(grpname, "Phi ");
1263   if(bPsi)
1264     strcat(grpname, "Psi ");
1265   if(bOmega)
1266     strcat(grpname, "Omega ");
1267   if(bChi){
1268     strcat(grpname, "Chi 1-") ;
1269     sprintf(grpname + strlen(grpname), "%i", maxchi);
1270   }
1271
1272
1273   low_ana_dih_trans(bDo_ot, opt2fn("-ot",NFILE,fnm),
1274                     bDo_oh, opt2fn("-oh",NFILE,fnm),maxchi,
1275                     dih, nlist, dlist, nf, nactdih, grpname, multiplicity,
1276                     time, FALSE, core_frac,oenv);
1277
1278   /* Order parameters */
1279   order_params(log,opt2fn("-o",NFILE,fnm),maxchi,nlist,dlist,
1280                ftp2fn_null(efPDB,NFILE,fnm),bfac_init,
1281                &atoms,x,ePBC,box,bPhi,bPsi,bChi,oenv);
1282
1283   /* Print ramachandran maps! */
1284   if (bRama)
1285     do_rama(nf,nlist,dlist,dih,bViol,bRamOmega,oenv);
1286
1287   if (bShift)
1288     do_pp2shifts(log,nf,nlist,dlist,dih);
1289
1290   /* rprint S^2, transitions, and rotamer occupancies to log */
1291   traj_t_ns = 0.001 * (time[nf-1]-time[0]) ;
1292   pr_dlist(log,nlist,dlist,traj_t_ns,edPrintST,bPhi,bPsi,bChi,bOmega,maxchi);
1293   pr_dlist(log,nlist,dlist,traj_t_ns,edPrintRO,bPhi,bPsi,bChi,bOmega,maxchi);
1294   ffclose(log);
1295   /* transitions to xvg */
1296   if (bDo_rt)
1297     print_transitions(opt2fn("-rt",NFILE,fnm),maxchi,nlist,dlist,
1298                       &atoms,x,box,bPhi,bPsi,bChi,traj_t_ns,oenv);
1299
1300   /* chi_product trajectories (ie one "rotamer number" for each residue) */
1301   if (bChiProduct && bChi ) {
1302     snew(chi_lookup,nlist) ;
1303     for (i=0;i<nlist;i++)
1304       snew(chi_lookup[i], maxchi) ;
1305     mk_chi_lookup(chi_lookup, maxchi, dih, nlist, dlist);
1306
1307     get_chi_product_traj(dih,nf,nactdih,nlist,
1308                          maxchi,dlist,time,chi_lookup,multiplicity,
1309                          FALSE,bNormHisto, core_frac,bAll,
1310                          opt2fn("-cp",NFILE,fnm),oenv);
1311
1312     for (i=0;i<nlist;i++)
1313       sfree(chi_lookup[i]);
1314   }
1315
1316   /* Correlation comes last because it fucks up the angles */
1317   if (bCorr)
1318     do_dihcorr(opt2fn("-corr",NFILE,fnm),nf,ndih,dih,dt,nlist,dlist,time,
1319                maxchi,bPhi,bPsi,bChi,bOmega,oenv);
1320
1321
1322   do_view(oenv,opt2fn("-o",NFILE,fnm),"-nxy");
1323   do_view(oenv,opt2fn("-jc",NFILE,fnm),"-nxy");
1324   if (bCorr)
1325     do_view(oenv,opt2fn("-corr",NFILE,fnm),"-nxy");
1326
1327   gmx_residuetype_destroy(rt);
1328
1329   thanx(stderr);
1330
1331   return 0;
1332 }