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