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