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