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.
38 /* This file is completely threadsafe - keep it that way! */
51 #include "gromacs/fileio/confio.h"
52 #include "gromacs/gmxpreprocess/gpp_nextnb.h"
53 #include "gromacs/gmxpreprocess/grompp_impl.h"
54 #include "gromacs/gmxpreprocess/notset.h"
55 #include "gromacs/gmxpreprocess/pgutil.h"
56 #include "gromacs/gmxpreprocess/topio.h"
57 #include "gromacs/gmxpreprocess/toputil.h"
58 #include "gromacs/math/vec.h"
59 #include "gromacs/topology/ifunc.h"
60 #include "gromacs/utility/cstringutil.h"
61 #include "gromacs/utility/enumerationhelpers.h"
62 #include "gromacs/utility/fatalerror.h"
63 #include "gromacs/utility/smalloc.h"
65 #include "hackblock.h"
68 #define DIHEDRAL_WAS_SET_IN_RTP 0
69 static bool was_dihedral_set_in_rtp(const InteractionOfType& dih)
71 // This is a bad way to check this, but I don't know how to make this better now.
72 gmx::ArrayRef<const real> forceParam = dih.forceParam();
73 return forceParam[MAXFORCEPARAM - 1] == DIHEDRAL_WAS_SET_IN_RTP;
76 typedef bool (*peq)(const InteractionOfType& p1, const InteractionOfType& p2);
78 static bool acomp(const InteractionOfType& a1, const InteractionOfType& a2)
82 if (((ac = (a1.aj() - a2.aj())) != 0) || ((ac = (a1.ai() - a2.ai())) != 0))
88 return (a1.ak() < a2.ak());
92 static bool pcomp(const InteractionOfType& a1, const InteractionOfType& a2)
96 if ((pc = (a1.ai() - a2.ai())) != 0)
102 return (a1.aj() < a2.aj());
106 static bool dcomp(const InteractionOfType& d1, const InteractionOfType& d2)
110 /* First sort by J & K (the two central) atoms */
111 if (((dc = (d1.aj() - d2.aj())) != 0) || ((dc = (d1.ak() - d2.ak())) != 0))
112 { // NOLINT bugprone-branch-clone
115 /* Then make sure to put rtp dihedrals before generated ones */
116 else if (was_dihedral_set_in_rtp(d1) && !was_dihedral_set_in_rtp(d2))
120 else if (!was_dihedral_set_in_rtp(d1) && was_dihedral_set_in_rtp(d2))
124 /* Then sort by I and J (two outer) atoms */
125 else if (((dc = (d1.ai() - d2.ai())) != 0) || ((dc = (d1.al() - d2.al())) != 0))
131 // AMBER force fields with type 9 dihedrals can reach here, where we sort on
132 // the contents of the string that names the macro for the parameters.
133 return std::lexicographical_compare(d1.interactionTypeName().begin(),
134 d1.interactionTypeName().end(),
135 d2.interactionTypeName().begin(),
136 d2.interactionTypeName().end());
141 static bool is_dihedral_on_same_bond(const InteractionOfType& p1, const InteractionOfType& p2)
143 return ((p1.aj() == p2.aj()) && (p1.ak() == p2.ak()))
144 || ((p1.aj() == p2.ak()) && (p1.ak() == p2.aj()));
148 static bool preq(const InteractionOfType& p1, const InteractionOfType& p2)
150 return (p1.ai() == p2.ai()) && (p1.aj() == p2.aj());
153 static void rm2par(std::vector<InteractionOfType>* p, peq eq)
160 for (auto param = p->begin() + 1; param != p->end();)
162 auto prev = param - 1;
163 if (eq(*param, *prev))
165 param = p->erase(param);
174 static void cppar(gmx::ArrayRef<const InteractionOfType> types, gmx::ArrayRef<InteractionsOfType> plist, int ftype)
177 for (const auto& type : types)
179 plist[ftype].interactionTypes.push_back(type);
183 static bool idcomp(const InteractionOfType& a, const InteractionOfType& b)
187 if (((d = (a.ai() - b.ai())) != 0) || ((d = (a.al() - b.al())) != 0) || ((d = (a.aj() - b.aj())) != 0))
193 return (a.ak() < b.ak());
197 static void sort_id(gmx::ArrayRef<InteractionOfType> ps)
201 for (auto& parm : ps)
205 std::sort(ps.begin(), ps.end(), idcomp);
209 static int n_hydro(gmx::ArrayRef<const int> a, char*** atomname)
213 for (auto atom = a.begin(); atom < a.end(); atom += 3)
215 const char* aname = *atomname[*atom];
216 char c0 = toupper(aname[0]);
221 else if ((static_cast<int>(strlen(aname)) > 1) && (c0 >= '0') && (c0 <= '9'))
223 char c1 = toupper(aname[1]);
233 /* Clean up the dihedrals (both generated and read from the .rtp
235 static std::vector<InteractionOfType> clean_dih(gmx::ArrayRef<const InteractionOfType> dih,
236 gmx::ArrayRef<const InteractionOfType> improper,
238 bool bKeepAllGeneratedDihedrals,
239 bool bRemoveDihedralIfWithImproper)
241 /* Construct the list of the indices of the dihedrals
242 * (i.e. generated or read) that might be kept. */
243 std::vector<std::pair<InteractionOfType, int>> newDihedrals;
244 if (bKeepAllGeneratedDihedrals)
246 fprintf(stderr, "Keeping all generated dihedrals\n");
248 for (const auto& dihedral : dih)
250 newDihedrals.emplace_back(std::pair<InteractionOfType, int>(dihedral, i++));
255 /* Check if generated dihedral i should be removed. The
256 * dihedrals have been sorted by dcomp() above, so all those
257 * on the same two central atoms are together, with those from
258 * the .rtp file preceding those that were automatically
259 * generated. We remove the latter if the former exist. */
261 for (auto dihedral = dih.begin(); dihedral != dih.end(); dihedral++)
263 /* Keep the dihedrals that were defined in the .rtp file,
264 * and the dihedrals that were generated and different
265 * from the last one (whether it was generated or not). */
266 if (was_dihedral_set_in_rtp(*dihedral) || dihedral == dih.begin()
267 || !is_dihedral_on_same_bond(*dihedral, *(dihedral - 1)))
269 newDihedrals.emplace_back(std::pair<InteractionOfType, int>(*dihedral, i++));
274 for (auto dihedral = newDihedrals.begin(); dihedral != newDihedrals.end();)
276 bool bWasSetInRTP = was_dihedral_set_in_rtp(dihedral->first);
278 if (!bWasSetInRTP && bRemoveDihedralIfWithImproper)
280 /* Remove the dihedral if there is an improper on the same
282 for (auto imp = improper.begin(); imp != improper.end() && bKeep; ++imp)
284 bKeep = !is_dihedral_on_same_bond(dihedral->first, *imp);
290 /* If we don't want all dihedrals, we want to select the
291 * ones with the fewest hydrogens. Note that any generated
292 * dihedrals on the same bond as an .rtp dihedral may have
293 * been already pruned above in the construction of
294 * index[]. However, their parameters are still present,
295 * and l is looping over this dihedral and all of its
296 * pruned siblings. */
297 int bestl = dihedral->second;
298 if (!bKeepAllGeneratedDihedrals && !bWasSetInRTP)
300 /* Minimum number of hydrogens for i and l atoms */
302 int next = dihedral->second + 1;
303 for (int l = dihedral->second;
304 (l < next && is_dihedral_on_same_bond(dihedral->first, dih[l]));
307 int nh = n_hydro(dih[l].atoms(), atoms->atomname);
318 dihedral->first = dih[bestl];
328 dihedral = newDihedrals.erase(dihedral);
331 std::vector<InteractionOfType> finalDihedrals;
332 finalDihedrals.reserve(newDihedrals.size());
333 for (const auto& param : newDihedrals)
335 finalDihedrals.emplace_back(param.first);
337 return finalDihedrals;
340 static std::vector<InteractionOfType> get_impropers(t_atoms* atoms,
341 gmx::ArrayRef<MoleculePatchDatabase> globalPatches,
343 gmx::ArrayRef<const int> cyclicBondsIndex)
345 std::vector<InteractionOfType> improper;
347 /* Add all the impropers from the residue database to the list. */
349 if (!globalPatches.empty())
351 for (int i = 0; (i < atoms->nres); i++)
353 BondedInteractionList* impropers = &globalPatches[i].rb[BondedTypes::ImproperDihedrals];
354 for (const auto& bondeds : impropers->b)
358 for (int k = 0; (k < 4) && !bStop; k++)
360 const int entry = search_atom(
361 bondeds.a[k].c_str(), start, atoms, "improper", bAllowMissing, cyclicBondsIndex);
365 ai.emplace_back(entry);
375 improper.emplace_back(InteractionOfType(ai, {}, bondeds.s));
378 while ((start < atoms->nr) && (atoms->atom[start].resind == i))
388 static int nb_dist(t_nextnb* nnb, int ai, int aj)
400 nrexcl = nnb->nrexcl[ai];
401 for (nre = 1; (nre < nnb->nrex); nre++)
404 for (nrx = 0; (nrx < nrexcl[nre]); nrx++)
406 if ((aj == a[nrx]) && (NRE == -1))
415 static bool is_hydro(t_atoms* atoms, int ai)
417 return ((*(atoms->atomname[ai]))[0] == 'H');
420 static void get_atomnames_min(int n,
421 gmx::ArrayRef<std::string> anm,
424 gmx::ArrayRef<const int> a)
426 /* Assume ascending residue numbering */
427 for (int m = 0; m < n; m++)
429 if (atoms->atom[a[m]].resind < resind)
433 else if (atoms->atom[a[m]].resind > resind)
441 anm[m].append(*(atoms->atomname[a[m]]));
445 static void gen_excls(t_atoms* atoms,
447 gmx::ArrayRef<MoleculePatchDatabase> globalPatches,
449 gmx::ArrayRef<const int> cyclicBondsIndex)
452 for (int a = 0; a < atoms->nr; a++)
454 int r = atoms->atom[a].resind;
455 if (a == atoms->nr - 1 || atoms->atom[a + 1].resind != r)
457 BondedInteractionList* hbexcl = &globalPatches[r].rb[BondedTypes::Exclusions];
459 for (const auto& bondeds : hbexcl->b)
461 const char* anm = bondeds.a[0].c_str();
462 int i1 = search_atom(anm, astart, atoms, "exclusion", bAllowMissing, cyclicBondsIndex);
463 anm = bondeds.a[1].c_str();
464 int i2 = search_atom(anm, astart, atoms, "exclusion", bAllowMissing, cyclicBondsIndex);
465 if (i1 != -1 && i2 != -1)
473 srenew(excls[i1].e, excls[i1].nr + 1);
474 excls[i1].e[excls[i1].nr] = i2;
483 for (int a = 0; a < atoms->nr; a++)
487 std::sort(excls[a].e, excls[a].e + excls[a].nr);
492 static void remove_excl(t_excls* excls, int remove)
496 for (i = remove + 1; i < excls->nr; i++)
498 excls->e[i - 1] = excls->e[i];
504 static void clean_excls(t_nextnb* nnb, int nrexcl, t_excls excls[])
506 int i, j, j1, k, k1, l, l1, e;
511 /* extract all i-j-k-l neighbours from nnb struct */
512 for (i = 0; (i < nnb->nr); i++)
514 /* For all particles */
517 for (j = 0; (j < nnb->nrexcl[i][1]); j++)
519 /* For all first neighbours */
520 j1 = nnb->a[i][1][j];
522 for (e = 0; e < excl->nr; e++)
524 if (excl->e[e] == j1)
526 remove_excl(excl, e);
532 for (k = 0; (k < nnb->nrexcl[j1][1]); k++)
534 /* For all first neighbours of j1 */
535 k1 = nnb->a[j1][1][k];
537 for (e = 0; e < excl->nr; e++)
539 if (excl->e[e] == k1)
541 remove_excl(excl, e);
547 for (l = 0; (l < nnb->nrexcl[k1][1]); l++)
549 /* For all first neighbours of k1 */
550 l1 = nnb->a[k1][1][l];
552 for (e = 0; e < excl->nr; e++)
554 if (excl->e[e] == l1)
556 remove_excl(excl, e);
569 * Generate pairs, angles and dihedrals from .rtp settings
571 * \param[in,out] atoms Global information about atoms in topology.
572 * \param[in] rtpFFDB Residue type database from force field.
573 * \param[in,out] plist Information about listed interactions.
574 * \param[in,out] excls Pair interaction exclusions.
575 * \param[in,out] globalPatches Information about possible residue modifications.
576 * \param[in] bAllowMissing True if missing interaction information is allowed.
577 * AKA allow cartoon physics
578 * \param[in] cyclicBondsIndex Information about bonds creating cyclic molecules.
579 * Empty if no such bonds exist.
581 void gen_pad(t_atoms* atoms,
582 gmx::ArrayRef<const PreprocessResidue> rtpFFDB,
583 gmx::ArrayRef<InteractionsOfType> plist,
585 gmx::ArrayRef<MoleculePatchDatabase> globalPatches,
587 gmx::ArrayRef<const int> cyclicBondsIndex)
590 init_nnb(&nnb, atoms->nr, 4);
591 gen_nnb(&nnb, plist);
592 print_nnb(&nnb, "NNB");
594 /* These are the angles, dihedrals and pairs that we generate
595 * from the bonds. The ones that are already there from the rtp file
598 std::vector<InteractionOfType> ang;
599 std::vector<InteractionOfType> dih;
600 std::vector<InteractionOfType> pai;
601 std::vector<InteractionOfType> improper;
603 std::array<std::string, 4> anm;
605 if (!globalPatches.empty())
607 gen_excls(atoms, excls, globalPatches, bAllowMissing, cyclicBondsIndex);
608 /* mark all entries as not matched yet */
609 for (int i = 0; i < atoms->nres; i++)
611 for (auto bondedsList : globalPatches[i].rb)
613 for (auto& bondeds : bondedsList.b)
615 bondeds.match = false;
621 /* Extract all i-j-k-l neighbours from nnb struct to generate all
622 * angles and dihedrals. */
623 for (int i = 0; (i < nnb.nr); i++)
625 /* For all particles */
626 for (int j = 0; (j < nnb.nrexcl[i][1]); j++)
628 /* For all first neighbours */
629 int j1 = nnb.a[i][1][j];
630 for (int k = 0; (k < nnb.nrexcl[j1][1]); k++)
632 /* For all first neighbours of j1 */
633 int k1 = nnb.a[j1][1][k];
636 /* Generate every angle only once */
639 std::vector<int> atomNumbers = { i, j1, k1 };
641 if (!globalPatches.empty())
643 int minres = atoms->atom[i].resind;
645 for (int m = 1; m < 3; m++)
647 minres = std::min(minres, atoms->atom[atomNumbers[m]].resind);
648 maxres = std::max(maxres, atoms->atom[atomNumbers[m]].resind);
650 int res = 2 * minres - maxres;
653 res += maxres - minres;
654 get_atomnames_min(3, anm, res, atoms, atomNumbers);
655 BondedInteractionList* hbang = &globalPatches[res].rb[BondedTypes::Angles];
656 for (auto& bondeds : hbang->b)
658 if (anm[1] == bondeds.aj())
661 for (int m = 0; m < 3; m += 2)
664 || ((anm[m] == bondeds.ai())
665 && (anm[2 - m] == bondeds.ak())));
670 /* Mark that we found a match for this entry */
671 bondeds.match = true;
675 } while (res < maxres);
677 ang.push_back(InteractionOfType(atomNumbers, {}, name));
679 /* Generate every dihedral, 1-4 exclusion and 1-4 interaction
683 for (int l = 0; (l < nnb.nrexcl[k1][1]); l++)
685 /* For all first neighbours of k1 */
686 int l1 = nnb.a[k1][1][l];
687 if ((l1 != i) && (l1 != j1))
689 std::vector<int> atomNumbers = { i, j1, k1, l1 };
692 if (!globalPatches.empty())
694 int minres = atoms->atom[i].resind;
696 for (int m = 1; m < 4; m++)
698 minres = std::min(minres, atoms->atom[atomNumbers[m]].resind);
699 maxres = std::max(maxres, atoms->atom[atomNumbers[m]].resind);
701 int res = 2 * minres - maxres;
704 res += maxres - minres;
705 get_atomnames_min(4, anm, res, atoms, atomNumbers);
706 BondedInteractionList* hbdih =
707 &globalPatches[res].rb[BondedTypes::ProperDihedrals];
708 for (auto& bondeds : hbdih->b)
711 for (int m = 0; m < 2; m++)
714 || ((anm[3 * m] == bondeds.ai())
715 && (anm[1 + m] == bondeds.aj())
716 && (anm[2 - m] == bondeds.ak())
717 && (anm[3 - 3 * m] == bondeds.al())));
722 /* Mark that we found a match for this entry */
723 bondeds.match = true;
725 /* Set the last parameter to be able to see
726 if the dihedral was in the rtp list.
729 dih.push_back(InteractionOfType(atomNumbers, {}, name));
730 dih.back().setForceParameter(
731 MAXFORCEPARAM - 1, DIHEDRAL_WAS_SET_IN_RTP);
734 } while (res < maxres);
738 std::vector<int> atoms = { i, j1, k1, l1 };
739 dih.push_back(InteractionOfType(atoms, {}, ""));
742 int nbd = nb_dist(&nnb, i, l1);
745 int i1 = std::min(i, l1);
746 int i2 = std::max(i, l1);
748 for (int m = 0; m < excls[i1].nr; m++)
750 bExcl = bExcl || excls[i1].e[m] == i2;
754 if (rtpFFDB[0].bGenerateHH14Interactions
755 || !(is_hydro(atoms, i1) && is_hydro(atoms, i2)))
757 std::vector<int> atoms = { i1, i2 };
758 pai.push_back(InteractionOfType(atoms, {}, ""));
770 if (!globalPatches.empty())
772 /* The above approach is great in that we double-check that e.g. an angle
773 * really corresponds to three atoms connected by bonds, but this is not
774 * generally true. Go through the angle and dihedral hackblocks to add
775 * entries that we have not yet marked as matched when going through bonds.
777 for (int i = 0; i < atoms->nres; i++)
779 /* Add remaining angles from hackblock */
780 BondedInteractionList* hbang = &globalPatches[i].rb[BondedTypes::Angles];
781 for (auto& bondeds : hbang->b)
785 /* We already used this entry, continue to the next */
788 /* Hm - entry not used, let's see if we can find all atoms */
789 std::vector<int> atomNumbers;
791 gmx::ArrayRef<const int>::iterator cyclicBondsIterator;
792 for (int k = 0; k < 3 && bFound; k++)
794 const char* p = bondeds.a[k].c_str();
799 cyclicBondsIterator =
800 std::find(cyclicBondsIndex.begin(), cyclicBondsIndex.end(), res--);
801 if (cyclicBondsIterator != cyclicBondsIndex.end()
802 && !((cyclicBondsIterator - cyclicBondsIndex.begin()) & 1))
804 res = *(++cyclicBondsIterator);
807 else if (p[0] == '+')
810 cyclicBondsIterator =
811 std::find(cyclicBondsIndex.begin(), cyclicBondsIndex.end(), res++);
812 if (cyclicBondsIterator != cyclicBondsIndex.end()
813 && ((cyclicBondsIterator - cyclicBondsIndex.begin()) & 1))
815 res = *(--cyclicBondsIterator);
818 atomNumbers.emplace_back(search_res_atom(p, res, atoms, "angle", TRUE));
819 bFound = (atomNumbers.back() != -1);
824 bondeds.match = true;
825 /* Incrementing nang means we save this angle */
826 ang.push_back(InteractionOfType(atomNumbers, {}, bondeds.s));
830 /* Add remaining dihedrals from hackblock */
831 BondedInteractionList* hbdih = &globalPatches[i].rb[BondedTypes::ProperDihedrals];
832 for (auto& bondeds : hbdih->b)
836 /* We already used this entry, continue to the next */
839 /* Hm - entry not used, let's see if we can find all atoms */
840 std::vector<int> atomNumbers;
842 gmx::ArrayRef<const int>::iterator cyclicBondsIterator;
843 for (int k = 0; k < 4 && bFound; k++)
845 const char* p = bondeds.a[k].c_str();
850 cyclicBondsIterator =
851 std::find(cyclicBondsIndex.begin(), cyclicBondsIndex.end(), res);
853 if (cyclicBondsIterator != cyclicBondsIndex.end()
854 && !((cyclicBondsIterator - cyclicBondsIndex.begin()) & 1))
856 res = *(++cyclicBondsIterator);
859 else if (p[0] == '+')
862 cyclicBondsIterator =
863 std::find(cyclicBondsIndex.begin(), cyclicBondsIndex.end(), res);
865 if (cyclicBondsIterator != cyclicBondsIndex.end()
866 && ((cyclicBondsIterator - cyclicBondsIndex.begin()) & 1))
868 res = *(--cyclicBondsIterator);
871 atomNumbers.emplace_back(search_res_atom(p, res, atoms, "dihedral", TRUE));
872 bFound = (atomNumbers.back() != -1);
877 bondeds.match = true;
878 /* Incrementing ndih means we save this dihedral */
879 dih.push_back(InteractionOfType(atomNumbers, {}, bondeds.s));
885 /* Sort angles with respect to j-i-k (middle atom first) */
888 std::sort(ang.begin(), ang.end(), acomp);
891 /* Sort dihedrals with respect to j-k-i-l (middle atoms first) */
894 std::sort(dih.begin(), dih.end(), dcomp);
900 std::sort(pai.begin(), pai.end(), pcomp);
904 /* Remove doubles, could occur in 6-rings, such as phenyls,
905 maybe one does not want this when fudgeQQ < 1.
907 fprintf(stderr, "Before cleaning: %zu pairs\n", pai.size());
911 /* Get the impropers from the database */
912 improper = get_impropers(atoms, globalPatches, bAllowMissing, cyclicBondsIndex);
914 /* Sort the impropers */
919 fprintf(stderr, "Before cleaning: %zu dihedrals\n", dih.size());
920 dih = clean_dih(dih, improper, atoms, rtpFFDB[0].bKeepAllGeneratedDihedrals, rtpFFDB[0].bRemoveDihedralIfWithImproper);
923 /* Now we have unique lists of angles and dihedrals
924 * Copy them into the destination struct
926 cppar(ang, plist, F_ANGLES);
927 cppar(dih, plist, F_PDIHS);
928 cppar(improper, plist, F_IDIHS);
929 cppar(pai, plist, F_LJ14);
931 /* Remove all exclusions which are within nrexcl */
932 clean_excls(&nnb, rtpFFDB[0].nrexcl, excls);