src/tools/dlist.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_dlist_c = "$Id$";
  30
  31 #include <stdlib.h>
  32 #include "string2.h"
  33 #include "pp2shift.h"
  34 #include "smalloc.h"
  35
  36 t_dlist *mk_dlist(FILE *log,
  37                   t_atoms *atoms, int *nlist,
  38                   bool bPhi, bool bPsi, bool bChi, int maxchi,
  39                   int r0,int naa,char **aa)
  40 {
  41   int     ires,i,j,k;
  42   t_dihatms atm,prev;
  43   int     nl=0,nc[edMax],ndih;
  44   bool    bDih;
  45   char    *thisres;
  46   t_dlist *dl;
  47
  48   snew(dl,atoms->nres+1);
  49   prev.C = -1;
  50   for(i=0; (i<edMax); i++)
  51     nc[i]=0;
  52   ires = -1;
  53   i    =  0;
  54   while (i<atoms->nr) {
  55     ires=atoms->atom[i].resnr;
  56
  57     /* Initiate all atom numbers to -1 */
  58     atm.minC=atm.H=atm.N=atm.C=atm.O=-1;
  59     for(j=0; (j<MAXCHI+3); j++)
  60       atm.Cn[j]=-1;
  61
  62     /* Look for atoms in this residue */
  63     while ((i<atoms->nr) && (atoms->atom[i].resnr == ires)) {
  64       if ((strcmp(*(atoms->atomname[i]),"H") == 0) ||
  65           (strcmp(*(atoms->atomname[i]),"H1") == 0) )
  66         atm.H=i;
  67       else if (strcmp(*(atoms->atomname[i]),"N") == 0)
  68         atm.N=i;
  69       else if (strcmp(*(atoms->atomname[i]),"C") == 0)
  70         atm.C=i;
  71       else if ((strcmp(*(atoms->atomname[i]),"O") == 0) ||
  72                (strcmp(*(atoms->atomname[i]),"O1") == 0))
  73         atm.O=i;
  74       else if (strcmp(*(atoms->atomname[i]),"CA") == 0)
  75         atm.Cn[1]=i;
  76       else if (strcmp(*(atoms->atomname[i]),"CB") == 0)
  77         atm.Cn[2]=i;
  78       else if ((strcmp(*(atoms->atomname[i]),"CG") == 0)  ||
  79                (strcmp(*(atoms->atomname[i]),"CG1") == 0) ||
  80                (strcmp(*(atoms->atomname[i]),"OG") == 0)  ||
  81                (strcmp(*(atoms->atomname[i]),"OG1") == 0) ||
  82                (strcmp(*(atoms->atomname[i]),"SG") == 0))
  83         atm.Cn[3]=i;
  84       else if ((strcmp(*(atoms->atomname[i]),"CD") == 0)  ||
  85                (strcmp(*(atoms->atomname[i]),"CD1") == 0) ||
  86                (strcmp(*(atoms->atomname[i]),"SD") == 0)  ||
  87                (strcmp(*(atoms->atomname[i]),"OD1") == 0) ||
  88                (strcmp(*(atoms->atomname[i]),"ND1") == 0) ||
  89                (strcmp(*(atoms->atomname[i]),"HG")  == 0) ||
  90                (strcmp(*(atoms->atomname[i]),"HG1")  == 0))
  91         atm.Cn[4]=i;
  92       else if ((strcmp(*(atoms->atomname[i]),"CE") == 0) ||
  93                (strcmp(*(atoms->atomname[i]),"CE1") == 0) ||
  94                (strcmp(*(atoms->atomname[i]),"OE1") == 0) ||
  95                (strcmp(*(atoms->atomname[i]),"NE") == 0))
  96         atm.Cn[5]=i;
  97       else if ((strcmp(*(atoms->atomname[i]),"CZ") == 0) ||
  98                (strcmp(*(atoms->atomname[i]),"NZ") == 0))
  99         atm.Cn[6]=i;
 100       else if (strcmp(*(atoms->atomname[i]),"NH1") == 0)
 101         atm.Cn[7]=i;
 102       i++;
 103     }
 104
 105     /* Special case for Pro, has no H */
 106     if (strcmp(*(atoms->resname[ires]),"PRO") == 0)
 107       atm.H=atm.Cn[4];
 108     /* Carbon from previous residue */
 109     if (prev.C != -1)
 110       atm.minC = prev.C;
 111     if (prev.O != -1)
 112       atm.minO = prev.O;
 113     prev = atm;
 114
 115     thisres=*(atoms->resname[ires]);
 116
 117     /* Check how many dihedrals we have */
 118     if ((atm.N != -1) && (atm.Cn[1] != -1) && (atm.C != -1) &&
 119       (atm.O != -1) && ((atm.H != -1) || (atm.minC != -1))) {
 120       dl[nl].resnr     = ires+1;
 121       dl[nl].atm       = atm;
 122       dl[nl].atm.Cn[0] = atm.N;
 123       if ((atm.Cn[3] != -1) && (atm.Cn[2] != -1) && (atm.Cn[1] != -1)) {
 124         nc[0]++;
 125         if (atm.Cn[4] != -1) {
 126           nc[1]++;
 127           if (atm.Cn[5] != -1) {
 128             nc[2]++;
 129             if (atm.Cn[6] != -1) {
 130               nc[3]++;
 131               if (atm.Cn[7] != -1) {
 132                 nc[4]++;
 133                 if (atm.Cn[8] != -1) {
 134                   nc[5]++;
 135                 }
 136               }
 137             }
 138           }
 139         }
 140       }
 141       if ((atm.minC != -1) && (atm.minO != -1))
 142         nc[6]++;
 143       for(k=0; (k<naa); k++) {
 144         if (strcasecmp(aa[k],thisres) == 0)
 145           break;
 146       }
 147       dl[nl].index=k;
 148
 149       sprintf(dl[nl].name,"%s%d",thisres,ires+r0);
 150       nl++;
 151     }
 152     else if (debug)
 153       fprintf(debug,"Could not find N atom but could find other atoms"
 154               " in residue %s%d\n",thisres,ires+r0);
 155   }
 156   fprintf(stderr,"\n");
 157   fprintf(log,"\n");
 158   fprintf(log,"There are %d residues with dihedrals\n",nl);
 159   j=0;
 160   if (bPhi) j+=nl;
 161   if (bPsi) j+=nl;
 162   if (bChi)
 163     for(i=0; (i<maxchi); i++)
 164       j+=nc[i];
 165   fprintf(log,"There are %d dihedrals\n",j);
 166   fprintf(log,"Dihedral: ");
 167   if (bPhi)
 168     fprintf(log," Phi  ");
 169   if (bPsi)
 170     fprintf(log," Psi  ");
 171   if (bChi)
 172     for(i=0; (i<maxchi); i++)
 173       fprintf(log,"Chi%d  ",i+1);
 174   fprintf(log,"\nNumber:   ");
 175   if (bPhi)
 176     fprintf(log,"%4d  ",nl);
 177   if (bPsi)
 178     fprintf(log,"%4d  ",nl);
 179   if (bChi)
 180     for(i=0; (i<maxchi); i++)
 181       fprintf(log,"%4d  ",nc[i]);
 182   fprintf(log,"\n");
 183
 184   *nlist=nl;
 185
 186   return dl;
 187 }
 188
 189 bool has_dihedral(int Dih,t_dlist *dl)
 190 {
 191   bool b = FALSE;
 192   int  ddd;
 193
 194 #define BBB(x) (dl->atm.##x != -1)
 195   switch (Dih) {
 196   case edPhi:
 197     b = (BBB(H) && BBB(N) && BBB(Cn[1]) && BBB(C));
 198     break;
 199   case edPsi:
 200     b = (BBB(N) && BBB(Cn[1]) && BBB(C) && BBB(O));
 201     break;
 202   case edOmega:
 203     b = (BBB(minO) && BBB(minC) && BBB(N) && BBB(Cn[1]));
 204     break;
 205   case edChi1:
 206   case edChi2:
 207   case edChi3:
 208   case edChi4:
 209   case edChi5:
 210   case edChi6:
 211     ddd = Dih - edChi1;
 212     b   = (BBB(Cn[ddd]) && BBB(Cn[ddd+1]) && BBB(Cn[ddd+2]) && BBB(Cn[ddd+3]));
 213     break;
 214   default:
 215     pr_dlist(stdout,1,dl,1);
 216     fatal_error(0,"Non existant dihedral %d in file %s, line %d",
 217                 Dih,__FILE__,__LINE__);
 218   }
 219   return b;
 220 }
 221
 222 static void pr_props(FILE *fp,t_dlist *dl,int nDih,real dt)
 223 {
 224   fprintf(fp,"  %6.2f  %6.2f\n",dl->ntr[nDih]/dt,dl->S2[nDih]);
 225 }
 226
 227 void pr_dlist(FILE *fp,int nl,t_dlist dl[],real dt)
 228 {
 229   int i, Xi;
 230
 231   for(i=0; (i<nl); i++) {
 232     fprintf(fp,"Residue %s\n",dl[i].name);
 233     fprintf(fp," Angle [  AI,  AJ,  AK,  AL]  #tr/ns  S^2D  \n"
 234                "--------------------------------------------\n");
 235     fprintf(fp,"   Phi [%4d,%4d,%4d,%4d]",
 236             (dl[i].atm.H == -1) ? dl[i].atm.minC : dl[i].atm.H,
 237             dl[i].atm.N,dl[i].atm.Cn[1],dl[i].atm.C);
 238     pr_props(fp,&dl[i],edPhi,dt);
 239     fprintf(fp,"   Psi [%4d,%4d,%4d,%4d]",dl[i].atm.N,dl[i].atm.Cn[1],
 240             dl[i].atm.C,dl[i].atm.O);
 241     pr_props(fp,&dl[i],edPsi,dt);
 242     fprintf(fp," Omega [%4d,%4d,%4d,%4d]",dl[i].atm.minO,dl[i].atm.minC,
 243             dl[i].atm.N,dl[i].atm.Cn[1]);
 244     pr_props(fp,&dl[i],edOmega,dt);
 245     for(Xi=0; Xi<MAXCHI; Xi++)
 246       if (dl[i].atm.Cn[Xi+3] != -1) {
 247         fprintf(fp,"   Chi%d[%4d,%4d,%4d,%4d]",Xi+1,dl[i].atm.Cn[Xi],
 248                 dl[i].atm.Cn[Xi+1],dl[i].atm.Cn[Xi+2],
 249                 dl[i].atm.Cn[Xi+3]);
 250         pr_props(fp,&dl[i],Xi+2,dt);
 251       }
 252     fprintf(fp,"\n");
 253   }
 254 }
 255
 256 int pr_trans(FILE *fp,int nl,t_dlist dl[],real dt,int Xi)
 257 {
 258   int  i,nn,nz;
 259
 260   nz=0;
 261   fprintf(fp,"\\begin{table}[h]\n");
 262   fprintf(fp,"\\caption{Number of dihedral transitions per nanosecond}\n");
 263   fprintf(fp,"\\begin{tabular}{|l|l|}\n");
 264   fprintf(fp,"\\hline\n");
 265   fprintf(fp,"Residue\t&$\\chi_%d$\t\\\\\n",Xi+1);
 266   for(i=0; (i<nl); i++) {
 267     nn=dl[i].ntr[Xi]/dt;
 268
 269     if (nn == 0) {
 270       fprintf(fp,"%s\t&\\HL{%d}\t\\\\\n",dl[i].name,nn);
 271       nz++;
 272     }
 273     else if (nn > 0)
 274       fprintf(fp,"%s\t&\\%d\t\\\\\n",dl[i].name,nn);
 275   }
 276   fprintf(fp,"\\hline\n");
 277   fprintf(fp,"\\end{tabular}\n");
 278   fprintf(fp,"\\end{table}\n\n");
 279
 280   return nz;
 281 }
 282