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