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