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,2018,2019, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
49 #include "gromacs/gmxpreprocess/fflibutil.h"
50 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
51 #include "gromacs/gmxpreprocess/grompp_impl.h"
52 #include "gromacs/gmxpreprocess/notset.h"
53 #include "gromacs/gmxpreprocess/pgutil.h"
54 #include "gromacs/topology/atoms.h"
55 #include "gromacs/topology/symtab.h"
56 #include "gromacs/utility/cstringutil.h"
57 #include "gromacs/utility/fatalerror.h"
58 #include "gromacs/utility/futil.h"
59 #include "gromacs/utility/smalloc.h"
60 #include "gromacs/utility/strdb.h"
62 #include "hackblock.h"
64 PreprocessingAtomTypes read_atype(const char *ffdir, t_symtab *tab)
67 char buf[STRLEN], name[STRLEN];
72 std::vector<std::string> files = fflib_search_file_end(ffdir, ".atp", TRUE);
74 PreprocessingAtomTypes at;
76 for (const auto &filename : files)
78 in = fflib_open(filename);
81 /* Skip blank or comment-only lines */
84 if (fgets2(buf, STRLEN, in) != nullptr)
90 while ((feof(in) == 0) && strlen(buf) == 0);
92 if (sscanf(buf, "%s%lf", name, &m) == 2)
95 at.addType(tab, *a, name, InteractionType({}, {}), 0, 0);
96 fprintf(stderr, "\rAtomtype %d", ++nratt);
101 fprintf(stderr, "\nInvalid format: %s\n", buf);
106 fprintf(stderr, "\n");
112 static void print_resatoms(FILE *out, const PreprocessingAtomTypes &atype, const PreprocessResidue &rtpDBEntry)
114 fprintf(out, "[ %s ]\n", rtpDBEntry.resname.c_str());
115 fprintf(out, " [ atoms ]\n");
117 for (int j = 0; (j < rtpDBEntry.natom()); j++)
119 int tp = rtpDBEntry.atom[j].type;
120 const char *tpnm = atype.atomNameFromAtomType(tp);
123 gmx_fatal(FARGS, "Incorrect atomtype (%d)", tp);
125 fprintf(out, "%6s %6s %8.3f %6d\n",
126 *(rtpDBEntry.atomname[j]), tpnm, rtpDBEntry.atom[j].q, rtpDBEntry.cgnr[j]);
130 static bool read_atoms(FILE *in, char *line,
131 PreprocessResidue *r0, t_symtab *tab, PreprocessingAtomTypes *atype)
134 char buf[256], buf1[256];
139 r0->atomname.clear();
141 while (get_a_line(in, line, STRLEN) && (strchr(line, '[') == nullptr))
143 if (sscanf(line, "%s%s%lf%d", buf, buf1, &q, &cg) != 4)
147 r0->atomname.push_back(put_symtab(tab, buf));
148 r0->atom.emplace_back();
149 r0->atom.back().q = q;
150 r0->cgnr.push_back(cg);
151 int j = atype->atomTypeFromName(buf1);
154 gmx_fatal(FARGS, "Atom type %s (residue %s) not found in atomtype "
155 "database", buf1, r0->resname.c_str());
157 r0->atom.back().type = j;
158 r0->atom.back().m = atype->atomMassAFromAtomType(j);
164 static bool read_bondeds(int bt, FILE *in, char *line, PreprocessResidue *rtpDBEntry)
168 while (get_a_line(in, line, STRLEN) && (strchr(line, '[') == nullptr))
172 rtpDBEntry->rb[bt].b.emplace_back();
173 BondedInteraction *newBond = &rtpDBEntry->rb[bt].b.back();
174 for (int j = 0; j < btsNiatoms[bt]; j++)
176 if (sscanf(line+n, "%s%n", str, &ni) == 1)
186 while (isspace(line[n]))
197 static void print_resbondeds(FILE *out, int bt, const PreprocessResidue &rtpDBEntry)
199 if (!rtpDBEntry.rb[bt].b.empty())
201 fprintf(out, " [ %s ]\n", btsNames[bt]);
203 for (const auto &b : rtpDBEntry.rb[bt].b)
205 for (int j = 0; j < btsNiatoms[bt]; j++)
207 fprintf(out, "%6s ", b.a[j].c_str());
211 fprintf(out, " %s", b.s.c_str());
218 static void check_rtp(gmx::ArrayRef<const PreprocessResidue> rtpDBEntry, const std::string &libfn)
220 /* check for double entries, assuming list is already sorted */
221 for (auto it = rtpDBEntry.begin() + 1; it != rtpDBEntry.end(); it++)
224 if (gmx::equalCaseInsensitive(prev->resname, it->resname))
226 fprintf(stderr, "WARNING double entry %s in file %s\n",
227 it->resname.c_str(), libfn.c_str());
232 static int get_bt(char* header)
236 for (i = 0; i < ebtsNR; i++)
238 if (gmx_strcasecmp(btsNames[i], header) == 0)
246 /* print all the ebtsNR type numbers */
247 static void print_resall_header(FILE *out, gmx::ArrayRef<const PreprocessResidue> rtpDBEntry)
249 fprintf(out, "[ bondedtypes ]\n");
250 fprintf(out, "; bonds angles dihedrals impropers all_dihedrals nr_exclusions HH14 remove_dih\n");
251 fprintf(out, " %5d %6d %9d %9d %14d %14d %14d %14d\n\n",
252 rtpDBEntry[0].rb[0].type,
253 rtpDBEntry[0].rb[1].type,
254 rtpDBEntry[0].rb[2].type,
255 rtpDBEntry[0].rb[3].type,
256 static_cast<int>(rtpDBEntry[0].bKeepAllGeneratedDihedrals),
257 rtpDBEntry[0].nrexcl,
258 static_cast<int>(rtpDBEntry[0].bGenerateHH14Interactions),
259 static_cast<int>(rtpDBEntry[0].bRemoveDihedralIfWithImproper));
262 void print_resall(FILE *out, gmx::ArrayRef<const PreprocessResidue> rtpDBEntry,
263 const PreprocessingAtomTypes &atype)
265 if (rtpDBEntry.empty())
270 print_resall_header(out, rtpDBEntry);
272 for (const auto &r : rtpDBEntry)
276 print_resatoms(out, atype, r);
277 for (int bt = 0; bt < ebtsNR; bt++)
279 print_resbondeds(out, bt, r);
285 void readResidueDatabase(const std::string &rrdb, std::vector<PreprocessResidue> *rtpDBEntry,
286 PreprocessingAtomTypes *atype, t_symtab *tab,
287 bool bAllowOverrideRTP)
290 char filebase[STRLEN], line[STRLEN], header[STRLEN];
292 int dum1, dum2, dum3;
293 bool bNextResidue, bError;
295 fflib_filename_base(rrdb.c_str(), filebase, STRLEN);
297 in = fflib_open(rrdb);
299 PreprocessResidue header_settings;
301 /* these bonded parameters will overwritten be when *
302 * there is a [ bondedtypes ] entry in the .rtp file */
303 header_settings.rb[ebtsBONDS].type = 1; /* normal bonds */
304 header_settings.rb[ebtsANGLES].type = 1; /* normal angles */
305 header_settings.rb[ebtsPDIHS].type = 1; /* normal dihedrals */
306 header_settings.rb[ebtsIDIHS].type = 2; /* normal impropers */
307 header_settings.rb[ebtsEXCLS].type = 1; /* normal exclusions */
308 header_settings.rb[ebtsCMAP].type = 1; /* normal cmap torsions */
310 header_settings.bKeepAllGeneratedDihedrals = FALSE;
311 header_settings.nrexcl = 3;
312 header_settings.bGenerateHH14Interactions = TRUE;
313 header_settings.bRemoveDihedralIfWithImproper = TRUE;
315 /* Column 5 & 6 aren't really bonded types, but we include
316 * them here to avoid introducing a new section:
317 * Column 5 : This controls the generation of dihedrals from the bonding.
318 * All possible dihedrals are generated automatically. A value of
319 * 1 here means that all these are retained. A value of
320 * 0 here requires generated dihedrals be removed if
321 * * there are any dihedrals on the same central atoms
322 * specified in the residue topology, or
323 * * there are other identical generated dihedrals
324 * sharing the same central atoms, or
325 * * there are other generated dihedrals sharing the
326 * same central bond that have fewer hydrogen atoms
327 * Column 6: Number of bonded neighbors to exclude.
328 * Column 7: Generate 1,4 interactions between two hydrogen atoms
329 * Column 8: Remove proper dihedrals if centered on the same bond
330 * as an improper dihedral
332 get_a_line(in, line, STRLEN);
333 if (!get_header(line, header))
335 gmx_fatal(FARGS, "in .rtp file at line:\n%s\n", line);
337 if (gmx_strncasecmp("bondedtypes", header, 5) == 0)
339 get_a_line(in, line, STRLEN);
340 if ((nparam = sscanf(line, "%d %d %d %d %d %d %d %d",
341 &header_settings.rb[ebtsBONDS].type, &header_settings.rb[ebtsANGLES].type,
342 &header_settings.rb[ebtsPDIHS].type, &header_settings.rb[ebtsIDIHS].type,
343 &dum1, &header_settings.nrexcl, &dum2, &dum3)) < 4)
345 gmx_fatal(FARGS, "need 4 to 8 parameters in the header of .rtp file %s at line:\n%s\n", rrdb.c_str(), line);
347 header_settings.bKeepAllGeneratedDihedrals = (dum1 != 0);
348 header_settings.bGenerateHH14Interactions = (dum2 != 0);
349 header_settings.bRemoveDihedralIfWithImproper = (dum3 != 0);
350 get_a_line(in, line, STRLEN);
353 fprintf(stderr, "Using default: not generating all possible dihedrals\n");
354 header_settings.bKeepAllGeneratedDihedrals = FALSE;
358 fprintf(stderr, "Using default: excluding 3 bonded neighbors\n");
359 header_settings.nrexcl = 3;
363 fprintf(stderr, "Using default: generating 1,4 H--H interactions\n");
364 header_settings.bGenerateHH14Interactions = TRUE;
368 fprintf(stderr, "Using default: removing proper dihedrals found on the same bond as a proper dihedral\n");
369 header_settings.bRemoveDihedralIfWithImproper = TRUE;
375 "Reading .rtp file without '[ bondedtypes ]' directive,\n"
376 "Will proceed as if the entry was:\n");
377 print_resall_header(stderr, gmx::arrayRefFromArray(&header_settings, 1));
379 /* We don't know the current size of rrtp, but simply realloc immediately */
380 auto oldArrayEnd = rtpDBEntry->end() - 1;
383 /* Initialise rtp entry structure */
384 rtpDBEntry->push_back(header_settings);
385 PreprocessResidue *res = &rtpDBEntry->back();
386 if (!get_header(line, header))
388 gmx_fatal(FARGS, "in .rtp file at line:\n%s\n", line);
390 res->resname = header;
391 res->filebase = filebase;
393 get_a_line(in, line, STRLEN);
395 bNextResidue = FALSE;
398 if (!get_header(line, header))
407 /* header is an bonded directive */
408 bError = !read_bondeds(bt, in, line, res);
410 else if (gmx_strncasecmp("atoms", header, 5) == 0)
412 /* header is the atoms directive */
413 bError = !read_atoms(in, line, res, tab, atype);
417 /* else header must be a residue name */
423 gmx_fatal(FARGS, "in .rtp file in residue %s at line:\n%s\n",
424 res->resname.c_str(), line);
427 while ((feof(in) == 0) && !bNextResidue);
429 if (res->natom() == 0)
431 gmx_fatal(FARGS, "No atoms found in .rtp file in residue %s\n",
432 res->resname.c_str());
435 auto found = std::find_if(rtpDBEntry->begin(), rtpDBEntry->end()-1,
436 [&res](const PreprocessResidue &entry)
437 { return gmx::equalCaseInsensitive(res->resname, entry.resname); });
439 if (found == rtpDBEntry->end() - 1)
441 int size = rtpDBEntry->size() - 1;
442 fprintf(stderr, "\rResidue %d", size);
447 if (found >= oldArrayEnd)
449 gmx_fatal(FARGS, "Found a second entry for '%s' in '%s'",
450 res->resname.c_str(), rrdb.c_str());
452 if (bAllowOverrideRTP)
454 fprintf(stderr, "\n");
455 fprintf(stderr, "Found another rtp entry for '%s' in '%s', ignoring this entry and keeping the one from '%s.rtp'\n",
456 res->resname.c_str(), rrdb.c_str(), found->filebase.c_str());
457 /* We should free all the data for this entry.
458 * The current code gives a lot of dangling pointers.
460 rtpDBEntry->erase(rtpDBEntry->end() - 1);
464 gmx_fatal(FARGS, "Found rtp entries for '%s' in both '%s' and '%s'. If you want the first definition to override the second one, set the -rtpo option of pdb2gmx.", res->resname.c_str(), found->filebase.c_str(), rrdb.c_str());
470 fprintf(stderr, "\nSorting it all out...\n");
471 std::sort(rtpDBEntry->begin(), rtpDBEntry->end(), []
472 (const PreprocessResidue &a,
473 const PreprocessResidue &b)
474 {return std::lexicographical_compare(
475 a.resname.begin(), a.resname.end(),
476 b.resname.begin(), b.resname.end(),
477 [](const char &c1, const char &c2)
478 { return std::toupper(c1) < std::toupper(c2); }); });
480 check_rtp(*rtpDBEntry, rrdb);
483 /************************************************************
487 ***********************************************************/
488 static bool is_sign(char c)
490 return (c == '+' || c == '-');
493 /* Compares if the strins match except for a sign at the end
494 * of (only) one of the two.
496 static int neq_str_sign(const char *a1, const char *a2)
500 l1 = static_cast<int>(strlen(a1));
501 l2 = static_cast<int>(strlen(a2));
502 lm = std::min(l1, l2);
505 ((l1 == l2+1 && is_sign(a1[l1-1])) ||
506 (l2 == l1+1 && is_sign(a2[l2-1]))) &&
507 gmx_strncasecmp(a1, a2, lm) == 0)
517 std::string searchResidueDatabase(const std::string &key, gmx::ArrayRef<const PreprocessResidue> rtpDBEntry)
519 int nbest, best, besti;
524 /* We want to match at least one character */
526 for (auto it = rtpDBEntry.begin(); it != rtpDBEntry.end(); it++)
528 if (gmx::equalCaseInsensitive(key, it->resname))
530 besti = std::distance(rtpDBEntry.begin(), it);
536 /* Allow a mismatch of at most a sign character (with warning) */
537 int n = neq_str_sign(key.c_str(), it->resname.c_str());
539 n+1 >= gmx::index(key.length()) &&
540 n+1 >= gmx::index(it->resname.length()))
546 bestbuf = rtpDBEntry[besti].resname;
550 bestbuf.append(" or ");
551 bestbuf.append(it->resname);
558 besti = std::distance(rtpDBEntry.begin(), it);
566 gmx_fatal(FARGS, "Residue '%s' not found in residue topology database, looks a bit like %s", key.c_str(), bestbuf.c_str());
568 else if (besti == -1)
570 gmx_fatal(FARGS, "Residue '%s' not found in residue topology database", key.c_str());
572 if (!gmx::equalCaseInsensitive(rtpDBEntry[besti].resname, key))
575 "\nWARNING: '%s' not found in residue topology database, "
576 "trying to use '%s'\n\n", key.c_str(), rtpDBEntry[besti].resname.c_str());
579 return rtpDBEntry[besti].resname;
582 gmx::ArrayRef<const PreprocessResidue>::const_iterator
583 getDatabaseEntry(const std::string &rtpname, gmx::ArrayRef<const PreprocessResidue> rtpDBEntry)
585 auto found = std::find_if(rtpDBEntry.begin(), rtpDBEntry.end(),
586 [&rtpname](const PreprocessResidue &entry)
587 { return gmx::equalCaseInsensitive(rtpname, entry.resname); });
588 if (found == rtpDBEntry.end())
590 /* This should never happen, since searchResidueDatabase should have been called
591 * before calling getDatabaseEntry.
593 gmx_fatal(FARGS, "Residue type '%s' not found in residue topology database", rtpname.c_str());