src/gromacs/gmxana/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, 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 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include <stdio.h>
  42 #include <string.h>
  43 #include <stdlib.h>
  44 #include <math.h>
  45 #include <stdarg.h>
  46 /* Modified DvdS */
  47 #include "pbc.h"
  48 #include "macros.h"
  49 #include "vec.h"
  50 #include "smalloc.h"
  51 #include "nsc.h"
  52
  53 #define TEST_NSC 0
  54
  55 #define TEST_ARC 0
  56 #define TEST_DOD 0
  57 #define TEST_CUBE 0
  58
  59 #define UNSP_ICO_DOD      9
  60 #define UNSP_ICO_ARC     10
  61
  62 real   *xpunsp = NULL;
  63 real    del_cube;
  64 int    *ico_wk = NULL, *ico_pt = NULL;
  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   ndot, ico_cube_cb, i, j, k, l, ijk, tn, tl, tl2;
 534     real *xus;
 535     int  *work;
 536     real  x, y, z;
 537
 538     if (xpunsp)
 539     {
 540         free(xpunsp);
 541     }
 542     if (ico_wk)
 543     {
 544         free(ico_wk);
 545     }
 546
 547     k = 1; if (mode < 0)
 548     {
 549         k = 0; mode = -mode;
 550     }
 551     if (mode == UNSP_ICO_ARC)
 552     {
 553         ndot = ico_dot_arc(densit);
 554     }
 555     else if (mode == UNSP_ICO_DOD)
 556     {
 557         ndot = ico_dot_dod(densit);
 558     }
 559     else
 560     {
 561         WARNING("make_unsp: mode %c%d not allowed", (k) ? '+' : '-', mode);
 562         return 1;
 563     }
 564
 565     last_n_dot = ndot; last_densit = densit; last_unsp = mode;
 566     *num_dot   = ndot; if (k)
 567     {
 568         return 0;
 569     }
 570
 571     /* in the following the dots of the unit sphere may be resorted */
 572     last_unsp = -last_unsp;
 573
 574     /* determine distribution of points in elementary cubes */
 575     if (cubus)
 576     {
 577         ico_cube = cubus;
 578     }
 579     else
 580     {
 581         last_cubus = 0;
 582         i          = 1;
 583         while (i*i*i*2 < ndot)
 584         {
 585             i++;
 586         }
 587         ico_cube = max(i-1, 0);
 588     }
 589     ico_cube_cb = ico_cube*ico_cube*ico_cube;
 590     del_cube    = 2./((real)ico_cube);
 591     snew(work, ndot);
 592     xus = xpunsp;
 593     for (l = 0; l < ndot; l++)
 594     {
 595         i = max((int) floor((1.+xus[3*l])/del_cube), 0);
 596         if (i >= ico_cube)
 597         {
 598             i = ico_cube-1;
 599         }
 600         j = max((int) floor((1.+xus[1+3*l])/del_cube), 0);
 601         if (j >= ico_cube)
 602         {
 603             j = ico_cube-1;
 604         }
 605         k = max((int) floor((1.+xus[2+3*l])/del_cube), 0);
 606         if (k >= ico_cube)
 607         {
 608             k = ico_cube-1;
 609         }
 610         ijk     = i+j*ico_cube+k*ico_cube*ico_cube;
 611         work[l] = ijk;
 612     }
 613
 614     snew(ico_wk, 2*ico_cube_cb+1);
 615
 616     ico_pt = ico_wk+ico_cube_cb;
 617     for (l = 0; l < ndot; l++)
 618     {
 619         ico_wk[work[l]]++; /* dots per elementary cube */
 620     }
 621
 622     /* reordering of the coordinate array in accordance with box number */
 623     tn = 0;
 624     for (i = 0; i < ico_cube; i++)
 625     {
 626         for (j = 0; j < ico_cube; j++)
 627         {
 628             for (k = 0; k < ico_cube; k++)
 629             {
 630                 tl            = 0;
 631                 tl2           = tn;
 632                 ijk           = i+ico_cube*j+ico_cube*ico_cube*k;
 633                 *(ico_pt+ijk) = tn;
 634                 for (l = tl2; l < ndot; l++)
 635                 {
 636                     if (ijk == work[l])
 637                     {
 638                         x          = xus[3*l]; y = xus[1+3*l]; z = xus[2+3*l];
 639                         xus[3*l]   = xus[3*tn];
 640                         xus[1+3*l] = xus[1+3*tn]; xus[2+3*l] = xus[2+3*tn];
 641                         xus[3*tn]  = x; xus[1+3*tn] = y; xus[2+3*tn] = z;
 642                         ijk        = work[l]; work[l] = work[tn]; work[tn] = ijk;
 643                         tn++; tl++;
 644                     }
 645                 }
 646                 *(ico_wk+ijk) = tl;
 647             } /* cycle k */
 648         }     /* cycle j */
 649     }         /* cycle i */
 650     free(work);
 651
 652     return 0;
 653 }
 654
 655
 656 typedef struct _stwknb {
 657     real x;
 658     real y;
 659     real z;
 660     real dot;
 661 } Neighb;
 662
 663 int nsc_dclm_pbc(rvec *coords, real *radius, int nat,
 664                  int  densit, int mode,
 665                  real *value_of_area, real **at_area,
 666                  real *value_of_vol,
 667                  real **lidots, int *nu_dots,
 668                  atom_id index[], int ePBC, matrix box)
 669 {
 670     int      iat, i, ii, iii, ix, iy, iz, ixe, ixs, iye, iys, ize, izs, i_ac;
 671     int      jat, j, jj, jjj, jx, jy, jz;
 672     int      distribution;
 673     int      l;
 674     int      maxnei, nnei, last, maxdots = 0;
 675     int     *wkdot = NULL, *wkbox = NULL, *wkat1 = NULL, *wkatm = NULL;
 676     Neighb  *wknb, *ctnb;
 677     int      iii1, iii2, iiat, lfnr = 0, i_at, j_at;
 678     real     dx, dy, dz, dd, ai, aisq, ajsq, aj, as, a;
 679     real     xi, yi, zi, xs = 0., ys = 0., zs = 0.;
 680     real     dotarea, area, vol = 0.;
 681     real    *xus, *dots = NULL, *atom_area = NULL;
 682     int      nxbox, nybox, nzbox, nxy, nxyz;
 683     real     xmin = 0, ymin = 0, zmin = 0, xmax, ymax, zmax, ra2max, d, *pco;
 684     /* Added DvdS 2006-07-19 */
 685     t_pbc    pbc;
 686     rvec     ddx, *x = NULL;
 687     int      iat_xx, jat_xx;
 688
 689     distribution = unsp_type(densit);
 690     if (distribution != -last_unsp || last_cubus != 4 ||
 691         (densit != last_densit && densit != last_n_dot))
 692     {
 693         if (make_unsp(densit, (-distribution), &n_dot, 4))
 694         {
 695             return 1;
 696         }
 697     }
 698     xus = xpunsp;
 699
 700     dotarea = FOURPI/(real) n_dot;
 701     area    = 0.;
 702
 703     if (debug)
 704     {
 705         fprintf(debug, "nsc_dclm: n_dot=%5d %9.3f\n", n_dot, dotarea);
 706     }
 707
 708     /* start with neighbour list */
 709     /* calculate neighbour list with the box algorithm */
 710     if (nat == 0)
 711     {
 712         WARNING("nsc_dclm: no surface atoms selected");
 713         return 1;
 714     }
 715     if (mode & FLAG_VOLUME)
 716     {
 717         vol = 0.;
 718     }
 719     if (mode & FLAG_DOTS)
 720     {
 721         maxdots = (3*n_dot*nat)/10;
 722         /* should be set to NULL on first user call */
 723         if (dots == NULL)
 724         {
 725             snew(dots, maxdots);
 726         }
 727         else
 728         {
 729             srenew(dots, maxdots);
 730         }
 731
 732         lfnr = 0;
 733     }
 734     if (mode & FLAG_ATOM_AREA)
 735     {
 736         /* should be set to NULL on first user call */
 737         if (atom_area == NULL)
 738         {
 739             snew(atom_area, nat);
 740         }
 741         else
 742         {
 743             srenew(atom_area, nat);
 744         }
 745     }
 746
 747     /* Compute minimum size for grid cells */
 748     ra2max = radius[index[0]];
 749     for (iat_xx = 1; (iat_xx < nat); iat_xx++)
 750     {
 751         iat    = index[iat_xx];
 752         ra2max = max(ra2max, radius[iat]);
 753     }
 754     ra2max = 2*ra2max;
 755
 756     /* Added DvdS 2006-07-19 */
 757     /* Updated 2008-10-09 */
 758     if (box)
 759     {
 760         set_pbc(&pbc, ePBC, box);
 761         snew(x, nat);
 762         for (i = 0; (i < nat); i++)
 763         {
 764             iat  = index[0];
 765             copy_rvec(coords[iat], x[i]);
 766         }
 767         put_atoms_in_triclinic_unitcell(ecenterTRIC, box, nat, x);
 768         nxbox = max(1, floor(norm(box[XX])/ra2max));
 769         nybox = max(1, floor(norm(box[YY])/ra2max));
 770         nzbox = max(1, floor(norm(box[ZZ])/ra2max));
 771         if (debug)
 772         {
 773             fprintf(debug, "nbox = %d, %d, %d\n", nxbox, nybox, nzbox);
 774         }
 775     }
 776     else
 777     {
 778         /* dimensions of atomic set, cell edge is 2*ra_max */
 779         iat    = index[0];
 780         xmin   = coords[iat][XX]; xmax = xmin; xs = xmin;
 781         ymin   = coords[iat][YY]; ymax = ymin; ys = ymin;
 782         zmin   = coords[iat][ZZ]; zmax = zmin; zs = zmin;
 783         ra2max = radius[iat];
 784
 785         for (iat_xx = 1; (iat_xx < nat); iat_xx++)
 786         {
 787             iat  = index[iat_xx];
 788             pco  = coords[iat];
 789             xmin = min(xmin, *pco);     xmax = max(xmax, *pco);
 790             ymin = min(ymin, *(pco+1)); ymax = max(ymax, *(pco+1));
 791             zmin = min(zmin, *(pco+2)); zmax = max(zmax, *(pco+2));
 792             xs   = xs+ *pco; ys = ys+ *(pco+1); zs = zs+ *(pco+2);
 793         }
 794         xs = xs/ (real) nat;
 795         ys = ys/ (real) nat;
 796         zs = zs/ (real) nat;
 797         if (debug)
 798         {
 799             fprintf(debug, "nsc_dclm: n_dot=%5d ra2max=%9.3f %9.3f\n", n_dot, ra2max, dotarea);
 800         }
 801
 802         d    = xmax-xmin; nxbox = (int) max(ceil(d/ra2max), 1.);
 803         d    = (((real)nxbox)*ra2max-d)/2.;
 804         xmin = xmin-d; xmax = xmax+d;
 805         d    = ymax-ymin; nybox = (int) max(ceil(d/ra2max), 1.);
 806         d    = (((real)nybox)*ra2max-d)/2.;
 807         ymin = ymin-d; ymax = ymax+d;
 808         d    = zmax-zmin; nzbox = (int) max(ceil(d/ra2max), 1.);
 809         d    = (((real)nzbox)*ra2max-d)/2.;
 810         zmin = zmin-d; zmax = zmax+d;
 811     }
 812     /* Help variables */
 813     nxy  = nxbox*nybox;
 814     nxyz = nxy*nzbox;
 815
 816     /* box number of atoms */
 817     snew(wkatm, nat);
 818     snew(wkat1, nat);
 819     snew(wkdot, n_dot);
 820     snew(wkbox, nxyz+1);
 821
 822     if (box)
 823     {
 824         matrix box_1;
 825         rvec   x_1;
 826         int    ix, iy, iz, m;
 827         m_inv(box, box_1);
 828         for (i = 0; (i < nat); i++)
 829         {
 830             mvmul(box_1, x[i], x_1);
 831             ix = ((int)floor(x_1[XX]*nxbox) + 2*nxbox) % nxbox;
 832             iy = ((int)floor(x_1[YY]*nybox) + 2*nybox) % nybox;
 833             iz = ((int)floor(x_1[ZZ]*nzbox) + 2*nzbox) % nzbox;
 834             j  =  ix + iy*nxbox + iz*nxbox*nybox;
 835             if (debug)
 836             {
 837                 fprintf(debug, "Atom %d cell index %d. x = (%8.3f,%8.3f,%8.3f) fc = (%8.3f,%8.3f,%8.3f)\n",
 838                         i, j, x[i][XX], x[i][YY], x[i][ZZ], x_1[XX], x_1[YY], x_1[ZZ]);
 839             }
 840             range_check(j, 0, nxyz);
 841             wkat1[i] = j;
 842             wkbox[j]++;
 843         }
 844     }
 845     else
 846     {
 847         /* Put the atoms in their boxes */
 848         for (iat_xx = 0; (iat_xx < nat); iat_xx++)
 849         {
 850             iat           = index[iat_xx];
 851             pco           = coords[iat];
 852             i             = (int) max(floor((pco[XX]-xmin)/ra2max), 0);
 853             i             = min(i, nxbox-1);
 854             j             = (int) max(floor((pco[YY]-ymin)/ra2max), 0);
 855             j             = min(j, nybox-1);
 856             l             = (int) max(floor((pco[ZZ]-zmin)/ra2max), 0);
 857             l             = min(l, nzbox-1);
 858             i             = i+j*nxbox+l*nxy;
 859             wkat1[iat_xx] = i;
 860             wkbox[i]++;
 861         }
 862     }
 863
 864     /* sorting of atoms in accordance with box numbers */
 865     j = wkbox[0];
 866     for (i = 1; i < nxyz; i++)
 867     {
 868         j = max(wkbox[i], j);
 869     }
 870     for (i = 1; i <= nxyz; i++)
 871     {
 872         wkbox[i] += wkbox[i-1];
 873     }
 874
 875     /* maxnei = (int) floor(ra2max*ra2max*ra2max*0.5); */
 876     maxnei = min(nat, 27*j);
 877     snew(wknb, maxnei);
 878     for (iat_xx = 0; iat_xx < nat; iat_xx++)
 879     {
 880         iat = index[iat_xx];
 881         range_check(wkat1[iat_xx], 0, nxyz);
 882         wkatm[--wkbox[wkat1[iat_xx]]] = iat_xx;
 883         if (debug)
 884         {
 885             fprintf(debug, "atom %5d on place %5d\n", iat, wkbox[wkat1[iat_xx]]);
 886         }
 887     }
 888
 889     if (debug)
 890     {
 891         fprintf(debug, "nsc_dclm: n_dot=%5d ra2max=%9.3f %9.3f\n",
 892                 n_dot, ra2max, dotarea);
 893         fprintf(debug, "neighbour list calculated/box(xyz):%d %d %d\n",
 894                 nxbox, nybox, nzbox);
 895
 896         for (i = 0; i < nxyz; i++)
 897         {
 898             fprintf(debug, "box %6d : atoms %4d-%4d    %5d\n",
 899                     i, wkbox[i], wkbox[i+1]-1, wkbox[i+1]-wkbox[i]);
 900         }
 901         for (i = 0; i < nat; i++)
 902         {
 903             fprintf(debug, "list place %5d by atom %7d\n", i, index[wkatm[i]]);
 904         }
 905     }
 906
 907     /* calculate surface for all atoms, step cube-wise */
 908     for (iz = 0; iz < nzbox; iz++)
 909     {
 910         iii = iz*nxy;
 911         if (box)
 912         {
 913             izs = iz-1;
 914             ize = min(iz+2, izs+nzbox);
 915         }
 916         else
 917         {
 918             izs = max(iz-1, 0);
 919             ize = min(iz+2, nzbox);
 920         }
 921         for (iy = 0; iy < nybox; iy++)
 922         {
 923             ii = iy*nxbox+iii;
 924             if (box)
 925             {
 926                 iys = iy-1;
 927                 iye = min(iy+2, iys+nybox);
 928             }
 929             else
 930             {
 931                 iys = max(iy-1, 0);
 932                 iye = min(iy+2, nybox);
 933             }
 934             for (ix = 0; ix < nxbox; ix++)
 935             {
 936                 i    = ii+ix;
 937                 iii1 = wkbox[i];
 938                 iii2 = wkbox[i+1];
 939                 if (iii1 >= iii2)
 940                 {
 941                     continue;
 942                 }
 943                 if (box)
 944                 {
 945                     ixs = ix-1;
 946                     ixe = min(ix+2, ixs+nxbox);
 947                 }
 948                 else
 949                 {
 950                     ixs = max(ix-1, 0);
 951                     ixe = min(ix+2, nxbox);
 952                 }
 953                 iiat = 0;
 954                 /* make intermediate atom list */
 955                 for (jz = izs; jz < ize; jz++)
 956                 {
 957                     jjj = ((jz+nzbox) % nzbox)*nxy;
 958                     for (jy = iys; jy < iye; jy++)
 959                     {
 960                         jj = ((jy+nybox) % nybox)*nxbox+jjj;
 961                         for (jx = ixs; jx < ixe; jx++)
 962                         {
 963                             j = jj+((jx+nxbox) % nxbox);
 964                             for (jat = wkbox[j]; jat < wkbox[j+1]; jat++)
 965                             {
 966                                 range_check(wkatm[jat], 0, nat);
 967                                 range_check(iiat, 0, nat);
 968                                 wkat1[iiat] = wkatm[jat];
 969                                 iiat++;
 970                             } /* end of cycle "jat" */
 971                         }     /* end of cycle "jx" */
 972                     }         /* end of cycle "jy" */
 973                 }             /* end of cycle "jz" */
 974                 for (iat = iii1; iat < iii2; iat++)
 975                 {
 976                     i_at = index[wkatm[iat]];
 977                     ai   = radius[i_at];
 978                     aisq = ai*ai;
 979                     pco  = coords[i_at];
 980                     xi   = pco[XX]; yi = pco[YY]; zi = pco[ZZ];
 981                     for (i = 0; i < n_dot; i++)
 982                     {
 983                         wkdot[i] = 0;
 984                     }
 985
 986                     ctnb = wknb; nnei = 0;
 987                     for (j = 0; j < iiat; j++)
 988                     {
 989                         j_at = index[wkat1[j]];
 990                         if (j_at == i_at)
 991                         {
 992                             continue;
 993                         }
 994
 995                         aj   = radius[j_at];
 996                         ajsq = aj*aj;
 997                         pco  = coords[j_at];
 998
 999                         /* Added DvdS 2006-07-19 */
1000                         if (box)
1001                         {
1002                             /*rvec xxi;
1003
1004                                xxi[XX] = xi;
1005                                xxi[YY] = yi;
1006                                xxi[ZZ] = zi;
1007                                pbc_dx(&pbc,pco,xxi,ddx);*/
1008                             pbc_dx(&pbc, coords[j_at], coords[i_at], ddx);
1009                             dx = ddx[XX];
1010                             dy = ddx[YY];
1011                             dz = ddx[ZZ];
1012                         }
1013                         else
1014                         {
1015                             dx = pco[XX]-xi;
1016                             dy = pco[YY]-yi;
1017                             dz = pco[ZZ]-zi;
1018                         }
1019                         dd = dx*dx+dy*dy+dz*dz;
1020                         as = ai+aj;
1021                         if (dd > as*as)
1022                         {
1023                             continue;
1024                         }
1025                         nnei++;
1026                         ctnb->x   = dx;
1027                         ctnb->y   = dy;
1028                         ctnb->z   = dz;
1029                         ctnb->dot = (dd+aisq-ajsq)/(2.*ai); /* reference dot product */
1030                         ctnb++;
1031                     }
1032
1033                     /* check points on accessibility */
1034                     if (nnei)
1035                     {
1036                         last = 0; i_ac = 0;
1037                         for (l = 0; l < n_dot; l++)
1038                         {
1039                             if (xus[3*l]*(wknb+last)->x+
1040                                 xus[1+3*l]*(wknb+last)->y+
1041                                 xus[2+3*l]*(wknb+last)->z <= (wknb+last)->dot)
1042                             {
1043                                 for (j = 0; j < nnei; j++)
1044                                 {
1045                                     if (xus[3*l]*(wknb+j)->x+xus[1+3*l]*(wknb+j)->y+
1046                                         xus[2+3*l]*(wknb+j)->z > (wknb+j)->dot)
1047                                     {
1048                                         last = j;
1049                                         break;
1050                                     }
1051                                 }
1052                                 if (j >= nnei)
1053                                 {
1054                                     i_ac++;
1055                                     wkdot[l] = 1;
1056                                 }
1057                             } /* end of cycle j */
1058                         }     /* end of cycle l */
1059                     }
1060                     else
1061                     {
1062                         i_ac  = n_dot;
1063                         for (l = 0; l < n_dot; l++)
1064                         {
1065                             wkdot[l] = 1;
1066                         }
1067                     }
1068
1069                     if (debug)
1070                     {
1071                         fprintf(debug, "i_ac=%d, dotarea=%8.3f, aisq=%8.3f\n",
1072                                 i_ac, dotarea, aisq);
1073                     }
1074
1075                     a    = aisq*dotarea* (real) i_ac;
1076                     area = area + a;
1077                     if (mode & FLAG_ATOM_AREA)
1078                     {
1079                         range_check(wkatm[iat], 0, nat);
1080                         atom_area[wkatm[iat]] = a;
1081                     }
1082                     if (mode & FLAG_DOTS)
1083                     {
1084                         for (l = 0; l < n_dot; l++)
1085                         {
1086                             if (wkdot[l])
1087                             {
1088                                 lfnr++;
1089                                 if (maxdots <= 3*lfnr+1)
1090                                 {
1091                                     maxdots = maxdots+n_dot*3;
1092                                     srenew(dots, maxdots);
1093                                 }
1094                                 dots[3*lfnr-3] = ai*xus[3*l]+xi;
1095                                 dots[3*lfnr-2] = ai*xus[1+3*l]+yi;
1096                                 dots[3*lfnr-1] = ai*xus[2+3*l]+zi;
1097                             }
1098                         }
1099                     }
1100                     if (mode & FLAG_VOLUME)
1101                     {
1102                         dx = 0.; dy = 0.; dz = 0.;
1103                         for (l = 0; l < n_dot; l++)
1104                         {
1105                             if (wkdot[l])
1106                             {
1107                                 dx = dx+xus[3*l];
1108                                 dy = dy+xus[1+3*l];
1109                                 dz = dz+xus[2+3*l];
1110                             }
1111                         }
1112                         vol = vol+aisq*(dx*(xi-xs)+dy*(yi-ys)+dz*(zi-zs)+ai* (real) i_ac);
1113                     }
1114
1115                 } /* end of cycle "iat" */
1116             }     /* end of cycle "ix" */
1117         }         /* end of cycle "iy" */
1118     }             /* end of cycle "iz" */
1119
1120     sfree(wkatm);
1121     sfree(wkat1);
1122     sfree(wkdot);
1123     sfree(wkbox);
1124     sfree(wknb);
1125     if (box)
1126     {
1127         sfree(x);
1128     }
1129     if (mode & FLAG_VOLUME)
1130     {
1131         vol           = vol*FOURPI/(3.* (real) n_dot);
1132         *value_of_vol = vol;
1133     }
1134     if (mode & FLAG_DOTS)
1135     {
1136         *nu_dots = lfnr;
1137         *lidots  = dots;
1138     }
1139     if (mode & FLAG_ATOM_AREA)
1140     {
1141         *at_area = atom_area;
1142     }
1143     *value_of_area = area;
1144
1145     if (debug)
1146     {
1147         fprintf(debug, "area=%8.3f\n", area);
1148     }
1149
1150     return 0;
1151 }
1152
1153
1154 #if TEST_NSC > 0
1155 #define NAT 2
1156 main () {
1157
1158     int    i, j, ndots;
1159     real   co[3*NAT], ra[NAT], area, volume, a, b, c;
1160     real * dots;
1161     real * at_area;
1162     FILE  *fp;
1163
1164
1165     a  = 1.; c = 0.1;
1166     fp = fopen("nsc.txt", "w+");
1167     for (i = 1; i <= NAT; i++)
1168     {
1169         j         = i-1;
1170         co[3*i-3] = j*1*c;
1171         co[3*i-2] = j*1*c;
1172         co[3*i-1] = j*1*c;
1173         /*
1174            co[3*i-3] = i*1.4;
1175            co[3*i-2] = 0.;
1176            co[3*i-1] = 0.;
1177          */
1178         /*
1179            co[3*i-2] = a*0.3;
1180            a = -a; b=0;
1181            if (i%3 == 0) b=0.5;
1182            co[3*i-1] = b;
1183            ra[i-1] = 2.0;
1184          */
1185         ra[i-1] = (1.+j*0.5)*c;
1186     }
1187     /*
1188        if (NSC(co, ra, NAT, 42, NULL, &area,
1189      */
1190     if (NSC(co, ra, NAT, 42, NULL, &area,
1191             NULL, NULL, NULL, NULL))
1192     {
1193         ERROR("error in NSC");
1194     }
1195     fprintf(fp, "\n");
1196     fprintf(fp, "area     : %8.3f\n", area);
1197     fprintf(fp, "\n");
1198     fprintf(fp, "\n");
1199     fprintf(fp, "next call\n");
1200     fprintf(fp, "\n");
1201     fprintf(fp, "\n");
1202
1203     if (NSC(co, ra, NAT, 42, FLAG_VOLUME | FLAG_ATOM_AREA | FLAG_DOTS, &area,
1204             &at_area, &volume,
1205             &dots, &ndots))
1206     {
1207         ERROR("error in NSC");
1208     }
1209
1210     fprintf(fp, "\n");
1211     fprintf(fp, "area     : %8.3f\n", area);
1212     printf("area     : %8.3f\n", area);
1213     fprintf(fp, "volume   : %8.3f\n", volume);
1214     printf("volume   : %8.3f\n", volume);
1215     fprintf(fp, "ndots    : %8d\n", ndots);
1216     printf("ndots    : %8d\n", ndots);
1217     fprintf(fp, "\n");
1218     for (i = 1; i <= NAT; i++)
1219     {
1220         fprintf(fp, "%4d ATOM %7.2f %7.2f %7.2f  ra=%4.1f  area=%8.3f\n",
1221                 i, co[3*i-3], co[3*i-2], co[3*i-1], ra[i-1], at_area[i-1]);
1222     }
1223     fprintf(fp, "\n");
1224     fprintf(fp, "DOTS : %8d\n", ndots);
1225     for (i = 1; i <= ndots; i++)
1226     {
1227         fprintf(fp, "%4d DOTS %8.2f %8.2f %8.2f\n",
1228                 i, dots[3*i-3], dots[3*i-2], dots[3*i-1]);
1229     }
1230 }
1231 #endif