src/kernel/ter_db.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, 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 "sysstuff.h"
  43 #include "smalloc.h"
  44 #include "typedefs.h"
  45 #include "symtab.h"
  46 #include "futil.h"
  47 #include "resall.h"
  48 #include "h_db.h"
  49 #include "string2.h"
  50 #include "strdb.h"
  51 #include "gmx_fatal.h"
  52 #include "ter_db.h"
  53 #include "toputil.h"
  54 #include "gmxfio.h"
  55 #include "fflibutil.h"
  56
  57
  58 /* use bonded types definitions in hackblock.h */
  59 #define ekwRepl ebtsNR+1
  60 #define ekwAdd  ebtsNR+2
  61 #define ekwDel  ebtsNR+3
  62 #define ekwNR   3
  63 const char *kw_names[ekwNR] = {
  64   "replace", "add", "delete"
  65 };
  66
  67 int find_kw(char *keyw)
  68 {
  69   int i;
  70
  71   for(i=0; i<ebtsNR; i++)
  72     if (gmx_strcasecmp(btsNames[i],keyw) == 0)
  73       return i;
  74   for(i=0; i<ekwNR; i++)
  75     if (gmx_strcasecmp(kw_names[i],keyw) == 0)
  76       return ebtsNR + 1 + i;
  77
  78   return NOTSET;
  79 }
  80
  81 #define FATAL() gmx_fatal(FARGS,"Reading Termini Database: not enough items on line\n%s",line)
  82
  83 static void read_atom(char *line, gmx_bool bAdd,
  84                       char **nname, t_atom *a, gpp_atomtype_t atype, int *cgnr)
  85 {
  86   int    nr, i;
  87   char   buf[5][30],type[12];
  88   double m, q;
  89
  90   /* This code is messy, because of support for different formats:
  91    * for replace: [new name] <atom type> <m> <q> [cgnr (old format)]
  92    * for add:                <atom type> <m> <q> [cgnr]
  93    */
  94   nr = sscanf(line,"%s %s %s %s %s",buf[0],buf[1],buf[2],buf[3],buf[4]);
  95
  96   /* Here there an ambiguity due to the old replace format with cgnr,
  97    * which was read for years, but ignored in the rest of the code.
  98    * We have to assume that the atom type does not start with a digit
  99    * to make a line with 4 entries uniquely interpretable.
 100    */
 101   if (!bAdd && nr == 4 && isdigit(buf[1][0])) {
 102     nr = 3;
 103   }
 104
 105   if (nr < 3 || nr > 4) {
 106     gmx_fatal(FARGS,"Reading Termini Database: expected %d or %d items of atom data in stead of %d on line\n%s", 3, 4, nr, line);
 107   }
 108   i = 0;
 109   if (!bAdd) {
 110     if (nr == 4) {
 111       *nname = strdup(buf[i++]);
 112     } else {
 113       *nname = NULL;
 114     }
 115   }
 116   a->type = get_atomtype_type(buf[i++],atype);
 117   sscanf(buf[i++],"%lf",&m);
 118   a->m = m;
 119   sscanf(buf[i++],"%lf",&q);
 120   a->q = q;
 121   if (bAdd && nr == 4) {
 122     sscanf(buf[i++],"%d", cgnr);
 123   } else {
 124     *cgnr = NOTSET;
 125   }
 126 }
 127
 128 static void print_atom(FILE *out,t_atom *a,gpp_atomtype_t atype,char *newnm)
 129 {
 130   fprintf(out,"\t%s\t%g\t%g\n",
 131           get_atomtype_name(a->type,atype),a->m,a->q);
 132 }
 133
 134 static void print_ter_db(const char *ff,char C,int nb,t_hackblock tb[],
 135                          gpp_atomtype_t atype)
 136 {
 137   FILE *out;
 138   int i,j,k,bt,nrepl,nadd,ndel;
 139   char buf[STRLEN],nname[STRLEN];
 140
 141   sprintf(buf,"%s-%c.tdb",ff,C);
 142   out = gmx_fio_fopen(buf,"w");
 143
 144   for(i=0; (i<nb); i++) {
 145     fprintf(out,"[ %s ]\n",tb[i].name);
 146
 147     /* first count: */
 148     nrepl=0;
 149     nadd=0;
 150     ndel=0;
 151     for(j=0; j<tb[i].nhack; j++)
 152       if ( tb[i].hack[j].oname!=NULL && tb[i].hack[j].nname!=NULL )
 153         nrepl++;
 154       else if ( tb[i].hack[j].oname==NULL && tb[i].hack[j].nname!=NULL )
 155         nadd++;
 156       else if ( tb[i].hack[j].oname!=NULL && tb[i].hack[j].nname==NULL )
 157         ndel++;
 158       else if ( tb[i].hack[j].oname==NULL && tb[i].hack[j].nname==NULL )
 159         gmx_fatal(FARGS,"invalid hack (%s) in termini database",tb[i].name);
 160     if (nrepl) {
 161       fprintf(out,"[ %s ]\n",kw_names[ekwRepl-ebtsNR-1]);
 162       for(j=0; j<tb[i].nhack; j++)
 163         if ( tb[i].hack[j].oname!=NULL && tb[i].hack[j].nname!=NULL ) {
 164           fprintf(out,"%s\t",tb[i].hack[j].oname);
 165           print_atom(out,tb[i].hack[j].atom,atype,tb[i].hack[j].nname);
 166         }
 167     }
 168     if (nadd) {
 169       fprintf(out,"[ %s ]\n",kw_names[ekwAdd-ebtsNR-1]);
 170       for(j=0; j<tb[i].nhack; j++)
 171         if ( tb[i].hack[j].oname==NULL && tb[i].hack[j].nname!=NULL ) {
 172           print_ab(out,&(tb[i].hack[j]),tb[i].hack[j].nname);
 173           print_atom(out,tb[i].hack[j].atom,atype,tb[i].hack[j].nname);
 174         }
 175     }
 176     if (ndel) {
 177       fprintf(out,"[ %s ]\n",kw_names[ekwDel-ebtsNR-1]);
 178       for(j=0; j<tb[i].nhack; j++)
 179         if ( tb[i].hack[j].oname!=NULL && tb[i].hack[j].nname==NULL )
 180           fprintf(out,"%s\n",tb[i].hack[j].oname);
 181     }
 182     for(bt=0; bt<ebtsNR; bt++)
 183       if (tb[i].rb[bt].nb) {
 184         fprintf(out,"[ %s ]\n", btsNames[bt]);
 185         for(j=0; j<tb[i].rb[bt].nb; j++) {
 186           for(k=0; k<btsNiatoms[bt]; k++)
 187             fprintf(out,"%s%s",k?"\t":"",tb[i].rb[bt].b[j].a[k]);
 188           if ( tb[i].rb[bt].b[j].s )
 189             fprintf(out,"\t%s",tb[i].rb[bt].b[j].s);
 190           fprintf(out,"\n");
 191         }
 192       }
 193     fprintf(out,"\n");
 194   }
 195   gmx_fio_fclose(out);
 196 }
 197
 198 static void read_ter_db_file(char *fn,
 199                              int *ntbptr,t_hackblock **tbptr,
 200                              gpp_atomtype_t atype)
 201 {
 202   char       filebase[STRLEN],*ptr;
 203   FILE       *in;
 204   char       header[STRLEN],buf[STRLEN],line[STRLEN];
 205   t_hackblock *tb;
 206   int        i,j,n,ni,kwnr,nb,maxnb,nh;
 207
 208   fflib_filename_base(fn,filebase,STRLEN);
 209   /* Remove the C/N termini extension */
 210   ptr = strrchr(filebase,'.');
 211   if (ptr != NULL) {
 212     ptr[0] = '\0';
 213   }
 214
 215   in = fflib_open(fn);
 216   if (debug)
 217     fprintf(debug,"Opened %s\n",fn);
 218
 219   tb = *tbptr;
 220   nb = *ntbptr - 1;
 221   maxnb=0;
 222   kwnr=NOTSET;
 223   get_a_line(in,line,STRLEN);
 224   while (!feof(in)) {
 225     if (get_header(line,header)) {
 226       /* this is a new block, or a new keyword */
 227       kwnr=find_kw(header);
 228
 229       if (kwnr == NOTSET) {
 230         nb++;
 231         /* here starts a new block */
 232         if ( nb >= maxnb ) {
 233           maxnb = nb + 100;
 234           srenew(tb,maxnb);
 235         }
 236         clear_t_hackblock(&tb[nb]);
 237         tb[nb].name     = strdup(header);
 238         tb[nb].filebase = strdup(filebase);
 239       }
 240     } else {
 241       if (nb < 0)
 242         gmx_fatal(FARGS,"reading termini database: "
 243                     "directive expected before line:\n%s\n"
 244                     "This might be a file in an old format.",line);
 245       /* this is not a header, so it must be data */
 246       if (kwnr >= ebtsNR) {
 247         /* this is a hack: add/rename/delete atoms */
 248         /* make space for hacks */
 249         if (tb[nb].nhack >= tb[nb].maxhack) {
 250           tb[nb].maxhack+=10;
 251           srenew(tb[nb].hack, tb[nb].maxhack);
 252         }
 253         nh=tb[nb].nhack;
 254         clear_t_hack(&(tb[nb].hack[nh]));
 255         for(i=0; i<4; i++)
 256           tb[nb].hack[nh].a[i]=NULL;
 257         tb[nb].nhack++;
 258
 259         /* get data */
 260         n=0;
 261         if ( kwnr==ekwRepl || kwnr==ekwDel ) {
 262           if (sscanf(line, "%s%n", buf, &n) != 1)
 263             gmx_fatal(FARGS,"Reading Termini Database '%s': "
 264                       "expected atom name on line\n%s",fn,line);
 265           tb[nb].hack[nh].oname = strdup(buf);
 266           /* we only replace or delete one atom at a time */
 267           tb[nb].hack[nh].nr = 1;
 268         } else if ( kwnr==ekwAdd ) {
 269           read_ab(line, fn, &(tb[nb].hack[nh]));
 270           get_a_line(in, line, STRLEN);
 271         } else
 272           gmx_fatal(FARGS,"unimplemented keyword number %d (%s:%d)",
 273                       kwnr,__FILE__,__LINE__);
 274         if ( kwnr==ekwRepl || kwnr==ekwAdd ) {
 275           snew(tb[nb].hack[nh].atom, 1);
 276           read_atom(line+n,kwnr==ekwAdd,
 277                     &tb[nb].hack[nh].nname, tb[nb].hack[nh].atom, atype,
 278                     &tb[nb].hack[nh].cgnr);
 279           if (tb[nb].hack[nh].nname == NULL) {
 280             if (tb[nb].hack[nh].oname != NULL) {
 281               tb[nb].hack[nh].nname = strdup(tb[nb].hack[nh].oname);
 282             } else {
 283               gmx_fatal(FARGS,"Reading Termini Database '%s': don't know which name the new atom should have on line\n%s",fn,line);
 284             }
 285           }
 286         }
 287       } else if (kwnr >= 0 && kwnr < ebtsNR) {
 288         /* this is bonded data: bonds, angles, dihedrals or impropers */
 289         srenew(tb[nb].rb[kwnr].b,tb[nb].rb[kwnr].nb+1);
 290         n=0;
 291         for(j=0; j<btsNiatoms[kwnr]; j++) {
 292           if ( sscanf(line+n, "%s%n", buf, &ni) == 1 )
 293             tb[nb].rb[kwnr].b[tb[nb].rb[kwnr].nb].a[j] = strdup(buf);
 294           else
 295             gmx_fatal(FARGS,"Reading Termini Database '%s': expected %d atom names (found %d) on line\n%s", fn, btsNiatoms[kwnr], j-1, line);
 296           n+=ni;
 297         }
 298         for(   ; j<MAXATOMLIST; j++)
 299           tb[nb].rb[kwnr].b[tb[nb].rb[kwnr].nb].a[j] = NULL;
 300         strcpy(buf, "");
 301         sscanf(line+n, "%s", buf);
 302         tb[nb].rb[kwnr].b[tb[nb].rb[kwnr].nb].s = strdup(buf);
 303         tb[nb].rb[kwnr].nb++;
 304       } else
 305         gmx_fatal(FARGS,"Reading Termini Database: Expecting a header at line\n"
 306                     "%s",line);
 307     }
 308     get_a_line(in,line,STRLEN);
 309   }
 310   nb++;
 311   srenew(tb,nb);
 312
 313   ffclose(in);
 314
 315   *ntbptr = nb;
 316   *tbptr  = tb;
 317 }
 318
 319 int read_ter_db(const char *ffdir,char ter,
 320                 t_hackblock **tbptr,gpp_atomtype_t atype)
 321 {
 322   char ext[STRLEN];
 323   int  ntdbf,f;
 324   char **tdbf;
 325   int  ntb;
 326
 327   sprintf(ext,".%c.tdb",ter);
 328
 329   /* Search for termini database files.
 330    * Do not generate an error when none are found.
 331    */
 332   ntdbf = fflib_search_file_end(ffdir,ext,FALSE,&tdbf);
 333   ntb    = 0;
 334   *tbptr = NULL;
 335   for(f=0; f<ntdbf; f++) {
 336     read_ter_db_file(tdbf[f],&ntb,tbptr,atype);
 337     sfree(tdbf[f]);
 338   }
 339   sfree(tdbf);
 340
 341   if (debug) {
 342     print_ter_db("new",ter,ntb,*tbptr,atype);
 343   }
 344
 345   return ntb;
 346 }
 347
 348 t_hackblock **filter_ter(int nrtp,t_restp rtp[],
 349                          int nb,t_hackblock tb[],
 350                          const char *resname,
 351                          const char *rtpname,
 352                          int *nret)
 353 {
 354     /* Since some force fields (e.g. OPLS) needs different
 355      * atomtypes for different residues there could be a lot
 356      * of entries in the databases for specific residues
 357      * (e.g. GLY-NH3+, SER-NH3+, ALA-NH3+).
 358      *
 359      * To reduce the database size, we assume that a terminus specifier liker
 360      *
 361      * [ GLY|SER|ALA-NH3+ ]
 362      *
 363      * would cover all of the three residue types above.
 364      * Wildcards (*,?) are not OK. Don't worry about e.g. GLU vs. GLUH since
 365      * pdb2gmx only uses the first 3 letters when calling this routine.
 366      *
 367      * To automate this, this routines scans a list of termini
 368      * for the residue name "resname" and returns an allocated list of
 369      * pointers to the termini that could be applied to the
 370      * residue in question. The variable pointed to by nret will
 371      * contain the number of valid pointers in the list.
 372      * Remember to free the list when you are done with it...
 373      */
 374
 375     t_restp *   restp;
 376     int         i,j,n,len,none_idx;
 377     gmx_bool        found;
 378     char        *rtpname_match,*s,*s2,*c;
 379     t_hackblock **list;
 380
 381     rtpname_match = search_rtp(rtpname,nrtp,rtp);
 382     restp = get_restp(rtpname_match,nrtp,rtp);
 383
 384     n=0;
 385     list=NULL;
 386
 387     for(i=0;i<nb;i++)
 388     {
 389         s=tb[i].name;
 390         found=FALSE;
 391         do
 392         {
 393             /* The residue name should appear in a tdb file with the same base name
 394              * as the file containing the rtp entry.
 395              * This makes termini selection for different molecule types
 396              * much cleaner.
 397              */
 398             if (gmx_strcasecmp(restp->filebase,tb[i].filebase) == 0 &&
 399                 gmx_strncasecmp(resname,s,3) == 0)
 400             {
 401                 found=TRUE;
 402                 srenew(list,n+1);
 403                 list[n]=&(tb[i]);
 404                 n++;
 405             }
 406             else
 407             {
 408                 /* advance to next |-separated field */
 409                 s=strchr(s,'|');
 410                 if(s!=NULL)
 411                 {
 412                     s++;
 413                 }
 414             }
 415         }
 416         while(!found && s!=NULL);
 417     }
 418
 419     /* All residue-specific termini have been added. See if there
 420      * are some generic ones by searching for the occurence of
 421      * '-' in the name prior to the last position (which indicates charge).
 422      * The [ None ] alternative is special since we don't want that
 423      * to be the default, so we put it last in the list we return.
 424      */
 425     none_idx=-1;
 426     for(i=0;i<nb;i++)
 427     {
 428         s=tb[i].name;
 429         /* The residue name should appear in a tdb file with the same base name
 430          * as the file containing the rtp entry.
 431          * This makes termini selection for different molecule types
 432          * much cleaner.
 433          */
 434         if(gmx_strcasecmp(restp->filebase,tb[i].filebase) == 0)
 435         {
 436             if(!gmx_strcasecmp("None",s))
 437             {
 438                 none_idx=i;
 439             }
 440             else
 441             {
 442                 c=strchr(s,'-');
 443                 if(c==NULL || ((c-s+1)==strlen(s)))
 444                 {
 445                     /* Check that we haven't already added a residue-specific version
 446                      * of this terminus.
 447                      */
 448                     for(j=0;j<n && strstr((*list[j]).name,s)==NULL;j++);
 449                     if(j==n)
 450                     {
 451                         srenew(list,n+1);
 452                         list[n]=&(tb[i]);
 453                         n++;
 454                     }
 455                 }
 456             }
 457         }
 458     }
 459     if(none_idx>=0)
 460     {
 461         srenew(list,n+1);
 462         list[n]=&(tb[none_idx]);
 463         n++;
 464     }
 465
 466     *nret=n;
 467     return list;
 468 }
 469
 470
 471 t_hackblock *choose_ter(int nb,t_hackblock **tb,const char *title)
 472 {
 473   int i,sel,ret;
 474
 475   printf("%s\n",title);
 476   for(i=0; (i<nb); i++)
 477   {
 478     char *advice_string = "";
 479     if (0 == gmx_wcmatch("*ZWITTERION*", (*tb[i]).name))
 480     {
 481       advice_string = " (only use with zwitterions containing exactly one residue)";
 482     }
 483     printf("%2d: %s%s\n",i,(*tb[i]).name, advice_string);
 484   }
 485   do {
 486     ret=fscanf(stdin,"%d",&sel);
 487   } while ((ret != 1) || (sel < 0) || (sel >= nb));
 488
 489   return tb[sel];
 490 }
 491