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