File: | gromacs/mdlib/nsgrid.c |
Location: | line 292, column 46 |
Description: | The left operand of '*' is a garbage value |
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-2004, 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 | /* This file is completely threadsafe - keep it that way! */ | |||
38 | #ifdef HAVE_CONFIG_H1 | |||
39 | #include <config.h> | |||
40 | #endif | |||
41 | ||||
42 | #include <stdlib.h> | |||
43 | ||||
44 | #include "typedefs.h" | |||
45 | #include "types/commrec.h" | |||
46 | #include "macros.h" | |||
47 | #include "gromacs/utility/smalloc.h" | |||
48 | #include "nsgrid.h" | |||
49 | #include "gromacs/utility/fatalerror.h" | |||
50 | #include "gromacs/math/vec.h" | |||
51 | #include "network.h" | |||
52 | #include "domdec.h" | |||
53 | #include "pbc.h" | |||
54 | #include <stdio.h> | |||
55 | #include "gromacs/utility/futil.h" | |||
56 | #include "gromacs/fileio/pdbio.h" | |||
57 | ||||
58 | /*********************************** | |||
59 | * Grid Routines | |||
60 | ***********************************/ | |||
61 | ||||
62 | const char *range_warn = | |||
63 | "Explanation: During neighborsearching, we assign each particle to a grid\n" | |||
64 | "based on its coordinates. If your system contains collisions or parameter\n" | |||
65 | "errors that give particles very high velocities you might end up with some\n" | |||
66 | "coordinates being +-Infinity or NaN (not-a-number). Obviously, we cannot\n" | |||
67 | "put these on a grid, so this is usually where we detect those errors.\n" | |||
68 | "Make sure your system is properly energy-minimized and that the potential\n" | |||
69 | "energy seems reasonable before trying again."; | |||
70 | ||||
71 | static void calc_x_av_stddev(int n, rvec *x, rvec av, rvec stddev) | |||
72 | { | |||
73 | dvec s1, s2; | |||
74 | int i, d; | |||
75 | ||||
76 | clear_dvec(s1); | |||
77 | clear_dvec(s2); | |||
78 | ||||
79 | for (i = 0; i < n; i++) | |||
80 | { | |||
81 | for (d = 0; d < DIM3; d++) | |||
82 | { | |||
83 | s1[d] += x[i][d]; | |||
84 | s2[d] += x[i][d]*x[i][d]; | |||
85 | } | |||
86 | } | |||
87 | ||||
88 | dsvmul(1.0/n, s1, s1); | |||
89 | dsvmul(1.0/n, s2, s2); | |||
90 | ||||
91 | for (d = 0; d < DIM3; d++) | |||
92 | { | |||
93 | av[d] = s1[d]; | |||
94 | stddev[d] = sqrt(s2[d] - s1[d]*s1[d]); | |||
95 | } | |||
96 | } | |||
97 | ||||
98 | void get_nsgrid_boundaries_vac(real av, real stddev, | |||
99 | real *bound0, real *bound1, | |||
100 | real *bdens0, real *bdens1) | |||
101 | { | |||
102 | /* Set the grid to 2 times the standard deviation of | |||
103 | * the charge group centers in both directions. | |||
104 | * For a uniform bounded distribution the width is sqrt(3)*stddev, | |||
105 | * so all charge groups fall within the width. | |||
106 | * For a sphere stddev is r/sqrt(5): 99.2% falls within the width. | |||
107 | * For a Gaussian distribution 98% fall within the width. | |||
108 | */ | |||
109 | *bound0 = av - NSGRID_STDDEV_FAC2.0*stddev; | |||
110 | *bound1 = av + NSGRID_STDDEV_FAC2.0*stddev; | |||
111 | ||||
112 | *bdens0 = av - GRID_STDDEV_FACsqrt(3)*stddev; | |||
113 | *bdens1 = av + GRID_STDDEV_FACsqrt(3)*stddev; | |||
114 | } | |||
115 | ||||
116 | static void dd_box_bounds_to_ns_bounds(real box0, real box_size, | |||
117 | real *gr0, real *gr1) | |||
118 | { | |||
119 | real av, stddev; | |||
120 | ||||
121 | /* Redetermine av and stddev from the DD box boundaries */ | |||
122 | av = box0 + 0.5*box_size; | |||
123 | stddev = 0.5*box_size/GRID_STDDEV_FACsqrt(3); | |||
124 | ||||
125 | *gr0 = av - NSGRID_STDDEV_FAC2.0*stddev; | |||
126 | *gr1 = av + NSGRID_STDDEV_FAC2.0*stddev; | |||
127 | } | |||
128 | ||||
129 | void get_nsgrid_boundaries(int nboundeddim, matrix box, | |||
130 | gmx_domdec_t *dd, | |||
131 | gmx_ddbox_t *ddbox, rvec *gr0, rvec *gr1, | |||
132 | int ncg, rvec *cgcm, | |||
133 | rvec grid_x0, rvec grid_x1, | |||
134 | real *grid_density) | |||
135 | { | |||
136 | rvec av, stddev; | |||
137 | real vol, bdens0, bdens1; | |||
138 | int d; | |||
139 | ||||
140 | if (nboundeddim < DIM3) | |||
141 | { | |||
142 | calc_x_av_stddev(ncg, cgcm, av, stddev); | |||
143 | } | |||
144 | ||||
145 | vol = 1; | |||
146 | for (d = 0; d < DIM3; d++) | |||
147 | { | |||
148 | if (d < nboundeddim) | |||
149 | { | |||
150 | grid_x0[d] = (gr0 != NULL((void*)0) ? (*gr0)[d] : 0); | |||
151 | grid_x1[d] = (gr1 != NULL((void*)0) ? (*gr1)[d] : box[d][d]); | |||
152 | vol *= (grid_x1[d] - grid_x0[d]); | |||
153 | } | |||
154 | else | |||
155 | { | |||
156 | if (ddbox == NULL((void*)0)) | |||
157 | { | |||
158 | get_nsgrid_boundaries_vac(av[d], stddev[d], | |||
159 | &grid_x0[d], &grid_x1[d], | |||
160 | &bdens0, &bdens1); | |||
161 | } | |||
162 | else | |||
163 | { | |||
164 | /* Temporary fix which uses global ddbox boundaries | |||
165 | * for unbounded dimensions. | |||
166 | * Should be replaced by local boundaries, which makes | |||
167 | * the ns grid smaller and does not require global comm. | |||
168 | */ | |||
169 | dd_box_bounds_to_ns_bounds(ddbox->box0[d], ddbox->box_size[d], | |||
170 | &grid_x0[d], &grid_x1[d]); | |||
171 | bdens0 = grid_x0[d]; | |||
172 | bdens1 = grid_x1[d]; | |||
173 | } | |||
174 | /* Check for a DD cell not at a lower edge */ | |||
175 | if (dd != NULL((void*)0) && gr0 != NULL((void*)0) && dd->ci[d] > 0) | |||
176 | { | |||
177 | grid_x0[d] = (*gr0)[d]; | |||
178 | bdens0 = (*gr0)[d]; | |||
179 | } | |||
180 | /* Check for a DD cell not at a higher edge */ | |||
181 | if (dd != NULL((void*)0) && gr1 != NULL((void*)0) && dd->ci[d] < dd->nc[d]-1) | |||
182 | { | |||
183 | grid_x1[d] = (*gr1)[d]; | |||
184 | bdens1 = (*gr1)[d]; | |||
185 | } | |||
186 | vol *= (bdens1 - bdens0); | |||
187 | } | |||
188 | ||||
189 | if (debug) | |||
190 | { | |||
191 | fprintf(debug, "Set grid boundaries dim %d: %f %f\n", | |||
192 | d, grid_x0[d], grid_x1[d]); | |||
193 | } | |||
194 | } | |||
195 | ||||
196 | *grid_density = ncg/vol; | |||
197 | } | |||
198 | ||||
199 | static void set_grid_sizes(matrix box, rvec izones_x0, rvec izones_x1, real rlist, | |||
200 | const gmx_domdec_t *dd, const gmx_ddbox_t *ddbox, | |||
201 | t_grid *grid, | |||
202 | real grid_density) | |||
203 | { | |||
204 | int i, j; | |||
205 | gmx_bool bDD, bDDRect; | |||
206 | rvec av, stddev; | |||
207 | rvec izones_size; | |||
208 | real inv_r_ideal, size, add_tric, radd; | |||
209 | ||||
210 | for (i = 0; (i < DIM3); i++) | |||
| ||||
211 | { | |||
212 | if (debug) | |||
213 | { | |||
214 | fprintf(debug, | |||
215 | "set_grid_sizes, i-zone bounds for dim %d: %6.3f %6.3f\n", | |||
216 | i, izones_x0[i], izones_x1[i]); | |||
217 | } | |||
218 | izones_size[i] = izones_x1[i] - izones_x0[i]; | |||
219 | } | |||
220 | ||||
221 | /* Use the ideal number of cg's per cell to set the ideal cell size */ | |||
222 | inv_r_ideal = pow(grid_density/grid->ncg_ideal, 1.0/3.0); | |||
223 | if (rlist > 0 && inv_r_ideal*rlist < 1) | |||
224 | { | |||
225 | inv_r_ideal = 1/rlist; | |||
226 | } | |||
227 | if (debug) | |||
228 | { | |||
229 | fprintf(debug, "CG density %f ideal ns cell size %f\n", | |||
230 | grid_density, 1/inv_r_ideal); | |||
231 | } | |||
232 | ||||
233 | clear_rvec(grid->cell_offset); | |||
234 | for (i = 0; (i < DIM3); i++) | |||
235 | { | |||
236 | /* Initial settings, for DD might change below */ | |||
237 | grid->cell_offset[i] = izones_x0[i]; | |||
238 | size = izones_size[i]; | |||
239 | ||||
240 | bDD = dd && (dd->nc[i] > 1); | |||
241 | if (!bDD) | |||
242 | { | |||
243 | bDDRect = FALSE0; | |||
244 | } | |||
245 | else | |||
246 | { | |||
247 | /* With DD grid cell jumps only the first decomposition | |||
248 | * direction has uniform DD cell boundaries. | |||
249 | */ | |||
250 | bDDRect = !(ddbox->tric_dir[i] || | |||
251 | (dd->bGridJump && i != dd->dim[0])); | |||
252 | ||||
253 | radd = rlist; | |||
254 | if (i >= ddbox->npbcdim && | |||
255 | (rlist == 0 || | |||
256 | izones_x1[i] + radd > ddbox->box0[i] + ddbox->box_size[i])) | |||
257 | { | |||
258 | radd = ddbox->box0[i] + ddbox->box_size[i] - izones_x1[i]; | |||
259 | if (radd < 0) | |||
260 | { | |||
261 | radd = 0; | |||
262 | } | |||
263 | } | |||
264 | ||||
265 | /* With DD we only need a grid of one DD cell size + rlist */ | |||
266 | if (bDDRect) | |||
267 | { | |||
268 | size += radd; | |||
269 | } | |||
270 | else | |||
271 | { | |||
272 | size += radd/ddbox->skew_fac[i]; | |||
273 | } | |||
274 | ||||
275 | /* Check if the cell boundary in this direction is | |||
276 | * perpendicular to the Cartesian axis. | |||
277 | */ | |||
278 | for (j = i+1; j < grid->npbcdim; j++) | |||
279 | { | |||
280 | if (box[j][i] != 0) | |||
281 | { | |||
282 | /* Correct the offset for the home cell location */ | |||
283 | grid->cell_offset[i] += izones_x0[j]*box[j][i]/box[j][j]; | |||
284 | ||||
285 | /* Correct the offset and size for the off-diagonal | |||
286 | * displacement of opposing DD cell corners. | |||
287 | */ | |||
288 | /* Without rouding we would need to add: | |||
289 | * box[j][i]*rlist/(dd->skew_fac[i]*box[j][j]) | |||
290 | */ | |||
291 | /* Determine the shift for the corners of the triclinic box */ | |||
292 | add_tric = izones_size[j]*box[j][i]/box[j][j]; | |||
| ||||
293 | if (dd && dd->ndim == 1 && j == ZZ2) | |||
294 | { | |||
295 | /* With 1D domain decomposition the cg's are not in | |||
296 | * the triclinic box, but trilinic x-y and rectangular y-z. | |||
297 | * Therefore we need to add the shift from the trilinic | |||
298 | * corner to the corner at y=0. | |||
299 | */ | |||
300 | add_tric += -box[YY1][XX0]*box[ZZ2][YY1]/box[YY1][YY1]; | |||
301 | } | |||
302 | if (box[j][i] < 0) | |||
303 | { | |||
304 | grid->cell_offset[i] += add_tric; | |||
305 | size -= add_tric; | |||
306 | } | |||
307 | else | |||
308 | { | |||
309 | size += add_tric; | |||
310 | } | |||
311 | } | |||
312 | } | |||
313 | } | |||
314 | if (!bDDRect) | |||
315 | { | |||
316 | /* No DD or the box is triclinic is this direction: | |||
317 | * we will use the normal grid ns that checks all cells | |||
318 | * that are within cut-off distance of the i-particle. | |||
319 | */ | |||
320 | grid->n[i] = (int)(size*inv_r_ideal + 0.5); | |||
321 | if (grid->n[i] < 2) | |||
322 | { | |||
323 | grid->n[i] = 2; | |||
324 | } | |||
325 | grid->cell_size[i] = size/grid->n[i]; | |||
326 | grid->ncpddc[i] = 0; | |||
327 | } | |||
328 | else | |||
329 | { | |||
330 | /* We use grid->ncpddc[i] such that all particles | |||
331 | * in one ns cell belong to a single DD cell only. | |||
332 | * We can then beforehand exclude certain ns grid cells | |||
333 | * for non-home i-particles. | |||
334 | */ | |||
335 | grid->ncpddc[i] = (int)(izones_size[i]*inv_r_ideal + 0.5); | |||
336 | if (grid->ncpddc[i] < 2) | |||
337 | { | |||
338 | grid->ncpddc[i] = 2; | |||
339 | } | |||
340 | grid->cell_size[i] = izones_size[i]/grid->ncpddc[i]; | |||
341 | grid->n[i] = grid->ncpddc[i] + (int)(radd/grid->cell_size[i]) + 1; | |||
342 | } | |||
343 | if (debug) | |||
344 | { | |||
345 | fprintf(debug, "grid dim %d size %d x %f: %f - %f\n", | |||
346 | i, grid->n[i], grid->cell_size[i], | |||
347 | grid->cell_offset[i], | |||
348 | grid->cell_offset[i]+grid->n[i]*grid->cell_size[i]); | |||
349 | } | |||
350 | } | |||
351 | ||||
352 | if (debug) | |||
353 | { | |||
354 | fprintf(debug, "CG ncg ideal %d, actual density %.1f\n", | |||
355 | grid->ncg_ideal, grid_density*grid->cell_size[XX0]*grid->cell_size[YY1]*grid->cell_size[ZZ2]); | |||
356 | } | |||
357 | } | |||
358 | ||||
359 | t_grid *init_grid(FILE *fplog, t_forcerec *fr) | |||
360 | { | |||
361 | int d, m; | |||
362 | char *ptr; | |||
363 | t_grid *grid; | |||
364 | ||||
365 | snew(grid, 1)(grid) = save_calloc("grid", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 365, (1), sizeof(*(grid))); | |||
366 | ||||
367 | grid->npbcdim = ePBC2npbcdim(fr->ePBC); | |||
368 | ||||
369 | if (fr->ePBC == epbcXY && fr->nwall == 2) | |||
370 | { | |||
371 | grid->nboundeddim = 3; | |||
372 | } | |||
373 | else | |||
374 | { | |||
375 | grid->nboundeddim = grid->npbcdim; | |||
376 | } | |||
377 | ||||
378 | if (debug) | |||
379 | { | |||
380 | fprintf(debug, "The coordinates are bounded in %d dimensions\n", | |||
381 | grid->nboundeddim); | |||
382 | } | |||
383 | ||||
384 | /* The ideal number of cg's per ns grid cell seems to be 10 */ | |||
385 | grid->ncg_ideal = 10; | |||
386 | ptr = getenv("GMX_NSCELL_NCG"); | |||
387 | if (ptr) | |||
388 | { | |||
389 | sscanf(ptr, "%d", &grid->ncg_ideal); | |||
390 | if (fplog) | |||
391 | { | |||
392 | fprintf(fplog, "Set ncg_ideal to %d\n", grid->ncg_ideal); | |||
393 | } | |||
394 | if (grid->ncg_ideal <= 0) | |||
395 | { | |||
396 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c", 396, "The number of cg's per cell should be > 0"); | |||
397 | } | |||
398 | } | |||
399 | if (debug) | |||
400 | { | |||
401 | fprintf(debug, "Set ncg_ideal to %d\n", grid->ncg_ideal); | |||
402 | } | |||
403 | ||||
404 | return grid; | |||
405 | } | |||
406 | ||||
407 | void done_grid(t_grid *grid) | |||
408 | { | |||
409 | grid->nr = 0; | |||
410 | clear_ivec(grid->n); | |||
411 | grid->ncells = 0; | |||
412 | sfree(grid->cell_index)save_free("grid->cell_index", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 412, (grid->cell_index)); | |||
413 | sfree(grid->a)save_free("grid->a", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 413, (grid->a)); | |||
414 | sfree(grid->index)save_free("grid->index", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 414, (grid->index)); | |||
415 | sfree(grid->nra)save_free("grid->nra", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 415, (grid->nra)); | |||
416 | grid->cells_nalloc = 0; | |||
417 | sfree(grid->dcx2)save_free("grid->dcx2", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 417, (grid->dcx2)); | |||
418 | sfree(grid->dcy2)save_free("grid->dcy2", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 418, (grid->dcy2)); | |||
419 | sfree(grid->dcz2)save_free("grid->dcz2", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 419, (grid->dcz2)); | |||
420 | grid->dc_nalloc = 0; | |||
421 | ||||
422 | if (debug) | |||
423 | { | |||
424 | fprintf(debug, "Successfully freed memory for grid pointers."); | |||
425 | } | |||
426 | } | |||
427 | ||||
428 | int xyz2ci_(int nry, int nrz, int x, int y, int z) | |||
429 | /* Return the cell index */ | |||
430 | { | |||
431 | return (nry*nrz*x+nrz*y+z); | |||
432 | } | |||
433 | ||||
434 | void ci2xyz(t_grid *grid, int i, int *x, int *y, int *z) | |||
435 | /* Return x,y and z from the cell index */ | |||
436 | { | |||
437 | int ci; | |||
438 | ||||
439 | range_check_mesg(i, 0, grid->nr, range_warn)_range_check(i, 0, grid->nr, range_warn,"i", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 439); | |||
440 | ||||
441 | ci = grid->cell_index[i]; | |||
442 | *x = ci / (grid->n[YY1]*grid->n[ZZ2]); | |||
443 | ci -= (*x)*grid->n[YY1]*grid->n[ZZ2]; | |||
444 | *y = ci / grid->n[ZZ2]; | |||
445 | ci -= (*y)*grid->n[ZZ2]; | |||
446 | *z = ci; | |||
447 | } | |||
448 | ||||
449 | static int ci_not_used(ivec n) | |||
450 | { | |||
451 | /* Return an improbable value */ | |||
452 | return xyz2ci(n[YY], n[ZZ], 3*n[XX], 3*n[YY], 3*n[ZZ])((n[1])*(n[2])*(3*n[0])+(n[2])*(3*n[1])+(3*n[2])); | |||
453 | } | |||
454 | ||||
455 | static void set_grid_ncg(t_grid *grid, int ncg) | |||
456 | { | |||
457 | int nr_old, i; | |||
458 | ||||
459 | grid->nr = ncg; | |||
460 | if (grid->nr+1 > grid->nr_alloc) | |||
461 | { | |||
462 | nr_old = grid->nr_alloc; | |||
463 | grid->nr_alloc = over_alloc_dd(grid->nr) + 1; | |||
464 | srenew(grid->cell_index, grid->nr_alloc)(grid->cell_index) = save_realloc("grid->cell_index", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 464, (grid->cell_index), (grid->nr_alloc), sizeof(*(grid ->cell_index))); | |||
465 | for (i = nr_old; i < grid->nr_alloc; i++) | |||
466 | { | |||
467 | grid->cell_index[i] = 0; | |||
468 | } | |||
469 | srenew(grid->a, grid->nr_alloc)(grid->a) = save_realloc("grid->a", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 469, (grid->a), (grid->nr_alloc), sizeof(*(grid->a ))); | |||
470 | } | |||
471 | } | |||
472 | ||||
473 | void grid_first(FILE *fplog, t_grid *grid, | |||
474 | gmx_domdec_t *dd, const gmx_ddbox_t *ddbox, | |||
475 | matrix box, rvec izones_x0, rvec izones_x1, | |||
476 | real rlistlong, real grid_density) | |||
477 | { | |||
478 | int i, m; | |||
479 | ivec cx; | |||
480 | ||||
481 | set_grid_sizes(box, izones_x0, izones_x1, rlistlong, dd, ddbox, grid, grid_density); | |||
482 | ||||
483 | grid->ncells = grid->n[XX0]*grid->n[YY1]*grid->n[ZZ2]; | |||
484 | ||||
485 | if (grid->ncells+1 > grid->cells_nalloc) | |||
486 | { | |||
487 | /* Allocate double the size so we have some headroom */ | |||
488 | grid->cells_nalloc = 2*grid->ncells; | |||
489 | srenew(grid->nra, grid->cells_nalloc+1)(grid->nra) = save_realloc("grid->nra", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 489, (grid->nra), (grid->cells_nalloc+1), sizeof(*(grid ->nra))); | |||
490 | srenew(grid->index, grid->cells_nalloc+1)(grid->index) = save_realloc("grid->index", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 490, (grid->index), (grid->cells_nalloc+1), sizeof(*( grid->index))); | |||
491 | ||||
492 | if (fplog) | |||
493 | { | |||
494 | fprintf(fplog, "Grid: %d x %d x %d cells\n", | |||
495 | grid->n[XX0], grid->n[YY1], grid->n[ZZ2]); | |||
496 | } | |||
497 | } | |||
498 | ||||
499 | m = max(grid->n[XX], max(grid->n[YY], grid->n[ZZ]))(((grid->n[0]) > ((((grid->n[1]) > (grid->n[2] )) ? (grid->n[1]) : (grid->n[2]) ))) ? (grid->n[0]) : ((((grid->n[1]) > (grid->n[2])) ? (grid->n[1]) : (grid->n[2]) )) ); | |||
500 | if (m > grid->dc_nalloc) | |||
501 | { | |||
502 | /* Allocate with double the initial size for box scaling */ | |||
503 | grid->dc_nalloc = 2*m; | |||
504 | srenew(grid->dcx2, grid->dc_nalloc)(grid->dcx2) = save_realloc("grid->dcx2", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 504, (grid->dcx2), (grid->dc_nalloc), sizeof(*(grid-> dcx2))); | |||
505 | srenew(grid->dcy2, grid->dc_nalloc)(grid->dcy2) = save_realloc("grid->dcy2", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 505, (grid->dcy2), (grid->dc_nalloc), sizeof(*(grid-> dcy2))); | |||
506 | srenew(grid->dcz2, grid->dc_nalloc)(grid->dcz2) = save_realloc("grid->dcz2", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 506, (grid->dcz2), (grid->dc_nalloc), sizeof(*(grid-> dcz2))); | |||
507 | } | |||
508 | ||||
509 | grid->nr = 0; | |||
510 | for (i = 0; (i < grid->ncells); i++) | |||
511 | { | |||
512 | grid->nra[i] = 0; | |||
513 | } | |||
514 | } | |||
515 | ||||
516 | static void calc_bor(int cg0, int cg1, int ncg, int CG0[2], int CG1[2]) | |||
517 | { | |||
518 | if (cg1 > ncg) | |||
519 | { | |||
520 | CG0[0] = cg0; | |||
521 | CG1[0] = ncg; | |||
522 | CG0[1] = 0; | |||
523 | CG1[1] = cg1-ncg; | |||
524 | } | |||
525 | else | |||
526 | { | |||
527 | CG0[0] = cg0; | |||
528 | CG1[0] = cg1; | |||
529 | CG0[1] = 0; | |||
530 | CG1[1] = 0; | |||
531 | } | |||
532 | if (debug) | |||
533 | { | |||
534 | int m; | |||
535 | ||||
536 | fprintf(debug, "calc_bor: cg0=%d, cg1=%d, ncg=%d\n", cg0, cg1, ncg); | |||
537 | for (m = 0; (m < 2); m++) | |||
538 | { | |||
539 | fprintf(debug, "CG0[%d]=%d, CG1[%d]=%d\n", m, CG0[m], m, CG1[m]); | |||
540 | } | |||
541 | } | |||
542 | ||||
543 | } | |||
544 | ||||
545 | void calc_elemnr(t_grid *grid, int cg0, int cg1, int ncg) | |||
546 | { | |||
547 | int CG0[2], CG1[2]; | |||
548 | int *cell_index = grid->cell_index; | |||
549 | int *nra = grid->nra; | |||
550 | int i, m, ncells; | |||
551 | int ci, not_used; | |||
552 | ||||
553 | ncells = grid->ncells; | |||
554 | if (ncells <= 0) | |||
555 | { | |||
556 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c", 556, "Number of grid cells is zero. Probably the system and box collapsed.\n"); | |||
557 | } | |||
558 | ||||
559 | not_used = ci_not_used(grid->n); | |||
560 | ||||
561 | calc_bor(cg0, cg1, ncg, CG0, CG1); | |||
562 | for (m = 0; (m < 2); m++) | |||
563 | { | |||
564 | for (i = CG0[m]; (i < CG1[m]); i++) | |||
565 | { | |||
566 | ci = cell_index[i]; | |||
567 | if (ci != not_used) | |||
568 | { | |||
569 | range_check_mesg(ci, 0, ncells, range_warn)_range_check(ci, 0, ncells, range_warn,"ci", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 569); | |||
570 | nra[ci]++; | |||
571 | } | |||
572 | } | |||
573 | } | |||
574 | } | |||
575 | ||||
576 | void calc_ptrs(t_grid *grid) | |||
577 | { | |||
578 | int *index = grid->index; | |||
579 | int *nra = grid->nra; | |||
580 | int ix, iy, iz, ci, nr; | |||
581 | int nnra, ncells; | |||
582 | ||||
583 | ncells = grid->ncells; | |||
584 | if (ncells <= 0) | |||
585 | { | |||
586 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c", 586, "Number of grid cells is zero. Probably the system and box collapsed.\n"); | |||
587 | } | |||
588 | ||||
589 | ci = nr = 0; | |||
590 | for (ix = 0; (ix < grid->n[XX0]); ix++) | |||
591 | { | |||
592 | for (iy = 0; (iy < grid->n[YY1]); iy++) | |||
593 | { | |||
594 | for (iz = 0; (iz < grid->n[ZZ2]); iz++, ci++) | |||
595 | { | |||
596 | range_check_mesg(ci, 0, ncells, range_warn)_range_check(ci, 0, ncells, range_warn,"ci", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 596); | |||
597 | index[ci] = nr; | |||
598 | nnra = nra[ci]; | |||
599 | nr += nnra; | |||
600 | nra[ci] = 0; | |||
601 | } | |||
602 | } | |||
603 | } | |||
604 | } | |||
605 | ||||
606 | void grid_last(t_grid *grid, int cg0, int cg1, int ncg) | |||
607 | { | |||
608 | int CG0[2], CG1[2]; | |||
609 | int i, m; | |||
610 | int ci, not_used, ind, ncells; | |||
611 | int *cell_index = grid->cell_index; | |||
612 | int *nra = grid->nra; | |||
613 | int *index = grid->index; | |||
614 | int *a = grid->a; | |||
615 | ||||
616 | ncells = grid->ncells; | |||
617 | if (ncells <= 0) | |||
618 | { | |||
619 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c", 619, "Number of grid cells is zero. Probably the system and box collapsed.\n"); | |||
620 | } | |||
621 | ||||
622 | not_used = ci_not_used(grid->n); | |||
623 | ||||
624 | calc_bor(cg0, cg1, ncg, CG0, CG1); | |||
625 | for (m = 0; (m < 2); m++) | |||
626 | { | |||
627 | for (i = CG0[m]; (i < CG1[m]); i++) | |||
628 | { | |||
629 | ci = cell_index[i]; | |||
630 | if (ci != not_used) | |||
631 | { | |||
632 | range_check_mesg(ci, 0, ncells, range_warn)_range_check(ci, 0, ncells, range_warn,"ci", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 632); | |||
633 | ind = index[ci]+nra[ci]++; | |||
634 | range_check_mesg(ind, 0, grid->nr, range_warn)_range_check(ind, 0, grid->nr, range_warn,"ind", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 634); | |||
635 | a[ind] = i; | |||
636 | } | |||
637 | } | |||
638 | } | |||
639 | } | |||
640 | ||||
641 | void fill_grid(gmx_domdec_zones_t *dd_zones, | |||
642 | t_grid *grid, int ncg_tot, | |||
643 | int cg0, int cg1, rvec cg_cm[]) | |||
644 | { | |||
645 | int *cell_index; | |||
646 | int nrx, nry, nrz; | |||
647 | rvec n_box, offset; | |||
648 | int zone, ccg0, ccg1, cg, d, not_used; | |||
649 | ivec shift0, useall, b0, b1, ind; | |||
650 | gmx_bool bUse; | |||
651 | ||||
652 | if (cg0 == -1) | |||
653 | { | |||
654 | /* We have already filled the grid up to grid->ncg, | |||
655 | * continue from there. | |||
656 | */ | |||
657 | cg0 = grid->nr; | |||
658 | } | |||
659 | ||||
660 | set_grid_ncg(grid, ncg_tot); | |||
661 | ||||
662 | cell_index = grid->cell_index; | |||
663 | ||||
664 | /* Initiate cell borders */ | |||
665 | nrx = grid->n[XX0]; | |||
666 | nry = grid->n[YY1]; | |||
667 | nrz = grid->n[ZZ2]; | |||
668 | for (d = 0; d < DIM3; d++) | |||
669 | { | |||
670 | if (grid->cell_size[d] > 0) | |||
671 | { | |||
672 | n_box[d] = 1/grid->cell_size[d]; | |||
673 | } | |||
674 | else | |||
675 | { | |||
676 | n_box[d] = 0; | |||
677 | } | |||
678 | } | |||
679 | copy_rvec(grid->cell_offset, offset); | |||
680 | ||||
681 | if (debug) | |||
682 | { | |||
683 | fprintf(debug, "Filling grid from %d to %d\n", cg0, cg1); | |||
684 | } | |||
685 | ||||
686 | debug_gmx(); | |||
687 | if (dd_zones == NULL((void*)0)) | |||
688 | { | |||
689 | for (cg = cg0; cg < cg1; cg++) | |||
690 | { | |||
691 | for (d = 0; d < DIM3; d++) | |||
692 | { | |||
693 | ind[d] = (cg_cm[cg][d] - offset[d])*n_box[d]; | |||
694 | /* With pbc we should be done here. | |||
695 | * Without pbc cg's outside the grid | |||
696 | * should be assigned to the closest grid cell. | |||
697 | */ | |||
698 | if (ind[d] < 0) | |||
699 | { | |||
700 | ind[d] = 0; | |||
701 | } | |||
702 | else if (ind[d] >= grid->n[d]) | |||
703 | { | |||
704 | ind[d] = grid->n[d] - 1; | |||
705 | } | |||
706 | } | |||
707 | cell_index[cg] = xyz2ci(nry, nrz, ind[XX], ind[YY], ind[ZZ])((nry)*(nrz)*(ind[0])+(nrz)*(ind[1])+(ind[2])); | |||
708 | } | |||
709 | } | |||
710 | else | |||
711 | { | |||
712 | for (zone = 0; zone < dd_zones->n; zone++) | |||
713 | { | |||
714 | ccg0 = dd_zones->cg_range[zone]; | |||
715 | ccg1 = dd_zones->cg_range[zone+1]; | |||
716 | if (ccg1 <= cg0 || ccg0 >= cg1) | |||
717 | { | |||
718 | continue; | |||
719 | } | |||
720 | ||||
721 | /* Determine the ns grid cell limits for this DD zone */ | |||
722 | for (d = 0; d < DIM3; d++) | |||
723 | { | |||
724 | shift0[d] = dd_zones->shift[zone][d]; | |||
725 | useall[d] = (shift0[d] == 0 || d >= grid->npbcdim); | |||
726 | /* Check if we need to do normal or optimized grid assignments. | |||
727 | * Normal is required for dims without DD or triclinic dims. | |||
728 | * DD edge cell on dims without pbc will be automatically | |||
729 | * be correct, since the shift=0 zones with have b0 and b1 | |||
730 | * set to the grid boundaries and there are no shift=1 zones. | |||
731 | */ | |||
732 | if (grid->ncpddc[d] == 0) | |||
733 | { | |||
734 | b0[d] = 0; | |||
735 | b1[d] = grid->n[d]; | |||
736 | } | |||
737 | else | |||
738 | { | |||
739 | if (shift0[d] == 0) | |||
740 | { | |||
741 | b0[d] = 0; | |||
742 | b1[d] = grid->ncpddc[d]; | |||
743 | } | |||
744 | else | |||
745 | { | |||
746 | /* shift = 1 */ | |||
747 | b0[d] = grid->ncpddc[d]; | |||
748 | b1[d] = grid->n[d]; | |||
749 | } | |||
750 | } | |||
751 | } | |||
752 | ||||
753 | not_used = ci_not_used(grid->n); | |||
754 | ||||
755 | /* Put all the charge groups of this DD zone on the grid */ | |||
756 | for (cg = ccg0; cg < ccg1; cg++) | |||
757 | { | |||
758 | if (cell_index[cg] == -1) | |||
759 | { | |||
760 | /* This cg has moved to another node */ | |||
761 | cell_index[cg] = NSGRID_SIGNAL_MOVED_FAC4*grid->ncells; | |||
762 | continue; | |||
763 | } | |||
764 | ||||
765 | bUse = TRUE1; | |||
766 | for (d = 0; d < DIM3; d++) | |||
767 | { | |||
768 | ind[d] = (cg_cm[cg][d] - offset[d])*n_box[d]; | |||
769 | /* Here we have to correct for rounding problems, | |||
770 | * as this cg_cm to cell index operation is not necessarily | |||
771 | * binary identical to the operation for the DD zone assignment | |||
772 | * and therefore a cg could end up in an unused grid cell. | |||
773 | * For dimensions without pbc we need to check | |||
774 | * for cells on the edge if charge groups are beyond | |||
775 | * the grid and if so, store them in the closest cell. | |||
776 | */ | |||
777 | if (ind[d] < b0[d]) | |||
778 | { | |||
779 | ind[d] = b0[d]; | |||
780 | } | |||
781 | else if (ind[d] >= b1[d]) | |||
782 | { | |||
783 | if (useall[d]) | |||
784 | { | |||
785 | ind[d] = b1[d] - 1; | |||
786 | } | |||
787 | else | |||
788 | { | |||
789 | /* Charge groups in this DD zone further away than the cut-off | |||
790 | * in direction do not participate in non-bonded interactions. | |||
791 | */ | |||
792 | bUse = FALSE0; | |||
793 | } | |||
794 | } | |||
795 | } | |||
796 | if (cg > grid->nr_alloc) | |||
797 | { | |||
798 | fprintf(stderrstderr, "WARNING: nra_alloc %d cg0 %d cg1 %d cg %d\n", | |||
799 | grid->nr_alloc, cg0, cg1, cg); | |||
800 | } | |||
801 | if (bUse) | |||
802 | { | |||
803 | cell_index[cg] = xyz2ci(nry, nrz, ind[XX], ind[YY], ind[ZZ])((nry)*(nrz)*(ind[0])+(nrz)*(ind[1])+(ind[2])); | |||
804 | } | |||
805 | else | |||
806 | { | |||
807 | cell_index[cg] = not_used; | |||
808 | } | |||
809 | } | |||
810 | } | |||
811 | } | |||
812 | debug_gmx(); | |||
813 | ||||
814 | } | |||
815 | ||||
816 | void check_grid(t_grid *grid) | |||
817 | { | |||
818 | int ix, iy, iz, ci, cci, nra; | |||
819 | ||||
820 | if (grid->ncells <= 0) | |||
821 | { | |||
822 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c", 822, "Number of grid cells is zero. Probably the system and box collapsed.\n"); | |||
823 | } | |||
824 | ||||
825 | ci = 0; | |||
826 | cci = 0; | |||
827 | for (ix = 0; (ix < grid->n[XX0]); ix++) | |||
828 | { | |||
829 | for (iy = 0; (iy < grid->n[YY1]); iy++) | |||
830 | { | |||
831 | for (iz = 0; (iz < grid->n[ZZ2]); iz++, ci++) | |||
832 | { | |||
833 | if (ci > 0) | |||
834 | { | |||
835 | nra = grid->index[ci]-grid->index[cci]; | |||
836 | if (nra != grid->nra[cci]) | |||
837 | { | |||
838 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c", 838, "nra=%d, grid->nra=%d, cci=%d", | |||
839 | nra, grid->nra[cci], cci); | |||
840 | } | |||
841 | } | |||
842 | cci = xyz2ci(grid->n[YY], grid->n[ZZ], ix, iy, iz)((grid->n[1])*(grid->n[2])*(ix)+(grid->n[2])*(iy)+(iz )); | |||
843 | range_check_mesg(cci, 0, grid->ncells, range_warn)_range_check(cci, 0, grid->ncells, range_warn,"cci", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c" , 843); | |||
844 | ||||
845 | if (cci != ci) | |||
846 | { | |||
847 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/nsgrid.c", 847, "ci = %d, cci = %d", ci, cci); | |||
848 | } | |||
849 | } | |||
850 | } | |||
851 | } | |||
852 | } | |||
853 | ||||
854 | void print_grid(FILE *log, t_grid *grid) | |||
855 | { | |||
856 | int i, nra, index; | |||
857 | int ix, iy, iz, ci; | |||
858 | ||||
859 | fprintf(log, "nr: %d\n", grid->nr); | |||
860 | fprintf(log, "nrx: %d\n", grid->n[XX0]); | |||
861 | fprintf(log, "nry: %d\n", grid->n[YY1]); | |||
862 | fprintf(log, "nrz: %d\n", grid->n[ZZ2]); | |||
863 | fprintf(log, "ncg_ideal: %d\n", grid->ncg_ideal); | |||
864 | fprintf(log, " i cell_index\n"); | |||
865 | for (i = 0; (i < grid->nr); i++) | |||
866 | { | |||
867 | fprintf(log, "%5d %5d\n", i, grid->cell_index[i]); | |||
868 | } | |||
869 | fprintf(log, "cells\n"); | |||
870 | fprintf(log, " ix iy iz nr index cgs...\n"); | |||
871 | ci = 0; | |||
872 | for (ix = 0; (ix < grid->n[XX0]); ix++) | |||
873 | { | |||
874 | for (iy = 0; (iy < grid->n[YY1]); iy++) | |||
875 | { | |||
876 | for (iz = 0; (iz < grid->n[ZZ2]); iz++, ci++) | |||
877 | { | |||
878 | index = grid->index[ci]; | |||
879 | nra = grid->nra[ci]; | |||
880 | fprintf(log, "%3d%3d%3d%5d%5d", ix, iy, iz, nra, index); | |||
881 | for (i = 0; (i < nra); i++) | |||
882 | { | |||
883 | fprintf(log, "%5d", grid->a[index+i]); | |||
884 | } | |||
885 | fprintf(log, "\n"); | |||
886 | } | |||
887 | } | |||
888 | } | |||
889 | fflush(log); | |||
890 | } |