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