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/fileio/pdbio.h"
52 #include "gromacs/gmxpreprocess/add_par.h"
53 #include "gromacs/gmxpreprocess/fflibutil.h"
54 #include "gromacs/gmxpreprocess/gen_ad.h"
55 #include "gromacs/gmxpreprocess/gen_vsite.h"
56 #include "gromacs/gmxpreprocess/grompp_impl.h"
57 #include "gromacs/gmxpreprocess/h_db.h"
58 #include "gromacs/gmxpreprocess/notset.h"
59 #include "gromacs/gmxpreprocess/pgutil.h"
60 #include "gromacs/gmxpreprocess/specbond.h"
61 #include "gromacs/gmxpreprocess/topdirs.h"
62 #include "gromacs/gmxpreprocess/topio.h"
63 #include "gromacs/gmxpreprocess/toputil.h"
64 #include "gromacs/math/functions.h"
65 #include "gromacs/math/vec.h"
66 #include "gromacs/topology/residuetypes.h"
67 #include "gromacs/topology/symtab.h"
68 #include "gromacs/utility/binaryinformation.h"
69 #include "gromacs/utility/cstringutil.h"
70 #include "gromacs/utility/dir_separator.h"
71 #include "gromacs/utility/exceptions.h"
72 #include "gromacs/utility/fatalerror.h"
73 #include "gromacs/utility/futil.h"
74 #include "gromacs/utility/logger.h"
75 #include "gromacs/utility/niceheader.h"
76 #include "gromacs/utility/path.h"
77 #include "gromacs/utility/programcontext.h"
78 #include "gromacs/utility/smalloc.h"
79 #include "gromacs/utility/strconvert.h"
80 #include "gromacs/utility/strdb.h"
81 #include "gromacs/utility/stringutil.h"
82 #include "gromacs/utility/textwriter.h"
84 #include "hackblock.h"
87 /* this must correspond to enum in pdb2top.h */
88 const char* hh[ehisNR] = { "HISD", "HISE", "HISH", "HIS1" };
90 static int missing_atoms(const PreprocessResidue* rp, int resind, t_atoms* at, int i0, int i, const gmx::MDLogger& logger)
93 for (int j = 0; j < rp->natom(); j++)
95 const char* name = *(rp->atomname[j]);
97 for (int k = i0; k < i; k++)
99 bFound = (bFound || (gmx_strcasecmp(*(at->atomname[k]), name) == 0));
104 GMX_LOG(logger.warning)
106 .appendTextFormatted("atom %s is missing in residue %s %d in the pdb file",
108 *(at->resinfo[resind].name),
109 at->resinfo[resind].nr);
110 if (name[0] == 'H' || name[0] == 'h')
112 GMX_LOG(logger.warning)
114 .appendTextFormatted(
115 "You might need to add atom %s to the hydrogen database of "
117 "in the file %s.hdb (see the manual)",
119 *(at->resinfo[resind].rtp),
120 rp->filebase.c_str());
128 bool is_int(double x)
130 const double tol = 1e-4;
137 ix = gmx::roundToInt(x);
139 return (fabs(x - ix) < tol);
142 static void choose_ff_impl(const char* ffsel,
147 const gmx::MDLogger& logger)
149 std::vector<gmx::DataFileInfo> ffdirs = fflib_enumerate_forcefields();
150 const int nff = ssize(ffdirs);
152 /* Replace with unix path separators */
153 #if DIR_SEPARATOR != '/'
154 for (int i = 0; i < nff; ++i)
156 std::replace(ffdirs[i].dir.begin(), ffdirs[i].dir.end(), DIR_SEPARATOR, '/');
160 /* Store the force field names in ffs */
161 std::vector<std::string> ffs;
162 ffs.reserve(ffdirs.size());
163 for (int i = 0; i < nff; ++i)
165 ffs.push_back(gmx::stripSuffixIfPresent(ffdirs[i].name, fflib_forcefield_dir_ext()));
169 if (ffsel != nullptr)
174 for (int i = 0; i < nff; ++i)
178 /* Matching ff name */
182 if (ffdirs[i].dir == ".")
198 GMX_LOG(logger.warning)
200 .appendTextFormatted(
201 "Force field '%s' occurs in %d places. pdb2gmx is using the one in "
202 "the current directory. Use interactive selection "
203 "(not the -ff option) if you would prefer a different one.",
209 std::string message = gmx::formatString(
210 "Force field '%s' occurs in %d places, but not in "
211 "the current directory.\n"
212 "Run without the -ff switch and select the force "
213 "field interactively.",
216 GMX_THROW(gmx::InconsistentInputError(message));
219 else if (nfound == 0)
221 std::string message = gmx::formatString(
222 "Could not find force field '%s' in current directory, "
223 "install tree or GMXLIB path.",
225 GMX_THROW(gmx::InconsistentInputError(message));
230 std::vector<std::string> desc;
231 desc.reserve(ffdirs.size());
232 for (int i = 0; i < nff; ++i)
234 std::string docFileName(gmx::Path::join(ffdirs[i].dir, ffdirs[i].name, fflib_forcefield_doc()));
235 // TODO: Just try to open the file with a method that does not
236 // throw/bail out with a fatal error instead of multiple checks.
237 if (gmx::File::exists(docFileName, gmx::File::returnFalseOnError))
239 // TODO: Use a C++ API without such an intermediate/fixed-length buffer.
241 /* We don't use fflib_open, because we don't want printf's */
242 FILE* fp = gmx_ffopen(docFileName, "r");
243 get_a_line(fp, buf, STRLEN);
245 desc.emplace_back(buf);
249 desc.push_back(ffs[i]);
252 /* Order force fields from the same dir alphabetically
253 * and put deprecated force fields at the end.
255 for (int i = 0; i < nff; ++i)
257 for (int j = i + 1; j < nff; ++j)
259 if (ffdirs[i].dir == ffdirs[j].dir
260 && ((desc[i][0] == '[' && desc[j][0] != '[')
261 || ((desc[i][0] == '[' || desc[j][0] != '[')
262 && gmx_strcasecmp(desc[i].c_str(), desc[j].c_str()) > 0)))
264 std::swap(ffdirs[i].name, ffdirs[j].name);
265 std::swap(ffs[i], ffs[j]);
266 std::swap(desc[i], desc[j]);
271 GMX_LOG(logger.info).asParagraph().appendTextFormatted("Select the Force Field:");
272 for (int i = 0; i < nff; ++i)
274 if (i == 0 || ffdirs[i - 1].dir != ffdirs[i].dir)
276 if (ffdirs[i].dir == ".")
280 .appendTextFormatted("From current directory:");
286 .appendTextFormatted("From '%s':", ffdirs[i].dir.c_str());
289 GMX_LOG(logger.info).asParagraph().appendTextFormatted("%2d: %s", i + 1, desc[i].c_str());
293 // TODO: Add a C++ API for this.
298 pret = fgets(buf, STRLEN, stdin);
302 sel = strtol(buf, nullptr, 10);
305 } while (pret == nullptr || (sel < 0) || (sel >= nff));
307 /* Check for a current limitation of the fflib code.
308 * It will always read from the first ff directory in the list.
309 * This check assumes that the order of ffs matches the order
310 * in which fflib_open searches ff library files.
312 for (int i = 0; i < sel; i++)
314 if (ffs[i] == ffs[sel])
316 std::string message = gmx::formatString(
317 "Can only select the first of multiple force "
318 "field entries with directory name '%s%s' in "
319 "the list. If you want to use the next entry, "
320 "run pdb2gmx in a different directory, set GMXLIB "
321 "to point to the desired force field first, and/or "
322 "rename or move the force field directory present "
323 "in the current working directory.",
325 fflib_forcefield_dir_ext());
326 GMX_THROW(gmx::NotImplementedError(message));
335 if (ffs[sel].length() >= static_cast<size_t>(ff_maxlen))
337 std::string message = gmx::formatString("Length of force field name (%d) >= maxlen (%d)",
338 static_cast<int>(ffs[sel].length()),
340 GMX_THROW(gmx::InvalidInputError(message));
342 strcpy(forcefield, ffs[sel].c_str());
345 if (ffdirs[sel].bFromDefaultDir)
347 ffpath = ffdirs[sel].name;
351 ffpath = gmx::Path::join(ffdirs[sel].dir, ffdirs[sel].name);
353 if (ffpath.length() >= static_cast<size_t>(ffdir_maxlen))
355 std::string message = gmx::formatString("Length of force field dir (%d) >= maxlen (%d)",
356 static_cast<int>(ffpath.length()),
358 GMX_THROW(gmx::InvalidInputError(message));
360 strcpy(ffdir, ffpath.c_str());
363 void choose_ff(const char* ffsel, char* forcefield, int ff_maxlen, char* ffdir, int ffdir_maxlen, const gmx::MDLogger& logger)
367 choose_ff_impl(ffsel, forcefield, ff_maxlen, ffdir, ffdir_maxlen, logger);
369 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
372 void choose_watermodel(const char* wmsel, const char* ffdir, char** watermodel, const gmx::MDLogger& logger)
374 const char* fn_watermodels = "watermodels.dat";
381 if (strcmp(wmsel, "none") == 0)
383 *watermodel = nullptr;
387 else if (strcmp(wmsel, "select") != 0)
389 *watermodel = gmx_strdup(wmsel);
394 std::string filename = gmx::Path::join(ffdir, fn_watermodels);
395 if (!fflib_fexist(filename))
397 GMX_LOG(logger.warning)
399 .appendTextFormatted("No file '%s' found, will not include a water model", fn_watermodels);
400 *watermodel = nullptr;
405 fp = fflib_open(filename);
406 GMX_LOG(logger.info).asParagraph().appendTextFormatted("Select the Water Model:");
409 while (get_a_line(fp, buf, STRLEN))
411 srenew(model, nwm + 1);
412 snew(model[nwm], STRLEN);
413 sscanf(buf, "%s%n", model[nwm], &i);
417 GMX_LOG(logger.info).asParagraph().appendTextFormatted("%2d: %s", nwm + 1, buf + i);
426 GMX_LOG(logger.info).asParagraph().appendTextFormatted("%2d: %s", nwm + 1, "None");
431 pret = fgets(buf, STRLEN, stdin);
435 sel = strtol(buf, nullptr, 10);
438 } while (pret == nullptr || sel < 0 || sel > nwm);
442 *watermodel = nullptr;
446 *watermodel = gmx_strdup(model[sel]);
449 for (i = 0; i < nwm; i++)
456 static int name2type(t_atoms* at,
458 gmx::ArrayRef<const PreprocessResidue> usedPpResidues,
460 const gmx::MDLogger& logger)
462 int i, j, prevresind, i0, prevcg, cg, curcg;
478 for (i = 0; (i < at->nr); i++)
481 if (at->atom[i].resind != resind)
483 resind = at->atom[i].resind;
484 bool bProt = rt->namedResidueHasType(*(at->resinfo[resind].name), "Protein");
485 bNterm = bProt && (resind == 0);
488 nmissat += missing_atoms(&usedPpResidues[prevresind], prevresind, at, i0, i, logger);
492 if (at->atom[i].m == 0)
496 name = *(at->atomname[i]);
497 j = search_jtype(usedPpResidues[resind], name, bNterm);
498 at->atom[i].type = usedPpResidues[resind].atom[j].type;
499 at->atom[i].q = usedPpResidues[resind].atom[j].q;
500 at->atom[i].m = usedPpResidues[resind].atom[j].m;
501 cg = usedPpResidues[resind].cgnr[j];
502 /* A charge group number -1 signals a separate charge group
505 if ((cg == -1) || (cg != prevcg) || (resind != prevresind))
521 at->atom[i].typeB = at->atom[i].type;
522 at->atom[i].qB = at->atom[i].q;
523 at->atom[i].mB = at->atom[i].m;
525 nmissat += missing_atoms(&usedPpResidues[resind], resind, at, i0, i, logger);
530 static void print_top_heavy_H(FILE* out, real mHmult)
534 fprintf(out, "; Using deuterium instead of hydrogen\n\n");
536 else if (mHmult == 4.0)
538 fprintf(out, "#define HEAVY_H\n\n");
540 else if (mHmult != 1.0)
543 "WARNING: unsupported proton mass multiplier (%g) "
549 void print_top_comment(FILE* out, const char* filename, const char* ffdir, bool bITP)
551 char ffdir_parent[STRLEN];
556 gmx::TextWriter writer(out);
557 gmx::niceHeader(&writer, filename, ';');
558 writer.writeLine(gmx::formatString(";\tThis is a %s topology file", bITP ? "include" : "standalone"));
559 writer.writeLine(";");
561 gmx::BinaryInformationSettings settings;
562 settings.generatedByHeader(true);
563 settings.linePrefix(";\t");
564 gmx::printBinaryInformation(&writer, gmx::getProgramContext(), settings);
566 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
568 if (strchr(ffdir, '/') == nullptr)
570 fprintf(out, ";\tForce field was read from the standard GROMACS share directory.\n;\n\n");
572 else if (ffdir[0] == '.')
575 ";\tForce field was read from current directory or a relative path - path "
580 strncpy(ffdir_parent, ffdir, STRLEN - 1);
581 ffdir_parent[STRLEN - 1] = '\0'; /*make sure it is 0-terminated even for long string*/
582 p = strrchr(ffdir_parent, '/');
587 ";\tForce field data was read from:\n"
591 ";\tThis might be a non-standard force field location. When you use this topology, "
593 ";\tforce field must either be present in the current directory, or the location\n"
594 ";\tspecified in the GMXLIB path variable or with the 'include' mdp file "
600 void print_top_header(FILE* out, const char* filename, bool bITP, const char* ffdir, real mHmult)
604 print_top_comment(out, filename, ffdir, bITP);
606 print_top_heavy_H(out, mHmult);
607 fprintf(out, "; Include forcefield parameters\n");
609 p = strrchr(ffdir, '/');
610 p = (ffdir[0] == '.' || p == nullptr) ? ffdir : p + 1;
612 fprintf(out, "#include \"%s/%s\"\n\n", p, fflib_forcefield_itp());
615 static void print_top_posre(FILE* out, const char* pr)
617 fprintf(out, "; Include Position restraint file\n");
618 fprintf(out, "#ifdef POSRES\n");
619 fprintf(out, "#include \"%s\"\n", pr);
620 fprintf(out, "#endif\n\n");
623 static void print_top_water(FILE* out, const char* ffdir, const char* water)
628 fprintf(out, "; Include water topology\n");
630 p = strrchr(ffdir, '/');
631 p = (ffdir[0] == '.' || p == nullptr) ? ffdir : p + 1;
632 fprintf(out, "#include \"%s/%s.itp\"\n", p, water);
635 fprintf(out, "#ifdef POSRES_WATER\n");
636 fprintf(out, "; Position restraint for each water oxygen\n");
637 fprintf(out, "[ position_restraints ]\n");
638 fprintf(out, ";%3s %5s %9s %10s %10s\n", "i", "funct", "fcx", "fcy", "fcz");
639 fprintf(out, "%4d %4d %10g %10g %10g\n", 1, 1, 1000.0, 1000.0, 1000.0);
640 fprintf(out, "#endif\n");
643 sprintf(buf, "%s/ions.itp", p);
645 if (fflib_fexist(buf))
647 fprintf(out, "; Include topology for ions\n");
648 fprintf(out, "#include \"%s\"\n", buf);
653 static void print_top_system(FILE* out, const char* title)
655 fprintf(out, "[ %s ]\n", enumValueToString(Directive::d_system));
656 fprintf(out, "; Name\n");
657 fprintf(out, "%s\n\n", title[0] ? title : "Protein");
660 void print_top_mols(FILE* out,
664 gmx::ArrayRef<const std::string> incls,
665 gmx::ArrayRef<const t_mols> mols)
670 fprintf(out, "; Include chain topologies\n");
671 for (const auto& incl : incls)
673 fprintf(out, "#include \"%s\"\n", gmx::Path::getFilename(incl).c_str());
680 print_top_water(out, ffdir, water);
682 print_top_system(out, title);
686 fprintf(out, "[ %s ]\n", enumValueToString(Directive::d_molecules));
687 fprintf(out, "; %-15s %5s\n", "Compound", "#mols");
688 for (const auto& mol : mols)
690 fprintf(out, "%-15s %5d\n", mol.name.c_str(), mol.nr);
695 void write_top(FILE* out,
700 gmx::ArrayRef<const int> bts,
701 gmx::ArrayRef<const InteractionsOfType> plist,
703 PreprocessingAtomTypes* atype,
706 /* NOTE: nrexcl is not the size of *excl! */
708 if (at && atype && cgnr)
710 fprintf(out, "[ %s ]\n", enumValueToString(Directive::d_moleculetype));
711 fprintf(out, "; %-15s %5s\n", "Name", "nrexcl");
712 fprintf(out, "%-15s %5d\n\n", molname ? molname : "Protein", nrexcl);
714 print_atoms(out, atype, at, cgnr, bRTPresname);
716 out, at->nr, Directive::d_bonds, F_BONDS, bts[static_cast<int>(BondedTypes::Bonds)], plist);
717 print_bondeds(out, at->nr, Directive::d_constraints, F_CONSTR, 0, plist);
718 print_bondeds(out, at->nr, Directive::d_constraints, F_CONSTRNC, 0, plist);
719 print_bondeds(out, at->nr, Directive::d_pairs, F_LJ14, 0, plist);
720 print_excl(out, at->nr, excls);
722 out, at->nr, Directive::d_angles, F_ANGLES, bts[static_cast<int>(BondedTypes::Angles)], plist);
725 Directive::d_dihedrals,
727 bts[static_cast<int>(BondedTypes::ProperDihedrals)],
731 Directive::d_dihedrals,
733 bts[static_cast<int>(BondedTypes::ImproperDihedrals)],
735 print_bondeds(out, at->nr, Directive::d_cmap, F_CMAP, bts[static_cast<int>(BondedTypes::Cmap)], plist);
736 print_bondeds(out, at->nr, Directive::d_polarization, F_POLARIZATION, 0, plist);
737 print_bondeds(out, at->nr, Directive::d_thole_polarization, F_THOLE_POL, 0, plist);
738 print_bondeds(out, at->nr, Directive::d_vsites2, F_VSITE2, 0, plist);
739 print_bondeds(out, at->nr, Directive::d_vsites3, F_VSITE3, 0, plist);
740 print_bondeds(out, at->nr, Directive::d_vsites3, F_VSITE3FD, 0, plist);
741 print_bondeds(out, at->nr, Directive::d_vsites3, F_VSITE3FAD, 0, plist);
742 print_bondeds(out, at->nr, Directive::d_vsites3, F_VSITE3OUT, 0, plist);
743 print_bondeds(out, at->nr, Directive::d_vsites4, F_VSITE4FD, 0, plist);
744 print_bondeds(out, at->nr, Directive::d_vsites4, F_VSITE4FDN, 0, plist);
748 print_top_posre(out, pr);
754 static void do_ssbonds(InteractionsOfType* ps,
756 gmx::ArrayRef<const DisulfideBond> ssbonds,
759 for (const auto& bond : ssbonds)
761 int ri = bond.firstResidue;
762 int rj = bond.secondResidue;
763 int ai = search_res_atom(bond.firstAtom.c_str(), ri, atoms, "special bond", bAllowMissing);
764 int aj = search_res_atom(bond.secondAtom.c_str(), rj, atoms, "special bond", bAllowMissing);
765 if ((ai == -1) || (aj == -1))
768 "Trying to make impossible special bond (%s-%s)!",
769 bond.firstAtom.c_str(),
770 bond.secondAtom.c_str());
772 add_param(ps, ai, aj, {}, nullptr);
776 static void at2bonds(InteractionsOfType* psb,
777 gmx::ArrayRef<MoleculePatchDatabase> globalPatches,
779 gmx::ArrayRef<const gmx::RVec> x,
781 real short_bond_dist,
782 gmx::ArrayRef<const int> cyclicBondsIndex,
783 const gmx::MDLogger& logger)
785 real long_bond_dist2, short_bond_dist2;
788 long_bond_dist2 = gmx::square(long_bond_dist);
789 short_bond_dist2 = gmx::square(short_bond_dist);
800 GMX_LOG(logger.info).asParagraph().appendTextFormatted("Making bonds...");
802 for (int resind = 0; (resind < atoms->nres) && (i < atoms->nr); resind++)
804 /* add bonds from list of bonded interactions */
805 for (const auto& patch : globalPatches[resind].rb[BondedTypes::Bonds].b)
807 /* Unfortunately we can not issue errors or warnings
808 * for missing atoms in bonds, as the hydrogens and terminal atoms
809 * have not been added yet.
811 int ai = search_atom(patch.ai().c_str(), i, atoms, ptr, TRUE, cyclicBondsIndex);
812 int aj = search_atom(patch.aj().c_str(), i, atoms, ptr, TRUE, cyclicBondsIndex);
813 if (ai != -1 && aj != -1)
815 real dist2 = distance2(x[ai], x[aj]);
816 if (dist2 > long_bond_dist2)
819 GMX_LOG(logger.warning)
821 .appendTextFormatted(
822 "Long Bond (%d-%d = %g nm)", ai + 1, aj + 1, std::sqrt(dist2));
824 else if (dist2 < short_bond_dist2)
826 GMX_LOG(logger.warning)
828 .appendTextFormatted(
829 "Short Bond (%d-%d = %g nm)", ai + 1, aj + 1, std::sqrt(dist2));
831 add_param(psb, ai, aj, {}, patch.s.c_str());
834 /* add bonds from list of hacks (each added atom gets a bond) */
835 while ((i < atoms->nr) && (atoms->atom[i].resind == resind))
837 for (const auto& patch : globalPatches[resind].hack)
839 if ((patch.tp > 0 || patch.type() == MoleculePatchType::Add)
840 && patch.a[0] == *(atoms->atomname[i]))
844 case 9: /* COOH terminus */
845 add_param(psb, i, i + 1, {}, nullptr); /* C-O */
846 add_param(psb, i, i + 2, {}, nullptr); /* C-OA */
847 add_param(psb, i + 2, i + 3, {}, nullptr); /* OA-H */
850 for (int k = 0; (k < patch.nr); k++)
852 add_param(psb, i, i + k + 1, {}, nullptr);
859 /* we're now at the start of the next residue */
863 static bool pcompar(const InteractionOfType& a, const InteractionOfType& b)
867 if (((d = a.ai() - b.ai()) != 0) || ((d = a.aj() - b.aj()) != 0))
873 return a.interactionTypeName().length() > b.interactionTypeName().length();
877 static void clean_bonds(InteractionsOfType* ps, const gmx::MDLogger& logger)
882 for (auto& bond : ps->interactionTypes)
886 std::sort(ps->interactionTypes.begin(), ps->interactionTypes.end(), pcompar);
888 /* remove doubles, keep the first one always. */
889 int oldNumber = ps->size();
890 for (auto parm = ps->interactionTypes.begin() + 1; parm != ps->interactionTypes.end();)
892 auto prev = parm - 1;
893 if (parm->ai() == prev->ai() && parm->aj() == prev->aj())
895 parm = ps->interactionTypes.erase(parm);
904 .appendTextFormatted("Number of bonds was %d, now %zu", oldNumber, ps->size());
908 GMX_LOG(logger.info).asParagraph().appendTextFormatted("No bonds");
912 void print_sums(const t_atoms* atoms, bool bSystem, const gmx::MDLogger& logger)
920 where = " in system";
929 for (i = 0; (i < atoms->nr); i++)
931 m += atoms->atom[i].m;
932 qtot += atoms->atom[i].q;
935 GMX_LOG(logger.info).asParagraph().appendTextFormatted("Total mass%s %.3f a.m.u.", where, m);
936 GMX_LOG(logger.info).asParagraph().appendTextFormatted("Total charge%s %.3f e", where, qtot);
939 static void check_restp_type(const char* name, int t1, int t2)
943 gmx_fatal(FARGS, "Residues in one molecule have a different '%s' type: %d and %d", name, t1, t2);
947 static void check_restp_types(const PreprocessResidue& r0, const PreprocessResidue& r1)
949 check_restp_type("all dihedrals",
950 static_cast<int>(r0.bKeepAllGeneratedDihedrals),
951 static_cast<int>(r1.bKeepAllGeneratedDihedrals));
952 check_restp_type("nrexcl", r0.nrexcl, r1.nrexcl);
953 check_restp_type("HH14",
954 static_cast<int>(r0.bGenerateHH14Interactions),
955 static_cast<int>(r1.bGenerateHH14Interactions));
956 check_restp_type("remove dihedrals",
957 static_cast<int>(r0.bRemoveDihedralIfWithImproper),
958 static_cast<int>(r1.bRemoveDihedralIfWithImproper));
960 for (auto i : gmx::EnumerationWrapper<BondedTypes>{})
962 check_restp_type(enumValueToString(i), r0.rb[i].type, r1.rb[i].type);
966 static void add_atom_to_restp(PreprocessResidue* usedPpResidues,
969 const MoleculePatch* patch)
972 for (int k = 0; k < patch->nr; k++)
974 /* set counter in atomname */
975 std::string buf = patch->nname;
978 buf.append(gmx::formatString("%d", k + 1));
980 usedPpResidues->atomname.insert(usedPpResidues->atomname.begin() + at_start + 1 + k,
981 put_symtab(symtab, buf.c_str()));
982 usedPpResidues->atom.insert(usedPpResidues->atom.begin() + at_start + 1 + k, patch->atom.back());
983 if (patch->cgnr != NOTSET)
985 usedPpResidues->cgnr.insert(usedPpResidues->cgnr.begin() + at_start + 1 + k, patch->cgnr);
989 usedPpResidues->cgnr.insert(usedPpResidues->cgnr.begin() + at_start + 1 + k,
990 usedPpResidues->cgnr[at_start]);
995 void get_hackblocks_rtp(std::vector<MoleculePatchDatabase>* globalPatches,
996 std::vector<PreprocessResidue>* usedPpResidues,
997 gmx::ArrayRef<const PreprocessResidue> rtpFFDB,
1002 gmx::ArrayRef<MoleculePatchDatabase*> ntdb,
1003 gmx::ArrayRef<MoleculePatchDatabase*> ctdb,
1004 gmx::ArrayRef<const int> rn,
1005 gmx::ArrayRef<const int> rc,
1007 const gmx::MDLogger& logger)
1012 globalPatches->resize(nres);
1013 usedPpResidues->clear();
1014 /* first the termini */
1015 for (int i = 0; i < nterpairs; i++)
1017 if (rn[i] >= 0 && ntdb[i] != nullptr)
1019 copyModificationBlocks(*ntdb[i], &globalPatches->at(rn[i]));
1021 if (rc[i] >= 0 && ctdb[i] != nullptr)
1023 mergeAtomAndBondModifications(*ctdb[i], &globalPatches->at(rc[i]));
1027 /* then the whole rtp */
1028 for (int i = 0; i < nres; i++)
1030 /* Here we allow a mismatch of one character when looking for the rtp entry.
1031 * For such a mismatch there should be only one mismatching name.
1032 * This is mainly useful for small molecules such as ions.
1033 * Note that this will usually not work for protein, DNA and RNA,
1034 * since there the residue names should be listed in residuetypes.dat
1035 * and an error will have been generated earlier in the process.
1037 key = *resinfo[i].rtp;
1039 resinfo[i].rtp = put_symtab(symtab, searchResidueDatabase(key, rtpFFDB, logger).c_str());
1040 auto res = getDatabaseEntry(*resinfo[i].rtp, rtpFFDB);
1041 usedPpResidues->push_back(PreprocessResidue());
1042 PreprocessResidue* newentry = &usedPpResidues->back();
1043 copyPreprocessResidues(*res, newentry, symtab);
1045 /* Check that we do not have different bonded types in one molecule */
1046 check_restp_types(usedPpResidues->front(), *newentry);
1049 for (int j = 0; j < nterpairs && tern == -1; j++)
1057 for (int j = 0; j < nterpairs && terc == -1; j++)
1064 bRM = mergeBondedInteractionList(res->rb, globalPatches->at(i).rb, tern >= 0, terc >= 0);
1066 if (bRM && ((tern >= 0 && ntdb[tern] == nullptr) || (terc >= 0 && ctdb[terc] == nullptr)))
1068 const char* errString =
1069 "There is a dangling bond at at least one of the terminal ends and the force "
1070 "field does not provide terminal entries or files. Fix your terminal "
1071 "residues so that they match the residue database (.rtp) entries, or provide "
1072 "terminal database entries (.tdb).";
1075 GMX_LOG(logger.warning).asParagraph().appendTextFormatted("%s", errString);
1079 gmx_fatal(FARGS, "%s", errString);
1082 else if (bRM && ((tern >= 0 && ntdb[tern]->nhack() == 0) || (terc >= 0 && ctdb[terc]->nhack() == 0)))
1084 const char* errString =
1085 "There is a dangling bond at at least one of the terminal ends. Fix your "
1086 "coordinate file, add a new terminal database entry (.tdb), or select the "
1087 "proper existing terminal entry.";
1090 GMX_LOG(logger.warning).asParagraph().appendTextFormatted("%s", errString);
1094 gmx_fatal(FARGS, "%s", errString);
1099 /* Apply patchs to t_restp entries
1100 i.e. add's and deletes from termini database will be
1101 added to/removed from residue topology
1102 we'll do this on one big dirty loop, so it won't make easy reading! */
1103 for (auto modifiedResidue = globalPatches->begin(); modifiedResidue != globalPatches->end();
1106 const int pos = std::distance(globalPatches->begin(), modifiedResidue);
1107 PreprocessResidue* posres = &usedPpResidues->at(pos);
1108 for (auto patch = modifiedResidue->hack.begin(); patch != modifiedResidue->hack.end(); patch++)
1112 /* find atom in restp */
1113 auto found = std::find_if(
1114 posres->atomname.begin(), posres->atomname.end(), [&patch](char** name) {
1115 return (patch->oname.empty() && patch->a[0] == *name)
1116 || (patch->oname == *name);
1119 if (found == posres->atomname.end())
1121 /* If we are doing an atom rename only, we don't need
1122 * to generate a fatal error if the old name is not found
1125 /* Deleting can happen also only on the input atoms,
1126 * not necessarily always on the rtp entry.
1128 if (patch->type() == MoleculePatchType::Add)
1131 "atom %s not found in buiding block %d%s "
1132 "while combining tdb and rtp",
1133 patch->oname.empty() ? patch->a[0].c_str() : patch->oname.c_str(),
1140 int l = std::distance(posres->atomname.begin(), found);
1141 switch (patch->type())
1143 case MoleculePatchType::Add:
1146 add_atom_to_restp(posres, symtab, l, &(*patch));
1149 case MoleculePatchType::Delete:
1150 { /* we're deleting */
1151 posres->atom.erase(posres->atom.begin() + l);
1152 posres->atomname.erase(posres->atomname.begin() + l);
1153 posres->cgnr.erase(posres->cgnr.begin() + l);
1156 case MoleculePatchType::Replace:
1158 /* we're replacing */
1159 posres->atom[l] = patch->atom.back();
1160 posres->atomname[l] = put_symtab(symtab, patch->nname.c_str());
1161 if (patch->cgnr != NOTSET)
1163 posres->cgnr[l] = patch->cgnr;
1174 static bool atomname_cmp_nr(const char* anm, const MoleculePatch* patch, int* nr)
1181 return (gmx::equalCaseInsensitive(anm, patch->nname));
1185 if (isdigit(anm[strlen(anm) - 1]))
1187 *nr = anm[strlen(anm) - 1] - '0';
1193 if (*nr <= 0 || *nr > patch->nr)
1199 return (strlen(anm) == patch->nname.length() + 1
1200 && gmx_strncasecmp(anm, patch->nname.c_str(), patch->nname.length()) == 0);
1205 static bool match_atomnames_with_rtp_atom(t_atoms* pdba,
1206 gmx::ArrayRef<gmx::RVec> x,
1209 PreprocessResidue* localPpResidue,
1210 const MoleculePatchDatabase& singlePatch,
1212 const gmx::MDLogger& logger)
1219 oldnm = *pdba->atomname[atind];
1220 resnr = pdba->resinfo[pdba->atom[atind].resind].nr;
1223 for (auto patch = singlePatch.hack.begin(); patch != singlePatch.hack.end(); patch++)
1225 if (patch->type() == MoleculePatchType::Replace && gmx::equalCaseInsensitive(oldnm, patch->oname))
1227 /* This is a replace entry. */
1228 /* Check if we are not replacing a replaced atom. */
1229 bool bReplaceReplace = false;
1230 for (auto selfPatch = singlePatch.hack.begin(); selfPatch != singlePatch.hack.end(); selfPatch++)
1232 if (patch != selfPatch && selfPatch->type() == MoleculePatchType::Replace
1233 && gmx::equalCaseInsensitive(selfPatch->nname, patch->oname))
1235 /* The replace in patch replaces an atom that
1236 * was already replaced in selfPatch, we do not want
1237 * second or higher level replaces at this stage.
1239 bReplaceReplace = true;
1242 if (bReplaceReplace)
1244 /* Skip this replace. */
1248 /* This atom still has the old name, rename it */
1249 std::string newnm = patch->nname;
1250 auto found = std::find_if(
1251 localPpResidue->atomname.begin(),
1252 localPpResidue->atomname.end(),
1253 [&newnm](char** name) { return gmx::equalCaseInsensitive(newnm, *name); });
1254 if (found == localPpResidue->atomname.end())
1256 /* The new name is not present in the rtp.
1257 * We need to apply the replace also to the rtp entry.
1260 /* We need to find the add hack that can add this atom
1261 * to find out after which atom it should be added.
1263 bool bFoundInAdd = false;
1264 for (auto rtpModification = singlePatch.hack.begin();
1265 rtpModification != singlePatch.hack.end();
1268 int k = std::distance(localPpResidue->atomname.begin(), found);
1269 std::string start_at;
1270 if (rtpModification->type() == MoleculePatchType::Add
1271 && atomname_cmp_nr(newnm.c_str(), &(*rtpModification), &anmnr))
1275 start_at = singlePatch.hack[k].a[0];
1279 start_at = gmx::formatString(
1280 "%s%d", singlePatch.hack[k].nname.c_str(), anmnr - 1);
1282 auto found2 = std::find_if(localPpResidue->atomname.begin(),
1283 localPpResidue->atomname.end(),
1284 [&start_at](char** name) {
1285 return gmx::equalCaseInsensitive(start_at, *name);
1287 if (found2 == localPpResidue->atomname.end())
1290 "Could not find atom '%s' in residue building block '%s' to "
1293 localPpResidue->resname.c_str(),
1296 /* We can add the atom after atom start_nr */
1297 add_atom_to_restp(localPpResidue,
1299 std::distance(localPpResidue->atomname.begin(), found2),
1309 "Could not find an 'add' entry for atom named '%s' corresponding to "
1310 "the 'replace' entry from atom name '%s' to '%s' for tdb or hdb "
1311 "database of residue type '%s'",
1313 patch->oname.c_str(),
1314 patch->nname.c_str(),
1315 localPpResidue->resname.c_str());
1321 GMX_LOG(logger.info)
1323 .appendTextFormatted("Renaming atom '%s' in residue '%s' %d to '%s'",
1325 localPpResidue->resname.c_str(),
1329 /* Rename the atom in pdba */
1330 pdba->atomname[atind] = put_symtab(symtab, newnm.c_str());
1332 else if (patch->type() == MoleculePatchType::Delete
1333 && gmx::equalCaseInsensitive(oldnm, patch->oname))
1335 /* This is a delete entry, check if this atom is present
1336 * in the rtp entry of this residue.
1338 auto found3 = std::find_if(
1339 localPpResidue->atomname.begin(),
1340 localPpResidue->atomname.end(),
1341 [&oldnm](char** name) { return gmx::equalCaseInsensitive(oldnm, *name); });
1342 if (found3 == localPpResidue->atomname.end())
1344 /* This atom is not present in the rtp entry,
1345 * delete is now from the input pdba.
1349 GMX_LOG(logger.info)
1351 .appendTextFormatted("Deleting atom '%s' in residue '%s' %d",
1353 localPpResidue->resname.c_str(),
1356 /* We should free the atom name,
1357 * but it might be used multiple times in the symtab.
1358 * sfree(pdba->atomname[atind]);
1360 for (int k = atind + 1; k < pdba->nr; k++)
1362 pdba->atom[k - 1] = pdba->atom[k];
1363 pdba->atomname[k - 1] = pdba->atomname[k];
1364 copy_rvec(x[k], x[k - 1]);
1375 void match_atomnames_with_rtp(gmx::ArrayRef<PreprocessResidue> usedPpResidues,
1376 gmx::ArrayRef<MoleculePatchDatabase> globalPatches,
1379 gmx::ArrayRef<gmx::RVec> x,
1381 const gmx::MDLogger& logger)
1383 for (int i = 0; i < pdba->nr; i++)
1385 const char* oldnm = *pdba->atomname[i];
1386 PreprocessResidue* localPpResidue = &usedPpResidues[pdba->atom[i].resind];
1387 auto found = std::find_if(
1388 localPpResidue->atomname.begin(), localPpResidue->atomname.end(), [&oldnm](char** name) {
1389 return gmx::equalCaseInsensitive(oldnm, *name);
1391 if (found == localPpResidue->atomname.end())
1393 /* Not found yet, check if we have to rename this atom */
1394 if (match_atomnames_with_rtp_atom(
1395 pdba, x, symtab, i, localPpResidue, globalPatches[pdba->atom[i].resind], bVerbose, logger))
1397 /* We deleted this atom, decrease the atom counter by 1. */
1404 #define NUM_CMAP_ATOMS 5
1405 static void gen_cmap(InteractionsOfType* psb,
1406 gmx::ArrayRef<const PreprocessResidue> usedPpResidues,
1408 gmx::ArrayRef<const int> cyclicBondsIndex,
1409 const gmx::MDLogger& logger)
1413 t_resinfo* resinfo = atoms->resinfo;
1414 int nres = atoms->nres;
1415 int cmap_atomid[NUM_CMAP_ATOMS];
1416 int cmap_chainnum = -1;
1427 GMX_LOG(logger.info).asParagraph().appendTextFormatted("Making cmap torsions...");
1429 /* Most cmap entries use the N atom from the next residue, so the last
1430 * residue should not have its CMAP entry in that case, but for things like
1431 * dipeptides we sometimes define a complete CMAP entry inside a residue,
1432 * and in this case we need to process everything through the last residue.
1434 for (residx = 0; residx < nres; residx++)
1436 /* Add CMAP terms from the list of CMAP interactions */
1437 for (const auto& b : usedPpResidues[residx].rb[BondedTypes::Cmap].b)
1439 bool bAddCMAP = true;
1440 /* Loop over atoms in a candidate CMAP interaction and
1441 * check that they exist, are from the same chain and are
1442 * from residues labelled as protein. */
1443 for (int k = 0; k < NUM_CMAP_ATOMS && bAddCMAP; k++)
1445 /* Assign the pointer to the name of the next reference atom.
1446 * This can use -/+ labels to refer to previous/next residue.
1448 const char* pname = b.a[k].c_str();
1449 /* Skip this CMAP entry if it refers to residues before the
1450 * first or after the last residue.
1452 if ((cyclicBondsIndex.empty() && ((strchr(pname, '-') != nullptr) && (residx == 0)))
1453 || ((strchr(pname, '+') != nullptr) && (residx == nres - 1)))
1459 cmap_atomid[k] = search_atom(pname, i, atoms, ptr, TRUE, cyclicBondsIndex);
1460 bAddCMAP = bAddCMAP && (cmap_atomid[k] != -1);
1463 /* This break is necessary, because cmap_atomid[k]
1464 * == -1 cannot be safely used as an index
1465 * into the atom array. */
1468 int this_residue_index = atoms->atom[cmap_atomid[k]].resind;
1471 cmap_chainnum = resinfo[this_residue_index].chainnum;
1475 /* Does the residue for this atom have the same
1476 * chain number as the residues for previous
1478 bAddCMAP = bAddCMAP && cmap_chainnum == resinfo[this_residue_index].chainnum;
1480 /* Here we used to check that the residuetype was protein and
1481 * disable bAddCMAP if that was not the case. However, some
1482 * special residues (say, alanine dipeptides) might not adhere
1483 * to standard naming, and if we start calling them normal
1484 * protein residues the user will be bugged to select termini.
1486 * Instead, I believe that the right course of action is to
1487 * keep the CMAP interaction if it is present in the RTP file
1488 * and we correctly identified all atoms (which is the case
1505 if (residx < nres - 1)
1507 while (atoms->atom[i].resind < residx + 1)
1513 /* Start the next residue */
1516 static void scrub_charge_groups(int* cgnr, int natoms)
1520 for (i = 0; i < natoms; i++)
1527 void pdb2top(FILE* top_file,
1528 const char* posre_fn,
1529 const char* molname,
1531 std::vector<gmx::RVec>* x,
1532 PreprocessingAtomTypes* atype,
1534 gmx::ArrayRef<const PreprocessResidue> rtpFFDB,
1535 gmx::ArrayRef<PreprocessResidue> usedPpResidues,
1536 gmx::ArrayRef<MoleculePatchDatabase> globalPatches,
1539 bool bVsiteAromatics,
1542 gmx::ArrayRef<const DisulfideBond> ssbonds,
1543 real long_bond_dist,
1544 real short_bond_dist,
1550 gmx::ArrayRef<const int> cyclicBondsIndex,
1551 const gmx::MDLogger& logger)
1553 std::array<InteractionsOfType, F_NRE> plist;
1558 gmx::EnumerationArray<BondedTypes, int> bts;
1563 at2bonds(&(plist[F_BONDS]), globalPatches, atoms, *x, long_bond_dist, short_bond_dist, cyclicBondsIndex, logger);
1565 /* specbonds: disulphide bonds & heme-his */
1566 do_ssbonds(&(plist[F_BONDS]), atoms, ssbonds, bAllowMissing);
1568 nmissat = name2type(atoms, &cgnr, usedPpResidues, &rt, logger);
1573 GMX_LOG(logger.warning)
1575 .appendTextFormatted("There were %d missing atoms in molecule %s", nmissat, molname);
1580 "There were %d missing atoms in molecule %s, if you want to use this "
1581 "incomplete topology anyhow, use the option -missing",
1587 /* Cleanup bonds (sort and rm doubles) */
1588 clean_bonds(&(plist[F_BONDS]), logger);
1590 snew(vsite_type, atoms->nr);
1591 for (i = 0; i < atoms->nr; i++)
1593 vsite_type[i] = NOTSET;
1597 if (bVsiteAromatics)
1599 GMX_LOG(logger.info)
1601 .appendTextFormatted(
1602 "The conversion of aromatic rings into virtual sites is deprecated "
1603 "and may be removed in a future version of GROMACS");
1605 /* determine which atoms will be vsites and add dummy masses
1606 also renumber atom numbers in plist[0..F_NRE]! */
1607 do_vsites(rtpFFDB, atype, atoms, tab, x, plist, &vsite_type, &cgnr, mHmult, bVsiteAromatics, ffdir);
1610 /* Make Angles and Dihedrals */
1611 GMX_LOG(logger.info)
1613 .appendTextFormatted("Generating angles, dihedrals and pairs...");
1614 snew(excls, atoms->nr);
1615 gen_pad(atoms, usedPpResidues, plist, excls, globalPatches, bAllowMissing, cyclicBondsIndex);
1620 gen_cmap(&(plist[F_CMAP]), usedPpResidues, atoms, cyclicBondsIndex, logger);
1621 if (plist[F_CMAP].size() > 0)
1623 GMX_LOG(logger.info)
1625 .appendTextFormatted("There are %4zu cmap torsion pairs", plist[F_CMAP].size());
1629 /* set mass of all remaining hydrogen atoms */
1632 do_h_mass(&(plist[F_BONDS]), vsite_type, atoms, mHmult, bDeuterate);
1636 /* Cleanup bonds (sort and rm doubles) */
1637 /* clean_bonds(&(plist[F_BONDS]));*/
1639 GMX_LOG(logger.info)
1641 .appendTextFormatted(
1642 "There are %4zu dihedrals, %4zu impropers, %4zu angles\n"
1643 " %4zu pairs, %4zu bonds and %4zu virtual sites",
1644 plist[F_PDIHS].size(),
1645 plist[F_IDIHS].size(),
1646 plist[F_ANGLES].size(),
1647 plist[F_LJ14].size(),
1648 plist[F_BONDS].size(),
1649 plist[F_VSITE2].size() + plist[F_VSITE3].size() + plist[F_VSITE3FD].size()
1650 + plist[F_VSITE3FAD].size() + plist[F_VSITE3OUT].size()
1651 + plist[F_VSITE4FD].size() + plist[F_VSITE4FDN].size());
1653 print_sums(atoms, FALSE, logger);
1657 scrub_charge_groups(cgnr, atoms->nr);
1662 for (i = 0; i < atoms->nres; i++)
1664 atoms->resinfo[i].nr = i + 1;
1665 atoms->resinfo[i].ic = ' ';
1671 GMX_LOG(logger.info).asParagraph().appendTextFormatted("Writing topology");
1672 /* We can copy the bonded types from the first restp,
1673 * since the types have to be identical for all residues in one molecule.
1675 for (auto i : gmx::EnumerationWrapper<BondedTypes>{})
1677 bts[i] = usedPpResidues[0].rb[i].type;
1689 usedPpResidues[0].nrexcl);
1693 /* we should clean up hb and restp here, but that is a *L*O*T* of work! */
1695 for (i = 0; i < atoms->nr; i++)