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