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