2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2007, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
50 #include "gromacs/utility/smalloc.h"
59 #define UNSP_ICO_DOD 9
60 #define UNSP_ICO_ARC 10
64 int n_dot, ico_cube, last_n_dot = 0, last_densit = 0, last_unsp = 0;
67 #define FOURPI (4.*M_PI)
68 #define TORAD(A) ((A)*0.017453293)
71 #define UPDATE_FL __file__ = __FILE__, __line__ = __LINE__
72 const char * __file__; /* declared versions of macros */
73 int __line__; /* __FILE__ and __LINE__ */
78 #define ERROR UPDATE_FL, error
79 void error(const char *fmt, ...)
83 "\n---> ERROR when executing line %i in file %s !\n",
86 vfprintf(stderr, fmt, args);
88 fprintf(stderr, "\n---> Execution stopped !\n\n");
91 #define WARNING UPDATE_FL, warning2
92 void warning2(const char *fmt, ...)
96 "\n---> WARNING : line %i in file %s\n",
99 vfprintf(stderr, fmt, args);
101 fprintf(stderr, " ...!\n\n");
106 #define ASIN safe_asin
107 real safe_asin(real f)
109 if ( (fabs(f) < 1.00) )
113 if ( (fabs(f) - 1.00) <= DP_TOL)
115 ERROR("ASIN : invalid argument %f", f);
123 /* routines for dot distributions on the surface of the unit sphere */
126 void icosaeder_vertices(real *xus)
128 rh = sqrt(1.-2.*cos(TORAD(72.)))/(1.-cos(TORAD(72.)));
129 rg = cos(TORAD(72.))/(1.-cos(TORAD(72.)));
130 /* icosaeder vertices */
131 xus[ 0] = 0.; xus[ 1] = 0.; xus[ 2] = 1.;
132 xus[ 3] = rh*cos(TORAD(72.)); xus[ 4] = rh*sin(TORAD(72.)); xus[ 5] = rg;
133 xus[ 6] = rh*cos(TORAD(144.)); xus[ 7] = rh*sin(TORAD(144.)); xus[ 8] = rg;
134 xus[ 9] = rh*cos(TORAD(216.)); xus[10] = rh*sin(TORAD(216.)); xus[11] = rg;
135 xus[12] = rh*cos(TORAD(288.)); xus[13] = rh*sin(TORAD(288.)); xus[14] = rg;
136 xus[15] = rh; xus[16] = 0; xus[17] = rg;
137 xus[18] = rh*cos(TORAD(36.)); xus[19] = rh*sin(TORAD(36.)); xus[20] = -rg;
138 xus[21] = rh*cos(TORAD(108.)); xus[22] = rh*sin(TORAD(108.)); xus[23] = -rg;
139 xus[24] = -rh; xus[25] = 0; xus[26] = -rg;
140 xus[27] = rh*cos(TORAD(252.)); xus[28] = rh*sin(TORAD(252.)); xus[29] = -rg;
141 xus[30] = rh*cos(TORAD(324.)); xus[31] = rh*sin(TORAD(324.)); xus[32] = -rg;
142 xus[33] = 0.; xus[34] = 0.; xus[35] = -1.;
146 void divarc(real x1, real y1, real z1,
147 real x2, real y2, real z2,
148 int div1, int div2, real *xr, real *yr, real *zr)
151 real xd, yd, zd, dd, d1, d2, s, x, y, z;
152 real phi, sphi, cphi;
157 dd = sqrt(xd*xd+yd*yd+zd*zd);
160 ERROR("divarc: rotation axis of length %f", dd);
163 d1 = x1*x1+y1*y1+z1*z1;
166 ERROR("divarc: vector 1 of sq.length %f", d1);
168 d2 = x2*x2+y2*y2+z2*z2;
171 ERROR("divarc: vector 2 of sq.length %f", d2);
174 phi = ASIN(dd/sqrt(d1*d2));
175 phi = phi*((real)div1)/((real)div2);
176 sphi = sin(phi); cphi = cos(phi);
177 s = (x1*xd+y1*yd+z1*zd)/dd;
179 x = xd*s*(1.-cphi)/dd + x1 * cphi + (yd*z1-y1*zd)*sphi/dd;
180 y = yd*s*(1.-cphi)/dd + y1 * cphi + (zd*x1-z1*xd)*sphi/dd;
181 z = zd*s*(1.-cphi)/dd + z1 * cphi + (xd*y1-x1*yd)*sphi/dd;
182 dd = sqrt(x*x+y*y+z*z);
183 *xr = x/dd; *yr = y/dd; *zr = z/dd;
186 int ico_dot_arc(int densit) /* densit...required dots per unit sphere */
188 /* dot distribution on a unit sphere based on an icosaeder *
189 * great circle average refining of icosahedral face */
191 int i, j, k, tl, tl2, tn, tess;
192 real a, d, x, y, z, x2, y2, z2, x3, y3, z3;
193 real xij, yij, zij, xji, yji, zji, xik, yik, zik, xki, yki, zki,
194 xjk, yjk, zjk, xkj, ykj, zkj;
197 /* calculate tessalation level */
198 a = sqrt((((real) densit)-2.)/10.);
199 tess = (int) ceil(a);
200 n_dot = 10*tess*tess+2;
203 ERROR("ico_dot_arc: error in formula for tessalation level (%d->%d, %d)",
204 tess, n_dot, densit);
209 icosaeder_vertices(xus);
214 a = rh*rh*2.*(1.-cos(TORAD(72.)));
215 /* calculate tessalation of icosaeder edges */
216 for (i = 0; i < 11; i++)
218 for (j = i+1; j < 12; j++)
220 x = xus[3*i]-xus[3*j];
221 y = xus[1+3*i]-xus[1+3*j]; z = xus[2+3*i]-xus[2+3*j];
223 if (fabs(a-d) > DP_TOL)
227 for (tl = 1; tl < tess; tl++)
231 ERROR("ico_dot: tn exceeds dimension of xus");
233 divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
234 xus[3*j], xus[1+3*j], xus[2+3*j],
235 tl, tess, &xus[3*tn], &xus[1+3*tn], &xus[2+3*tn]);
240 /* calculate tessalation of icosaeder faces */
241 for (i = 0; i < 10; i++)
243 for (j = i+1; j < 11; j++)
245 x = xus[3*i]-xus[3*j];
246 y = xus[1+3*i]-xus[1+3*j]; z = xus[2+3*i]-xus[2+3*j];
248 if (fabs(a-d) > DP_TOL)
253 for (k = j+1; k < 12; k++)
255 x = xus[3*i]-xus[3*k];
256 y = xus[1+3*i]-xus[1+3*k]; z = xus[2+3*i]-xus[2+3*k];
258 if (fabs(a-d) > DP_TOL)
262 x = xus[3*j]-xus[3*k];
263 y = xus[1+3*j]-xus[1+3*k]; z = xus[2+3*j]-xus[2+3*k];
265 if (fabs(a-d) > DP_TOL)
269 for (tl = 1; tl < tess-1; tl++)
271 divarc(xus[3*j], xus[1+3*j], xus[2+3*j],
272 xus[3*i], xus[1+3*i], xus[2+3*i],
273 tl, tess, &xji, &yji, &zji);
274 divarc(xus[3*k], xus[1+3*k], xus[2+3*k],
275 xus[3*i], xus[1+3*i], xus[2+3*i],
276 tl, tess, &xki, &yki, &zki);
278 for (tl2 = 1; tl2 < tess-tl; tl2++)
280 divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
281 xus[3*j], xus[1+3*j], xus[2+3*j],
282 tl2, tess, &xij, &yij, &zij);
283 divarc(xus[3*k], xus[1+3*k], xus[2+3*k],
284 xus[3*j], xus[1+3*j], xus[2+3*j],
285 tl2, tess, &xkj, &ykj, &zkj);
286 divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
287 xus[3*k], xus[1+3*k], xus[2+3*k],
288 tess-tl-tl2, tess, &xik, &yik, &zik);
289 divarc(xus[3*j], xus[1+3*j], xus[2+3*j],
290 xus[3*k], xus[1+3*k], xus[2+3*k],
291 tess-tl-tl2, tess, &xjk, &yjk, &zjk);
294 ERROR("ico_dot: tn exceeds dimension of xus");
296 divarc(xki, yki, zki, xji, yji, zji, tl2, tess-tl,
298 divarc(xkj, ykj, zkj, xij, yij, zij, tl, tess-tl2,
300 divarc(xjk, yjk, zjk, xik, yik, zik, tl, tl+tl2,
302 x = x+x2+x3; y = y+y2+y3; z = z+z2+z3;
303 d = sqrt(x*x+y*y+z*z);
315 ERROR("ico_dot: n_dot(%d) and tn(%d) differ", n_dot, tn);
317 } /* end of if (tess > 1) */
320 } /* end of routine ico_dot_arc */
322 int ico_dot_dod(int densit) /* densit...required dots per unit sphere */
324 /* dot distribution on a unit sphere based on an icosaeder *
325 * great circle average refining of icosahedral face */
327 int i, j, k, tl, tl2, tn, tess, j1, j2;
328 real a, d, x, y, z, x2, y2, z2, x3, y3, z3, ai_d, adod;
329 real xij, yij, zij, xji, yji, zji, xik, yik, zik, xki, yki, zki,
330 xjk, yjk, zjk, xkj, ykj, zkj;
332 /* calculate tesselation level */
333 a = sqrt((((real) densit)-2.)/30.);
334 tess = max((int) ceil(a), 1);
335 n_dot = 30*tess*tess+2;
338 ERROR("ico_dot_dod: error in formula for tessalation level (%d->%d, %d)",
339 tess, n_dot, densit);
344 icosaeder_vertices(xus);
347 /* square of the edge of an icosaeder */
348 a = rh*rh*2.*(1.-cos(TORAD(72.)));
349 /* dodecaeder vertices */
350 for (i = 0; i < 10; i++)
352 for (j = i+1; j < 11; j++)
354 x = xus[3*i]-xus[3*j];
355 y = xus[1+3*i]-xus[1+3*j]; z = xus[2+3*i]-xus[2+3*j];
357 if (fabs(a-d) > DP_TOL)
361 for (k = j+1; k < 12; k++)
363 x = xus[3*i]-xus[3*k];
364 y = xus[1+3*i]-xus[1+3*k]; z = xus[2+3*i]-xus[2+3*k];
366 if (fabs(a-d) > DP_TOL)
370 x = xus[3*j]-xus[3*k];
371 y = xus[1+3*j]-xus[1+3*k]; z = xus[2+3*j]-xus[2+3*k];
373 if (fabs(a-d) > DP_TOL)
377 x = xus[ 3*i]+xus[ 3*j]+xus[ 3*k];
378 y = xus[1+3*i]+xus[1+3*j]+xus[1+3*k];
379 z = xus[2+3*i]+xus[2+3*j]+xus[2+3*k];
380 d = sqrt(x*x+y*y+z*z);
381 xus[3*tn] = x/d; xus[1+3*tn] = y/d; xus[2+3*tn] = z/d;
390 /* square of the edge of an dodecaeder */
391 adod = 4.*(cos(TORAD(108.))-cos(TORAD(120.)))/(1.-cos(TORAD(120.)));
392 /* square of the distance of two adjacent vertices of ico- and dodecaeder */
393 ai_d = 2.*(1.-sqrt(1.-a/3.));
395 /* calculate tessalation of mixed edges */
396 for (i = 0; i < 31; i++)
398 j1 = 12; j2 = 32; a = ai_d;
403 for (j = j1; j < j2; j++)
405 x = xus[3*i]-xus[3*j];
406 y = xus[1+3*i]-xus[1+3*j]; z = xus[2+3*i]-xus[2+3*j];
408 if (fabs(a-d) > DP_TOL)
412 for (tl = 1; tl < tess; tl++)
416 ERROR("ico_dot: tn exceeds dimension of xus");
418 divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
419 xus[3*j], xus[1+3*j], xus[2+3*j],
420 tl, tess, &xus[3*tn], &xus[1+3*tn], &xus[2+3*tn]);
425 /* calculate tessalation of pentakisdodecahedron faces */
426 for (i = 0; i < 12; i++)
428 for (j = 12; j < 31; j++)
430 x = xus[3*i]-xus[3*j];
431 y = xus[1+3*i]-xus[1+3*j]; z = xus[2+3*i]-xus[2+3*j];
433 if (fabs(ai_d-d) > DP_TOL)
438 for (k = j+1; k < 32; k++)
440 x = xus[3*i]-xus[3*k];
441 y = xus[1+3*i]-xus[1+3*k]; z = xus[2+3*i]-xus[2+3*k];
443 if (fabs(ai_d-d) > DP_TOL)
447 x = xus[3*j]-xus[3*k];
448 y = xus[1+3*j]-xus[1+3*k]; z = xus[2+3*j]-xus[2+3*k];
450 if (fabs(adod-d) > DP_TOL)
454 for (tl = 1; tl < tess-1; tl++)
456 divarc(xus[3*j], xus[1+3*j], xus[2+3*j],
457 xus[3*i], xus[1+3*i], xus[2+3*i],
458 tl, tess, &xji, &yji, &zji);
459 divarc(xus[3*k], xus[1+3*k], xus[2+3*k],
460 xus[3*i], xus[1+3*i], xus[2+3*i],
461 tl, tess, &xki, &yki, &zki);
463 for (tl2 = 1; tl2 < tess-tl; tl2++)
465 divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
466 xus[3*j], xus[1+3*j], xus[2+3*j],
467 tl2, tess, &xij, &yij, &zij);
468 divarc(xus[3*k], xus[1+3*k], xus[2+3*k],
469 xus[3*j], xus[1+3*j], xus[2+3*j],
470 tl2, tess, &xkj, &ykj, &zkj);
471 divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
472 xus[3*k], xus[1+3*k], xus[2+3*k],
473 tess-tl-tl2, tess, &xik, &yik, &zik);
474 divarc(xus[3*j], xus[1+3*j], xus[2+3*j],
475 xus[3*k], xus[1+3*k], xus[2+3*k],
476 tess-tl-tl2, tess, &xjk, &yjk, &zjk);
479 ERROR("ico_dot: tn exceeds dimension of xus");
481 divarc(xki, yki, zki, xji, yji, zji, tl2, tess-tl,
483 divarc(xkj, ykj, zkj, xij, yij, zij, tl, tess-tl2,
485 divarc(xjk, yjk, zjk, xik, yik, zik, tl, tl+tl2,
487 x = x+x2+x3; y = y+y2+y3; z = z+z2+z3;
488 d = sqrt(x*x+y*y+z*z);
500 ERROR("ico_dot: n_dot(%d) and tn(%d) differ", n_dot, tn);
502 } /* end of if (tess > 1) */
505 } /* end of routine ico_dot_dod */
507 int unsp_type(int densit)
511 while (10*i1*i1+2 < densit)
516 while (30*i2*i2+2 < densit)
520 if (10*i1*i1-2 < 30*i2*i2-2)
530 int make_unsp(int densit, int mode, int * num_dot, int cubus)
532 int *ico_wk, *ico_pt;
533 int ndot, ico_cube_cb, i, j, k, l, ijk, tn, tl, tl2;
547 if (mode == UNSP_ICO_ARC)
549 ndot = ico_dot_arc(densit);
551 else if (mode == UNSP_ICO_DOD)
553 ndot = ico_dot_dod(densit);
557 WARNING("make_unsp: mode %c%d not allowed", (k) ? '+' : '-', mode);
561 last_n_dot = ndot; last_densit = densit; last_unsp = mode;
562 *num_dot = ndot; if (k)
567 /* in the following the dots of the unit sphere may be resorted */
568 last_unsp = -last_unsp;
570 /* determine distribution of points in elementary cubes */
579 while (i*i*i*2 < ndot)
583 ico_cube = max(i-1, 0);
585 ico_cube_cb = ico_cube*ico_cube*ico_cube;
586 del_cube = 2./((real)ico_cube);
589 for (l = 0; l < ndot; l++)
591 i = max((int) floor((1.+xus[3*l])/del_cube), 0);
596 j = max((int) floor((1.+xus[1+3*l])/del_cube), 0);
601 k = max((int) floor((1.+xus[2+3*l])/del_cube), 0);
606 ijk = i+j*ico_cube+k*ico_cube*ico_cube;
610 snew(ico_wk, 2*ico_cube_cb+1);
612 ico_pt = ico_wk+ico_cube_cb;
613 for (l = 0; l < ndot; l++)
615 ico_wk[work[l]]++; /* dots per elementary cube */
618 /* reordering of the coordinate array in accordance with box number */
620 for (i = 0; i < ico_cube; i++)
622 for (j = 0; j < ico_cube; j++)
624 for (k = 0; k < ico_cube; k++)
628 ijk = i+ico_cube*j+ico_cube*ico_cube*k;
630 for (l = tl2; l < ndot; l++)
634 x = xus[3*l]; y = xus[1+3*l]; z = xus[2+3*l];
635 xus[3*l] = xus[3*tn];
636 xus[1+3*l] = xus[1+3*tn]; xus[2+3*l] = xus[2+3*tn];
637 xus[3*tn] = x; xus[1+3*tn] = y; xus[2+3*tn] = z;
638 ijk = work[l]; work[l] = work[tn]; work[tn] = ijk;
653 typedef struct _stwknb {
660 int nsc_dclm_pbc(const rvec *coords, real *radius, int nat,
661 int densit, int mode,
662 real *value_of_area, real **at_area,
664 real **lidots, int *nu_dots,
665 atom_id index[], int ePBC, matrix box)
667 int iat, i, ii, iii, ix, iy, iz, ixe, ixs, iye, iys, ize, izs, i_ac;
668 int jat, j, jj, jjj, jx, jy, jz;
671 int maxnei, nnei, last, maxdots = 0;
672 int *wkdot = NULL, *wkbox = NULL, *wkat1 = NULL, *wkatm = NULL;
674 int iii1, iii2, iiat, lfnr = 0, i_at, j_at;
675 real dx, dy, dz, dd, ai, aisq, ajsq, aj, as, a;
676 real xi, yi, zi, xs = 0., ys = 0., zs = 0.;
677 real dotarea, area, vol = 0.;
678 real *xus, *dots = NULL, *atom_area = NULL;
679 int nxbox, nybox, nzbox, nxy, nxyz;
680 real xmin = 0, ymin = 0, zmin = 0, xmax, ymax, zmax, ra2max, d;
682 /* Added DvdS 2006-07-19 */
687 distribution = unsp_type(densit);
688 if (distribution != -last_unsp || last_cubus != 4 ||
689 (densit != last_densit && densit != last_n_dot))
691 if (make_unsp(densit, (-distribution), &n_dot, 4))
698 dotarea = FOURPI/(real) n_dot;
703 fprintf(debug, "nsc_dclm: n_dot=%5d %9.3f\n", n_dot, dotarea);
706 /* start with neighbour list */
707 /* calculate neighbour list with the box algorithm */
710 WARNING("nsc_dclm: no surface atoms selected");
713 if (mode & FLAG_VOLUME)
717 if (mode & FLAG_DOTS)
719 maxdots = (3*n_dot*nat)/10;
720 /* should be set to NULL on first user call */
727 srenew(dots, maxdots);
732 if (mode & FLAG_ATOM_AREA)
734 /* should be set to NULL on first user call */
735 if (atom_area == NULL)
737 snew(atom_area, nat);
741 srenew(atom_area, nat);
745 /* Compute minimum size for grid cells */
746 ra2max = radius[index[0]];
747 for (iat_xx = 1; (iat_xx < nat); iat_xx++)
750 ra2max = max(ra2max, radius[iat]);
754 /* Added DvdS 2006-07-19 */
755 /* Updated 2008-10-09 */
758 set_pbc(&pbc, ePBC, box);
760 for (i = 0; (i < nat); i++)
763 copy_rvec(coords[iat], x[i]);
765 put_atoms_in_triclinic_unitcell(ecenterTRIC, box, nat, x);
766 nxbox = max(1, floor(norm(box[XX])/ra2max));
767 nybox = max(1, floor(norm(box[YY])/ra2max));
768 nzbox = max(1, floor(norm(box[ZZ])/ra2max));
771 fprintf(debug, "nbox = %d, %d, %d\n", nxbox, nybox, nzbox);
776 /* dimensions of atomic set, cell edge is 2*ra_max */
778 xmin = coords[iat][XX]; xmax = xmin; xs = xmin;
779 ymin = coords[iat][YY]; ymax = ymin; ys = ymin;
780 zmin = coords[iat][ZZ]; zmax = zmin; zs = zmin;
782 for (iat_xx = 1; (iat_xx < nat); iat_xx++)
786 xmin = min(xmin, *pco); xmax = max(xmax, *pco);
787 ymin = min(ymin, *(pco+1)); ymax = max(ymax, *(pco+1));
788 zmin = min(zmin, *(pco+2)); zmax = max(zmax, *(pco+2));
789 xs = xs+ *pco; ys = ys+ *(pco+1); zs = zs+ *(pco+2);
796 fprintf(debug, "nsc_dclm: n_dot=%5d ra2max=%9.3f %9.3f\n", n_dot, ra2max, dotarea);
799 d = xmax-xmin; nxbox = (int) max(ceil(d/ra2max), 1.);
800 d = (((real)nxbox)*ra2max-d)/2.;
801 xmin = xmin-d; xmax = xmax+d;
802 d = ymax-ymin; nybox = (int) max(ceil(d/ra2max), 1.);
803 d = (((real)nybox)*ra2max-d)/2.;
804 ymin = ymin-d; ymax = ymax+d;
805 d = zmax-zmin; nzbox = (int) max(ceil(d/ra2max), 1.);
806 d = (((real)nzbox)*ra2max-d)/2.;
807 zmin = zmin-d; zmax = zmax+d;
813 /* box number of atoms */
825 for (i = 0; (i < nat); i++)
827 mvmul(box_1, x[i], x_1);
828 ix = ((int)floor(x_1[XX]*nxbox) + 2*nxbox) % nxbox;
829 iy = ((int)floor(x_1[YY]*nybox) + 2*nybox) % nybox;
830 iz = ((int)floor(x_1[ZZ]*nzbox) + 2*nzbox) % nzbox;
831 j = ix + iy*nxbox + iz*nxbox*nybox;
834 fprintf(debug, "Atom %d cell index %d. x = (%8.3f,%8.3f,%8.3f) fc = (%8.3f,%8.3f,%8.3f)\n",
835 i, j, x[i][XX], x[i][YY], x[i][ZZ], x_1[XX], x_1[YY], x_1[ZZ]);
837 range_check(j, 0, nxyz);
844 /* Put the atoms in their boxes */
845 for (iat_xx = 0; (iat_xx < nat); iat_xx++)
849 i = (int) max(floor((pco[XX]-xmin)/ra2max), 0);
851 j = (int) max(floor((pco[YY]-ymin)/ra2max), 0);
853 l = (int) max(floor((pco[ZZ]-zmin)/ra2max), 0);
861 /* sorting of atoms in accordance with box numbers */
863 for (i = 1; i < nxyz; i++)
865 j = max(wkbox[i], j);
867 for (i = 1; i <= nxyz; i++)
869 wkbox[i] += wkbox[i-1];
872 /* maxnei = (int) floor(ra2max*ra2max*ra2max*0.5); */
873 maxnei = min(nat, 27*j);
875 for (iat_xx = 0; iat_xx < nat; iat_xx++)
878 range_check(wkat1[iat_xx], 0, nxyz);
879 wkatm[--wkbox[wkat1[iat_xx]]] = iat_xx;
882 fprintf(debug, "atom %5d on place %5d\n", iat, wkbox[wkat1[iat_xx]]);
888 fprintf(debug, "nsc_dclm: n_dot=%5d ra2max=%9.3f %9.3f\n",
889 n_dot, ra2max, dotarea);
890 fprintf(debug, "neighbour list calculated/box(xyz):%d %d %d\n",
891 nxbox, nybox, nzbox);
893 for (i = 0; i < nxyz; i++)
895 fprintf(debug, "box %6d : atoms %4d-%4d %5d\n",
896 i, wkbox[i], wkbox[i+1]-1, wkbox[i+1]-wkbox[i]);
898 for (i = 0; i < nat; i++)
900 fprintf(debug, "list place %5d by atom %7d\n", i, index[wkatm[i]]);
904 /* calculate surface for all atoms, step cube-wise */
905 for (iz = 0; iz < nzbox; iz++)
911 ize = min(iz+2, izs+nzbox);
916 ize = min(iz+2, nzbox);
918 for (iy = 0; iy < nybox; iy++)
924 iye = min(iy+2, iys+nybox);
929 iye = min(iy+2, nybox);
931 for (ix = 0; ix < nxbox; ix++)
943 ixe = min(ix+2, ixs+nxbox);
948 ixe = min(ix+2, nxbox);
951 /* make intermediate atom list */
952 for (jz = izs; jz < ize; jz++)
954 jjj = ((jz+nzbox) % nzbox)*nxy;
955 for (jy = iys; jy < iye; jy++)
957 jj = ((jy+nybox) % nybox)*nxbox+jjj;
958 for (jx = ixs; jx < ixe; jx++)
960 j = jj+((jx+nxbox) % nxbox);
961 for (jat = wkbox[j]; jat < wkbox[j+1]; jat++)
963 range_check(wkatm[jat], 0, nat);
964 range_check(iiat, 0, nat);
965 wkat1[iiat] = wkatm[jat];
967 } /* end of cycle "jat" */
968 } /* end of cycle "jx" */
969 } /* end of cycle "jy" */
970 } /* end of cycle "jz" */
971 for (iat = iii1; iat < iii2; iat++)
973 i_at = index[wkatm[iat]];
977 xi = pco[XX]; yi = pco[YY]; zi = pco[ZZ];
978 for (i = 0; i < n_dot; i++)
983 ctnb = wknb; nnei = 0;
984 for (j = 0; j < iiat; j++)
986 j_at = index[wkat1[j]];
996 /* Added DvdS 2006-07-19 */
1004 pbc_dx(&pbc,pco,xxi,ddx);*/
1005 pbc_dx(&pbc, coords[j_at], coords[i_at], ddx);
1016 dd = dx*dx+dy*dy+dz*dz;
1026 ctnb->dot = (dd+aisq-ajsq)/(2.*ai); /* reference dot product */
1030 /* check points on accessibility */
1034 for (l = 0; l < n_dot; l++)
1036 if (xus[3*l]*(wknb+last)->x+
1037 xus[1+3*l]*(wknb+last)->y+
1038 xus[2+3*l]*(wknb+last)->z <= (wknb+last)->dot)
1040 for (j = 0; j < nnei; j++)
1042 if (xus[3*l]*(wknb+j)->x+xus[1+3*l]*(wknb+j)->y+
1043 xus[2+3*l]*(wknb+j)->z > (wknb+j)->dot)
1054 } /* end of cycle j */
1055 } /* end of cycle l */
1060 for (l = 0; l < n_dot; l++)
1068 fprintf(debug, "i_ac=%d, dotarea=%8.3f, aisq=%8.3f\n",
1069 i_ac, dotarea, aisq);
1072 a = aisq*dotarea* (real) i_ac;
1074 if (mode & FLAG_ATOM_AREA)
1076 range_check(wkatm[iat], 0, nat);
1077 atom_area[wkatm[iat]] = a;
1079 if (mode & FLAG_DOTS)
1081 for (l = 0; l < n_dot; l++)
1086 if (maxdots <= 3*lfnr+1)
1088 maxdots = maxdots+n_dot*3;
1089 srenew(dots, maxdots);
1091 dots[3*lfnr-3] = ai*xus[3*l]+xi;
1092 dots[3*lfnr-2] = ai*xus[1+3*l]+yi;
1093 dots[3*lfnr-1] = ai*xus[2+3*l]+zi;
1097 if (mode & FLAG_VOLUME)
1099 dx = 0.; dy = 0.; dz = 0.;
1100 for (l = 0; l < n_dot; l++)
1109 vol = vol+aisq*(dx*(xi-xs)+dy*(yi-ys)+dz*(zi-zs)+ai* (real) i_ac);
1112 } /* end of cycle "iat" */
1113 } /* end of cycle "ix" */
1114 } /* end of cycle "iy" */
1115 } /* end of cycle "iz" */
1126 if (mode & FLAG_VOLUME)
1128 vol = vol*FOURPI/(3.* (real) n_dot);
1129 *value_of_vol = vol;
1131 if (mode & FLAG_DOTS)
1136 if (mode & FLAG_ATOM_AREA)
1138 *at_area = atom_area;
1140 *value_of_area = area;
1144 fprintf(debug, "area=%8.3f\n", area);
1156 real co[3*NAT], ra[NAT], area, volume, a, b, c;
1163 fp = fopen("nsc.txt", "w+");
1164 for (i = 1; i <= NAT; i++)
1178 if (i%3 == 0) b=0.5;
1182 ra[i-1] = (1.+j*0.5)*c;
1185 if (NSC(co, ra, NAT, 42, NULL, &area,
1187 if (NSC(co, ra, NAT, 42, NULL, &area,
1188 NULL, NULL, NULL, NULL))
1190 ERROR("error in NSC");
1193 fprintf(fp, "area : %8.3f\n", area);
1196 fprintf(fp, "next call\n");
1200 if (NSC(co, ra, NAT, 42, FLAG_VOLUME | FLAG_ATOM_AREA | FLAG_DOTS, &area,
1204 ERROR("error in NSC");
1208 fprintf(fp, "area : %8.3f\n", area);
1209 printf("area : %8.3f\n", area);
1210 fprintf(fp, "volume : %8.3f\n", volume);
1211 printf("volume : %8.3f\n", volume);
1212 fprintf(fp, "ndots : %8d\n", ndots);
1213 printf("ndots : %8d\n", ndots);
1215 for (i = 1; i <= NAT; i++)
1217 fprintf(fp, "%4d ATOM %7.2f %7.2f %7.2f ra=%4.1f area=%8.3f\n",
1218 i, co[3*i-3], co[3*i-2], co[3*i-1], ra[i-1], at_area[i-1]);
1221 fprintf(fp, "DOTS : %8d\n", ndots);
1222 for (i = 1; i <= ndots; i++)
1224 fprintf(fp, "%4d DOTS %8.2f %8.2f %8.2f\n",
1225 i, dots[3*i-3], dots[3*i-2], dots[3*i-1]);