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