src/kernel/specbond.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 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <ctype.h>
  40 #include <math.h>
  41 #include "typedefs.h"
  42 #include "pdbio.h"
  43 #include "strdb.h"
  44 #include "string2.h"
  45 #include "smalloc.h"
  46 #include "specbond.h"
  47 #include "pdb2top.h"
  48 #include "vec.h"
  49
  50 gmx_bool yesno(void)
  51 {
  52   char c;
  53
  54   do {
  55     c=toupper(fgetc(stdin));
  56   } while ((c != 'Y') && (c != 'N'));
  57
  58   return (c == 'Y');
  59 }
  60
  61 t_specbond *get_specbonds(int *nspecbond)
  62 {
  63   const 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].newres1= strdup(nr1buf);
  85       sb[n].newres2= 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 successfully\n",
  98           n,nlines,sbfile);
  99
 100   *nspecbond = n;
 101
 102   return sb;
 103 }
 104
 105 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].newres1);
 115     sfree(sb[i].newres2);
 116   }
 117 }
 118
 119 static gmx_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 (((strncmp(sb[i].res1,res,3) == 0) &&
 125          (gmx_strcasecmp(sb[i].atom1,atom) == 0)) ||
 126         ((strncmp(sb[i].res2,res,3) == 0) &&
 127          (gmx_strcasecmp(sb[i].atom2,atom) == 0)))
 128       return TRUE;
 129   }
 130   return FALSE;
 131 }
 132
 133 static gmx_bool is_bond(int nsb,t_specbond sb[],t_atoms *pdba,int a1,int a2,
 134                     real d,int *index_sb,gmx_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->resinfo[pdba->atom[a1].resind].name;
 142   res2=*pdba->resinfo[pdba->atom[a2].resind].name;
 143
 144   if (debug)
 145     fprintf(stderr,"Checking %s-%d %s-%d and %s-%d %s-%d: %g ",
 146             res1, pdba->resinfo[pdba->atom[a1].resind].nr, at1, a1+1,
 147             res2, pdba->resinfo[pdba->atom[a2].resind].nr, at2, a2+1, d);
 148
 149   for(i=0; (i<nsb); i++) {
 150     *index_sb = i;
 151     if (((strncmp(sb[i].res1,res1,3) == 0)  &&
 152          (gmx_strcasecmp(sb[i].atom1,at1) == 0) &&
 153          (strncmp(sb[i].res2,res2,3) == 0)  &&
 154          (gmx_strcasecmp(sb[i].atom2,at2) == 0))) {
 155       *bSwap = FALSE;
 156       if ((0.9*sb[i].length < d) && (1.1*sb[i].length > d)) {
 157         if (debug) fprintf(stderr,"%g\n", sb[i].length);
 158         return TRUE;
 159       }
 160     }
 161     if (((strncmp(sb[i].res1,res2,3) == 0)  &&
 162          (gmx_strcasecmp(sb[i].atom1,at2) == 0) &&
 163          (strncmp(sb[i].res2,res1,3) == 0)  &&
 164          (gmx_strcasecmp(sb[i].atom2,at1) == 0))) {
 165       *bSwap = TRUE;
 166       if ((0.9*sb[i].length < d) && (1.1*sb[i].length > d)) {
 167         if (debug) fprintf(stderr,"%g\n", sb[i].length);
 168         return TRUE;
 169       }
 170     }
 171   }
 172   if (debug) fprintf(stderr,"\n");
 173   return FALSE;
 174 }
 175
 176 static void rename_1res(t_atoms *pdba,int resind,char *newres,gmx_bool bVerbose)
 177 {
 178   if (bVerbose) {
 179     printf("Using rtp entry %s for %s %d\n",
 180            newres,
 181            *pdba->resinfo[resind].name,
 182            pdba->resinfo[resind].nr);
 183   }
 184   /* this used to free *resname, which fucks up the symtab! */
 185   snew(pdba->resinfo[resind].rtp,1);
 186   *pdba->resinfo[resind].rtp = strdup(newres);
 187 }
 188
 189 int mk_specbonds(t_atoms *pdba,rvec x[],gmx_bool bInteractive,
 190                  t_ssbond **specbonds,gmx_bool bVerbose)
 191 {
 192   t_specbond *sb=NULL;
 193   t_ssbond   *bonds=NULL;
 194   int  nsb;
 195   int  nspec,nbonds;
 196   int  *specp,*sgp;
 197   gmx_bool bDoit,bSwap;
 198   int  i,j,b,e,e2;
 199   int  ai,aj,index_sb;
 200   real **d;
 201   char buf[10];
 202
 203   nbonds = 0;
 204   sb     = get_specbonds(&nsb);
 205
 206   if (nsb > 0) {
 207     snew(specp,pdba->nr);
 208     snew(sgp,pdba->nr);
 209
 210     nspec = 0;
 211     for(i=0;(i<pdba->nr);i++) {
 212       /* Check if this atom is special and if it is not a double atom
 213        * in the input that still needs to be removed.
 214        */
 215       if (is_special(nsb,sb,*pdba->resinfo[pdba->atom[i].resind].name,
 216                      *pdba->atomname[i]) &&
 217           !(nspec > 0 &&
 218             pdba->atom[sgp[nspec-1]].resind == pdba->atom[i].resind &&
 219             gmx_strcasecmp(*pdba->atomname[sgp[nspec-1]],
 220                        *pdba->atomname[i]) == 0)) {
 221         specp[nspec] = pdba->atom[i].resind;
 222         sgp[nspec] = i;
 223         nspec++;
 224       }
 225     }
 226     /* distance matrix d[nspec][nspec] */
 227     snew(d,nspec);
 228     for(i=0; (i<nspec); i++)
 229       snew(d[i],nspec);
 230
 231     for(i=0; (i<nspec); i++) {
 232       ai=sgp[i];
 233       for(j=0; (j<nspec); j++) {
 234         aj=sgp[j];
 235         d[i][j]=sqrt(distance2(x[ai],x[aj]));
 236       }
 237     }
 238     if (nspec > 1) {
 239 #define MAXCOL 7
 240       fprintf(stderr,"Special Atom Distance matrix:\n");
 241       for(b=0; (b<nspec); b+=MAXCOL) {
 242         /* print resname/number column headings */
 243         fprintf(stderr,"%8s%8s","","");
 244         e=min(b+MAXCOL, nspec-1);
 245         for(i=b; (i<e); i++) {
 246           sprintf(buf,"%s%d",*pdba->resinfo[pdba->atom[sgp[i]].resind].name,
 247                   pdba->resinfo[specp[i]].nr);
 248           fprintf(stderr,"%8s",buf);
 249         }
 250         fprintf(stderr,"\n");
 251         /* print atomname/number column headings */
 252         fprintf(stderr,"%8s%8s","","");
 253         e=min(b+MAXCOL, nspec-1);
 254         for(i=b; (i<e); i++) {
 255           sprintf(buf,"%s%d",*pdba->atomname[sgp[i]], sgp[i]+1);
 256           fprintf(stderr,"%8s",buf);
 257         }
 258         fprintf(stderr,"\n");
 259         /* print matrix */
 260         e=min(b+MAXCOL, nspec);
 261         for(i=b+1; (i<nspec); i++) {
 262           sprintf(buf,"%s%d",*pdba->resinfo[pdba->atom[sgp[i]].resind].name,
 263                   pdba->resinfo[specp[i]].nr);
 264           fprintf(stderr,"%8s",buf);
 265           sprintf(buf,"%s%d", *pdba->atomname[sgp[i]], sgp[i]+1);
 266           fprintf(stderr,"%8s",buf);
 267           e2=min(i,e);
 268           for(j=b; (j<e2); j++)
 269             fprintf(stderr," %7.3f",d[i][j]);
 270           fprintf(stderr,"\n");
 271         }
 272       }
 273     }
 274
 275     snew(bonds,nspec);
 276
 277     for(i=0; (i<nspec); i++) {
 278       ai = sgp[i];
 279       for(j=i+1; (j<nspec); j++) {
 280         aj = sgp[j];
 281         if (is_bond(nsb,sb,pdba,ai,aj,d[i][j],&index_sb,&bSwap)) {
 282           fprintf(stderr,"%s %s-%d %s-%d and %s-%d %s-%d%s",
 283                   bInteractive ? "Link" : "Linking",
 284                   *pdba->resinfo[pdba->atom[ai].resind].name,
 285                   pdba->resinfo[specp[i]].nr,
 286                   *pdba->atomname[ai], ai+1,
 287                   *pdba->resinfo[pdba->atom[aj].resind].name,
 288                   pdba->resinfo[specp[j]].nr,
 289                   *pdba->atomname[aj], aj+1,
 290                   bInteractive ? " (y/n) ?" : "...\n");
 291           bDoit=bInteractive ? yesno() : TRUE;
 292
 293           if (bDoit) {
 294             /* Store the residue numbers in the bonds array */
 295             bonds[nbonds].res1 = specp[i];
 296             bonds[nbonds].res2 = specp[j];
 297             bonds[nbonds].a1   = strdup(*pdba->atomname[ai]);
 298             bonds[nbonds].a2   = strdup(*pdba->atomname[aj]);
 299             /* rename residues */
 300             if (bSwap) {
 301               rename_1res(pdba,specp[i],sb[index_sb].newres2,bVerbose);
 302               rename_1res(pdba,specp[j],sb[index_sb].newres1,bVerbose);
 303             }
 304             else {
 305               rename_1res(pdba,specp[i],sb[index_sb].newres1,bVerbose);
 306               rename_1res(pdba,specp[j],sb[index_sb].newres2,bVerbose);
 307             }
 308             nbonds++;
 309           }
 310         }
 311       }
 312     }
 313
 314     for(i=0; (i<nspec); i++)
 315       sfree(d[i]);
 316     sfree(d);
 317     sfree(sgp);
 318     sfree(specp);
 319
 320     done_specbonds(nsb,sb);
 321     sfree(sb);
 322   }
 323
 324   *specbonds=bonds;
 325
 326   return nbonds;
 327 }
 328