src/gromacs/gmxpreprocess/gen_ad.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,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.
  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 /* This file is completely threadsafe - keep it that way! */
  39 #include "gmxpre.h"
  40
  41 #include "gen_ad.h"
  42
  43 #include <cctype>
  44 #include <cmath>
  45 #include <cstdlib>
  46 #include <cstring>
  47
  48 #include <algorithm>
  49 #include <numeric>
  50
  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/fatalerror.h"
  62 #include "gromacs/utility/smalloc.h"
  63
  64 #include "hackblock.h"
  65 #include "resall.h"
  66
  67 #define DIHEDRAL_WAS_SET_IN_RTP 0
  68 static bool was_dihedral_set_in_rtp(const InteractionOfType& dih)
  69 {
  70     // This is a bad way to check this, but I don't know how to make this better now.
  71     gmx::ArrayRef<const real> forceParam = dih.forceParam();
  72     return forceParam[MAXFORCEPARAM - 1] == DIHEDRAL_WAS_SET_IN_RTP;
  73 }
  74
  75 typedef bool (*peq)(const InteractionOfType& p1, const InteractionOfType& p2);
  76
  77 static bool acomp(const InteractionOfType& a1, const InteractionOfType& a2)
  78 {
  79     int ac;
  80
  81     if (((ac = (a1.aj() - a2.aj())) != 0) || ((ac = (a1.ai() - a2.ai())) != 0))
  82     {
  83         return ac < 0;
  84     }
  85     else
  86     {
  87         return (a1.ak() < a2.ak());
  88     }
  89 }
  90
  91 static bool pcomp(const InteractionOfType& a1, const InteractionOfType& a2)
  92 {
  93     int pc;
  94
  95     if ((pc = (a1.ai() - a2.ai())) != 0)
  96     {
  97         return pc < 0;
  98     }
  99     else
 100     {
 101         return (a1.aj() < a2.aj());
 102     }
 103 }
 104
 105 static bool dcomp(const InteractionOfType& d1, const InteractionOfType& d2)
 106 {
 107     int dc;
 108
 109     /* First sort by J & K (the two central) atoms */
 110     if (((dc = (d1.aj() - d2.aj())) != 0) || ((dc = (d1.ak() - d2.ak())) != 0))
 111     { // NOLINT bugprone-branch-clone
 112         return dc < 0;
 113     }
 114     /* Then make sure to put rtp dihedrals before generated ones */
 115     else if (was_dihedral_set_in_rtp(d1) && !was_dihedral_set_in_rtp(d2))
 116     {
 117         return true;
 118     }
 119     else if (!was_dihedral_set_in_rtp(d1) && was_dihedral_set_in_rtp(d2))
 120     {
 121         return false;
 122     }
 123     /* Then sort by I and J (two outer) atoms */
 124     else if (((dc = (d1.ai() - d2.ai())) != 0) || ((dc = (d1.al() - d2.al())) != 0))
 125     {
 126         return dc < 0;
 127     }
 128     else
 129     {
 130         // AMBER force fields with type 9 dihedrals can reach here, where we sort on
 131         // the contents of the string that names the macro for the parameters.
 132         return std::lexicographical_compare(d1.interactionTypeName().begin(),
 133                                             d1.interactionTypeName().end(),
 134                                             d2.interactionTypeName().begin(),
 135                                             d2.interactionTypeName().end());
 136     }
 137 }
 138
 139
 140 static bool is_dihedral_on_same_bond(const InteractionOfType& p1, const InteractionOfType& p2)
 141 {
 142     return ((p1.aj() == p2.aj()) && (p1.ak() == p2.ak()))
 143            || ((p1.aj() == p2.ak()) && (p1.ak() == p2.aj()));
 144 }
 145
 146
 147 static bool preq(const InteractionOfType& p1, const InteractionOfType& p2)
 148 {
 149     return (p1.ai() == p2.ai()) && (p1.aj() == p2.aj());
 150 }
 151
 152 static void rm2par(std::vector<InteractionOfType>* p, peq eq)
 153 {
 154     if (p->empty())
 155     {
 156         return;
 157     }
 158
 159     for (auto param = p->begin() + 1; param != p->end();)
 160     {
 161         auto prev = param - 1;
 162         if (eq(*param, *prev))
 163         {
 164             param = p->erase(param);
 165         }
 166         else
 167         {
 168             ++param;
 169         }
 170     }
 171 }
 172
 173 static void cppar(gmx::ArrayRef<const InteractionOfType> types, gmx::ArrayRef<InteractionsOfType> plist, int ftype)
 174 {
 175     /* Keep old stuff */
 176     for (const auto& type : types)
 177     {
 178         plist[ftype].interactionTypes.push_back(type);
 179     }
 180 }
 181
 182 static bool idcomp(const InteractionOfType& a, const InteractionOfType& b)
 183 {
 184     int d;
 185
 186     if (((d = (a.ai() - b.ai())) != 0) || ((d = (a.al() - b.al())) != 0) || ((d = (a.aj() - b.aj())) != 0))
 187     {
 188         return d < 0;
 189     }
 190     else
 191     {
 192         return (a.ak() < b.ak());
 193     }
 194 }
 195
 196 static void sort_id(gmx::ArrayRef<InteractionOfType> ps)
 197 {
 198     if (ps.size() > 1)
 199     {
 200         for (auto& parm : ps)
 201         {
 202             parm.sortAtomIds();
 203         }
 204         std::sort(ps.begin(), ps.end(), idcomp);
 205     }
 206 }
 207
 208 static int n_hydro(gmx::ArrayRef<const int> a, char*** atomname)
 209 {
 210     int nh = 0;
 211
 212     for (auto atom = a.begin(); atom < a.end(); atom += 3)
 213     {
 214         const char* aname = *atomname[*atom];
 215         char        c0    = toupper(aname[0]);
 216         if (c0 == 'H')
 217         {
 218             nh++;
 219         }
 220         else if ((static_cast<int>(strlen(aname)) > 1) && (c0 >= '0') && (c0 <= '9'))
 221         {
 222             char c1 = toupper(aname[1]);
 223             if (c1 == 'H')
 224             {
 225                 nh++;
 226             }
 227         }
 228     }
 229     return nh;
 230 }
 231
 232 /* Clean up the dihedrals (both generated and read from the .rtp
 233  * file). */
 234 static std::vector<InteractionOfType> clean_dih(gmx::ArrayRef<const InteractionOfType> dih,
 235                                                 gmx::ArrayRef<const InteractionOfType> improper,
 236                                                 t_atoms*                               atoms,
 237                                                 bool bKeepAllGeneratedDihedrals,
 238                                                 bool bRemoveDihedralIfWithImproper)
 239 {
 240     /* Construct the list of the indices of the dihedrals
 241      * (i.e. generated or read) that might be kept. */
 242     std::vector<std::pair<InteractionOfType, int>> newDihedrals;
 243     if (bKeepAllGeneratedDihedrals)
 244     {
 245         fprintf(stderr, "Keeping all generated dihedrals\n");
 246         int i = 0;
 247         for (const auto& dihedral : dih)
 248         {
 249             newDihedrals.emplace_back(std::pair<InteractionOfType, int>(dihedral, i++));
 250         }
 251     }
 252     else
 253     {
 254         /* Check if generated dihedral i should be removed. The
 255          * dihedrals have been sorted by dcomp() above, so all those
 256          * on the same two central atoms are together, with those from
 257          * the .rtp file preceding those that were automatically
 258          * generated. We remove the latter if the former exist. */
 259         int i = 0;
 260         for (auto dihedral = dih.begin(); dihedral != dih.end(); dihedral++)
 261         {
 262             /* Keep the dihedrals that were defined in the .rtp file,
 263              * and the dihedrals that were generated and different
 264              * from the last one (whether it was generated or not). */
 265             if (was_dihedral_set_in_rtp(*dihedral) || dihedral == dih.begin()
 266                 || !is_dihedral_on_same_bond(*dihedral, *(dihedral - 1)))
 267             {
 268                 newDihedrals.emplace_back(std::pair<InteractionOfType, int>(*dihedral, i++));
 269             }
 270         }
 271     }
 272     int k = 0;
 273     for (auto dihedral = newDihedrals.begin(); dihedral != newDihedrals.end();)
 274     {
 275         bool bWasSetInRTP = was_dihedral_set_in_rtp(dihedral->first);
 276         bool bKeep        = true;
 277         if (!bWasSetInRTP && bRemoveDihedralIfWithImproper)
 278         {
 279             /* Remove the dihedral if there is an improper on the same
 280              * bond. */
 281             for (auto imp = improper.begin(); imp != improper.end() && bKeep; ++imp)
 282             {
 283                 bKeep = !is_dihedral_on_same_bond(dihedral->first, *imp);
 284             }
 285         }
 286
 287         if (bKeep)
 288         {
 289             /* If we don't want all dihedrals, we want to select the
 290              * ones with the fewest hydrogens. Note that any generated
 291              * dihedrals on the same bond as an .rtp dihedral may have
 292              * been already pruned above in the construction of
 293              * index[]. However, their parameters are still present,
 294              * and l is looping over this dihedral and all of its
 295              * pruned siblings. */
 296             int bestl = dihedral->second;
 297             if (!bKeepAllGeneratedDihedrals && !bWasSetInRTP)
 298             {
 299                 /* Minimum number of hydrogens for i and l atoms */
 300                 int minh = 2;
 301                 int next = dihedral->second + 1;
 302                 for (int l = dihedral->second;
 303                      (l < next && is_dihedral_on_same_bond(dihedral->first, dih[l]));
 304                      l++)
 305                 {
 306                     int nh = n_hydro(dih[l].atoms(), atoms->atomname);
 307                     if (nh < minh)
 308                     {
 309                         minh  = nh;
 310                         bestl = l;
 311                     }
 312                     if (0 == minh)
 313                     {
 314                         break;
 315                     }
 316                 }
 317                 dihedral->first = dih[bestl];
 318             }
 319             if (k == bestl)
 320             {
 321                 ++dihedral;
 322             }
 323             k++;
 324         }
 325         else
 326         {
 327             dihedral = newDihedrals.erase(dihedral);
 328         }
 329     }
 330     std::vector<InteractionOfType> finalDihedrals;
 331     finalDihedrals.reserve(newDihedrals.size());
 332     for (const auto& param : newDihedrals)
 333     {
 334         finalDihedrals.emplace_back(param.first);
 335     }
 336     return finalDihedrals;
 337 }
 338
 339 static std::vector<InteractionOfType> get_impropers(t_atoms*                             atoms,
 340                                                     gmx::ArrayRef<MoleculePatchDatabase> globalPatches,
 341                                                     bool                     bAllowMissing,
 342                                                     gmx::ArrayRef<const int> cyclicBondsIndex)
 343 {
 344     std::vector<InteractionOfType> improper;
 345
 346     /* Add all the impropers from the residue database to the list. */
 347     int start = 0;
 348     if (!globalPatches.empty())
 349     {
 350         for (int i = 0; (i < atoms->nres); i++)
 351         {
 352             BondedInteractionList* impropers = &globalPatches[i].rb[BondedTypes::ImproperDihedrals];
 353             for (const auto& bondeds : impropers->b)
 354             {
 355                 bool             bStop = false;
 356                 std::vector<int> ai;
 357                 for (int k = 0; (k < 4) && !bStop; k++)
 358                 {
 359                     const int entry = search_atom(
 360                             bondeds.a[k].c_str(), start, atoms, "improper", bAllowMissing, cyclicBondsIndex);
 361
 362                     if (entry != -1)
 363                     {
 364                         ai.emplace_back(entry);
 365                     }
 366                     else
 367                     {
 368                         bStop = true;
 369                     }
 370                 }
 371                 if (!bStop)
 372                 {
 373                     /* Not broken out */
 374                     improper.emplace_back(InteractionOfType(ai, {}, bondeds.s));
 375                 }
 376             }
 377             while ((start < atoms->nr) && (atoms->atom[start].resind == i))
 378             {
 379                 start++;
 380             }
 381         }
 382     }
 383
 384     return improper;
 385 }
 386
 387 static int nb_dist(t_nextnb* nnb, int ai, int aj)
 388 {
 389     int  nre, nrx, NRE;
 390     int* nrexcl;
 391     int* a;
 392
 393     if (ai == aj)
 394     {
 395         return 0;
 396     }
 397
 398     NRE    = -1;
 399     nrexcl = nnb->nrexcl[ai];
 400     for (nre = 1; (nre < nnb->nrex); nre++)
 401     {
 402         a = nnb->a[ai][nre];
 403         for (nrx = 0; (nrx < nrexcl[nre]); nrx++)
 404         {
 405             if ((aj == a[nrx]) && (NRE == -1))
 406             {
 407                 NRE = nre;
 408             }
 409         }
 410     }
 411     return NRE;
 412 }
 413
 414 static bool is_hydro(t_atoms* atoms, int ai)
 415 {
 416     return ((*(atoms->atomname[ai]))[0] == 'H');
 417 }
 418
 419 static void get_atomnames_min(int                        n,
 420                               gmx::ArrayRef<std::string> anm,
 421                               int                        resind,
 422                               t_atoms*                   atoms,
 423                               gmx::ArrayRef<const int>   a)
 424 {
 425     /* Assume ascending residue numbering */
 426     for (int m = 0; m < n; m++)
 427     {
 428         if (atoms->atom[a[m]].resind < resind)
 429         {
 430             anm[m] = "-";
 431         }
 432         else if (atoms->atom[a[m]].resind > resind)
 433         {
 434             anm[m] = "+";
 435         }
 436         else
 437         {
 438             anm[m] = "";
 439         }
 440         anm[m].append(*(atoms->atomname[a[m]]));
 441     }
 442 }
 443
 444 static void gen_excls(t_atoms*                             atoms,
 445                       t_excls*                             excls,
 446                       gmx::ArrayRef<MoleculePatchDatabase> globalPatches,
 447                       bool                                 bAllowMissing,
 448                       gmx::ArrayRef<const int>             cyclicBondsIndex)
 449 {
 450     int astart = 0;
 451     for (int a = 0; a < atoms->nr; a++)
 452     {
 453         int r = atoms->atom[a].resind;
 454         if (a == atoms->nr - 1 || atoms->atom[a + 1].resind != r)
 455         {
 456             BondedInteractionList* hbexcl = &globalPatches[r].rb[BondedTypes::Exclusions];
 457
 458             for (const auto& bondeds : hbexcl->b)
 459             {
 460                 const char* anm = bondeds.a[0].c_str();
 461                 int i1 = search_atom(anm, astart, atoms, "exclusion", bAllowMissing, cyclicBondsIndex);
 462                 anm    = bondeds.a[1].c_str();
 463                 int i2 = search_atom(anm, astart, atoms, "exclusion", bAllowMissing, cyclicBondsIndex);
 464                 if (i1 != -1 && i2 != -1)
 465                 {
 466                     if (i1 > i2)
 467                     {
 468                         int itmp = i1;
 469                         i1       = i2;
 470                         i2       = itmp;
 471                     }
 472                     srenew(excls[i1].e, excls[i1].nr + 1);
 473                     excls[i1].e[excls[i1].nr] = i2;
 474                     excls[i1].nr++;
 475                 }
 476             }
 477
 478             astart = a + 1;
 479         }
 480     }
 481
 482     for (int a = 0; a < atoms->nr; a++)
 483     {
 484         if (excls[a].nr > 1)
 485         {
 486             std::sort(excls[a].e, excls[a].e + excls[a].nr);
 487         }
 488     }
 489 }
 490
 491 static void remove_excl(t_excls* excls, int remove)
 492 {
 493     int i;
 494
 495     for (i = remove + 1; i < excls->nr; i++)
 496     {
 497         excls->e[i - 1] = excls->e[i];
 498     }
 499
 500     excls->nr--;
 501 }
 502
 503 static void clean_excls(t_nextnb* nnb, int nrexcl, t_excls excls[])
 504 {
 505     int      i, j, j1, k, k1, l, l1, e;
 506     t_excls* excl;
 507
 508     if (nrexcl >= 1)
 509     {
 510         /* extract all i-j-k-l neighbours from nnb struct */
 511         for (i = 0; (i < nnb->nr); i++)
 512         {
 513             /* For all particles */
 514             excl = &excls[i];
 515
 516             for (j = 0; (j < nnb->nrexcl[i][1]); j++)
 517             {
 518                 /* For all first neighbours */
 519                 j1 = nnb->a[i][1][j];
 520
 521                 for (e = 0; e < excl->nr; e++)
 522                 {
 523                     if (excl->e[e] == j1)
 524                     {
 525                         remove_excl(excl, e);
 526                     }
 527                 }
 528
 529                 if (nrexcl >= 2)
 530                 {
 531                     for (k = 0; (k < nnb->nrexcl[j1][1]); k++)
 532                     {
 533                         /* For all first neighbours of j1 */
 534                         k1 = nnb->a[j1][1][k];
 535
 536                         for (e = 0; e < excl->nr; e++)
 537                         {
 538                             if (excl->e[e] == k1)
 539                             {
 540                                 remove_excl(excl, e);
 541                             }
 542                         }
 543
 544                         if (nrexcl >= 3)
 545                         {
 546                             for (l = 0; (l < nnb->nrexcl[k1][1]); l++)
 547                             {
 548                                 /* For all first neighbours of k1 */
 549                                 l1 = nnb->a[k1][1][l];
 550
 551                                 for (e = 0; e < excl->nr; e++)
 552                                 {
 553                                     if (excl->e[e] == l1)
 554                                     {
 555                                         remove_excl(excl, e);
 556                                     }
 557                                 }
 558                             }
 559                         }
 560                     }
 561                 }
 562             }
 563         }
 564     }
 565 }
 566
 567 /*! \brief
 568  * Generate pairs, angles and dihedrals from .rtp settings
 569  *
 570  * \param[in,out] atoms            Global information about atoms in topology.
 571  * \param[in]     rtpFFDB          Residue type database from force field.
 572  * \param[in,out] plist            Information about listed interactions.
 573  * \param[in,out] excls            Pair interaction exclusions.
 574  * \param[in,out] globalPatches    Information about possible residue modifications.
 575  * \param[in]     bAllowMissing    True if missing interaction information is allowed.
 576  *                                 AKA allow cartoon physics
 577  * \param[in]     cyclicBondsIndex Information about bonds creating cyclic molecules.
 578  *                                 Empty if no such bonds exist.
 579  */
 580 void gen_pad(t_atoms*                               atoms,
 581              gmx::ArrayRef<const PreprocessResidue> rtpFFDB,
 582              gmx::ArrayRef<InteractionsOfType>      plist,
 583              t_excls                                excls[],
 584              gmx::ArrayRef<MoleculePatchDatabase>   globalPatches,
 585              bool                                   bAllowMissing,
 586              gmx::ArrayRef<const int>               cyclicBondsIndex)
 587 {
 588     t_nextnb nnb;
 589     init_nnb(&nnb, atoms->nr, 4);
 590     gen_nnb(&nnb, plist);
 591     print_nnb(&nnb, "NNB");
 592
 593     /* These are the angles, dihedrals and pairs that we generate
 594      * from the bonds. The ones that are already there from the rtp file
 595      * will be retained.
 596      */
 597     std::vector<InteractionOfType> ang;
 598     std::vector<InteractionOfType> dih;
 599     std::vector<InteractionOfType> pai;
 600     std::vector<InteractionOfType> improper;
 601
 602     std::array<std::string, 4> anm;
 603
 604     if (!globalPatches.empty())
 605     {
 606         gen_excls(atoms, excls, globalPatches, bAllowMissing, cyclicBondsIndex);
 607         /* mark all entries as not matched yet */
 608         for (int i = 0; i < atoms->nres; i++)
 609         {
 610             for (auto bondedsList : globalPatches[i].rb)
 611             {
 612                 for (auto& bondeds : bondedsList.b)
 613                 {
 614                     bondeds.match = false;
 615                 }
 616             }
 617         }
 618     }
 619
 620     /* Extract all i-j-k-l neighbours from nnb struct to generate all
 621      * angles and dihedrals. */
 622     for (int i = 0; (i < nnb.nr); i++)
 623     {
 624         /* For all particles */
 625         for (int j = 0; (j < nnb.nrexcl[i][1]); j++)
 626         {
 627             /* For all first neighbours */
 628             int j1 = nnb.a[i][1][j];
 629             for (int k = 0; (k < nnb.nrexcl[j1][1]); k++)
 630             {
 631                 /* For all first neighbours of j1 */
 632                 int k1 = nnb.a[j1][1][k];
 633                 if (k1 != i)
 634                 {
 635                     /* Generate every angle only once */
 636                     if (i < k1)
 637                     {
 638                         std::vector<int> atomNumbers = { i, j1, k1 };
 639                         std::string      name;
 640                         if (!globalPatches.empty())
 641                         {
 642                             int minres = atoms->atom[i].resind;
 643                             int maxres = minres;
 644                             for (int m = 1; m < 3; m++)
 645                             {
 646                                 minres = std::min(minres, atoms->atom[atomNumbers[m]].resind);
 647                                 maxres = std::max(maxres, atoms->atom[atomNumbers[m]].resind);
 648                             }
 649                             int res = 2 * minres - maxres;
 650                             do
 651                             {
 652                                 res += maxres - minres;
 653                                 get_atomnames_min(3, anm, res, atoms, atomNumbers);
 654                                 BondedInteractionList* hbang = &globalPatches[res].rb[BondedTypes::Angles];
 655                                 for (auto& bondeds : hbang->b)
 656                                 {
 657                                     if (anm[1] == bondeds.aj())
 658                                     {
 659                                         bool bFound = false;
 660                                         for (int m = 0; m < 3; m += 2)
 661                                         {
 662                                             bFound = (bFound
 663                                                       || ((anm[m] == bondeds.ai())
 664                                                           && (anm[2 - m] == bondeds.ak())));
 665                                         }
 666                                         if (bFound)
 667                                         {
 668                                             name = bondeds.s;
 669                                             /* Mark that we found a match for this entry */
 670                                             bondeds.match = true;
 671                                         }
 672                                     }
 673                                 }
 674                             } while (res < maxres);
 675                         }
 676                         ang.push_back(InteractionOfType(atomNumbers, {}, name));
 677                     }
 678                     /* Generate every dihedral, 1-4 exclusion and 1-4 interaction
 679                        only once */
 680                     if (j1 < k1)
 681                     {
 682                         for (int l = 0; (l < nnb.nrexcl[k1][1]); l++)
 683                         {
 684                             /* For all first neighbours of k1 */
 685                             int l1 = nnb.a[k1][1][l];
 686                             if ((l1 != i) && (l1 != j1))
 687                             {
 688                                 std::vector<int> atomNumbers = { i, j1, k1, l1 };
 689                                 std::string      name;
 690                                 int              nFound = 0;
 691                                 if (!globalPatches.empty())
 692                                 {
 693                                     int minres = atoms->atom[i].resind;
 694                                     int maxres = minres;
 695                                     for (int m = 1; m < 4; m++)
 696                                     {
 697                                         minres = std::min(minres, atoms->atom[atomNumbers[m]].resind);
 698                                         maxres = std::max(maxres, atoms->atom[atomNumbers[m]].resind);
 699                                     }
 700                                     int res = 2 * minres - maxres;
 701                                     do
 702                                     {
 703                                         res += maxres - minres;
 704                                         get_atomnames_min(4, anm, res, atoms, atomNumbers);
 705                                         BondedInteractionList* hbdih =
 706                                                 &globalPatches[res].rb[BondedTypes::ProperDihedrals];
 707                                         for (auto& bondeds : hbdih->b)
 708                                         {
 709                                             bool bFound = false;
 710                                             for (int m = 0; m < 2; m++)
 711                                             {
 712                                                 bFound = (bFound
 713                                                           || ((anm[3 * m] == bondeds.ai())
 714                                                               && (anm[1 + m] == bondeds.aj())
 715                                                               && (anm[2 - m] == bondeds.ak())
 716                                                               && (anm[3 - 3 * m] == bondeds.al())));
 717                                             }
 718                                             if (bFound)
 719                                             {
 720                                                 name = bondeds.s;
 721                                                 /* Mark that we found a match for this entry */
 722                                                 bondeds.match = true;
 723
 724                                                 /* Set the last parameter to be able to see
 725                                                    if the dihedral was in the rtp list.
 726                                                  */
 727                                                 nFound++;
 728                                                 dih.push_back(InteractionOfType(atomNumbers, {}, name));
 729                                                 dih.back().setForceParameter(
 730                                                         MAXFORCEPARAM - 1, DIHEDRAL_WAS_SET_IN_RTP);
 731                                             }
 732                                         }
 733                                     } while (res < maxres);
 734                                 }
 735                                 if (nFound == 0)
 736                                 {
 737                                     std::vector<int> atoms = { i, j1, k1, l1 };
 738                                     dih.push_back(InteractionOfType(atoms, {}, ""));
 739                                 }
 740
 741                                 int nbd = nb_dist(&nnb, i, l1);
 742                                 if (nbd == 3)
 743                                 {
 744                                     int  i1    = std::min(i, l1);
 745                                     int  i2    = std::max(i, l1);
 746                                     bool bExcl = false;
 747                                     for (int m = 0; m < excls[i1].nr; m++)
 748                                     {
 749                                         bExcl = bExcl || excls[i1].e[m] == i2;
 750                                     }
 751                                     if (!bExcl)
 752                                     {
 753                                         if (rtpFFDB[0].bGenerateHH14Interactions
 754                                             || !(is_hydro(atoms, i1) && is_hydro(atoms, i2)))
 755                                         {
 756                                             std::vector<int> atoms = { i1, i2 };
 757                                             pai.push_back(InteractionOfType(atoms, {}, ""));
 758                                         }
 759                                     }
 760                                 }
 761                             }
 762                         }
 763                     }
 764                 }
 765             }
 766         }
 767     }
 768
 769     if (!globalPatches.empty())
 770     {
 771         /* The above approach is great in that we double-check that e.g. an angle
 772          * really corresponds to three atoms connected by bonds, but this is not
 773          * generally true. Go through the angle and dihedral hackblocks to add
 774          * entries that we have not yet marked as matched when going through bonds.
 775          */
 776         for (int i = 0; i < atoms->nres; i++)
 777         {
 778             /* Add remaining angles from hackblock */
 779             BondedInteractionList* hbang = &globalPatches[i].rb[BondedTypes::Angles];
 780             for (auto& bondeds : hbang->b)
 781             {
 782                 if (bondeds.match)
 783                 {
 784                     /* We already used this entry, continue to the next */
 785                     continue;
 786                 }
 787                 /* Hm - entry not used, let's see if we can find all atoms */
 788                 std::vector<int>                   atomNumbers;
 789                 bool                               bFound = true;
 790                 gmx::ArrayRef<const int>::iterator cyclicBondsIterator;
 791                 for (int k = 0; k < 3 && bFound; k++)
 792                 {
 793                     const char* p   = bondeds.a[k].c_str();
 794                     int         res = i;
 795                     if (p[0] == '-')
 796                     {
 797                         p++;
 798                         cyclicBondsIterator =
 799                                 std::find(cyclicBondsIndex.begin(), cyclicBondsIndex.end(), res--);
 800                         if (cyclicBondsIterator != cyclicBondsIndex.end()
 801                             && !((cyclicBondsIterator - cyclicBondsIndex.begin()) & 1))
 802                         {
 803                             res = *(++cyclicBondsIterator);
 804                         }
 805                     }
 806                     else if (p[0] == '+')
 807                     {
 808                         p++;
 809                         cyclicBondsIterator =
 810                                 std::find(cyclicBondsIndex.begin(), cyclicBondsIndex.end(), res++);
 811                         if (cyclicBondsIterator != cyclicBondsIndex.end()
 812                             && ((cyclicBondsIterator - cyclicBondsIndex.begin()) & 1))
 813                         {
 814                             res = *(--cyclicBondsIterator);
 815                         }
 816                     }
 817                     atomNumbers.emplace_back(search_res_atom(p, res, atoms, "angle", TRUE));
 818                     bFound = (atomNumbers.back() != -1);
 819                 }
 820
 821                 if (bFound)
 822                 {
 823                     bondeds.match = true;
 824                     /* Incrementing nang means we save this angle */
 825                     ang.push_back(InteractionOfType(atomNumbers, {}, bondeds.s));
 826                 }
 827             }
 828
 829             /* Add remaining dihedrals from hackblock */
 830             BondedInteractionList* hbdih = &globalPatches[i].rb[BondedTypes::ProperDihedrals];
 831             for (auto& bondeds : hbdih->b)
 832             {
 833                 if (bondeds.match)
 834                 {
 835                     /* We already used this entry, continue to the next */
 836                     continue;
 837                 }
 838                 /* Hm - entry not used, let's see if we can find all atoms */
 839                 std::vector<int>                   atomNumbers;
 840                 bool                               bFound = true;
 841                 gmx::ArrayRef<const int>::iterator cyclicBondsIterator;
 842                 for (int k = 0; k < 4 && bFound; k++)
 843                 {
 844                     const char* p   = bondeds.a[k].c_str();
 845                     int         res = i;
 846                     if (p[0] == '-')
 847                     {
 848                         p++;
 849                         cyclicBondsIterator =
 850                                 std::find(cyclicBondsIndex.begin(), cyclicBondsIndex.end(), res);
 851                         res--;
 852                         if (cyclicBondsIterator != cyclicBondsIndex.end()
 853                             && !((cyclicBondsIterator - cyclicBondsIndex.begin()) & 1))
 854                         {
 855                             res = *(++cyclicBondsIterator);
 856                         }
 857                     }
 858                     else if (p[0] == '+')
 859                     {
 860                         p++;
 861                         cyclicBondsIterator =
 862                                 std::find(cyclicBondsIndex.begin(), cyclicBondsIndex.end(), res);
 863                         res++;
 864                         if (cyclicBondsIterator != cyclicBondsIndex.end()
 865                             && ((cyclicBondsIterator - cyclicBondsIndex.begin()) & 1))
 866                         {
 867                             res = *(--cyclicBondsIterator);
 868                         }
 869                     }
 870                     atomNumbers.emplace_back(search_res_atom(p, res, atoms, "dihedral", TRUE));
 871                     bFound = (atomNumbers.back() != -1);
 872                 }
 873
 874                 if (bFound)
 875                 {
 876                     bondeds.match = true;
 877                     /* Incrementing ndih means we save this dihedral */
 878                     dih.push_back(InteractionOfType(atomNumbers, {}, bondeds.s));
 879                 }
 880             }
 881         }
 882     }
 883
 884     /* Sort angles with respect to j-i-k (middle atom first) */
 885     if (ang.size() > 1)
 886     {
 887         std::sort(ang.begin(), ang.end(), acomp);
 888     }
 889
 890     /* Sort dihedrals with respect to j-k-i-l (middle atoms first) */
 891     if (dih.size() > 1)
 892     {
 893         std::sort(dih.begin(), dih.end(), dcomp);
 894     }
 895
 896     /* Sort the pairs */
 897     if (pai.size() > 1)
 898     {
 899         std::sort(pai.begin(), pai.end(), pcomp);
 900     }
 901     if (!pai.empty())
 902     {
 903         /* Remove doubles, could occur in 6-rings, such as phenyls,
 904            maybe one does not want this when fudgeQQ < 1.
 905          */
 906         fprintf(stderr, "Before cleaning: %zu pairs\n", pai.size());
 907         rm2par(&pai, preq);
 908     }
 909
 910     /* Get the impropers from the database */
 911     improper = get_impropers(atoms, globalPatches, bAllowMissing, cyclicBondsIndex);
 912
 913     /* Sort the impropers */
 914     sort_id(improper);
 915
 916     if (!dih.empty())
 917     {
 918         fprintf(stderr, "Before cleaning: %zu dihedrals\n", dih.size());
 919         dih = clean_dih(dih, improper, atoms, rtpFFDB[0].bKeepAllGeneratedDihedrals, rtpFFDB[0].bRemoveDihedralIfWithImproper);
 920     }
 921
 922     /* Now we have unique lists of angles and dihedrals
 923      * Copy them into the destination struct
 924      */
 925     cppar(ang, plist, F_ANGLES);
 926     cppar(dih, plist, F_PDIHS);
 927     cppar(improper, plist, F_IDIHS);
 928     cppar(pai, plist, F_LJ14);
 929
 930     /* Remove all exclusions which are within nrexcl */
 931     clean_excls(&nnb, rtpFFDB[0].nrexcl, excls);
 932     done_nnb(&nnb);
 933 }