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