src/programs/pdb2gmx/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 gmx_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                       gmx_bool donor[],gmx_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 gmx_bool chk_hbonds(int i,t_atoms *pdba, rvec x[],
 121                        gmx_bool ad[],gmx_bool hbond[],rvec xh,
 122                        real angle,real dist)
 123 {
 124   gmx_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   gmx_bool *donor,*acceptor;
 188   gmx_bool *hbond,bHaveH=FALSE;
 189   gmx_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
 200   i = 0;
 201   while (i < natom &&
 202          gmx_strcasecmp(*pdba->resinfo[pdba->atom[i].resind].name,"HIS") != 0)
 203   {
 204     i++;
 205   }
 206   if (natom == i)
 207   {
 208     return;
 209   }
 210
 211   /* A histidine residue exists that requires automated assignment, so
 212    * doing the analysis of donors and acceptors is worthwhile. */
 213   fprintf(stderr,
 214           "Analysing hydrogen-bonding network for automated assigment of histidine\n"
 215           " protonation.");
 216
 217   snew(donor,natom);
 218   snew(acceptor,natom);
 219   snew(hbond,natom);
 220   snew(hb,1);
 221
 222   nd=na=0;
 223   for(j=0; (j<natom); j++) {
 224     if (in_strings(*pdba->atomname[j],NPA,prot_acc) != -1) {
 225       acceptor[j] = TRUE;
 226       na++;
 227     }
 228     if (in_strings(*pdba->atomname[j],NPD,prot_don) != -1) {
 229       donor[j] = TRUE;
 230       nd++;
 231     }
 232   }
 233   fprintf(stderr," %d donors and %d acceptors were found.\n",nd,na);
 234   chk_allhb(pdba,x,hb,donor,acceptor,dist);
 235   if (debug)
 236     pr_hbonds(debug,hb,pdba);
 237   fprintf(stderr,"There are %d hydrogen bonds\n",hb->nra);
 238
 239   /* Now do the HIS stuff */
 240   hisind=-1;
 241   while (i < natom)
 242   {
 243     if (gmx_strcasecmp(*pdba->resinfo[pdba->atom[i].resind].name,"HIS") != 0)
 244     {
 245       i++;
 246     }
 247     else
 248     {
 249       if (pdba->atom[i].resind != hisind) {
 250         hisind=pdba->atom[i].resind;
 251
 252         /* Find the  atoms in the ring */
 253         nd1=ne2=cg=cd2=ce1=-1;
 254         while (i<natom && pdba->atom[i].resind==hisind) {
 255           atomnm = *pdba->atomname[i];
 256           if (strcmp(atomnm,"CD2") == 0)
 257             cd2 = i;
 258           else if (strcmp(atomnm,"CG") == 0)
 259             cg  = i;
 260           else if (strcmp(atomnm,"CE1") == 0)
 261             ce1 = i;
 262           else if (strcmp(atomnm,"ND1") == 0)
 263             nd1 = i;
 264           else if (strcmp(atomnm,"NE2") == 0)
 265             ne2 = i;
 266
 267           i++;
 268         }
 269
 270         if (!((cg == -1 ) || (cd2 == -1) || (ce1 == -1) ||
 271               (nd1 == -1) || (ne2 == -1))) {
 272           calc_ringh(x[nd1],x[cg],x[ce1],xh1);
 273           calc_ringh(x[ne2],x[ce1],x[cd2],xh2);
 274
 275           bHDd = chk_hbonds(nd1,pdba,x,acceptor,hbond,xh1,angle,dist);
 276           chk_hbonds(nd1,pdba,x,donor,hbond,xh1,angle,dist);
 277           bHEd = chk_hbonds(ne2,pdba,x,acceptor,hbond,xh2,angle,dist);
 278           chk_hbonds(ne2,pdba,x,donor,hbond,xh2,angle,dist);
 279
 280           if (bHDd) {
 281             if (bHEd)
 282               type = ehisH;
 283             else
 284               type = ehisA;
 285           }
 286           else
 287             type = ehisB;
 288           fprintf(stderr,"Will use %s for residue %d\n",
 289                   hh[type],pdba->resinfo[hisind].nr);
 290         }
 291         else
 292           gmx_fatal(FARGS,"Incomplete ring in HIS%d",
 293                     pdba->resinfo[hisind].nr);
 294
 295         snew(pdba->resinfo[hisind].rtp,1);
 296         *pdba->resinfo[hisind].rtp = strdup(hh[type]);
 297       }
 298     }
 299   }
 300   done_blocka(hb);
 301   sfree(hb);
 302   sfree(donor);
 303   sfree(acceptor);
 304   sfree(hbond);
 305 }