File: | gromacs/trajectoryanalysis/modules/nsc.c |
Location: | line 807, column 24 |
Description: | Value stored to 'zmax' is never read |
1 | /* |
2 | * This file is part of the GROMACS molecular simulation package. |
3 | * |
4 | * Copyright (c) 1991-2000, University of Groningen, The Netherlands. |
5 | * Copyright (c) 2001-2007, The GROMACS development team. |
6 | * 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. |
10 | * |
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. |
15 | * |
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. |
20 | * |
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. |
25 | * |
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. |
33 | * |
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. |
36 | */ |
37 | #ifdef HAVE_CONFIG_H1 |
38 | #include <config.h> |
39 | #endif |
40 | |
41 | #include <stdio.h> |
42 | #include <string.h> |
43 | #include <stdlib.h> |
44 | #include <math.h> |
45 | #include <stdarg.h> |
46 | /* Modified DvdS */ |
47 | #include "pbc.h" |
48 | #include "macros.h" |
49 | #include "gromacs/math/vec.h" |
50 | #include "gromacs/utility/smalloc.h" |
51 | #include "nsc.h" |
52 | |
53 | #define TEST_NSC0 0 |
54 | |
55 | #define TEST_ARC0 0 |
56 | #define TEST_DOD0 0 |
57 | #define TEST_CUBE0 0 |
58 | |
59 | #define UNSP_ICO_DOD9 9 |
60 | #define UNSP_ICO_ARC10 10 |
61 | |
62 | real *xpunsp = NULL((void*)0); |
63 | real del_cube; |
64 | int n_dot, ico_cube, last_n_dot = 0, last_densit = 0, last_unsp = 0; |
65 | int last_cubus = 0; |
66 | |
67 | #define FOURPI(4.*3.14159265358979323846) (4.*M_PI3.14159265358979323846) |
68 | #define TORAD(A)((A)*0.017453293) ((A)*0.017453293) |
69 | #define DP_TOL0.001 0.001 |
70 | |
71 | #define UPDATE_FL__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 71 __file__ = __FILE__"/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c", __line__ = __LINE__71 |
72 | const char * __file__; /* declared versions of macros */ |
73 | int __line__; /* __FILE__ and __LINE__ */ |
74 | |
75 | #ifdef ERROR__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 75, error |
76 | #undef ERROR__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 76, error |
77 | #endif |
78 | #define ERROR__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 78, error UPDATE_FL__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 78, error |
79 | void error(const char *fmt, ...) |
80 | { |
81 | va_list args; |
82 | fprintf(stderrstderr, |
83 | "\n---> ERROR when executing line %i in file %s !\n", |
84 | __line__, __file__); |
85 | va_start(args, fmt)__builtin_va_start(args, fmt); |
86 | vfprintf(stderrstderr, fmt, args); |
87 | va_end(args)__builtin_va_end(args); |
88 | fprintf(stderrstderr, "\n---> Execution stopped !\n\n"); |
89 | } |
90 | |
91 | #define WARNING__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 91, warning2 UPDATE_FL__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 91, warning2 |
92 | void warning2(const char *fmt, ...) |
93 | { |
94 | va_list args; |
95 | fprintf(stderrstderr, |
96 | "\n---> WARNING : line %i in file %s\n", |
97 | __line__, __file__); |
98 | va_start(args, fmt)__builtin_va_start(args, fmt); |
99 | vfprintf(stderrstderr, fmt, args); |
100 | va_end(args)__builtin_va_end(args); |
101 | fprintf(stderrstderr, " ...!\n\n"); |
102 | fflush(stderrstderr); |
103 | fflush(stdoutstdout); |
104 | } |
105 | |
106 | #define ASINsafe_asin safe_asin |
107 | real safe_asin(real f) |
108 | { |
109 | if ( (fabs(f) < 1.00) ) |
110 | { |
111 | return( asin(f) ); |
112 | } |
113 | if ( (fabs(f) - 1.00) <= DP_TOL0.001) |
114 | { |
115 | ERROR__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 115, error("ASIN : invalid argument %f", f); |
116 | } |
117 | return(M_PI_21.57079632679489661923); |
118 | } |
119 | |
120 | |
121 | |
122 | |
123 | /* routines for dot distributions on the surface of the unit sphere */ |
124 | real rg, rh; |
125 | |
126 | void icosaeder_vertices(real *xus) |
127 | { |
128 | rh = sqrt(1.-2.*cos(TORAD(72.)((72.)*0.017453293)))/(1.-cos(TORAD(72.)((72.)*0.017453293))); |
129 | rg = cos(TORAD(72.)((72.)*0.017453293))/(1.-cos(TORAD(72.)((72.)*0.017453293))); |
130 | /* icosaeder vertices */ |
131 | xus[ 0] = 0.; xus[ 1] = 0.; xus[ 2] = 1.; |
132 | xus[ 3] = rh*cos(TORAD(72.)((72.)*0.017453293)); xus[ 4] = rh*sin(TORAD(72.)((72.)*0.017453293)); xus[ 5] = rg; |
133 | xus[ 6] = rh*cos(TORAD(144.)((144.)*0.017453293)); xus[ 7] = rh*sin(TORAD(144.)((144.)*0.017453293)); xus[ 8] = rg; |
134 | xus[ 9] = rh*cos(TORAD(216.)((216.)*0.017453293)); xus[10] = rh*sin(TORAD(216.)((216.)*0.017453293)); xus[11] = rg; |
135 | xus[12] = rh*cos(TORAD(288.)((288.)*0.017453293)); xus[13] = rh*sin(TORAD(288.)((288.)*0.017453293)); xus[14] = rg; |
136 | xus[15] = rh; xus[16] = 0; xus[17] = rg; |
137 | xus[18] = rh*cos(TORAD(36.)((36.)*0.017453293)); xus[19] = rh*sin(TORAD(36.)((36.)*0.017453293)); xus[20] = -rg; |
138 | xus[21] = rh*cos(TORAD(108.)((108.)*0.017453293)); xus[22] = rh*sin(TORAD(108.)((108.)*0.017453293)); xus[23] = -rg; |
139 | xus[24] = -rh; xus[25] = 0; xus[26] = -rg; |
140 | xus[27] = rh*cos(TORAD(252.)((252.)*0.017453293)); xus[28] = rh*sin(TORAD(252.)((252.)*0.017453293)); xus[29] = -rg; |
141 | xus[30] = rh*cos(TORAD(324.)((324.)*0.017453293)); xus[31] = rh*sin(TORAD(324.)((324.)*0.017453293)); xus[32] = -rg; |
142 | xus[33] = 0.; xus[34] = 0.; xus[35] = -1.; |
143 | } |
144 | |
145 | |
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) |
149 | { |
150 | |
151 | real xd, yd, zd, dd, d1, d2, s, x, y, z; |
152 | real phi, sphi, cphi; |
153 | |
154 | xd = y1*z2-y2*z1; |
155 | yd = z1*x2-z2*x1; |
156 | zd = x1*y2-x2*y1; |
157 | dd = sqrt(xd*xd+yd*yd+zd*zd); |
158 | if (dd < DP_TOL0.001) |
159 | { |
160 | ERROR__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 160, error("divarc: rotation axis of length %f", dd); |
161 | } |
162 | |
163 | d1 = x1*x1+y1*y1+z1*z1; |
164 | if (d1 < 0.5) |
165 | { |
166 | ERROR__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 166, error("divarc: vector 1 of sq.length %f", d1); |
167 | } |
168 | d2 = x2*x2+y2*y2+z2*z2; |
169 | if (d2 < 0.5) |
170 | { |
171 | ERROR__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 171, error("divarc: vector 2 of sq.length %f", d2); |
172 | } |
173 | |
174 | phi = ASINsafe_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; |
178 | |
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; |
184 | } |
185 | |
186 | int ico_dot_arc(int densit) /* densit...required dots per unit sphere */ |
187 | { |
188 | /* dot distribution on a unit sphere based on an icosaeder * |
189 | * great circle average refining of icosahedral face */ |
190 | |
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; |
195 | real *xus = NULL((void*)0); |
196 | |
197 | /* calculate tessalation level */ |
198 | a = sqrt((((real) densit)-2.)/10.); |
199 | tess = (int) ceil(a); |
200 | n_dot = 10*tess*tess+2; |
201 | if (n_dot < densit) |
202 | { |
203 | ERROR__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 203, error("ico_dot_arc: error in formula for tessalation level (%d->%d, %d)", |
204 | tess, n_dot, densit); |
205 | } |
206 | |
207 | snew(xus, 3*n_dot)(xus) = save_calloc("xus", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 207, (3*n_dot), sizeof(*(xus))); |
208 | xpunsp = xus; |
209 | icosaeder_vertices(xus); |
210 | |
211 | if (tess > 1) |
212 | { |
213 | tn = 12; |
214 | a = rh*rh*2.*(1.-cos(TORAD(72.)((72.)*0.017453293))); |
215 | /* calculate tessalation of icosaeder edges */ |
216 | for (i = 0; i < 11; i++) |
217 | { |
218 | for (j = i+1; j < 12; j++) |
219 | { |
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]; |
222 | d = x*x+y*y+z*z; |
223 | if (fabs(a-d) > DP_TOL0.001) |
224 | { |
225 | continue; |
226 | } |
227 | for (tl = 1; tl < tess; tl++) |
228 | { |
229 | if (tn >= n_dot) |
230 | { |
231 | ERROR__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 231, error("ico_dot: tn exceeds dimension of xus"); |
232 | } |
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]); |
236 | tn++; |
237 | } |
238 | } |
239 | } |
240 | /* calculate tessalation of icosaeder faces */ |
241 | for (i = 0; i < 10; i++) |
242 | { |
243 | for (j = i+1; j < 11; j++) |
244 | { |
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]; |
247 | d = x*x+y*y+z*z; |
248 | if (fabs(a-d) > DP_TOL0.001) |
249 | { |
250 | continue; |
251 | } |
252 | |
253 | for (k = j+1; k < 12; k++) |
254 | { |
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]; |
257 | d = x*x+y*y+z*z; |
258 | if (fabs(a-d) > DP_TOL0.001) |
259 | { |
260 | continue; |
261 | } |
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]; |
264 | d = x*x+y*y+z*z; |
265 | if (fabs(a-d) > DP_TOL0.001) |
266 | { |
267 | continue; |
268 | } |
269 | for (tl = 1; tl < tess-1; tl++) |
270 | { |
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); |
277 | |
278 | for (tl2 = 1; tl2 < tess-tl; tl2++) |
279 | { |
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); |
292 | if (tn >= n_dot) |
293 | { |
294 | ERROR__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 294, error("ico_dot: tn exceeds dimension of xus"); |
295 | } |
296 | divarc(xki, yki, zki, xji, yji, zji, tl2, tess-tl, |
297 | &x, &y, &z); |
298 | divarc(xkj, ykj, zkj, xij, yij, zij, tl, tess-tl2, |
299 | &x2, &y2, &z2); |
300 | divarc(xjk, yjk, zjk, xik, yik, zik, tl, tl+tl2, |
301 | &x3, &y3, &z3); |
302 | x = x+x2+x3; y = y+y2+y3; z = z+z2+z3; |
303 | d = sqrt(x*x+y*y+z*z); |
304 | xus[3*tn] = x/d; |
305 | xus[1+3*tn] = y/d; |
306 | xus[2+3*tn] = z/d; |
307 | tn++; |
308 | } /* cycle tl2 */ |
309 | } /* cycle tl */ |
310 | } /* cycle k */ |
311 | } /* cycle j */ |
312 | } /* cycle i */ |
313 | if (n_dot != tn) |
314 | { |
315 | ERROR__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 315, error("ico_dot: n_dot(%d) and tn(%d) differ", n_dot, tn); |
316 | } |
317 | } /* end of if (tess > 1) */ |
318 | |
319 | return n_dot; |
320 | } /* end of routine ico_dot_arc */ |
321 | |
322 | int ico_dot_dod(int densit) /* densit...required dots per unit sphere */ |
323 | { |
324 | /* dot distribution on a unit sphere based on an icosaeder * |
325 | * great circle average refining of icosahedral face */ |
326 | |
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; |
331 | real *xus = NULL((void*)0); |
332 | /* calculate tesselation level */ |
333 | a = sqrt((((real) densit)-2.)/30.); |
334 | tess = max((int) ceil(a), 1)((((int) ceil(a)) > (1)) ? ((int) ceil(a)) : (1) ); |
335 | n_dot = 30*tess*tess+2; |
336 | if (n_dot < densit) |
337 | { |
338 | ERROR__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 338, error("ico_dot_dod: error in formula for tessalation level (%d->%d, %d)", |
339 | tess, n_dot, densit); |
340 | } |
341 | |
342 | snew(xus, 3*n_dot)(xus) = save_calloc("xus", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 342, (3*n_dot), sizeof(*(xus))); |
343 | xpunsp = xus; |
344 | icosaeder_vertices(xus); |
345 | |
346 | tn = 12; |
347 | /* square of the edge of an icosaeder */ |
348 | a = rh*rh*2.*(1.-cos(TORAD(72.)((72.)*0.017453293))); |
349 | /* dodecaeder vertices */ |
350 | for (i = 0; i < 10; i++) |
351 | { |
352 | for (j = i+1; j < 11; j++) |
353 | { |
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]; |
356 | d = x*x+y*y+z*z; |
357 | if (fabs(a-d) > DP_TOL0.001) |
358 | { |
359 | continue; |
360 | } |
361 | for (k = j+1; k < 12; k++) |
362 | { |
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]; |
365 | d = x*x+y*y+z*z; |
366 | if (fabs(a-d) > DP_TOL0.001) |
367 | { |
368 | continue; |
369 | } |
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]; |
372 | d = x*x+y*y+z*z; |
373 | if (fabs(a-d) > DP_TOL0.001) |
374 | { |
375 | continue; |
376 | } |
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; |
382 | tn++; |
383 | } |
384 | } |
385 | } |
386 | |
387 | if (tess > 1) |
388 | { |
389 | tn = 32; |
390 | /* square of the edge of an dodecaeder */ |
391 | adod = 4.*(cos(TORAD(108.)((108.)*0.017453293))-cos(TORAD(120.)((120.)*0.017453293)))/(1.-cos(TORAD(120.)((120.)*0.017453293))); |
392 | /* square of the distance of two adjacent vertices of ico- and dodecaeder */ |
393 | ai_d = 2.*(1.-sqrt(1.-a/3.)); |
394 | |
395 | /* calculate tessalation of mixed edges */ |
396 | for (i = 0; i < 31; i++) |
397 | { |
398 | j1 = 12; j2 = 32; a = ai_d; |
399 | if (i >= 12) |
400 | { |
401 | j1 = i+1; a = adod; |
402 | } |
403 | for (j = j1; j < j2; j++) |
404 | { |
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]; |
407 | d = x*x+y*y+z*z; |
408 | if (fabs(a-d) > DP_TOL0.001) |
409 | { |
410 | continue; |
411 | } |
412 | for (tl = 1; tl < tess; tl++) |
413 | { |
414 | if (tn >= n_dot) |
415 | { |
416 | ERROR__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 416, error("ico_dot: tn exceeds dimension of xus"); |
417 | } |
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]); |
421 | tn++; |
422 | } |
423 | } |
424 | } |
425 | /* calculate tessalation of pentakisdodecahedron faces */ |
426 | for (i = 0; i < 12; i++) |
427 | { |
428 | for (j = 12; j < 31; j++) |
429 | { |
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]; |
432 | d = x*x+y*y+z*z; |
433 | if (fabs(ai_d-d) > DP_TOL0.001) |
434 | { |
435 | continue; |
436 | } |
437 | |
438 | for (k = j+1; k < 32; k++) |
439 | { |
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]; |
442 | d = x*x+y*y+z*z; |
443 | if (fabs(ai_d-d) > DP_TOL0.001) |
444 | { |
445 | continue; |
446 | } |
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]; |
449 | d = x*x+y*y+z*z; |
450 | if (fabs(adod-d) > DP_TOL0.001) |
451 | { |
452 | continue; |
453 | } |
454 | for (tl = 1; tl < tess-1; tl++) |
455 | { |
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); |
462 | |
463 | for (tl2 = 1; tl2 < tess-tl; tl2++) |
464 | { |
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); |
477 | if (tn >= n_dot) |
478 | { |
479 | ERROR__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 479, error("ico_dot: tn exceeds dimension of xus"); |
480 | } |
481 | divarc(xki, yki, zki, xji, yji, zji, tl2, tess-tl, |
482 | &x, &y, &z); |
483 | divarc(xkj, ykj, zkj, xij, yij, zij, tl, tess-tl2, |
484 | &x2, &y2, &z2); |
485 | divarc(xjk, yjk, zjk, xik, yik, zik, tl, tl+tl2, |
486 | &x3, &y3, &z3); |
487 | x = x+x2+x3; y = y+y2+y3; z = z+z2+z3; |
488 | d = sqrt(x*x+y*y+z*z); |
489 | xus[3*tn] = x/d; |
490 | xus[1+3*tn] = y/d; |
491 | xus[2+3*tn] = z/d; |
492 | tn++; |
493 | } /* cycle tl2 */ |
494 | } /* cycle tl */ |
495 | } /* cycle k */ |
496 | } /* cycle j */ |
497 | } /* cycle i */ |
498 | if (n_dot != tn) |
499 | { |
500 | ERROR__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 500, error("ico_dot: n_dot(%d) and tn(%d) differ", n_dot, tn); |
501 | } |
502 | } /* end of if (tess > 1) */ |
503 | |
504 | return n_dot; |
505 | } /* end of routine ico_dot_dod */ |
506 | |
507 | int unsp_type(int densit) |
508 | { |
509 | int i1, i2; |
510 | i1 = 1; |
511 | while (10*i1*i1+2 < densit) |
512 | { |
513 | i1++; |
514 | } |
515 | i2 = 1; |
516 | while (30*i2*i2+2 < densit) |
517 | { |
518 | i2++; |
519 | } |
520 | if (10*i1*i1-2 < 30*i2*i2-2) |
521 | { |
522 | return UNSP_ICO_ARC10; |
523 | } |
524 | else |
525 | { |
526 | return UNSP_ICO_DOD9; |
527 | } |
528 | } |
529 | |
530 | int make_unsp(int densit, int mode, int * num_dot, int cubus) |
531 | { |
532 | int *ico_wk, *ico_pt; |
533 | int ndot, ico_cube_cb, i, j, k, l, ijk, tn, tl, tl2; |
534 | real *xus; |
535 | int *work; |
536 | real x, y, z; |
537 | |
538 | if (xpunsp) |
539 | { |
540 | free(xpunsp); |
541 | } |
542 | |
543 | k = 1; if (mode < 0) |
544 | { |
545 | k = 0; mode = -mode; |
546 | } |
547 | if (mode == UNSP_ICO_ARC10) |
548 | { |
549 | ndot = ico_dot_arc(densit); |
550 | } |
551 | else if (mode == UNSP_ICO_DOD9) |
552 | { |
553 | ndot = ico_dot_dod(densit); |
554 | } |
555 | else |
556 | { |
557 | WARNING__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 557, warning2("make_unsp: mode %c%d not allowed", (k) ? '+' : '-', mode); |
558 | return 1; |
559 | } |
560 | |
561 | last_n_dot = ndot; last_densit = densit; last_unsp = mode; |
562 | *num_dot = ndot; if (k) |
563 | { |
564 | return 0; |
565 | } |
566 | |
567 | /* in the following the dots of the unit sphere may be resorted */ |
568 | last_unsp = -last_unsp; |
569 | |
570 | /* determine distribution of points in elementary cubes */ |
571 | if (cubus) |
572 | { |
573 | ico_cube = cubus; |
574 | } |
575 | else |
576 | { |
577 | last_cubus = 0; |
578 | i = 1; |
579 | while (i*i*i*2 < ndot) |
580 | { |
581 | i++; |
582 | } |
583 | ico_cube = max(i-1, 0)(((i-1) > (0)) ? (i-1) : (0) ); |
584 | } |
585 | ico_cube_cb = ico_cube*ico_cube*ico_cube; |
586 | del_cube = 2./((real)ico_cube); |
587 | snew(work, ndot)(work) = save_calloc("work", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 587, (ndot), sizeof(*(work))); |
588 | xus = xpunsp; |
589 | for (l = 0; l < ndot; l++) |
590 | { |
591 | i = max((int) floor((1.+xus[3*l])/del_cube), 0)((((int) floor((1.+xus[3*l])/del_cube)) > (0)) ? ((int) floor ((1.+xus[3*l])/del_cube)) : (0) ); |
592 | if (i >= ico_cube) |
593 | { |
594 | i = ico_cube-1; |
595 | } |
596 | j = max((int) floor((1.+xus[1+3*l])/del_cube), 0)((((int) floor((1.+xus[1+3*l])/del_cube)) > (0)) ? ((int) floor ((1.+xus[1+3*l])/del_cube)) : (0) ); |
597 | if (j >= ico_cube) |
598 | { |
599 | j = ico_cube-1; |
600 | } |
601 | k = max((int) floor((1.+xus[2+3*l])/del_cube), 0)((((int) floor((1.+xus[2+3*l])/del_cube)) > (0)) ? ((int) floor ((1.+xus[2+3*l])/del_cube)) : (0) ); |
602 | if (k >= ico_cube) |
603 | { |
604 | k = ico_cube-1; |
605 | } |
606 | ijk = i+j*ico_cube+k*ico_cube*ico_cube; |
607 | work[l] = ijk; |
608 | } |
609 | |
610 | snew(ico_wk, 2*ico_cube_cb+1)(ico_wk) = save_calloc("ico_wk", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 610, (2*ico_cube_cb+1), sizeof(*(ico_wk))); |
611 | |
612 | ico_pt = ico_wk+ico_cube_cb; |
613 | for (l = 0; l < ndot; l++) |
614 | { |
615 | ico_wk[work[l]]++; /* dots per elementary cube */ |
616 | } |
617 | |
618 | /* reordering of the coordinate array in accordance with box number */ |
619 | tn = 0; |
620 | for (i = 0; i < ico_cube; i++) |
621 | { |
622 | for (j = 0; j < ico_cube; j++) |
623 | { |
624 | for (k = 0; k < ico_cube; k++) |
625 | { |
626 | tl = 0; |
627 | tl2 = tn; |
628 | ijk = i+ico_cube*j+ico_cube*ico_cube*k; |
629 | *(ico_pt+ijk) = tn; |
630 | for (l = tl2; l < ndot; l++) |
631 | { |
632 | if (ijk == work[l]) |
633 | { |
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; |
639 | tn++; tl++; |
640 | } |
641 | } |
642 | *(ico_wk+ijk) = tl; |
643 | } /* cycle k */ |
644 | } /* cycle j */ |
645 | } /* cycle i */ |
646 | sfree(ico_wk)save_free("ico_wk", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 646, (ico_wk)); |
647 | sfree(work)save_free("work", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 647, (work)); |
648 | |
649 | return 0; |
650 | } |
651 | |
652 | |
653 | typedef struct _stwknb { |
654 | real x; |
655 | real y; |
656 | real z; |
657 | real dot; |
658 | } Neighb; |
659 | |
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, |
663 | real *value_of_vol, |
664 | real **lidots, int *nu_dots, |
665 | atom_id index[], int ePBC, matrix box) |
666 | { |
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; |
669 | int distribution; |
670 | int l; |
671 | int maxnei, nnei, last, maxdots = 0; |
672 | int *wkdot = NULL((void*)0), *wkbox = NULL((void*)0), *wkat1 = NULL((void*)0), *wkatm = NULL((void*)0); |
673 | Neighb *wknb, *ctnb; |
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((void*)0), *atom_area = NULL((void*)0); |
679 | int nxbox, nybox, nzbox, nxy, nxyz; |
680 | real xmin = 0, ymin = 0, zmin = 0, xmax, ymax, zmax, ra2max, d; |
681 | const real *pco; |
682 | /* Added DvdS 2006-07-19 */ |
683 | t_pbc pbc; |
684 | rvec ddx, *x = NULL((void*)0); |
685 | int iat_xx, jat_xx; |
686 | |
687 | distribution = unsp_type(densit); |
688 | if (distribution != -last_unsp || last_cubus != 4 || |
689 | (densit != last_densit && densit != last_n_dot)) |
690 | { |
691 | if (make_unsp(densit, (-distribution), &n_dot, 4)) |
692 | { |
693 | return 1; |
694 | } |
695 | } |
696 | xus = xpunsp; |
697 | |
698 | dotarea = FOURPI(4.*3.14159265358979323846)/(real) n_dot; |
699 | area = 0.; |
700 | |
701 | if (debug) |
702 | { |
703 | fprintf(debug, "nsc_dclm: n_dot=%5d %9.3f\n", n_dot, dotarea); |
704 | } |
705 | |
706 | /* start with neighbour list */ |
707 | /* calculate neighbour list with the box algorithm */ |
708 | if (nat == 0) |
709 | { |
710 | WARNING__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 710, warning2("nsc_dclm: no surface atoms selected"); |
711 | return 1; |
712 | } |
713 | if (mode & FLAG_VOLUME02) |
714 | { |
715 | vol = 0.; |
716 | } |
717 | if (mode & FLAG_DOTS01) |
718 | { |
719 | maxdots = (3*n_dot*nat)/10; |
720 | /* should be set to NULL on first user call */ |
721 | if (dots == NULL((void*)0)) |
722 | { |
723 | snew(dots, maxdots)(dots) = save_calloc("dots", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 723, (maxdots), sizeof(*(dots))); |
724 | } |
725 | else |
726 | { |
727 | srenew(dots, maxdots)(dots) = save_realloc("dots", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 727, (dots), (maxdots), sizeof(*(dots))); |
728 | } |
729 | |
730 | lfnr = 0; |
731 | } |
732 | if (mode & FLAG_ATOM_AREA04) |
733 | { |
734 | /* should be set to NULL on first user call */ |
735 | if (atom_area == NULL((void*)0)) |
736 | { |
737 | snew(atom_area, nat)(atom_area) = save_calloc("atom_area", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 737, (nat), sizeof(*(atom_area))); |
738 | } |
739 | else |
740 | { |
741 | srenew(atom_area, nat)(atom_area) = save_realloc("atom_area", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 741, (atom_area), (nat), sizeof(*(atom_area))); |
742 | } |
743 | } |
744 | |
745 | /* Compute minimum size for grid cells */ |
746 | ra2max = radius[index[0]]; |
747 | for (iat_xx = 1; (iat_xx < nat); iat_xx++) |
748 | { |
749 | iat = index[iat_xx]; |
750 | ra2max = max(ra2max, radius[iat])(((ra2max) > (radius[iat])) ? (ra2max) : (radius[iat]) ); |
751 | } |
752 | ra2max = 2*ra2max; |
753 | |
754 | /* Added DvdS 2006-07-19 */ |
755 | /* Updated 2008-10-09 */ |
756 | if (box) |
757 | { |
758 | set_pbc(&pbc, ePBC, box); |
759 | snew(x, nat)(x) = save_calloc("x", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 759, (nat), sizeof(*(x))); |
760 | for (i = 0; (i < nat); i++) |
761 | { |
762 | iat = index[0]; |
763 | copy_rvec(coords[iat], x[i]); |
764 | } |
765 | put_atoms_in_triclinic_unitcell(ecenterTRIC, box, nat, x); |
766 | nxbox = max(1, floor(norm(box[XX])/ra2max))(((1) > (floor(norm(box[0])/ra2max))) ? (1) : (floor(norm( box[0])/ra2max)) ); |
767 | nybox = max(1, floor(norm(box[YY])/ra2max))(((1) > (floor(norm(box[1])/ra2max))) ? (1) : (floor(norm( box[1])/ra2max)) ); |
768 | nzbox = max(1, floor(norm(box[ZZ])/ra2max))(((1) > (floor(norm(box[2])/ra2max))) ? (1) : (floor(norm( box[2])/ra2max)) ); |
769 | if (debug) |
770 | { |
771 | fprintf(debug, "nbox = %d, %d, %d\n", nxbox, nybox, nzbox); |
772 | } |
773 | } |
774 | else |
775 | { |
776 | /* dimensions of atomic set, cell edge is 2*ra_max */ |
777 | iat = index[0]; |
778 | xmin = coords[iat][XX0]; xmax = xmin; xs = xmin; |
779 | ymin = coords[iat][YY1]; ymax = ymin; ys = ymin; |
780 | zmin = coords[iat][ZZ2]; zmax = zmin; zs = zmin; |
781 | |
782 | for (iat_xx = 1; (iat_xx < nat); iat_xx++) |
783 | { |
784 | iat = index[iat_xx]; |
785 | pco = coords[iat]; |
786 | xmin = min(xmin, *pco)(((xmin) < (*pco)) ? (xmin) : (*pco) ); xmax = max(xmax, *pco)(((xmax) > (*pco)) ? (xmax) : (*pco) ); |
787 | ymin = min(ymin, *(pco+1))(((ymin) < (*(pco+1))) ? (ymin) : (*(pco+1)) ); ymax = max(ymax, *(pco+1))(((ymax) > (*(pco+1))) ? (ymax) : (*(pco+1)) ); |
788 | zmin = min(zmin, *(pco+2))(((zmin) < (*(pco+2))) ? (zmin) : (*(pco+2)) ); zmax = max(zmax, *(pco+2))(((zmax) > (*(pco+2))) ? (zmax) : (*(pco+2)) ); |
789 | xs = xs+ *pco; ys = ys+ *(pco+1); zs = zs+ *(pco+2); |
790 | } |
791 | xs = xs/ (real) nat; |
792 | ys = ys/ (real) nat; |
793 | zs = zs/ (real) nat; |
794 | if (debug) |
795 | { |
796 | fprintf(debug, "nsc_dclm: n_dot=%5d ra2max=%9.3f %9.3f\n", n_dot, ra2max, dotarea); |
797 | } |
798 | |
799 | d = xmax-xmin; nxbox = (int) max(ceil(d/ra2max), 1.)(((ceil(d/ra2max)) > (1.)) ? (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.)(((ceil(d/ra2max)) > (1.)) ? (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.)(((ceil(d/ra2max)) > (1.)) ? (ceil(d/ra2max)) : (1.) ); |
806 | d = (((real)nzbox)*ra2max-d)/2.; |
807 | zmin = zmin-d; zmax = zmax+d; |
Value stored to 'zmax' is never read | |
808 | } |
809 | /* Help variables */ |
810 | nxy = nxbox*nybox; |
811 | nxyz = nxy*nzbox; |
812 | |
813 | /* box number of atoms */ |
814 | snew(wkatm, nat)(wkatm) = save_calloc("wkatm", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 814, (nat), sizeof(*(wkatm))); |
815 | snew(wkat1, nat)(wkat1) = save_calloc("wkat1", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 815, (nat), sizeof(*(wkat1))); |
816 | snew(wkdot, n_dot)(wkdot) = save_calloc("wkdot", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 816, (n_dot), sizeof(*(wkdot))); |
817 | snew(wkbox, nxyz+1)(wkbox) = save_calloc("wkbox", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 817, (nxyz+1), sizeof(*(wkbox))); |
818 | |
819 | if (box) |
820 | { |
821 | matrix box_1; |
822 | rvec x_1; |
823 | int ix, iy, iz, m; |
824 | m_inv(box, box_1); |
825 | for (i = 0; (i < nat); i++) |
826 | { |
827 | mvmul(box_1, x[i], x_1); |
828 | ix = ((int)floor(x_1[XX0]*nxbox) + 2*nxbox) % nxbox; |
829 | iy = ((int)floor(x_1[YY1]*nybox) + 2*nybox) % nybox; |
830 | iz = ((int)floor(x_1[ZZ2]*nzbox) + 2*nzbox) % nzbox; |
831 | j = ix + iy*nxbox + iz*nxbox*nybox; |
832 | if (debug) |
833 | { |
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][XX0], x[i][YY1], x[i][ZZ2], x_1[XX0], x_1[YY1], x_1[ZZ2]); |
836 | } |
837 | range_check(j, 0, nxyz)_range_check(j, 0, nxyz, ((void*)0),"j", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 837); |
838 | wkat1[i] = j; |
839 | wkbox[j]++; |
840 | } |
841 | } |
842 | else |
843 | { |
844 | /* Put the atoms in their boxes */ |
845 | for (iat_xx = 0; (iat_xx < nat); iat_xx++) |
846 | { |
847 | iat = index[iat_xx]; |
848 | pco = coords[iat]; |
849 | i = (int) max(floor((pco[XX]-xmin)/ra2max), 0)(((floor((pco[0]-xmin)/ra2max)) > (0)) ? (floor((pco[0]-xmin )/ra2max)) : (0) ); |
850 | i = min(i, nxbox-1)(((i) < (nxbox-1)) ? (i) : (nxbox-1) ); |
851 | j = (int) max(floor((pco[YY]-ymin)/ra2max), 0)(((floor((pco[1]-ymin)/ra2max)) > (0)) ? (floor((pco[1]-ymin )/ra2max)) : (0) ); |
852 | j = min(j, nybox-1)(((j) < (nybox-1)) ? (j) : (nybox-1) ); |
853 | l = (int) max(floor((pco[ZZ]-zmin)/ra2max), 0)(((floor((pco[2]-zmin)/ra2max)) > (0)) ? (floor((pco[2]-zmin )/ra2max)) : (0) ); |
854 | l = min(l, nzbox-1)(((l) < (nzbox-1)) ? (l) : (nzbox-1) ); |
855 | i = i+j*nxbox+l*nxy; |
856 | wkat1[iat_xx] = i; |
857 | wkbox[i]++; |
858 | } |
859 | } |
860 | |
861 | /* sorting of atoms in accordance with box numbers */ |
862 | j = wkbox[0]; |
863 | for (i = 1; i < nxyz; i++) |
864 | { |
865 | j = max(wkbox[i], j)(((wkbox[i]) > (j)) ? (wkbox[i]) : (j) ); |
866 | } |
867 | for (i = 1; i <= nxyz; i++) |
868 | { |
869 | wkbox[i] += wkbox[i-1]; |
870 | } |
871 | |
872 | /* maxnei = (int) floor(ra2max*ra2max*ra2max*0.5); */ |
873 | maxnei = min(nat, 27*j)(((nat) < (27*j)) ? (nat) : (27*j) ); |
874 | snew(wknb, maxnei)(wknb) = save_calloc("wknb", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 874, (maxnei), sizeof(*(wknb))); |
875 | for (iat_xx = 0; iat_xx < nat; iat_xx++) |
876 | { |
877 | iat = index[iat_xx]; |
878 | range_check(wkat1[iat_xx], 0, nxyz)_range_check(wkat1[iat_xx], 0, nxyz, ((void*)0),"wkat1[iat_xx]" , "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 878); |
879 | wkatm[--wkbox[wkat1[iat_xx]]] = iat_xx; |
880 | if (debug) |
881 | { |
882 | fprintf(debug, "atom %5d on place %5d\n", iat, wkbox[wkat1[iat_xx]]); |
883 | } |
884 | } |
885 | |
886 | if (debug) |
887 | { |
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); |
892 | |
893 | for (i = 0; i < nxyz; i++) |
894 | { |
895 | fprintf(debug, "box %6d : atoms %4d-%4d %5d\n", |
896 | i, wkbox[i], wkbox[i+1]-1, wkbox[i+1]-wkbox[i]); |
897 | } |
898 | for (i = 0; i < nat; i++) |
899 | { |
900 | fprintf(debug, "list place %5d by atom %7d\n", i, index[wkatm[i]]); |
901 | } |
902 | } |
903 | |
904 | /* calculate surface for all atoms, step cube-wise */ |
905 | for (iz = 0; iz < nzbox; iz++) |
906 | { |
907 | iii = iz*nxy; |
908 | if (box) |
909 | { |
910 | izs = iz-1; |
911 | ize = min(iz+2, izs+nzbox)(((iz+2) < (izs+nzbox)) ? (iz+2) : (izs+nzbox) ); |
912 | } |
913 | else |
914 | { |
915 | izs = max(iz-1, 0)(((iz-1) > (0)) ? (iz-1) : (0) ); |
916 | ize = min(iz+2, nzbox)(((iz+2) < (nzbox)) ? (iz+2) : (nzbox) ); |
917 | } |
918 | for (iy = 0; iy < nybox; iy++) |
919 | { |
920 | ii = iy*nxbox+iii; |
921 | if (box) |
922 | { |
923 | iys = iy-1; |
924 | iye = min(iy+2, iys+nybox)(((iy+2) < (iys+nybox)) ? (iy+2) : (iys+nybox) ); |
925 | } |
926 | else |
927 | { |
928 | iys = max(iy-1, 0)(((iy-1) > (0)) ? (iy-1) : (0) ); |
929 | iye = min(iy+2, nybox)(((iy+2) < (nybox)) ? (iy+2) : (nybox) ); |
930 | } |
931 | for (ix = 0; ix < nxbox; ix++) |
932 | { |
933 | i = ii+ix; |
934 | iii1 = wkbox[i]; |
935 | iii2 = wkbox[i+1]; |
936 | if (iii1 >= iii2) |
937 | { |
938 | continue; |
939 | } |
940 | if (box) |
941 | { |
942 | ixs = ix-1; |
943 | ixe = min(ix+2, ixs+nxbox)(((ix+2) < (ixs+nxbox)) ? (ix+2) : (ixs+nxbox) ); |
944 | } |
945 | else |
946 | { |
947 | ixs = max(ix-1, 0)(((ix-1) > (0)) ? (ix-1) : (0) ); |
948 | ixe = min(ix+2, nxbox)(((ix+2) < (nxbox)) ? (ix+2) : (nxbox) ); |
949 | } |
950 | iiat = 0; |
951 | /* make intermediate atom list */ |
952 | for (jz = izs; jz < ize; jz++) |
953 | { |
954 | jjj = ((jz+nzbox) % nzbox)*nxy; |
955 | for (jy = iys; jy < iye; jy++) |
956 | { |
957 | jj = ((jy+nybox) % nybox)*nxbox+jjj; |
958 | for (jx = ixs; jx < ixe; jx++) |
959 | { |
960 | j = jj+((jx+nxbox) % nxbox); |
961 | for (jat = wkbox[j]; jat < wkbox[j+1]; jat++) |
962 | { |
963 | range_check(wkatm[jat], 0, nat)_range_check(wkatm[jat], 0, nat, ((void*)0),"wkatm[jat]", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 963); |
964 | range_check(iiat, 0, nat)_range_check(iiat, 0, nat, ((void*)0),"iiat", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 964); |
965 | wkat1[iiat] = wkatm[jat]; |
966 | iiat++; |
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++) |
972 | { |
973 | i_at = index[wkatm[iat]]; |
974 | ai = radius[i_at]; |
975 | aisq = ai*ai; |
976 | pco = coords[i_at]; |
977 | xi = pco[XX0]; yi = pco[YY1]; zi = pco[ZZ2]; |
978 | for (i = 0; i < n_dot; i++) |
979 | { |
980 | wkdot[i] = 0; |
981 | } |
982 | |
983 | ctnb = wknb; nnei = 0; |
984 | for (j = 0; j < iiat; j++) |
985 | { |
986 | j_at = index[wkat1[j]]; |
987 | if (j_at == i_at) |
988 | { |
989 | continue; |
990 | } |
991 | |
992 | aj = radius[j_at]; |
993 | ajsq = aj*aj; |
994 | pco = coords[j_at]; |
995 | |
996 | /* Added DvdS 2006-07-19 */ |
997 | if (box) |
998 | { |
999 | /*rvec xxi; |
1000 | |
1001 | xxi[XX] = xi; |
1002 | xxi[YY] = yi; |
1003 | xxi[ZZ] = zi; |
1004 | pbc_dx(&pbc,pco,xxi,ddx);*/ |
1005 | pbc_dx(&pbc, coords[j_at], coords[i_at], ddx); |
1006 | dx = ddx[XX0]; |
1007 | dy = ddx[YY1]; |
1008 | dz = ddx[ZZ2]; |
1009 | } |
1010 | else |
1011 | { |
1012 | dx = pco[XX0]-xi; |
1013 | dy = pco[YY1]-yi; |
1014 | dz = pco[ZZ2]-zi; |
1015 | } |
1016 | dd = dx*dx+dy*dy+dz*dz; |
1017 | as = ai+aj; |
1018 | if (dd > as*as) |
1019 | { |
1020 | continue; |
1021 | } |
1022 | nnei++; |
1023 | ctnb->x = dx; |
1024 | ctnb->y = dy; |
1025 | ctnb->z = dz; |
1026 | ctnb->dot = (dd+aisq-ajsq)/(2.*ai); /* reference dot product */ |
1027 | ctnb++; |
1028 | } |
1029 | |
1030 | /* check points on accessibility */ |
1031 | if (nnei) |
1032 | { |
1033 | last = 0; i_ac = 0; |
1034 | for (l = 0; l < n_dot; l++) |
1035 | { |
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) |
1039 | { |
1040 | for (j = 0; j < nnei; j++) |
1041 | { |
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) |
1044 | { |
1045 | last = j; |
1046 | break; |
1047 | } |
1048 | } |
1049 | if (j >= nnei) |
1050 | { |
1051 | i_ac++; |
1052 | wkdot[l] = 1; |
1053 | } |
1054 | } /* end of cycle j */ |
1055 | } /* end of cycle l */ |
1056 | } |
1057 | else |
1058 | { |
1059 | i_ac = n_dot; |
1060 | for (l = 0; l < n_dot; l++) |
1061 | { |
1062 | wkdot[l] = 1; |
1063 | } |
1064 | } |
1065 | |
1066 | if (debug) |
1067 | { |
1068 | fprintf(debug, "i_ac=%d, dotarea=%8.3f, aisq=%8.3f\n", |
1069 | i_ac, dotarea, aisq); |
1070 | } |
1071 | |
1072 | a = aisq*dotarea* (real) i_ac; |
1073 | area = area + a; |
1074 | if (mode & FLAG_ATOM_AREA04) |
1075 | { |
1076 | range_check(wkatm[iat], 0, nat)_range_check(wkatm[iat], 0, nat, ((void*)0),"wkatm[iat]", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 1076); |
1077 | atom_area[wkatm[iat]] = a; |
1078 | } |
1079 | if (mode & FLAG_DOTS01) |
1080 | { |
1081 | for (l = 0; l < n_dot; l++) |
1082 | { |
1083 | if (wkdot[l]) |
1084 | { |
1085 | lfnr++; |
1086 | if (maxdots <= 3*lfnr+1) |
1087 | { |
1088 | maxdots = maxdots+n_dot*3; |
1089 | srenew(dots, maxdots)(dots) = save_realloc("dots", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 1089, (dots), (maxdots), sizeof(*(dots))); |
1090 | } |
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; |
1094 | } |
1095 | } |
1096 | } |
1097 | if (mode & FLAG_VOLUME02) |
1098 | { |
1099 | dx = 0.; dy = 0.; dz = 0.; |
1100 | for (l = 0; l < n_dot; l++) |
1101 | { |
1102 | if (wkdot[l]) |
1103 | { |
1104 | dx = dx+xus[3*l]; |
1105 | dy = dy+xus[1+3*l]; |
1106 | dz = dz+xus[2+3*l]; |
1107 | } |
1108 | } |
1109 | vol = vol+aisq*(dx*(xi-xs)+dy*(yi-ys)+dz*(zi-zs)+ai* (real) i_ac); |
1110 | } |
1111 | |
1112 | } /* end of cycle "iat" */ |
1113 | } /* end of cycle "ix" */ |
1114 | } /* end of cycle "iy" */ |
1115 | } /* end of cycle "iz" */ |
1116 | |
1117 | sfree(wkatm)save_free("wkatm", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 1117, (wkatm)); |
1118 | sfree(wkat1)save_free("wkat1", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 1118, (wkat1)); |
1119 | sfree(wkdot)save_free("wkdot", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 1119, (wkdot)); |
1120 | sfree(wkbox)save_free("wkbox", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 1120, (wkbox)); |
1121 | sfree(wknb)save_free("wknb", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 1121, (wknb)); |
1122 | if (box) |
1123 | { |
1124 | sfree(x)save_free("x", "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , 1124, (x)); |
1125 | } |
1126 | if (mode & FLAG_VOLUME02) |
1127 | { |
1128 | vol = vol*FOURPI(4.*3.14159265358979323846)/(3.* (real) n_dot); |
1129 | *value_of_vol = vol; |
1130 | } |
1131 | if (mode & FLAG_DOTS01) |
1132 | { |
1133 | *nu_dots = lfnr; |
1134 | *lidots = dots; |
1135 | } |
1136 | if (mode & FLAG_ATOM_AREA04) |
1137 | { |
1138 | *at_area = atom_area; |
1139 | } |
1140 | *value_of_area = area; |
1141 | |
1142 | if (debug) |
1143 | { |
1144 | fprintf(debug, "area=%8.3f\n", area); |
1145 | } |
1146 | |
1147 | return 0; |
1148 | } |
1149 | |
1150 | |
1151 | #if TEST_NSC0 > 0 |
1152 | #define NAT 2 |
1153 | main () { |
1154 | |
1155 | int i, j, ndots; |
1156 | real co[3*NAT], ra[NAT], area, volume, a, b, c; |
1157 | real * dots; |
1158 | real * at_area; |
1159 | FILE *fp; |
1160 | |
1161 | |
1162 | a = 1.; c = 0.1; |
1163 | fp = fopen("nsc.txt", "w+"); |
1164 | for (i = 1; i <= NAT; i++) |
1165 | { |
1166 | j = i-1; |
1167 | co[3*i-3] = j*1*c; |
1168 | co[3*i-2] = j*1*c; |
1169 | co[3*i-1] = j*1*c; |
1170 | /* |
1171 | co[3*i-3] = i*1.4; |
1172 | co[3*i-2] = 0.; |
1173 | co[3*i-1] = 0.; |
1174 | */ |
1175 | /* |
1176 | co[3*i-2] = a*0.3; |
1177 | a = -a; b=0; |
1178 | if (i%3 == 0) b=0.5; |
1179 | co[3*i-1] = b; |
1180 | ra[i-1] = 2.0; |
1181 | */ |
1182 | ra[i-1] = (1.+j*0.5)*c; |
1183 | } |
1184 | /* |
1185 | if (NSC(co, ra, NAT, 42, NULL, &area, |
1186 | */ |
1187 | if (NSC(co, ra, NAT, 42, NULL((void*)0), &area, |
1188 | NULL((void*)0), NULL((void*)0), NULL((void*)0), NULL((void*)0))) |
1189 | { |
1190 | ERROR__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 1190, error("error in NSC"); |
1191 | } |
1192 | fprintf(fp, "\n"); |
1193 | fprintf(fp, "area : %8.3f\n", area); |
1194 | fprintf(fp, "\n"); |
1195 | fprintf(fp, "\n"); |
1196 | fprintf(fp, "next call\n"); |
1197 | fprintf(fp, "\n"); |
1198 | fprintf(fp, "\n"); |
1199 | |
1200 | if (NSC(co, ra, NAT, 42, FLAG_VOLUME02 | FLAG_ATOM_AREA04 | FLAG_DOTS01, &area, |
1201 | &at_area, &volume, |
1202 | &dots, &ndots)) |
1203 | { |
1204 | ERROR__file__ = "/home/alexxy/Develop/gromacs/src/gromacs/trajectoryanalysis/modules/nsc.c" , __line__ = 1204, error("error in NSC"); |
1205 | } |
1206 | |
1207 | fprintf(fp, "\n"); |
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); |
1214 | fprintf(fp, "\n"); |
1215 | for (i = 1; i <= NAT; i++) |
1216 | { |
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]); |
1219 | } |
1220 | fprintf(fp, "\n"); |
1221 | fprintf(fp, "DOTS : %8d\n", ndots); |
1222 | for (i = 1; i <= ndots; i++) |
1223 | { |
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]); |
1226 | } |
1227 | } |
1228 | #endif |