src/kernel/hizzie.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 /* This file is completely threadsafe - keep it that way! */
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  42
  43 #include <stdio.h>
  44 #include <string.h>
  45 #include "typedefs.h"
  46 #include "pdbio.h"
  47 #include "smalloc.h"
  48 #include "vec.h"
  49 #include "physics.h"
  50 #include "toputil.h"
  51 #include "pdb2top.h"
  52 #include "string2.h"
  53 #include "hizzie.h"
  54
  55 static int in_strings(char *key,int nstr,const char **str)
  56 {
  57   int j;
  58
  59   for(j=0; (j<nstr); j++)
  60     if (strcmp(str[j],key) == 0)
  61       return j;
  62
  63   return -1;
  64 }
  65
  66 static gmx_bool hbond(rvec x[],int i,int j,real distance)
  67 {
  68   real tol = distance*distance;
  69   rvec   tmp;
  70
  71   rvec_sub(x[i],x[j],tmp);
  72
  73   return (iprod(tmp,tmp) < tol);
  74 }
  75
  76 static void chk_allhb(t_atoms *pdba,rvec x[],t_blocka *hb,
  77                       gmx_bool donor[],gmx_bool accept[],real dist)
  78 {
  79   int i,j,k,ii,natom;
  80
  81   natom=pdba->nr;
  82   snew(hb->index,natom+1);
  83   snew(hb->a,6*natom);
  84   hb->nr  = natom;
  85   hb->nra = 6*natom;
  86
  87   k = ii = 0;
  88   hb->index[ii++] = 0;
  89   for(i=0; (i<natom); i++) {
  90     if (donor[i]) {
  91       for(j=i+1; (j<natom); j++)
  92         if ((accept[j]) && (hbond(x,i,j,dist)))
  93           hb->a[k++] = j;
  94     }
  95     else if (accept[i]) {
  96       for(j=i+1; (j<natom); j++)
  97         if ((donor[j]) && (hbond(x,i,j,dist)))
  98           hb->a[k++] = j;
  99     }
 100     hb->index[ii++] = k;
 101   }
 102   hb->nra = k;
 103 }
 104
 105 static void pr_hbonds(FILE *fp,t_blocka *hb,t_atoms *pdba)
 106 {
 107   int i,j,k,j0,j1;
 108
 109   fprintf(fp,"Dumping all hydrogen bonds!\n");
 110   for(i=0; (i<hb->nr); i++) {
 111     j0=hb->index[i];
 112     j1=hb->index[i+1];
 113     for(j=j0; (j<j1); j++) {
 114       k=hb->a[j];
 115       fprintf(fp,"%5s%4d%5s - %5s%4d%5s\n",
 116               *pdba->resinfo[pdba->atom[i].resind].name,
 117               pdba->resinfo[pdba->atom[i].resind].nr,*pdba->atomname[i],
 118               *pdba->resinfo[pdba->atom[k].resind].name,
 119               pdba->resinfo[pdba->atom[k].resind].nr,*pdba->atomname[k]);
 120     }
 121   }
 122 }
 123
 124 static gmx_bool chk_hbonds(int i,t_atoms *pdba, rvec x[],
 125                        gmx_bool ad[],gmx_bool hbond[],rvec xh,
 126                        real angle,real dist)
 127 {
 128   gmx_bool bHB;
 129   int  j,aj,ri,natom;
 130   real d2,dist2,a;
 131   rvec nh,oh;
 132
 133   natom=pdba->nr;
 134   bHB = FALSE;
 135   ri = pdba->atom[i].resind;
 136   dist2=sqr(dist);
 137   for(j=0; (j<natom); j++) {
 138     /* Check whether the other atom is a donor/acceptor and not i */
 139     if ((ad[j]) && (j != i)) {
 140       /* Check whether the other atom is on the same ring as well */
 141       if ((pdba->atom[j].resind != ri) ||
 142           ((strcmp(*pdba->atomname[j],"ND1") != 0) &&
 143            (strcmp(*pdba->atomname[j],"NE2") != 0))) {
 144         aj = j;
 145         d2  = distance2(x[i],x[j]);
 146         rvec_sub(x[i],xh,nh);
 147         rvec_sub(x[aj],xh,oh);
 148         a  = RAD2DEG * acos(cos_angle(nh,oh));
 149         if ((d2 < dist2) && (a > angle)) {
 150           if (debug)
 151             fprintf(debug,
 152                     "HBOND between %s%d-%s and %s%d-%s is %g nm, %g deg\n",
 153                     *pdba->resinfo[pdba->atom[i].resind].name,
 154                     pdba->resinfo[pdba->atom[i].resind].nr,*pdba->atomname[i],
 155                     *pdba->resinfo[pdba->atom[aj].resind].name,
 156                     pdba->resinfo[pdba->atom[aj].resind].nr,*pdba->atomname[aj],
 157                     sqrt(d2),a);
 158           hbond[i] = TRUE;
 159           bHB      = TRUE;
 160         }
 161       }
 162     }
 163   }
 164   return bHB;
 165 }
 166
 167 static void calc_ringh(rvec xattach,rvec xb,rvec xc,rvec xh)
 168 {
 169   rvec tab,tac;
 170   real n;
 171
 172   /* Add a proton on a ring to atom attach at distance 0.1 nm */
 173   rvec_sub(xattach,xb,tab);
 174   rvec_sub(xattach,xc,tac);
 175   rvec_add(tab,tac,xh);
 176   n=0.1/norm(xh);
 177   svmul(n,xh,xh);
 178   rvec_inc(xh,xattach);
 179 }
 180
 181 void set_histp(t_atoms *pdba,rvec *x,real angle,real dist){
 182   static const char *prot_acc[] = {
 183     "O", "OD1", "OD2", "OE1", "OE2", "OG", "OG1", "OH", "OW"
 184   };
 185 #define NPA asize(prot_acc)
 186   static const char *prot_don[] = {
 187     "N", "NH1", "NH2", "NE", "ND1", "ND2", "NE2", "NZ", "OG", "OG1", "OH", "NE1", "OW"
 188   };
 189 #define NPD asize(prot_don)
 190
 191   gmx_bool *donor,*acceptor;
 192   gmx_bool *hbond,bHaveH=FALSE;
 193   gmx_bool bHDd,bHEd;
 194   rvec xh1,xh2;
 195   int  natom;
 196   int  i,j,nd,na,aj,hisind,his0,type=-1;
 197   int  nd1,ne2,cg,cd2,ce1;
 198   t_blocka *hb;
 199   real d;
 200   char *atomnm;
 201
 202   natom=pdba->nr;
 203
 204   i = 0;
 205   while (i < natom &&
 206          gmx_strcasecmp(*pdba->resinfo[pdba->atom[i].resind].name,"HIS") != 0)
 207   {
 208     i++;
 209   }
 210   if (natom == i)
 211   {
 212     return;
 213   }
 214
 215   /* A histidine residue exists that requires automated assignment, so
 216    * doing the analysis of donors and acceptors is worthwhile. */
 217   fprintf(stderr,
 218           "Analysing hydrogen-bonding network for automated assignment of histidine\n"
 219           " protonation.");
 220
 221   snew(donor,natom);
 222   snew(acceptor,natom);
 223   snew(hbond,natom);
 224   snew(hb,1);
 225
 226   nd=na=0;
 227   for(j=0; (j<natom); j++) {
 228     if (in_strings(*pdba->atomname[j],NPA,prot_acc) != -1) {
 229       acceptor[j] = TRUE;
 230       na++;
 231     }
 232     if (in_strings(*pdba->atomname[j],NPD,prot_don) != -1) {
 233       donor[j] = TRUE;
 234       nd++;
 235     }
 236   }
 237   fprintf(stderr," %d donors and %d acceptors were found.\n",nd,na);
 238   chk_allhb(pdba,x,hb,donor,acceptor,dist);
 239   if (debug)
 240     pr_hbonds(debug,hb,pdba);
 241   fprintf(stderr,"There are %d hydrogen bonds\n",hb->nra);
 242
 243   /* Now do the HIS stuff */
 244   hisind=-1;
 245   while (i < natom)
 246   {
 247     if (gmx_strcasecmp(*pdba->resinfo[pdba->atom[i].resind].name,"HIS") != 0)
 248     {
 249       i++;
 250     }
 251     else
 252     {
 253       if (pdba->atom[i].resind != hisind) {
 254         hisind=pdba->atom[i].resind;
 255
 256         /* Find the  atoms in the ring */
 257         nd1=ne2=cg=cd2=ce1=-1;
 258         while (i<natom && pdba->atom[i].resind==hisind) {
 259           atomnm = *pdba->atomname[i];
 260           if (strcmp(atomnm,"CD2") == 0)
 261             cd2 = i;
 262           else if (strcmp(atomnm,"CG") == 0)
 263             cg  = i;
 264           else if (strcmp(atomnm,"CE1") == 0)
 265             ce1 = i;
 266           else if (strcmp(atomnm,"ND1") == 0)
 267             nd1 = i;
 268           else if (strcmp(atomnm,"NE2") == 0)
 269             ne2 = i;
 270
 271           i++;
 272         }
 273
 274         if (!((cg == -1 ) || (cd2 == -1) || (ce1 == -1) ||
 275               (nd1 == -1) || (ne2 == -1))) {
 276           calc_ringh(x[nd1],x[cg],x[ce1],xh1);
 277           calc_ringh(x[ne2],x[ce1],x[cd2],xh2);
 278
 279           bHDd = chk_hbonds(nd1,pdba,x,acceptor,hbond,xh1,angle,dist);
 280           chk_hbonds(nd1,pdba,x,donor,hbond,xh1,angle,dist);
 281           bHEd = chk_hbonds(ne2,pdba,x,acceptor,hbond,xh2,angle,dist);
 282           chk_hbonds(ne2,pdba,x,donor,hbond,xh2,angle,dist);
 283
 284           if (bHDd) {
 285             if (bHEd)
 286               type = ehisH;
 287             else
 288               type = ehisA;
 289           }
 290           else
 291             type = ehisB;
 292           fprintf(stderr,"Will use %s for residue %d\n",
 293                   hh[type],pdba->resinfo[hisind].nr);
 294         }
 295         else
 296           gmx_fatal(FARGS,"Incomplete ring in HIS%d",
 297                     pdba->resinfo[hisind].nr);
 298
 299         snew(pdba->resinfo[hisind].rtp,1);
 300         *pdba->resinfo[hisind].rtp = strdup(hh[type]);
 301       }
 302     }
 303   }
 304   done_blocka(hb);
 305   sfree(hb);
 306   sfree(donor);
 307   sfree(acceptor);
 308   sfree(hbond);
 309 }