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