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