src/gromacs/gmxpreprocess/gen_ad.cpp

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