Project

General

Profile

specbond.c

specbond.c - David van der Spoel, 04/16/2006 10:52 PM

 
1
/*
2
 * $Id: specbond.c,v 1.22 2004/01/23 17:23:16 lindahl Exp $
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 3.2.0
11
 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12
 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13
 * Copyright (c) 2001-2004, The GROMACS development team,
14
 * check out http://www.gromacs.org for more information.
15

16
 * This program is free software; you can redistribute it and/or
17
 * modify it under the terms of the GNU General Public License
18
 * as published by the Free Software Foundation; either version 2
19
 * of the License, or (at your option) any later version.
20
 * 
21
 * If you want to redistribute modifications, please consider that
22
 * scientific software is very special. Version control is crucial -
23
 * bugs must be traceable. We will be happy to consider code for
24
 * inclusion in the official distribution, but derived work must not
25
 * be called official GROMACS. Details are found in the README & COPYING
26
 * files - if they are missing, get the official version at www.gromacs.org.
27
 * 
28
 * To help us fund GROMACS development, we humbly ask that you cite
29
 * the papers on the package - you can find them in the top README file.
30
 * 
31
 * For more info, check our website at http://www.gromacs.org
32
 * 
33
 * And Hey:
34
 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
35
 */
36
#ifdef HAVE_CONFIG_H
37
#include <config.h>
38
#endif
39

    
40
#include <ctype.h>
41
#include <math.h>
42
#include "typedefs.h"
43
#include "pdbio.h"
44
#include "strdb.h"
45
#include "string2.h"
46
#include "smalloc.h"
47
#include "specbond.h"
48
#include "pdb2top.h"
49
#include "vec.h"
50

    
51
bool yesno(void)
52
{
53
  char c;
54

    
55
  do {
56
    c=toupper(fgetc(stdin));
57
  } while ((c != 'Y') && (c != 'N'));
58
  
59
  return (c == 'Y');
60
}
61

    
62
t_specbond *get_specbonds(int *nspecbond)
63
{
64
  const char  *sbfile="specbond.dat";
65
  
66
  t_specbond *sb=NULL;
67
  char   r1buf[32],r2buf[32],a1buf[32],a2buf[32],nr1buf[32],nr2buf[32];
68
  double length;
69
  int    nb1,nb2;
70
  char   **lines;
71
  int    nlines,i,n;
72
  
73
  nlines = get_lines(sbfile,&lines);
74
  if (nlines > 0)
75
    snew(sb,nlines);
76
  
77
  n = 0;
78
  for(i=0; (i<nlines); i++) {
79
    if (sscanf(lines[i],"%s%s%d%s%s%d%lf%s%s",
80
               r1buf,a1buf,&nb1,r2buf,a2buf,&nb2,&length,nr1buf,nr2buf) != 9) 
81
      fprintf(stderr,"Invalid line '%s' in %s\n",lines[i],sbfile);
82
    else {
83
      sb[n].res1   = strdup(r1buf);
84
      sb[n].res2   = strdup(r2buf);
85
      sb[n].newres1= strdup(nr1buf);
86
      sb[n].newres2= strdup(nr2buf);
87
      sb[n].atom1  = strdup(a1buf);
88
      sb[n].atom2  = strdup(a2buf);
89
      sb[n].nbond1 = nb1;
90
      sb[n].nbond2 = nb2;
91
      sb[n].length = length;
92
      n++;
93
    }
94
    sfree(lines[i]);
95
  }
96
  if (nlines > 0)
97
    sfree(lines);
98
  fprintf(stderr,"%d out of %d lines of %s converted succesfully\n",
99
          n,nlines,sbfile);
100
          
101
  *nspecbond = n;
102
  
103
  return sb;
104
}
105

    
106
void done_specbonds(int nsb,t_specbond sb[])
107
{
108
  int i;
109
  
110
  for(i=0; (i<nsb); i++) {
111
    sfree(sb[i].res1);
112
    sfree(sb[i].res2);
113
    sfree(sb[i].atom1);
114
    sfree(sb[i].atom2);
115
    sfree(sb[i].newres1);
116
    sfree(sb[i].newres2);
117
  }
118
}
119

    
120
static bool is_special(int nsb,t_specbond sb[],char *res,char *atom)
121
{
122
  int i;
123
  
124
  for(i=0; (i<nsb); i++) {
125
    if (((strncmp(sb[i].res1,res,3) == 0) && 
126
         (strcasecmp(sb[i].atom1,atom) == 0)) ||
127
        ((strncmp(sb[i].res2,res,3) == 0) && 
128
         (strcasecmp(sb[i].atom2,atom) == 0)))
129
      return TRUE;
130
  }
131
  return FALSE;
132
}
133

    
134
static bool is_bond(int nsb,t_specbond sb[],t_atoms *pdba,int a1,int a2,
135
                    real d,int *index_sb,bool *bSwap)
