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