src/gromacs/trajectoryanalysis/modules/nsc.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-2007, The GROMACS development team.
   6  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "config.h"
  40
  41 #include <stdio.h>
  42 #include <string.h>
  43 #include <stdlib.h>
  44 #include <math.h>
  45 #include <stdarg.h>
  46
  47 /* Modified DvdS */
  48 #include "gromacs/pbcutil/pbc.h"
  49 #include "gromacs/legacyheaders/macros.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/utility/smalloc.h"
  52 #include "nsc.h"
  53
  54 #define TEST_NSC 0
  55
  56 #define TEST_ARC 0
  57 #define TEST_DOD 0
  58 #define TEST_CUBE 0
  59
  60 #define UNSP_ICO_DOD      9
  61 #define UNSP_ICO_ARC     10
  62
  63 real   *xpunsp = NULL;
  64 real    del_cube;
  65 int     n_dot, ico_cube, last_n_dot = 0, last_densit = 0, last_unsp = 0;
  66 int     last_cubus = 0;
  67
  68 #define FOURPI (4.*M_PI)
  69 #define TORAD(A)     ((A)*0.017453293)
  70 #define DP_TOL     0.001
  71
  72 #define UPDATE_FL  __file__ = __FILE__, __line__ = __LINE__
  73 const char * __file__;   /* declared versions of macros */
  74 int          __line__;   /* __FILE__  and __LINE__ */
  75
  76 #ifdef ERROR
  77 #undef ERROR
  78 #endif
  79 #define ERROR UPDATE_FL, error
  80 void error(const char *fmt, ...)
  81 {
  82     va_list args;
  83     fprintf(stderr,
  84             "\n---> ERROR when executing line %i in file %s !\n",
  85             __line__, __file__);
  86     va_start(args, fmt);
  87     vfprintf(stderr, fmt, args);
  88     va_end(args);
  89     fprintf(stderr, "\n---> Execution stopped !\n\n");
  90 }
  91
  92 #define WARNING UPDATE_FL, warning2
  93 void warning2(const char *fmt, ...)
  94 {
  95     va_list args;
  96     fprintf(stderr,
  97             "\n---> WARNING : line %i in file %s\n",
  98             __line__, __file__);
  99     va_start(args, fmt);
 100     vfprintf(stderr, fmt, args);
 101     va_end(args);
 102     fprintf(stderr, " ...!\n\n");
 103     fflush(stderr);
 104     fflush(stdout);
 105 }
 106
 107 #define ASIN safe_asin
 108 real safe_asin(real f)
 109 {
 110     if ( (fabs(f) < 1.00) )
 111     {
 112         return( asin(f) );
 113     }
 114     if ( (fabs(f) - 1.00)  <= DP_TOL)
 115     {
 116         ERROR("ASIN : invalid argument %f", f);
 117     }
 118     return(M_PI_2);
 119 }
 120
 121
 122
 123
 124 /* routines for dot distributions on the surface of the unit sphere */
 125 real rg, rh;
 126
 127 void icosaeder_vertices(real *xus)
 128 {
 129     rh = sqrt(1.-2.*cos(TORAD(72.)))/(1.-cos(TORAD(72.)));
 130     rg = cos(TORAD(72.))/(1.-cos(TORAD(72.)));
 131     /* icosaeder vertices */
 132     xus[ 0] = 0.;                  xus[ 1] = 0.;                  xus[ 2] = 1.;
 133     xus[ 3] = rh*cos(TORAD(72.));  xus[ 4] = rh*sin(TORAD(72.));  xus[ 5] = rg;
 134     xus[ 6] = rh*cos(TORAD(144.)); xus[ 7] = rh*sin(TORAD(144.)); xus[ 8] = rg;
 135     xus[ 9] = rh*cos(TORAD(216.)); xus[10] = rh*sin(TORAD(216.)); xus[11] = rg;
 136     xus[12] = rh*cos(TORAD(288.)); xus[13] = rh*sin(TORAD(288.)); xus[14] = rg;
 137     xus[15] = rh;                  xus[16] = 0;                   xus[17] = rg;
 138     xus[18] = rh*cos(TORAD(36.));  xus[19] = rh*sin(TORAD(36.));  xus[20] = -rg;
 139     xus[21] = rh*cos(TORAD(108.)); xus[22] = rh*sin(TORAD(108.)); xus[23] = -rg;
 140     xus[24] = -rh;                 xus[25] = 0;                   xus[26] = -rg;
 141     xus[27] = rh*cos(TORAD(252.)); xus[28] = rh*sin(TORAD(252.)); xus[29] = -rg;
 142     xus[30] = rh*cos(TORAD(324.)); xus[31] = rh*sin(TORAD(324.)); xus[32] = -rg;
 143     xus[33] = 0.;                  xus[34] = 0.;                  xus[35] = -1.;
 144 }
 145
 146
 147 void divarc(real x1, real y1, real z1,
 148             real x2, real y2, real z2,
 149             int div1, int div2, real *xr, real *yr, real *zr)
 150 {
 151
 152     real xd, yd, zd, dd, d1, d2, s, x, y, z;
 153     real phi, sphi, cphi;
 154
 155     xd = y1*z2-y2*z1;
 156     yd = z1*x2-z2*x1;
 157     zd = x1*y2-x2*y1;
 158     dd = sqrt(xd*xd+yd*yd+zd*zd);
 159     if (dd < DP_TOL)
 160     {
 161         ERROR("divarc: rotation axis of length %f", dd);
 162     }
 163
 164     d1 = x1*x1+y1*y1+z1*z1;
 165     if (d1 < 0.5)
 166     {
 167         ERROR("divarc: vector 1 of sq.length %f", d1);
 168     }
 169     d2 = x2*x2+y2*y2+z2*z2;
 170     if (d2 < 0.5)
 171     {
 172         ERROR("divarc: vector 2 of sq.length %f", d2);
 173     }
 174
 175     phi  = ASIN(dd/sqrt(d1*d2));
 176     phi  = phi*((real)div1)/((real)div2);
 177     sphi = sin(phi); cphi = cos(phi);
 178     s    = (x1*xd+y1*yd+z1*zd)/dd;
 179
 180     x   = xd*s*(1.-cphi)/dd + x1 * cphi + (yd*z1-y1*zd)*sphi/dd;
 181     y   = yd*s*(1.-cphi)/dd + y1 * cphi + (zd*x1-z1*xd)*sphi/dd;
 182     z   = zd*s*(1.-cphi)/dd + z1 * cphi + (xd*y1-x1*yd)*sphi/dd;
 183     dd  = sqrt(x*x+y*y+z*z);
 184     *xr = x/dd; *yr = y/dd; *zr = z/dd;
 185 }
 186
 187 int ico_dot_arc(int densit)   /* densit...required dots per unit sphere */
 188 {
 189     /* dot distribution on a unit sphere based on an icosaeder *
 190      * great circle average refining of icosahedral face       */
 191
 192     int   i, j, k, tl, tl2, tn, tess;
 193     real  a, d, x, y, z, x2, y2, z2, x3, y3, z3;
 194     real  xij, yij, zij, xji, yji, zji, xik, yik, zik, xki, yki, zki,
 195           xjk, yjk, zjk, xkj, ykj, zkj;
 196     real *xus = NULL;
 197
 198     /* calculate tessalation level */
 199     a     = sqrt((((real) densit)-2.)/10.);
 200     tess  = (int) ceil(a);
 201     n_dot = 10*tess*tess+2;
 202     if (n_dot < densit)
 203     {
 204         ERROR("ico_dot_arc: error in formula for tessalation level (%d->%d, %d)",
 205               tess, n_dot, densit);
 206     }
 207
 208     snew(xus, 3*n_dot);
 209     xpunsp = xus;
 210     icosaeder_vertices(xus);
 211
 212     if (tess > 1)
 213     {
 214         tn = 12;
 215         a  = rh*rh*2.*(1.-cos(TORAD(72.)));
 216         /* calculate tessalation of icosaeder edges */
 217         for (i = 0; i < 11; i++)
 218         {
 219             for (j = i+1; j < 12; j++)
 220             {
 221                 x = xus[3*i]-xus[3*j];
 222                 y = xus[1+3*i]-xus[1+3*j]; z = xus[2+3*i]-xus[2+3*j];
 223                 d = x*x+y*y+z*z;
 224                 if (fabs(a-d) > DP_TOL)
 225                 {
 226                     continue;
 227                 }
 228                 for (tl = 1; tl < tess; tl++)
 229                 {
 230                     if (tn >= n_dot)
 231                     {
 232                         ERROR("ico_dot: tn exceeds dimension of xus");
 233                     }
 234                     divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
 235                            xus[3*j], xus[1+3*j], xus[2+3*j],
 236                            tl, tess, &xus[3*tn], &xus[1+3*tn], &xus[2+3*tn]);
 237                     tn++;
 238                 }
 239             }
 240         }
 241         /* calculate tessalation of icosaeder faces */
 242         for (i = 0; i < 10; i++)
 243         {
 244             for (j = i+1; j < 11; j++)
 245             {
 246                 x = xus[3*i]-xus[3*j];
 247                 y = xus[1+3*i]-xus[1+3*j]; z = xus[2+3*i]-xus[2+3*j];
 248                 d = x*x+y*y+z*z;
 249                 if (fabs(a-d) > DP_TOL)
 250                 {
 251                     continue;
 252                 }
 253
 254                 for (k = j+1; k < 12; k++)
 255                 {
 256                     x = xus[3*i]-xus[3*k];
 257                     y = xus[1+3*i]-xus[1+3*k]; z = xus[2+3*i]-xus[2+3*k];
 258                     d = x*x+y*y+z*z;
 259                     if (fabs(a-d) > DP_TOL)
 260                     {
 261                         continue;
 262                     }
 263                     x = xus[3*j]-xus[3*k];
 264                     y = xus[1+3*j]-xus[1+3*k]; z = xus[2+3*j]-xus[2+3*k];
 265                     d = x*x+y*y+z*z;
 266                     if (fabs(a-d) > DP_TOL)
 267                     {
 268                         continue;
 269                     }
 270                     for (tl = 1; tl < tess-1; tl++)
 271                     {
 272                         divarc(xus[3*j], xus[1+3*j], xus[2+3*j],
 273                                xus[3*i], xus[1+3*i], xus[2+3*i],
 274                                tl, tess, &xji, &yji, &zji);
 275                         divarc(xus[3*k], xus[1+3*k], xus[2+3*k],
 276                                xus[3*i], xus[1+3*i], xus[2+3*i],
 277                                tl, tess, &xki, &yki, &zki);
 278
 279                         for (tl2 = 1; tl2 < tess-tl; tl2++)
 280                         {
 281                             divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
 282                                    xus[3*j], xus[1+3*j], xus[2+3*j],
 283                                    tl2, tess, &xij, &yij, &zij);
 284                             divarc(xus[3*k], xus[1+3*k], xus[2+3*k],
 285                                    xus[3*j], xus[1+3*j], xus[2+3*j],
 286                                    tl2, tess, &xkj, &ykj, &zkj);
 287                             divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
 288                                    xus[3*k], xus[1+3*k], xus[2+3*k],
 289                                    tess-tl-tl2, tess, &xik, &yik, &zik);
 290                             divarc(xus[3*j], xus[1+3*j], xus[2+3*j],
 291                                    xus[3*k], xus[1+3*k], xus[2+3*k],
 292                                    tess-tl-tl2, tess, &xjk, &yjk, &zjk);
 293                             if (tn >= n_dot)
 294                             {
 295                                 ERROR("ico_dot: tn exceeds dimension of xus");
 296                             }
 297                             divarc(xki, yki, zki, xji, yji, zji, tl2, tess-tl,
 298                                    &x, &y, &z);
 299                             divarc(xkj, ykj, zkj, xij, yij, zij, tl, tess-tl2,
 300                                    &x2, &y2, &z2);
 301                             divarc(xjk, yjk, zjk, xik, yik, zik, tl, tl+tl2,
 302                                    &x3, &y3, &z3);
 303                             x           = x+x2+x3; y = y+y2+y3; z = z+z2+z3;
 304                             d           = sqrt(x*x+y*y+z*z);
 305                             xus[3*tn]   = x/d;
 306                             xus[1+3*tn] = y/d;
 307                             xus[2+3*tn] = z/d;
 308                             tn++;
 309                         } /* cycle tl2 */
 310                     }     /* cycle tl */
 311                 }         /* cycle k */
 312             }             /* cycle j */
 313         }                 /* cycle i */
 314         if (n_dot != tn)
 315         {
 316             ERROR("ico_dot: n_dot(%d) and tn(%d) differ", n_dot, tn);
 317         }
 318     }   /* end of if (tess > 1) */
 319
 320     return n_dot;
 321 }                           /* end of routine ico_dot_arc */
 322
 323 int ico_dot_dod(int densit) /* densit...required dots per unit sphere */
 324 {
 325     /* dot distribution on a unit sphere based on an icosaeder *
 326      * great circle average refining of icosahedral face       */
 327
 328     int   i, j, k, tl, tl2, tn, tess, j1, j2;
 329     real  a, d, x, y, z, x2, y2, z2, x3, y3, z3, ai_d, adod;
 330     real  xij, yij, zij, xji, yji, zji, xik, yik, zik, xki, yki, zki,
 331           xjk, yjk, zjk, xkj, ykj, zkj;
 332     real *xus = NULL;
 333     /* calculate tesselation level */
 334     a     = sqrt((((real) densit)-2.)/30.);
 335     tess  = max((int) ceil(a), 1);
 336     n_dot = 30*tess*tess+2;
 337     if (n_dot < densit)
 338     {
 339         ERROR("ico_dot_dod: error in formula for tessalation level (%d->%d, %d)",
 340               tess, n_dot, densit);
 341     }
 342
 343     snew(xus, 3*n_dot);
 344     xpunsp = xus;
 345     icosaeder_vertices(xus);
 346
 347     tn = 12;
 348     /* square of the edge of an icosaeder */
 349     a = rh*rh*2.*(1.-cos(TORAD(72.)));
 350     /* dodecaeder vertices */
 351     for (i = 0; i < 10; i++)
 352     {
 353         for (j = i+1; j < 11; j++)
 354         {
 355             x = xus[3*i]-xus[3*j];
 356             y = xus[1+3*i]-xus[1+3*j]; z = xus[2+3*i]-xus[2+3*j];
 357             d = x*x+y*y+z*z;
 358             if (fabs(a-d) > DP_TOL)
 359             {
 360                 continue;
 361             }
 362             for (k = j+1; k < 12; k++)
 363             {
 364                 x = xus[3*i]-xus[3*k];
 365                 y = xus[1+3*i]-xus[1+3*k]; z = xus[2+3*i]-xus[2+3*k];
 366                 d = x*x+y*y+z*z;
 367                 if (fabs(a-d) > DP_TOL)
 368                 {
 369                     continue;
 370                 }
 371                 x = xus[3*j]-xus[3*k];
 372                 y = xus[1+3*j]-xus[1+3*k]; z = xus[2+3*j]-xus[2+3*k];
 373                 d = x*x+y*y+z*z;
 374                 if (fabs(a-d) > DP_TOL)
 375                 {
 376                     continue;
 377                 }
 378                 x         = xus[  3*i]+xus[  3*j]+xus[  3*k];
 379                 y         = xus[1+3*i]+xus[1+3*j]+xus[1+3*k];
 380                 z         = xus[2+3*i]+xus[2+3*j]+xus[2+3*k];
 381                 d         = sqrt(x*x+y*y+z*z);
 382                 xus[3*tn] = x/d; xus[1+3*tn] = y/d; xus[2+3*tn] = z/d;
 383                 tn++;
 384             }
 385         }
 386     }
 387
 388     if (tess > 1)
 389     {
 390         tn = 32;
 391         /* square of the edge of an dodecaeder */
 392         adod = 4.*(cos(TORAD(108.))-cos(TORAD(120.)))/(1.-cos(TORAD(120.)));
 393         /* square of the distance of two adjacent vertices of ico- and dodecaeder */
 394         ai_d = 2.*(1.-sqrt(1.-a/3.));
 395
 396         /* calculate tessalation of mixed edges */
 397         for (i = 0; i < 31; i++)
 398         {
 399             j1 = 12; j2 = 32; a = ai_d;
 400             if (i >= 12)
 401             {
 402                 j1 = i+1; a = adod;
 403             }
 404             for (j = j1; j < j2; j++)
 405             {
 406                 x = xus[3*i]-xus[3*j];
 407                 y = xus[1+3*i]-xus[1+3*j]; z = xus[2+3*i]-xus[2+3*j];
 408                 d = x*x+y*y+z*z;
 409                 if (fabs(a-d) > DP_TOL)
 410                 {
 411                     continue;
 412                 }
 413                 for (tl = 1; tl < tess; tl++)
 414                 {
 415                     if (tn >= n_dot)
 416                     {
 417                         ERROR("ico_dot: tn exceeds dimension of xus");
 418                     }
 419                     divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
 420                            xus[3*j], xus[1+3*j], xus[2+3*j],
 421                            tl, tess, &xus[3*tn], &xus[1+3*tn], &xus[2+3*tn]);
 422                     tn++;
 423                 }
 424             }
 425         }
 426         /* calculate tessalation of pentakisdodecahedron faces */
 427         for (i = 0; i < 12; i++)
 428         {
 429             for (j = 12; j < 31; j++)
 430             {
 431                 x = xus[3*i]-xus[3*j];
 432                 y = xus[1+3*i]-xus[1+3*j]; z = xus[2+3*i]-xus[2+3*j];
 433                 d = x*x+y*y+z*z;
 434                 if (fabs(ai_d-d) > DP_TOL)
 435                 {
 436                     continue;
 437                 }
 438
 439                 for (k = j+1; k < 32; k++)
 440                 {
 441                     x = xus[3*i]-xus[3*k];
 442                     y = xus[1+3*i]-xus[1+3*k]; z = xus[2+3*i]-xus[2+3*k];
 443                     d = x*x+y*y+z*z;
 444                     if (fabs(ai_d-d) > DP_TOL)
 445                     {
 446                         continue;
 447                     }
 448                     x = xus[3*j]-xus[3*k];
 449                     y = xus[1+3*j]-xus[1+3*k]; z = xus[2+3*j]-xus[2+3*k];
 450                     d = x*x+y*y+z*z;
 451                     if (fabs(adod-d) > DP_TOL)
 452                     {
 453                         continue;
 454                     }
 455                     for (tl = 1; tl < tess-1; tl++)
 456                     {
 457                         divarc(xus[3*j], xus[1+3*j], xus[2+3*j],
 458                                xus[3*i], xus[1+3*i], xus[2+3*i],
 459                                tl, tess, &xji, &yji, &zji);
 460                         divarc(xus[3*k], xus[1+3*k], xus[2+3*k],
 461                                xus[3*i], xus[1+3*i], xus[2+3*i],
 462                                tl, tess, &xki, &yki, &zki);
 463
 464                         for (tl2 = 1; tl2 < tess-tl; tl2++)
 465                         {
 466                             divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
 467                                    xus[3*j], xus[1+3*j], xus[2+3*j],
 468                                    tl2, tess, &xij, &yij, &zij);
 469                             divarc(xus[3*k], xus[1+3*k], xus[2+3*k],
 470                                    xus[3*j], xus[1+3*j], xus[2+3*j],
 471                                    tl2, tess, &xkj, &ykj, &zkj);
 472                             divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
 473                                    xus[3*k], xus[1+3*k], xus[2+3*k],
 474                                    tess-tl-tl2, tess, &xik, &yik, &zik);
 475                             divarc(xus[3*j], xus[1+3*j], xus[2+3*j],
 476                                    xus[3*k], xus[1+3*k], xus[2+3*k],
 477                                    tess-tl-tl2, tess, &xjk, &yjk, &zjk);
 478                             if (tn >= n_dot)
 479                             {
 480                                 ERROR("ico_dot: tn exceeds dimension of xus");
 481                             }
 482                             divarc(xki, yki, zki, xji, yji, zji, tl2, tess-tl,
 483                                    &x, &y, &z);
 484                             divarc(xkj, ykj, zkj, xij, yij, zij, tl, tess-tl2,
 485                                    &x2, &y2, &z2);
 486                             divarc(xjk, yjk, zjk, xik, yik, zik, tl, tl+tl2,
 487                                    &x3, &y3, &z3);
 488                             x           = x+x2+x3; y = y+y2+y3; z = z+z2+z3;
 489                             d           = sqrt(x*x+y*y+z*z);
 490                             xus[3*tn]   = x/d;
 491                             xus[1+3*tn] = y/d;
 492                             xus[2+3*tn] = z/d;
 493                             tn++;
 494                         } /* cycle tl2 */
 495                     }     /* cycle tl */
 496                 }         /* cycle k */
 497             }             /* cycle j */
 498         }                 /* cycle i */
 499         if (n_dot != tn)
 500         {
 501             ERROR("ico_dot: n_dot(%d) and tn(%d) differ", n_dot, tn);
 502         }
 503     }   /* end of if (tess > 1) */
 504
 505     return n_dot;
 506 }       /* end of routine ico_dot_dod */
 507
 508 int unsp_type(int densit)
 509 {
 510     int i1, i2;
 511     i1 = 1;
 512     while (10*i1*i1+2 < densit)
 513     {
 514         i1++;
 515     }
 516     i2 = 1;
 517     while (30*i2*i2+2 < densit)
 518     {
 519         i2++;
 520     }
 521     if (10*i1*i1-2 < 30*i2*i2-2)
 522     {
 523         return UNSP_ICO_ARC;
 524     }
 525     else
 526     {
 527         return UNSP_ICO_DOD;
 528     }
 529 }
 530
 531 int make_unsp(int densit, int mode, int * num_dot, int cubus)
 532 {
 533     int  *ico_wk, *ico_pt;
 534     int   ndot, ico_cube_cb, i, j, k, l, ijk, tn, tl, tl2;
 535     real *xus;
 536     int  *work;
 537     real  x, y, z;
 538
 539     if (xpunsp)
 540     {
 541         free(xpunsp);
 542     }
 543
 544     k = 1; if (mode < 0)
 545     {
 546         k = 0; mode = -mode;
 547     }
 548     if (mode == UNSP_ICO_ARC)
 549     {
 550         ndot = ico_dot_arc(densit);
 551     }
 552     else if (mode == UNSP_ICO_DOD)
 553     {
 554         ndot = ico_dot_dod(densit);
 555     }
 556     else
 557     {
 558         WARNING("make_unsp: mode %c%d not allowed", (k) ? '+' : '-', mode);
 559         return 1;
 560     }
 561
 562     last_n_dot = ndot; last_densit = densit; last_unsp = mode;
 563     *num_dot   = ndot; if (k)
 564     {
 565         return 0;
 566     }
 567
 568     /* in the following the dots of the unit sphere may be resorted */
 569     last_unsp = -last_unsp;
 570
 571     /* determine distribution of points in elementary cubes */
 572     if (cubus)
 573     {
 574         ico_cube = cubus;
 575     }
 576     else
 577     {
 578         last_cubus = 0;
 579         i          = 1;
 580         while (i*i*i*2 < ndot)
 581         {
 582             i++;
 583         }
 584         ico_cube = max(i-1, 0);
 585     }
 586     ico_cube_cb = ico_cube*ico_cube*ico_cube;
 587     del_cube    = 2./((real)ico_cube);
 588     snew(work, ndot);
 589     xus = xpunsp;
 590     for (l = 0; l < ndot; l++)
 591     {
 592         i = max((int) floor((1.+xus[3*l])/del_cube), 0);
 593         if (i >= ico_cube)
 594         {
 595             i = ico_cube-1;
 596         }
 597         j = max((int) floor((1.+xus[1+3*l])/del_cube), 0);
 598         if (j >= ico_cube)
 599         {
 600             j = ico_cube-1;
 601         }
 602         k = max((int) floor((1.+xus[2+3*l])/del_cube), 0);
 603         if (k >= ico_cube)
 604         {
 605             k = ico_cube-1;
 606         }
 607         ijk     = i+j*ico_cube+k*ico_cube*ico_cube;
 608         work[l] = ijk;
 609     }
 610
 611     snew(ico_wk, 2*ico_cube_cb+1);
 612
 613     ico_pt = ico_wk+ico_cube_cb;
 614     for (l = 0; l < ndot; l++)
 615     {
 616         ico_wk[work[l]]++; /* dots per elementary cube */
 617     }
 618
 619     /* reordering of the coordinate array in accordance with box number */
 620     tn = 0;
 621     for (i = 0; i < ico_cube; i++)
 622     {
 623         for (j = 0; j < ico_cube; j++)
 624         {
 625             for (k = 0; k < ico_cube; k++)
 626             {
 627                 tl            = 0;
 628                 tl2           = tn;
 629                 ijk           = i+ico_cube*j+ico_cube*ico_cube*k;
 630                 *(ico_pt+ijk) = tn;
 631                 for (l = tl2; l < ndot; l++)
 632                 {
 633                     if (ijk == work[l])
 634                     {
 635                         x          = xus[3*l]; y = xus[1+3*l]; z = xus[2+3*l];
 636                         xus[3*l]   = xus[3*tn];
 637                         xus[1+3*l] = xus[1+3*tn]; xus[2+3*l] = xus[2+3*tn];
 638                         xus[3*tn]  = x; xus[1+3*tn] = y; xus[2+3*tn] = z;
 639                         ijk        = work[l]; work[l] = work[tn]; work[tn] = ijk;
 640                         tn++; tl++;
 641                     }
 642                 }
 643                 *(ico_wk+ijk) = tl;
 644             } /* cycle k */
 645         }     /* cycle j */
 646     }         /* cycle i */
 647     sfree(ico_wk);
 648     sfree(work);
 649
 650     return 0;
 651 }
 652
 653
 654 typedef struct _stwknb {
 655     real x;
 656     real y;
 657     real z;
 658     real dot;
 659 } Neighb;
 660
 661 int nsc_dclm_pbc(const rvec *coords, real *radius, int nat,
 662                  int  densit, int mode,
 663                  real *value_of_area, real **at_area,
 664                  real *value_of_vol,
 665                  real **lidots, int *nu_dots,
 666                  atom_id index[], int ePBC, matrix box)
 667 {
 668     int         iat, i, ii, iii, ix, iy, iz, ixe, ixs, iye, iys, ize, izs, i_ac;
 669     int         jat, j, jj, jjj, jx, jy, jz;
 670     int         distribution;
 671     int         l;
 672     int         maxnei, nnei, last, maxdots = 0;
 673     int        *wkdot = NULL, *wkbox = NULL, *wkat1 = NULL, *wkatm = NULL;
 674     Neighb     *wknb, *ctnb;
 675     int         iii1, iii2, iiat, lfnr = 0, i_at, j_at;
 676     real        dx, dy, dz, dd, ai, aisq, ajsq, aj, as, a;
 677     real        xi, yi, zi, xs = 0., ys = 0., zs = 0.;
 678     real        dotarea, area, vol = 0.;
 679     real       *xus, *dots = NULL, *atom_area = NULL;
 680     int         nxbox, nybox, nzbox, nxy, nxyz;
 681     real        xmin = 0, ymin = 0, zmin = 0, xmax, ymax, zmax, ra2max, d;
 682     const real *pco;
 683     /* Added DvdS 2006-07-19 */
 684     t_pbc       pbc;
 685     rvec        ddx, *x = NULL;
 686     int         iat_xx, jat_xx;
 687
 688     distribution = unsp_type(densit);
 689     if (distribution != -last_unsp || last_cubus != 4 ||
 690         (densit != last_densit && densit != last_n_dot))
 691     {
 692         if (make_unsp(densit, (-distribution), &n_dot, 4))
 693         {
 694             return 1;
 695         }
 696     }
 697     xus = xpunsp;
 698
 699     dotarea = FOURPI/(real) n_dot;
 700     area    = 0.;
 701
 702     if (debug)
 703     {
 704         fprintf(debug, "nsc_dclm: n_dot=%5d %9.3f\n", n_dot, dotarea);
 705     }
 706
 707     /* start with neighbour list */
 708     /* calculate neighbour list with the box algorithm */
 709     if (nat == 0)
 710     {
 711         WARNING("nsc_dclm: no surface atoms selected");
 712         return 1;
 713     }
 714     if (mode & FLAG_VOLUME)
 715     {
 716         vol = 0.;
 717     }
 718     if (mode & FLAG_DOTS)
 719     {
 720         maxdots = (3*n_dot*nat)/10;
 721         /* should be set to NULL on first user call */
 722         if (dots == NULL)
 723         {
 724             snew(dots, maxdots);
 725         }
 726         else
 727         {
 728             srenew(dots, maxdots);
 729         }
 730
 731         lfnr = 0;
 732     }
 733     if (mode & FLAG_ATOM_AREA)
 734     {
 735         /* should be set to NULL on first user call */
 736         if (atom_area == NULL)
 737         {
 738             snew(atom_area, nat);
 739         }
 740         else
 741         {
 742             srenew(atom_area, nat);
 743         }
 744     }
 745
 746     /* Compute minimum size for grid cells */
 747     ra2max = radius[index[0]];
 748     for (iat_xx = 1; (iat_xx < nat); iat_xx++)
 749     {
 750         iat    = index[iat_xx];
 751         ra2max = max(ra2max, radius[iat]);
 752     }
 753     ra2max = 2*ra2max;
 754
 755     /* Added DvdS 2006-07-19 */
 756     /* Updated 2008-10-09 */
 757     if (box)
 758     {
 759         set_pbc(&pbc, ePBC, box);
 760         snew(x, nat);
 761         for (i = 0; (i < nat); i++)
 762         {
 763             iat  = index[0];
 764             copy_rvec(coords[iat], x[i]);
 765         }
 766         put_atoms_in_triclinic_unitcell(ecenterTRIC, box, nat, x);
 767         nxbox = max(1, floor(norm(box[XX])/ra2max));
 768         nybox = max(1, floor(norm(box[YY])/ra2max));
 769         nzbox = max(1, floor(norm(box[ZZ])/ra2max));
 770         if (debug)
 771         {
 772             fprintf(debug, "nbox = %d, %d, %d\n", nxbox, nybox, nzbox);
 773         }
 774     }
 775     else
 776     {
 777         /* dimensions of atomic set, cell edge is 2*ra_max */
 778         iat    = index[0];
 779         xmin   = coords[iat][XX]; xmax = xmin; xs = xmin;
 780         ymin   = coords[iat][YY]; ymax = ymin; ys = ymin;
 781         zmin   = coords[iat][ZZ]; zmax = zmin; zs = zmin;
 782
 783         for (iat_xx = 1; (iat_xx < nat); iat_xx++)
 784         {
 785             iat  = index[iat_xx];
 786             pco  = coords[iat];
 787             xmin = min(xmin, *pco);     xmax = max(xmax, *pco);
 788             ymin = min(ymin, *(pco+1)); ymax = max(ymax, *(pco+1));
 789             zmin = min(zmin, *(pco+2)); zmax = max(zmax, *(pco+2));
 790             xs   = xs+ *pco; ys = ys+ *(pco+1); zs = zs+ *(pco+2);
 791         }
 792         xs = xs/ (real) nat;
 793         ys = ys/ (real) nat;
 794         zs = zs/ (real) nat;
 795         if (debug)
 796         {
 797             fprintf(debug, "nsc_dclm: n_dot=%5d ra2max=%9.3f %9.3f\n", n_dot, ra2max, dotarea);
 798         }
 799
 800         d    = xmax-xmin; nxbox = (int) max(ceil(d/ra2max), 1.);
 801         d    = (((real)nxbox)*ra2max-d)/2.;
 802         xmin = xmin-d; xmax = xmax+d;
 803         d    = ymax-ymin; nybox = (int) max(ceil(d/ra2max), 1.);
 804         d    = (((real)nybox)*ra2max-d)/2.;
 805         ymin = ymin-d; ymax = ymax+d;
 806         d    = zmax-zmin; nzbox = (int) max(ceil(d/ra2max), 1.);
 807         d    = (((real)nzbox)*ra2max-d)/2.;
 808         zmin = zmin-d; zmax = zmax+d;
 809     }
 810     /* Help variables */
 811     nxy  = nxbox*nybox;
 812     nxyz = nxy*nzbox;
 813
 814     /* box number of atoms */
 815     snew(wkatm, nat);
 816     snew(wkat1, nat);
 817     snew(wkdot, n_dot);
 818     snew(wkbox, nxyz+1);
 819
 820     if (box)
 821     {
 822         matrix box_1;
 823         rvec   x_1;
 824         int    ix, iy, iz, m;
 825         m_inv(box, box_1);
 826         for (i = 0; (i < nat); i++)
 827         {
 828             mvmul(box_1, x[i], x_1);
 829             ix = ((int)floor(x_1[XX]*nxbox) + 2*nxbox) % nxbox;
 830             iy = ((int)floor(x_1[YY]*nybox) + 2*nybox) % nybox;
 831             iz = ((int)floor(x_1[ZZ]*nzbox) + 2*nzbox) % nzbox;
 832             j  =  ix + iy*nxbox + iz*nxbox*nybox;
 833             if (debug)
 834             {
 835                 fprintf(debug, "Atom %d cell index %d. x = (%8.3f,%8.3f,%8.3f) fc = (%8.3f,%8.3f,%8.3f)\n",
 836                         i, j, x[i][XX], x[i][YY], x[i][ZZ], x_1[XX], x_1[YY], x_1[ZZ]);
 837             }
 838             range_check(j, 0, nxyz);
 839             wkat1[i] = j;
 840             wkbox[j]++;
 841         }
 842     }
 843     else
 844     {
 845         /* Put the atoms in their boxes */
 846         for (iat_xx = 0; (iat_xx < nat); iat_xx++)
 847         {
 848             iat           = index[iat_xx];
 849             pco           = coords[iat];
 850             i             = (int) max(floor((pco[XX]-xmin)/ra2max), 0);
 851             i             = min(i, nxbox-1);
 852             j             = (int) max(floor((pco[YY]-ymin)/ra2max), 0);
 853             j             = min(j, nybox-1);
 854             l             = (int) max(floor((pco[ZZ]-zmin)/ra2max), 0);
 855             l             = min(l, nzbox-1);
 856             i             = i+j*nxbox+l*nxy;
 857             wkat1[iat_xx] = i;
 858             wkbox[i]++;
 859         }
 860     }
 861
 862     /* sorting of atoms in accordance with box numbers */
 863     j = wkbox[0];
 864     for (i = 1; i < nxyz; i++)
 865     {
 866         j = max(wkbox[i], j);
 867     }
 868     for (i = 1; i <= nxyz; i++)
 869     {
 870         wkbox[i] += wkbox[i-1];
 871     }
 872
 873     /* maxnei = (int) floor(ra2max*ra2max*ra2max*0.5); */
 874     maxnei = min(nat, 27*j);
 875     snew(wknb, maxnei);
 876     for (iat_xx = 0; iat_xx < nat; iat_xx++)
 877     {
 878         iat = index[iat_xx];
 879         range_check(wkat1[iat_xx], 0, nxyz);
 880         wkatm[--wkbox[wkat1[iat_xx]]] = iat_xx;
 881         if (debug)
 882         {
 883             fprintf(debug, "atom %5d on place %5d\n", iat, wkbox[wkat1[iat_xx]]);
 884         }
 885     }
 886
 887     if (debug)
 888     {
 889         fprintf(debug, "nsc_dclm: n_dot=%5d ra2max=%9.3f %9.3f\n",
 890                 n_dot, ra2max, dotarea);
 891         fprintf(debug, "neighbour list calculated/box(xyz):%d %d %d\n",
 892                 nxbox, nybox, nzbox);
 893
 894         for (i = 0; i < nxyz; i++)
 895         {
 896             fprintf(debug, "box %6d : atoms %4d-%4d    %5d\n",
 897                     i, wkbox[i], wkbox[i+1]-1, wkbox[i+1]-wkbox[i]);
 898         }
 899         for (i = 0; i < nat; i++)
 900         {
 901             fprintf(debug, "list place %5d by atom %7d\n", i, index[wkatm[i]]);
 902         }
 903     }
 904
 905     /* calculate surface for all atoms, step cube-wise */
 906     for (iz = 0; iz < nzbox; iz++)
 907     {
 908         iii = iz*nxy;
 909         if (box)
 910         {
 911             izs = iz-1;
 912             ize = min(iz+2, izs+nzbox);
 913         }
 914         else
 915         {
 916             izs = max(iz-1, 0);
 917             ize = min(iz+2, nzbox);
 918         }
 919         for (iy = 0; iy < nybox; iy++)
 920         {
 921             ii = iy*nxbox+iii;
 922             if (box)
 923             {
 924                 iys = iy-1;
 925                 iye = min(iy+2, iys+nybox);
 926             }
 927             else
 928             {
 929                 iys = max(iy-1, 0);
 930                 iye = min(iy+2, nybox);
 931             }
 932             for (ix = 0; ix < nxbox; ix++)
 933             {
 934                 i    = ii+ix;
 935                 iii1 = wkbox[i];
 936                 iii2 = wkbox[i+1];
 937                 if (iii1 >= iii2)
 938                 {
 939                     continue;
 940                 }
 941                 if (box)
 942                 {
 943                     ixs = ix-1;
 944                     ixe = min(ix+2, ixs+nxbox);
 945                 }
 946                 else
 947                 {
 948                     ixs = max(ix-1, 0);
 949                     ixe = min(ix+2, nxbox);
 950                 }
 951                 iiat = 0;
 952                 /* make intermediate atom list */
 953                 for (jz = izs; jz < ize; jz++)
 954                 {
 955                     jjj = ((jz+nzbox) % nzbox)*nxy;
 956                     for (jy = iys; jy < iye; jy++)
 957                     {
 958                         jj = ((jy+nybox) % nybox)*nxbox+jjj;
 959                         for (jx = ixs; jx < ixe; jx++)
 960                         {
 961                             j = jj+((jx+nxbox) % nxbox);
 962                             for (jat = wkbox[j]; jat < wkbox[j+1]; jat++)
 963                             {
 964                                 range_check(wkatm[jat], 0, nat);
 965                                 range_check(iiat, 0, nat);
 966                                 wkat1[iiat] = wkatm[jat];
 967                                 iiat++;
 968                             } /* end of cycle "jat" */
 969                         }     /* end of cycle "jx" */
 970                     }         /* end of cycle "jy" */
 971                 }             /* end of cycle "jz" */
 972                 for (iat = iii1; iat < iii2; iat++)
 973                 {
 974                     i_at = index[wkatm[iat]];
 975                     ai   = radius[i_at];
 976                     aisq = ai*ai;
 977                     pco  = coords[i_at];
 978                     xi   = pco[XX]; yi = pco[YY]; zi = pco[ZZ];
 979                     for (i = 0; i < n_dot; i++)
 980                     {
 981                         wkdot[i] = 0;
 982                     }
 983
 984                     ctnb = wknb; nnei = 0;
 985                     for (j = 0; j < iiat; j++)
 986                     {
 987                         j_at = index[wkat1[j]];
 988                         if (j_at == i_at)
 989                         {
 990                             continue;
 991                         }
 992
 993                         aj   = radius[j_at];
 994                         ajsq = aj*aj;
 995                         pco  = coords[j_at];
 996
 997                         /* Added DvdS 2006-07-19 */
 998                         if (box)
 999                         {
1000                             /*rvec xxi;
1001
1002                                xxi[XX] = xi;
1003                                xxi[YY] = yi;
1004                                xxi[ZZ] = zi;
1005                                pbc_dx(&pbc,pco,xxi,ddx);*/
1006                             pbc_dx(&pbc, coords[j_at], coords[i_at], ddx);
1007                             dx = ddx[XX];
1008                             dy = ddx[YY];
1009                             dz = ddx[ZZ];
1010                         }
1011                         else
1012                         {
1013                             dx = pco[XX]-xi;
1014                             dy = pco[YY]-yi;
1015                             dz = pco[ZZ]-zi;
1016                         }
1017                         dd = dx*dx+dy*dy+dz*dz;
1018                         as = ai+aj;
1019                         if (dd > as*as)
1020                         {
1021                             continue;
1022                         }
1023                         nnei++;
1024                         ctnb->x   = dx;
1025                         ctnb->y   = dy;
1026                         ctnb->z   = dz;
1027                         ctnb->dot = (dd+aisq-ajsq)/(2.*ai); /* reference dot product */
1028                         ctnb++;
1029                     }
1030
1031                     /* check points on accessibility */
1032                     if (nnei)
1033                     {
1034                         last = 0; i_ac = 0;
1035                         for (l = 0; l < n_dot; l++)
1036                         {
1037                             if (xus[3*l]*(wknb+last)->x+
1038                                 xus[1+3*l]*(wknb+last)->y+
1039                                 xus[2+3*l]*(wknb+last)->z <= (wknb+last)->dot)
1040                             {
1041                                 for (j = 0; j < nnei; j++)
1042                                 {
1043                                     if (xus[3*l]*(wknb+j)->x+xus[1+3*l]*(wknb+j)->y+
1044                                         xus[2+3*l]*(wknb+j)->z > (wknb+j)->dot)
1045                                     {
1046                                         last = j;
1047                                         break;
1048                                     }
1049                                 }
1050                                 if (j >= nnei)
1051                                 {
1052                                     i_ac++;
1053                                     wkdot[l] = 1;
1054                                 }
1055                             } /* end of cycle j */
1056                         }     /* end of cycle l */
1057                     }
1058                     else
1059                     {
1060                         i_ac  = n_dot;
1061                         for (l = 0; l < n_dot; l++)
1062                         {
1063                             wkdot[l] = 1;
1064                         }
1065                     }
1066
1067                     if (debug)
1068                     {
1069                         fprintf(debug, "i_ac=%d, dotarea=%8.3f, aisq=%8.3f\n",
1070                                 i_ac, dotarea, aisq);
1071                     }
1072
1073                     a    = aisq*dotarea* (real) i_ac;
1074                     area = area + a;
1075                     if (mode & FLAG_ATOM_AREA)
1076                     {
1077                         range_check(wkatm[iat], 0, nat);
1078                         atom_area[wkatm[iat]] = a;
1079                     }
1080                     if (mode & FLAG_DOTS)
1081                     {
1082                         for (l = 0; l < n_dot; l++)
1083                         {
1084                             if (wkdot[l])
1085                             {
1086                                 lfnr++;
1087                                 if (maxdots <= 3*lfnr+1)
1088                                 {
1089                                     maxdots = maxdots+n_dot*3;
1090                                     srenew(dots, maxdots);
1091                                 }
1092                                 dots[3*lfnr-3] = ai*xus[3*l]+xi;
1093                                 dots[3*lfnr-2] = ai*xus[1+3*l]+yi;
1094                                 dots[3*lfnr-1] = ai*xus[2+3*l]+zi;
1095                             }
1096                         }
1097                     }
1098                     if (mode & FLAG_VOLUME)
1099                     {
1100                         dx = 0.; dy = 0.; dz = 0.;
1101                         for (l = 0; l < n_dot; l++)
1102                         {
1103                             if (wkdot[l])
1104                             {
1105                                 dx = dx+xus[3*l];
1106                                 dy = dy+xus[1+3*l];
1107                                 dz = dz+xus[2+3*l];
1108                             }
1109                         }
1110                         vol = vol+aisq*(dx*(xi-xs)+dy*(yi-ys)+dz*(zi-zs)+ai* (real) i_ac);
1111                     }
1112
1113                 } /* end of cycle "iat" */
1114             }     /* end of cycle "ix" */
1115         }         /* end of cycle "iy" */
1116     }             /* end of cycle "iz" */
1117
1118     sfree(wkatm);
1119     sfree(wkat1);
1120     sfree(wkdot);
1121     sfree(wkbox);
1122     sfree(wknb);
1123     if (box)
1124     {
1125         sfree(x);
1126     }
1127     if (mode & FLAG_VOLUME)
1128     {
1129         vol           = vol*FOURPI/(3.* (real) n_dot);
1130         *value_of_vol = vol;
1131     }
1132     if (mode & FLAG_DOTS)
1133     {
1134         *nu_dots = lfnr;
1135         *lidots  = dots;
1136     }
1137     if (mode & FLAG_ATOM_AREA)
1138     {
1139         *at_area = atom_area;
1140     }
1141     *value_of_area = area;
1142
1143     if (debug)
1144     {
1145         fprintf(debug, "area=%8.3f\n", area);
1146     }
1147
1148     return 0;
1149 }
1150
1151
1152 #if TEST_NSC > 0
1153 #define NAT 2
1154 main () {
1155
1156     int    i, j, ndots;
1157     real   co[3*NAT], ra[NAT], area, volume, a, b, c;
1158     real * dots;
1159     real * at_area;
1160     FILE  *fp;
1161
1162
1163     a  = 1.; c = 0.1;
1164     fp = fopen("nsc.txt", "w+");
1165     for (i = 1; i <= NAT; i++)
1166     {
1167         j         = i-1;
1168         co[3*i-3] = j*1*c;
1169         co[3*i-2] = j*1*c;
1170         co[3*i-1] = j*1*c;
1171         /*
1172            co[3*i-3] = i*1.4;
1173            co[3*i-2] = 0.;
1174            co[3*i-1] = 0.;
1175          */
1176         /*
1177            co[3*i-2] = a*0.3;
1178            a = -a; b=0;
1179            if (i%3 == 0) b=0.5;
1180            co[3*i-1] = b;
1181            ra[i-1] = 2.0;
1182          */
1183         ra[i-1] = (1.+j*0.5)*c;
1184     }
1185     /*
1186        if (NSC(co, ra, NAT, 42, NULL, &area,
1187      */
1188     if (NSC(co, ra, NAT, 42, NULL, &area,
1189             NULL, NULL, NULL, NULL))
1190     {
1191         ERROR("error in NSC");
1192     }
1193     fprintf(fp, "\n");
1194     fprintf(fp, "area     : %8.3f\n", area);
1195     fprintf(fp, "\n");
1196     fprintf(fp, "\n");
1197     fprintf(fp, "next call\n");
1198     fprintf(fp, "\n");
1199     fprintf(fp, "\n");
1200
1201     if (NSC(co, ra, NAT, 42, FLAG_VOLUME | FLAG_ATOM_AREA | FLAG_DOTS, &area,
1202             &at_area, &volume,
1203             &dots, &ndots))
1204     {
1205         ERROR("error in NSC");
1206     }
1207
1208     fprintf(fp, "\n");
1209     fprintf(fp, "area     : %8.3f\n", area);
1210     printf("area     : %8.3f\n", area);
1211     fprintf(fp, "volume   : %8.3f\n", volume);
1212     printf("volume   : %8.3f\n", volume);
1213     fprintf(fp, "ndots    : %8d\n", ndots);
1214     printf("ndots    : %8d\n", ndots);
1215     fprintf(fp, "\n");
1216     for (i = 1; i <= NAT; i++)
1217     {
1218         fprintf(fp, "%4d ATOM %7.2f %7.2f %7.2f  ra=%4.1f  area=%8.3f\n",
1219                 i, co[3*i-3], co[3*i-2], co[3*i-1], ra[i-1], at_area[i-1]);
1220     }
1221     fprintf(fp, "\n");
1222     fprintf(fp, "DOTS : %8d\n", ndots);
1223     for (i = 1; i <= ndots; i++)
1224     {
1225         fprintf(fp, "%4d DOTS %8.2f %8.2f %8.2f\n",
1226                 i, dots[3*i-3], dots[3*i-2], dots[3*i-1]);
1227     }
1228 }
1229 #endif