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