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.
   7  * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "ter_db.h"
  41
  42 #include <cctype>
  43 #include <cstring>
  44
  45 #include <algorithm>
  46 #include <string>
  47 #include <vector>
  48
  49 #include "gromacs/fileio/gmxfio.h"
  50 #include "gromacs/gmxpreprocess/fflibutil.h"
  51 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
  52 #include "gromacs/gmxpreprocess/grompp_impl.h"
  53 #include "gromacs/gmxpreprocess/h_db.h"
  54 #include "gromacs/gmxpreprocess/notset.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 #include "resall.h"
  64
  65 /* use bonded types definitions in hackblock.h */
  66 #define ekwRepl (ebtsNR + 1)
  67 #define ekwAdd (ebtsNR + 2)
  68 #define ekwDel (ebtsNR + 3)
  69 #define ekwNR 3
  70 static const char* kw_names[ekwNR] = { "replace", "add", "delete" };
  71
  72 static int find_kw(char* keyw)
  73 {
  74     int i;
  75
  76     for (i = 0; i < ebtsNR; i++)
  77     {
  78         if (gmx_strcasecmp(btsNames[i], keyw) == 0)
  79         {
  80             return i;
  81         }
  82     }
  83     for (i = 0; i < ekwNR; i++)
  84     {
  85         if (gmx_strcasecmp(kw_names[i], keyw) == 0)
  86         {
  87             return ebtsNR + 1 + i;
  88         }
  89     }
  90
  91     return NOTSET;
  92 }
  93
  94 #define FATAL() gmx_fatal(FARGS, "Reading Termini Database: not enough items on line\n%s", line)
  95
  96 static void read_atom(char* line, bool bAdd, std::string* nname, t_atom* a, PreprocessingAtomTypes* 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,
 121                   "Reading Termini Database: expected %d or %d items of atom data in stead of %d "
 122                   "on line\n%s",
 123                   3,
 124                   4,
 125                   nr,
 126                   line);
 127     }
 128     i = 0;
 129     if (!bAdd)
 130     {
 131         if (nr == 4)
 132         {
 133             *nname = buf[i++];
 134         }
 135         else
 136         {
 137             *nname = "";
 138         }
 139     }
 140     a->type = *atype->atomTypeFromName(buf[i++]);
 141     sscanf(buf[i++], "%lf", &m);
 142     a->m = m;
 143     sscanf(buf[i++], "%lf", &q);
 144     a->q = q;
 145     if (bAdd && nr == 4)
 146     {
 147         sscanf(buf[i++], "%d", cgnr);
 148     }
 149     else
 150     {
 151         *cgnr = NOTSET;
 152     }
 153 }
 154
 155 static void print_atom(FILE* out, const t_atom& a, PreprocessingAtomTypes* atype)
 156 {
 157     fprintf(out, "\t%s\t%g\t%g\n", *atype->atomNameFromAtomType(a.type), a.m, a.q);
 158 }
 159
 160 static void print_ter_db(const char*                                ff,
 161                          char                                       C,
 162                          gmx::ArrayRef<const MoleculePatchDatabase> tb,
 163                          PreprocessingAtomTypes*                    atype)
 164 {
 165     std::string buf = gmx::formatString("%s-%c.tdb", ff, C);
 166     FILE*       out = gmx_fio_fopen(buf.c_str(), "w");
 167
 168     for (const auto& modification : tb)
 169     {
 170         fprintf(out, "[ %s ]\n", modification.name.c_str());
 171
 172         if (std::any_of(modification.hack.begin(), modification.hack.end(), [](const auto& mod) {
 173                 return mod.type() == MoleculePatchType::Replace;
 174             }))
 175         {
 176             fprintf(out, "[ %s ]\n", kw_names[ekwRepl - ebtsNR - 1]);
 177             for (const auto& hack : modification.hack)
 178             {
 179                 if (hack.type() == MoleculePatchType::Replace)
 180                 {
 181                     fprintf(out, "%s\t", hack.oname.c_str());
 182                     print_atom(out, hack.atom.back(), atype);
 183                 }
 184             }
 185         }
 186         if (std::any_of(modification.hack.begin(), modification.hack.end(), [](const auto& mod) {
 187                 return mod.type() == MoleculePatchType::Add;
 188             }))
 189         {
 190             fprintf(out, "[ %s ]\n", kw_names[ekwAdd - ebtsNR - 1]);
 191             for (const auto& hack : modification.hack)
 192             {
 193                 if (hack.type() == MoleculePatchType::Add)
 194                 {
 195                     print_ab(out, hack, hack.nname.c_str());
 196                     print_atom(out, hack.atom.back(), atype);
 197                 }
 198             }
 199         }
 200         if (std::any_of(modification.hack.begin(), modification.hack.end(), [](const auto& mod) {
 201                 return mod.type() == MoleculePatchType::Delete;
 202             }))
 203         {
 204             fprintf(out, "[ %s ]\n", kw_names[ekwDel - ebtsNR - 1]);
 205             for (const auto& hack : modification.hack)
 206             {
 207                 if (hack.type() == MoleculePatchType::Delete)
 208                 {
 209                     fprintf(out, "%s\n", hack.oname.c_str());
 210                 }
 211             }
 212         }
 213         for (int bt = 0; bt < ebtsNR; bt++)
 214         {
 215             if (!modification.rb[bt].b.empty())
 216             {
 217                 fprintf(out, "[ %s ]\n", btsNames[bt]);
 218                 for (const auto& b : modification.rb[bt].b)
 219                 {
 220                     for (int k = 0; k < btsNiatoms[bt]; k++)
 221                     {
 222                         fprintf(out, "%s%s", k ? "\t" : "", b.a[k].c_str());
 223                     }
 224                     if (!b.s.empty())
 225                     {
 226                         fprintf(out, "\t%s", b.s.c_str());
 227                     }
 228                     fprintf(out, "\n");
 229                 }
 230             }
 231         }
 232         fprintf(out, "\n");
 233     }
 234     gmx_fio_fclose(out);
 235 }
 236
 237 static void read_ter_db_file(const char*                         fn,
 238                              std::vector<MoleculePatchDatabase>* tbptr,
 239                              PreprocessingAtomTypes*             atype)
 240 {
 241     char filebase[STRLEN];
 242     char header[STRLEN], buf[STRLEN], line[STRLEN];
 243
 244     fflib_filename_base(fn, filebase, STRLEN);
 245     /* Remove the C/N termini extension */
 246     char* ptr = strrchr(filebase, '.');
 247     if (ptr != nullptr)
 248     {
 249         ptr[0] = '\0';
 250     }
 251
 252     FILE* in = fflib_open(fn);
 253
 254     int kwnr = NOTSET;
 255     get_a_line(in, line, STRLEN);
 256     MoleculePatchDatabase* block = nullptr;
 257     while (!feof(in))
 258     {
 259         if (get_header(line, header))
 260         {
 261             /* this is a new block, or a new keyword */
 262             kwnr = find_kw(header);
 263
 264             if (kwnr == NOTSET)
 265             {
 266                 tbptr->emplace_back(MoleculePatchDatabase());
 267                 block = &tbptr->back();
 268                 clearModificationBlock(block);
 269                 block->name     = header;
 270                 block->filebase = filebase;
 271             }
 272         }
 273         else
 274         {
 275             if (block == nullptr)
 276             {
 277                 gmx_fatal(FARGS,
 278                           "reading termini database: "
 279                           "directive expected before line:\n%s\n"
 280                           "This might be a file in an old format.",
 281                           line);
 282             }
 283             /* this is not a header, so it must be data */
 284             if (kwnr >= ebtsNR)
 285             {
 286                 /* this is a hack: add/rename/delete atoms */
 287                 /* make space for hacks */
 288                 block->hack.emplace_back(MoleculePatch());
 289                 MoleculePatch* hack = &block->hack.back();
 290
 291                 /* get data */
 292                 int n = 0;
 293                 if (kwnr == ekwRepl || kwnr == ekwDel)
 294                 {
 295                     if (sscanf(line, "%s%n", buf, &n) != 1)
 296                     {
 297                         gmx_fatal(FARGS,
 298                                   "Reading Termini Database '%s': "
 299                                   "expected atom name on line\n%s",
 300                                   fn,
 301                                   line);
 302                     }
 303                     hack->oname = buf;
 304                     /* we only replace or delete one atom at a time */
 305                     hack->nr = 1;
 306                 }
 307                 else if (kwnr == ekwAdd)
 308                 {
 309                     read_ab(line, fn, hack);
 310                     get_a_line(in, line, STRLEN);
 311                 }
 312                 else
 313                 {
 314                     gmx_fatal(FARGS, "unimplemented keyword number %d (%s:%d)", kwnr, __FILE__, __LINE__);
 315                 }
 316                 if (kwnr == ekwRepl || kwnr == ekwAdd)
 317                 {
 318                     hack->atom.emplace_back();
 319                     read_atom(line + n, kwnr == ekwAdd, &hack->nname, &hack->atom.back(), atype, &hack->cgnr);
 320                     if (hack->nname.empty())
 321                     {
 322                         if (!hack->oname.empty())
 323                         {
 324                             hack->nname = hack->oname;
 325                         }
 326                         else
 327                         {
 328                             gmx_fatal(FARGS,
 329                                       "Reading Termini Database '%s': don't know which name the "
 330                                       "new atom should have on line\n%s",
 331                                       fn,
 332                                       line);
 333                         }
 334                     }
 335                 }
 336             }
 337             else if (kwnr >= 0 && kwnr < ebtsNR)
 338             {
 339                 /* this is bonded data: bonds, angles, dihedrals or impropers */
 340                 int n = 0;
 341                 block->rb[kwnr].b.emplace_back();
 342                 BondedInteraction* newBond = &block->rb[kwnr].b.back();
 343                 for (int j = 0; j < btsNiatoms[kwnr]; j++)
 344                 {
 345                     int ni;
 346                     if (sscanf(line + n, "%s%n", buf, &ni) == 1)
 347                     {
 348                         newBond->a[j] = buf;
 349                     }
 350                     else
 351                     {
 352                         gmx_fatal(FARGS,
 353                                   "Reading Termini Database '%s': expected %d atom names (found "
 354                                   "%d) on line\n%s",
 355                                   fn,
 356                                   btsNiatoms[kwnr],
 357                                   j - 1,
 358                                   line);
 359                     }
 360                     n += ni;
 361                 }
 362                 strcpy(buf, "");
 363                 sscanf(line + n, "%s", buf);
 364                 newBond->s = buf;
 365             }
 366             else
 367             {
 368                 gmx_fatal(FARGS,
 369                           "Reading Termini Database: Expecting a header at line\n"
 370                           "%s",
 371                           line);
 372             }
 373         }
 374         get_a_line(in, line, STRLEN);
 375     }
 376
 377     gmx_ffclose(in);
 378 }
 379
 380 int read_ter_db(const char* ffdir, char ter, std::vector<MoleculePatchDatabase>* tbptr, PreprocessingAtomTypes* atype)
 381 {
 382     std::string ext = gmx::formatString(".%c.tdb", ter);
 383
 384     /* Search for termini database files.
 385      * Do not generate an error when none are found.
 386      */
 387     std::vector<std::string> tdbf = fflib_search_file_end(ffdir, ext.c_str(), FALSE);
 388     tbptr->clear();
 389     for (const auto& filename : tdbf)
 390     {
 391         read_ter_db_file(filename.c_str(), tbptr, atype);
 392     }
 393
 394     if (debug)
 395     {
 396         print_ter_db("new", ter, *tbptr, atype);
 397     }
 398
 399     return tbptr->size();
 400 }
 401
 402 std::vector<MoleculePatchDatabase*> filter_ter(gmx::ArrayRef<MoleculePatchDatabase> tb, const char* resname)
 403 {
 404     // TODO Four years later, no force fields have ever used this, so decide status of this feature
 405     /* Since some force fields (e.g. OPLS) needs different
 406      * atomtypes for different residues there could be a lot
 407      * of entries in the databases for specific residues
 408      * (e.g. GLY-NH3+, SER-NH3+, ALA-NH3+).
 409      *
 410      * To reduce the database size, we assume that a terminus specifier liker
 411      *
 412      * [ GLY|SER|ALA-NH3+ ]
 413      *
 414      * would cover all of the three residue types above.
 415      * Wildcards (*,?) are not OK. Don't worry about e.g. GLU vs. GLUH since
 416      * pdb2gmx only uses the first 3 letters when calling this routine.
 417      *
 418      * To automate this, this routines scans a list of termini
 419      * for the residue name "resname" and returns an allocated list of
 420      * pointers to the termini that could be applied to the
 421      * residue in question. The variable pointed to by nret will
 422      * contain the number of valid pointers in the list.
 423      * Remember to free the list when you are done with it...
 424      */
 425
 426     auto                                none_idx = tb.end();
 427     std::vector<MoleculePatchDatabase*> list;
 428
 429     for (auto it = tb.begin(); it != tb.end(); it++)
 430     {
 431         const char* s     = it->name.c_str();
 432         bool        found = false;
 433         do
 434         {
 435             if (gmx::equalCaseInsensitive(resname, s, 3))
 436             {
 437                 found = true;
 438                 list.push_back(&*it);
 439             }
 440             else
 441             {
 442                 /* advance to next |-separated field */
 443                 s = strchr(s, '|');
 444                 if (s != nullptr)
 445                 {
 446                     s++;
 447                 }
 448             }
 449         } while (!found && s != nullptr);
 450     }
 451
 452     /* All residue-specific termini have been added. We might have to fall
 453      * back on generic termini, which are characterized by not having
 454      * '-' in the name prior to the last position (which indicates charge).
 455      * The [ None ] alternative is special since we don't want that
 456      * to be the default, so we put it last in the list we return.
 457      */
 458     for (auto it = tb.begin(); it != tb.end(); it++)
 459     {
 460         const char* s = it->name.c_str();
 461         if (gmx::equalCaseInsensitive("None", it->name))
 462         {
 463             none_idx = it;
 464         }
 465         else
 466         {
 467             /* Time to see if there's a generic terminus that matches.
 468                Is there a hyphen? */
 469             const char* c = strchr(s, '-');
 470
 471             /* A conjunction hyphen normally indicates a residue-specific
 472                terminus, which is named like "GLY-COOH". A generic terminus
 473                won't have a hyphen. */
 474             bool bFoundAnyHyphen = (c != nullptr);
 475             /* '-' as the last character indicates charge, so if that's
 476                the only one found e.g. "COO-", then it was not a conjunction
 477                hyphen, so this is a generic terminus */
 478             bool bOnlyFoundChargeHyphen = (bFoundAnyHyphen && *(c + 1) == '\0');
 479             /* Thus, "GLY-COO-" is not recognized as a generic terminus. */
 480             bool bFoundGenericTerminus = !bFoundAnyHyphen || bOnlyFoundChargeHyphen;
 481             if (bFoundGenericTerminus)
 482             {
 483                 /* Check that we haven't already added a residue-specific version
 484                  * of this terminus.
 485                  */
 486                 auto found = std::find_if(list.begin(), list.end(), [&s](const MoleculePatchDatabase* b) {
 487                     return strstr(b->name.c_str(), s) != nullptr;
 488                 });
 489                 if (found == list.end())
 490                 {
 491                     list.push_back(&*it);
 492                 }
 493             }
 494         }
 495     }
 496     if (none_idx != tb.end())
 497     {
 498         list.push_back(&*none_idx);
 499     }
 500
 501     return list;
 502 }
 503
 504
 505 MoleculePatchDatabase* choose_ter(gmx::ArrayRef<MoleculePatchDatabase*> tb, const char* title)
 506 {
 507     int sel, ret;
 508
 509     printf("%s\n", title);
 510     int i = 0;
 511     for (const auto& modification : tb)
 512     {
 513         bool bIsZwitterion = (0 == gmx_wcmatch("*ZWITTERION*", modification->name.c_str()));
 514         printf("%2d: %s%s\n",
 515                i,
 516                modification->name.c_str(),
 517                bIsZwitterion ? " (only use with zwitterions containing exactly one residue)" : "");
 518         i++;
 519     }
 520     do
 521     {
 522         ret = fscanf(stdin, "%d", &sel);
 523     } while ((ret != 1) || (sel < 0) || (sel >= tb.ssize()));
 524
 525     return tb[sel];
 526 }