src/tools/gmx_chi.c

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