src/tools/dlist.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
  39 #include <stdlib.h>
  40
  41 #include "string2.h"
  42 #include "smalloc.h"
  43 #include "gstat.h"
  44 #include "gmx_fatal.h"
  45
  46 t_dlist *mk_dlist(FILE *log,
  47                   t_atoms *atoms, int *nlist,
  48                   bool bPhi, bool bPsi, bool bChi, bool bHChi,
  49                   int maxchi,int r0,int naa,char **aa)
  50 {
  51   int     ires,i,j,k,ii;
  52   t_dihatms atm,prev;
  53   int     nl=0,nc[edMax];
  54   char    *thisres;
  55   t_dlist *dl;
  56
  57   snew(dl,atoms->nres+1);
  58   prev.C = prev.O = -1;
  59   for(i=0; (i<edMax); i++)
  60     nc[i]=0;
  61   ires = -1;
  62   i    =  0;
  63   while (i<atoms->nr) {
  64     ires = atoms->atom[i].resind;
  65
  66     /* Initiate all atom numbers to -1 */
  67     atm.minC=atm.H=atm.N=atm.C=atm.O=atm.minO=-1;
  68     for(j=0; (j<MAXCHI+3); j++)
  69       atm.Cn[j]=-1;
  70
  71     /* Look for atoms in this residue */
  72     /* maybe should allow for chis to hydrogens? */
  73     while ((i<atoms->nr) && (atoms->atom[i].resind == ires)) {
  74       if ((strcmp(*(atoms->atomname[i]),"H") == 0) ||
  75           (strcmp(*(atoms->atomname[i]),"H1") == 0) )
  76         atm.H=i;
  77       else if (strcmp(*(atoms->atomname[i]),"N") == 0)
  78         atm.N=i;
  79       else if (strcmp(*(atoms->atomname[i]),"C") == 0)
  80         atm.C=i;
  81       else if ((strcmp(*(atoms->atomname[i]),"O") == 0) ||
  82                (strcmp(*(atoms->atomname[i]),"O1") == 0))
  83         atm.O=i;
  84       else if (strcmp(*(atoms->atomname[i]),"CA") == 0)
  85         atm.Cn[1]=i;
  86       else if (strcmp(*(atoms->atomname[i]),"CB") == 0)
  87         atm.Cn[2]=i;
  88       else if ((strcmp(*(atoms->atomname[i]),"CG") == 0)  ||
  89                (strcmp(*(atoms->atomname[i]),"CG1") == 0) ||
  90                (strcmp(*(atoms->atomname[i]),"OG") == 0)  ||
  91                (strcmp(*(atoms->atomname[i]),"OG1") == 0) ||
  92                (strcmp(*(atoms->atomname[i]),"SG") == 0))
  93         atm.Cn[3]=i;
  94       else if ((strcmp(*(atoms->atomname[i]),"CD") == 0)  ||
  95                (strcmp(*(atoms->atomname[i]),"CD1") == 0) ||
  96                (strcmp(*(atoms->atomname[i]),"SD") == 0)  ||
  97                (strcmp(*(atoms->atomname[i]),"OD1") == 0) ||
  98                (strcmp(*(atoms->atomname[i]),"ND1") == 0))
  99         atm.Cn[4]=i;
 100       /* by grs - split the Cn[4] into 2 bits to check allowing dih to H */
 101       else if (bHChi && ((strcmp(*(atoms->atomname[i]),"HG")  == 0) ||
 102                (strcmp(*(atoms->atomname[i]),"HG1")  == 0)) )
 103         atm.Cn[4]=i;
 104       else if ((strcmp(*(atoms->atomname[i]),"CE") == 0) ||
 105                (strcmp(*(atoms->atomname[i]),"CE1") == 0) ||
 106                (strcmp(*(atoms->atomname[i]),"OE1") == 0) ||
 107                (strcmp(*(atoms->atomname[i]),"NE") == 0))
 108         atm.Cn[5]=i;
 109       else if ((strcmp(*(atoms->atomname[i]),"CZ") == 0) ||
 110                (strcmp(*(atoms->atomname[i]),"NZ") == 0))
 111         atm.Cn[6]=i;
 112       /* HChi flag here too */
 113       else if (bHChi && (strcmp(*(atoms->atomname[i]),"NH1") == 0))
 114         atm.Cn[7]=i;
 115       i++;
 116     }
 117
 118     thisres = *(atoms->resinfo[ires].name);
 119
 120     /* added by grs - special case for aromatics, whose chis above 2 are
 121        not real and produce rubbish output - so set back to -1 */
 122     if (strcmp(thisres,"PHE") == 0 ||
 123         strcmp(thisres,"TYR") == 0 ||
 124         strcmp(thisres,"PTR") == 0 ||
 125         strcmp(thisres,"TRP") == 0 ||
 126         strcmp(thisres,"HIS") == 0 ||
 127         strcmp(thisres,"HISA") == 0 ||
 128         strcmp(thisres,"HISB") == 0 )  {
 129       for (ii=5 ; ii<=7 ; ii++)
 130         atm.Cn[ii]=-1;
 131     }
 132     /* end fixing aromatics */
 133
 134     /* Special case for Pro, has no H */
 135     if (strcmp(thisres,"PRO") == 0)
 136       atm.H=atm.Cn[4];
 137     /* Carbon from previous residue */
 138     if (prev.C != -1)
 139       atm.minC = prev.C;
 140     if (prev.O != -1)
 141       atm.minO = prev.O;
 142     prev = atm;
 143
 144     /* Check how many dihedrals we have */
 145     if ((atm.N != -1) && (atm.Cn[1] != -1) && (atm.C != -1) &&
 146       (atm.O != -1) && ((atm.H != -1) || (atm.minC != -1))) {
 147       dl[nl].resnr     = ires+1;
 148       dl[nl].atm       = atm;
 149       dl[nl].atm.Cn[0] = atm.N;
 150       if ((atm.Cn[3] != -1) && (atm.Cn[2] != -1) && (atm.Cn[1] != -1)) {
 151         nc[0]++;
 152         if (atm.Cn[4] != -1) {
 153           nc[1]++;
 154           if (atm.Cn[5] != -1) {
 155             nc[2]++;
 156             if (atm.Cn[6] != -1) {
 157               nc[3]++;
 158               if (atm.Cn[7] != -1) {
 159                 nc[4]++;
 160                 if (atm.Cn[8] != -1) {
 161                   nc[5]++;
 162                 }
 163               }
 164             }
 165           }
 166         }
 167       }
 168       if ((atm.minC != -1) && (atm.minO != -1))
 169         nc[6]++;
 170       for(k=0; (k<naa); k++) {
 171         if (gmx_strcasecmp(aa[k],thisres) == 0)
 172           break;
 173       }
 174       dl[nl].index=k;
 175
 176       sprintf(dl[nl].name,"%s%d",thisres,ires+r0);
 177       nl++;
 178     }
 179     else if (debug)
 180       fprintf(debug,"Could not find N atom but could find other atoms"
 181               " in residue %s%d\n",thisres,ires+r0);
 182   }
 183   fprintf(stderr,"\n");
 184   fprintf(log,"\n");
 185   fprintf(log,"There are %d residues with dihedrals\n",nl);
 186   j=0;
 187   if (bPhi) j+=nl;
 188   if (bPsi) j+=nl;
 189   if (bChi)
 190     for(i=0; (i<maxchi); i++)
 191       j+=nc[i];
 192   fprintf(log,"There are %d dihedrals\n",j);
 193   fprintf(log,"Dihedral: ");
 194   if (bPhi)
 195     fprintf(log," Phi  ");
 196   if (bPsi)
 197     fprintf(log," Psi  ");
 198   if (bChi)
 199     for(i=0; (i<maxchi); i++)
 200       fprintf(log,"Chi%d  ",i+1);
 201   fprintf(log,"\nNumber:   ");
 202   if (bPhi)
 203     fprintf(log,"%4d  ",nl);
 204   if (bPsi)
 205     fprintf(log,"%4d  ",nl);
 206   if (bChi)
 207     for(i=0; (i<maxchi); i++)
 208       fprintf(log,"%4d  ",nc[i]);
 209   fprintf(log,"\n");
 210
 211   *nlist=nl;
 212
 213   return dl;
 214 }
 215
 216 bool has_dihedral(int Dih,t_dlist *dl)
 217 {
 218   bool b = FALSE;
 219   int  ddd;
 220
 221   switch (Dih) {
 222   case edPhi:
 223     b = ((dl->atm.H!=-1) && (dl->atm.N!=-1) && (dl->atm.Cn[1]!=-1) && (dl->atm.C!=-1));
 224     break;
 225   case edPsi:
 226     b = ((dl->atm.N!=-1) && (dl->atm.Cn[1]!=-1) && (dl->atm.C!=-1) && (dl->atm.O!=-1));
 227     break;
 228   case edOmega:
 229     b = ((dl->atm.minO!=-1) && (dl->atm.minC!=-1) && (dl->atm.N!=-1) && (dl->atm.Cn[1]!=-1));
 230     break;
 231   case edChi1:
 232   case edChi2:
 233   case edChi3:
 234   case edChi4:
 235   case edChi5:
 236   case edChi6:
 237     ddd = Dih - edChi1;
 238     b   = ((dl->atm.Cn[ddd]!=-1) &&  (dl->atm.Cn[ddd+1]!=-1)&&
 239            (dl->atm.Cn[ddd+2]!=-1) && (dl->atm.Cn[ddd+3]!=-1));
 240     break;
 241   default:
 242     pr_dlist(stdout,1,dl,1,0,TRUE,TRUE,TRUE,TRUE,MAXCHI);
 243     gmx_fatal(FARGS,"Non existant dihedral %d in file %s, line %d",
 244                 Dih,__FILE__,__LINE__);
 245   }
 246   return b;
 247 }
 248
 249 static void pr_one_ro(FILE *fp,t_dlist *dl,int nDih,real dt)
 250 {
 251   int k ;
 252   for(k=0;k<NROT;k++)
 253     fprintf(fp,"  %6.2f",dl->rot_occ[nDih][k]);
 254   fprintf(fp,"\n");
 255 }
 256
 257 static void pr_ntr_s2(FILE *fp,t_dlist *dl,int nDih,real dt)
 258 {
 259   fprintf(fp,"  %6.2f  %6.2f\n",(dt == 0) ? 0 : dl->ntr[nDih]/dt,dl->S2[nDih]);
 260 }
 261
 262 void pr_dlist(FILE *fp,int nl,t_dlist dl[],real dt, int printtype,
 263 bool bPhi, bool bPsi,bool bChi,bool bOmega, int maxchi)
 264 {
 265   int i, Xi;
 266
 267   void  (*pr_props)(FILE *, t_dlist *, int, real);
 268
 269   /* Analysis of dihedral transitions etc */
 270
 271   if (printtype == edPrintST){
 272     pr_props=pr_ntr_s2 ;
 273     fprintf(stderr,"Now printing out transitions and OPs...\n");
 274   }else{
 275     pr_props=pr_one_ro ;
 276     fprintf(stderr,"Now printing out rotamer occupancies...\n");
 277     fprintf(fp,"\nXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX\n\n");
 278   }
 279
 280   /* change atom numbers from 0 based to 1 based */
 281   for(i=0; (i<nl); i++) {
 282     fprintf(fp,"Residue %s\n",dl[i].name);
 283     if (printtype == edPrintST){
 284       fprintf(fp," Angle [   AI,   AJ,   AK,   AL]  #tr/ns  S^2D  \n"
 285                  "--------------------------------------------\n");
 286     } else {
 287       fprintf(fp," Angle [   AI,   AJ,   AK,   AL]  rotamers  0  g(-)  t  g(+)\n"
 288                  "--------------------------------------------\n");
 289     }
 290     if (bPhi) {
 291       fprintf(fp,"   Phi [%5d,%5d,%5d,%5d]",
 292             (dl[i].atm.H == -1) ? 1+dl[i].atm.minC : 1+dl[i].atm.H,
 293             1+dl[i].atm.N, 1+dl[i].atm.Cn[1], 1+dl[i].atm.C);
 294       pr_props(fp,&dl[i],edPhi,dt);
 295     }
 296     if (bPsi) {
 297       fprintf(fp,"   Psi [%5d,%5d,%5d,%5d]",1+dl[i].atm.N, 1+dl[i].atm.Cn[1],
 298             1+dl[i].atm.C, 1+dl[i].atm.O);
 299       pr_props(fp,&dl[i],edPsi,dt);
 300     }
 301     if (bOmega && has_dihedral(edOmega,&(dl[i]))) {
 302       fprintf(fp," Omega [%5d,%5d,%5d,%5d]",1+dl[i].atm.minO, 1+dl[i].atm.minC,
 303               1+dl[i].atm.N, 1+dl[i].atm.Cn[1]);
 304       pr_props(fp,&dl[i],edOmega,dt);
 305     }
 306     for(Xi=0; Xi<MAXCHI; Xi++)
 307       if (bChi && (Xi < maxchi) && (dl[i].atm.Cn[Xi+3] != -1) ) {
 308         fprintf(fp,"   Chi%d[%5d,%5d,%5d,%5d]",Xi+1, 1+dl[i].atm.Cn[Xi],
 309                 1+dl[i].atm.Cn[Xi+1], 1+dl[i].atm.Cn[Xi+2],
 310                 1+dl[i].atm.Cn[Xi+3]);
 311         pr_props(fp,&dl[i],Xi+edChi1,dt); /* Xi+2 was wrong here */
 312       }
 313     fprintf(fp,"\n");
 314   }
 315 }
 316
 317
 318
 319 int pr_trans(FILE *fp,int nl,t_dlist dl[],real dt,int Xi)
 320 {
 321   /* never called at the moment */
 322
 323   int  i,nn,nz;
 324
 325   nz=0;
 326   fprintf(fp,"\\begin{table}[h]\n");
 327   fprintf(fp,"\\caption{Number of dihedral transitions per nanosecond}\n");
 328   fprintf(fp,"\\begin{tabular}{|l|l|}\n");
 329   fprintf(fp,"\\hline\n");
 330   fprintf(fp,"Residue\t&$\\chi_%d$\t\\\\\n",Xi+1);
 331   for(i=0; (i<nl); i++) {
 332     nn=dl[i].ntr[Xi]/dt;
 333
 334     if (nn == 0) {
 335       fprintf(fp,"%s\t&\\HL{%d}\t\\\\\n",dl[i].name,nn);
 336       nz++;
 337     }
 338     else if (nn > 0)
 339       fprintf(fp,"%s\t&\\%d\t\\\\\n",dl[i].name,nn);
 340   }
 341   fprintf(fp,"\\hline\n");
 342   fprintf(fp,"\\end{tabular}\n");
 343   fprintf(fp,"\\end{table}\n\n");
 344
 345   return nz;
 346 }
 347