src/gromacs/gmxpreprocess/gen_ad.c

   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, 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 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <ctype.h>
  43 #include <math.h>
  44 #include <string.h>
  45
  46 #include "sysstuff.h"
  47 #include "macros.h"
  48 #include "gromacs/utility/smalloc.h"
  49 #include "gromacs/utility/cstringutil.h"
  50 #include "gromacs/fileio/confio.h"
  51 #include "vec.h"
  52 #include "pbc.h"
  53 #include "toputil.h"
  54 #include "topio.h"
  55 #include "gpp_nextnb.h"
  56 #include "symtab.h"
  57 #include "macros.h"
  58 #include "gmx_fatal.h"
  59 #include "pgutil.h"
  60 #include "resall.h"
  61 #include "gen_ad.h"
  62
  63 #define DIHEDRAL_WAS_SET_IN_RTP 0
  64 static gmx_bool was_dihedral_set_in_rtp(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     t_param *p1, *p2;
  74     int      ac;
  75
  76     p1 = (t_param *)a1;
  77     p2 = (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     t_param *p1, *p2;
  95     int      pc;
  96
  97     p1 = (t_param *)a1;
  98     p2 = (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     t_param *p1, *p2;
 112     int      dc;
 113
 114     p1 = (t_param *)d1;
 115     p2 = (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] %u a[1] %u a[2] %u a[3] %u\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, atom_id 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 atom_id_comp(const void *a, const void *b)
 308 {
 309     return (*(atom_id *)a) - (*(atom_id *)b);
 310 }
 311
 312 static int eq_imp(atom_id a1[], atom_id a2[])
 313 {
 314     int b1[4], b2[4];
 315     int j;
 316
 317     for (j = 0; (j < 4); j++)
 318     {
 319         b1[j] = a1[j];
 320         b2[j] = a2[j];
 321     }
 322     qsort(b1, 4, (size_t)sizeof(b1[0]), int_comp);
 323     qsort(b2, 4, (size_t)sizeof(b2[0]), int_comp);
 324
 325     for (j = 0; (j < 4); j++)
 326     {
 327         if (b1[j] != b2[j])
 328         {
 329             return FALSE;
 330         }
 331     }
 332
 333     return TRUE;
 334 }
 335
 336 static int idcomp(const void *a, const void *b)
 337 {
 338     t_param *pa, *pb;
 339     int      d;
 340
 341     pa = (t_param *)a;
 342     pb = (t_param *)b;
 343     if ((d = (pa->a[0]-pb->a[0])) != 0)
 344     {
 345         return d;
 346     }
 347     else if ((d = (pa->a[3]-pb->a[3])) != 0)
 348     {
 349         return d;
 350     }
 351     else if ((d = (pa->a[1]-pb->a[1])) != 0)
 352     {
 353         return d;
 354     }
 355     else
 356     {
 357         return (int) (pa->a[2]-pb->a[2]);
 358     }
 359 }
 360
 361 static void sort_id(int nr, t_param ps[])
 362 {
 363     int i, tmp;
 364
 365     /* First swap order of atoms around if necessary */
 366     for (i = 0; (i < nr); i++)
 367     {
 368         if (ps[i].a[3] < ps[i].a[0])
 369         {
 370             tmp = ps[i].a[3]; ps[i].a[3] = ps[i].a[0]; ps[i].a[0] = tmp;
 371             tmp = ps[i].a[2]; ps[i].a[2] = ps[i].a[1]; ps[i].a[1] = tmp;
 372         }
 373     }
 374     /* Now sort it */
 375     if (nr > 1)
 376     {
 377         qsort(ps, nr, (size_t)sizeof(ps[0]), idcomp);
 378     }
 379 }
 380
 381 static int n_hydro(atom_id a[], char ***atomname)
 382 {
 383     int  i, nh = 0;
 384     char c0, c1, *aname;
 385
 386     for (i = 0; (i < 4); i += 3)
 387     {
 388         aname = *atomname[a[i]];
 389         c0    = toupper(aname[0]);
 390         if (c0 == 'H')
 391         {
 392             nh++;
 393         }
 394         else if (((int)strlen(aname) > 1) && (c0 >= '0') && (c0 <= '9'))
 395         {
 396             c1 = toupper(aname[1]);
 397             if (c1 == 'H')
 398             {
 399                 nh++;
 400             }
 401         }
 402     }
 403     return nh;
 404 }
 405
 406 /* Clean up the dihedrals (both generated and read from the .rtp
 407  * file). */
 408 static void clean_dih(t_param *dih, int *ndih, t_param improper[], int nimproper,
 409                       t_atoms *atoms, gmx_bool bKeepAllGeneratedDihedrals,
 410                       gmx_bool bRemoveDihedralIfWithImproper)
 411 {
 412     int   i, j, k, l;
 413     int  *index, nind;
 414
 415     /* Construct the list of the indices of the dihedrals
 416      * (i.e. generated or read) that might be kept. */
 417     snew(index, *ndih+1);
 418     if (bKeepAllGeneratedDihedrals)
 419     {
 420         fprintf(stderr, "Keeping all generated dihedrals\n");
 421         nind = *ndih;
 422         for (i = 0; i < nind; i++)
 423         {
 424             index[i] = i;
 425         }
 426         index[nind] = *ndih;
 427     }
 428     else
 429     {
 430         nind = 0;
 431         /* Check if generated dihedral i should be removed. The
 432          * dihedrals have been sorted by dcomp() above, so all those
 433          * on the same two central atoms are together, with those from
 434          * the .rtp file preceding those that were automatically
 435          * generated. We remove the latter if the former exist. */
 436         for (i = 0; i < *ndih; i++)
 437         {
 438             /* Keep the dihedrals that were defined in the .rtp file,
 439              * and the dihedrals that were generated and different
 440              * from the last one (whether it was generated or not). */
 441             if (was_dihedral_set_in_rtp(&dih[i]) ||
 442                 0 == i ||
 443                 !is_dihedral_on_same_bond(&dih[i], &dih[i-1]))
 444             {
 445                 index[nind++] = i;
 446             }
 447         }
 448         index[nind] = *ndih;
 449     }
 450
 451     k = 0;
 452     for (i = 0; i < nind; i++)
 453     {
 454         gmx_bool bWasSetInRTP = was_dihedral_set_in_rtp(&dih[index[i]]);
 455         gmx_bool bKeep        = TRUE;
 456         if (!bWasSetInRTP && bRemoveDihedralIfWithImproper)
 457         {
 458             /* Remove the dihedral if there is an improper on the same
 459              * bond. */
 460             for (j = 0; j < nimproper && bKeep; j++)
 461             {
 462                 bKeep = !is_dihedral_on_same_bond(&dih[index[i]], &improper[j]);
 463             }
 464         }
 465
 466         if (bKeep)
 467         {
 468             /* If we don't want all dihedrals, we want to select the
 469              * ones with the fewest hydrogens. Note that any generated
 470              * dihedrals on the same bond as an .rtp dihedral may have
 471              * been already pruned above in the construction of
 472              * index[]. However, their parameters are still present,
 473              * and l is looping over this dihedral and all of its
 474              * pruned siblings. */
 475             int bestl = index[i];
 476             if (!bKeepAllGeneratedDihedrals && !bWasSetInRTP)
 477             {
 478                 /* Minimum number of hydrogens for i and l atoms */
 479                 int minh = 2;
 480                 for (l = index[i];
 481                      (l < index[i+1] &&
 482                       is_dihedral_on_same_bond(&dih[index[i]], &dih[l]));
 483                      l++)
 484                 {
 485                     int nh = n_hydro(dih[l].a, atoms->atomname);
 486                     if (nh < minh)
 487                     {
 488                         minh  = nh;
 489                         bestl = l;
 490                     }
 491                     if (0 == minh)
 492                     {
 493                         break;
 494                     }
 495                 }
 496             }
 497             if (k != bestl)
 498             {
 499                 cpparam(&dih[k], &dih[bestl]);
 500             }
 501             k++;
 502         }
 503     }
 504
 505     for (i = k; i < *ndih; i++)
 506     {
 507         strcpy(dih[i].s, "");
 508     }
 509     *ndih = k;
 510
 511     sfree(index);
 512 }
 513
 514 static int get_impropers(t_atoms *atoms, t_hackblock hb[], t_param **improper,
 515                          gmx_bool bAllowMissing)
 516 {
 517     char         *a0;
 518     t_rbondeds   *impropers;
 519     t_rbonded    *hbimproper;
 520     int           nimproper, i, j, k, r, start, ninc, nalloc;
 521     atom_id       ai[MAXATOMLIST];
 522     gmx_bool      bStop;
 523
 524     ninc   = 500;
 525     nalloc = ninc;
 526     snew(*improper, nalloc);
 527
 528     /* Add all the impropers from the residue database to the list. */
 529     nimproper = 0;
 530     start     = 0;
 531     if (hb != NULL)
 532     {
 533         for (i = 0; (i < atoms->nres); i++)
 534         {
 535             impropers = &hb[i].rb[ebtsIDIHS];
 536             for (j = 0; (j < impropers->nb); j++)
 537             {
 538                 bStop = FALSE;
 539                 for (k = 0; (k < 4) && !bStop; k++)
 540                 {
 541                     ai[k] = search_atom(impropers->b[j].a[k], start,
 542                                         atoms,
 543                                         "improper", bAllowMissing);
 544                     if (ai[k] == NO_ATID)
 545                     {
 546                         bStop = TRUE;
 547                     }
 548                 }
 549                 if (!bStop)
 550                 {
 551                     if (nimproper == nalloc)
 552                     {
 553                         nalloc += ninc;
 554                         srenew(*improper, nalloc);
 555                     }
 556                     /* Not broken out */
 557                     set_p(&((*improper)[nimproper]), ai, NULL, impropers->b[j].s);
 558                     nimproper++;
 559                 }
 560             }
 561             while ((start < atoms->nr) && (atoms->atom[start].resind == i))
 562             {
 563                 start++;
 564             }
 565         }
 566     }
 567
 568     return nimproper;
 569 }
 570
 571 static int nb_dist(t_nextnb *nnb, int ai, int aj)
 572 {
 573     int  nre, nrx, NRE;
 574     int *nrexcl;
 575     int *a;
 576
 577     if (ai == aj)
 578     {
 579         return 0;
 580     }
 581
 582     NRE    = -1;
 583     nrexcl = nnb->nrexcl[ai];
 584     for (nre = 1; (nre < nnb->nrex); nre++)
 585     {
 586         a = nnb->a[ai][nre];
 587         for (nrx = 0; (nrx < nrexcl[nre]); nrx++)
 588         {
 589             if ((aj == a[nrx]) && (NRE == -1))
 590             {
 591                 NRE = nre;
 592             }
 593         }
 594     }
 595     return NRE;
 596 }
 597
 598 gmx_bool is_hydro(t_atoms *atoms, int ai)
 599 {
 600     return ((*(atoms->atomname[ai]))[0] == 'H');
 601 }
 602
 603 static void get_atomnames_min(int n, char **anm,
 604                               int resind, t_atoms *atoms, atom_id *a)
 605 {
 606     int m;
 607
 608     /* Assume ascending residue numbering */
 609     for (m = 0; m < n; m++)
 610     {
 611         if (atoms->atom[a[m]].resind < resind)
 612         {
 613             strcpy(anm[m], "-");
 614         }
 615         else if (atoms->atom[a[m]].resind > resind)
 616         {
 617             strcpy(anm[m], "+");
 618         }
 619         else
 620         {
 621             strcpy(anm[m], "");
 622         }
 623         strcat(anm[m], *(atoms->atomname[a[m]]));
 624     }
 625 }
 626
 627 static void gen_excls(t_atoms *atoms, t_excls *excls, t_hackblock hb[],
 628                       gmx_bool bAllowMissing)
 629 {
 630     int         r;
 631     atom_id     a, astart, i1, i2, itmp;
 632     t_rbondeds *hbexcl;
 633     int         e;
 634     char       *anm;
 635
 636     astart = 0;
 637     for (a = 0; a < atoms->nr; a++)
 638     {
 639         r = atoms->atom[a].resind;
 640         if (a == atoms->nr-1 || atoms->atom[a+1].resind != r)
 641         {
 642             hbexcl = &hb[r].rb[ebtsEXCLS];
 643
 644             for (e = 0; e < hbexcl->nb; e++)
 645             {
 646                 anm = hbexcl->b[e].a[0];
 647                 i1  = search_atom(anm, astart, atoms,
 648                                   "exclusion", bAllowMissing);
 649                 anm = hbexcl->b[e].a[1];
 650                 i2  = search_atom(anm, astart, atoms,
 651                                   "exclusion", bAllowMissing);
 652                 if (i1 != NO_ATID && i2 != NO_ATID)
 653                 {
 654                     if (i1 > i2)
 655                     {
 656                         itmp = i1;
 657                         i1   = i2;
 658                         i2   = itmp;
 659                     }
 660                     srenew(excls[i1].e, excls[i1].nr+1);
 661                     excls[i1].e[excls[i1].nr] = i2;
 662                     excls[i1].nr++;
 663                 }
 664             }
 665
 666             astart = a+1;
 667         }
 668     }
 669
 670     for (a = 0; a < atoms->nr; a++)
 671     {
 672         if (excls[a].nr > 1)
 673         {
 674             qsort(excls[a].e, excls[a].nr, (size_t)sizeof(atom_id), atom_id_comp);
 675         }
 676     }
 677 }
 678
 679 static void remove_excl(t_excls *excls, int remove)
 680 {
 681     int i;
 682
 683     for (i = remove+1; i < excls->nr; i++)
 684     {
 685         excls->e[i-1] = excls->e[i];
 686     }
 687
 688     excls->nr--;
 689 }
 690
 691 void clean_excls(t_nextnb *nnb, int nrexcl, t_excls excls[])
 692 {
 693     int      i, j, j1, k, k1, l, l1, m, n, e;
 694     t_excls *excl;
 695
 696     if (nrexcl >= 1)
 697     {
 698         /* extract all i-j-k-l neighbours from nnb struct */
 699         for (i = 0; (i < nnb->nr); i++)
 700         {
 701             /* For all particles */
 702             excl = &excls[i];
 703
 704             for (j = 0; (j < nnb->nrexcl[i][1]); j++)
 705             {
 706                 /* For all first neighbours */
 707                 j1 = nnb->a[i][1][j];
 708
 709                 for (e = 0; e < excl->nr; e++)
 710                 {
 711                     if (excl->e[e] == j1)
 712                     {
 713                         remove_excl(excl, e);
 714                     }
 715                 }
 716
 717                 if (nrexcl >= 2)
 718                 {
 719                     for (k = 0; (k < nnb->nrexcl[j1][1]); k++)
 720                     {
 721                         /* For all first neighbours of j1 */
 722                         k1 = nnb->a[j1][1][k];
 723
 724                         for (e = 0; e < excl->nr; e++)
 725                         {
 726                             if (excl->e[e] == k1)
 727                             {
 728                                 remove_excl(excl, e);
 729                             }
 730                         }
 731
 732                         if (nrexcl >= 3)
 733                         {
 734                             for (l = 0; (l < nnb->nrexcl[k1][1]); l++)
 735                             {
 736                                 /* For all first neighbours of k1 */
 737                                 l1 = nnb->a[k1][1][l];
 738
 739                                 for (e = 0; e < excl->nr; e++)
 740                                 {
 741                                     if (excl->e[e] == l1)
 742                                     {
 743                                         remove_excl(excl, e);
 744                                     }
 745                                 }
 746                             }
 747                         }
 748                     }
 749                 }
 750             }
 751         }
 752     }
 753 }
 754
 755 void generate_excls(t_nextnb *nnb, int nrexcl, t_excls excls[])
 756 {
 757     int      i, j, j1, k, k1, l, l1, m, n, e, N;
 758     t_excls *excl;
 759
 760     for (N = 1; (N < min(nrexcl, nnb->nrex)); N++)
 761     {
 762         /* extract all i-j-k-l neighbours from nnb struct */
 763         for (i = 0; (i < nnb->nr); i++)
 764         {
 765             /* For all particles */
 766             excl      = &excls[i];
 767             n         = excl->nr;
 768             excl->nr += nnb->nrexcl[i][N];
 769             srenew(excl->e, excl->nr);
 770             for (j = 0; (j < nnb->nrexcl[i][N]); j++)
 771             {
 772                 /* For all first neighbours */
 773                 if (nnb->a[i][N][j] != i)
 774                 {
 775                     excl->e[n++] = nnb->a[i][N][j];
 776                 }
 777             }
 778         }
 779     }
 780 }
 781
 782 /* Generate pairs, angles and dihedrals from .rtp settings */
 783 void gen_pad(t_nextnb *nnb, t_atoms *atoms, t_restp rtp[],
 784              t_params plist[], t_excls excls[], t_hackblock hb[],
 785              gmx_bool bAllowMissing)
 786 {
 787     t_param    *ang, *dih, *pai, *improper;
 788     t_rbondeds *hbang, *hbdih;
 789     char      **anm;
 790     const char *p;
 791     int         res, minres, maxres;
 792     int         i, j, j1, k, k1, l, l1, m, n, i1, i2;
 793     int         ninc, maxang, maxdih, maxpai;
 794     int         nang, ndih, npai, nimproper, nbd;
 795     int         nFound;
 796     gmx_bool    bFound, bExcl;
 797
 798     /* These are the angles, dihedrals and pairs that we generate
 799      * from the bonds. The ones that are already there from the rtp file
 800      * will be retained.
 801      */
 802     nang   = 0;
 803     npai   = 0;
 804     ndih   = 0;
 805     ninc   = 500;
 806     maxang = maxdih = maxpai = ninc;
 807     snew(ang, maxang);
 808     snew(dih, maxdih);
 809     snew(pai, maxpai);
 810
 811     snew(anm, 4);
 812     for (i = 0; i < 4; i++)
 813     {
 814         snew(anm[i], 12);
 815     }
 816
 817     if (hb)
 818     {
 819         gen_excls(atoms, excls, hb, bAllowMissing);
 820         /* mark all entries as not matched yet */
 821         for (i = 0; i < atoms->nres; i++)
 822         {
 823             for (j = 0; j < ebtsNR; j++)
 824             {
 825                 for (k = 0; k < hb[i].rb[j].nb; k++)
 826                 {
 827                     hb[i].rb[j].b[k].match = FALSE;
 828                 }
 829             }
 830         }
 831     }
 832
 833     /* Extract all i-j-k-l neighbours from nnb struct to generate all
 834      * angles and dihedrals. */
 835     for (i = 0; (i < nnb->nr); i++)
 836     {
 837         /* For all particles */
 838         for (j = 0; (j < nnb->nrexcl[i][1]); j++)
 839         {
 840             /* For all first neighbours */
 841             j1 = nnb->a[i][1][j];
 842             for (k = 0; (k < nnb->nrexcl[j1][1]); k++)
 843             {
 844                 /* For all first neighbours of j1 */
 845                 k1 = nnb->a[j1][1][k];
 846                 if (k1 != i)
 847                 {
 848                     /* Generate every angle only once */
 849                     if (i < k1)
 850                     {
 851                         if (nang == maxang)
 852                         {
 853                             maxang += ninc;
 854                             srenew(ang, maxang);
 855                         }
 856                         ang[nang].AI = i;
 857                         ang[nang].AJ = j1;
 858                         ang[nang].AK = k1;
 859                         ang[nang].C0 = NOTSET;
 860                         ang[nang].C1 = NOTSET;
 861                         set_p_string(&(ang[nang]), "");
 862                         if (hb)
 863                         {
 864                             minres = atoms->atom[ang[nang].a[0]].resind;
 865                             maxres = minres;
 866                             for (m = 1; m < 3; m++)
 867                             {
 868                                 minres = min(minres, atoms->atom[ang[nang].a[m]].resind);
 869                                 maxres = max(maxres, atoms->atom[ang[nang].a[m]].resind);
 870                             }
 871                             res = 2*minres-maxres;
 872                             do
 873                             {
 874                                 res += maxres-minres;
 875                                 get_atomnames_min(3, anm, res, atoms, ang[nang].a);
 876                                 hbang = &hb[res].rb[ebtsANGLES];
 877                                 for (l = 0; (l < hbang->nb); l++)
 878                                 {
 879                                     if (strcmp(anm[1], hbang->b[l].AJ) == 0)
 880                                     {
 881                                         bFound = FALSE;
 882                                         for (m = 0; m < 3; m += 2)
 883                                         {
 884                                             bFound = (bFound ||
 885                                                       ((strcmp(anm[m], hbang->b[l].AI) == 0) &&
 886                                                        (strcmp(anm[2-m], hbang->b[l].AK) == 0)));
 887                                         }
 888                                         if (bFound)
 889                                         {
 890                                             set_p_string(&(ang[nang]), hbang->b[l].s);
 891                                             /* Mark that we found a match for this entry */
 892                                             hbang->b[l].match = TRUE;
 893                                         }
 894                                     }
 895                                 }
 896                             }
 897                             while (res < maxres);
 898                         }
 899                         nang++;
 900                     }
 901                     /* Generate every dihedral, 1-4 exclusion and 1-4 interaction
 902                        only once */
 903                     if (j1 < k1)
 904                     {
 905                         for (l = 0; (l < nnb->nrexcl[k1][1]); l++)
 906                         {
 907                             /* For all first neighbours of k1 */
 908                             l1 = nnb->a[k1][1][l];
 909                             if ((l1 != i) && (l1 != j1))
 910                             {
 911                                 if (ndih == maxdih)
 912                                 {
 913                                     maxdih += ninc;
 914                                     srenew(dih, maxdih);
 915                                 }
 916                                 dih[ndih].AI = i;
 917                                 dih[ndih].AJ = j1;
 918                                 dih[ndih].AK = k1;
 919                                 dih[ndih].AL = l1;
 920                                 for (m = 0; m < MAXFORCEPARAM; m++)
 921                                 {
 922                                     dih[ndih].c[m] = NOTSET;
 923                                 }
 924                                 set_p_string(&(dih[ndih]), "");
 925                                 nFound = 0;
 926                                 if (hb)
 927                                 {
 928                                     minres = atoms->atom[dih[ndih].a[0]].resind;
 929                                     maxres = minres;
 930                                     for (m = 1; m < 4; m++)
 931                                     {
 932                                         minres = min(minres, atoms->atom[dih[ndih].a[m]].resind);
 933                                         maxres = max(maxres, atoms->atom[dih[ndih].a[m]].resind);
 934                                     }
 935                                     res = 2*minres-maxres;
 936                                     do
 937                                     {
 938                                         res += maxres-minres;
 939                                         get_atomnames_min(4, anm, res, atoms, dih[ndih].a);
 940                                         hbdih = &hb[res].rb[ebtsPDIHS];
 941                                         for (n = 0; (n < hbdih->nb); n++)
 942                                         {
 943                                             bFound = FALSE;
 944                                             for (m = 0; m < 2; m++)
 945                                             {
 946                                                 bFound = (bFound ||
 947                                                           ((strcmp(anm[3*m],  hbdih->b[n].AI) == 0) &&
 948                                                            (strcmp(anm[1+m],  hbdih->b[n].AJ) == 0) &&
 949                                                            (strcmp(anm[2-m],  hbdih->b[n].AK) == 0) &&
 950                                                            (strcmp(anm[3-3*m], hbdih->b[n].AL) == 0)));
 951                                             }
 952                                             if (bFound)
 953                                             {
 954                                                 set_p_string(&dih[ndih], hbdih->b[n].s);
 955                                                 /* Mark that we found a match for this entry */
 956                                                 hbdih->b[n].match = TRUE;
 957
 958                                                 /* Set the last parameter to be able to see
 959                                                    if the dihedral was in the rtp list.
 960                                                  */
 961                                                 dih[ndih].c[MAXFORCEPARAM-1] = DIHEDRAL_WAS_SET_IN_RTP;
 962                                                 nFound++;
 963                                                 ndih++;
 964                                                 /* Set the next direct in case the rtp contains
 965                                                    multiple entries for this dihedral.
 966                                                  */
 967                                                 if (ndih == maxdih)
 968                                                 {
 969                                                     maxdih += ninc;
 970                                                     srenew(dih, maxdih);
 971                                                 }
 972                                                 dih[ndih].AI = i;
 973                                                 dih[ndih].AJ = j1;
 974                                                 dih[ndih].AK = k1;
 975                                                 dih[ndih].AL = l1;
 976                                                 for (m = 0; m < MAXFORCEPARAM; m++)
 977                                                 {
 978                                                     dih[ndih].c[m] = NOTSET;
 979                                                 }
 980                                             }
 981                                         }
 982                                     }
 983                                     while (res < maxres);
 984                                 }
 985                                 if (nFound == 0)
 986                                 {
 987                                     if (ndih == maxdih)
 988                                     {
 989                                         maxdih += ninc;
 990                                         srenew(dih, maxdih);
 991                                     }
 992                                     dih[ndih].AI = i;
 993                                     dih[ndih].AJ = j1;
 994                                     dih[ndih].AK = k1;
 995                                     dih[ndih].AL = l1;
 996                                     for (m = 0; m < MAXFORCEPARAM; m++)
 997                                     {
 998                                         dih[ndih].c[m] = NOTSET;
 999                                     }
1000                                     set_p_string(&(dih[ndih]), "");
1001                                     ndih++;
1002                                 }
1003
1004                                 nbd = nb_dist(nnb, i, l1);
1005                                 if (debug)
1006                                 {
1007                                     fprintf(debug, "Distance (%d-%d) = %d\n", i+1, l1+1, nbd);
1008                                 }
1009                                 if (nbd == 3)
1010                                 {
1011                                     i1    = min(i, l1);
1012                                     i2    = max(i, l1);
1013                                     bExcl = FALSE;
1014                                     for (m = 0; m < excls[i1].nr; m++)
1015                                     {
1016                                         bExcl = bExcl || excls[i1].e[m] == i2;
1017                                     }
1018                                     if (!bExcl)
1019                                     {
1020                                         if (rtp[0].bGenerateHH14Interactions ||
1021                                             !(is_hydro(atoms, i1) && is_hydro(atoms, i2)))
1022                                         {
1023                                             if (npai == maxpai)
1024                                             {
1025                                                 maxpai += ninc;
1026                                                 srenew(pai, maxpai);
1027                                             }
1028                                             pai[npai].AI = i1;
1029                                             pai[npai].AJ = i2;
1030                                             pai[npai].C0 = NOTSET;
1031                                             pai[npai].C1 = NOTSET;
1032                                             set_p_string(&(pai[npai]), "");
1033                                             npai++;
1034                                         }
1035                                     }
1036                                 }
1037                             }
1038                         }
1039                     }
1040                 }
1041             }
1042         }
1043     }
1044
1045     /* The above approach is great in that we double-check that e.g. an angle
1046      * really corresponds to three atoms connected by bonds, but this is not
1047      * generally true. Go through the angle and dihedral hackblocks to add
1048      * entries that we have not yet marked as matched when going through bonds.
1049      */
1050     for (i = 0; i < atoms->nres; i++)
1051     {
1052         /* Add remaining angles from hackblock */
1053         hbang = &hb[i].rb[ebtsANGLES];
1054         for (j = 0; j < hbang->nb; j++)
1055         {
1056             if (hbang->b[j].match == TRUE)
1057             {
1058                 /* We already used this entry, continue to the next */
1059                 continue;
1060             }
1061             /* Hm - entry not used, let's see if we can find all atoms */
1062             if (nang == maxang)
1063             {
1064                 maxang += ninc;
1065                 srenew(ang, maxang);
1066             }
1067             bFound = TRUE;
1068             for (k = 0; k < 3 && bFound; k++)
1069             {
1070                 p   = hbang->b[j].a[k];
1071                 res = i;
1072                 if (p[0] == '-')
1073                 {
1074                     p++;
1075                     res--;
1076                 }
1077                 else if (p[0] == '+')
1078                 {
1079                     p++;
1080                     res++;
1081                 }
1082                 ang[nang].a[k] = search_res_atom(p, res, atoms, "angle", TRUE);
1083                 bFound         = (ang[nang].a[k] != NO_ATID);
1084             }
1085             ang[nang].C0 = NOTSET;
1086             ang[nang].C1 = NOTSET;
1087
1088             if (bFound)
1089             {
1090                 set_p_string(&(ang[nang]), hbang->b[j].s);
1091                 hbang->b[j].match = TRUE;
1092                 /* Incrementing nang means we save this angle */
1093                 nang++;
1094             }
1095         }
1096
1097         /* Add remaining dihedrals from hackblock */
1098         hbdih = &hb[i].rb[ebtsPDIHS];
1099         for (j = 0; j < hbdih->nb; j++)
1100         {
1101             if (hbdih->b[j].match == TRUE)
1102             {
1103                 /* We already used this entry, continue to the next */
1104                 continue;
1105             }
1106             /* Hm - entry not used, let's see if we can find all atoms */
1107             if (ndih == maxdih)
1108             {
1109                 maxdih += ninc;
1110                 srenew(dih, maxdih);
1111             }
1112             bFound = TRUE;
1113             for (k = 0; k < 4 && bFound; k++)
1114             {
1115                 p   = hbdih->b[j].a[k];
1116                 res = i;
1117                 if (p[0] == '-')
1118                 {
1119                     p++;
1120                     res--;
1121                 }
1122                 else if (p[0] == '+')
1123                 {
1124                     p++;
1125                     res++;
1126                 }
1127                 dih[ndih].a[k] = search_res_atom(p, res, atoms, "dihedral", TRUE);
1128                 bFound         = (dih[ndih].a[k] != NO_ATID);
1129             }
1130             for (m = 0; m < MAXFORCEPARAM; m++)
1131             {
1132                 dih[ndih].c[m] = NOTSET;
1133             }
1134
1135             if (bFound)
1136             {
1137                 set_p_string(&(dih[ndih]), hbdih->b[j].s);
1138                 hbdih->b[j].match = TRUE;
1139                 /* Incrementing ndih means we save this dihedral */
1140                 ndih++;
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 }