Bug Summary

File:gromacs/mdlib/nsgrid.c
Location:line 665, column 5
Description:Value stored to 'nrx' 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-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
62const 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
71static 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
98void 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
116static 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
129void 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
199static 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
359t_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
407void 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
428int 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
434void 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
449static 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
455static 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
473void 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
516static 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
545void 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
576void 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
606void 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
641void 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];
Value stored to 'nrx' is never read
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
816void 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
854void 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}