1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
31 * For more info, check our website at http://www.gromacs.org
34 * GROwing Monsters And Cloning Shrimps
36 /* This file is completely threadsafe - keep it that way! */
45 #include "types/commrec.h"
49 #include "gmx_fatal.h"
59 /***********************************
61 ***********************************/
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.";
72 static void calc_x_av_stddev(int n,rvec *x,rvec av,rvec stddev)
85 s2[d] += x[i][d]*x[i][d];
95 stddev[d] = sqrt(s2[d] - s1[d]*s1[d]);
99 void get_nsgrid_boundaries_vac(real av,real stddev,
100 real *bound0,real *bound1,
101 real *bdens0,real *bdens1)
103 /* Set the grid to 2 times the standard deviation of
104 * the charge group centers in both directions.
105 * For a uniform bounded distribution the width is sqrt(3)*stddev,
106 * so all charge groups fall within the width.
107 * For a sphere stddev is r/sqrt(5): 99.2% falls within the width.
108 * For a Gaussian distribution 98% fall within the width.
110 *bound0 = av - NSGRID_STDDEV_FAC*stddev;
111 *bound1 = av + NSGRID_STDDEV_FAC*stddev;
113 *bdens0 = av - GRID_STDDEV_FAC*stddev;
114 *bdens1 = av + GRID_STDDEV_FAC*stddev;
117 static void dd_box_bounds_to_ns_bounds(real box0,real box_size,
122 /* Redetermine av and stddev from the DD box boundaries */
123 av = box0 + 0.5*box_size;
124 stddev = 0.5*box_size/GRID_STDDEV_FAC;
126 *gr0 = av - NSGRID_STDDEV_FAC*stddev;
127 *gr1 = av + NSGRID_STDDEV_FAC*stddev;
130 void get_nsgrid_boundaries(int nboundeddim,matrix box,
132 gmx_ddbox_t *ddbox,rvec *gr0,rvec *gr1,
134 rvec grid_x0,rvec grid_x1,
138 real vol,bdens0,bdens1;
141 if (nboundeddim < DIM)
143 calc_x_av_stddev(ncg,cgcm,av,stddev);
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]);
159 get_nsgrid_boundaries_vac(av[d],stddev[d],
160 &grid_x0[d],&grid_x1[d],
165 /* Temporary fix which uses global ddbox boundaries
166 * for unbounded dimensions.
167 * Should be replaced by local boundaries, which makes
168 * the ns grid smaller and does not require global comm.
170 dd_box_bounds_to_ns_bounds(ddbox->box0[d],ddbox->box_size[d],
171 &grid_x0[d],&grid_x1[d]);
175 /* Check for a DD cell not at a lower edge */
176 if (dd != NULL && gr0 != NULL && dd->ci[d] > 0)
178 grid_x0[d] = (*gr0)[d];
181 /* Check for a DD cell not at a higher edge */
182 if (dd != NULL && gr1 != NULL && dd->ci[d] < dd->nc[d]-1)
184 grid_x1[d] = (*gr1)[d];
187 vol *= (bdens1 - bdens0);
192 fprintf(debug,"Set grid boundaries dim %d: %f %f\n",
193 d,grid_x0[d],grid_x1[d]);
197 *grid_density = ncg/vol;
200 static void set_grid_sizes(matrix box,rvec izones_x0,rvec izones_x1,real rlist,
201 const gmx_domdec_t *dd,const gmx_ddbox_t *ddbox,
206 gmx_bool bDD,bDDRect;
209 real inv_r_ideal,size,add_tric,radd;
211 for(i=0; (i<DIM); i++)
216 "set_grid_sizes, i-zone bounds for dim %d: %6.3f %6.3f\n",
217 i,izones_x0[i],izones_x1[i]);
219 izones_size[i] = izones_x1[i] - izones_x0[i];
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)
226 inv_r_ideal = 1/rlist;
230 fprintf(debug,"CG density %f ideal ns cell size %f\n",
231 grid_density,1/inv_r_ideal);
234 clear_rvec(grid->cell_offset);
235 for(i=0; (i<DIM); i++)
237 /* Initial settings, for DD might change below */
238 grid->cell_offset[i] = izones_x0[i];
239 size = izones_size[i];
241 bDD = dd && (dd->nc[i] > 1);
248 /* With DD grid cell jumps only the first decomposition
249 * direction has uniform DD cell boundaries.
251 bDDRect = !(ddbox->tric_dir[i] ||
252 (dd->bGridJump && i != dd->dim[0]));
255 if (i >= ddbox->npbcdim &&
257 izones_x1[i] + radd > ddbox->box0[i] + ddbox->box_size[i]))
259 radd = ddbox->box0[i] + ddbox->box_size[i] - izones_x1[i];
266 /* With DD we only need a grid of one DD cell size + rlist */
273 size += radd/ddbox->skew_fac[i];
276 /* Check if the cell boundary in this direction is
277 * perpendicular to the Cartesian axis.
279 for(j=i+1; j<grid->npbcdim; j++)
283 /* Correct the offset for the home cell location */
284 grid->cell_offset[i] += izones_x0[j]*box[j][i]/box[j][j];
286 /* Correct the offset and size for the off-diagonal
287 * displacement of opposing DD cell corners.
289 /* Without rouding we would need to add:
290 * box[j][i]*rlist/(dd->skew_fac[i]*box[j][j])
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)
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.
301 add_tric += -box[YY][XX]*box[ZZ][YY]/box[YY][YY];
305 grid->cell_offset[i] += add_tric;
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.
321 grid->n[i] = (int)(size*inv_r_ideal + 0.5);
326 grid->cell_size[i] = size/grid->n[i];
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.
336 grid->ncpddc[i] = (int)(izones_size[i]*inv_r_ideal + 0.5);
337 if (grid->ncpddc[i] < 2)
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;
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]);
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]);
360 t_grid *init_grid(FILE *fplog,t_forcerec *fr)
368 grid->npbcdim = ePBC2npbcdim(fr->ePBC);
370 if (fr->ePBC == epbcXY && fr->nwall == 2)
372 grid->nboundeddim = 3;
376 grid->nboundeddim = grid->npbcdim;
381 fprintf(debug,"The coordinates are bounded in %d dimensions\n",
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");
390 sscanf(ptr,"%d",&grid->ncg_ideal);
393 fprintf(fplog,"Set ncg_ideal to %d\n",grid->ncg_ideal);
395 if (grid->ncg_ideal <= 0)
397 gmx_fatal(FARGS,"The number of cg's per cell should be > 0");
402 fprintf(debug,"Set ncg_ideal to %d\n",grid->ncg_ideal);
408 void done_grid(t_grid *grid)
413 sfree(grid->cell_index);
417 grid->cells_nalloc = 0;
424 fprintf(debug,"Successfully freed memory for grid pointers.");
427 int xyz2ci_(int nry,int nrz,int x,int y,int z)
428 /* Return the cell index */
430 return (nry*nrz*x+nrz*y+z);
433 void ci2xyz(t_grid *grid, int i, int *x, int *y, int *z)
434 /* Return x,y and z from the cell index */
438 range_check_mesg(i,0,grid->nr,range_warn);
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];
448 static int ci_not_used(ivec n)
450 /* Return an improbable value */
451 return xyz2ci(n[YY],n[ZZ],3*n[XX],3*n[YY],3*n[ZZ]);
454 static void set_grid_ncg(t_grid *grid,int 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);
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)
477 set_grid_sizes(box,izones_x0,izones_x1,rlistlong,dd,ddbox,grid,grid_density);
479 grid->ncells = grid->n[XX]*grid->n[YY]*grid->n[ZZ];
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);
488 fprintf(fplog,"Grid: %d x %d x %d cells\n",
489 grid->n[XX],grid->n[YY],grid->n[ZZ]);
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);
502 for(i=0; (i<grid->ncells); i++) {
507 static void calc_bor(int cg0,int cg1,int ncg,int CG0[2],int CG1[2])
524 fprintf(debug,"calc_bor: cg0=%d, cg1=%d, ncg=%d\n",cg0,cg1,ncg);
526 fprintf(debug,"CG0[%d]=%d, CG1[%d]=%d\n",m,CG0[m],m,CG1[m]);
531 void calc_elemnr(FILE *fplog,t_grid *grid,int cg0,int cg1,int ncg)
534 int *cell_index=grid->cell_index;
541 gmx_fatal(FARGS,"Number of grid cells is zero. Probably the system and box collapsed.\n");
543 not_used = ci_not_used(grid->n);
545 calc_bor(cg0,cg1,ncg,CG0,CG1);
547 for(i=CG0[m]; (i<CG1[m]); i++) {
549 if (ci != not_used) {
550 range_check_mesg(ci,0,ncells,range_warn);
556 void calc_ptrs(t_grid *grid)
558 int *index = grid->index;
559 int *nra = grid->nra;
565 gmx_fatal(FARGS,"Number of grid cells is zero. Probably the system and box collapsed.\n");
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);
579 void grid_last(FILE *log,t_grid *grid,int cg0,int cg1,int ncg)
583 int ci,not_used,ind,ncells;
584 int *cell_index = grid->cell_index;
585 int *nra = grid->nra;
586 int *index = grid->index;
591 gmx_fatal(FARGS,"Number of grid cells is zero. Probably the system and box collapsed.\n");
593 not_used = ci_not_used(grid->n);
595 calc_bor(cg0,cg1,ncg,CG0,CG1);
597 for(i=CG0[m]; (i<CG1[m]); 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);
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[])
616 int zone,ccg0,ccg1,cg,d,not_used;
617 ivec shift0,useall,b0,b1,ind;
622 /* We have already filled the grid up to grid->ncg,
623 * continue from there.
628 set_grid_ncg(grid,ncg_tot);
630 cell_index = grid->cell_index;
632 /* Initiate cell borders */
638 if (grid->cell_size[d] > 0)
640 n_box[d] = 1/grid->cell_size[d];
647 copy_rvec(grid->cell_offset,offset);
651 fprintf(debug,"Filling grid from %d to %d\n",cg0,cg1);
655 if (dd_zones == NULL)
657 for (cg=cg0; cg<cg1; cg++)
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.
670 else if (ind[d] >= grid->n[d])
672 ind[d] = grid->n[d] - 1;
675 cell_index[cg] = xyz2ci(nry,nrz,ind[XX],ind[YY],ind[ZZ]);
680 for(zone=0; zone<dd_zones->n; zone++)
682 ccg0 = dd_zones->cg_range[zone];
683 ccg1 = dd_zones->cg_range[zone+1];
684 if (ccg1 <= cg0 || ccg0 >= cg1)
689 /* Determine the ns grid cell limits for this DD zone */
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.
700 if (grid->ncpddc[d] == 0)
710 b1[d] = grid->ncpddc[d];
715 b0[d] = grid->ncpddc[d];
721 not_used = ci_not_used(grid->n);
723 /* Put all the charge groups of this DD zone on the grid */
724 for(cg=ccg0; cg<ccg1; cg++)
726 if (cell_index[cg] == -1)
728 /* This cg has moved to another node */
729 cell_index[cg] = NSGRID_SIGNAL_MOVED_FAC*grid->ncells;
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.
745 if (ind[d] < b0[d]) {
748 else if (ind[d] >= b1[d])
756 /* Charge groups in this DD zone further away than the cut-off
757 * in direction do not participate in non-bonded interactions.
763 if (cg > grid->nr_alloc)
765 fprintf(stderr,"WARNING: nra_alloc %d cg0 %d cg1 %d cg %d\n",
766 grid->nr_alloc,cg0,cg1,cg);
770 cell_index[cg] = xyz2ci(nry,nrz,ind[XX],ind[YY],ind[ZZ]);
774 cell_index[cg] = not_used;
783 void check_grid(FILE *log,t_grid *grid)
785 int ix,iy,iz,ci,cci,nra;
788 gmx_fatal(FARGS,"Number of grid cells is zero. Probably the system and box collapsed.\n");
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++) {
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);
801 cci=xyz2ci(grid->n[YY],grid->n[ZZ],ix,iy,iz);
802 range_check_mesg(cci,0,grid->ncells,range_warn);
805 gmx_fatal(FARGS,"ci = %d, cci = %d",ci,cci);
809 void print_grid(FILE *log,t_grid *grid)
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");
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];
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]);
838 void mv_grid(t_commrec *cr,t_grid *grid)
843 #define next ((cur+1) % (cr->nnodes-cr->npmenodes))
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));
852 start = cgindex[next];
853 nr = cgindex[next+1] - start;
854 gmx_rx(cr,GMX_RIGHT,&(ci[start]),nr*sizeof(*ci));
856 gmx_tx_wait(cr, GMX_LEFT);
857 gmx_rx_wait(cr, GMX_RIGHT);