src/kernel/hizzie.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_hizzie_c = "$Id$";
  30
  31 #include <stdio.h>
  32 #include <string.h>
  33 #include "typedefs.h"
  34 #include "pdbio.h"
  35 #include "smalloc.h"
  36 #include "vec.h"
  37 #include "physics.h"
  38 #include "toputil.h"
  39 #include "pdb2top.h"
  40
  41 static int in_strings(char *key,int nstr,char **str)
  42 {
  43   int j;
  44
  45   for(j=0; (j<nstr); j++)
  46     if (strcmp(str[j],key) == 0)
  47       return j;
  48
  49   return -1;
  50 }
  51
  52 static bool hbond(rvec x[],int i,int j,real distance)
  53 {
  54   real tol = distance*distance;
  55   rvec   tmp;
  56
  57   rvec_sub(x[i],x[j],tmp);
  58
  59   return (iprod(tmp,tmp) < tol);
  60 }
  61
  62 static void chk_allhb(t_atoms *pdba,rvec x[],t_block *hb,
  63                       bool donor[],bool accept[],real dist)
  64 {
  65   int i,j,k,ii,natom;
  66
  67   natom=pdba->nr;
  68   snew(hb->index,natom+1);
  69   snew(hb->a,6*natom);
  70   hb->nr  = natom;
  71   hb->nra = 6*natom;
  72
  73   k = ii = 0;
  74   hb->index[ii++] = 0;
  75   for(i=0; (i<natom); i++) {
  76     if (donor[i]) {
  77       for(j=i+1; (j<natom); j++)
  78         if ((accept[j]) && (hbond(x,i,j,dist)))
  79           hb->a[k++] = j;
  80     }
  81     else if (accept[i]) {
  82       for(j=i+1; (j<natom); j++)
  83         if ((donor[j]) && (hbond(x,i,j,dist)))
  84           hb->a[k++] = j;
  85     }
  86     hb->index[ii++] = k;
  87   }
  88   hb->nra = k;
  89 }
  90
  91 static void pr_hbonds(FILE *fp,t_block *hb,t_atoms *pdba)
  92 {
  93   int i,j,k,j0,j1;
  94
  95   fprintf(fp,"Dumping all hydrogen bonds!\n");
  96   for(i=0; (i<hb->nr); i++) {
  97     j0=hb->index[i];
  98     j1=hb->index[i+1];
  99     for(j=j0; (j<j1); j++) {
 100       k=hb->a[j];
 101       fprintf(fp,"%5s%4d%5s - %5s%4d%5s\n",
 102               *pdba->resname[pdba->atom[i].resnr],
 103               pdba->atom[i].resnr+1,*pdba->atomname[i],
 104               *pdba->resname[pdba->atom[k].resnr],
 105               pdba->atom[k].resnr+1,*pdba->atomname[k]);
 106     }
 107   }
 108 }
 109
 110 static bool chk_hbonds(int i,t_atoms *pdba, rvec x[],
 111                        bool ad[],bool hbond[],rvec xh,
 112                        real angle,real dist)
 113 {
 114   bool bHB;
 115   int  j,aj,ri,natom;
 116   real d2,dist2,a;
 117   rvec nh,oh;
 118
 119   natom=pdba->nr;
 120   bHB = FALSE;
 121   ri = pdba->atom[i].resnr;
 122   dist2=sqr(dist);
 123   for(j=0; (j<natom); j++) {
 124     /* Check whether the other atom is a donor/acceptor and not i */
 125     if ((ad[j]) && (j != i)) {
 126       /* Check whether the other atom is on the same ring as well */
 127       if ((pdba->atom[j].resnr != ri) ||
 128           ((strcmp(*pdba->atomname[j],"ND1") != 0) &&
 129            (strcmp(*pdba->atomname[j],"NE2") != 0))) {
 130         aj = j;
 131         d2  = distance2(x[i],x[j]);
 132         rvec_sub(x[i],xh,nh);
 133         rvec_sub(x[aj],xh,oh);
 134         a  = RAD2DEG * acos(cos_angle(nh,oh));
 135         if ((d2 < dist2) && (a > angle)) {
 136           if (debug)
 137             fprintf(debug,
 138                     "HBOND between %s%d-%s and %s%d-%s is %g nm, %g deg\n",
 139                     *pdba->resname[pdba->atom[i].resnr],
 140                     pdba->atom[i].resnr+1, *pdba->atomname[i],
 141                     *pdba->resname[pdba->atom[aj].resnr],
 142                     pdba->atom[aj].resnr+1,*pdba->atomname[aj],sqrt(d2),a);
 143           hbond[i] = TRUE;
 144           bHB      = TRUE;
 145         }
 146       }
 147     }
 148   }
 149   return bHB;
 150 }
 151
 152 static void calc_ringh(rvec xattach,rvec xb,rvec xc,rvec xh)
 153 {
 154   rvec tab,tac;
 155   real n;
 156
 157   /* Add a proton on a ring to atom attach at distance 0.1 nm */
 158   rvec_sub(xattach,xb,tab);
 159   rvec_sub(xattach,xc,tac);
 160   rvec_add(tab,tac,xh);
 161   n=0.1/norm(xh);
 162   svmul(n,xh,xh);
 163   rvec_inc(xh,xattach);
 164 }
 165
 166 void set_histp(t_atoms *pdba,rvec *x,real angle,real dist){
 167   static char *prot_acc[] = {
 168     "O", "OD1", "OD2", "OE1", "OE2", "OG", "OG1", "OH", "OW"
 169   };
 170 #define NPA asize(prot_acc)
 171   static char *prot_don[] = {
 172     "N", "NH1", "NH2", "NE", "ND1", "ND2", "NE2", "NZ", "OG", "OG1", "OH", "NE1", "OW"
 173   };
 174 #define NPD asize(prot_don)
 175
 176   bool *donor,*acceptor;
 177   bool *hbond,bHaveH=FALSE;
 178   bool bHDd,bHEd;
 179   rvec xh1,xh2;
 180   int  natom;
 181   int  i,j,nd,na,aj,hisnr,his0,type=-1;
 182   int  nd1,ne2,cg,cd2,ce1;
 183   t_block *hb;
 184   real d;
 185   char *atomnm;
 186
 187   natom=pdba->nr;
 188   snew(donor,natom);
 189   snew(acceptor,natom);
 190   snew(hbond,natom);
 191   snew(hb,1);
 192
 193   nd=na=0;
 194   for(i=0; (i<natom); i++) {
 195     if (in_strings(*pdba->atomname[i],NPA,prot_acc) != -1) {
 196       acceptor[i] = TRUE;
 197       na++;
 198     }
 199     if (in_strings(*pdba->atomname[i],NPD,prot_don) != -1) {
 200       donor[i] = TRUE;
 201       nd++;
 202     }
 203   }
 204   fprintf(stderr,"There are %d donors and %d acceptors\n",nd,na);
 205   chk_allhb(pdba,x,hb,donor,acceptor,dist);
 206   if (debug)
 207     pr_hbonds(debug,hb,pdba);
 208   fprintf(stderr,"There are %d hydrogen bonds\n",hb->nra);
 209
 210   /* Now do the HIS stuff */
 211   hisnr=-1;
 212   for(i=0; (i<natom); ) {
 213     if (strcasecmp(*pdba->resname[pdba->atom[i].resnr],"HIS") != 0)
 214       i++;
 215     else {
 216       if (pdba->atom[i].resnr != hisnr) {
 217         hisnr=pdba->atom[i].resnr;
 218
 219         /* Find the  atoms in the ring */
 220         nd1=ne2=cg=cd2=ce1=-1;
 221         for(j=i; (pdba->atom[j].resnr==hisnr) && (j<natom); j++) {
 222           atomnm=*pdba->atomname[j];
 223           if (strcmp(atomnm,"CD2") == 0)
 224             cd2=j;
 225           else if (strcmp(atomnm,"CG") == 0)
 226             cg=j;
 227           else if (strcmp(atomnm,"CE1") == 0)
 228             ce1=j;
 229           else if (strcmp(atomnm,"ND1") == 0)
 230             nd1=j;
 231           else if (strcmp(atomnm,"NE2") == 0)
 232             ne2=j;
 233         }
 234
 235         if (!((cg == -1 ) || (cd2 == -1) || (ce1 == -1) ||
 236               (nd1 == -1) || (ne2 == -1))) {
 237           calc_ringh(x[nd1],x[cg],x[ce1],xh1);
 238           calc_ringh(x[ne2],x[ce1],x[cd2],xh2);
 239
 240           bHDd = chk_hbonds(nd1,pdba,x,acceptor,hbond,xh1,angle,dist);
 241           chk_hbonds(nd1,pdba,x,donor,hbond,xh1,angle,dist);
 242           bHEd = chk_hbonds(ne2,pdba,x,acceptor,hbond,xh2,angle,dist);
 243           chk_hbonds(ne2,pdba,x,donor,hbond,xh2,angle,dist);
 244
 245           if (bHDd) {
 246             if (bHEd)
 247               type = ehisH;
 248             else
 249               type = ehisA;
 250           }
 251           else
 252             type = ehisB;
 253           fprintf(stderr,"Will use %s for residue %d\n",hh[type],hisnr+1);
 254         }
 255         else
 256           fatal_error(0,"Incomplete ring in HIS%d",hisnr+1);
 257
 258         sfree(*pdba->resname[hisnr]);
 259         *pdba->resname[hisnr]=strdup(hh[type]);
 260       }
 261     }
 262   }
 263   done_block(hb);
 264   sfree(hb);
 265   sfree(donor);
 266   sfree(acceptor);
 267   sfree(hbond);
 268 }