src/gromacs/gmxpreprocess/ter_db.cpp

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