src/gromacs/mdlib/nsgrid.cpp

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