Bug Summary

File:gromacs/trajectoryanalysis/modules/nsc.c
Location:line 804, column 24
Description:Value stored to 'ymax' is never read

Annotated Source Code

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
62real *xpunsp = NULL((void*)0);
63real del_cube;
64int n_dot, ico_cube, last_n_dot = 0, last_densit = 0, last_unsp = 0;
65int 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
72const char * __file__; /* declared versions of macros */
73int __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
79void 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
92void 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
107real 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 */
124real rg, rh;
125
126void 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
146void 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
186int 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
322int 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
507int 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
530int 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
653typedef struct _stwknb {
654 real x;
655 real y;
656 real z;
657 real dot;
658} Neighb;
659
660int 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;
Value stored to 'ymax' is never read
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;
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
1153main () {
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