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