src/kernel/gen_ad.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_gen_ad_c = "$Id$";
  30
  31 #include <math.h>
  32 #include <ctype.h>
  33 #include "sysstuff.h"
  34 #include "macros.h"
  35 #include "smalloc.h"
  36 #include "string2.h"
  37 #include "confio.h"
  38 #include "vec.h"
  39 #include "pbc.h"
  40 #include "toputil.h"
  41 #include "topio.h"
  42 #include "topexcl.h"
  43 #include "symtab.h"
  44 #include "macros.h"
  45 #include "fatal.h"
  46 #include "pgutil.h"
  47 #include "resall.h"
  48
  49 typedef bool (*peq)(t_param *p1, t_param *p2);
  50
  51 static int acomp(const void *a1, const void *a2)
  52 {
  53   t_param *p1,*p2;
  54   int     ac;
  55
  56   p1=(t_param *)a1;
  57   p2=(t_param *)a2;
  58   if ((ac=(p1->AJ-p2->AJ))!=0)
  59     return ac;
  60   else if ((ac=(p1->AI-p2->AI))!=0)
  61     return ac;
  62   else
  63     return (p1->AK-p2->AK);
  64 }
  65
  66 static int pcomp(const void *a1, const void *a2)
  67 {
  68   t_param *p1,*p2;
  69   int     pc;
  70
  71   p1=(t_param *)a1;
  72   p2=(t_param *)a2;
  73   if ((pc=(p1->AI-p2->AI))!=0)
  74     return pc;
  75   else
  76     return (p1->AJ-p2->AJ);
  77 }
  78
  79 static int dcomp(const void *d1, const void *d2)
  80 {
  81   t_param *p1,*p2;
  82   int     dc;
  83
  84   p1=(t_param *)d1;
  85   p2=(t_param *)d2;
  86   if ((dc=(p1->AJ-p2->AJ))!=0)
  87     return dc;
  88   else if ((dc=(p1->AK-p2->AK))!=0)
  89     return dc;
  90   else if ((dc=(p1->AI-p2->AI))!=0)
  91     return dc;
  92   else
  93     return (p1->AL-p2->AL);
  94 }
  95
  96 static bool aeq(t_param *p1, t_param *p2)
  97 {
  98   if (p1->AJ!=p2->AJ)
  99     return FALSE;
 100   else if (((p1->AI==p2->AI) && (p1->AK==p2->AK)) ||
 101            ((p1->AI==p2->AK) && (p1->AK==p2->AI)))
 102     return TRUE;
 103   else
 104     return FALSE;
 105 }
 106
 107 static bool deq2(t_param *p1, t_param *p2)
 108 {
 109   /* if bAlldih is true, dihedrals are only equal when
 110      ijkl = ijkl or ijkl =lkji*/
 111   if (((p1->AI==p2->AI) && (p1->AJ==p2->AJ) &&
 112        (p1->AK==p2->AK) && (p1->AL==p2->AL)) ||
 113       ((p1->AI==p2->AL) && (p1->AJ==p2->AK) &&
 114        (p1->AK==p2->AJ) && (p1->AL==p2->AI)))
 115     return TRUE;
 116   else
 117     return FALSE;
 118 }
 119
 120 static bool deq(t_param *p1, t_param *p2)
 121 {
 122   if (((p1->AJ==p2->AJ) && (p1->AK==p2->AK)) ||
 123       ((p1->AJ==p2->AK) && (p1->AK==p2->AJ)))
 124     return TRUE;
 125   else
 126     return FALSE;
 127 }
 128
 129 static bool remove_dih(t_param *p, int i, int np)
 130      /* check if dihedral p[i] should be removed */
 131 {
 132   bool bMidEq,bRem;
 133   int j;
 134
 135   if (i>0)
 136     bMidEq = deq(&p[i],&p[i-1]);
 137   else
 138     bMidEq = FALSE;
 139
 140   if (p[i].c[MAXFORCEPARAM-1]==NOTSET) {
 141     /* also remove p[i] if there is a dihedral on the same bond
 142        which has parameters set */
 143     bRem = bMidEq;
 144     j=i+1;
 145     while (!bRem && (j<np) && deq(&p[i],&p[j])) {
 146       bRem = (p[j].c[MAXFORCEPARAM-1] != NOTSET);
 147       j++;
 148     }
 149   } else
 150     bRem = bMidEq && (((p[i].AI==p[i-1].AI) && (p[i].AL==p[i-1].AL)) ||
 151                       ((p[i].AI==p[i-1].AL) && (p[i].AL==p[i-1].AI)));
 152
 153   return bRem;
 154 }
 155
 156 static bool preq(t_param *p1, t_param *p2)
 157 {
 158   if ((p1->AI==p2->AI) && (p1->AJ==p2->AJ))
 159     return TRUE;
 160   else
 161     return FALSE;
 162 }
 163
 164 static void rm2par(t_param p[], int *np, peq eq)
 165 {
 166   int *index,nind;
 167   int i,j;
 168
 169   if ((*np)==0)
 170     return;
 171
 172   snew(index,*np);
 173   nind=0;
 174     index[nind++]=0;
 175   for(i=1; (i<(*np)); i++)
 176     if (!eq(&p[i],&p[i-1]))
 177       index[nind++]=i;
 178   /* Index now holds pointers to all the non-equal params,
 179    * this only works when p is sorted of course
 180    */
 181   for(i=0; (i<nind); i++) {
 182     for(j=0; (j<MAXATOMLIST); j++)
 183       p[i].a[j]=p[index[i]].a[j];
 184     for(j=0; (j<MAXFORCEPARAM); j++)
 185       p[i].c[j]=p[index[i]].c[j];
 186     if (p[index[i]].a[0] == p[index[i]].a[1]) {
 187       if (debug)
 188         fprintf(debug,
 189                 "Something VERY strange is going on in rm2par (gen_ad.c)\n"
 190                 "a[0] %u a[1] %u a[2] %u a[3] %u\n",
 191                 p[i].a[0],p[i].a[1],p[i].a[2],p[i].a[3]);
 192       p[i].s = strdup("");
 193     } else if (index[i] > i) {
 194       /* Copy the string only if it comes from somewhere else
 195        * otherwise we will end up copying a random (newly freed) pointer.
 196        * Since the index is sorted we only have to test for index[i] > i.
 197        */
 198       sfree(p[i].s);
 199       p[i].s = strdup(p[index[i]].s);
 200     }
 201   }
 202   (*np)=nind;
 203
 204   sfree(index);
 205 }
 206
 207 static void cppar(t_param p[], int np, t_params plist[], int ftype)
 208 {
 209   int      i,j,nral,nrfp;
 210   t_params *ps;
 211
 212   ps   = &plist[ftype];
 213   nral = NRAL(ftype);
 214   nrfp = NRFP(ftype);
 215
 216   /* Keep old stuff */
 217   pr_alloc(np,ps);
 218   for(i=0; (i<np); i++) {
 219     for(j=0; (j<nral); j++)
 220       ps->param[ps->nr].a[j] = p[i].a[j];
 221     for(j=0; (j<nrfp); j++)
 222       ps->param[ps->nr].c[j] = p[i].c[j];
 223     ps->param[ps->nr].s=strdup(p[i].s);
 224     ps->nr++;
 225   }
 226 }
 227
 228 static void cpparam(t_param *dest, t_param *src)
 229 {
 230   int j;
 231
 232   for(j=0; (j<MAXATOMLIST); j++)
 233     dest->a[j]=src->a[j];
 234   for(j=0; (j<MAXFORCEPARAM); j++)
 235     dest->c[j]=src->c[j];
 236   dest->s=strdup(src->s);
 237 }
 238
 239 static void set_p(t_param *p, atom_id ai[4], real *c, char *s)
 240 {
 241   int j;
 242
 243   for(j=0; (j<4); j++)
 244     p->a[j]=ai[j];
 245   for(j=0; (j<MAXFORCEPARAM); j++)
 246     if (c)
 247       p->c[j]=c[j];
 248     else
 249       p->c[j]=NOTSET;
 250   sfree(p->s);
 251   if (s)
 252     p->s=strdup(s);
 253   else
 254     p->s=strdup("");
 255 }
 256
 257 static int int_comp(const void *a,const void *b)
 258 {
 259   return (*(int *)a) - (*(int *)b);
 260 }
 261
 262 static int eq_imp(atom_id a1[],atom_id a2[])
 263 {
 264   int b1[4],b2[4];
 265   int j;
 266
 267   for(j=0; (j<4); j++) {
 268     b1[j]=a1[j];
 269     b2[j]=a2[j];
 270   }
 271   qsort(b1,4,(size_t)sizeof(b1[0]),int_comp);
 272   qsort(b2,4,(size_t)sizeof(b2[0]),int_comp);
 273
 274   for(j=0; (j<4); j++)
 275     if (b1[j] != b2[j])
 276       return FALSE;
 277
 278   return TRUE;
 279 }
 280
 281 static bool ideq(t_param *p1, t_param *p2)
 282 {
 283   return eq_imp(p1->a,p2->a);
 284 }
 285
 286 static int idcomp(const void *a,const void *b)
 287 {
 288   t_param *pa,*pb;
 289   int     d;
 290
 291   pa=(t_param *)a;
 292   pb=(t_param *)b;
 293   if ((d=(pa->a[0]-pb->a[0])) != 0)
 294     return d;
 295   else if ((d=(pa->a[3]-pb->a[3])) != 0)
 296     return d;
 297   else if ((d=(pa->a[1]-pb->a[1])) != 0)
 298     return d;
 299   else
 300     return (int) (pa->a[2]-pb->a[2]);
 301 }
 302
 303 static void sort_id(int nr,t_param ps[])
 304 {
 305   qsort(ps,nr,(size_t)sizeof(ps[0]),idcomp);
 306 }
 307
 308 static t_rbonded *is_imp(t_param *p,t_atoms *atoms,t_hackblock hb[])
 309 {
 310   int        j,n,maxresnr,start;
 311   atom_id    aa0,a0[MAXATOMLIST];
 312   char      *atom;
 313   t_rbondeds *idihs;
 314
 315   if (!hb)
 316     return NULL;
 317   /* Find the max residue number in this dihedral */
 318   maxresnr=0;
 319   for(j=0; j<4; j++)
 320     maxresnr=max(maxresnr,atoms->atom[p->a[j]].resnr);
 321
 322   /* Now find the start of this residue */
 323   for(start=0; start < atoms->nr; start++)
 324     if (atoms->atom[start].resnr == maxresnr)
 325       break;
 326
 327   /* See if there are any impropers defined for this residue */
 328   idihs = &hb[atoms->atom[start].resnr].rb[ebtsIDIHS];
 329   for(n=0; n < idihs->nb; n++) {
 330     for(j=0; (j<4); j++) {
 331       atom=idihs->b[n].a[j];
 332       aa0=search_atom(atom,start,atoms->nr,atoms->atom,atoms->atomname);
 333       if (aa0 == NO_ATID) {
 334         if (debug)
 335           fprintf(debug,"Atom %s not found in res %d (maxresnr=%d) "
 336                   "in is_imp\n",atom,atoms->atom[start].resnr,maxresnr);
 337         break;
 338       } else
 339         a0[j] = aa0;
 340     }
 341     if (j==4) /* Not broken out */
 342       if (eq_imp(p->a,a0))
 343         return &idihs->b[n];
 344   }
 345   return NULL;
 346 }
 347
 348 static int n_hydro(atom_id a[],char ***atomname)
 349 {
 350   int i,nh=0;
 351   char c0,c1,*aname;
 352
 353   for(i=0; (i<4); i+=3) {
 354     aname=*atomname[a[i]];
 355     c0=toupper(aname[0]);
 356     if (c0 == 'H')
 357       nh++;
 358     else if (((int)strlen(aname) > 1) && (c0 >= '0') && (c0 <= '9')) {
 359       c1=toupper(aname[1]);
 360       if (c1 == 'H')
 361         nh++;
 362     }
 363   }
 364   return nh;
 365 }
 366
 367 static void pdih2idih(t_param *alldih, int *nalldih,t_param idih[],int *nidih,
 368                       t_atoms *atoms, t_hackblock hb[],bool bAlldih)
 369 {
 370   t_param   *dih;
 371   int       ndih;
 372   char      *a0;
 373   t_rbondeds *idihs;
 374   t_rbonded  *hbidih;
 375   int       i,j,k,l,start,aa0;
 376   int       *index,nind;
 377   atom_id   ai[MAXATOMLIST];
 378   bool      bStop,bIsSet,bKeep;
 379   int bestl,nh,minh;
 380
 381   /* First add all the impropers from the residue database
 382    * to the list.
 383    */
 384   start=0;
 385   if (hb) {
 386     for(i=0; (i<atoms->nres); i++) {
 387       idihs=&hb[i].rb[ebtsIDIHS];
 388       for(j=0; (j<idihs->nb); j++) {
 389         bStop=FALSE;
 390         for(k=0; (k<4) && !bStop; k++) {
 391           ai[k]=search_atom(idihs->b[j].a[k],start,
 392                             atoms->nr,atoms->atom,atoms->atomname);
 393           if (ai[k] == NO_ATID) {
 394             if (debug)
 395               fprintf(debug,"Atom %s (%d) not found in res %d in pdih2idih\n",
 396                       idihs->b[j].a[k], k, atoms->atom[start].resnr);
 397             bStop=TRUE;
 398           }
 399         }
 400         if (!bStop) {
 401           /* Not broken out */
 402           set_p(&idih[*nidih],ai,NULL,idihs->b[j].s);
 403         (*nidih)++;
 404         }
 405       }
 406       while ((start<atoms->nr) && (atoms->atom[start].resnr==i))
 407         start++;
 408     }
 409   }
 410   if (*nalldih == 0)
 411     return;
 412
 413   /* Copy the impropers and dihedrals to seperate arrays. */
 414   snew(dih,*nalldih);
 415   ndih = 0;
 416   for(i=0; i<*nalldih; i++) {
 417     hbidih = is_imp(&alldih[i],atoms,hb);
 418     if ( hbidih ) {
 419       set_p(&idih[*nidih],alldih[i].a,NULL,hbidih->s);
 420       (*nidih)++;
 421     } else {
 422       cpparam(&dih[ndih],&alldih[i]);
 423       ndih++;
 424     }
 425   }
 426
 427   /* Now, because this list still contains the double entries,
 428    * keep the dihedral with parameters or the first one.
 429    */
 430
 431   snew(index,ndih);
 432   nind=0;
 433   if (bAlldih) {
 434     fprintf(stderr,"bAlldih = true\n");
 435     for(i=0; i<ndih; i++)
 436       if ((i==0) || !deq2(&dih[i],&dih[i-1]))
 437         index[nind++]=i;
 438   } else {
 439     for(i=0; i<ndih; i++)
 440       if (!remove_dih(dih,i,ndih))
 441         index[nind++]=i;
 442   }
 443   index[nind]=ndih;
 444
 445   /* if we don't want all dihedrals, we need to select the ones with the
 446    *  fewest hydrogens
 447    */
 448
 449   k=0;
 450   for(i=0; i<nind; i++) {
 451     bIsSet = (dih[index[i]].c[MAXFORCEPARAM-1] != NOTSET);
 452     bKeep = TRUE;
 453     if (!bIsSet)
 454       /* remove the dihedral if there is an improper on the same bond */
 455       for(j=0; (j<(*nidih)) && bKeep; j++)
 456         bKeep = !deq(&dih[index[i]],&idih[j]);
 457
 458     if (bKeep) {
 459       /* Now select the "fittest" dihedral:
 460        * the one with as few hydrogens as possible
 461        */
 462
 463       /* Best choice to get dihedral from */
 464       bestl=index[i];
 465       if (!bAlldih && !bIsSet) {
 466         /* Minimum number of hydrogens for i and l atoms */
 467         minh=2;
 468         for(l=index[i]; (l<index[i+1]) && deq(&dih[index[i]],&dih[l]); l++) {
 469           if ((nh=n_hydro(dih[l].a,atoms->atomname)) < minh) {
 470             minh=nh;
 471             bestl=l;
 472           }
 473           if (minh == 0)
 474             break;
 475         }
 476       }
 477       for(j=0; (j<MAXATOMLIST); j++)
 478         alldih[k].a[j] = dih[bestl].a[j];
 479       for(j=0; (j<MAXFORCEPARAM); j++)
 480         alldih[k].c[j] = dih[bestl].c[j];
 481       sfree(alldih[k].s);
 482       alldih[k].s = strdup(dih[bestl].s);
 483       k++;
 484     }
 485   }
 486   for (i=k; i < *nalldih; i++)
 487     sfree(alldih[i].s);
 488   *nalldih = k;
 489
 490   for(i=0; i<ndih; i++)
 491     sfree(dih[i].s);
 492   sfree(dih);
 493   sfree(index);
 494 }
 495
 496 static int nb_dist(t_nextnb *nnb,int ai,int aj)
 497 {
 498   int nre,nrx,NRE;
 499   int *nrexcl;
 500   int *a;
 501
 502   if (ai == aj)
 503     return 0;
 504
 505   NRE=-1;
 506   nrexcl=nnb->nrexcl[ai];
 507   for(nre=1; (nre < nnb->nrex); nre++) {
 508     a=nnb->a[ai][nre];
 509     for(nrx=0; (nrx < nrexcl[nre]); nrx++) {
 510       if ((aj == a[nrx]) && (NRE == -1))
 511         NRE=nre;
 512     }
 513   }
 514   return NRE;
 515 }
 516
 517 bool is_hydro(t_atoms *atoms,int ai)
 518 {
 519   return ((*(atoms->atomname[ai]))[0] == 'H');
 520 }
 521
 522 static void get_atomnames_min(int n,char anm[4][12],
 523                               int res,t_atoms *atoms,atom_id *a)
 524 {
 525   int m;
 526
 527   /* Assume ascending residue numbering */
 528   for(m=0; m<n; m++) {
 529     if (atoms->atom[a[m]].resnr < res)
 530       strcpy(anm[m],"-");
 531     else if (atoms->atom[a[m]].resnr > res)
 532       strcpy(anm[m],"+");
 533     else
 534       strcpy(anm[m],"");
 535     strcat(anm[m],*(atoms->atomname[a[m]]));
 536   }
 537 }
 538
 539 void gen_pad(t_nextnb *nnb, t_atoms *atoms, bool bH14,
 540              t_params plist[], t_hackblock hb[], bool bAlldih)
 541 {
 542   t_param *ang,*dih,*pai,*idih;
 543   t_rbondeds *hbang, *hbdih;
 544   char    anm[4][12];
 545   int     res,minres,maxres;
 546   int     i,j,j1,k,k1,l,l1,m,n;
 547   int     maxang,maxdih,maxidih,maxpai;
 548   int     nang,ndih,npai,nidih,nbd;
 549   int     dang,ddih;
 550   bool    bFound;
 551
 552   /* These are the angles, pairs, impropers and dihedrals that we generate
 553    * from the bonds. The ones that are already there from the rtp file
 554    * will be retained.
 555    */
 556   nang    = 0;
 557   nidih   = 0;
 558   npai    = 0;
 559   ndih    = 0;
 560   dang    = 6*nnb->nr;
 561   ddih    = 24*nnb->nr;
 562   maxang  = dang;
 563   maxdih  = maxpai = maxidih = ddih;
 564   snew(ang, maxang);
 565   snew(dih, maxdih);
 566   snew(pai, maxpai);
 567   snew(idih,maxidih);
 568
 569   /* extract all i-j-k-l neighbours from nnb struct */
 570   for(i=0; (i<nnb->nr); i++)
 571     /* For all particles */
 572     for(j=0; (j<nnb->nrexcl[i][1]); j++) {
 573       /* For all first neighbours */
 574       j1=nnb->a[i][1][j];
 575       for(k=0; (k<nnb->nrexcl[j1][1]); k++) {
 576         /* For all first neighbours of j1 */
 577         k1=nnb->a[j1][1][k];
 578         if (k1 != i) {
 579           if (nang == maxang) {
 580             srenew(ang,maxang+dang);
 581             memset(ang+maxang,0,sizeof(ang[0]));
 582             maxang += dang;
 583           }
 584           ang[nang].AJ=j1;
 585           if (i < k1) {
 586             ang[nang].AI=i;
 587             ang[nang].AK=k1;
 588           }
 589           else {
 590             ang[nang].AI=k1;
 591             ang[nang].AK=i;
 592           }
 593           ang[nang].C0=NOTSET;
 594           ang[nang].C1=NOTSET;
 595           ang[nang].s=strdup("");
 596           if (hb) {
 597             minres = atoms->atom[ang[nang].a[0]].resnr;
 598             maxres = minres;
 599             for(m=1; m<3; m++) {
 600               minres = min(minres,atoms->atom[ang[nang].a[m]].resnr);
 601               maxres = max(maxres,atoms->atom[ang[nang].a[m]].resnr);
 602             }
 603             res = 2*minres-maxres;
 604             do {
 605               res += maxres-minres;
 606               hbang=&hb[res].rb[ebtsANGLES];
 607               for(l=0; (l<hbang->nb); l++) {
 608                 get_atomnames_min(3,anm,res,atoms,ang[nang].a);
 609                 if (strcmp(anm[1],hbang->b[l].AJ)==0) {
 610                   bFound=FALSE;
 611                   for (m=0; m<3; m+=2)
 612                     bFound=(bFound ||
 613                             ((strcmp(anm[m],hbang->b[l].AI)==0) &&
 614                              (strcmp(anm[2-m],hbang->b[l].AK)==0)));
 615                   if (bFound) {
 616                     sfree(ang[nang].s);
 617                     if (hbang->b[l].s)
 618                       ang[nang].s = strdup(hbang->b[l].s);
 619                     else
 620                       ang[nang].s = strdup("");
 621                   }
 622                 }
 623               }
 624             } while (res < maxres);
 625           }
 626           nang++;
 627           for(l=0; (l<nnb->nrexcl[k1][1]); l++) {
 628             /* For all first neighbours of k1 */
 629             l1=nnb->a[k1][1][l];
 630             if ((l1 != i) && (l1 != j1)) {
 631               if (ndih == maxdih) {
 632                 srenew(dih,maxdih+ddih);
 633                 memset(dih+maxdih,0,sizeof(dih[0]));
 634                 srenew(idih,maxdih+ddih);
 635                 memset(idih+maxdih,0,sizeof(idih[0]));
 636                 srenew(pai,maxdih+ddih);
 637                 memset(pai+maxdih,0,sizeof(pai[0]));
 638                 maxdih += ddih;
 639               }
 640               if (j1 < k1) {
 641                 dih[ndih].AI=i;
 642                 dih[ndih].AJ=j1;
 643                 dih[ndih].AK=k1;
 644                 dih[ndih].AL=l1;
 645               }
 646               else {
 647                 dih[ndih].AI=l1;
 648                 dih[ndih].AJ=k1;
 649                 dih[ndih].AK=j1;
 650                 dih[ndih].AL=i;
 651               }
 652               for (m=0; m<MAXFORCEPARAM; m++)
 653                 dih[ndih].c[m]=NOTSET;
 654               dih[ndih].s=strdup("");
 655               if (hb) {
 656                 minres = atoms->atom[dih[ndih].a[0]].resnr;
 657                 maxres = minres;
 658                 for(m=1; m<4; m++) {
 659                   minres = min(minres,atoms->atom[dih[ndih].a[m]].resnr);
 660                   maxres = max(maxres,atoms->atom[dih[ndih].a[m]].resnr);
 661                 }
 662                 res = 2*minres-maxres;
 663                 do {
 664                   res += maxres-minres;
 665                   hbdih=&hb[res].rb[ebtsPDIHS];
 666                   for(n=0; (n<hbdih->nb); n++) {
 667                     get_atomnames_min(4,anm,res,atoms,dih[ndih].a);
 668                     bFound=FALSE;
 669                     for (m=0; m<2; m++)
 670                       bFound=(bFound ||
 671                               ((strcmp(anm[3*m],  hbdih->b[n].AI)==0) &&
 672                                (strcmp(anm[1+m],  hbdih->b[n].AJ)==0) &&
 673                                (strcmp(anm[2-m],  hbdih->b[n].AK)==0) &&
 674                                (strcmp(anm[3-3*m],hbdih->b[n].AL)==0)));
 675                     if (bFound) {
 676                       sfree(dih[ndih].s);
 677                       if (hbdih->b[n].s)
 678                       dih[ndih].s = strdup(hbdih->b[n].s);
 679                       else
 680                         dih[ndih].s = strdup("");
 681                       /* Set the last parameter to be able to see
 682                          if the dihedral was in the rtp list */
 683                       dih[ndih].c[MAXFORCEPARAM-1] = 0;
 684                     }
 685                   }
 686                 } while (res < maxres);
 687               }
 688               nbd=nb_dist(nnb,i,l1);
 689               if (debug)
 690                 fprintf(debug,"Distance (%d-%d) = %d\n",i+1,l1+1,nbd);
 691               if (nbd == 3) {
 692                 pai[npai].AI=min(i,l1);
 693                 pai[npai].AJ=max(i,l1);
 694                 pai[npai].C0=NOTSET;
 695                 pai[npai].C1=NOTSET;
 696                 sfree(pai[npai].s);
 697                 pai[npai].s =strdup("");
 698                 if (bH14 || !(is_hydro(atoms,pai[npai].AI) &&
 699                               is_hydro(atoms,pai[npai].AJ)))
 700                   npai++;
 701               }
 702
 703               ndih++;
 704             }
 705           }
 706         }
 707       }
 708     }
 709
 710   /* We now have a params list with double entries for each angle,
 711    * and even more for each dihedral. We will remove these now.
 712    */
 713   /* Sort angles with respect to j-i-k (middle atom first) */
 714   qsort(ang,nang,(size_t)sizeof(ang[0]),acomp);
 715   rm2par(ang,&nang,aeq);
 716
 717   /* Sort dihedrals with respect to j-k-i-l (middle atoms first) */
 718   fprintf(stderr,"before sorting: %d dihedrals\n",ndih);
 719   qsort(dih,ndih,(size_t)sizeof(dih[0]),dcomp);
 720   pdih2idih(dih,&ndih,idih,&nidih,atoms,hb,bAlldih);
 721   fprintf(stderr,"after sorting: %d dihedrals\n",ndih);
 722
 723   /* Now the dihedrals are sorted and doubles removed, this has to be done
 724    * for impropers too
 725    */
 726   sort_id(nidih,idih);
 727   rm2par(idih,&nidih,ideq);
 728
 729   /* And for the pairs */
 730   fprintf(stderr,"There are %d pairs before sorting\n",npai);
 731   qsort(pai,npai,(size_t)sizeof(pai[0]),pcomp);
 732   rm2par(pai,&npai,preq);
 733
 734   /* Now we have unique lists of angles and dihedrals
 735    * Copy them into the destination struct
 736    */
 737   cppar(ang, nang, plist,F_ANGLES);
 738   cppar(dih, ndih, plist,F_PDIHS);
 739   cppar(idih,nidih,plist,F_IDIHS);
 740   cppar(pai, npai, plist,F_LJ14);
 741
 742
 743   for(i=0; i<nang; i++)
 744     sfree(ang[i].s);
 745   sfree(ang);
 746
 747   for(i=0; i<ndih; i++)
 748     sfree(dih[i].s);
 749   sfree(dih);
 750
 751   for(i=0; i<nidih; i++)
 752     sfree(idih[i].s);
 753   sfree(idih);
 754
 755   for(i=0; i<npai; i++)
 756     sfree(pai[i].s);
 757   sfree(pai);
 758 }
 759