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