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