src/kernel/specbond.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_specbond_c = "$Id$";
  30
  31 #include <ctype.h>
  32 #include <math.h>
  33 #include "typedefs.h"
  34 #include "pdbio.h"
  35 #include "strdb.h"
  36 #include "string2.h"
  37 #include "smalloc.h"
  38 #include "specbond.h"
  39 #include "pdb2top.h"
  40 #include "vec.h"
  41
  42 bool yesno(void)
  43 {
  44   char c;
  45
  46   do {
  47     c=toupper(fgetc(stdin));
  48   } while ((c != 'Y') && (c != 'N'));
  49
  50   return (c == 'Y');
  51 }
  52
  53 typedef struct {
  54   char *res1, *res2;
  55   char *atom1,*atom2;
  56   char *nres1,*nres2;
  57   int  nbond1,nbond2;
  58   real length;
  59 } t_specbond;
  60
  61 static t_specbond *get_specbonds(int *nspecbond)
  62 {
  63   static char  *sbfile="specbond.dat";
  64
  65   t_specbond *sb=NULL;
  66   char   r1buf[32],r2buf[32],a1buf[32],a2buf[32],nr1buf[32],nr2buf[32];
  67   double length;
  68   int    nb1,nb2;
  69   char   **lines;
  70   int    nlines,i,n;
  71
  72   nlines = get_lines(sbfile,&lines);
  73   if (nlines > 0)
  74     snew(sb,nlines);
  75
  76   n = 0;
  77   for(i=0; (i<nlines); i++) {
  78     if (sscanf(lines[i],"%s%s%d%s%s%d%lf%s%s",
  79                r1buf,a1buf,&nb1,r2buf,a2buf,&nb2,&length,nr1buf,nr2buf) != 9)
  80       fprintf(stderr,"Invalid line '%s' in %s\n",lines[i],sbfile);
  81     else {
  82       sb[n].res1   = strdup(r1buf);
  83       sb[n].res2   = strdup(r2buf);
  84       sb[n].nres1  = strdup(nr1buf);
  85       sb[n].nres2  = strdup(nr2buf);
  86       sb[n].atom1  = strdup(a1buf);
  87       sb[n].atom2  = strdup(a2buf);
  88       sb[n].nbond1 = nb1;
  89       sb[n].nbond2 = nb2;
  90       sb[n].length = length;
  91       n++;
  92     }
  93     sfree(lines[i]);
  94   }
  95   if (nlines > 0)
  96     sfree(lines);
  97   fprintf(stderr,"%d out of %d lines of %s converted succesfully\n",
  98           n,nlines,sbfile);
  99
 100   *nspecbond = n;
 101
 102   return sb;
 103 }
 104
 105 static void done_specbonds(int nsb,t_specbond sb[])
 106 {
 107   int i;
 108
 109   for(i=0; (i<nsb); i++) {
 110     sfree(sb[i].res1);
 111     sfree(sb[i].res2);
 112     sfree(sb[i].atom1);
 113     sfree(sb[i].atom2);
 114     sfree(sb[i].nres1);
 115     sfree(sb[i].nres2);
 116   }
 117 }
 118
 119 static bool is_special(int nsb,t_specbond sb[],char *res,char *atom)
 120 {
 121   int i;
 122
 123   for(i=0; (i<nsb); i++) {
 124     if (((strcasecmp(sb[i].res1,res) == 0) &&
 125          (strcasecmp(sb[i].atom1,atom) == 0)) ||
 126         ((strcasecmp(sb[i].res2,res) == 0) &&
 127          (strcasecmp(sb[i].atom2,atom) == 0)))
 128       return TRUE;
 129   }
 130   return FALSE;
 131 }
 132
 133 static bool is_bond(int nsb,t_specbond sb[],t_atoms *pdba,int a1,int a2,
 134                     real d,int *nb,bool *bSwap)
 135 {
 136   int i;
 137   char *at1,*at2,*res1,*res2;
 138
 139   at1=*pdba->atomname[a1];
 140   at2=*pdba->atomname[a2];
 141   res1=*pdba->resname[pdba->atom[a1].resnr];
 142   res2=*pdba->resname[pdba->atom[a2].resnr];
 143
 144   for(i=0; (i<nsb); i++) {
 145     *nb = i;
 146     if (((strcasecmp(sb[i].res1,res1) == 0)  &&
 147          (strcasecmp(sb[i].atom1,at1) == 0) &&
 148          (strcasecmp(sb[i].res2,res2) == 0)  &&
 149          (strcasecmp(sb[i].atom2,at2) == 0))) {
 150       *bSwap = FALSE;
 151       if ((0.9*sb[i].length < d) && (1.1*sb[i].length > d))
 152         return TRUE;
 153     }
 154     if (((strcasecmp(sb[i].res1,res2) == 0)  &&
 155          (strcasecmp(sb[i].atom1,at2) == 0) &&
 156          (strcasecmp(sb[i].res2,res1) == 0)  &&
 157          (strcasecmp(sb[i].atom2,at1) == 0))) {
 158       *bSwap = TRUE;
 159       if ((0.9*sb[i].length < d) && (1.1*sb[i].length > d))
 160         return TRUE;
 161     }
 162   }
 163   return FALSE;
 164 }
 165
 166 static void rename_1res(t_atoms *pdba,int resnr,char *nres)
 167 {
 168   sfree(*pdba->resname[resnr]);
 169   *pdba->resname[resnr]=strdup(nres);
 170 }
 171
 172 int mk_specbonds(t_atoms *pdba,rvec x[],bool bInteractive,
 173                  t_ssbond **specbonds)
 174 {
 175   t_specbond *sb=NULL;
 176   t_ssbond   *bonds=NULL;
 177   int  nsb,natoms;
 178   int  nspec,nbonds;
 179   int  *specp,*sgp;
 180   bool bDoit,bSwap;
 181   int  i,j,b,e,e2;
 182   int  ai,aj,index_sb;
 183   real **d;
 184   char buf[10];
 185
 186   natoms=pdba->nr;
 187   /* Initiate variables that will be exported */
 188   nbonds = 0;
 189   sb=get_specbonds(&nsb);
 190
 191   if (nsb > 0) {
 192     snew(specp,natoms);
 193     snew(sgp,natoms);
 194
 195     nspec = 0;
 196     for(i=0;(i<natoms);i++) {
 197       if (is_special(nsb,sb,*pdba->resname[pdba->atom[i].resnr],
 198                      *pdba->atomname[i])) {
 199         specp[nspec] = pdba->atom[i].resnr;
 200         sgp[nspec] = i;
 201         nspec++;
 202       }
 203     }
 204     srenew(specp,nspec);
 205     srenew(sgp,nspec);
 206
 207     /* distance matrix d[nspec][nspec] */
 208     snew(d,nspec);
 209     for(i=0; (i<nspec); i++)
 210       snew(d[i],nspec);
 211
 212     for(i=0; (i<nspec); i++) {
 213       ai=sgp[i];
 214       for(j=0; (j<nspec); j++) {
 215         aj=sgp[j];
 216         d[i][j]=sqrt(distance2(x[ai],x[aj]));
 217       }
 218     }
 219     if (nspec > 1) {
 220 #define MAXCOL 8
 221       fprintf(stderr,"Special Atom Distance matrix:\n");
 222       for(b=0; (b<nspec); b+=MAXCOL) {
 223         fprintf(stderr,"%8s","");
 224         e=min(b+MAXCOL, nspec-1);
 225         for(i=b; (i<e); i++) {
 226           sprintf(buf,"%s%d",*pdba->resname[pdba->atom[sgp[i]].resnr],
 227                   specp[i]+1);
 228           fprintf(stderr,"%8s",buf);
 229         }
 230         fprintf(stderr,"\n");
 231         e=min(b+MAXCOL, nspec);
 232         for(i=b+1; (i<nspec); i++) {
 233           sprintf(buf,"%s%d",*pdba->resname[pdba->atom[sgp[i]].resnr],
 234                   specp[i]+1);
 235           fprintf(stderr,"%8s",buf);
 236           e2=min(i,e);
 237           for(j=b; (j<e2); j++)
 238           fprintf(stderr," %7.3f",d[i][j]);
 239           fprintf(stderr,"\n");
 240         }
 241       }
 242     }
 243     i=j=0;
 244
 245     snew(bonds,nspec/2);
 246
 247     for(i=0; (i<nspec); i++) {
 248       ai = sgp[i];
 249       for(j=i+1; (j<nspec); j++) {
 250         aj = sgp[j];
 251         if (is_bond(nsb,sb,pdba,ai,aj,d[i][j],
 252                     &index_sb,&bSwap)) {
 253           if (bInteractive) {
 254             fprintf(stderr,"Link Res%4d and Res%4d (y/n) ?",specp[i]+1,specp[j]+1);
 255             bDoit=yesno();
 256           }
 257           else {
 258             fprintf(stderr,"Linking Res%4d and Res%4d...\n",specp[i]+1,specp[j]+1);
 259             bDoit=TRUE;
 260           }
 261           if (bDoit) {
 262             /* Store the residue numbers in the bonds array */
 263             bonds[nbonds].res1 = specp[i];
 264             bonds[nbonds].res2 = specp[j];
 265             bonds[nbonds].a1   = strdup(*pdba->atomname[ai]);
 266             bonds[nbonds].a2   = strdup(*pdba->atomname[aj]);
 267             if (bSwap) {
 268               rename_1res(pdba,specp[i],sb[index_sb].nres2);
 269               rename_1res(pdba,specp[j],sb[index_sb].nres1);
 270             }
 271             else {
 272               rename_1res(pdba,specp[i],sb[index_sb].nres1);
 273               rename_1res(pdba,specp[j],sb[index_sb].nres2);
 274             }
 275             nbonds++;
 276           }
 277         }
 278       }
 279     }
 280
 281     for(i=0; (i<nspec); i++)
 282       sfree(d[i]);
 283     sfree(d);
 284     sfree(sgp);
 285     sfree(specp);
 286
 287     done_specbonds(nsb,sb);
 288     sfree(sb);
 289   }
 290
 291   *specbonds=bonds;
 292
 293   return nbonds;
 294 }
 295