136
{
137
  int i;
138
  char *at1,*at2,*res1,*res2;
139
  
140
  at1=*pdba->atomname[a1];
141
  at2=*pdba->atomname[a2];
142
  res1=*pdba->resname[pdba->atom[a1].resnr];
143
  res2=*pdba->resname[pdba->atom[a2].resnr];
144
  
145
  if (debug) 
146
    fprintf(stderr,"Checking %s-%d %s-%d and %s-%d %s-%d: %g ",
147
            res1, pdba->atom[a1].resnr+1, at1, a1+1,
148
            res2, pdba->atom[a2].resnr+1, at2, a2+1, d);
149
                    
150
  for(i=0; (i<nsb); i++) {
151
    *index_sb = i;
152
    if (((strncmp(sb[i].res1,res1,3) == 0)  && 
153
         (strcasecmp(sb[i].atom1,at1) == 0) &&
154
         (strncmp(sb[i].res2,res2,3) == 0)  && 
155
         (strcasecmp(sb[i].atom2,at2) == 0))) {
156
      *bSwap = FALSE;
157
      if ((0.9*sb[i].length < d) && (1.1*sb[i].length > d)) {
158
        if (debug) fprintf(stderr,"%g\n", sb[i].length);
159
        return TRUE;
160
      }
161
    }
162
    if (((strncmp(sb[i].res1,res2,3) == 0)  && 
163
         (strcasecmp(sb[i].atom1,at2) == 0) &&
164
         (strncmp(sb[i].res2,res1,3) == 0)  && 
165
         (strcasecmp(sb[i].atom2,at1) == 0))) {
166
      *bSwap = TRUE;
167
      if ((0.9*sb[i].length < d) && (1.1*sb[i].length > d)) {
168
        if (debug) fprintf(stderr,"%g\n", sb[i].length);
169
        return TRUE;
170
      }
171
    }
172
  }
173
  if (debug) fprintf(stderr,"\n");
174
  return FALSE;
175
}
176

    
177
static void rename_1res(t_atoms *pdba,int resnr,char *newres)
178
{
179
  if(debug) fprintf(stderr,"Renaming %s-%d to %s\n", 
180
                    *pdba->resname[resnr], resnr+1, newres);
181
  /* this used to free *resname, which fucks up the symtab! */
182
  snew(pdba->resname[resnr],1);
183
  *pdba->resname[resnr]=strdup(newres);
184
}
185

    
186
int mk_specbonds(t_atoms *pdba,rvec x[],bool bInteractive,
187
                 t_ssbond **specbonds)
