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