src/tools/nsc.c

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