188
{
189
  t_specbond *sb=NULL;
190
  t_ssbond   *bonds=NULL;
191
  int  nsb;
192
  int  nspec,nbonds;
193
  int  *specp,*sgp;
194
  bool bDoit,bSwap;
195
  int  i,j,b,e,e2;
196
  int  ai,aj,index_sb;
197
  real **d;
198
  char buf[10];
199

    
200
  nbonds = 0;
201
  sb     = get_specbonds(&nsb);
202
  
203
  if (nsb > 0) {
204
    snew(specp,pdba->nr);
205
    snew(sgp,pdba->nr);
206
    
207
    nspec = 0;
208
    for(i=0;(i<pdba->nr);i++) {
209
      if (is_special(nsb,sb,*pdba->resname[pdba->atom[i].resnr],
210
                     *pdba->atomname[i])) {
211
        specp[nspec] = pdba->atom[i].resnr;
212
        sgp[nspec] = i;
213
        nspec++;
214
      }
215
    }
216
    /* distance matrix d[nspec][nspec] */
217
    snew(d,nspec);
218
    for(i=0; (i<nspec); i++)
219
      snew(d[i],nspec);
220
    
221
    for(i=0; (i<nspec); i++) {
222
      ai=sgp[i];
223
      for(j=0; (j<nspec); j++) {
224
        aj=sgp[j];
225
        d[i][j]=sqrt(distance2(x[ai],x[aj]));
226
      }
227
    }
228
    if (nspec > 1) {
229
#define MAXCOL 7
230
      fprintf(stderr,"Special Atom Distance matrix:\n");
231
      for(b=0; (b<nspec); b+=MAXCOL) {
232
        /* print resname/number column headings */
233
        fprintf(stderr,"%8s%8s","","");
234
        e=min(b+MAXCOL, nspec-1);
235
        for(i=b; (i<e); i++) {
236
          sprintf(buf,"%s%d",*pdba->resname[pdba->atom[sgp[i]].resnr],
237
                  specp[i]+1);
238
          fprintf(stderr,"%8s",buf);
239
        }
240
        fprintf(stderr,"\n");
241
        /* print atomname/number column headings */
242
        fprintf(stderr,"%8s%8s","","");
243
        e=min(b+MAXCOL, nspec-1);
244
        for(i=b; (i<e); i++) {
245
          sprintf(buf,"%s%d",*pdba->atomname[sgp[i]], sgp[i]+1);
246
          fprintf(stderr,"%8s",buf);
247
        }
248
        fprintf(stderr,"\n");
249
        /* print matrix */
250
        e=min(b+MAXCOL, nspec);
251
        for(i=b+1; (i<nspec); i++) {
252
          sprintf(buf,"%s%d",*pdba->resname[pdba->atom[sgp[i]].resnr],
253
                  specp[i]+1);
254
          fprintf(stderr,"%8s",buf);
255
          sprintf(buf,"%s%d", *pdba->atomname[sgp[i]], sgp[i]+1);
256
          fprintf(stderr,"%8s",buf);
257
          e2=min(i,e);
258
          for(j=b; (j<e2); j++)
259
            fprintf(stderr," %7.3f",d[i][j]);
260
          fprintf(stderr,"\n");
261
        }
262
      }
263
    }
264
    
265
    snew(bonds,nspec);
266

    
267
    for(i=0; (i<nspec); i++) {
268
      ai = sgp[i];
269
      for(j=i+1; (j<nspec); j++) {
270
        aj = sgp[j];
271
        if (is_bond(nsb,sb,pdba,ai,aj,d[i][j],&index_sb,&bSwap)) {
272
          fprintf(stderr,"%s %s-%d %s-%d and %s-%d %s-%d%s",
273
                  bInteractive ? "Link" : "Linking",
274
                  *pdba->resname[pdba->atom[ai].resnr], specp[i]+1, 
275
                  *pdba->atomname[ai], ai+1,
276
                  *pdba->resname[pdba->atom[aj].resnr], specp[j]+1, 
277
                  *pdba->atomname[aj], aj+1,
278
                  bInteractive ? " (y/n) ?" : "...\n");
279
          bDoit=bInteractive ? yesno() : TRUE;
280
          
281
          if (bDoit) {
282
            /* Store the residue numbers in the bonds array */
283
            bonds[nbonds].res1 = specp[i];
284
            bonds[nbonds].res2 = specp[j];
285
            bonds[nbonds].a1   = strdup(*pdba->atomname[ai]);
286
            bonds[nbonds].a2   = strdup(*pdba->atomname[aj]);
287
            /* rename residues */
288
            if (bSwap) {
289
              rename_1res(pdba,specp[i],sb[index_sb].newres2);
290
              rename_1res(pdba,specp[j],sb[index_sb].newres1);
291
            }
292
            else {
293
              rename_1res(pdba,specp[i],sb[index_sb].newres1);
294
              rename_1res(pdba,specp[j],sb[index_sb].newres2);
295
            }
296
            nbonds++;
297
          }
298
        }
299
      }
300
    }
301
    
302
    for(i=0; (i<nspec); i++)
303
      sfree(d[i]);
304
    sfree(d);
305
    sfree(sgp);
306
    sfree(specp);
307
 
308
    done_specbonds(nsb,sb);
309
    sfree(sb);
310
  }
311
  
312
  *specbonds=bonds;
313
  
314
  return nbonds;
315
}
316