src/gromacs/mdlib/nsgrid.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROwing Monsters And Cloning Shrimps
  35  */
  36 /* This file is completely threadsafe - keep it that way! */
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include <stdlib.h>
  42
  43 #include "sysstuff.h"
  44 #include "typedefs.h"
  45 #include "types/commrec.h"
  46 #include "macros.h"
  47 #include "smalloc.h"
  48 #include "nsgrid.h"
  49 #include "gmx_fatal.h"
  50 #include "vec.h"
  51 #include "network.h"
  52 #include "domdec.h"
  53 #include "partdec.h"
  54 #include "pbc.h"
  55 #include <stdio.h>
  56 #include "futil.h"
  57 #include "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     fprintf(debug,"Successfully freed memory for grid pointers.");
 425 }
 426
 427 int xyz2ci_(int nry,int nrz,int x,int y,int z)
 428 /* Return the cell index */
 429 {
 430   return (nry*nrz*x+nrz*y+z);
 431 }
 432
 433 void ci2xyz(t_grid *grid, int i, int *x, int *y, int *z)
 434 /* Return x,y and z from the cell index */
 435 {
 436   int ci;
 437
 438   range_check_mesg(i,0,grid->nr,range_warn);
 439
 440   ci = grid->cell_index[i];
 441   *x  = ci / (grid->n[YY]*grid->n[ZZ]);
 442   ci -= (*x)*grid->n[YY]*grid->n[ZZ];
 443   *y  = ci / grid->n[ZZ];
 444   ci -= (*y)*grid->n[ZZ];
 445   *z  = ci;
 446 }
 447
 448 static int ci_not_used(ivec n)
 449 {
 450   /* Return an improbable value */
 451   return xyz2ci(n[YY],n[ZZ],3*n[XX],3*n[YY],3*n[ZZ]);
 452 }
 453
 454 static void set_grid_ncg(t_grid *grid,int ncg)
 455 {
 456   int nr_old,i;
 457
 458   grid->nr = ncg;
 459   if (grid->nr+1 > grid->nr_alloc) {
 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       grid->cell_index[i] = 0;
 465     srenew(grid->a,grid->nr_alloc);
 466   }
 467 }
 468
 469 void grid_first(FILE *fplog,t_grid *grid,
 470                 gmx_domdec_t *dd,const gmx_ddbox_t *ddbox,
 471                 int ePBC,matrix box,rvec izones_x0,rvec izones_x1,
 472                 real rlistlong,real grid_density)
 473 {
 474   int    i,m;
 475   ivec   cx;
 476
 477   set_grid_sizes(box,izones_x0,izones_x1,rlistlong,dd,ddbox,grid,grid_density);
 478
 479   grid->ncells = grid->n[XX]*grid->n[YY]*grid->n[ZZ];
 480
 481   if (grid->ncells+1 > grid->cells_nalloc) {
 482     /* Allocate double the size so we have some headroom */
 483     grid->cells_nalloc = 2*grid->ncells;
 484     srenew(grid->nra,  grid->cells_nalloc+1);
 485     srenew(grid->index,grid->cells_nalloc+1);
 486
 487     if (fplog)
 488       fprintf(fplog,"Grid: %d x %d x %d cells\n",
 489               grid->n[XX],grid->n[YY],grid->n[ZZ]);
 490   }
 491
 492   m = max(grid->n[XX],max(grid->n[YY],grid->n[ZZ]));
 493   if (m > grid->dc_nalloc) {
 494     /* Allocate with double the initial size for box scaling */
 495     grid->dc_nalloc = 2*m;
 496     srenew(grid->dcx2,grid->dc_nalloc);
 497     srenew(grid->dcy2,grid->dc_nalloc);
 498     srenew(grid->dcz2,grid->dc_nalloc);
 499   }
 500
 501   grid->nr = 0;
 502   for(i=0; (i<grid->ncells); i++) {
 503     grid->nra[i] = 0;
 504   }
 505 }
 506
 507 static void calc_bor(int cg0,int cg1,int ncg,int CG0[2],int CG1[2])
 508 {
 509   if (cg1 > ncg) {
 510     CG0[0]=cg0;
 511     CG1[0]=ncg;
 512     CG0[1]=0;
 513     CG1[1]=cg1-ncg;
 514   }
 515   else {
 516     CG0[0]=cg0;
 517     CG1[0]=cg1;
 518     CG0[1]=0;
 519     CG1[1]=0;
 520   }
 521   if (debug) {
 522     int m;
 523
 524     fprintf(debug,"calc_bor: cg0=%d, cg1=%d, ncg=%d\n",cg0,cg1,ncg);
 525     for(m=0; (m<2); m++)
 526       fprintf(debug,"CG0[%d]=%d, CG1[%d]=%d\n",m,CG0[m],m,CG1[m]);
 527   }
 528
 529 }
 530
 531 void calc_elemnr(FILE *fplog,t_grid *grid,int cg0,int cg1,int ncg)
 532 {
 533   int    CG0[2],CG1[2];
 534   int    *cell_index=grid->cell_index;
 535   int    *nra=grid->nra;
 536   int    i,m,ncells;
 537   int    ci,not_used;
 538
 539   ncells=grid->ncells;
 540   if(ncells<=0)
 541     gmx_fatal(FARGS,"Number of grid cells is zero. Probably the system and box collapsed.\n");
 542
 543   not_used = ci_not_used(grid->n);
 544
 545   calc_bor(cg0,cg1,ncg,CG0,CG1);
 546   for(m=0; (m<2); m++)
 547     for(i=CG0[m]; (i<CG1[m]); i++) {
 548       ci = cell_index[i];
 549       if (ci != not_used) {
 550           range_check_mesg(ci,0,ncells,range_warn);
 551         nra[ci]++;
 552       }
 553     }
 554 }
 555
 556 void calc_ptrs(t_grid *grid)
 557 {
 558   int *index = grid->index;
 559   int *nra   = grid->nra;
 560   int ix,iy,iz,ci,nr;
 561   int nnra,ncells;
 562
 563   ncells=grid->ncells;
 564   if(ncells<=0)
 565     gmx_fatal(FARGS,"Number of grid cells is zero. Probably the system and box collapsed.\n");
 566
 567   ci=nr=0;
 568   for(ix=0; (ix < grid->n[XX]); ix++)
 569     for(iy=0; (iy < grid->n[YY]); iy++)
 570       for(iz=0; (iz < grid->n[ZZ]); iz++,ci++) {
 571           range_check_mesg(ci,0,ncells,range_warn);
 572         index[ci] = nr;
 573         nnra      = nra[ci];
 574         nr       += nnra;
 575         nra[ci]   = 0;
 576       }
 577 }
 578
 579 void grid_last(FILE *log,t_grid *grid,int cg0,int cg1,int ncg)
 580 {
 581   int    CG0[2],CG1[2];
 582   int    i,m;
 583   int    ci,not_used,ind,ncells;
 584   int    *cell_index = grid->cell_index;
 585   int    *nra        = grid->nra;
 586   int    *index      = grid->index;
 587   int    *a          = grid->a;
 588
 589   ncells=grid->ncells;
 590   if (ncells <= 0)
 591     gmx_fatal(FARGS,"Number of grid cells is zero. Probably the system and box collapsed.\n");
 592
 593   not_used = ci_not_used(grid->n);
 594
 595   calc_bor(cg0,cg1,ncg,CG0,CG1);
 596   for(m=0; (m<2); m++)
 597     for(i=CG0[m]; (i<CG1[m]); i++) {
 598       ci     = cell_index[i];
 599       if (ci != not_used) {
 600           range_check_mesg(ci,0,ncells,range_warn);
 601         ind    = index[ci]+nra[ci]++;
 602         range_check_mesg(ind,0,grid->nr,range_warn);
 603         a[ind] = i;
 604       }
 605     }
 606 }
 607
 608 void fill_grid(FILE *log,
 609                gmx_domdec_zones_t *dd_zones,
 610                t_grid *grid,int ncg_tot,
 611                int cg0,int cg1,rvec cg_cm[])
 612 {
 613     int    *cell_index;
 614     int    nrx,nry,nrz;
 615     rvec   n_box,offset;
 616     int    zone,ccg0,ccg1,cg,d,not_used;
 617     ivec   shift0,useall,b0,b1,ind;
 618     gmx_bool   bUse;
 619
 620     if (cg0 == -1)
 621     {
 622         /* We have already filled the grid up to grid->ncg,
 623          * continue from there.
 624          */
 625         cg0 = grid->nr;
 626     }
 627
 628     set_grid_ncg(grid,ncg_tot);
 629
 630     cell_index = grid->cell_index;
 631
 632     /* Initiate cell borders */
 633     nrx = grid->n[XX];
 634     nry = grid->n[YY];
 635     nrz = grid->n[ZZ];
 636     for(d=0; d<DIM; d++)
 637     {
 638         if (grid->cell_size[d] > 0)
 639         {
 640             n_box[d] = 1/grid->cell_size[d];
 641         }
 642         else
 643         {
 644             n_box[d] = 0;
 645         }
 646     }
 647     copy_rvec(grid->cell_offset,offset);
 648
 649     if (debug)
 650     {
 651         fprintf(debug,"Filling grid from %d to %d\n",cg0,cg1);
 652     }
 653
 654     debug_gmx();
 655     if (dd_zones == NULL)
 656     {
 657         for (cg=cg0; cg<cg1; cg++)
 658         {
 659             for(d=0; d<DIM; d++)
 660             {
 661                 ind[d] = (cg_cm[cg][d] - offset[d])*n_box[d];
 662                 /* With pbc we should be done here.
 663                  * Without pbc cg's outside the grid
 664                  * should be assigned to the closest grid cell.
 665                  */
 666                 if (ind[d] < 0)
 667                 {
 668                     ind[d] = 0;
 669                 }
 670                 else if (ind[d] >= grid->n[d])
 671                 {
 672                     ind[d] = grid->n[d] - 1;
 673                 }
 674             }
 675             cell_index[cg] = xyz2ci(nry,nrz,ind[XX],ind[YY],ind[ZZ]);
 676         }
 677     }
 678     else
 679     {
 680         for(zone=0; zone<dd_zones->n; zone++)
 681         {
 682             ccg0 = dd_zones->cg_range[zone];
 683             ccg1 = dd_zones->cg_range[zone+1];
 684             if (ccg1 <= cg0 || ccg0 >= cg1)
 685             {
 686                 continue;
 687             }
 688
 689             /* Determine the ns grid cell limits for this DD zone */
 690             for(d=0; d<DIM; d++)
 691             {
 692                 shift0[d] = dd_zones->shift[zone][d];
 693                 useall[d] = (shift0[d] == 0 || d >= grid->npbcdim);
 694                 /* Check if we need to do normal or optimized grid assignments.
 695                  * Normal is required for dims without DD or triclinic dims.
 696                  * DD edge cell on dims without pbc will be automatically
 697                  * be correct, since the shift=0 zones with have b0 and b1
 698                  * set to the grid boundaries and there are no shift=1 zones.
 699                  */
 700                 if (grid->ncpddc[d] == 0)
 701                 {
 702                     b0[d] = 0;
 703                     b1[d] = grid->n[d];
 704                 }
 705                 else
 706                 {
 707                     if (shift0[d] == 0)
 708                     {
 709                         b0[d] = 0;
 710                         b1[d] = grid->ncpddc[d];
 711                     }
 712                     else
 713                     {
 714                         /* shift = 1 */
 715                         b0[d] = grid->ncpddc[d];
 716                         b1[d] = grid->n[d];
 717                     }
 718                 }
 719             }
 720
 721             not_used = ci_not_used(grid->n);
 722
 723             /* Put all the charge groups of this DD zone on the grid */
 724             for(cg=ccg0; cg<ccg1; cg++)
 725             {
 726                 if (cell_index[cg] == -1)
 727                 {
 728                     /* This cg has moved to another node */
 729                     cell_index[cg] = NSGRID_SIGNAL_MOVED_FAC*grid->ncells;
 730                     continue;
 731                 }
 732
 733                 bUse = TRUE;
 734                 for(d=0; d<DIM; d++)
 735                 {
 736                     ind[d] = (cg_cm[cg][d] - offset[d])*n_box[d];
 737                     /* Here we have to correct for rounding problems,
 738                      * as this cg_cm to cell index operation is not necessarily
 739                      * binary identical to the operation for the DD zone assignment
 740                      * and therefore a cg could end up in an unused grid cell.
 741                      * For dimensions without pbc we need to check
 742                      * for cells on the edge if charge groups are beyond
 743                      * the grid and if so, store them in the closest cell.
 744                      */
 745                     if (ind[d] < b0[d]) {
 746                         ind[d] = b0[d];
 747                     }
 748                     else if (ind[d] >= b1[d])
 749                     {
 750                         if (useall[d])
 751                         {
 752                             ind[d] = b1[d] - 1;
 753                         }
 754                         else
 755                         {
 756                             /* Charge groups in this DD zone further away than the cut-off
 757                              * in direction do not participate in non-bonded interactions.
 758                              */
 759                             bUse = FALSE;
 760                         }
 761                     }
 762                 }
 763                 if (cg > grid->nr_alloc)
 764                 {
 765                     fprintf(stderr,"WARNING: nra_alloc %d cg0 %d cg1 %d cg %d\n",
 766                             grid->nr_alloc,cg0,cg1,cg);
 767                 }
 768                 if (bUse)
 769                 {
 770                     cell_index[cg] = xyz2ci(nry,nrz,ind[XX],ind[YY],ind[ZZ]);
 771                 }
 772                 else
 773                 {
 774                     cell_index[cg] = not_used;
 775                 }
 776             }
 777         }
 778     }
 779     debug_gmx();
 780
 781 }
 782
 783 void check_grid(FILE *log,t_grid *grid)
 784 {
 785   int ix,iy,iz,ci,cci,nra;
 786
 787   if(grid->ncells<=0)
 788     gmx_fatal(FARGS,"Number of grid cells is zero. Probably the system and box collapsed.\n");
 789
 790   ci=0;
 791   cci=0;
 792   for(ix=0; (ix<grid->n[XX]); ix++)
 793     for(iy=0; (iy<grid->n[YY]); iy++)
 794       for(iz=0; (iz<grid->n[ZZ]); iz++,ci++) {
 795         if (ci > 0) {
 796           nra=grid->index[ci]-grid->index[cci];
 797           if (nra != grid->nra[cci])
 798             gmx_fatal(FARGS,"nra=%d, grid->nra=%d, cci=%d",
 799                         nra,grid->nra[cci],cci);
 800         }
 801         cci=xyz2ci(grid->n[YY],grid->n[ZZ],ix,iy,iz);
 802         range_check_mesg(cci,0,grid->ncells,range_warn);
 803
 804         if (cci != ci)
 805           gmx_fatal(FARGS,"ci = %d, cci = %d",ci,cci);
 806       }
 807 }
 808
 809 void print_grid(FILE *log,t_grid *grid)
 810 {
 811   int i,nra,index;
 812   int ix,iy,iz,ci;
 813
 814   fprintf(log,"nr:        %d\n",grid->nr);
 815   fprintf(log,"nrx:       %d\n",grid->n[XX]);
 816   fprintf(log,"nry:       %d\n",grid->n[YY]);
 817   fprintf(log,"nrz:       %d\n",grid->n[ZZ]);
 818   fprintf(log,"ncg_ideal: %d\n",grid->ncg_ideal);
 819   fprintf(log,"    i  cell_index\n");
 820   for(i=0; (i<grid->nr); i++)
 821     fprintf(log,"%5d  %5d\n",i,grid->cell_index[i]);
 822   fprintf(log,"cells\n");
 823   fprintf(log," ix iy iz   nr  index  cgs...\n");
 824   ci=0;
 825   for(ix=0; (ix<grid->n[XX]); ix++)
 826     for(iy=0; (iy<grid->n[YY]); iy++)
 827       for(iz=0; (iz<grid->n[ZZ]); iz++,ci++) {
 828         index=grid->index[ci];
 829         nra=grid->nra[ci];
 830         fprintf(log,"%3d%3d%3d%5d%5d",ix,iy,iz,nra,index);
 831         for(i=0; (i<nra); i++)
 832           fprintf(log,"%5d",grid->a[index+i]);
 833         fprintf(log,"\n");
 834       }
 835   fflush(log);
 836 }
 837
 838 void mv_grid(t_commrec *cr,t_grid *grid)
 839 {
 840   int i,start,nr;
 841   int cur=cr->nodeid;
 842   int *ci,*cgindex;
 843 #define next ((cur+1) % (cr->nnodes-cr->npmenodes))
 844
 845   ci = grid->cell_index;
 846   cgindex = pd_cgindex(cr);
 847   for(i=0; (i<cr->nnodes-1); i++) {
 848     start = cgindex[cur];
 849     nr    = cgindex[cur+1] - start;
 850     gmx_tx(cr,GMX_LEFT,&(ci[start]),nr*sizeof(*ci));
 851
 852     start = cgindex[next];
 853     nr    = cgindex[next+1] - start;
 854     gmx_rx(cr,GMX_RIGHT,&(ci[start]),nr*sizeof(*ci));
 855
 856     gmx_tx_wait(cr, GMX_LEFT);
 857     gmx_rx_wait(cr, GMX_RIGHT);
 858
 859     cur=next;
 860   }
 861 }
 862