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