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