src/gromacs/mdlib/nsgrid.c

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