3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Groningen Machine for Chemical Simulation
57 #define UNSP_ICO_DOD 9
58 #define UNSP_ICO_ARC 10
62 int *ico_wk=NULL, *ico_pt=NULL;
63 int n_dot, ico_cube, last_n_dot=0, last_densit=0, last_unsp=0;
66 #define FOURPI (4.*M_PI)
67 #define TORAD(A) ((A)*0.017453293)
70 #define UPDATE_FL __file__=__FILE__,__line__=__LINE__
71 const char * __file__; /* declared versions of macros */
72 int __line__; /* __FILE__ and __LINE__ */
74 #define ERROR UPDATE_FL,error
75 void error(const char *fmt, ...) {
78 "\n---> ERROR when executing line %i in file %s !\n",
81 vfprintf(stderr, fmt, args);
83 fprintf(stderr, "\n---> Execution stopped !\n\n");
86 #define WARNING UPDATE_FL,warning2
87 void warning2(const char *fmt, ...) {
90 "\n---> WARNING : line %i in file %s\n",
93 vfprintf(stderr, fmt, args);
95 fprintf(stderr, " ...!\n\n");
100 #define ASIN safe_asin
101 real safe_asin(real f) {
102 if ( (fabs(f) < 1.00) ) return( asin(f) );
103 if ( (fabs(f) - 1.00) <= DP_TOL )
104 ERROR("ASIN : invalid argument %f", f);
108 #define CALLOC(n, size) mycalloc(__FILE__,__LINE__, n, size)
109 void * mycalloc(const char * filename, const int linenr,
110 size_t nelem, size_t elsize) {
112 ip = (int *) calloc(nelem, elsize);
114 ERROR("CALLOC : failed in file %s at line %d", filename, linenr);
118 #define REALLOC(ptr, size) myrealloc(__FILE__,__LINE__, ptr, size)
119 void * myrealloc(const char * filename, const int linenr,
120 void * ptr, size_t size) {
122 ip = (int *) realloc(ptr, size);
124 ERROR("REALLOC : failed in file %s at line %d", filename, linenr);
129 /* routines for dot distributions on the surface of the unit sphere */
132 void icosaeder_vertices(real *xus) {
133 rh = sqrt(1.-2.*cos(TORAD(72.)))/(1.-cos(TORAD(72.)));
134 rg = cos(TORAD(72.))/(1.-cos(TORAD(72.)));
135 /* icosaeder vertices */
136 xus[ 0] = 0.; xus[ 1] = 0.; xus[ 2] = 1.;
137 xus[ 3] = rh*cos(TORAD(72.)); xus[ 4] = rh*sin(TORAD(72.)); xus[ 5] = rg;
138 xus[ 6] = rh*cos(TORAD(144.)); xus[ 7] = rh*sin(TORAD(144.)); xus[ 8] = rg;
139 xus[ 9] = rh*cos(TORAD(216.)); xus[10] = rh*sin(TORAD(216.)); xus[11] = rg;
140 xus[12] = rh*cos(TORAD(288.)); xus[13] = rh*sin(TORAD(288.)); xus[14] = rg;
141 xus[15] = rh; xus[16] = 0; xus[17] = rg;
142 xus[18] = rh*cos(TORAD(36.)); xus[19] = rh*sin(TORAD(36.)); xus[20] = -rg;
143 xus[21] = rh*cos(TORAD(108.)); xus[22] = rh*sin(TORAD(108.)); xus[23] = -rg;
144 xus[24] = -rh; xus[25] = 0; xus[26] = -rg;
145 xus[27] = rh*cos(TORAD(252.)); xus[28] = rh*sin(TORAD(252.)); xus[29] = -rg;
146 xus[30] = rh*cos(TORAD(324.)); xus[31] = rh*sin(TORAD(324.)); xus[32] = -rg;
147 xus[33] = 0.; xus[34] = 0.; xus[35] = -1.;
151 void divarc(real x1, real y1, real z1,
152 real x2, real y2, real z2,
153 int div1, int div2, real *xr, real *yr, real *zr) {
155 real xd, yd, zd, dd, d1, d2, s, x, y, z;
156 real phi, sphi, cphi;
161 dd = sqrt(xd*xd+yd*yd+zd*zd);
162 if (dd < DP_TOL) ERROR("divarc: rotation axis of length %f", dd);
164 d1 = x1*x1+y1*y1+z1*z1;
165 if (d1 < 0.5) ERROR("divarc: vector 1 of sq.length %f", d1);
166 d2 = x2*x2+y2*y2+z2*z2;
167 if (d2 < 0.5) ERROR("divarc: vector 2 of sq.length %f", d2);
169 phi = ASIN(dd/sqrt(d1*d2));
170 phi = phi*((real)div1)/((real)div2);
171 sphi = sin(phi); cphi = cos(phi);
172 s = (x1*xd+y1*yd+z1*zd)/dd;
174 x = xd*s*(1.-cphi)/dd + x1 * cphi + (yd*z1-y1*zd)*sphi/dd;
175 y = yd*s*(1.-cphi)/dd + y1 * cphi + (zd*x1-z1*xd)*sphi/dd;
176 z = zd*s*(1.-cphi)/dd + z1 * cphi + (xd*y1-x1*yd)*sphi/dd;
177 dd = sqrt(x*x+y*y+z*z);
178 *xr = x/dd; *yr = y/dd; *zr = z/dd;
181 int ico_dot_arc(int densit) { /* densit...required dots per unit sphere */
182 /* dot distribution on a unit sphere based on an icosaeder *
183 * great circle average refining of icosahedral face */
185 int i, j, k, tl, tl2, tn, tess;
186 real a, d, x, y, z, x2, y2, z2, x3, y3, z3;
187 real xij, yij, zij, xji, yji, zji, xik, yik, zik, xki, yki, zki,
188 xjk, yjk, zjk, xkj, ykj, zkj;
191 /* calculate tessalation level */
192 a = sqrt((((real) densit)-2.)/10.);
193 tess = (int) ceil(a);
194 n_dot = 10*tess*tess+2;
195 if (n_dot < densit) {
196 ERROR("ico_dot_arc: error in formula for tessalation level (%d->%d, %d)",
197 tess, n_dot, densit);
200 xus = (real *) CALLOC(3*n_dot, sizeof(real));
202 icosaeder_vertices(xus);
206 a = rh*rh*2.*(1.-cos(TORAD(72.)));
207 /* calculate tessalation of icosaeder edges */
208 for (i=0; i<11; i++) {
209 for (j=i+1; j<12; j++) {
210 x = xus[3*i]-xus[3*j];
211 y = xus[1+3*i]-xus[1+3*j]; z = xus[2+3*i]-xus[2+3*j];
213 if (fabs(a-d) > DP_TOL) continue;
214 for (tl=1; tl<tess; tl++) {
215 if (tn >= n_dot) { ERROR("ico_dot: tn exceeds dimension of xus"); }
216 divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
217 xus[3*j], xus[1+3*j], xus[2+3*j],
218 tl, tess, &xus[3*tn], &xus[1+3*tn], &xus[2+3*tn]);
223 /* calculate tessalation of icosaeder faces */
224 for (i=0; i<10; i++) {
225 for (j=i+1; j<11; j++) {
226 x = xus[3*i]-xus[3*j];
227 y = xus[1+3*i]-xus[1+3*j]; z = xus[2+3*i]-xus[2+3*j];
229 if (fabs(a-d) > DP_TOL) continue;
231 for (k=j+1; k<12; k++) {
232 x = xus[3*i]-xus[3*k];
233 y = xus[1+3*i]-xus[1+3*k]; z = xus[2+3*i]-xus[2+3*k];
235 if (fabs(a-d) > DP_TOL) continue;
236 x = xus[3*j]-xus[3*k];
237 y = xus[1+3*j]-xus[1+3*k]; z = xus[2+3*j]-xus[2+3*k];
239 if (fabs(a-d) > DP_TOL) continue;
240 for (tl=1; tl<tess-1; tl++) {
241 divarc(xus[3*j], xus[1+3*j], xus[2+3*j],
242 xus[3*i], xus[1+3*i], xus[2+3*i],
243 tl, tess, &xji, &yji, &zji);
244 divarc(xus[3*k], xus[1+3*k], xus[2+3*k],
245 xus[3*i], xus[1+3*i], xus[2+3*i],
246 tl, tess, &xki, &yki, &zki);
248 for (tl2=1; tl2<tess-tl; tl2++) {
249 divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
250 xus[3*j], xus[1+3*j], xus[2+3*j],
251 tl2, tess, &xij, &yij, &zij);
252 divarc(xus[3*k], xus[1+3*k], xus[2+3*k],
253 xus[3*j], xus[1+3*j], xus[2+3*j],
254 tl2, tess, &xkj, &ykj, &zkj);
255 divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
256 xus[3*k], xus[1+3*k], xus[2+3*k],
257 tess-tl-tl2, tess, &xik, &yik, &zik);
258 divarc(xus[3*j], xus[1+3*j], xus[2+3*j],
259 xus[3*k], xus[1+3*k], xus[2+3*k],
260 tess-tl-tl2, tess, &xjk, &yjk, &zjk);
261 if (tn >= n_dot) ERROR("ico_dot: tn exceeds dimension of xus");
262 divarc(xki, yki, zki, xji, yji, zji, tl2, tess-tl,
264 divarc(xkj, ykj, zkj, xij, yij, zij, tl, tess-tl2,
266 divarc(xjk, yjk, zjk, xik, yik, zik, tl, tl+tl2,
268 x = x+x2+x3; y = y+y2+y3; z = z+z2+z3;
269 d = sqrt(x*x+y*y+z*z);
280 ERROR("ico_dot: n_dot(%d) and tn(%d) differ", n_dot, tn);
282 } /* end of if (tess > 1) */
284 } /* end of routine ico_dot_arc */
286 int ico_dot_dod(int densit) { /* densit...required dots per unit sphere */
287 /* dot distribution on a unit sphere based on an icosaeder *
288 * great circle average refining of icosahedral face */
290 int i, j, k, tl, tl2, tn, tess, j1, j2;
291 real a, d, x, y, z, x2, y2, z2, x3, y3, z3, ai_d, adod;
292 real xij, yij, zij, xji, yji, zji, xik, yik, zik, xki, yki, zki,
293 xjk, yjk, zjk, xkj, ykj, zkj;
295 /* calculate tesselation level */
296 a = sqrt((((real) densit)-2.)/30.);
297 tess = max((int) ceil(a), 1);
298 n_dot = 30*tess*tess+2;
299 if (n_dot < densit) {
300 ERROR("ico_dot_dod: error in formula for tessalation level (%d->%d, %d)",
301 tess, n_dot, densit);
304 xus = (real *) CALLOC(3*n_dot, sizeof(real));
306 icosaeder_vertices(xus);
309 /* square of the edge of an icosaeder */
310 a = rh*rh*2.*(1.-cos(TORAD(72.)));
311 /* dodecaeder vertices */
312 for (i=0; i<10; i++) {
313 for (j=i+1; j<11; j++) {
314 x = xus[3*i]-xus[3*j];
315 y = xus[1+3*i]-xus[1+3*j]; z = xus[2+3*i]-xus[2+3*j];
317 if (fabs(a-d) > DP_TOL) continue;
318 for (k=j+1; k<12; k++) {
319 x = xus[3*i]-xus[3*k];
320 y = xus[1+3*i]-xus[1+3*k]; z = xus[2+3*i]-xus[2+3*k];
322 if (fabs(a-d) > DP_TOL) continue;
323 x = xus[3*j]-xus[3*k];
324 y = xus[1+3*j]-xus[1+3*k]; z = xus[2+3*j]-xus[2+3*k];
326 if (fabs(a-d) > DP_TOL) continue;
327 x = xus[ 3*i]+xus[ 3*j]+xus[ 3*k];
328 y = xus[1+3*i]+xus[1+3*j]+xus[1+3*k];
329 z = xus[2+3*i]+xus[2+3*j]+xus[2+3*k];
330 d = sqrt(x*x+y*y+z*z);
331 xus[3*tn]=x/d; xus[1+3*tn]=y/d; xus[2+3*tn]=z/d;
339 /* square of the edge of an dodecaeder */
340 adod = 4.*(cos(TORAD(108.))-cos(TORAD(120.)))/(1.-cos(TORAD(120.)));
341 /* square of the distance of two adjacent vertices of ico- and dodecaeder */
342 ai_d = 2.*(1.-sqrt(1.-a/3.));
344 /* calculate tessalation of mixed edges */
345 for (i=0; i<31; i++) {
346 j1 = 12; j2 = 32; a = ai_d;
347 if (i>=12) { j1=i+1; a = adod; }
348 for (j=j1; j<j2; j++) {
349 x = xus[3*i]-xus[3*j];
350 y = xus[1+3*i]-xus[1+3*j]; z = xus[2+3*i]-xus[2+3*j];
352 if (fabs(a-d) > DP_TOL) continue;
353 for (tl=1; tl<tess; tl++) {
355 ERROR("ico_dot: tn exceeds dimension of xus");
357 divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
358 xus[3*j], xus[1+3*j], xus[2+3*j],
359 tl, tess, &xus[3*tn], &xus[1+3*tn], &xus[2+3*tn]);
364 /* calculate tessalation of pentakisdodecahedron faces */
365 for (i=0; i<12; i++) {
366 for (j=12; j<31; j++) {
367 x = xus[3*i]-xus[3*j];
368 y = xus[1+3*i]-xus[1+3*j]; z = xus[2+3*i]-xus[2+3*j];
370 if (fabs(ai_d-d) > DP_TOL) continue;
372 for (k=j+1; k<32; k++) {
373 x = xus[3*i]-xus[3*k];
374 y = xus[1+3*i]-xus[1+3*k]; z = xus[2+3*i]-xus[2+3*k];
376 if (fabs(ai_d-d) > DP_TOL) continue;
377 x = xus[3*j]-xus[3*k];
378 y = xus[1+3*j]-xus[1+3*k]; z = xus[2+3*j]-xus[2+3*k];
380 if (fabs(adod-d) > DP_TOL) continue;
381 for (tl=1; tl<tess-1; tl++) {
382 divarc(xus[3*j], xus[1+3*j], xus[2+3*j],
383 xus[3*i], xus[1+3*i], xus[2+3*i],
384 tl, tess, &xji, &yji, &zji);
385 divarc(xus[3*k], xus[1+3*k], xus[2+3*k],
386 xus[3*i], xus[1+3*i], xus[2+3*i],
387 tl, tess, &xki, &yki, &zki);
389 for (tl2=1; tl2<tess-tl; tl2++) {
390 divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
391 xus[3*j], xus[1+3*j], xus[2+3*j],
392 tl2, tess, &xij, &yij, &zij);
393 divarc(xus[3*k], xus[1+3*k], xus[2+3*k],
394 xus[3*j], xus[1+3*j], xus[2+3*j],
395 tl2, tess, &xkj, &ykj, &zkj);
396 divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
397 xus[3*k], xus[1+3*k], xus[2+3*k],
398 tess-tl-tl2, tess, &xik, &yik, &zik);
399 divarc(xus[3*j], xus[1+3*j], xus[2+3*j],
400 xus[3*k], xus[1+3*k], xus[2+3*k],
401 tess-tl-tl2, tess, &xjk, &yjk, &zjk);
403 ERROR("ico_dot: tn exceeds dimension of xus");
405 divarc(xki, yki, zki, xji, yji, zji, tl2, tess-tl,
407 divarc(xkj, ykj, zkj, xij, yij, zij, tl, tess-tl2,
409 divarc(xjk, yjk, zjk, xik, yik, zik, tl, tl+tl2,
411 x = x+x2+x3; y = y+y2+y3; z = z+z2+z3;
412 d = sqrt(x*x+y*y+z*z);
423 ERROR("ico_dot: n_dot(%d) and tn(%d) differ", n_dot, tn);
425 } /* end of if (tess > 1) */
427 } /* end of routine ico_dot_dod */
429 int unsp_type(int densit) {
432 while (10*i1*i1+2 < densit) i1++;
434 while (30*i2*i2+2 < densit) i2++;
435 if (10*i1*i1-2 < 30*i2*i2-2) return UNSP_ICO_ARC;
436 else return UNSP_ICO_DOD;
439 int make_unsp(int densit, int mode, int * num_dot, int cubus) {
440 int ndot, ico_cube_cb, i, j, k, l, ijk, tn, tl, tl2;
445 if (xpunsp) free(xpunsp); if (ico_wk) free(ico_wk);
447 k=1; if (mode < 0) { k=0; mode = -mode; }
448 if (mode == UNSP_ICO_ARC) { ndot = ico_dot_arc(densit); }
449 else if (mode == UNSP_ICO_DOD) { ndot = ico_dot_dod(densit); }
451 WARNING("make_unsp: mode %c%d not allowed", (k) ? '+':'-',mode);
455 last_n_dot = ndot; last_densit = densit; last_unsp = mode;
456 *num_dot=ndot; if (k) return 0;
458 /* in the following the dots of the unit sphere may be resorted */
459 last_unsp = -last_unsp;
461 /* determine distribution of points in elementary cubes */
468 while (i*i*i*2 < ndot) i++;
469 ico_cube = max(i-1, 0);
471 ico_cube_cb = ico_cube*ico_cube*ico_cube;
472 del_cube=2./((real)ico_cube);
473 work = (int *) CALLOC(ndot, sizeof(int));
475 for (l=0; l<ndot; l++) {
476 i = max((int) floor((1.+xus[3*l])/del_cube), 0);
477 if (i>=ico_cube) i = ico_cube-1;
478 j = max((int) floor((1.+xus[1+3*l])/del_cube), 0);
479 if (j>=ico_cube) j = ico_cube-1;
480 k = max((int) floor((1.+xus[2+3*l])/del_cube), 0);
481 if (k>=ico_cube) k = ico_cube-1;
482 ijk = i+j*ico_cube+k*ico_cube*ico_cube;
486 ico_wk = (int *) CALLOC(2*ico_cube_cb+1, sizeof(int));
487 ico_pt = ico_wk+ico_cube_cb;
488 for (l=0; l<ndot; l++) {
489 ico_wk[work[l]]++; /* dots per elementary cube */
492 /* reordering of the coordinate array in accordance with box number */
494 for (i=0; i<ico_cube; i++) {
495 for (j=0; j<ico_cube; j++) {
496 for (k=0; k<ico_cube; k++) {
499 ijk = i+ico_cube*j+ico_cube*ico_cube*k;
501 for (l=tl2; l<ndot; l++) {
502 if (ijk == work[l]) {
503 x = xus[3*l]; y = xus[1+3*l]; z = xus[2+3*l];
504 xus[3*l] = xus[3*tn];
505 xus[1+3*l] = xus[1+3*tn]; xus[2+3*l] = xus[2+3*tn];
506 xus[3*tn] = x; xus[1+3*tn] = y; xus[2+3*tn] = z;
507 ijk = work[l]; work[l]=work[tn]; work[tn]=ijk;
515 free(work); return 0;
519 typedef struct _stwknb {
526 int nsc_dclm_pbc(rvec *coords, real *radius, int nat,
527 int densit, int mode,
528 real *value_of_area, real **at_area,
530 real **lidots, int *nu_dots,
531 atom_id index[],int ePBC,matrix box)
533 int iat, i, ii, iii, ix, iy, iz, ixe, ixs, iye, iys, ize, izs, i_ac;
534 int jat, j, jj, jjj, jx, jy, jz;
537 int maxnei, nnei, last, maxdots=0;
538 int *wkdot=NULL, *wkbox=NULL, *wkat1=NULL, *wkatm=NULL;
540 int iii1, iii2, iiat, lfnr=0, i_at, j_at;
541 real dx, dy, dz, dd, ai, aisq, ajsq, aj, as, a;
542 real xi, yi, zi, xs=0., ys=0., zs=0.;
543 real dotarea, area, vol=0.;
544 real *xus, *dots=NULL, *atom_area=NULL;
545 int nxbox, nybox, nzbox, nxy, nxyz;
546 real xmin=0, ymin=0, zmin=0, xmax, ymax, zmax, ra2max, d, *pco;
547 /* Added DvdS 2006-07-19 */
552 distribution = unsp_type(densit);
553 if (distribution != -last_unsp || last_cubus != 4 ||
554 (densit != last_densit && densit != last_n_dot)) {
555 if (make_unsp(densit, (-distribution), &n_dot, 4))
560 dotarea = FOURPI/(real) n_dot;
564 fprintf(debug,"nsc_dclm: n_dot=%5d %9.3f\n", n_dot, dotarea);
566 /* start with neighbour list */
567 /* calculate neighbour list with the box algorithm */
569 WARNING("nsc_dclm: no surface atoms selected");
572 if (mode & FLAG_VOLUME)
574 if (mode & FLAG_DOTS) {
575 maxdots = (3*n_dot*nat)/10;
579 if (mode & FLAG_ATOM_AREA)
582 /* Compute minimum size for grid cells */
583 ra2max = radius[index[0]];
584 for (iat_xx=1; (iat_xx<nat); iat_xx++) {
586 ra2max = max(ra2max, radius[iat]);
590 /* Added DvdS 2006-07-19 */
591 /* Updated 2008-10-09 */
593 set_pbc(&pbc,ePBC,box);
595 for(i=0; (i<nat); i++) {
597 copy_rvec(coords[iat],x[i]);
599 put_atoms_in_triclinic_unitcell(ecenterTRIC,box,nat,x);
600 nxbox = max(1,floor(norm(box[XX])/ra2max));
601 nybox = max(1,floor(norm(box[YY])/ra2max));
602 nzbox = max(1,floor(norm(box[ZZ])/ra2max));
604 fprintf(debug,"nbox = %d, %d, %d\n",nxbox,nybox,nzbox);
607 /* dimensions of atomic set, cell edge is 2*ra_max */
609 xmin = coords[iat][XX]; xmax = xmin; xs=xmin;
610 ymin = coords[iat][YY]; ymax = ymin; ys=ymin;
611 zmin = coords[iat][ZZ]; zmax = zmin; zs=zmin;
612 ra2max = radius[iat];
614 for (iat_xx=1; (iat_xx<nat); iat_xx++) {
617 xmin = min(xmin, *pco); xmax = max(xmax, *pco);
618 ymin = min(ymin, *(pco+1)); ymax = max(ymax, *(pco+1));
619 zmin = min(zmin, *(pco+2)); zmax = max(zmax, *(pco+2));
620 xs= xs+ *pco; ys = ys+ *(pco+1); zs= zs+ *(pco+2);
626 fprintf(debug,"nsc_dclm: n_dot=%5d ra2max=%9.3f %9.3f\n", n_dot, ra2max, dotarea);
628 d = xmax-xmin; nxbox = (int) max(ceil(d/ra2max), 1.);
629 d = (((real)nxbox)*ra2max-d)/2.;
630 xmin = xmin-d; xmax = xmax+d;
631 d = ymax-ymin; nybox = (int) max(ceil(d/ra2max), 1.);
632 d = (((real)nybox)*ra2max-d)/2.;
633 ymin = ymin-d; ymax = ymax+d;
634 d = zmax-zmin; nzbox = (int) max(ceil(d/ra2max), 1.);
635 d = (((real)nzbox)*ra2max-d)/2.;
636 zmin = zmin-d; zmax = zmax+d;
642 /* box number of atoms */
653 for(i=0; (i<nat); i++) {
654 mvmul(box_1,x[i],x_1);
655 ix = ((int)floor(x_1[XX]*nxbox) + nxbox) % nxbox;
656 iy = ((int)floor(x_1[YY]*nybox) + nybox) % nybox;
657 iz = ((int)floor(x_1[ZZ]*nzbox) + nzbox) % nzbox;
658 j = ix + iy*nxbox + iz*nxbox*nybox;
660 fprintf(debug,"Atom %d cell index %d. x = (%8.3f,%8.3f,%8.3f) fc = (%8.3f,%8.3f,%8.3f)\n",
661 i,j,x[i][XX],x[i][YY],x[i][ZZ],x_1[XX],x_1[YY],x_1[ZZ]);
662 range_check(j,0,nxyz);
668 /* Put the atoms in their boxes */
669 for (iat_xx=0; (iat_xx<nat); iat_xx++) {
672 i = (int) max(floor((pco[XX]-xmin)/ra2max), 0);
674 j = (int) max(floor((pco[YY]-ymin)/ra2max), 0);
676 l = (int) max(floor((pco[ZZ]-zmin)/ra2max), 0);
684 /* sorting of atoms in accordance with box numbers */
686 for (i=1; i<nxyz; i++)
688 for (i=1; i<=nxyz; i++)
689 wkbox[i] += wkbox[i-1];
691 /* maxnei = (int) floor(ra2max*ra2max*ra2max*0.5); */
692 maxnei = min(nat, 27*j);
694 for (iat_xx=0; iat_xx<nat; iat_xx++) {
696 range_check(wkat1[iat_xx],0,nxyz);
697 wkatm[--wkbox[wkat1[iat_xx]]] = iat_xx;
699 fprintf(debug,"atom %5d on place %5d\n", iat, wkbox[wkat1[iat_xx]]);
703 fprintf(debug,"nsc_dclm: n_dot=%5d ra2max=%9.3f %9.3f\n",
704 n_dot, ra2max, dotarea);
705 fprintf(debug,"neighbour list calculated/box(xyz):%d %d %d\n",
706 nxbox, nybox, nzbox);
708 for (i=0; i<nxyz; i++)
709 fprintf(debug,"box %6d : atoms %4d-%4d %5d\n",
710 i, wkbox[i], wkbox[i+1]-1, wkbox[i+1]-wkbox[i]);
711 for (i=0; i<nat; i++) {
712 fprintf(debug,"list place %5d by atom %7d\n", i, index[wkatm[i]]);
716 /* calculate surface for all atoms, step cube-wise */
717 for (iz=0; iz<nzbox; iz++) {
721 ize = min(iz+2,izs+nzbox);
725 ize = min(iz+2, nzbox);
727 for (iy=0; iy<nybox; iy++) {
731 iye = min(iy+2,iys+nybox);
735 iye = min(iy+2, nybox);
737 for (ix=0; ix<nxbox; ix++) {
745 ixe = min(ix+2,ixs+nxbox);
749 ixe = min(ix+2, nxbox);
752 /* make intermediate atom list */
753 for (jz=izs; jz<ize; jz++) {
754 jjj = ((jz+nzbox) % nzbox)*nxy;
755 for (jy=iys; jy<iye; jy++) {
756 jj = ((jy+nybox) % nybox)*nxbox+jjj;
757 for (jx=ixs; jx<ixe; jx++) {
758 j = jj+((jx+nxbox) % nxbox);
759 for (jat=wkbox[j]; jat<wkbox[j+1]; jat++) {
760 range_check(wkatm[jat],0,nat);
761 range_check(iiat,0,nat);
762 wkat1[iiat] = wkatm[jat];
764 } /* end of cycle "jat" */
765 } /* end of cycle "jx" */
766 } /* end of cycle "jy" */
767 } /* end of cycle "jz" */
768 for (iat=iii1; iat<iii2; iat++) {
769 i_at = index[wkatm[iat]];
773 xi = pco[XX]; yi = pco[YY]; zi = pco[ZZ];
774 for (i=0; i<n_dot; i++)
777 ctnb = wknb; nnei = 0;
778 for (j=0; j<iiat; j++) {
779 j_at = index[wkat1[j]];
787 /* Added DvdS 2006-07-19 */
794 pbc_dx(&pbc,pco,xxi,ddx);*/
795 pbc_dx(&pbc,coords[i_at],coords[j_at],ddx);
805 dd = dx*dx+dy*dy+dz*dz;
813 ctnb->dot = (dd+aisq-ajsq)/(2.*ai); /* reference dot product */
817 /* check points on accessibility */
820 for (l=0; l<n_dot; l++) {
821 if (xus[3*l]*(wknb+last)->x+
822 xus[1+3*l]*(wknb+last)->y+
823 xus[2+3*l]*(wknb+last)->z <= (wknb+last)->dot) {
824 for (j=0; j<nnei; j++) {
825 if (xus[3*l]*(wknb+j)->x+xus[1+3*l]*(wknb+j)->y+
826 xus[2+3*l]*(wknb+j)->z > (wknb+j)->dot) {
835 } /* end of cycle j */
836 } /* end of cycle l */
840 for (l=0; l < n_dot; l++)
845 fprintf(debug,"i_ac=%d, dotarea=%8.3f, aisq=%8.3f\n",
846 i_ac, dotarea, aisq);
848 a = aisq*dotarea* (real) i_ac;
850 if (mode & FLAG_ATOM_AREA) {
851 range_check(wkatm[iat],0,nat);
852 atom_area[wkatm[iat]] = a;
854 if (mode & FLAG_DOTS) {
855 for (l=0; l<n_dot; l++) {
858 if (maxdots <= 3*lfnr+1) {
859 maxdots = maxdots+n_dot*3;
860 srenew(dots,maxdots);
862 dots[3*lfnr-3] = ai*xus[3*l]+xi;
863 dots[3*lfnr-2] = ai*xus[1+3*l]+yi;
864 dots[3*lfnr-1] = ai*xus[2+3*l]+zi;
868 if (mode & FLAG_VOLUME) {
870 for (l=0; l<n_dot; l++) {
877 vol = vol+aisq*(dx*(xi-xs)+dy*(yi-ys)+dz*(zi-zs)+ai* (real) i_ac);
880 } /* end of cycle "iat" */
881 } /* end of cycle "ix" */
882 } /* end of cycle "iy" */
883 } /* end of cycle "iz" */
893 if (mode & FLAG_VOLUME) {
894 vol = vol*FOURPI/(3.* (real) n_dot);
897 if (mode & FLAG_DOTS) {
901 if (mode & FLAG_ATOM_AREA) {
902 *at_area = atom_area;
904 *value_of_area = area;
907 fprintf(debug,"area=%8.3f\n", area);
918 real co[3*NAT], ra[NAT], area, volume, a, b, c;
925 fp = fopen("nsc.txt", "w+");
926 for (i=1; i<=NAT; i++) {
943 ra[i-1] = (1.+j*0.5)*c;
946 if (NSC(co, ra, NAT, 42, NULL, &area,
948 if (NSC(co, ra, NAT, 42, NULL, &area,
949 NULL,NULL,NULL,NULL)) ERROR("error in NSC");
951 fprintf(fp, "area : %8.3f\n", area);
954 fprintf(fp, "next call\n");
958 if (NSC(co, ra, NAT, 42,FLAG_VOLUME | FLAG_ATOM_AREA | FLAG_DOTS, &area,
960 &dots, &ndots)) ERROR("error in NSC");
963 fprintf(fp, "area : %8.3f\n", area);
964 printf("area : %8.3f\n", area);
965 fprintf(fp, "volume : %8.3f\n", volume);
966 printf("volume : %8.3f\n", volume);
967 fprintf(fp, "ndots : %8d\n", ndots);
968 printf("ndots : %8d\n", ndots);
970 for (i=1; i<=NAT; i++) {
971 fprintf(fp, "%4d ATOM %7.2f %7.2f %7.2f ra=%4.1f area=%8.3f\n",
972 i, co[3*i-3], co[3*i-2], co[3*i-1], ra[i-1], at_area[i-1]);
975 fprintf(fp, "DOTS : %8d\n", ndots);
976 for (i=1; i<=ndots; i++) {
977 fprintf(fp, "%4d DOTS %8.2f %8.2f %8.2f\n",
978 i, dots[3*i-3], dots[3*i-2], dots[3*i-1]);