2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
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"
67 #include "hackblock.h"
69 PreprocessingAtomTypes read_atype(const char* ffdir, t_symtab* tab)
72 char buf[STRLEN], name[STRLEN];
76 std::vector<std::string> files = fflib_search_file_end(ffdir, ".atp", TRUE);
78 PreprocessingAtomTypes at;
80 for (const auto& filename : files)
82 in = fflib_open(filename);
85 /* Skip blank or comment-only lines */
88 if (fgets2(buf, STRLEN, in) != nullptr)
93 } while ((feof(in) == 0) && strlen(buf) == 0);
95 if (sscanf(buf, "%s%lf", name, &m) == 2)
98 at.addType(tab, *a, name, InteractionOfType({}, {}), 0, 0);
102 gmx_fatal(FARGS, "Invalid atomtype format: '%s'", buf);
111 static void print_resatoms(FILE* out, const PreprocessingAtomTypes& atype, const PreprocessResidue& rtpDBEntry)
113 fprintf(out, "[ %s ]\n", rtpDBEntry.resname.c_str());
114 fprintf(out, " [ atoms ]\n");
116 for (int j = 0; (j < rtpDBEntry.natom()); j++)
118 int tp = rtpDBEntry.atom[j].type;
119 auto tpnm = atype.atomNameFromAtomType(tp);
120 if (!tpnm.has_value())
122 gmx_fatal(FARGS, "Incorrect atomtype (%d)", tp);
125 "%6s %6s %8.3f %6d\n",
126 *(rtpDBEntry.atomname[j]),
128 rtpDBEntry.atom[j].q,
133 static bool read_atoms(FILE* in, char* line, PreprocessResidue* r0, t_symtab* tab, PreprocessingAtomTypes* atype)
136 char buf[256], buf1[256];
141 r0->atomname.clear();
143 while (get_a_line(in, line, STRLEN) && (strchr(line, '[') == nullptr))
145 if (sscanf(line, "%s%s%lf%d", buf, buf1, &q, &cg) != 4)
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);
157 "Atom type %s (residue %s) not found in atomtype "
160 r0->resname.c_str());
162 r0->atom.back().type = *j;
163 r0->atom.back().m = *atype->atomMassFromAtomType(*j);
169 static bool read_bondeds(BondedTypes bt, FILE* in, char* line, PreprocessResidue* rtpDBEntry)
173 while (get_a_line(in, line, STRLEN) && (strchr(line, '[') == nullptr))
177 rtpDBEntry->rb[bt].b.emplace_back();
178 BondedInteraction* newBond = &rtpDBEntry->rb[bt].b.back();
179 for (int j = 0; j < enumValueToNumIAtoms(bt); j++)
181 if (sscanf(line + n, "%s%n", str, &ni) == 1)
191 while (isspace(line[n]))
196 newBond->s = line + n;
202 static void print_resbondeds(FILE* out, BondedTypes bt, const PreprocessResidue& rtpDBEntry)
204 if (!rtpDBEntry.rb[bt].b.empty())
206 fprintf(out, " [ %s ]\n", enumValueToString(bt));
208 for (const auto& b : rtpDBEntry.rb[bt].b)
210 for (int j = 0; j < enumValueToNumIAtoms(bt); j++)
212 fprintf(out, "%6s ", b.a[j].c_str());
216 fprintf(out, " %s", b.s.c_str());
223 static void check_rtp(gmx::ArrayRef<const PreprocessResidue> rtpDBEntry,
224 const std::string& libfn,
225 const gmx::MDLogger& logger)
227 /* check for double entries, assuming list is already sorted */
228 for (auto it = rtpDBEntry.begin() + 1; it != rtpDBEntry.end(); it++)
231 if (gmx::equalCaseInsensitive(prev->resname, it->resname))
233 GMX_LOG(logger.warning)
235 .appendTextFormatted("Double entry %s in file %s", it->resname.c_str(), libfn.c_str());
240 static std::optional<BondedTypes> get_bt(char* header)
242 gmx::StringToEnumValueConverter<BondedTypes, enumValueToString> converter;
243 return converter.valueFrom(header);
246 /* print all the BondedTypes type numbers */
247 static void print_resall_header(FILE* out, gmx::ArrayRef<const PreprocessResidue> rtpDBEntry)
249 fprintf(out, "[ bondedtypes ]\n");
251 "; bonds angles dihedrals impropers all_dihedrals nr_exclusions HH14 "
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));
266 static void print_resall_log(const gmx::MDLogger& logger, gmx::ArrayRef<const PreprocessResidue> rtpDBEntry)
268 GMX_LOG(logger.info).asParagraph().appendTextFormatted("[ bondedtypes ]");
271 .appendTextFormatted(
272 "; bonds angles dihedrals impropers all_dihedrals nr_exclusions HH14 "
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));
288 void print_resall(FILE* out, gmx::ArrayRef<const PreprocessResidue> rtpDBEntry, const PreprocessingAtomTypes& atype)
290 if (rtpDBEntry.empty())
295 print_resall_header(out, rtpDBEntry);
297 for (const auto& r : rtpDBEntry)
301 print_resatoms(out, atype, r);
302 for (auto bt : gmx::EnumerationWrapper<BondedTypes>{})
304 print_resbondeds(out, bt, r);
310 void readResidueDatabase(const std::string& rrdb,
311 std::vector<PreprocessResidue>* rtpDBEntry,
312 PreprocessingAtomTypes* atype,
314 const gmx::MDLogger& logger,
315 bool bAllowOverrideRTP)
318 char filebase[STRLEN], line[STRLEN], header[STRLEN];
320 int dum1, dum2, dum3;
321 bool bNextResidue, bError;
323 fflib_filename_base(rrdb.c_str(), filebase, STRLEN);
325 in = fflib_open(rrdb);
327 PreprocessResidue header_settings;
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 */
338 header_settings.bKeepAllGeneratedDihedrals = FALSE;
339 header_settings.nrexcl = 3;
340 header_settings.bGenerateHH14Interactions = TRUE;
341 header_settings.bRemoveDihedralIfWithImproper = TRUE;
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
360 get_a_line(in, line, STRLEN);
361 if (!get_header(line, header))
363 gmx_fatal(FARGS, "in .rtp file at line:\n%s\n", line);
365 if (gmx::equalCaseInsensitive("bondedtypes", header, 5))
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,
375 &header_settings.nrexcl,
381 "need 4 to 8 parameters in the header of .rtp file %s at line:\n%s\n",
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);
393 .appendTextFormatted("Using default: not generating all possible dihedrals");
394 header_settings.bKeepAllGeneratedDihedrals = FALSE;
400 .appendTextFormatted("Using default: excluding 3 bonded neighbors");
401 header_settings.nrexcl = 3;
407 .appendTextFormatted("Using default: generating 1,4 H--H interactions");
408 header_settings.bGenerateHH14Interactions = TRUE;
412 GMX_LOG(logger.warning)
414 .appendTextFormatted(
415 "Using default: removing proper dihedrals found on the same bond as a "
417 header_settings.bRemoveDihedralIfWithImproper = TRUE;
422 GMX_LOG(logger.warning)
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));
429 /* We don't know the current size of rrtp, but simply realloc immediately */
430 auto oldArrayEnd = rtpDBEntry->end();
433 /* Initialise rtp entry structure */
434 rtpDBEntry->push_back(header_settings);
435 PreprocessResidue* res = &rtpDBEntry->back();
436 if (!get_header(line, header))
438 gmx_fatal(FARGS, "in .rtp file at line:\n%s\n", line);
440 res->resname = header;
441 res->filebase = filebase;
443 get_a_line(in, line, STRLEN);
445 bNextResidue = FALSE;
448 if (!get_header(line, header))
454 auto bt = get_bt(header);
457 /* header is an bonded directive */
458 bError = !read_bondeds(*bt, in, line, res);
460 else if (gmx::equalCaseInsensitive("atoms", header, 5))
462 /* header is the atoms directive */
463 bError = !read_atoms(in, line, res, tab, atype);
467 /* else header must be a residue name */
473 gmx_fatal(FARGS, "in .rtp file in residue %s at line:\n%s\n", res->resname.c_str(), line);
475 } while ((feof(in) == 0) && !bNextResidue);
477 if (res->natom() == 0)
479 gmx_fatal(FARGS, "No atoms found in .rtp file in residue %s\n", res->resname.c_str());
482 auto found = std::find_if(
483 rtpDBEntry->begin(), rtpDBEntry->end() - 1, [&res](const PreprocessResidue& entry) {
484 return gmx::equalCaseInsensitive(entry.resname, res->resname);
487 if (found != rtpDBEntry->end() - 1)
489 if (found >= oldArrayEnd)
491 gmx_fatal(FARGS, "Found a second entry for '%s' in '%s'", res->resname.c_str(), rrdb.c_str());
493 if (bAllowOverrideRTP)
495 GMX_LOG(logger.warning)
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(),
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.
506 rtpDBEntry->erase(rtpDBEntry->end() - 1);
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(),
521 std::sort(rtpDBEntry->begin(), rtpDBEntry->end(), [](const PreprocessResidue& a, const PreprocessResidue& b) {
522 return std::lexicographical_compare(
527 [](const char& c1, const char& c2) { return std::toupper(c1) < std::toupper(c2); });
530 check_rtp(*rtpDBEntry, rrdb, logger);
533 /************************************************************
537 ***********************************************************/
538 static bool is_sign(char c)
540 return (c == '+' || c == '-');
543 /* Compares if the strins match except for a sign at the end
544 * of (only) one of the two.
546 static int neq_str_sign(const char* a1, const char* a2)
550 l1 = static_cast<int>(strlen(a1));
551 l2 = static_cast<int>(strlen(a2));
552 lm = std::min(l1, l2);
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))
565 std::string searchResidueDatabase(const std::string& key,
566 gmx::ArrayRef<const PreprocessResidue> rtpDBEntry,
567 const gmx::MDLogger& logger)
569 int nbest, best, besti;
574 /* We want to match at least one character */
576 for (auto it = rtpDBEntry.begin(); it != rtpDBEntry.end(); it++)
578 if (gmx::equalCaseInsensitive(key, it->resname))
580 besti = std::distance(rtpDBEntry.begin(), it);
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()))
594 bestbuf = rtpDBEntry[besti].resname;
598 bestbuf.append(" or ");
599 bestbuf.append(it->resname);
606 besti = std::distance(rtpDBEntry.begin(), it);
615 "Residue '%s' not found in residue topology database, looks a bit like %s",
619 else if (besti == -1)
621 gmx_fatal(FARGS, "Residue '%s' not found in residue topology database", key.c_str());
623 if (!gmx::equalCaseInsensitive(rtpDBEntry[besti].resname, key))
625 GMX_LOG(logger.warning)
627 .appendTextFormatted(
628 "'%s' not found in residue topology database, "
629 "trying to use '%s'",
631 rtpDBEntry[besti].resname.c_str());
634 return rtpDBEntry[besti].resname;
637 gmx::ArrayRef<const PreprocessResidue>::const_iterator
638 getDatabaseEntry(const std::string& rtpname, gmx::ArrayRef<const PreprocessResidue> rtpDBEntry)
640 auto found = std::find_if(
641 rtpDBEntry.begin(), rtpDBEntry.end(), [&rtpname](const PreprocessResidue& entry) {
642 return gmx::equalCaseInsensitive(rtpname, entry.resname);
644 if (found == rtpDBEntry.end())
646 /* This should never happen, since searchResidueDatabase should have been called
647 * before calling getDatabaseEntry.
649 gmx_fatal(FARGS, "Residue type '%s' not found in residue topology database", rtpname.c_str());