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