src/kernel/specbond.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <ctype.h>
  43 #include <math.h>
  44 #include "typedefs.h"
  45 #include "pdbio.h"
  46 #include "strdb.h"
  47 #include "string2.h"
  48 #include "smalloc.h"
  49 #include "specbond.h"
  50 #include "pdb2top.h"
  51 #include "vec.h"
  52
  53 gmx_bool yesno(void)
  54 {
  55   char c;
  56
  57   do {
  58     c=toupper(fgetc(stdin));
  59   } while ((c != 'Y') && (c != 'N'));
  60
  61   return (c == 'Y');
  62 }
  63
  64 t_specbond *get_specbonds(int *nspecbond)
  65 {
  66   const char  *sbfile="specbond.dat";
  67
  68   t_specbond *sb=NULL;
  69   char   r1buf[32],r2buf[32],a1buf[32],a2buf[32],nr1buf[32],nr2buf[32];
  70   double length;
  71   int    nb1,nb2;
  72   char   **lines;
  73   int    nlines,i,n;
  74
  75   nlines = get_lines(sbfile,&lines);
  76   if (nlines > 0)
  77     snew(sb,nlines);
  78
  79   n = 0;
  80   for(i=0; (i<nlines); i++) {
  81     if (sscanf(lines[i],"%s%s%d%s%s%d%lf%s%s",
  82                r1buf,a1buf,&nb1,r2buf,a2buf,&nb2,&length,nr1buf,nr2buf) != 9)
  83       fprintf(stderr,"Invalid line '%s' in %s\n",lines[i],sbfile);
  84     else {
  85       sb[n].res1   = strdup(r1buf);
  86       sb[n].res2   = strdup(r2buf);
  87       sb[n].newres1= strdup(nr1buf);
  88       sb[n].newres2= strdup(nr2buf);
  89       sb[n].atom1  = strdup(a1buf);
  90       sb[n].atom2  = strdup(a2buf);
  91       sb[n].nbond1 = nb1;
  92       sb[n].nbond2 = nb2;
  93       sb[n].length = length;
  94       n++;
  95     }
  96     sfree(lines[i]);
  97   }
  98   if (nlines > 0)
  99     sfree(lines);
 100   fprintf(stderr,"%d out of %d lines of %s converted successfully\n",
 101           n,nlines,sbfile);
 102
 103   *nspecbond = n;
 104
 105   return sb;
 106 }
 107
 108 void done_specbonds(int nsb,t_specbond sb[])
 109 {
 110   int i;
 111
 112   for(i=0; (i<nsb); i++) {
 113     sfree(sb[i].res1);
 114     sfree(sb[i].res2);
 115     sfree(sb[i].atom1);
 116     sfree(sb[i].atom2);
 117     sfree(sb[i].newres1);
 118     sfree(sb[i].newres2);
 119   }
 120 }
 121
 122 static gmx_bool is_special(int nsb,t_specbond sb[],char *res,char *atom)
 123 {
 124   int i;
 125
 126   for(i=0; (i<nsb); i++) {
 127     if (((strncmp(sb[i].res1,res,3) == 0) &&
 128          (gmx_strcasecmp(sb[i].atom1,atom) == 0)) ||
 129         ((strncmp(sb[i].res2,res,3) == 0) &&
 130          (gmx_strcasecmp(sb[i].atom2,atom) == 0)))
 131       return TRUE;
 132   }
 133   return FALSE;
 134 }
 135
 136 static gmx_bool is_bond(int nsb,t_specbond sb[],t_atoms *pdba,int a1,int a2,
 137                     real d,int *index_sb,gmx_bool *bSwap)
 138 {
 139   int i;
 140   char *at1,*at2,*res1,*res2;
 141
 142   at1=*pdba->atomname[a1];
 143   at2=*pdba->atomname[a2];
 144   res1=*pdba->resinfo[pdba->atom[a1].resind].name;
 145   res2=*pdba->resinfo[pdba->atom[a2].resind].name;
 146
 147   if (debug)
 148     fprintf(stderr,"Checking %s-%d %s-%d and %s-%d %s-%d: %g ",
 149             res1, pdba->resinfo[pdba->atom[a1].resind].nr, at1, a1+1,
 150             res2, pdba->resinfo[pdba->atom[a2].resind].nr, at2, a2+1, d);
 151
 152   for(i=0; (i<nsb); i++) {
 153     *index_sb = i;
 154     if (((strncmp(sb[i].res1,res1,3) == 0)  &&
 155          (gmx_strcasecmp(sb[i].atom1,at1) == 0) &&
 156          (strncmp(sb[i].res2,res2,3) == 0)  &&
 157          (gmx_strcasecmp(sb[i].atom2,at2) == 0))) {
 158       *bSwap = FALSE;
 159       if ((0.9*sb[i].length < d) && (1.1*sb[i].length > d)) {
 160         if (debug) fprintf(stderr,"%g\n", sb[i].length);
 161         return TRUE;
 162       }
 163     }
 164     if (((strncmp(sb[i].res1,res2,3) == 0)  &&
 165          (gmx_strcasecmp(sb[i].atom1,at2) == 0) &&
 166          (strncmp(sb[i].res2,res1,3) == 0)  &&
 167          (gmx_strcasecmp(sb[i].atom2,at1) == 0))) {
 168       *bSwap = TRUE;
 169       if ((0.9*sb[i].length < d) && (1.1*sb[i].length > d)) {
 170         if (debug) fprintf(stderr,"%g\n", sb[i].length);
 171         return TRUE;
 172       }
 173     }
 174   }
 175   if (debug) fprintf(stderr,"\n");
 176   return FALSE;
 177 }
 178
 179 static void rename_1res(t_atoms *pdba,int resind,char *newres,gmx_bool bVerbose)
 180 {
 181   if (bVerbose) {
 182     printf("Using rtp entry %s for %s %d\n",
 183            newres,
 184            *pdba->resinfo[resind].name,
 185            pdba->resinfo[resind].nr);
 186   }
 187   /* this used to free *resname, which fucks up the symtab! */
 188   snew(pdba->resinfo[resind].rtp,1);
 189   *pdba->resinfo[resind].rtp = strdup(newres);
 190 }
 191
 192 int mk_specbonds(t_atoms *pdba,rvec x[],gmx_bool bInteractive,
 193                  t_ssbond **specbonds,gmx_bool bVerbose)
 194 {
 195   t_specbond *sb=NULL;
 196   t_ssbond   *bonds=NULL;
 197   int  nsb;
 198   int  nspec,nbonds;
 199   int  *specp,*sgp;
 200   gmx_bool bDoit,bSwap;
 201   int  i,j,b,e,e2;
 202   int  ai,aj,index_sb;
 203   real **d;
 204   char buf[10];
 205
 206   nbonds = 0;
 207   sb     = get_specbonds(&nsb);
 208
 209   if (nsb > 0) {
 210     snew(specp,pdba->nr);
 211     snew(sgp,pdba->nr);
 212
 213     nspec = 0;
 214     for(i=0;(i<pdba->nr);i++) {
 215       /* Check if this atom is special and if it is not a double atom
 216        * in the input that still needs to be removed.
 217        */
 218       if (is_special(nsb,sb,*pdba->resinfo[pdba->atom[i].resind].name,
 219                      *pdba->atomname[i]) &&
 220           !(nspec > 0 &&
 221             pdba->atom[sgp[nspec-1]].resind == pdba->atom[i].resind &&
 222             gmx_strcasecmp(*pdba->atomname[sgp[nspec-1]],
 223                        *pdba->atomname[i]) == 0)) {
 224         specp[nspec] = pdba->atom[i].resind;
 225         sgp[nspec] = i;
 226         nspec++;
 227       }
 228     }
 229     /* distance matrix d[nspec][nspec] */
 230     snew(d,nspec);
 231     for(i=0; (i<nspec); i++)
 232       snew(d[i],nspec);
 233
 234     for(i=0; (i<nspec); i++) {
 235       ai=sgp[i];
 236       for(j=0; (j<nspec); j++) {
 237         aj=sgp[j];
 238         d[i][j]=sqrt(distance2(x[ai],x[aj]));
 239       }
 240     }
 241     if (nspec > 1) {
 242 #define MAXCOL 7
 243       fprintf(stderr,"Special Atom Distance matrix:\n");
 244       for(b=0; (b<nspec); b+=MAXCOL) {
 245         /* print resname/number column headings */
 246         fprintf(stderr,"%8s%8s","","");
 247         e=min(b+MAXCOL, nspec-1);
 248         for(i=b; (i<e); i++) {
 249           sprintf(buf,"%s%d",*pdba->resinfo[pdba->atom[sgp[i]].resind].name,
 250                   pdba->resinfo[specp[i]].nr);
 251           fprintf(stderr,"%8s",buf);
 252         }
 253         fprintf(stderr,"\n");
 254         /* print atomname/number column headings */
 255         fprintf(stderr,"%8s%8s","","");
 256         e=min(b+MAXCOL, nspec-1);
 257         for(i=b; (i<e); i++) {
 258           sprintf(buf,"%s%d",*pdba->atomname[sgp[i]], sgp[i]+1);
 259           fprintf(stderr,"%8s",buf);
 260         }
 261         fprintf(stderr,"\n");
 262         /* print matrix */
 263         e=min(b+MAXCOL, nspec);
 264         for(i=b+1; (i<nspec); i++) {
 265           sprintf(buf,"%s%d",*pdba->resinfo[pdba->atom[sgp[i]].resind].name,
 266                   pdba->resinfo[specp[i]].nr);
 267           fprintf(stderr,"%8s",buf);
 268           sprintf(buf,"%s%d", *pdba->atomname[sgp[i]], sgp[i]+1);
 269           fprintf(stderr,"%8s",buf);
 270           e2=min(i,e);
 271           for(j=b; (j<e2); j++)
 272             fprintf(stderr," %7.3f",d[i][j]);
 273           fprintf(stderr,"\n");
 274         }
 275       }
 276     }
 277
 278     snew(bonds,nspec);
 279
 280     for(i=0; (i<nspec); i++) {
 281       ai = sgp[i];
 282       for(j=i+1; (j<nspec); j++) {
 283         aj = sgp[j];
 284         if (is_bond(nsb,sb,pdba,ai,aj,d[i][j],&index_sb,&bSwap)) {
 285           fprintf(stderr,"%s %s-%d %s-%d and %s-%d %s-%d%s",
 286                   bInteractive ? "Link" : "Linking",
 287                   *pdba->resinfo[pdba->atom[ai].resind].name,
 288                   pdba->resinfo[specp[i]].nr,
 289                   *pdba->atomname[ai], ai+1,
 290                   *pdba->resinfo[pdba->atom[aj].resind].name,
 291                   pdba->resinfo[specp[j]].nr,
 292                   *pdba->atomname[aj], aj+1,
 293                   bInteractive ? " (y/n) ?" : "...\n");
 294           bDoit=bInteractive ? yesno() : TRUE;
 295
 296           if (bDoit) {
 297             /* Store the residue numbers in the bonds array */
 298             bonds[nbonds].res1 = specp[i];
 299             bonds[nbonds].res2 = specp[j];
 300             bonds[nbonds].a1   = strdup(*pdba->atomname[ai]);
 301             bonds[nbonds].a2   = strdup(*pdba->atomname[aj]);
 302             /* rename residues */
 303             if (bSwap) {
 304               rename_1res(pdba,specp[i],sb[index_sb].newres2,bVerbose);
 305               rename_1res(pdba,specp[j],sb[index_sb].newres1,bVerbose);
 306             }
 307             else {
 308               rename_1res(pdba,specp[i],sb[index_sb].newres1,bVerbose);
 309               rename_1res(pdba,specp[j],sb[index_sb].newres2,bVerbose);
 310             }
 311             nbonds++;
 312           }
 313         }
 314       }
 315     }
 316
 317     for(i=0; (i<nspec); i++)
 318       sfree(d[i]);
 319     sfree(d);
 320     sfree(sgp);
 321     sfree(specp);
 322
 323     done_specbonds(nsb,sb);
 324     sfree(sb);
 325   }
 326
 327   *specbonds=bonds;
 328
 329   return nbonds;
 330 }
 331