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