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