src/gromacs/gmxpreprocess/gen_ad.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  34  */
  35 /* This file is completely threadsafe - keep it that way! */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <math.h>
  41 #include <ctype.h>
  42 #include "sysstuff.h"
  43 #include "macros.h"
  44 #include "smalloc.h"
  45 #include "string2.h"
  46 #include "confio.h"
  47 #include "vec.h"
  48 #include "pbc.h"
  49 #include "toputil.h"
  50 #include "topio.h"
  51 #include "gpp_nextnb.h"
  52 #include "symtab.h"
  53 #include "macros.h"
  54 #include "gmx_fatal.h"
  55 #include "pgutil.h"
  56 #include "resall.h"
  57 #include "gen_ad.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     return ac;
  76   else if ((ac=(p1->AI-p2->AI))!=0)
  77     return ac;
  78   else
  79     return (p1->AK-p2->AK);
  80 }
  81
  82 static int pcomp(const void *a1, const void *a2)
  83 {
  84   t_param *p1,*p2;
  85   int     pc;
  86
  87   p1=(t_param *)a1;
  88   p2=(t_param *)a2;
  89   if ((pc=(p1->AI-p2->AI))!=0)
  90     return pc;
  91   else
  92     return (p1->AJ-p2->AJ);
  93 }
  94
  95 static int dcomp(const void *d1, const void *d2)
  96 {
  97     t_param *p1,*p2;
  98     int     dc;
  99
 100     p1=(t_param *)d1;
 101     p2=(t_param *)d2;
 102     /* First sort by J & K (the two central) atoms */
 103     if ((dc=(p1->AJ-p2->AJ))!=0)
 104     {
 105         return dc;
 106     }
 107     else if ((dc=(p1->AK-p2->AK))!=0)
 108     {
 109         return dc;
 110     }
 111     /* Then make sure to put rtp dihedrals before generated ones */
 112     else if (was_dihedral_set_in_rtp(p1) &&
 113              !was_dihedral_set_in_rtp(p2))
 114     {
 115         return -1;
 116     }
 117     else if (!was_dihedral_set_in_rtp(p1) &&
 118              was_dihedral_set_in_rtp(p2))
 119     {
 120         return 1;
 121     }
 122     /* Finally, sort by I and J (two outer) atoms */
 123     else if ((dc=(p1->AI-p2->AI))!=0)
 124     {
 125         return dc;
 126     }
 127     else
 128     {
 129         return (p1->AL-p2->AL);
 130     }
 131 }
 132
 133
 134 static gmx_bool is_dihedral_on_same_bond(t_param *p1, t_param *p2)
 135 {
 136   if (((p1->AJ==p2->AJ) && (p1->AK==p2->AK)) ||
 137       ((p1->AJ==p2->AK) && (p1->AK==p2->AJ)))
 138     return TRUE;
 139   else
 140     return FALSE;
 141 }
 142
 143
 144 static gmx_bool preq(t_param *p1, t_param *p2)
 145 {
 146   if ((p1->AI==p2->AI) && (p1->AJ==p2->AJ))
 147     return TRUE;
 148   else
 149     return FALSE;
 150 }
 151
 152 static void rm2par(t_param p[], int *np, peq eq)
 153 {
 154   int *index,nind;
 155   int i,j;
 156
 157   if ((*np)==0)
 158     return;
 159
 160   snew(index,*np);
 161   nind=0;
 162     index[nind++]=0;
 163   for(i=1; (i<(*np)); i++)
 164     if (!eq(&p[i],&p[i-1]))
 165       index[nind++]=i;
 166   /* Index now holds pointers to all the non-equal params,
 167    * this only works when p is sorted of course
 168    */
 169   for(i=0; (i<nind); i++) {
 170     for(j=0; (j<MAXATOMLIST); j++)
 171       p[i].a[j]=p[index[i]].a[j];
 172     for(j=0; (j<MAXFORCEPARAM); j++)
 173       p[i].c[j]=p[index[i]].c[j];
 174     if (p[index[i]].a[0] == p[index[i]].a[1]) {
 175       if (debug)
 176         fprintf(debug,
 177                 "Something VERY strange is going on in rm2par (gen_ad.c)\n"
 178                 "a[0] %u a[1] %u a[2] %u a[3] %u\n",
 179                 p[i].a[0],p[i].a[1],p[i].a[2],p[i].a[3]);
 180       strcpy(p[i].s,"");
 181     } else if (index[i] > i) {
 182       /* Copy the string only if it comes from somewhere else
 183        * otherwise we will end up copying a random (newly freed) pointer.
 184        * Since the index is sorted we only have to test for index[i] > i.
 185        */
 186       strcpy(p[i].s,p[index[i]].s);
 187     }
 188   }
 189   (*np)=nind;
 190
 191   sfree(index);
 192 }
 193
 194 static void cppar(t_param p[], int np, t_params plist[], int ftype)
 195 {
 196   int      i,j,nral,nrfp;
 197   t_params *ps;
 198
 199   ps   = &plist[ftype];
 200   nral = NRAL(ftype);
 201   nrfp = NRFP(ftype);
 202
 203   /* Keep old stuff */
 204   pr_alloc(np,ps);
 205   for(i=0; (i<np); i++) {
 206     for(j=0; (j<nral); j++)
 207       ps->param[ps->nr].a[j] = p[i].a[j];
 208     for(j=0; (j<nrfp); j++)
 209       ps->param[ps->nr].c[j] = p[i].c[j];
 210     for(j=0; (j<MAXSLEN); j++)
 211       ps->param[ps->nr].s[j] = p[i].s[j];
 212     ps->nr++;
 213   }
 214 }
 215
 216 static void cpparam(t_param *dest, t_param *src)
 217 {
 218   int j;
 219
 220   for(j=0; (j<MAXATOMLIST); j++)
 221     dest->a[j] = src->a[j];
 222   for(j=0; (j<MAXFORCEPARAM); j++)
 223     dest->c[j] = src->c[j];
 224   for(j=0; (j<MAXSLEN); j++)
 225     dest->s[j] = src->s[j];
 226 }
 227
 228 static void set_p(t_param *p, atom_id ai[4], real *c, char *s)
 229 {
 230   int j;
 231
 232   for(j=0; (j<4); j++)
 233     p->a[j]=ai[j];
 234   for(j=0; (j<MAXFORCEPARAM); j++)
 235     if (c)
 236       p->c[j]=c[j];
 237     else
 238       p->c[j]=NOTSET;
 239
 240   set_p_string(p,s);
 241 }
 242
 243 static int int_comp(const void *a,const void *b)
 244 {
 245   return (*(int *)a) - (*(int *)b);
 246 }
 247
 248 static int atom_id_comp(const void *a,const void *b)
 249 {
 250   return (*(atom_id *)a) - (*(atom_id *)b);
 251 }
 252
 253 static int eq_imp(atom_id a1[],atom_id a2[])
 254 {
 255   int b1[4],b2[4];
 256   int j;
 257
 258   for(j=0; (j<4); j++) {
 259     b1[j]=a1[j];
 260     b2[j]=a2[j];
 261   }
 262   qsort(b1,4,(size_t)sizeof(b1[0]),int_comp);
 263   qsort(b2,4,(size_t)sizeof(b2[0]),int_comp);
 264
 265   for(j=0; (j<4); j++)
 266     if (b1[j] != b2[j])
 267       return FALSE;
 268
 269   return TRUE;
 270 }
 271
 272 static int idcomp(const void *a,const void *b)
 273 {
 274   t_param *pa,*pb;
 275   int     d;
 276
 277   pa=(t_param *)a;
 278   pb=(t_param *)b;
 279   if ((d=(pa->a[0]-pb->a[0])) != 0)
 280     return d;
 281   else if ((d=(pa->a[3]-pb->a[3])) != 0)
 282     return d;
 283   else if ((d=(pa->a[1]-pb->a[1])) != 0)
 284     return d;
 285   else
 286     return (int) (pa->a[2]-pb->a[2]);
 287 }
 288
 289 static void sort_id(int nr,t_param ps[])
 290 {
 291   int i,tmp;
 292
 293   /* First swap order of atoms around if necessary */
 294   for(i=0; (i<nr); i++) {
 295     if (ps[i].a[3] < ps[i].a[0]) {
 296       tmp = ps[i].a[3]; ps[i].a[3] = ps[i].a[0]; ps[i].a[0] = tmp;
 297       tmp = ps[i].a[2]; ps[i].a[2] = ps[i].a[1]; ps[i].a[1] = tmp;
 298     }
 299   }
 300   /* Now sort it */
 301   if (nr > 1)
 302     qsort(ps,nr,(size_t)sizeof(ps[0]),idcomp);
 303 }
 304
 305 static int n_hydro(atom_id a[],char ***atomname)
 306 {
 307   int i,nh=0;
 308   char c0,c1,*aname;
 309
 310   for(i=0; (i<4); i+=3) {
 311     aname=*atomname[a[i]];
 312     c0=toupper(aname[0]);
 313     if (c0 == 'H')
 314       nh++;
 315     else if (((int)strlen(aname) > 1) && (c0 >= '0') && (c0 <= '9')) {
 316       c1=toupper(aname[1]);
 317       if (c1 == 'H')
 318         nh++;
 319     }
 320   }
 321   return nh;
 322 }
 323
 324 /* Clean up the dihedrals (both generated and read from the .rtp
 325  * file). */
 326 static void clean_dih(t_param *dih, int *ndih,t_param improper[],int nimproper,
 327                       t_atoms *atoms,gmx_bool bKeepAllGeneratedDihedrals,
 328                       gmx_bool bRemoveDihedralIfWithImproper)
 329 {
 330     int  i,j,k,l;
 331     int  *index,nind;
 332
 333     /* Construct the list of the indices of the dihedrals
 334      * (i.e. generated or read) that might be kept. */
 335     snew(index, *ndih+1);
 336     if (bKeepAllGeneratedDihedrals)
 337     {
 338         fprintf(stderr,"Keeping all generated dihedrals\n");
 339         nind = *ndih;
 340         for(i = 0; i < nind; i++)
 341         {
 342             index[i] = i;
 343         }
 344         index[nind] = *ndih;
 345     }
 346     else
 347     {
 348         nind = 0;
 349         /* Check if generated dihedral i should be removed. The
 350          * dihedrals have been sorted by dcomp() above, so all those
 351          * on the same two central atoms are together, with those from
 352          * the .rtp file preceding those that were automatically
 353          * generated. We remove the latter if the former exist. */
 354         for(i = 0; i < *ndih; i++)
 355         {
 356             /* Keep the dihedrals that were defined in the .rtp file,
 357              * and the dihedrals that were generated and different
 358              * from the last one (whether it was generated or not). */
 359             if (was_dihedral_set_in_rtp(&dih[i]) ||
 360                 0 == i ||
 361                 !is_dihedral_on_same_bond(&dih[i],&dih[i-1]))
 362             {
 363                 index[nind++] = i;
 364             }
 365         }
 366         index[nind] = *ndih;
 367     }
 368
 369     k=0;
 370     for(i=0; i<nind; i++)
 371     {
 372         gmx_bool bWasSetInRTP = was_dihedral_set_in_rtp(&dih[index[i]]);
 373         gmx_bool bKeep = TRUE;
 374         if (!bWasSetInRTP && bRemoveDihedralIfWithImproper)
 375         {
 376             /* Remove the dihedral if there is an improper on the same
 377              * bond. */
 378             for(j = 0; j < nimproper && bKeep; j++)
 379             {
 380                 bKeep = !is_dihedral_on_same_bond(&dih[index[i]],&improper[j]);
 381             }
 382         }
 383
 384         if (bKeep)
 385         {
 386             /* If we don't want all dihedrals, we want to select the
 387              * ones with the fewest hydrogens. Note that any generated
 388              * dihedrals on the same bond as an .rtp dihedral may have
 389              * been already pruned above in the construction of
 390              * index[]. However, their parameters are still present,
 391              * and l is looping over this dihedral and all of its
 392              * pruned siblings. */
 393             int bestl = index[i];
 394             if (!bKeepAllGeneratedDihedrals && !bWasSetInRTP)
 395             {
 396                 /* Minimum number of hydrogens for i and l atoms */
 397                 int minh = 2;
 398                 for(l = index[i];
 399                     (l < index[i+1] &&
 400                      is_dihedral_on_same_bond(&dih[index[i]],&dih[l]));
 401                     l++)
 402                 {
 403                     int nh = n_hydro(dih[l].a,atoms->atomname);
 404                     if (nh < minh)
 405                     {
 406                         minh=nh;
 407                         bestl=l;
 408                     }
 409                     if (0 == minh)
 410                     {
 411                         break;
 412                     }
 413                 }
 414             }
 415             if (k != bestl)
 416             {
 417                 cpparam(&dih[k],&dih[bestl]);
 418             }
 419             k++;
 420         }
 421     }
 422
 423     for (i = k; i < *ndih; i++)
 424     {
 425         strcpy(dih[i].s,"");
 426     }
 427     *ndih = k;
 428
 429     sfree(index);
 430 }
 431
 432 static int get_impropers(t_atoms *atoms,t_hackblock hb[],t_param **improper,
 433                          gmx_bool bAllowMissing)
 434 {
 435   char      *a0;
 436   t_rbondeds *impropers;
 437   t_rbonded  *hbimproper;
 438   int       nimproper,i,j,k,r,start,ninc,nalloc;
 439   atom_id   ai[MAXATOMLIST];
 440   gmx_bool      bStop;
 441
 442   ninc = 500;
 443   nalloc = ninc;
 444   snew(*improper,nalloc);
 445
 446   /* Add all the impropers from the residue database to the list. */
 447   nimproper = 0;
 448   start = 0;
 449   if (hb != NULL) {
 450     for(i=0; (i<atoms->nres); i++) {
 451       impropers=&hb[i].rb[ebtsIDIHS];
 452       for(j=0; (j<impropers->nb); j++) {
 453         bStop=FALSE;
 454         for(k=0; (k<4) && !bStop; k++) {
 455           ai[k] = search_atom(impropers->b[j].a[k],start,
 456                           atoms,
 457                               "improper",bAllowMissing);
 458           if (ai[k] == NO_ATID)
 459             bStop = TRUE;
 460         }
 461         if (!bStop) {
 462           if (nimproper == nalloc) {
 463             nalloc += ninc;
 464             srenew(*improper,nalloc);
 465           }
 466           /* Not broken out */
 467           set_p(&((*improper)[nimproper]),ai,NULL,impropers->b[j].s);
 468           nimproper++;
 469         }
 470       }
 471       while ((start<atoms->nr) && (atoms->atom[start].resind == i))
 472         start++;
 473     }
 474   }
 475
 476   return nimproper;
 477 }
 478
 479 static int nb_dist(t_nextnb *nnb,int ai,int aj)
 480 {
 481   int nre,nrx,NRE;
 482   int *nrexcl;
 483   int *a;
 484
 485   if (ai == aj)
 486     return 0;
 487
 488   NRE=-1;
 489   nrexcl=nnb->nrexcl[ai];
 490   for(nre=1; (nre < nnb->nrex); nre++) {
 491     a=nnb->a[ai][nre];
 492     for(nrx=0; (nrx < nrexcl[nre]); nrx++) {
 493       if ((aj == a[nrx]) && (NRE == -1))
 494         NRE=nre;
 495     }
 496   }
 497   return NRE;
 498 }
 499
 500 gmx_bool is_hydro(t_atoms *atoms,int ai)
 501 {
 502   return ((*(atoms->atomname[ai]))[0] == 'H');
 503 }
 504
 505 static void get_atomnames_min(int n,char **anm,
 506                               int resind,t_atoms *atoms,atom_id *a)
 507 {
 508   int m;
 509
 510   /* Assume ascending residue numbering */
 511   for(m=0; m<n; m++) {
 512     if (atoms->atom[a[m]].resind < resind)
 513       strcpy(anm[m],"-");
 514     else if (atoms->atom[a[m]].resind > resind)
 515       strcpy(anm[m],"+");
 516     else
 517       strcpy(anm[m],"");
 518     strcat(anm[m],*(atoms->atomname[a[m]]));
 519   }
 520 }
 521
 522 static void gen_excls(t_atoms *atoms, t_excls *excls, t_hackblock hb[],
 523                       gmx_bool bAllowMissing)
 524 {
 525   int        r;
 526   atom_id    a,astart,i1,i2,itmp;
 527   t_rbondeds *hbexcl;
 528   int        e;
 529   char       *anm;
 530
 531   astart = 0;
 532   for(a=0; a<atoms->nr; a++) {
 533     r = atoms->atom[a].resind;
 534     if (a==atoms->nr-1 || atoms->atom[a+1].resind != r) {
 535       hbexcl = &hb[r].rb[ebtsEXCLS];
 536
 537       for(e=0; e<hbexcl->nb; e++) {
 538         anm = hbexcl->b[e].a[0];
 539         i1 = search_atom(anm,astart,atoms,
 540                          "exclusion",bAllowMissing);
 541         anm = hbexcl->b[e].a[1];
 542         i2 = search_atom(anm,astart,atoms,
 543                          "exclusion",bAllowMissing);
 544         if (i1!=NO_ATID && i2!=NO_ATID) {
 545           if (i1 > i2) {
 546             itmp = i1;
 547             i1 = i2;
 548             i2 = itmp;
 549           }
 550           srenew(excls[i1].e,excls[i1].nr+1);
 551           excls[i1].e[excls[i1].nr] = i2;
 552           excls[i1].nr++;
 553         }
 554       }
 555
 556       astart = a+1;
 557     }
 558   }
 559
 560   for(a=0; a<atoms->nr; a++)
 561     if (excls[a].nr > 1)
 562       qsort(excls[a].e,excls[a].nr,(size_t)sizeof(atom_id),atom_id_comp);
 563 }
 564
 565 static void remove_excl(t_excls *excls, int remove)
 566 {
 567   int i;
 568
 569   for(i=remove+1; i<excls->nr; i++)
 570     excls->e[i-1] = excls->e[i];
 571
 572   excls->nr--;
 573 }
 574
 575 void clean_excls(t_nextnb *nnb, int nrexcl, t_excls excls[])
 576 {
 577   int i,j,j1,k,k1,l,l1,m,n,e;
 578   t_excls *excl;
 579
 580   if (nrexcl >= 1)
 581     /* extract all i-j-k-l neighbours from nnb struct */
 582     for(i=0; (i<nnb->nr); i++) {
 583       /* For all particles */
 584       excl = &excls[i];
 585
 586       for(j=0; (j<nnb->nrexcl[i][1]); j++) {
 587         /* For all first neighbours */
 588         j1=nnb->a[i][1][j];
 589
 590         for(e=0; e<excl->nr; e++)
 591           if (excl->e[e] == j1)
 592             remove_excl(excl,e);
 593
 594         if (nrexcl >= 2)
 595           for(k=0; (k<nnb->nrexcl[j1][1]); k++) {
 596             /* For all first neighbours of j1 */
 597             k1=nnb->a[j1][1][k];
 598
 599             for(e=0; e<excl->nr; e++)
 600               if (excl->e[e] == k1)
 601                 remove_excl(excl,e);
 602
 603             if (nrexcl >= 3)
 604               for(l=0; (l<nnb->nrexcl[k1][1]); l++) {
 605                 /* For all first neighbours of k1 */
 606                 l1=nnb->a[k1][1][l];
 607
 608                 for(e=0; e<excl->nr; e++)
 609                   if (excl->e[e] == l1)
 610                     remove_excl(excl,e);
 611               }
 612           }
 613       }
 614     }
 615 }
 616
 617 void generate_excls(t_nextnb *nnb, int nrexcl, t_excls excls[])
 618 {
 619   int i,j,j1,k,k1,l,l1,m,n,e,N;
 620   t_excls *excl;
 621
 622   for(N=1; (N<min(nrexcl,nnb->nrex)); N++) {
 623     /* extract all i-j-k-l neighbours from nnb struct */
 624     for(i=0; (i<nnb->nr); i++) {
 625       /* For all particles */
 626       excl = &excls[i];
 627       n = excl->nr;
 628       excl->nr += nnb->nrexcl[i][N];
 629       srenew(excl->e,excl->nr);
 630       for(j=0; (j<nnb->nrexcl[i][N]); j++)
 631         /* For all first neighbours */
 632         if (nnb->a[i][N][j] != i)
 633           excl->e[n++] = nnb->a[i][N][j];
 634     }
 635   }
 636 }
 637
 638 /* Generate pairs, angles and dihedrals from .rtp settings */
 639 void gen_pad(t_nextnb *nnb, t_atoms *atoms, t_restp rtp[],
 640              t_params plist[], t_excls excls[], t_hackblock hb[],
 641              gmx_bool bAllowMissing)
 642 {
 643   t_param *ang,*dih,*pai,*improper;
 644   t_rbondeds *hbang, *hbdih;
 645   char    **anm;
 646   int     res,minres,maxres;
 647   int     i,j,j1,k,k1,l,l1,m,n,i1,i2;
 648   int     ninc,maxang,maxdih,maxpai;
 649   int     nang,ndih,npai,nimproper,nbd;
 650   int     nFound;
 651   gmx_bool    bFound,bExcl;
 652
 653
 654   /* These are the angles, dihedrals and pairs that we generate
 655    * from the bonds. The ones that are already there from the rtp file
 656    * will be retained.
 657    */
 658   nang   = 0;
 659   npai   = 0;
 660   ndih   = 0;
 661   ninc   = 500;
 662   maxang = maxdih = maxpai = ninc;
 663   snew(ang, maxang);
 664   snew(dih, maxdih);
 665   snew(pai, maxpai);
 666
 667   snew(anm,4);
 668   for(i=0;i<4;i++)
 669     snew(anm[i],12);
 670
 671   if (hb)
 672     gen_excls(atoms,excls,hb,bAllowMissing);
 673
 674   /* Extract all i-j-k-l neighbours from nnb struct to generate all
 675    * angles and dihedrals. */
 676   for(i=0; (i<nnb->nr); i++)
 677     /* For all particles */
 678     for(j=0; (j<nnb->nrexcl[i][1]); j++) {
 679       /* For all first neighbours */
 680       j1=nnb->a[i][1][j];
 681       for(k=0; (k<nnb->nrexcl[j1][1]); k++) {
 682         /* For all first neighbours of j1 */
 683         k1=nnb->a[j1][1][k];
 684         if (k1 != i) {
 685           /* Generate every angle only once */
 686           if (i < k1) {
 687             if (nang == maxang) {
 688               maxang += ninc;
 689               srenew(ang,maxang);
 690             }
 691             ang[nang].AI=i;
 692             ang[nang].AJ=j1;
 693             ang[nang].AK=k1;
 694             ang[nang].C0=NOTSET;
 695             ang[nang].C1=NOTSET;
 696             set_p_string(&(ang[nang]),"");
 697             if (hb) {
 698               minres = atoms->atom[ang[nang].a[0]].resind;
 699               maxres = minres;
 700               for(m=1; m<3; m++) {
 701                 minres = min(minres,atoms->atom[ang[nang].a[m]].resind);
 702                 maxres = max(maxres,atoms->atom[ang[nang].a[m]].resind);
 703               }
 704               res = 2*minres-maxres;
 705               do {
 706                 res += maxres-minres;
 707                 get_atomnames_min(3,anm,res,atoms,ang[nang].a);
 708                 hbang=&hb[res].rb[ebtsANGLES];
 709                 for(l=0; (l<hbang->nb); l++) {
 710                   if (strcmp(anm[1],hbang->b[l].AJ)==0) {
 711                     bFound=FALSE;
 712                     for (m=0; m<3; m+=2)
 713                       bFound=(bFound ||
 714                               ((strcmp(anm[m],hbang->b[l].AI)==0) &&
 715                                (strcmp(anm[2-m],hbang->b[l].AK)==0)));
 716                     if (bFound) {
 717                       set_p_string(&(ang[nang]),hbang->b[l].s);
 718                     }
 719                   }
 720                 }
 721               } while (res < maxres);
 722             }
 723             nang++;
 724           }
 725           /* Generate every dihedral, 1-4 exclusion and 1-4 interaction
 726              only once */
 727           if (j1 < k1) {
 728             for(l=0; (l<nnb->nrexcl[k1][1]); l++) {
 729               /* For all first neighbours of k1 */
 730               l1=nnb->a[k1][1][l];
 731               if ((l1 != i) && (l1 != j1)) {
 732                 if (ndih == maxdih) {
 733                   maxdih += ninc;
 734                   srenew(dih,maxdih);
 735                 }
 736                 dih[ndih].AI=i;
 737                 dih[ndih].AJ=j1;
 738                 dih[ndih].AK=k1;
 739                 dih[ndih].AL=l1;
 740                 for (m=0; m<MAXFORCEPARAM; m++)
 741                   dih[ndih].c[m]=NOTSET;
 742                 set_p_string(&(dih[ndih]),"");
 743                 nFound = 0;
 744                 if (hb) {
 745                   minres = atoms->atom[dih[ndih].a[0]].resind;
 746                   maxres = minres;
 747                   for(m=1; m<4; m++) {
 748                     minres = min(minres,atoms->atom[dih[ndih].a[m]].resind);
 749                     maxres = max(maxres,atoms->atom[dih[ndih].a[m]].resind);
 750                   }
 751                   res = 2*minres-maxres;
 752                   do {
 753                     res += maxres-minres;
 754                     get_atomnames_min(4,anm,res,atoms,dih[ndih].a);
 755                     hbdih=&hb[res].rb[ebtsPDIHS];
 756                     for(n=0; (n<hbdih->nb); n++) {
 757                       bFound=FALSE;
 758                       for (m=0; m<2; m++)
 759                         bFound=(bFound ||
 760                                 ((strcmp(anm[3*m],  hbdih->b[n].AI)==0) &&
 761                                  (strcmp(anm[1+m],  hbdih->b[n].AJ)==0) &&
 762                                  (strcmp(anm[2-m],  hbdih->b[n].AK)==0) &&
 763                                  (strcmp(anm[3-3*m],hbdih->b[n].AL)==0)));
 764                       if (bFound) {
 765                         set_p_string(&dih[ndih],hbdih->b[n].s);
 766
 767                         /* Set the last parameter to be able to see
 768                            if the dihedral was in the rtp list.
 769                            */
 770                         dih[ndih].c[MAXFORCEPARAM-1] = DIHEDRAL_WAS_SET_IN_RTP;
 771                         nFound++;
 772                         ndih++;
 773                         /* Set the next direct in case the rtp contains
 774                            multiple entries for this dihedral.
 775                            */
 776                         if (ndih == maxdih) {
 777                           maxdih += ninc;
 778                           srenew(dih,maxdih);
 779                         }
 780                         dih[ndih].AI=i;
 781                         dih[ndih].AJ=j1;
 782                         dih[ndih].AK=k1;
 783                         dih[ndih].AL=l1;
 784                         for (m=0; m<MAXFORCEPARAM; m++)
 785                           dih[ndih].c[m]=NOTSET;
 786                       }
 787                     }
 788                   } while (res < maxres);
 789                 }
 790                 if (nFound == 0) {
 791                   if (ndih == maxdih) {
 792                     maxdih += ninc;
 793                     srenew(dih,maxdih);
 794                   }
 795                   dih[ndih].AI=i;
 796                   dih[ndih].AJ=j1;
 797                   dih[ndih].AK=k1;
 798                   dih[ndih].AL=l1;
 799                   for (m=0; m<MAXFORCEPARAM; m++)
 800                     dih[ndih].c[m]=NOTSET;
 801                   set_p_string(&(dih[ndih]),"");
 802                   ndih++;
 803                 }
 804
 805                 nbd=nb_dist(nnb,i,l1);
 806                 if (debug)
 807                   fprintf(debug,"Distance (%d-%d) = %d\n",i+1,l1+1,nbd);
 808                 if (nbd == 3) {
 809                   i1 = min(i,l1);
 810                   i2 = max(i,l1);
 811                   bExcl = FALSE;
 812                   for(m=0; m<excls[i1].nr; m++)
 813                     bExcl = bExcl || excls[i1].e[m]==i2;
 814                   if (!bExcl) {
 815                     if (rtp[0].bGenerateHH14Interactions ||
 816                         !(is_hydro(atoms,i1) && is_hydro(atoms,i2))) {
 817                       if (npai == maxpai) {
 818                         maxpai += ninc;
 819                         srenew(pai,maxpai);
 820                       }
 821                       pai[npai].AI=i1;
 822                       pai[npai].AJ=i2;
 823                       pai[npai].C0=NOTSET;
 824                       pai[npai].C1=NOTSET;
 825                       set_p_string(&(pai[npai]),"");
 826                       npai++;
 827                     }
 828                   }
 829                 }
 830               }
 831             }
 832           }
 833         }
 834       }
 835     }
 836
 837   /* Sort angles with respect to j-i-k (middle atom first) */
 838   if (nang > 1)
 839     qsort(ang,nang,(size_t)sizeof(ang[0]),acomp);
 840
 841   /* Sort dihedrals with respect to j-k-i-l (middle atoms first) */
 842   if (ndih > 1)
 843     qsort(dih,ndih,(size_t)sizeof(dih[0]),dcomp);
 844
 845   /* Sort the pairs */
 846   if (npai > 1)
 847     qsort(pai,npai,(size_t)sizeof(pai[0]),pcomp);
 848   if (npai > 0) {
 849     /* Remove doubles, could occur in 6-rings, such as phenyls,
 850        maybe one does not want this when fudgeQQ < 1.
 851        */
 852     fprintf(stderr,"Before cleaning: %d pairs\n",npai);
 853     rm2par(pai,&npai,preq);
 854   }
 855
 856   /* Get the impropers from the database */
 857   nimproper = get_impropers(atoms,hb,&improper,bAllowMissing);
 858
 859   /* Sort the impropers */
 860   sort_id(nimproper,improper);
 861
 862   if (ndih > 0) {
 863     fprintf(stderr,"Before cleaning: %d dihedrals\n",ndih);
 864     clean_dih(dih,&ndih,improper,nimproper,atoms,
 865               rtp[0].bKeepAllGeneratedDihedrals,
 866               rtp[0].bRemoveDihedralIfWithImproper);
 867   }
 868
 869   /* Now we have unique lists of angles and dihedrals
 870    * Copy them into the destination struct
 871    */
 872   cppar(ang, nang, plist,F_ANGLES);
 873   cppar(dih, ndih, plist,F_PDIHS);
 874   cppar(improper,nimproper,plist,F_IDIHS);
 875   cppar(pai, npai, plist,F_LJ14);
 876
 877   /* Remove all exclusions which are within nrexcl */
 878   clean_excls(nnb,rtp[0].nrexcl,excls);
 879
 880   sfree(ang);
 881   sfree(dih);
 882   sfree(improper);
 883   sfree(pai);
 884 }
 885