src/tools/smooth.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  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_smooth_c = "$Id$";
  30
  31 #define NMRLEN  99.9
  32 #define OVERLAP(a_lb,a_ub,b_lb,b_ub) !((a_ub)<(b_lb)) || ((b_ub)<(a_lb))
  33 /* Bound smoothing for cdist*/
  34 /* Adam Kirrander 990119 */
  35
  36 real _mysqrt(real x,char *fn,int line)
  37 {
  38   if (x < 0) {
  39     if (debug)
  40       fprintf(debug,"MYSQRT: negative argument %g called from %s, line %d\n",
  41               x,fn,line);
  42     return 0;
  43   }
  44   else
  45     return sqrt(x);
  46 }
  47
  48 #define mysqrt(x) _mysqrt(x,__FILE__,__LINE__)
  49
  50 void check_bounds(int ai,int aj,real lb,real ub)
  51 {
  52   if (lb > ub)
  53     fatal_error(0,"Strange bounds for %d-%d: lb=%f, ub=%f",
  54                 ai+1,aj+1,lb,ub);
  55 }
  56
  57 void check_len_(int ai,int aj,real lb,real ub,real len,char *file,int line)
  58 {
  59   if (len != NMRLEN)
  60     if (((len > ub) || (len < lb)) && len > 0.0)
  61       fatal_error(0,"%s, line %d: Ideal len. outside bounds for %d-%d:"
  62                   " lb=%f, ub=%f, len=%f",
  63                   file,line,ai+1,aj+1,lb,ub,len);
  64 }
  65
  66 #define check_len(ai,aj,lb,ub,len) check_len_(ai,aj,lb,ub,len,__FILE__,__LINE__)
  67
  68 void check_tri(t_dist *d,int natoms,int i,int j,int k)
  69 {
  70   real a_lb,a_ub,a_len,b_lb,b_ub,b_len,c_lb,c_ub,c_len;
  71
  72   a_lb  = d_lb(d,natoms,i,j);
  73   a_ub  = d_ub(d,natoms,i,j);
  74   a_len = d_len(d,natoms,i,j);
  75
  76   b_lb  = d_lb(d,natoms,i,k);
  77   b_ub  = d_ub(d,natoms,i,k);
  78   b_len = d_len(d,natoms,i,k);
  79
  80   c_lb  = d_lb(d,natoms,j,k);
  81   c_ub  = d_ub(d,natoms,j,k);
  82   c_len = d_len(d,natoms,j,k);
  83
  84   check_len(i,j,a_lb,a_ub,a_len);
  85   check_len(i,k,b_lb,b_ub,b_len);
  86   check_len(j,k,c_lb,c_ub,c_len);
  87
  88   check_bounds(i,j,a_lb,a_ub);
  89   check_bounds(i,k,b_lb,b_ub);
  90   check_bounds(j,k,c_lb,c_ub);
  91
  92   /* Check if triangle inequality is fulfilled */
  93   if ((a_lb > (b_ub+c_ub)) ||
  94       (b_lb > (a_ub+c_ub)) ||
  95       (c_lb > (a_ub+b_ub))) {
  96     fprintf(stderr,"Impossible triangle (%d,%d,%d):\n"\
  97             "a_lb=%10.5f, a_ub=%10.5f, a_len=%10.5f\n"\
  98             "b_lb=%10.5f, b_ub=%10.5f, b_len=%10.5f\n"
  99             "c_lb=%10.5f, c_ub=%10.5f, c_len=%10.5f\n",\
 100             i,j,k,a_lb,a_ub,a_len,b_lb,b_ub,b_len,c_lb,c_ub,c_len);
 101     exit(1);
 102   }
 103 }
 104
 105 int triangle_upper_bound (t_dist *d,int natoms,real tol)
 106 {
 107   double possible;
 108   int    i,j,k,nsmooth,innerloop,count=0,nloops=0,ntotal=0;
 109   real   a_lb,a_ub,a_len,b_lb,b_ub,b_len,c_lb,c_ub,c_len;
 110
 111   if (debug)
 112     fprintf(debug,"Just entered triangle_upper_bound!\n");
 113
 114   /* Loop over all triangles until no smoothing occurs */
 115   do {
 116     nsmooth=0;
 117     nloops=nloops+1;
 118     for (i=0; (i < natoms); i++) {
 119       for (j=i+1; (j < natoms); j++) {
 120         if (!dist_set(d,natoms,i,j)) continue;
 121         for (k=j+1; (k < natoms); k++) {
 122
 123           count=count+1;
 124
 125           /* Should make use of distances between zero-weighted if that can
 126              help us to define better distances for e.g. loops. But make sure
 127              we don't waste time smoothing zero-weighted tripples...*/
 128
 129           /* Are all distances defined? */
 130           if (!dist_set(d,natoms,i,k) || !dist_set(d,natoms,j,k)) continue;
 131
 132           /* smooth the triangle */
 133           do {
 134             /* Reset inner counter */
 135             innerloop=0;
 136
 137             check_tri(d,natoms,i,j,k);
 138
 139             /* Read bounds */
 140             a_lb  = d_lb(d,natoms,i,j);
 141             a_ub  = d_ub(d,natoms,i,j);
 142             a_len = d_len(d,natoms,i,j);
 143
 144             b_lb  = d_lb(d,natoms,i,k);
 145             b_ub  = d_ub(d,natoms,i,k);
 146             b_len = d_len(d,natoms,i,k);
 147
 148             c_lb  = d_lb(d,natoms,j,k);
 149             c_ub  = d_ub(d,natoms,j,k);
 150             c_len = d_len(d,natoms,j,k);
 151
 152             /* Check if triangle inequality is fulfilled */
 153             if (a_ub > (b_ub+c_ub+tol)) {
 154               set_dist(d,natoms,i,j,a_lb,b_ub+c_ub,a_len);
 155               a_ub = d_ub(d,natoms,i,j);
 156               innerloop++;
 157             }
 158             if (b_ub > (a_ub+c_ub+tol)) {
 159               set_dist(d,natoms,i,k,b_lb,a_ub+c_ub,b_len);
 160               b_ub = d_ub(d,natoms,i,k);
 161               innerloop++;
 162             }
 163             if (c_ub > (a_ub+b_ub+tol)) {
 164               set_dist(d,natoms,j,k,c_lb,a_ub+b_ub,c_len);
 165               c_ub = d_ub(d,natoms,j,k);
 166               innerloop++;
 167             }
 168
 169             nsmooth += innerloop;
 170           }
 171           while (innerloop > 0);
 172         }
 173       }
 174     }
 175     ntotal += nsmooth;
 176   }
 177   while (nsmooth>0);
 178
 179   possible = natoms*(natoms-1.0)*(natoms-2.0)/6;
 180   fprintf(stderr,"Checked %d triples (of %.0f) in %d rounds of ub"
 181           " triangle smoothing.\n",
 182           count/nloops,possible,nloops);
 183   fprintf(stderr,"Smoothed %d upper bounds with triagonal ineq.\n",ntotal);
 184
 185   return ntotal;
 186 }
 187
 188
 189 int triangle_lower_bound (t_dist *d,int natoms,real tol)
 190 {
 191
 192   int i,j,k,nsmooth,innerloop,ntotal=0;
 193   real a_lb,a_ub,a_len,b_lb,b_ub,b_len,c_lb,c_ub,c_len,new_lb;
 194
 195   /*fprintf(stderr,"Just entered triangle_lower_bound!\n");*/
 196
 197   /* Loop over all triangles until no smoothing occurs */
 198   do {
 199     nsmooth=0;
 200     for (i=0; (i < natoms); i++) {
 201       for (j=i+1; (j < natoms); j++) {
 202         if (!dist_set(d,natoms,i,j)) continue;
 203         for (k=j+1; (k < natoms); k++) {
 204
 205           /* Are all distances defined? If not, continue */
 206           if (!dist_set(d,natoms,i,k) || !dist_set(d,natoms,j,k)) continue;
 207
 208
 209           /* smooth the triangle */
 210           do {
 211             /* Reset inner counter */
 212             innerloop=0;
 213
 214             /* Read bounds */
 215             a_lb  = d_lb(d,natoms,i,j);
 216             a_ub  = d_ub(d,natoms,i,j);
 217             a_len = d_len(d,natoms,i,j);
 218
 219             b_lb = d_lb(d,natoms,i,k);
 220             b_ub = d_ub(d,natoms,i,k);
 221             b_len = d_len(d,natoms,i,k);
 222
 223             c_lb  = d_lb(d,natoms,j,k);
 224             c_ub  = d_ub(d,natoms,j,k);
 225             c_len = d_len(d,natoms,j,k);
 226
 227             /* Smooth lower bound */
 228             if (!OVERLAP(b_lb,b_ub,c_lb,c_ub)) {
 229               new_lb = min(abs(b_lb-c_ub),abs(c_lb-b_ub));
 230               if ( tol < new_lb - a_lb) {
 231                 set_dist(d,natoms,i,j,new_lb,a_ub,a_len);
 232                 nsmooth++;
 233                 innerloop++;
 234                 /*fprintf(stderr,"Smoothed %d-%d",i,j);*/
 235               }
 236             }
 237             if (!OVERLAP(a_lb,a_ub,c_lb,c_ub)) {
 238               new_lb = min(abs(a_lb-c_ub),abs(c_lb-a_ub));
 239               if ( tol < new_lb - b_lb) {
 240                 set_dist(d,natoms,i,k,new_lb,b_ub,b_len);
 241                 nsmooth++;
 242                 innerloop++;
 243                 /*fprintf(stderr,"Smoothed %d-%d",i,k);*/
 244               }
 245             }
 246             if (!OVERLAP(a_lb,a_ub,b_lb,b_ub)) {
 247               new_lb = min(abs(a_lb-b_ub),abs(b_lb-a_ub));
 248               /*Programming error? */
 249               /* if ( (a_ub+b_ub) < c_ub ) { */
 250               /* New if-statement 990609 Adam */
 251               if ( tol < new_lb - c_lb) {
 252                 set_dist(d,natoms,j,k,new_lb,c_ub,c_len);
 253                 nsmooth++;
 254                 innerloop++;
 255                 /*fprintf(stderr,"Smoothed %d-%d",j,k);*/
 256               }
 257             }
 258
 259             /* Check that bounds make sense */
 260             check_bounds(i,j,a_lb,a_ub);
 261             check_bounds(i,k,b_lb,b_ub);
 262             check_bounds(j,k,c_lb,c_ub);
 263           }
 264           while (innerloop > 0);
 265         }
 266       }
 267     }
 268     ntotal += nsmooth;
 269   }
 270   while (nsmooth>0);
 271   fprintf(stderr,"Smoothed %d lower bounds with triagonal ineq.\n",ntotal);
 272   return ntotal;
 273 }
 274
 275
 276 int do_triangle (t_dist *d,t_atoms *atoms,real tol)
 277 {
 278   int ntot=0;
 279   /* Do triangular smoothing only */
 280   int natoms = atoms->nr;
 281
 282   fprintf(stderr,"Starting triangle smoothing\n");
 283   ntot += triangle_upper_bound (d,natoms,tol);
 284   ntot += triangle_lower_bound (d,natoms,tol);
 285
 286   return ntot;
 287 }
 288
 289 /* Tetrangle-bound smoothing for cdist.
 290    Adam Kirrander 990209 */
 291
 292 /* According Havel,Kuntz and Crippen 1983 Bull.Math.Biol */
 293
 294
 295
 296
 297 /* FUNCTIONS ASSOCIATED WITH TRIANGLE_LIMITED */
 298 /* Routines to check if triangle inequalty limits are attainable between
 299    points 1 and 4, given that the 5 other distances are defined.
 300    Sets logical TUL=true if the upper limit of distance 1-4 is defined by
 301    triangle inequalty, and TLL=true if corresponding lower limit is set by
 302    the triangle ineq. Written to be used in conjunction with tetragonal
 303    bound smoothing in cdist.
 304
 305    The main routine is trianglelimited, and it is the one that should be
 306    called. The rest are supporting routines.
 307    Ref.: Bull.Math.Biol.Vol45,Nr5,665-720,Havel,Crippen,Kuntz,1983 */
 308
 309 /* Adam Kirrander    990209  */
 310
 311 /*#define OVERLAP(a_lb,a_ub,b_lb,b_ub) !((a_ub<b_lb) || (b_ub<a_lb))*/
 312 /* ((a_ub >= b_lb) && (b_ub >= a_lb)) */
 313 #define apex(AC,AD,BC,BD) mysqrt((BC*(AD*AD-BD*BD)+BD*(AC*AC-BC*BC))/(BC+BD))
 314
 315
 316 bool _triangle_bound_attainable(real BC,real BD,real AClb,real ACub,
 317                                 real ADlb,real ADub,real ABlb,real ABub,
 318                                 char *fn,int line)
 319 {
 320   /* Assume triangle with A,C and D at corners and B at C-D.
 321      Given the distances AC,AD,BC,BD it calculates the distance AB
 322      at its min and max, which corresponds to AC and AD being min and
 323      max respectively.
 324      If this intervall overlaps with the set bounds for AB, it follows that
 325      the points C,B and D can be made colinear and thus are restricted by the
 326      triangle inequality. */
 327
 328   real ABmin,ABmax;
 329
 330   if (debug)
 331     fprintf(debug,"%s, line %d: ABlb: %g, ABub: %g, AClb: %g, ACub: %g, ADlb: %g, ADub: %g, BC: %g, BD: %g\n",
 332             fn,line,ABlb,ABub,AClb,ACub,ADlb,ADub,BC,BD);
 333   ABmax = apex(ACub,ADub,BC,BD);
 334     /* Triangle can't be formed, i.e. there is a risk all 4 points can
 335      become colinear */
 336   if ( (AClb+ADlb) < (BC+BD) )
 337     return ( ABlb <= ABmax );
 338
 339   /* I am not sure the above truly is the optimal solution, but it should
 340    * be the safe solution
 341    */
 342
 343   ABmin = apex(AClb,ADlb,BC,BD);
 344
 345   return (OVERLAP(ABlb,ABub,ABmin,ABmax));
 346 }
 347
 348 #define triangle_bound_attainable(BC,BD,AClb,ACub,ADlb,ADub,ABlb,ABub) \
 349        _triangle_bound_attainable(BC,BD,AClb,ACub,ADlb,ADub,ABlb,ABub, \
 350                                   __FILE__,__LINE__)
 351
 352
 353 void triangle_limited(int *ATMS,int natoms,t_dist *d,bool *TUL, bool *TLL)
 354 {
 355   /* Given 4 atoms it checks if the triangle inequality lower or upper bounds
 356      for the distance 1-4 are attainable. */
 357   /* Situation looks something like the following: */
 358   /* Can:                                          */
 359   /*   1                      4                    */
 360   /*    \\\               ////                     */
 361   /*        2------------3                         */
 362   /* become:                                       */
 363   /*   1---------------3------4                    */
 364   /*    \\\\\\\ ///////                            */
 365   /*           2                                   */
 366
 367
 368   int N1=ATMS[0],N2=ATMS[1],N3=ATMS[2],N4=ATMS[3];
 369   real d12[2],d13[2],d14[2],d23[2],d24[2],d34[2];
 370
 371
 372   /* Initiate and read bounds */
 373   *TUL = FALSE;
 374   *TLL = FALSE;
 375
 376   d12[0] = d_lb(d,natoms,N1,N2);
 377   d12[1] = d_ub(d,natoms,N1,N2);
 378
 379   d13[0] = d_lb(d,natoms,N1,N3);
 380   d13[1] = d_ub(d,natoms,N1,N3);
 381
 382   d14[0] = d_lb(d,natoms,N1,N4);
 383   d14[1] = d_ub(d,natoms,N1,N4);
 384
 385   d23[0] = d_lb(d,natoms,N2,N3);
 386   d23[1] = d_ub(d,natoms,N2,N3);
 387
 388   d24[0] = d_lb(d,natoms,N2,N4);
 389   d24[1] = d_ub(d,natoms,N2,N4);
 390
 391   d34[0] = d_lb(d,natoms,N3,N4);
 392   d34[1] = d_ub(d,natoms,N3,N4);
 393
 394   /* Check if UPPER triangle inequality limit is attainable. */
 395   if ( d12[1] + d24[1] < d13[1] + d34[1] ) {
 396     /* Corresponds to N1,N2 and N4 colinear.
 397        N1=D,N2=B,N3=A,N4=C ;in notation of subroutines */
 398     *TUL = triangle_bound_attainable(d24[1],d12[1],d34[0],d34[1],
 399                                      d13[0],d13[1],d23[0],d23[1]);
 400   }
 401   else if ( d12[1] + d24[1] > d13[1] + d34[1] ) {
 402     /* Corresponds to N1,N3 and N4 colinear.
 403        N1=D,N2=A,N3=B,N4=C ;in notation of subroutines */
 404     *TUL = triangle_bound_attainable(d34[1],d13[1],d24[0],d24[1],
 405                                      d12[0],d12[1],d23[0],d23[1]);
 406   }
 407   /* if N2 and N3 can superimpose TUL is true by necessity (AB=0) */
 408   else if (d23[0] == 0) *TUL = TRUE;
 409
 410   /* Check if LOWER triangle inequality limit is attainable */
 411   if ( (d13[0] - d34[1] <= 0) && (d34[0] - d13[1] <= 0) &&
 412        (d24[0] - d12[1] <= 0) && (d12[0] - d24[1] <= 0) ) {
 413     /* if all inverse triangle ineq. limits are zero then */
 414     *TLL = TRUE;
 415   }
 416   else if (d12[0] - d24[1] > 0) {
 417     /* Corresponds to N1,N4 and N2 colinear.
 418        A=N3,B=N4,C=N1,D=N2 */
 419     *TLL = triangle_bound_attainable(d24[1],d12[0]-d24[1],d23[0],d23[1],
 420                                      d13[0],d13[1],d34[0],d34[1]);
 421   }
 422   else if (d24[0] - d12[1] > 0) {
 423     /* Corresponds to N2,N1 and N4 colinear.
 424        A=N3,B=N1,C=N2,D=N4 */
 425     *TLL = triangle_bound_attainable(d12[1],d24[0]-d12[1],d23[0],d23[1],
 426                                      d34[0],d34[1],d13[0],d13[1]);
 427   }
 428   else if (d13[0] - d34[1] > 0) {
 429     /* Corresponds to N1,N4 and N3 colinear.
 430        A=N2,B=N4,C=N3,D=N1 */
 431     *TLL = triangle_bound_attainable(d34[1],d13[0]-d34[1],d23[0],d23[1],
 432                                      d12[0],d12[1],d24[0],d24[1]);
 433   }
 434   else if (d34[0] - d13[1] > 0) {
 435     /* Corresponds to N3,N1 and N4 colinear.
 436        A=N2,B=N1,C=N3,D=N4 */
 437     *TLL = triangle_bound_attainable(d13[1],d34[0]-d13[1],d23[0],d23[1],
 438                                      d24[0],d24[1],d12[0],d12[1]);
 439   }
 440
 441 }
 442 /* END OF FUNCTIONS >TRIANGLE_LIMITED< */
 443
 444
 445
 446 /* FUNCTIONS ASSOCIATED WITH TETRAGONAL-SMOOTHING */
 447
 448 void minmax(real array[],int SIZE,int *min,int *max)
 449 {
 450   /* Finds min and max (indices thereof) in array.
 451      Early version that only accepts real.
 452      Fix it!
 453      Adam Kirrander 990211                         */
 454
 455   int i;
 456   *min = 0;
 457   *max = 0;
 458
 459   for (i=1; (i<SIZE); i++) {
 460     if (array[i] < array[*min])
 461       *min = i;
 462     else if (array[i] > array[*max])
 463       *max = i;
 464   }
 465 }
 466
 467 void swap (int *x,int *y)
 468 {
 469   /* Swaps x and y.
 470      Early version that only can swap integers.
 471      Fix it!
 472      Adam Kirrander 990211                         */
 473
 474   int temp;
 475   temp = *x;
 476   *x = *y;
 477   *y = temp;
 478 }
 479
 480
 481 bool alldistset (int i,int j,int k,int l,int natoms,t_dist *d)
 482 {
 483   /* Returns FALSE if all distances are not set */
 484   if (!dist_set(d,natoms,i,j)) return FALSE;
 485   if (!dist_set(d,natoms,i,k)) return FALSE;
 486   if (!dist_set(d,natoms,i,l)) return FALSE;
 487   if (!dist_set(d,natoms,j,k)) return FALSE;
 488   if (!dist_set(d,natoms,j,l)) return FALSE;
 489   if (!dist_set(d,natoms,k,l)) return FALSE;
 490
 491   return TRUE;
 492 }
 493
 494 int nr_nb (int i,int j,int k,int l,int natoms,t_dist *d)
 495 {
 496   /* Counts the number of nb's */
 497   int nnb=0;
 498
 499   if (d_len(d,natoms,i,j) == 0.0) nnb++;
 500   if (d_len(d,natoms,i,k) == 0.0) nnb++;
 501   if (d_len(d,natoms,i,l) == 0.0) nnb++;
 502   if (d_len(d,natoms,j,k) == 0.0) nnb++;
 503   if (d_len(d,natoms,j,l) == 0.0) nnb++;
 504   if (d_len(d,natoms,k,l) == 0.0) nnb++;
 505
 506   return nnb;
 507 }
 508
 509 void sort_atoms(int natoms,t_dist *d,int *ATMS)
 510 {
 511   /* Sorts atoms. The nb-atoms end up in ATMS[0] and ATMS[3]. */
 512
 513   if (d_len(d,natoms,ATMS[0],ATMS[1]) == 0.0)
 514     swap(&ATMS[1],&ATMS[3]);
 515   else if (d_len(d,natoms,ATMS[0],ATMS[2]) == 0.0)
 516     swap(&ATMS[2],&ATMS[3]);
 517   else if (d_len(d,natoms,ATMS[1],ATMS[2]) == 0.0) {
 518     swap(&ATMS[0],&ATMS[1]);
 519     swap(&ATMS[2],&ATMS[3]);
 520   }
 521   else if (d_len(d,natoms,ATMS[1],ATMS[3]) == 0.0)
 522     swap(&ATMS[1],&ATMS[0]);
 523   else if (d_len(d,natoms,ATMS[2],ATMS[3]) == 0.0)
 524     swap(&ATMS[2],&ATMS[0]);
 525
 526   /* Put the two middle ones in order (ambivalent procedure, necessary?) */
 527   if (d_len(d,natoms,ATMS[0],ATMS[1]) > d_len(d,natoms,ATMS[0],ATMS[2]))
 528     swap(&ATMS[1],&ATMS[2]);
 529 }
 530
 531 real solve_tetrangle(real dAB,real dAC,real dBC,real dBD,real dCD,int cosphi)
 532 {
 533   /* Solves tetrangle eq. for distance AD.
 534      cosphi=cos(phi), where phi is the dihedral angle
 535      cosphi=1 corresponds to min, cosphi=-1 to max
 536      eq. solved and optimized by Maple, watch out for bugs */
 537
 538   real t1,t2,t3,t4,t5,t7,t8,t9,t11,t12,t13,t14,t15,t17,t20;
 539
 540   t1 = dAB*dAB;
 541   t2 = dBC*dBC;
 542   t3 = t1*t2;
 543   t4 = t1*t1;
 544   t5 = dAC*dAC;
 545   t7 = t2*t2;
 546   t8 = t5*t2;
 547   t9 = t5*t5;
 548   t11 = dCD*dCD;
 549   t12 = t2*t11;
 550   t13 = dBD*dBD;
 551   t14 = t2*t13;
 552   t15 = t11*t11;
 553   t17 = t13*t13;
 554   t20 = mysqrt((-2.0*t3+t4-2.0*t1*t5+t7-2.0*t8+t9)*(-2.0*t12+t7-2.0*t14+t15
 555                                                     -2.0*t11*t13+t17));
 556   return mysqrt(-(cosphi*t20-t8-t14+t7-t3-t1*t11+t1*t13-t12+t5*t11-t5*t13)/(2*t2));
 557   /* t20 = mysqrt((-2.0*(t3+t1*t5+t8)     + t4 + t7  + t9)*
 558      (-2.0*(t12+t14+t11*t13) + t7 + t15 + t17));
 559      return mysqrt(-0.5*(cosphi*t20-t8-t14+t7-t3-t1*t11+t1*t13-t12+t5*t11-t5*t13))/dBC;*/
 560 }
 561
 562
 563 int tetrangle_bound(int *ATMS,int cosphi,int natoms,t_dist *d,real tol)
 564 {
 565   int sol=0,D1,D2,D3,D4,D5,nsmooth=0,dmin,dmax;
 566   real dAD[32],dAB,dAC,dBC,dBD,dCD,present_lb,present_ub,present_len;
 567
 568   /*fprintf(stderr,"entering tetrangle_bound\n");*/
 569
 570   /* Try all 32 combinations of bounds */
 571   for (D1=0; (D1<2); D1++) {
 572     for (D2=0; (D2<2); D2++) {
 573       for (D3=0; (D3<2); D3++) {
 574         for (D4=0; (D4<2); D4++) {
 575           for (D5=0; (D5<2); D5++) {
 576             dAB = (D1) ? d_lb(d,natoms,ATMS[0],ATMS[1]):
 577               d_ub(d,natoms,ATMS[0],ATMS[1]);
 578             dAC = (D2) ? d_lb(d,natoms,ATMS[0],ATMS[2]):
 579               d_ub(d,natoms,ATMS[0],ATMS[2]);
 580             dBC = (D3) ? d_lb(d,natoms,ATMS[1],ATMS[2]):
 581               d_ub(d,natoms,ATMS[1],ATMS[2]);
 582             dBD = (D4) ? d_lb(d,natoms,ATMS[1],ATMS[3]):
 583               d_ub(d,natoms,ATMS[1],ATMS[3]);
 584             dCD = (D5) ? d_lb(d,natoms,ATMS[2],ATMS[3]):
 585               d_ub(d,natoms,ATMS[2],ATMS[3]);
 586
 587             dAD[sol] = solve_tetrangle(dAB,dAC,dBC,dBD,dCD,cosphi);
 588             if (debug)
 589               fprintf(debug,"dAD[%d]=%10.5f\n",sol,dAD[sol]);
 590             sol++;
 591           }
 592         }
 593       }
 594     }
 595   }
 596
 597   /* we need these to re-set one of the bounds for the distance */
 598   present_len = d_len(d,natoms,ATMS[0],ATMS[3]);
 599   present_lb  = d_lb(d,natoms,ATMS[0],ATMS[3]);
 600   present_ub  = d_ub(d,natoms,ATMS[0],ATMS[3]);
 601
 602   /* Set the new bound(s) */
 603   minmax(dAD,32,&dmin,&dmax);
 604   if (debug)
 605     fprintf(debug,"Max=dAD[%d] (%12g), Min=dAD[%d] (%12g)\n",
 606             dmax,dAD[dmax],dmin,dAD[dmin]);
 607   /* Set new upper limit  */
 608   if ((cosphi == -1)  &&  ((present_ub - dAD[dmax]) > tol) ) {
 609     set_dist(d,natoms,ATMS[0],ATMS[3],present_lb,dAD[dmax],present_len);
 610     if (debug)
 611       fprintf(debug,"Corrected %d-%d-%d-%d ub to %10.5f\n",ATMS[0]+1,
 612               ATMS[1]+1,ATMS[2]+1,ATMS[3]+1,dAD[dmax]);
 613     nsmooth++;
 614   }
 615   /* Set new lower limit  */
 616   else if ((cosphi == 1)  &&  ( tol < (dAD[dmin] - present_lb))) {
 617     set_dist(d,natoms,ATMS[0],ATMS[3],dAD[dmin],present_ub,present_len);
 618     if (debug)
 619       fprintf(debug,"Corrected %d-%d-%d-%d lb to %10.5f\n",ATMS[0]+1,
 620               ATMS[1]+1,ATMS[2]+1,ATMS[3]+1,dAD[dmin]);
 621     nsmooth++;
 622   }
 623
 624   /* Check the new bounds */
 625   if (d_ub(d,natoms,ATMS[0],ATMS[3]) < d_lb(d,natoms,ATMS[0],ATMS[3]))
 626     fatal_error(0,"Fatal error for tetr. smooth of (%d,%d) ub=%f, lb=%f",
 627                 ATMS[0],ATMS[3],d_ub(d,natoms,ATMS[0],ATMS[3]),
 628                 d_lb(d,natoms,ATMS[0],ATMS[3]));
 629
 630   return nsmooth;
 631 }
 632
 633
 634
 635 int tetrangle (t_dist *d,int natoms,real tol)
 636 {
 637
 638   int i,j,k,l,AT[4],nsmooth=0,nnonbonds1=0,nnonbonds2=0;
 639   bool TLL=FALSE,TUL=FALSE;
 640
 641   fprintf(stderr,"Starting tetrangle smoothing\n");
 642   /* Loops over all sets of four atoms */
 643   for (i=0;i<(natoms-3);i++) {
 644     for (j=i+1;j<(natoms-2);j++) {
 645       if (!dist_set(d,natoms,i,j)) continue;
 646       for (k=j+1;k<(natoms-1);k++) {
 647         if (!dist_set(d,natoms,i,k)) continue;
 648         if (!dist_set(d,natoms,j,k)) continue;
 649         nnonbonds1 = 0;
 650         if (d_len(d,natoms,i,j) == 0.0) nnonbonds1++;
 651         if (d_len(d,natoms,i,k) == 0.0) nnonbonds1++;
 652         if (d_len(d,natoms,j,k) == 0.0) nnonbonds1++;
 653         if (nnonbonds1 > 1) continue;
 654         for (l=k+1;l<natoms;l++) {
 655           if (!dist_set(d,natoms,i,l)) continue;
 656           if (!dist_set(d,natoms,j,l)) continue;
 657           if (!dist_set(d,natoms,k,l)) continue;
 658           nnonbonds2 = 0;
 659           if (d_len(d,natoms,i,l) == 0.0) nnonbonds2++;
 660           if (d_len(d,natoms,j,l) == 0.0) nnonbonds2++;
 661           if (d_len(d,natoms,k,l) == 0.0) nnonbonds2++;
 662           if ( (nnonbonds1+nnonbonds2) != 1) continue;
 663
 664           /*fprintf(stderr,"Trying %d-%d-%d-%d\n",i+1,j+1,k+1,l+1);*/
 665
 666           /* The two following subroutines have been nested into the loops
 667              above. Doesn't look as good, but should be substantially
 668              faster. */
 669           /*      if (!alldistset(i,j,k,l,natoms,d)) continue; */
 670           /*       if (nr_nb(i,j,k,l,natoms,d) != 1)  continue; */
 671
 672           /* Copy indices to array AT for easier handling & sort them */
 673           AT[0] = i;
 674           AT[1] = j;
 675           AT[2] = k;
 676           AT[3] = l;
 677           /*fprintf(stderr,"OK'd atoms    (%d,%d,%d,%d)\n",i+1,j+1,k+1,l+1);*/
 678           sort_atoms(natoms,d,AT);
 679           /*fprintf(stderr,"After sorting (%d,%d,%d,%d)\n",AT[0]+1,AT[1]+1,AT[2]+1,AT[3]+1);*/
 680
 681           /* Are the distances limited by the triangle-ineq.? */
 682           triangle_limited(AT,natoms,d,&TUL,&TLL);
 683           /*if (TUL) fprintf(stderr,"TUL=true\n");
 684             else fprintf(stderr,"TUL=false\n");
 685             if (TLL) fprintf(stderr,"TLL=true\n");
 686             else fprintf(stderr,"TLL=false\n");   */
 687
 688           /* If not, correct bounds according to tetrangle-ineq. */
 689           /* upper limit */
 690           if (!TUL) nsmooth += tetrangle_bound(AT,-1,natoms,d,tol);
 691           /* lower limit */
 692           if (!TLL) nsmooth += tetrangle_bound(AT,1,natoms,d,tol);
 693
 694           /* Comment: If we were to calculate upper and lower bound
 695              simultaneously we could make the code a bit faster     */
 696         }
 697       }
 698     }
 699   }
 700   return nsmooth;
 701 }
 702
 703 /* END OF FUNCTIONS >TETRANGLE-SMOOTH< */
 704
 705
 706
 707 void do_smooth (t_dist *d,t_atoms *atoms,real tol)
 708 {
 709   /* Triangular and tetragonal bound smoothing.
 710      Do we need to nestle triangular and tetragonal
 711      smoothing or can we do them separately? THINK! */
 712
 713   /* In this case I think we can get away with this, because
 714      we are not using distances that we correct with tetragonal
 715      ineq. to smooth other distances. The triangle smoothing
 716      should then spread the tighter bounds to other distances,
 717      not corrected by the teragonal ineq.
 718      (Of course, if you correct one distance by tetragonal ineq.
 719      it will become shorter, and other distances that are involved
 720      in triagonal ineq. with it will be affected.)
 721   */
 722
 723   /* It should only do 2 loops, and in the second no tetragonal
 724      correction will be made. This is because (in the current version)
 725      we only correct nb's that are connected by bonded distances, and
 726      these are unlikely to change during triagonal smoothing. */
 727
 728   /* If this turns out to be true => SKIP DO-LOOP AND JUST HAVE
 729      TRIAGONAL-TETRAGONAL-TRIAGONAL */
 730
 731   int ntri;
 732   int nloops = 0,nsmooth,natoms = atoms->nr;
 733
 734   if (debug)
 735     fprintf(debug,"Starting tetragonal smoothing routine\n");
 736
 737   /* ntri = do_triangle(d,atoms,tol);
 738      do {
 739      nloops += 1;
 740
 741      nsmooth = tetrangle (d,natoms,tol);
 742      fprintf(stderr,"Smoothed %d distances with tetr. ineq.\n",nsmooth);
 743      if (nsmooth > 0)
 744      ntri = do_triangle(d,atoms,tol);
 745      else
 746      ntri = 0;
 747
 748      fprintf(stderr,"Finished %d loop(s) of smoothing.\n",nloops);
 749      } while (ntri > 0);
 750   */
 751   ntri = do_triangle(d,atoms,tol);
 752   nsmooth = tetrangle (d,natoms,tol);
 753   fprintf(stderr,"Smoothed %d distances with tetr. ineq.\n",nsmooth);
 754   ntri = do_triangle(d,atoms,tol);
 755
 756   fprintf(stderr,"Finished smoothing.\n");
 757 }