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