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