src/gromacs/gmxpreprocess/resall.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,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,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 "resall.h"
  41
  42 #include <cctype>
  43 #include <cstdlib>
  44 #include <cstring>
  45
  46 #include <algorithm>
  47 #include <optional>
  48 #include <string>
  49 #include <vector>
  50
  51 #include "gromacs/gmxpreprocess/fflibutil.h"
  52 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
  53 #include "gromacs/gmxpreprocess/grompp_impl.h"
  54 #include "gromacs/gmxpreprocess/notset.h"
  55 #include "gromacs/gmxpreprocess/pgutil.h"
  56 #include "gromacs/topology/atoms.h"
  57 #include "gromacs/topology/symtab.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/logger.h"
  63 #include "gromacs/utility/smalloc.h"
  64 #include "gromacs/utility/strdb.h"
  65 #include "gromacs/utility/stringtoenumvalueconverter.h"
  66
  67 #include "hackblock.h"
  68
  69 PreprocessingAtomTypes read_atype(const char* ffdir, t_symtab* tab)
  70 {
  71     FILE*   in;
  72     char    buf[STRLEN], name[STRLEN];
  73     double  m;
  74     t_atom* a;
  75
  76     std::vector<std::string> files = fflib_search_file_end(ffdir, ".atp", TRUE);
  77     snew(a, 1);
  78     PreprocessingAtomTypes at;
  79
  80     for (const auto& filename : files)
  81     {
  82         in = fflib_open(filename);
  83         while (!feof(in))
  84         {
  85             /* Skip blank or comment-only lines */
  86             do
  87             {
  88                 if (fgets2(buf, STRLEN, in) != nullptr)
  89                 {
  90                     strip_comment(buf);
  91                     trim(buf);
  92                 }
  93             } while ((feof(in) == 0) && strlen(buf) == 0);
  94
  95             if (sscanf(buf, "%s%lf", name, &m) == 2)
  96             {
  97                 a->m = m;
  98                 at.addType(tab, *a, name, InteractionOfType({}, {}), 0, 0);
  99             }
 100             else
 101             {
 102                 gmx_fatal(FARGS, "Invalid atomtype format: '%s'", buf);
 103             }
 104         }
 105         gmx_ffclose(in);
 106     }
 107     sfree(a);
 108     return at;
 109 }
 110
 111 static void print_resatoms(FILE* out, const PreprocessingAtomTypes& atype, const PreprocessResidue& rtpDBEntry)
 112 {
 113     fprintf(out, "[ %s ]\n", rtpDBEntry.resname.c_str());
 114     fprintf(out, " [ atoms ]\n");
 115
 116     for (int j = 0; (j < rtpDBEntry.natom()); j++)
 117     {
 118         int  tp   = rtpDBEntry.atom[j].type;
 119         auto tpnm = atype.atomNameFromAtomType(tp);
 120         if (!tpnm.has_value())
 121         {
 122             gmx_fatal(FARGS, "Incorrect atomtype (%d)", tp);
 123         }
 124         fprintf(out,
 125                 "%6s  %6s  %8.3f  %6d\n",
 126                 *(rtpDBEntry.atomname[j]),
 127                 *tpnm,
 128                 rtpDBEntry.atom[j].q,
 129                 rtpDBEntry.cgnr[j]);
 130     }
 131 }
 132
 133 static bool read_atoms(FILE* in, char* line, PreprocessResidue* r0, t_symtab* tab, PreprocessingAtomTypes* atype)
 134 {
 135     int    cg;
 136     char   buf[256], buf1[256];
 137     double q;
 138
 139     /* Read Atoms */
 140     r0->atom.clear();
 141     r0->atomname.clear();
 142     r0->cgnr.clear();
 143     while (get_a_line(in, line, STRLEN) && (strchr(line, '[') == nullptr))
 144     {
 145         if (sscanf(line, "%s%s%lf%d", buf, buf1, &q, &cg) != 4)
 146         {
 147             return FALSE;
 148         }
 149         r0->atomname.push_back(put_symtab(tab, buf));
 150         r0->atom.emplace_back();
 151         r0->atom.back().q = q;
 152         r0->cgnr.push_back(cg);
 153         auto j = atype->atomTypeFromName(buf1);
 154         if (!j.has_value())
 155         {
 156             gmx_fatal(FARGS,
 157                       "Atom type %s (residue %s) not found in atomtype "
 158                       "database",
 159                       buf1,
 160                       r0->resname.c_str());
 161         }
 162         r0->atom.back().type = *j;
 163         r0->atom.back().m    = *atype->atomMassFromAtomType(*j);
 164     }
 165
 166     return TRUE;
 167 }
 168
 169 static bool read_bondeds(BondedTypes bt, FILE* in, char* line, PreprocessResidue* rtpDBEntry)
 170 {
 171     char str[STRLEN];
 172
 173     while (get_a_line(in, line, STRLEN) && (strchr(line, '[') == nullptr))
 174     {
 175         int n = 0;
 176         int ni;
 177         rtpDBEntry->rb[bt].b.emplace_back();
 178         BondedInteraction* newBond = &rtpDBEntry->rb[bt].b.back();
 179         for (int j = 0; j < enumValueToNumIAtoms(bt); j++)
 180         {
 181             if (sscanf(line + n, "%s%n", str, &ni) == 1)
 182             {
 183                 newBond->a[j] = str;
 184             }
 185             else
 186             {
 187                 return FALSE;
 188             }
 189             n += ni;
 190         }
 191         while (isspace(line[n]))
 192         {
 193             n++;
 194         }
 195         rtrim(line + n);
 196         newBond->s = line + n;
 197     }
 198
 199     return TRUE;
 200 }
 201
 202 static void print_resbondeds(FILE* out, BondedTypes bt, const PreprocessResidue& rtpDBEntry)
 203 {
 204     if (!rtpDBEntry.rb[bt].b.empty())
 205     {
 206         fprintf(out, " [ %s ]\n", enumValueToString(bt));
 207
 208         for (const auto& b : rtpDBEntry.rb[bt].b)
 209         {
 210             for (int j = 0; j < enumValueToNumIAtoms(bt); j++)
 211             {
 212                 fprintf(out, "%6s ", b.a[j].c_str());
 213             }
 214             if (!b.s.empty())
 215             {
 216                 fprintf(out, "    %s", b.s.c_str());
 217             }
 218             fprintf(out, "\n");
 219         }
 220     }
 221 }
 222
 223 static void check_rtp(gmx::ArrayRef<const PreprocessResidue> rtpDBEntry,
 224                       const std::string&                     libfn,
 225                       const gmx::MDLogger&                   logger)
 226 {
 227     /* check for double entries, assuming list is already sorted */
 228     for (auto it = rtpDBEntry.begin() + 1; it != rtpDBEntry.end(); it++)
 229     {
 230         auto prev = it - 1;
 231         if (gmx::equalCaseInsensitive(prev->resname, it->resname))
 232         {
 233             GMX_LOG(logger.warning)
 234                     .asParagraph()
 235                     .appendTextFormatted("Double entry %s in file %s", it->resname.c_str(), libfn.c_str());
 236         }
 237     }
 238 }
 239
 240 static std::optional<BondedTypes> get_bt(char* header)
 241 {
 242     gmx::StringToEnumValueConverter<BondedTypes, enumValueToString> converter;
 243     return converter.valueFrom(header);
 244 }
 245
 246 /* print all the BondedTypes type numbers */
 247 static void print_resall_header(FILE* out, gmx::ArrayRef<const PreprocessResidue> rtpDBEntry)
 248 {
 249     fprintf(out, "[ bondedtypes ]\n");
 250     fprintf(out,
 251             "; bonds  angles  dihedrals  impropers all_dihedrals nr_exclusions  HH14  "
 252             "remove_dih\n");
 253     fprintf(out,
 254             " %5d  %6d  %9d  %9d  %14d  %14d %14d %14d\n\n",
 255             rtpDBEntry[0].rb[0].type,
 256             rtpDBEntry[0].rb[1].type,
 257             rtpDBEntry[0].rb[2].type,
 258             rtpDBEntry[0].rb[3].type,
 259             static_cast<int>(rtpDBEntry[0].bKeepAllGeneratedDihedrals),
 260             rtpDBEntry[0].nrexcl,
 261             static_cast<int>(rtpDBEntry[0].bGenerateHH14Interactions),
 262             static_cast<int>(rtpDBEntry[0].bRemoveDihedralIfWithImproper));
 263 }
 264
 265
 266 static void print_resall_log(const gmx::MDLogger& logger, gmx::ArrayRef<const PreprocessResidue> rtpDBEntry)
 267 {
 268     GMX_LOG(logger.info).asParagraph().appendTextFormatted("[ bondedtypes ]");
 269     GMX_LOG(logger.info)
 270             .asParagraph()
 271             .appendTextFormatted(
 272                     "; bonds  angles  dihedrals  impropers all_dihedrals nr_exclusions  HH14  "
 273                     "remove_dih");
 274     GMX_LOG(logger.info)
 275             .asParagraph()
 276             .appendTextFormatted(" %5d  %6d  %9d  %9d  %14d  %14d %14d %14d",
 277                                  rtpDBEntry[0].rb[0].type,
 278                                  rtpDBEntry[0].rb[1].type,
 279                                  rtpDBEntry[0].rb[2].type,
 280                                  rtpDBEntry[0].rb[3].type,
 281                                  static_cast<int>(rtpDBEntry[0].bKeepAllGeneratedDihedrals),
 282                                  rtpDBEntry[0].nrexcl,
 283                                  static_cast<int>(rtpDBEntry[0].bGenerateHH14Interactions),
 284                                  static_cast<int>(rtpDBEntry[0].bRemoveDihedralIfWithImproper));
 285 }
 286
 287
 288 void print_resall(FILE* out, gmx::ArrayRef<const PreprocessResidue> rtpDBEntry, const PreprocessingAtomTypes& atype)
 289 {
 290     if (rtpDBEntry.empty())
 291     {
 292         return;
 293     }
 294
 295     print_resall_header(out, rtpDBEntry);
 296
 297     for (const auto& r : rtpDBEntry)
 298     {
 299         if (r.natom() > 0)
 300         {
 301             print_resatoms(out, atype, r);
 302             for (auto bt : gmx::EnumerationWrapper<BondedTypes>{})
 303             {
 304                 print_resbondeds(out, bt, r);
 305             }
 306         }
 307     }
 308 }
 309
 310 void readResidueDatabase(const std::string&              rrdb,
 311                          std::vector<PreprocessResidue>* rtpDBEntry,
 312                          PreprocessingAtomTypes*         atype,
 313                          t_symtab*                       tab,
 314                          const gmx::MDLogger&            logger,
 315                          bool                            bAllowOverrideRTP)
 316 {
 317     FILE* in;
 318     char  filebase[STRLEN], line[STRLEN], header[STRLEN];
 319     int   nparam;
 320     int   dum1, dum2, dum3;
 321     bool  bNextResidue, bError;
 322
 323     fflib_filename_base(rrdb.c_str(), filebase, STRLEN);
 324
 325     in = fflib_open(rrdb);
 326
 327     PreprocessResidue header_settings;
 328
 329     /* these bonded parameters will overwritten be when  *
 330      * there is a [ bondedtypes ] entry in the .rtp file */
 331     header_settings.rb[BondedTypes::Bonds].type             = 1; /* normal bonds     */
 332     header_settings.rb[BondedTypes::Angles].type            = 1; /* normal angles    */
 333     header_settings.rb[BondedTypes::ProperDihedrals].type   = 1; /* normal dihedrals */
 334     header_settings.rb[BondedTypes::ImproperDihedrals].type = 2; /* normal impropers */
 335     header_settings.rb[BondedTypes::Exclusions].type        = 1; /* normal exclusions */
 336     header_settings.rb[BondedTypes::Cmap].type              = 1; /* normal cmap torsions */
 337
 338     header_settings.bKeepAllGeneratedDihedrals    = FALSE;
 339     header_settings.nrexcl                        = 3;
 340     header_settings.bGenerateHH14Interactions     = TRUE;
 341     header_settings.bRemoveDihedralIfWithImproper = TRUE;
 342
 343     /* Column 5 & 6 aren't really bonded types, but we include
 344      * them here to avoid introducing a new section:
 345      * Column 5 : This controls the generation of dihedrals from the bonding.
 346      *            All possible dihedrals are generated automatically. A value of
 347      *            1 here means that all these are retained. A value of
 348      *            0 here requires generated dihedrals be removed if
 349      *              * there are any dihedrals on the same central atoms
 350      *                specified in the residue topology, or
 351      *              * there are other identical generated dihedrals
 352      *                sharing the same central atoms, or
 353      *              * there are other generated dihedrals sharing the
 354      *                same central bond that have fewer hydrogen atoms
 355      * Column 6: Number of bonded neighbors to exclude.
 356      * Column 7: Generate 1,4 interactions between two hydrogen atoms
 357      * Column 8: Remove proper dihedrals if centered on the same bond
 358      *           as an improper dihedral
 359      */
 360     get_a_line(in, line, STRLEN);
 361     if (!get_header(line, header))
 362     {
 363         gmx_fatal(FARGS, "in .rtp file at line:\n%s\n", line);
 364     }
 365     if (gmx::equalCaseInsensitive("bondedtypes", header, 5))
 366     {
 367         get_a_line(in, line, STRLEN);
 368         if ((nparam = sscanf(line,
 369                              "%d %d %d %d %d %d %d %d",
 370                              &header_settings.rb[BondedTypes::Bonds].type,
 371                              &header_settings.rb[BondedTypes::Angles].type,
 372                              &header_settings.rb[BondedTypes::ProperDihedrals].type,
 373                              &header_settings.rb[BondedTypes::ImproperDihedrals].type,
 374                              &dum1,
 375                              &header_settings.nrexcl,
 376                              &dum2,
 377                              &dum3))
 378             < 4)
 379         {
 380             gmx_fatal(FARGS,
 381                       "need 4 to 8 parameters in the header of .rtp file %s at line:\n%s\n",
 382                       rrdb.c_str(),
 383                       line);
 384         }
 385         header_settings.bKeepAllGeneratedDihedrals    = (dum1 != 0);
 386         header_settings.bGenerateHH14Interactions     = (dum2 != 0);
 387         header_settings.bRemoveDihedralIfWithImproper = (dum3 != 0);
 388         get_a_line(in, line, STRLEN);
 389         if (nparam < 5)
 390         {
 391             GMX_LOG(logger.info)
 392                     .asParagraph()
 393                     .appendTextFormatted("Using default: not generating all possible dihedrals");
 394             header_settings.bKeepAllGeneratedDihedrals = FALSE;
 395         }
 396         if (nparam < 6)
 397         {
 398             GMX_LOG(logger.info)
 399                     .asParagraph()
 400                     .appendTextFormatted("Using default: excluding 3 bonded neighbors");
 401             header_settings.nrexcl = 3;
 402         }
 403         if (nparam < 7)
 404         {
 405             GMX_LOG(logger.info)
 406                     .asParagraph()
 407                     .appendTextFormatted("Using default: generating 1,4 H--H interactions");
 408             header_settings.bGenerateHH14Interactions = TRUE;
 409         }
 410         if (nparam < 8)
 411         {
 412             GMX_LOG(logger.warning)
 413                     .asParagraph()
 414                     .appendTextFormatted(
 415                             "Using default: removing proper dihedrals found on the same bond as a "
 416                             "proper dihedral");
 417             header_settings.bRemoveDihedralIfWithImproper = TRUE;
 418         }
 419     }
 420     else
 421     {
 422         GMX_LOG(logger.warning)
 423                 .asParagraph()
 424                 .appendTextFormatted(
 425                         "Reading .rtp file without '[ bondedtypes ]' directive, "
 426                         "Will proceed as if the entry was:");
 427         print_resall_log(logger, gmx::arrayRefFromArray(&header_settings, 1));
 428     }
 429     /* We don't know the current size of rrtp, but simply realloc immediately */
 430     auto oldArrayEnd = rtpDBEntry->end();
 431     while (!feof(in))
 432     {
 433         /* Initialise rtp entry structure */
 434         rtpDBEntry->push_back(header_settings);
 435         PreprocessResidue* res = &rtpDBEntry->back();
 436         if (!get_header(line, header))
 437         {
 438             gmx_fatal(FARGS, "in .rtp file at line:\n%s\n", line);
 439         }
 440         res->resname  = header;
 441         res->filebase = filebase;
 442
 443         get_a_line(in, line, STRLEN);
 444         bError       = FALSE;
 445         bNextResidue = FALSE;
 446         do
 447         {
 448             if (!get_header(line, header))
 449             {
 450                 bError = TRUE;
 451             }
 452             else
 453             {
 454                 auto bt = get_bt(header);
 455                 if (bt.has_value())
 456                 {
 457                     /* header is an bonded directive */
 458                     bError = !read_bondeds(*bt, in, line, res);
 459                 }
 460                 else if (gmx::equalCaseInsensitive("atoms", header, 5))
 461                 {
 462                     /* header is the atoms directive */
 463                     bError = !read_atoms(in, line, res, tab, atype);
 464                 }
 465                 else
 466                 {
 467                     /* else header must be a residue name */
 468                     bNextResidue = TRUE;
 469                 }
 470             }
 471             if (bError)
 472             {
 473                 gmx_fatal(FARGS, "in .rtp file in residue %s at line:\n%s\n", res->resname.c_str(), line);
 474             }
 475         } while ((feof(in) == 0) && !bNextResidue);
 476
 477         if (res->natom() == 0)
 478         {
 479             gmx_fatal(FARGS, "No atoms found in .rtp file in residue %s\n", res->resname.c_str());
 480         }
 481
 482         auto found = std::find_if(
 483                 rtpDBEntry->begin(), rtpDBEntry->end() - 1, [&res](const PreprocessResidue& entry) {
 484                     return gmx::equalCaseInsensitive(entry.resname, res->resname);
 485                 });
 486
 487         if (found != rtpDBEntry->end() - 1)
 488         {
 489             if (found >= oldArrayEnd)
 490             {
 491                 gmx_fatal(FARGS, "Found a second entry for '%s' in '%s'", res->resname.c_str(), rrdb.c_str());
 492             }
 493             if (bAllowOverrideRTP)
 494             {
 495                 GMX_LOG(logger.warning)
 496                         .asParagraph()
 497                         .appendTextFormatted(
 498                                 "Found another rtp entry for '%s' in '%s',"
 499                                 " ignoring this entry and keeping the one from '%s.rtp'",
 500                                 res->resname.c_str(),
 501                                 rrdb.c_str(),
 502                                 found->filebase.c_str());
 503                 /* We should free all the data for this entry.
 504                  * The current code gives a lot of dangling pointers.
 505                  */
 506                 rtpDBEntry->erase(rtpDBEntry->end() - 1);
 507             }
 508             else
 509             {
 510                 gmx_fatal(FARGS,
 511                           "Found rtp entries for '%s' in both '%s' and '%s'. If you want the first "
 512                           "definition to override the second one, set the -rtpo option of pdb2gmx.",
 513                           res->resname.c_str(),
 514                           found->filebase.c_str(),
 515                           rrdb.c_str());
 516             }
 517         }
 518     }
 519     gmx_ffclose(in);
 520
 521     std::sort(rtpDBEntry->begin(), rtpDBEntry->end(), [](const PreprocessResidue& a, const PreprocessResidue& b) {
 522         return std::lexicographical_compare(
 523                 a.resname.begin(),
 524                 a.resname.end(),
 525                 b.resname.begin(),
 526                 b.resname.end(),
 527                 [](const char& c1, const char& c2) { return std::toupper(c1) < std::toupper(c2); });
 528     });
 529
 530     check_rtp(*rtpDBEntry, rrdb, logger);
 531 }
 532
 533 /************************************************************
 534  *
 535  *                  SEARCH   ROUTINES
 536  *
 537  ***********************************************************/
 538 static bool is_sign(char c)
 539 {
 540     return (c == '+' || c == '-');
 541 }
 542
 543 /* Compares if the strins match except for a sign at the end
 544  * of (only) one of the two.
 545  */
 546 static int neq_str_sign(const char* a1, const char* a2)
 547 {
 548     int l1, l2, lm;
 549
 550     l1 = static_cast<int>(strlen(a1));
 551     l2 = static_cast<int>(strlen(a2));
 552     lm = std::min(l1, l2);
 553
 554     if (lm >= 1 && ((l1 == l2 + 1 && is_sign(a1[l1 - 1])) || (l2 == l1 + 1 && is_sign(a2[l2 - 1])))
 555         && gmx::equalCaseInsensitive(a1, a2, lm))
 556     {
 557         return lm;
 558     }
 559     else
 560     {
 561         return 0;
 562     }
 563 }
 564
 565 std::string searchResidueDatabase(const std::string&                     key,
 566                                   gmx::ArrayRef<const PreprocessResidue> rtpDBEntry,
 567                                   const gmx::MDLogger&                   logger)
 568 {
 569     int         nbest, best, besti;
 570     std::string bestbuf;
 571
 572     nbest = 0;
 573     besti = -1;
 574     /* We want to match at least one character */
 575     best = 1;
 576     for (auto it = rtpDBEntry.begin(); it != rtpDBEntry.end(); it++)
 577     {
 578         if (gmx::equalCaseInsensitive(key, it->resname))
 579         {
 580             besti = std::distance(rtpDBEntry.begin(), it);
 581             nbest = 1;
 582             break;
 583         }
 584         else
 585         {
 586             /* Allow a mismatch of at most a sign character (with warning) */
 587             int n = neq_str_sign(key.c_str(), it->resname.c_str());
 588             if (n >= best && n + 1 >= gmx::index(key.length()) && n + 1 >= gmx::index(it->resname.length()))
 589             {
 590                 if (n == best)
 591                 {
 592                     if (nbest == 1)
 593                     {
 594                         bestbuf = rtpDBEntry[besti].resname;
 595                     }
 596                     if (nbest >= 1)
 597                     {
 598                         bestbuf.append(" or ");
 599                         bestbuf.append(it->resname);
 600                     }
 601                 }
 602                 else
 603                 {
 604                     nbest = 0;
 605                 }
 606                 besti = std::distance(rtpDBEntry.begin(), it);
 607                 best  = n;
 608                 nbest++;
 609             }
 610         }
 611     }
 612     if (nbest > 1)
 613     {
 614         gmx_fatal(FARGS,
 615                   "Residue '%s' not found in residue topology database, looks a bit like %s",
 616                   key.c_str(),
 617                   bestbuf.c_str());
 618     }
 619     else if (besti == -1)
 620     {
 621         gmx_fatal(FARGS, "Residue '%s' not found in residue topology database", key.c_str());
 622     }
 623     if (!gmx::equalCaseInsensitive(rtpDBEntry[besti].resname, key))
 624     {
 625         GMX_LOG(logger.warning)
 626                 .asParagraph()
 627                 .appendTextFormatted(
 628                         "'%s' not found in residue topology database, "
 629                         "trying to use '%s'",
 630                         key.c_str(),
 631                         rtpDBEntry[besti].resname.c_str());
 632     }
 633
 634     return rtpDBEntry[besti].resname;
 635 }
 636
 637 gmx::ArrayRef<const PreprocessResidue>::const_iterator
 638 getDatabaseEntry(const std::string& rtpname, gmx::ArrayRef<const PreprocessResidue> rtpDBEntry)
 639 {
 640     auto found = std::find_if(
 641             rtpDBEntry.begin(), rtpDBEntry.end(), [&rtpname](const PreprocessResidue& entry) {
 642                 return gmx::equalCaseInsensitive(rtpname, entry.resname);
 643             });
 644     if (found == rtpDBEntry.end())
 645     {
 646         /* This should never happen, since searchResidueDatabase should have been called
 647          * before calling getDatabaseEntry.
 648          */
 649         gmx_fatal(FARGS, "Residue type '%s' not found in residue topology database", rtpname.c_str());
 650     }
 651
 652     return found;
 653 }