src/gromacs/gmxpreprocess/gen_ad.cpp

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