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, 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 <ctype.h>
  43 #include <stdlib.h>
  44 #include <string.h>
  45
  46 #include <cmath>
  47
  48 #include <algorithm>
  49
  50 #include "gromacs/fileio/confio.h"
  51 #include "gromacs/gmxpreprocess/gpp_nextnb.h"
  52 #include "gromacs/gmxpreprocess/notset.h"
  53 #include "gromacs/gmxpreprocess/pgutil.h"
  54 #include "gromacs/gmxpreprocess/resall.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 #define DIHEDRAL_WAS_SET_IN_RTP 0
  64 static bool was_dihedral_set_in_rtp(const t_param *dih)
  65 {
  66     return dih->c[MAXFORCEPARAM-1] == DIHEDRAL_WAS_SET_IN_RTP;
  67 }
  68
  69 typedef bool (*peq)(t_param *p1, t_param *p2);
  70
  71 static int acomp(const void *a1, const void *a2)
  72 {
  73     const t_param *p1, *p2;
  74     int            ac;
  75
  76     p1 = static_cast<const t_param *>(a1);
  77     p2 = static_cast<const t_param *>(a2);
  78     if ((ac = (p1->aj()-p2->aj())) != 0)
  79     {
  80         return ac;
  81     }
  82     else if ((ac = (p1->ai()-p2->ai())) != 0)
  83     {
  84         return ac;
  85     }
  86     else
  87     {
  88         return (p1->ak()-p2->ak());
  89     }
  90 }
  91
  92 static int pcomp(const void *a1, const void *a2)
  93 {
  94     const t_param *p1, *p2;
  95     int            pc;
  96
  97     p1 = static_cast<const t_param *>(a1);
  98     p2 = static_cast<const t_param *>(a2);
  99     if ((pc = (p1->ai()-p2->ai())) != 0)
 100     {
 101         return pc;
 102     }
 103     else
 104     {
 105         return (p1->aj()-p2->aj());
 106     }
 107 }
 108
 109 static int dcomp(const void *d1, const void *d2)
 110 {
 111     const t_param *p1, *p2;
 112     int            dc;
 113
 114     p1 = static_cast<const t_param *>(d1);
 115     p2 = static_cast<const t_param *>(d2);
 116     /* First sort by J & K (the two central) atoms */
 117     if ((dc = (p1->aj()-p2->aj())) != 0)
 118     {
 119         return dc;
 120     }
 121     else if ((dc = (p1->ak()-p2->ak())) != 0)
 122     {
 123         return dc;
 124     }
 125     /* Then make sure to put rtp dihedrals before generated ones */
 126     else if (was_dihedral_set_in_rtp(p1) &&
 127              !was_dihedral_set_in_rtp(p2))
 128     {
 129         return -1;
 130     }
 131     else if (!was_dihedral_set_in_rtp(p1) &&
 132              was_dihedral_set_in_rtp(p2))
 133     {
 134         return 1;
 135     }
 136     /* Then sort by I and J (two outer) atoms */
 137     else if ((dc = (p1->ai()-p2->ai())) != 0)
 138     {
 139         return dc;
 140     }
 141     else if ((dc = (p1->al()-p2->al())) != 0)
 142     {
 143         return dc;
 144     }
 145     else
 146     {
 147         // AMBER force fields with type 9 dihedrals can reach here, where we sort on
 148         // the contents of the string that names the macro for the parameters.
 149         return strcmp(p1->s, p2->s);
 150     }
 151 }
 152
 153
 154 static bool is_dihedral_on_same_bond(t_param *p1, t_param *p2)
 155 {
 156     return ((p1->aj() == p2->aj()) && (p1->ak() == p2->ak())) ||
 157            ((p1->aj() == p2->ak()) && (p1->ak() == p2->aj()));
 158 }
 159
 160
 161 static bool preq(t_param *p1, t_param *p2)
 162 {
 163     return (p1->ai() == p2->ai()) && (p1->aj() == p2->aj());
 164 }
 165
 166 static void rm2par(t_param p[], int *np, peq eq)
 167 {
 168     int *index, nind;
 169     int  i, j;
 170
 171     if ((*np) == 0)
 172     {
 173         return;
 174     }
 175
 176     snew(index, *np);
 177     nind          = 0;
 178     index[nind++] = 0;
 179     for (i = 1; (i < (*np)); i++)
 180     {
 181         if (!eq(&p[i], &p[i-1]))
 182         {
 183             index[nind++] = i;
 184         }
 185     }
 186     /* Index now holds pointers to all the non-equal params,
 187      * this only works when p is sorted of course
 188      */
 189     for (i = 0; (i < nind); i++)
 190     {
 191         for (j = 0; (j < MAXATOMLIST); j++)
 192         {
 193             p[i].a[j] = p[index[i]].a[j];
 194         }
 195         for (j = 0; (j < MAXFORCEPARAM); j++)
 196         {
 197             p[i].c[j] = p[index[i]].c[j];
 198         }
 199         if (p[index[i]].a[0] == p[index[i]].a[1])
 200         {
 201             strcpy(p[i].s, "");
 202         }
 203         else if (index[i] > i)
 204         {
 205             /* Copy the string only if it comes from somewhere else
 206              * otherwise we will end up copying a random (newly freed) pointer.
 207              * Since the index is sorted we only have to test for index[i] > i.
 208              */
 209             strcpy(p[i].s, p[index[i]].s);
 210         }
 211     }
 212     (*np) = nind;
 213
 214     sfree(index);
 215 }
 216
 217 static void cppar(t_param p[], int np, t_params plist[], int ftype)
 218 {
 219     int       i, j, nral, nrfp;
 220     t_params *ps;
 221
 222     ps   = &plist[ftype];
 223     nral = NRAL(ftype);
 224     nrfp = NRFP(ftype);
 225
 226     /* Keep old stuff */
 227     pr_alloc(np, ps);
 228     for (i = 0; (i < np); i++)
 229     {
 230         for (j = 0; (j < nral); j++)
 231         {
 232             ps->param[ps->nr].a[j] = p[i].a[j];
 233         }
 234         for (j = 0; (j < nrfp); j++)
 235         {
 236             ps->param[ps->nr].c[j] = p[i].c[j];
 237         }
 238         for (j = 0; (j < MAXSLEN); j++)
 239         {
 240             ps->param[ps->nr].s[j] = p[i].s[j];
 241         }
 242         ps->nr++;
 243     }
 244 }
 245
 246 static void cpparam(t_param *dest, t_param *src)
 247 {
 248     int j;
 249
 250     for (j = 0; (j < MAXATOMLIST); j++)
 251     {
 252         dest->a[j] = src->a[j];
 253     }
 254     for (j = 0; (j < MAXFORCEPARAM); j++)
 255     {
 256         dest->c[j] = src->c[j];
 257     }
 258     for (j = 0; (j < MAXSLEN); j++)
 259     {
 260         dest->s[j] = src->s[j];
 261     }
 262 }
 263
 264 static void set_p(t_param *p, const int ai[4], const real *c, char *s)
 265 {
 266     int j;
 267
 268     for (j = 0; (j < 4); j++)
 269     {
 270         p->a[j] = ai[j];
 271     }
 272     for (j = 0; (j < MAXFORCEPARAM); j++)
 273     {
 274         if (c)
 275         {
 276             p->c[j] = c[j];
 277         }
 278         else
 279         {
 280             p->c[j] = NOTSET;
 281         }
 282     }
 283
 284     set_p_string(p, s);
 285 }
 286
 287 static int int_comp(const void *a, const void *b)
 288 {
 289     return (*(int *)a) - (*(int *)b);
 290 }
 291
 292 static int idcomp(const void *a, const void *b)
 293 {
 294     const t_param *pa, *pb;
 295     int            d;
 296
 297     pa = static_cast<const t_param *>(a);
 298     pb = static_cast<const t_param *>(b);
 299     if ((d = (pa->a[0]-pb->a[0])) != 0)
 300     {
 301         return d;
 302     }
 303     else if ((d = (pa->a[3]-pb->a[3])) != 0)
 304     {
 305         return d;
 306     }
 307     else if ((d = (pa->a[1]-pb->a[1])) != 0)
 308     {
 309         return d;
 310     }
 311     else
 312     {
 313         return (pa->a[2]-pb->a[2]);
 314     }
 315 }
 316
 317 static void sort_id(int nr, t_param ps[])
 318 {
 319     int i, tmp;
 320
 321     /* First swap order of atoms around if necessary */
 322     for (i = 0; (i < nr); i++)
 323     {
 324         if (ps[i].a[3] < ps[i].a[0])
 325         {
 326             tmp = ps[i].a[3]; ps[i].a[3] = ps[i].a[0]; ps[i].a[0] = tmp;
 327             tmp = ps[i].a[2]; ps[i].a[2] = ps[i].a[1]; ps[i].a[1] = tmp;
 328         }
 329     }
 330     /* Now sort it */
 331     if (nr > 1)
 332     {
 333         qsort(ps, nr, static_cast<size_t>(sizeof(ps[0])), idcomp);
 334     }
 335 }
 336
 337 static int n_hydro(const int a[], char ***atomname)
 338 {
 339     int  i, nh = 0;
 340     char c0, c1, *aname;
 341
 342     for (i = 0; (i < 4); i += 3)
 343     {
 344         aname = *atomname[a[i]];
 345         c0    = toupper(aname[0]);
 346         if (c0 == 'H')
 347         {
 348             nh++;
 349         }
 350         else if ((static_cast<int>(strlen(aname)) > 1) && (c0 >= '0') && (c0 <= '9'))
 351         {
 352             c1 = toupper(aname[1]);
 353             if (c1 == 'H')
 354             {
 355                 nh++;
 356             }
 357         }
 358     }
 359     return nh;
 360 }
 361
 362 /* Clean up the dihedrals (both generated and read from the .rtp
 363  * file). */
 364 static void clean_dih(t_param *dih, int *ndih, t_param improper[], int nimproper,
 365                       t_atoms *atoms, bool bKeepAllGeneratedDihedrals,
 366                       bool bRemoveDihedralIfWithImproper)
 367 {
 368     int   i, j, k, l;
 369     int  *index, nind;
 370
 371     /* Construct the list of the indices of the dihedrals
 372      * (i.e. generated or read) that might be kept. */
 373     snew(index, *ndih+1);
 374     if (bKeepAllGeneratedDihedrals)
 375     {
 376         fprintf(stderr, "Keeping all generated dihedrals\n");
 377         nind = *ndih;
 378         for (i = 0; i < nind; i++)
 379         {
 380             index[i] = i;
 381         }
 382         index[nind] = *ndih;
 383     }
 384     else
 385     {
 386         nind = 0;
 387         /* Check if generated dihedral i should be removed. The
 388          * dihedrals have been sorted by dcomp() above, so all those
 389          * on the same two central atoms are together, with those from
 390          * the .rtp file preceding those that were automatically
 391          * generated. We remove the latter if the former exist. */
 392         for (i = 0; i < *ndih; i++)
 393         {
 394             /* Keep the dihedrals that were defined in the .rtp file,
 395              * and the dihedrals that were generated and different
 396              * from the last one (whether it was generated or not). */
 397             if (was_dihedral_set_in_rtp(&dih[i]) ||
 398                 0 == i ||
 399                 !is_dihedral_on_same_bond(&dih[i], &dih[i-1]))
 400             {
 401                 index[nind++] = i;
 402             }
 403         }
 404         index[nind] = *ndih;
 405     }
 406
 407     k = 0;
 408     for (i = 0; i < nind; i++)
 409     {
 410         bool bWasSetInRTP = was_dihedral_set_in_rtp(&dih[index[i]]);
 411         bool bKeep        = TRUE;
 412         if (!bWasSetInRTP && bRemoveDihedralIfWithImproper)
 413         {
 414             /* Remove the dihedral if there is an improper on the same
 415              * bond. */
 416             for (j = 0; j < nimproper && bKeep; j++)
 417             {
 418                 bKeep = !is_dihedral_on_same_bond(&dih[index[i]], &improper[j]);
 419             }
 420         }
 421
 422         if (bKeep)
 423         {
 424             /* If we don't want all dihedrals, we want to select the
 425              * ones with the fewest hydrogens. Note that any generated
 426              * dihedrals on the same bond as an .rtp dihedral may have
 427              * been already pruned above in the construction of
 428              * index[]. However, their parameters are still present,
 429              * and l is looping over this dihedral and all of its
 430              * pruned siblings. */
 431             int bestl = index[i];
 432             if (!bKeepAllGeneratedDihedrals && !bWasSetInRTP)
 433             {
 434                 /* Minimum number of hydrogens for i and l atoms */
 435                 int minh = 2;
 436                 for (l = index[i];
 437                      (l < index[i+1] &&
 438                       is_dihedral_on_same_bond(&dih[index[i]], &dih[l]));
 439                      l++)
 440                 {
 441                     int nh = n_hydro(dih[l].a, atoms->atomname);
 442                     if (nh < minh)
 443                     {
 444                         minh  = nh;
 445                         bestl = l;
 446                     }
 447                     if (0 == minh)
 448                     {
 449                         break;
 450                     }
 451                 }
 452             }
 453             if (k != bestl)
 454             {
 455                 cpparam(&dih[k], &dih[bestl]);
 456             }
 457             k++;
 458         }
 459     }
 460
 461     for (i = k; i < *ndih; i++)
 462     {
 463         strcpy(dih[i].s, "");
 464     }
 465     *ndih = k;
 466
 467     sfree(index);
 468 }
 469
 470 static int get_impropers(t_atoms *atoms, t_hackblock hb[], t_param **improper,
 471                          bool bAllowMissing)
 472 {
 473     t_rbondeds   *impropers;
 474     int           nimproper, i, j, k, start, ninc, nalloc;
 475     int           ai[MAXATOMLIST];
 476     bool          bStop;
 477
 478     ninc   = 500;
 479     nalloc = ninc;
 480     snew(*improper, nalloc);
 481
 482     /* Add all the impropers from the residue database to the list. */
 483     nimproper = 0;
 484     start     = 0;
 485     if (hb != nullptr)
 486     {
 487         for (i = 0; (i < atoms->nres); i++)
 488         {
 489             impropers = &hb[i].rb[ebtsIDIHS];
 490             for (j = 0; (j < impropers->nb); j++)
 491             {
 492                 bStop = FALSE;
 493                 for (k = 0; (k < 4) && !bStop; k++)
 494                 {
 495                     ai[k] = search_atom(impropers->b[j].a[k], start,
 496                                         atoms,
 497                                         "improper", bAllowMissing);
 498                     if (ai[k] == -1)
 499                     {
 500                         bStop = TRUE;
 501                     }
 502                 }
 503                 if (!bStop)
 504                 {
 505                     if (nimproper == nalloc)
 506                     {
 507                         nalloc += ninc;
 508                         srenew(*improper, nalloc);
 509                     }
 510                     /* Not broken out */
 511                     set_p(&((*improper)[nimproper]), ai, nullptr, impropers->b[j].s);
 512                     nimproper++;
 513                 }
 514             }
 515             while ((start < atoms->nr) && (atoms->atom[start].resind == i))
 516             {
 517                 start++;
 518             }
 519         }
 520     }
 521
 522     return nimproper;
 523 }
 524
 525 static int nb_dist(t_nextnb *nnb, int ai, int aj)
 526 {
 527     int  nre, nrx, NRE;
 528     int *nrexcl;
 529     int *a;
 530
 531     if (ai == aj)
 532     {
 533         return 0;
 534     }
 535
 536     NRE    = -1;
 537     nrexcl = nnb->nrexcl[ai];
 538     for (nre = 1; (nre < nnb->nrex); nre++)
 539     {
 540         a = nnb->a[ai][nre];
 541         for (nrx = 0; (nrx < nrexcl[nre]); nrx++)
 542         {
 543             if ((aj == a[nrx]) && (NRE == -1))
 544             {
 545                 NRE = nre;
 546             }
 547         }
 548     }
 549     return NRE;
 550 }
 551
 552 static bool is_hydro(t_atoms *atoms, int ai)
 553 {
 554     return ((*(atoms->atomname[ai]))[0] == 'H');
 555 }
 556
 557 static void get_atomnames_min(int n, char **anm,
 558                               int resind, t_atoms *atoms, const int *a)
 559 {
 560     int m;
 561
 562     /* Assume ascending residue numbering */
 563     for (m = 0; m < n; m++)
 564     {
 565         if (atoms->atom[a[m]].resind < resind)
 566         {
 567             strcpy(anm[m], "-");
 568         }
 569         else if (atoms->atom[a[m]].resind > resind)
 570         {
 571             strcpy(anm[m], "+");
 572         }
 573         else
 574         {
 575             strcpy(anm[m], "");
 576         }
 577         strcat(anm[m], *(atoms->atomname[a[m]]));
 578     }
 579 }
 580
 581 static void gen_excls(t_atoms *atoms, t_excls *excls, t_hackblock hb[],
 582                       bool bAllowMissing)
 583 {
 584     int         r;
 585     int         a, astart, i1, i2, itmp;
 586     t_rbondeds *hbexcl;
 587     int         e;
 588     char       *anm;
 589
 590     astart = 0;
 591     for (a = 0; a < atoms->nr; a++)
 592     {
 593         r = atoms->atom[a].resind;
 594         if (a == atoms->nr-1 || atoms->atom[a+1].resind != r)
 595         {
 596             hbexcl = &hb[r].rb[ebtsEXCLS];
 597
 598             for (e = 0; e < hbexcl->nb; e++)
 599             {
 600                 anm = hbexcl->b[e].a[0];
 601                 i1  = search_atom(anm, astart, atoms,
 602                                   "exclusion", bAllowMissing);
 603                 anm = hbexcl->b[e].a[1];
 604                 i2  = search_atom(anm, astart, atoms,
 605                                   "exclusion", bAllowMissing);
 606                 if (i1 != -1 && i2 != -1)
 607                 {
 608                     if (i1 > i2)
 609                     {
 610                         itmp = i1;
 611                         i1   = i2;
 612                         i2   = itmp;
 613                     }
 614                     srenew(excls[i1].e, excls[i1].nr+1);
 615                     excls[i1].e[excls[i1].nr] = i2;
 616                     excls[i1].nr++;
 617                 }
 618             }
 619
 620             astart = a+1;
 621         }
 622     }
 623
 624     for (a = 0; a < atoms->nr; a++)
 625     {
 626         if (excls[a].nr > 1)
 627         {
 628             qsort(excls[a].e, excls[a].nr, static_cast<size_t>(sizeof(int)), int_comp);
 629         }
 630     }
 631 }
 632
 633 static void remove_excl(t_excls *excls, int remove)
 634 {
 635     int i;
 636
 637     for (i = remove+1; i < excls->nr; i++)
 638     {
 639         excls->e[i-1] = excls->e[i];
 640     }
 641
 642     excls->nr--;
 643 }
 644
 645 void clean_excls(t_nextnb *nnb, int nrexcl, t_excls excls[])
 646 {
 647     int      i, j, j1, k, k1, l, l1, e;
 648     t_excls *excl;
 649
 650     if (nrexcl >= 1)
 651     {
 652         /* extract all i-j-k-l neighbours from nnb struct */
 653         for (i = 0; (i < nnb->nr); i++)
 654         {
 655             /* For all particles */
 656             excl = &excls[i];
 657
 658             for (j = 0; (j < nnb->nrexcl[i][1]); j++)
 659             {
 660                 /* For all first neighbours */
 661                 j1 = nnb->a[i][1][j];
 662
 663                 for (e = 0; e < excl->nr; e++)
 664                 {
 665                     if (excl->e[e] == j1)
 666                     {
 667                         remove_excl(excl, e);
 668                     }
 669                 }
 670
 671                 if (nrexcl >= 2)
 672                 {
 673                     for (k = 0; (k < nnb->nrexcl[j1][1]); k++)
 674                     {
 675                         /* For all first neighbours of j1 */
 676                         k1 = nnb->a[j1][1][k];
 677
 678                         for (e = 0; e < excl->nr; e++)
 679                         {
 680                             if (excl->e[e] == k1)
 681                             {
 682                                 remove_excl(excl, e);
 683                             }
 684                         }
 685
 686                         if (nrexcl >= 3)
 687                         {
 688                             for (l = 0; (l < nnb->nrexcl[k1][1]); l++)
 689                             {
 690                                 /* For all first neighbours of k1 */
 691                                 l1 = nnb->a[k1][1][l];
 692
 693                                 for (e = 0; e < excl->nr; e++)
 694                                 {
 695                                     if (excl->e[e] == l1)
 696                                     {
 697                                         remove_excl(excl, e);
 698                                     }
 699                                 }
 700                             }
 701                         }
 702                     }
 703                 }
 704             }
 705         }
 706     }
 707 }
 708
 709 void generate_excls(t_nextnb *nnb, int nrexcl, t_excls excls[])
 710 {
 711     int      i, j, n, N;
 712     t_excls *excl;
 713
 714     for (N = 1; (N < std::min(nrexcl, nnb->nrex)); N++)
 715     {
 716         /* extract all i-j-k-l neighbours from nnb struct */
 717         for (i = 0; (i < nnb->nr); i++)
 718         {
 719             /* For all particles */
 720             excl      = &excls[i];
 721             n         = excl->nr;
 722             excl->nr += nnb->nrexcl[i][N];
 723             srenew(excl->e, excl->nr);
 724             for (j = 0; (j < nnb->nrexcl[i][N]); j++)
 725             {
 726                 /* For all first neighbours */
 727                 if (nnb->a[i][N][j] != i)
 728                 {
 729                     excl->e[n++] = nnb->a[i][N][j];
 730                 }
 731             }
 732         }
 733     }
 734 }
 735
 736 /* Generate pairs, angles and dihedrals from .rtp settings */
 737 void gen_pad(t_nextnb *nnb, t_atoms *atoms, t_restp rtp[],
 738              t_params plist[], t_excls excls[], t_hackblock hb[],
 739              bool bAllowMissing)
 740 {
 741     t_param    *ang, *dih, *pai, *improper;
 742     t_rbondeds *hbang, *hbdih;
 743     char      **anm;
 744     const char *p;
 745     int         res, minres, maxres;
 746     int         i, j, j1, k, k1, l, l1, m, n, i1, i2;
 747     int         ninc, maxang, maxdih, maxpai;
 748     int         nang, ndih, npai, nimproper, nbd;
 749     int         nFound;
 750     bool        bFound, bExcl;
 751
 752     /* These are the angles, dihedrals and pairs that we generate
 753      * from the bonds. The ones that are already there from the rtp file
 754      * will be retained.
 755      */
 756     nang   = 0;
 757     npai   = 0;
 758     ndih   = 0;
 759     ninc   = 500;
 760     maxang = maxdih = maxpai = ninc;
 761     snew(ang, maxang);
 762     snew(dih, maxdih);
 763     snew(pai, maxpai);
 764
 765     snew(anm, 4);
 766     for (i = 0; i < 4; i++)
 767     {
 768         snew(anm[i], 12);
 769     }
 770
 771     if (hb)
 772     {
 773         gen_excls(atoms, excls, hb, bAllowMissing);
 774         /* mark all entries as not matched yet */
 775         for (i = 0; i < atoms->nres; i++)
 776         {
 777             for (j = 0; j < ebtsNR; j++)
 778             {
 779                 for (k = 0; k < hb[i].rb[j].nb; k++)
 780                 {
 781                     hb[i].rb[j].b[k].match = FALSE;
 782                 }
 783             }
 784         }
 785     }
 786
 787     /* Extract all i-j-k-l neighbours from nnb struct to generate all
 788      * angles and dihedrals. */
 789     for (i = 0; (i < nnb->nr); i++)
 790     {
 791         /* For all particles */
 792         for (j = 0; (j < nnb->nrexcl[i][1]); j++)
 793         {
 794             /* For all first neighbours */
 795             j1 = nnb->a[i][1][j];
 796             for (k = 0; (k < nnb->nrexcl[j1][1]); k++)
 797             {
 798                 /* For all first neighbours of j1 */
 799                 k1 = nnb->a[j1][1][k];
 800                 if (k1 != i)
 801                 {
 802                     /* Generate every angle only once */
 803                     if (i < k1)
 804                     {
 805                         if (nang == maxang)
 806                         {
 807                             maxang += ninc;
 808                             srenew(ang, maxang);
 809                         }
 810                         ang[nang].ai() = i;
 811                         ang[nang].aj() = j1;
 812                         ang[nang].ak() = k1;
 813                         ang[nang].c0() = NOTSET;
 814                         ang[nang].c1() = NOTSET;
 815                         set_p_string(&(ang[nang]), "");
 816                         if (hb)
 817                         {
 818                             minres = atoms->atom[ang[nang].a[0]].resind;
 819                             maxres = minres;
 820                             for (m = 1; m < 3; m++)
 821                             {
 822                                 minres = std::min(minres, atoms->atom[ang[nang].a[m]].resind);
 823                                 maxres = std::max(maxres, atoms->atom[ang[nang].a[m]].resind);
 824                             }
 825                             res = 2*minres-maxres;
 826                             do
 827                             {
 828                                 res += maxres-minres;
 829                                 get_atomnames_min(3, anm, res, atoms, ang[nang].a);
 830                                 hbang = &hb[res].rb[ebtsANGLES];
 831                                 for (l = 0; (l < hbang->nb); l++)
 832                                 {
 833                                     if (strcmp(anm[1], hbang->b[l].aj()) == 0)
 834                                     {
 835                                         bFound = FALSE;
 836                                         for (m = 0; m < 3; m += 2)
 837                                         {
 838                                             bFound = (bFound ||
 839                                                       ((strcmp(anm[m], hbang->b[l].ai()) == 0) &&
 840                                                        (strcmp(anm[2-m], hbang->b[l].ak()) == 0)));
 841                                         }
 842                                         if (bFound)
 843                                         {
 844                                             set_p_string(&(ang[nang]), hbang->b[l].s);
 845                                             /* Mark that we found a match for this entry */
 846                                             hbang->b[l].match = TRUE;
 847                                         }
 848                                     }
 849                                 }
 850                             }
 851                             while (res < maxres);
 852                         }
 853                         nang++;
 854                     }
 855                     /* Generate every dihedral, 1-4 exclusion and 1-4 interaction
 856                        only once */
 857                     if (j1 < k1)
 858                     {
 859                         for (l = 0; (l < nnb->nrexcl[k1][1]); l++)
 860                         {
 861                             /* For all first neighbours of k1 */
 862                             l1 = nnb->a[k1][1][l];
 863                             if ((l1 != i) && (l1 != j1))
 864                             {
 865                                 if (ndih == maxdih)
 866                                 {
 867                                     maxdih += ninc;
 868                                     srenew(dih, maxdih);
 869                                 }
 870                                 dih[ndih].ai() = i;
 871                                 dih[ndih].aj() = j1;
 872                                 dih[ndih].ak() = k1;
 873                                 dih[ndih].al() = l1;
 874                                 for (m = 0; m < MAXFORCEPARAM; m++)
 875                                 {
 876                                     dih[ndih].c[m] = NOTSET;
 877                                 }
 878                                 set_p_string(&(dih[ndih]), "");
 879                                 nFound = 0;
 880                                 if (hb)
 881                                 {
 882                                     minres = atoms->atom[dih[ndih].a[0]].resind;
 883                                     maxres = minres;
 884                                     for (m = 1; m < 4; m++)
 885                                     {
 886                                         minres = std::min(minres, atoms->atom[dih[ndih].a[m]].resind);
 887                                         maxres = std::max(maxres, atoms->atom[dih[ndih].a[m]].resind);
 888                                     }
 889                                     res = 2*minres-maxres;
 890                                     do
 891                                     {
 892                                         res += maxres-minres;
 893                                         get_atomnames_min(4, anm, res, atoms, dih[ndih].a);
 894                                         hbdih = &hb[res].rb[ebtsPDIHS];
 895                                         for (n = 0; (n < hbdih->nb); n++)
 896                                         {
 897                                             bFound = FALSE;
 898                                             for (m = 0; m < 2; m++)
 899                                             {
 900                                                 bFound = (bFound ||
 901                                                           ((strcmp(anm[3*m],  hbdih->b[n].ai()) == 0) &&
 902                                                            (strcmp(anm[1+m],  hbdih->b[n].aj()) == 0) &&
 903                                                            (strcmp(anm[2-m],  hbdih->b[n].ak()) == 0) &&
 904                                                            (strcmp(anm[3-3*m], hbdih->b[n].al()) == 0)));
 905                                             }
 906                                             if (bFound)
 907                                             {
 908                                                 set_p_string(&dih[ndih], hbdih->b[n].s);
 909                                                 /* Mark that we found a match for this entry */
 910                                                 hbdih->b[n].match = TRUE;
 911
 912                                                 /* Set the last parameter to be able to see
 913                                                    if the dihedral was in the rtp list.
 914                                                  */
 915                                                 dih[ndih].c[MAXFORCEPARAM-1] = DIHEDRAL_WAS_SET_IN_RTP;
 916                                                 nFound++;
 917                                                 ndih++;
 918                                                 /* Set the next direct in case the rtp contains
 919                                                    multiple entries for this dihedral.
 920                                                  */
 921                                                 if (ndih == maxdih)
 922                                                 {
 923                                                     maxdih += ninc;
 924                                                     srenew(dih, maxdih);
 925                                                 }
 926                                                 dih[ndih].ai() = i;
 927                                                 dih[ndih].aj() = j1;
 928                                                 dih[ndih].ak() = k1;
 929                                                 dih[ndih].al() = l1;
 930                                                 for (m = 0; m < MAXFORCEPARAM; m++)
 931                                                 {
 932                                                     dih[ndih].c[m] = NOTSET;
 933                                                 }
 934                                             }
 935                                         }
 936                                     }
 937                                     while (res < maxres);
 938                                 }
 939                                 if (nFound == 0)
 940                                 {
 941                                     if (ndih == maxdih)
 942                                     {
 943                                         maxdih += ninc;
 944                                         srenew(dih, maxdih);
 945                                     }
 946                                     dih[ndih].ai() = i;
 947                                     dih[ndih].aj() = j1;
 948                                     dih[ndih].ak() = k1;
 949                                     dih[ndih].al() = l1;
 950                                     for (m = 0; m < MAXFORCEPARAM; m++)
 951                                     {
 952                                         dih[ndih].c[m] = NOTSET;
 953                                     }
 954                                     set_p_string(&(dih[ndih]), "");
 955                                     ndih++;
 956                                 }
 957
 958                                 nbd = nb_dist(nnb, i, l1);
 959                                 if (nbd == 3)
 960                                 {
 961                                     i1    = std::min(i, l1);
 962                                     i2    = std::max(i, l1);
 963                                     bExcl = FALSE;
 964                                     for (m = 0; m < excls[i1].nr; m++)
 965                                     {
 966                                         bExcl = bExcl || excls[i1].e[m] == i2;
 967                                     }
 968                                     if (!bExcl)
 969                                     {
 970                                         if (rtp[0].bGenerateHH14Interactions ||
 971                                             !(is_hydro(atoms, i1) && is_hydro(atoms, i2)))
 972                                         {
 973                                             if (npai == maxpai)
 974                                             {
 975                                                 maxpai += ninc;
 976                                                 srenew(pai, maxpai);
 977                                             }
 978                                             pai[npai].ai() = i1;
 979                                             pai[npai].aj() = i2;
 980                                             pai[npai].c0() = NOTSET;
 981                                             pai[npai].c1() = NOTSET;
 982                                             set_p_string(&(pai[npai]), "");
 983                                             npai++;
 984                                         }
 985                                     }
 986                                 }
 987                             }
 988                         }
 989                     }
 990                 }
 991             }
 992         }
 993     }
 994
 995     if (hb)
 996     {
 997         /* The above approach is great in that we double-check that e.g. an angle
 998          * really corresponds to three atoms connected by bonds, but this is not
 999          * generally true. Go through the angle and dihedral hackblocks to add
1000          * entries that we have not yet marked as matched when going through bonds.
1001          */
1002         for (i = 0; i < atoms->nres; i++)
1003         {
1004             /* Add remaining angles from hackblock */
1005             hbang = &hb[i].rb[ebtsANGLES];
1006             for (j = 0; j < hbang->nb; j++)
1007             {
1008                 if (hbang->b[j].match == TRUE)
1009                 {
1010                     /* We already used this entry, continue to the next */
1011                     continue;
1012                 }
1013                 /* Hm - entry not used, let's see if we can find all atoms */
1014                 if (nang == maxang)
1015                 {
1016                     maxang += ninc;
1017                     srenew(ang, maxang);
1018                 }
1019                 bFound = TRUE;
1020                 for (k = 0; k < 3 && bFound; k++)
1021                 {
1022                     p   = hbang->b[j].a[k];
1023                     res = i;
1024                     if (p[0] == '-')
1025                     {
1026                         p++;
1027                         res--;
1028                     }
1029                     else if (p[0] == '+')
1030                     {
1031                         p++;
1032                         res++;
1033                     }
1034                     ang[nang].a[k] = search_res_atom(p, res, atoms, "angle", TRUE);
1035                     bFound         = (ang[nang].a[k] != -1);
1036                 }
1037                 ang[nang].c0() = NOTSET;
1038                 ang[nang].c1() = NOTSET;
1039
1040                 if (bFound)
1041                 {
1042                     set_p_string(&(ang[nang]), hbang->b[j].s);
1043                     hbang->b[j].match = TRUE;
1044                     /* Incrementing nang means we save this angle */
1045                     nang++;
1046                 }
1047             }
1048
1049             /* Add remaining dihedrals from hackblock */
1050             hbdih = &hb[i].rb[ebtsPDIHS];
1051             for (j = 0; j < hbdih->nb; j++)
1052             {
1053                 if (hbdih->b[j].match == TRUE)
1054                 {
1055                     /* We already used this entry, continue to the next */
1056                     continue;
1057                 }
1058                 /* Hm - entry not used, let's see if we can find all atoms */
1059                 if (ndih == maxdih)
1060                 {
1061                     maxdih += ninc;
1062                     srenew(dih, maxdih);
1063                 }
1064                 bFound = TRUE;
1065                 for (k = 0; k < 4 && bFound; k++)
1066                 {
1067                     p   = hbdih->b[j].a[k];
1068                     res = i;
1069                     if (p[0] == '-')
1070                     {
1071                         p++;
1072                         res--;
1073                     }
1074                     else if (p[0] == '+')
1075                     {
1076                         p++;
1077                         res++;
1078                     }
1079                     dih[ndih].a[k] = search_res_atom(p, res, atoms, "dihedral", TRUE);
1080                     bFound         = (dih[ndih].a[k] != -1);
1081                 }
1082                 for (m = 0; m < MAXFORCEPARAM; m++)
1083                 {
1084                     dih[ndih].c[m] = NOTSET;
1085                 }
1086
1087                 if (bFound)
1088                 {
1089                     set_p_string(&(dih[ndih]), hbdih->b[j].s);
1090                     hbdih->b[j].match = TRUE;
1091                     /* Incrementing ndih means we save this dihedral */
1092                     ndih++;
1093                 }
1094             }
1095         }
1096     }
1097
1098     /* Sort angles with respect to j-i-k (middle atom first) */
1099     if (nang > 1)
1100     {
1101         qsort(ang, nang, static_cast<size_t>(sizeof(ang[0])), acomp);
1102     }
1103
1104     /* Sort dihedrals with respect to j-k-i-l (middle atoms first) */
1105     if (ndih > 1)
1106     {
1107         qsort(dih, ndih, static_cast<size_t>(sizeof(dih[0])), dcomp);
1108     }
1109
1110     /* Sort the pairs */
1111     if (npai > 1)
1112     {
1113         qsort(pai, npai, static_cast<size_t>(sizeof(pai[0])), pcomp);
1114     }
1115     if (npai > 0)
1116     {
1117         /* Remove doubles, could occur in 6-rings, such as phenyls,
1118            maybe one does not want this when fudgeQQ < 1.
1119          */
1120         fprintf(stderr, "Before cleaning: %d pairs\n", npai);
1121         rm2par(pai, &npai, preq);
1122     }
1123
1124     /* Get the impropers from the database */
1125     nimproper = get_impropers(atoms, hb, &improper, bAllowMissing);
1126
1127     /* Sort the impropers */
1128     sort_id(nimproper, improper);
1129
1130     if (ndih > 0)
1131     {
1132         fprintf(stderr, "Before cleaning: %d dihedrals\n", ndih);
1133         clean_dih(dih, &ndih, improper, nimproper, atoms,
1134                   rtp[0].bKeepAllGeneratedDihedrals,
1135                   rtp[0].bRemoveDihedralIfWithImproper);
1136     }
1137
1138     /* Now we have unique lists of angles and dihedrals
1139      * Copy them into the destination struct
1140      */
1141     cppar(ang, nang, plist, F_ANGLES);
1142     cppar(dih, ndih, plist, F_PDIHS);
1143     cppar(improper, nimproper, plist, F_IDIHS);
1144     cppar(pai, npai, plist, F_LJ14);
1145
1146     /* Remove all exclusions which are within nrexcl */
1147     clean_excls(nnb, rtp[0].nrexcl, excls);
1148
1149     sfree(ang);
1150     sfree(dih);
1151     sfree(improper);
1152     sfree(pai);
1153 }