2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2017,2018, 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.
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.
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.
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.
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.
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.
37 /* This file is completely threadsafe - keep it that way! */
48 #include "gromacs/domdec/domdec.h"
49 #include "gromacs/domdec/domdec_struct.h"
50 #include "gromacs/fileio/pdbio.h"
51 #include "gromacs/gmxlib/network.h"
52 #include "gromacs/math/vec.h"
53 #include "gromacs/mdtypes/forcerec.h"
54 #include "gromacs/pbcutil/pbc.h"
55 #include "gromacs/utility/fatalerror.h"
56 #include "gromacs/utility/futil.h"
57 #include "gromacs/utility/smalloc.h"
59 /***********************************
61 ***********************************/
63 static const char *range_warn =
64 "Explanation: During neighborsearching, we assign each particle to a grid\n"
65 "based on its coordinates. If your system contains collisions or parameter\n"
66 "errors that give particles very high velocities you might end up with some\n"
67 "coordinates being +-Infinity or NaN (not-a-number). Obviously, we cannot\n"
68 "put these on a grid, so this is usually where we detect those errors.\n"
69 "Make sure your system is properly energy-minimized and that the potential\n"
70 "energy seems reasonable before trying again.";
72 static void calc_x_av_stddev(int n, rvec *x, rvec av, rvec stddev)
80 for (i = 0; i < n; i++)
82 for (d = 0; d < DIM; d++)
85 s2[d] += x[i][d]*x[i][d];
89 dsvmul(1.0/n, s1, s1);
90 dsvmul(1.0/n, s2, s2);
92 for (d = 0; d < DIM; d++)
95 stddev[d] = std::sqrt(s2[d] - s1[d]*s1[d]);
99 static void get_nsgrid_boundaries_vac(real av, real stddev,
100 real *bound0, real *bound1,
101 real *bdens0, real *bdens1)
103 /* Set the grid to 2 times the standard deviation of
104 * the charge group centers in both directions.
105 * For a uniform bounded distribution the width is sqrt(3)*stddev,
106 * so all charge groups fall within the width.
107 * For a sphere stddev is r/sqrt(5): 99.2% falls within the width.
108 * For a Gaussian distribution 98% fall within the width.
110 *bound0 = av - NSGRID_STDDEV_FAC*stddev;
111 *bound1 = av + NSGRID_STDDEV_FAC*stddev;
113 *bdens0 = av - GRID_STDDEV_FAC*stddev;
114 *bdens1 = av + GRID_STDDEV_FAC*stddev;
117 static void dd_box_bounds_to_ns_bounds(real box0, real box_size,
118 real *gr0, real *gr1)
122 /* Redetermine av and stddev from the DD box boundaries */
123 av = box0 + 0.5*box_size;
124 stddev = 0.5*box_size/GRID_STDDEV_FAC;
126 *gr0 = av - NSGRID_STDDEV_FAC*stddev;
127 *gr1 = av + NSGRID_STDDEV_FAC*stddev;
130 void get_nsgrid_boundaries(int nboundeddim, matrix box,
132 gmx_ddbox_t *ddbox, rvec *gr0, rvec *gr1,
134 rvec grid_x0, rvec grid_x1,
138 real vol, bdens0, bdens1;
141 if (nboundeddim < DIM)
143 calc_x_av_stddev(ncg, cgcm, av, stddev);
147 for (d = 0; d < DIM; d++)
151 grid_x0[d] = (gr0 != nullptr ? (*gr0)[d] : 0);
152 grid_x1[d] = (gr1 != nullptr ? (*gr1)[d] : box[d][d]);
153 vol *= (grid_x1[d] - grid_x0[d]);
157 if (ddbox == nullptr)
159 get_nsgrid_boundaries_vac(av[d], stddev[d],
160 &grid_x0[d], &grid_x1[d],
165 /* Temporary fix which uses global ddbox boundaries
166 * for unbounded dimensions.
167 * Should be replaced by local boundaries, which makes
168 * the ns grid smaller and does not require global comm.
170 dd_box_bounds_to_ns_bounds(ddbox->box0[d], ddbox->box_size[d],
171 &grid_x0[d], &grid_x1[d]);
175 /* Check for a DD cell not at a lower edge */
176 if (dd != nullptr && gr0 != nullptr && dd->ci[d] > 0)
178 grid_x0[d] = (*gr0)[d];
181 /* Check for a DD cell not at a higher edge */
182 if (dd != nullptr && gr1 != nullptr && dd->ci[d] < dd->nc[d]-1)
184 grid_x1[d] = (*gr1)[d];
187 vol *= (bdens1 - bdens0);
192 fprintf(debug, "Set grid boundaries dim %d: %f %f\n",
193 d, grid_x0[d], grid_x1[d]);
197 *grid_density = ncg/vol;
200 static void set_grid_sizes(matrix box, rvec izones_x0, rvec izones_x1, real rlist,
201 const gmx_domdec_t *dd, const gmx_ddbox_t *ddbox,
206 gmx_bool bDD, bDDRect;
208 real inv_r_ideal, size, add_tric, radd;
210 for (i = 0; (i < DIM); i++)
215 "set_grid_sizes, i-zone bounds for dim %d: %6.3f %6.3f\n",
216 i, izones_x0[i], izones_x1[i]);
218 izones_size[i] = izones_x1[i] - izones_x0[i];
221 /* Use the ideal number of cg's per cell to set the ideal cell size */
222 inv_r_ideal = std::cbrt(grid_density/grid->ncg_ideal);
223 if (rlist > 0 && inv_r_ideal*rlist < 1)
225 inv_r_ideal = 1/rlist;
229 fprintf(debug, "CG density %f ideal ns cell size %f\n",
230 grid_density, 1/inv_r_ideal);
233 clear_rvec(grid->cell_offset);
234 for (i = 0; (i < DIM); i++)
236 /* Initial settings, for DD might change below */
237 grid->cell_offset[i] = izones_x0[i];
238 size = izones_size[i];
240 bDD = dd && (dd->nc[i] > 1);
247 /* With DD grid cell jumps only the first decomposition
248 * direction has uniform DD cell boundaries.
250 bDDRect = !(ddbox->tric_dir[i] ||
251 (dd_dlb_is_on(dd) && i != dd->dim[0]));
254 if (i >= ddbox->npbcdim &&
256 izones_x1[i] + radd > ddbox->box0[i] + ddbox->box_size[i]))
258 radd = ddbox->box0[i] + ddbox->box_size[i] - izones_x1[i];
265 /* With DD we only need a grid of one DD cell size + rlist */
272 size += radd/ddbox->skew_fac[i];
275 /* Check if the cell boundary in this direction is
276 * perpendicular to the Cartesian axis.
277 * Since grid->npbcdim isan integer that in principle can take
278 * any value, we help the compiler avoid warnings and potentially
279 * optimize by ensuring that j < DIM here.
281 for (j = i+1; j < grid->npbcdim && j < DIM; j++)
285 /* Correct the offset for the home cell location */
286 grid->cell_offset[i] += izones_x0[j]*box[j][i]/box[j][j];
288 /* Correct the offset and size for the off-diagonal
289 * displacement of opposing DD cell corners.
291 /* Without rouding we would need to add:
292 * box[j][i]*rlist/(dd->skew_fac[i]*box[j][j])
294 /* Determine the shift for the corners of the triclinic box */
295 add_tric = izones_size[j]*box[j][i]/box[j][j];
296 if (dd->ndim == 1 && j == ZZ)
298 /* With 1D domain decomposition the cg's are not in
299 * the triclinic box, but trilinic x-y and rectangular y-z.
300 * Therefore we need to add the shift from the trilinic
301 * corner to the corner at y=0.
303 add_tric += -box[YY][XX]*box[ZZ][YY]/box[YY][YY];
307 grid->cell_offset[i] += add_tric;
319 /* No DD or the box is triclinic is this direction:
320 * we will use the normal grid ns that checks all cells
321 * that are within cut-off distance of the i-particle.
323 grid->n[i] = static_cast<int>(size*inv_r_ideal + 0.5);
328 grid->cell_size[i] = size/grid->n[i];
333 /* We use grid->ncpddc[i] such that all particles
334 * in one ns cell belong to a single DD cell only.
335 * We can then beforehand exclude certain ns grid cells
336 * for non-home i-particles.
338 grid->ncpddc[i] = static_cast<int>(izones_size[i]*inv_r_ideal + 0.5);
339 if (grid->ncpddc[i] < 2)
343 grid->cell_size[i] = izones_size[i]/grid->ncpddc[i];
344 grid->n[i] = grid->ncpddc[i] + static_cast<int>(radd/grid->cell_size[i]) + 1;
348 fprintf(debug, "grid dim %d size %d x %f: %f - %f\n",
349 i, grid->n[i], grid->cell_size[i],
350 grid->cell_offset[i],
351 grid->cell_offset[i]+grid->n[i]*grid->cell_size[i]);
357 fprintf(debug, "CG ncg ideal %d, actual density %.1f\n",
358 grid->ncg_ideal, grid_density*grid->cell_size[XX]*grid->cell_size[YY]*grid->cell_size[ZZ]);
362 t_grid *init_grid(FILE *fplog, t_forcerec *fr)
369 grid->npbcdim = ePBC2npbcdim(fr->ePBC);
371 if (fr->ePBC == epbcXY && fr->nwall == 2)
373 grid->nboundeddim = 3;
377 grid->nboundeddim = grid->npbcdim;
382 fprintf(debug, "The coordinates are bounded in %d dimensions\n",
386 /* The ideal number of cg's per ns grid cell seems to be 10 */
387 grid->ncg_ideal = 10;
388 ptr = getenv("GMX_NSCELL_NCG");
391 sscanf(ptr, "%d", &grid->ncg_ideal);
394 fprintf(fplog, "Set ncg_ideal to %d\n", grid->ncg_ideal);
396 if (grid->ncg_ideal <= 0)
398 gmx_fatal(FARGS, "The number of cg's per cell should be > 0");
403 fprintf(debug, "Set ncg_ideal to %d\n", grid->ncg_ideal);
409 void done_grid(t_grid *grid)
418 sfree(grid->cell_index);
422 grid->cells_nalloc = 0;
430 fprintf(debug, "Successfully freed memory for grid pointers.");
435 int xyz2ci_(int nry, int nrz, int x, int y, int z)
436 /* Return the cell index */
438 return (nry*nrz*x+nrz*y+z);
441 void ci2xyz(t_grid *grid, int i, int *x, int *y, int *z)
442 /* Return x,y and z from the cell index */
446 range_check_mesg(i, 0, grid->nr, range_warn);
448 ci = grid->cell_index[i];
449 *x = ci / (grid->n[YY]*grid->n[ZZ]);
450 ci -= (*x)*grid->n[YY]*grid->n[ZZ];
451 *y = ci / grid->n[ZZ];
452 ci -= (*y)*grid->n[ZZ];
456 static int ci_not_used(const ivec n)
458 /* Return an improbable value */
459 return xyz2ci(n[YY], n[ZZ], 3*n[XX], 3*n[YY], 3*n[ZZ]);
462 static void set_grid_ncg(t_grid *grid, int ncg)
467 if (grid->nr+1 > grid->nr_alloc)
469 nr_old = grid->nr_alloc;
470 grid->nr_alloc = over_alloc_dd(grid->nr) + 1;
471 srenew(grid->cell_index, grid->nr_alloc);
472 for (i = nr_old; i < grid->nr_alloc; i++)
474 grid->cell_index[i] = 0;
476 srenew(grid->a, grid->nr_alloc);
480 void grid_first(FILE *fplog, t_grid *grid,
481 gmx_domdec_t *dd, const gmx_ddbox_t *ddbox,
482 matrix box, rvec izones_x0, rvec izones_x1,
483 real rlist, real grid_density)
487 set_grid_sizes(box, izones_x0, izones_x1, rlist, dd, ddbox, grid, grid_density);
489 grid->ncells = grid->n[XX]*grid->n[YY]*grid->n[ZZ];
491 if (grid->ncells+1 > grid->cells_nalloc)
493 /* Allocate double the size so we have some headroom */
494 grid->cells_nalloc = 2*grid->ncells;
495 srenew(grid->nra, grid->cells_nalloc+1);
496 srenew(grid->index, grid->cells_nalloc+1);
500 fprintf(fplog, "Grid: %d x %d x %d cells\n",
501 grid->n[XX], grid->n[YY], grid->n[ZZ]);
505 m = std::max(grid->n[XX], std::max(grid->n[YY], grid->n[ZZ]));
506 if (m > grid->dc_nalloc)
508 /* Allocate with double the initial size for box scaling */
509 grid->dc_nalloc = 2*m;
510 srenew(grid->dcx2, grid->dc_nalloc);
511 srenew(grid->dcy2, grid->dc_nalloc);
512 srenew(grid->dcz2, grid->dc_nalloc);
516 for (i = 0; (i < grid->ncells); i++)
522 static void calc_bor(int cg0, int cg1, int ncg, int CG0[2], int CG1[2])
542 fprintf(debug, "calc_bor: cg0=%d, cg1=%d, ncg=%d\n", cg0, cg1, ncg);
543 for (m = 0; (m < 2); m++)
545 fprintf(debug, "CG0[%d]=%d, CG1[%d]=%d\n", m, CG0[m], m, CG1[m]);
551 void calc_elemnr(t_grid *grid, int cg0, int cg1, int ncg)
554 int *cell_index = grid->cell_index;
555 int *nra = grid->nra;
559 ncells = grid->ncells;
562 gmx_fatal(FARGS, "Number of grid cells is zero. Probably the system and box collapsed.\n");
565 not_used = ci_not_used(grid->n);
567 calc_bor(cg0, cg1, ncg, CG0, CG1);
568 for (m = 0; (m < 2); m++)
570 for (i = CG0[m]; (i < CG1[m]); i++)
575 range_check_mesg(ci, 0, ncells, range_warn);
582 void calc_ptrs(t_grid *grid)
584 int *index = grid->index;
585 int *nra = grid->nra;
586 int ix, iy, iz, ci, nr;
589 ncells = grid->ncells;
592 gmx_fatal(FARGS, "Number of grid cells is zero. Probably the system and box collapsed.\n");
596 for (ix = 0; (ix < grid->n[XX]); ix++)
598 for (iy = 0; (iy < grid->n[YY]); iy++)
600 for (iz = 0; (iz < grid->n[ZZ]); iz++, ci++)
602 range_check_mesg(ci, 0, ncells, range_warn);
612 void grid_last(t_grid *grid, int cg0, int cg1, int ncg)
616 int ci, not_used, ind, ncells;
617 int *cell_index = grid->cell_index;
618 int *nra = grid->nra;
619 int *index = grid->index;
622 ncells = grid->ncells;
625 gmx_fatal(FARGS, "Number of grid cells is zero. Probably the system and box collapsed.\n");
628 not_used = ci_not_used(grid->n);
630 calc_bor(cg0, cg1, ncg, CG0, CG1);
631 for (m = 0; (m < 2); m++)
633 for (i = CG0[m]; (i < CG1[m]); i++)
638 range_check_mesg(ci, 0, ncells, range_warn);
639 ind = index[ci]+nra[ci]++;
640 range_check_mesg(ind, 0, grid->nr, range_warn);
647 void fill_grid(gmx_domdec_zones_t *dd_zones,
648 t_grid *grid, int ncg_tot,
649 int cg0, int cg1, rvec cg_cm[])
654 int zone, ccg0, ccg1, cg, d, not_used;
655 ivec shift0, useall, b0, b1, ind;
660 /* We have already filled the grid up to grid->ncg,
661 * continue from there.
666 set_grid_ncg(grid, ncg_tot);
668 cell_index = grid->cell_index;
670 /* Initiate cell borders */
673 for (d = 0; d < DIM; d++)
675 if (grid->cell_size[d] > 0)
677 n_box[d] = 1/grid->cell_size[d];
684 copy_rvec(grid->cell_offset, offset);
688 fprintf(debug, "Filling grid from %d to %d\n", cg0, cg1);
691 if (dd_zones == nullptr)
693 for (cg = cg0; cg < cg1; cg++)
695 for (d = 0; d < DIM; d++)
697 ind[d] = static_cast<int>((cg_cm[cg][d] - offset[d])*n_box[d]);
698 /* With pbc we should be done here.
699 * Without pbc cg's outside the grid
700 * should be assigned to the closest grid cell.
706 else if (ind[d] >= grid->n[d])
708 ind[d] = grid->n[d] - 1;
711 cell_index[cg] = xyz2ci(nry, nrz, ind[XX], ind[YY], ind[ZZ]);
716 for (zone = 0; zone < dd_zones->n; zone++)
718 ccg0 = dd_zones->cg_range[zone];
719 ccg1 = dd_zones->cg_range[zone+1];
720 if (ccg1 <= cg0 || ccg0 >= cg1)
725 /* Determine the ns grid cell limits for this DD zone */
726 for (d = 0; d < DIM; d++)
728 shift0[d] = dd_zones->shift[zone][d];
729 useall[d] = (shift0[d] == 0 || d >= grid->npbcdim);
730 /* Check if we need to do normal or optimized grid assignments.
731 * Normal is required for dims without DD or triclinic dims.
732 * DD edge cell on dims without pbc will be automatically
733 * be correct, since the shift=0 zones with have b0 and b1
734 * set to the grid boundaries and there are no shift=1 zones.
736 if (grid->ncpddc[d] == 0)
746 b1[d] = grid->ncpddc[d];
751 b0[d] = grid->ncpddc[d];
757 not_used = ci_not_used(grid->n);
759 /* Put all the charge groups of this DD zone on the grid */
760 for (cg = ccg0; cg < ccg1; cg++)
762 if (cell_index[cg] == -1)
764 /* This cg has moved to another node */
765 cell_index[cg] = NSGRID_SIGNAL_MOVED_FAC*grid->ncells;
770 for (d = 0; d < DIM; d++)
772 ind[d] = static_cast<int>((cg_cm[cg][d] - offset[d])*n_box[d]);
773 /* Here we have to correct for rounding problems,
774 * as this cg_cm to cell index operation is not necessarily
775 * binary identical to the operation for the DD zone assignment
776 * and therefore a cg could end up in an unused grid cell.
777 * For dimensions without pbc we need to check
778 * for cells on the edge if charge groups are beyond
779 * the grid and if so, store them in the closest cell.
785 else if (ind[d] >= b1[d])
793 /* Charge groups in this DD zone further away than the cut-off
794 * in direction do not participate in non-bonded interactions.
800 if (cg > grid->nr_alloc)
802 fprintf(stderr, "WARNING: nra_alloc %d cg0 %d cg1 %d cg %d\n",
803 grid->nr_alloc, cg0, cg1, cg);
807 cell_index[cg] = xyz2ci(nry, nrz, ind[XX], ind[YY], ind[ZZ]);
811 cell_index[cg] = not_used;
819 void check_grid(t_grid *grid)
821 int ix, iy, iz, ci, cci, nra;
823 if (grid->ncells <= 0)
825 gmx_fatal(FARGS, "Number of grid cells is zero. Probably the system and box collapsed.\n");
830 for (ix = 0; (ix < grid->n[XX]); ix++)
832 for (iy = 0; (iy < grid->n[YY]); iy++)
834 for (iz = 0; (iz < grid->n[ZZ]); iz++, ci++)
838 nra = grid->index[ci]-grid->index[cci];
839 if (nra != grid->nra[cci])
841 gmx_fatal(FARGS, "nra=%d, grid->nra=%d, cci=%d",
842 nra, grid->nra[cci], cci);
845 cci = xyz2ci(grid->n[YY], grid->n[ZZ], ix, iy, iz);
846 range_check_mesg(cci, 0, grid->ncells, range_warn);
850 gmx_fatal(FARGS, "ci = %d, cci = %d", ci, cci);
857 void print_grid(FILE *log, t_grid *grid)
862 fprintf(log, "nr: %d\n", grid->nr);
863 fprintf(log, "nrx: %d\n", grid->n[XX]);
864 fprintf(log, "nry: %d\n", grid->n[YY]);
865 fprintf(log, "nrz: %d\n", grid->n[ZZ]);
866 fprintf(log, "ncg_ideal: %d\n", grid->ncg_ideal);
867 fprintf(log, " i cell_index\n");
868 for (i = 0; (i < grid->nr); i++)
870 fprintf(log, "%5d %5d\n", i, grid->cell_index[i]);
872 fprintf(log, "cells\n");
873 fprintf(log, " ix iy iz nr index cgs...\n");
875 for (ix = 0; (ix < grid->n[XX]); ix++)
877 for (iy = 0; (iy < grid->n[YY]); iy++)
879 for (iz = 0; (iz < grid->n[ZZ]); iz++, ci++)
881 index = grid->index[ci];
883 fprintf(log, "%3d%3d%3d%5d%5d", ix, iy, iz, nra, index);
884 for (i = 0; (i < nra); i++)
886 fprintf(log, "%5d", grid->a[index+i]);