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