2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
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.
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.
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.
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.
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.
37 /* This file is completely threadsafe - keep it that way! */
49 #include "gromacs/domdec/domdec.h"
50 #include "gromacs/fileio/pdbio.h"
51 #include "gromacs/legacyheaders/network.h"
52 #include "gromacs/legacyheaders/typedefs.h"
53 #include "gromacs/legacyheaders/types/commrec.h"
54 #include "gromacs/math/vec.h"
55 #include "gromacs/pbcutil/pbc.h"
56 #include "gromacs/utility/fatalerror.h"
57 #include "gromacs/utility/futil.h"
58 #include "gromacs/utility/smalloc.h"
60 /***********************************
62 ***********************************/
64 const char *range_warn
=
65 "Explanation: During neighborsearching, we assign each particle to a grid\n"
66 "based on its coordinates. If your system contains collisions or parameter\n"
67 "errors that give particles very high velocities you might end up with some\n"
68 "coordinates being +-Infinity or NaN (not-a-number). Obviously, we cannot\n"
69 "put these on a grid, so this is usually where we detect those errors.\n"
70 "Make sure your system is properly energy-minimized and that the potential\n"
71 "energy seems reasonable before trying again.";
73 static void calc_x_av_stddev(int n
, rvec
*x
, rvec av
, rvec stddev
)
81 for (i
= 0; i
< n
; i
++)
83 for (d
= 0; d
< DIM
; d
++)
86 s2
[d
] += x
[i
][d
]*x
[i
][d
];
90 dsvmul(1.0/n
, s1
, s1
);
91 dsvmul(1.0/n
, s2
, s2
);
93 for (d
= 0; d
< DIM
; d
++)
96 stddev
[d
] = std::sqrt(s2
[d
] - s1
[d
]*s1
[d
]);
100 void get_nsgrid_boundaries_vac(real av
, real stddev
,
101 real
*bound0
, real
*bound1
,
102 real
*bdens0
, real
*bdens1
)
104 /* Set the grid to 2 times the standard deviation of
105 * the charge group centers in both directions.
106 * For a uniform bounded distribution the width is sqrt(3)*stddev,
107 * so all charge groups fall within the width.
108 * For a sphere stddev is r/sqrt(5): 99.2% falls within the width.
109 * For a Gaussian distribution 98% fall within the width.
111 *bound0
= av
- NSGRID_STDDEV_FAC
*stddev
;
112 *bound1
= av
+ NSGRID_STDDEV_FAC
*stddev
;
114 *bdens0
= av
- GRID_STDDEV_FAC
*stddev
;
115 *bdens1
= av
+ GRID_STDDEV_FAC
*stddev
;
118 static void dd_box_bounds_to_ns_bounds(real box0
, real box_size
,
119 real
*gr0
, real
*gr1
)
123 /* Redetermine av and stddev from the DD box boundaries */
124 av
= box0
+ 0.5*box_size
;
125 stddev
= 0.5*box_size
/GRID_STDDEV_FAC
;
127 *gr0
= av
- NSGRID_STDDEV_FAC
*stddev
;
128 *gr1
= av
+ NSGRID_STDDEV_FAC
*stddev
;
131 void get_nsgrid_boundaries(int nboundeddim
, matrix box
,
133 gmx_ddbox_t
*ddbox
, rvec
*gr0
, rvec
*gr1
,
135 rvec grid_x0
, rvec grid_x1
,
139 real vol
, bdens0
, bdens1
;
142 if (nboundeddim
< DIM
)
144 calc_x_av_stddev(ncg
, cgcm
, av
, stddev
);
148 for (d
= 0; d
< DIM
; d
++)
152 grid_x0
[d
] = (gr0
!= NULL
? (*gr0
)[d
] : 0);
153 grid_x1
[d
] = (gr1
!= NULL
? (*gr1
)[d
] : box
[d
][d
]);
154 vol
*= (grid_x1
[d
] - grid_x0
[d
]);
160 get_nsgrid_boundaries_vac(av
[d
], stddev
[d
],
161 &grid_x0
[d
], &grid_x1
[d
],
166 /* Temporary fix which uses global ddbox boundaries
167 * for unbounded dimensions.
168 * Should be replaced by local boundaries, which makes
169 * the ns grid smaller and does not require global comm.
171 dd_box_bounds_to_ns_bounds(ddbox
->box0
[d
], ddbox
->box_size
[d
],
172 &grid_x0
[d
], &grid_x1
[d
]);
176 /* Check for a DD cell not at a lower edge */
177 if (dd
!= NULL
&& gr0
!= NULL
&& dd
->ci
[d
] > 0)
179 grid_x0
[d
] = (*gr0
)[d
];
182 /* Check for a DD cell not at a higher edge */
183 if (dd
!= NULL
&& gr1
!= NULL
&& dd
->ci
[d
] < dd
->nc
[d
]-1)
185 grid_x1
[d
] = (*gr1
)[d
];
188 vol
*= (bdens1
- bdens0
);
193 fprintf(debug
, "Set grid boundaries dim %d: %f %f\n",
194 d
, grid_x0
[d
], grid_x1
[d
]);
198 *grid_density
= ncg
/vol
;
201 static void set_grid_sizes(matrix box
, rvec izones_x0
, rvec izones_x1
, real rlist
,
202 const gmx_domdec_t
*dd
, const gmx_ddbox_t
*ddbox
,
207 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
= std::pow((real
)(grid_density
/grid
->ncg_ideal
), (real
)(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_dlb_is_on(dd
) && 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.
278 * Since grid->npbcdim isan integer that in principle can take
279 * any value, we help the compiler avoid warnings and potentially
280 * optimize by ensuring that j < DIM here.
282 for (j
= i
+1; j
< grid
->npbcdim
&& j
< DIM
; j
++)
286 /* Correct the offset for the home cell location */
287 grid
->cell_offset
[i
] += izones_x0
[j
]*box
[j
][i
]/box
[j
][j
];
289 /* Correct the offset and size for the off-diagonal
290 * displacement of opposing DD cell corners.
292 /* Without rouding we would need to add:
293 * box[j][i]*rlist/(dd->skew_fac[i]*box[j][j])
295 /* Determine the shift for the corners of the triclinic box */
296 add_tric
= izones_size
[j
]*box
[j
][i
]/box
[j
][j
];
297 if (dd
->ndim
== 1 && j
== ZZ
)
299 /* With 1D domain decomposition the cg's are not in
300 * the triclinic box, but trilinic x-y and rectangular y-z.
301 * Therefore we need to add the shift from the trilinic
302 * corner to the corner at y=0.
304 add_tric
+= -box
[YY
][XX
]*box
[ZZ
][YY
]/box
[YY
][YY
];
308 grid
->cell_offset
[i
] += add_tric
;
320 /* No DD or the box is triclinic is this direction:
321 * we will use the normal grid ns that checks all cells
322 * that are within cut-off distance of the i-particle.
324 grid
->n
[i
] = (int)(size
*inv_r_ideal
+ 0.5);
329 grid
->cell_size
[i
] = size
/grid
->n
[i
];
334 /* We use grid->ncpddc[i] such that all particles
335 * in one ns cell belong to a single DD cell only.
336 * We can then beforehand exclude certain ns grid cells
337 * for non-home i-particles.
339 grid
->ncpddc
[i
] = (int)(izones_size
[i
]*inv_r_ideal
+ 0.5);
340 if (grid
->ncpddc
[i
] < 2)
344 grid
->cell_size
[i
] = izones_size
[i
]/grid
->ncpddc
[i
];
345 grid
->n
[i
] = grid
->ncpddc
[i
] + (int)(radd
/grid
->cell_size
[i
]) + 1;
349 fprintf(debug
, "grid dim %d size %d x %f: %f - %f\n",
350 i
, grid
->n
[i
], grid
->cell_size
[i
],
351 grid
->cell_offset
[i
],
352 grid
->cell_offset
[i
]+grid
->n
[i
]*grid
->cell_size
[i
]);
358 fprintf(debug
, "CG ncg ideal %d, actual density %.1f\n",
359 grid
->ncg_ideal
, grid_density
*grid
->cell_size
[XX
]*grid
->cell_size
[YY
]*grid
->cell_size
[ZZ
]);
363 t_grid
*init_grid(FILE *fplog
, t_forcerec
*fr
)
370 grid
->npbcdim
= ePBC2npbcdim(fr
->ePBC
);
372 if (fr
->ePBC
== epbcXY
&& fr
->nwall
== 2)
374 grid
->nboundeddim
= 3;
378 grid
->nboundeddim
= grid
->npbcdim
;
383 fprintf(debug
, "The coordinates are bounded in %d dimensions\n",
387 /* The ideal number of cg's per ns grid cell seems to be 10 */
388 grid
->ncg_ideal
= 10;
389 ptr
= getenv("GMX_NSCELL_NCG");
392 sscanf(ptr
, "%d", &grid
->ncg_ideal
);
395 fprintf(fplog
, "Set ncg_ideal to %d\n", grid
->ncg_ideal
);
397 if (grid
->ncg_ideal
<= 0)
399 gmx_fatal(FARGS
, "The number of cg's per cell should be > 0");
404 fprintf(debug
, "Set ncg_ideal to %d\n", grid
->ncg_ideal
);
410 void done_grid(t_grid
*grid
)
415 sfree(grid
->cell_index
);
419 grid
->cells_nalloc
= 0;
427 fprintf(debug
, "Successfully freed memory for grid pointers.");
431 int xyz2ci_(int nry
, int nrz
, int x
, int y
, int z
)
432 /* Return the cell index */
434 return (nry
*nrz
*x
+nrz
*y
+z
);
437 void ci2xyz(t_grid
*grid
, int i
, int *x
, int *y
, int *z
)
438 /* Return x,y and z from the cell index */
442 range_check_mesg(i
, 0, grid
->nr
, range_warn
);
444 ci
= grid
->cell_index
[i
];
445 *x
= ci
/ (grid
->n
[YY
]*grid
->n
[ZZ
]);
446 ci
-= (*x
)*grid
->n
[YY
]*grid
->n
[ZZ
];
447 *y
= ci
/ grid
->n
[ZZ
];
448 ci
-= (*y
)*grid
->n
[ZZ
];
452 static int ci_not_used(ivec n
)
454 /* Return an improbable value */
455 return xyz2ci(n
[YY
], n
[ZZ
], 3*n
[XX
], 3*n
[YY
], 3*n
[ZZ
]);
458 static void set_grid_ncg(t_grid
*grid
, int ncg
)
463 if (grid
->nr
+1 > grid
->nr_alloc
)
465 nr_old
= grid
->nr_alloc
;
466 grid
->nr_alloc
= over_alloc_dd(grid
->nr
) + 1;
467 srenew(grid
->cell_index
, grid
->nr_alloc
);
468 for (i
= nr_old
; i
< grid
->nr_alloc
; i
++)
470 grid
->cell_index
[i
] = 0;
472 srenew(grid
->a
, grid
->nr_alloc
);
476 void grid_first(FILE *fplog
, t_grid
*grid
,
477 gmx_domdec_t
*dd
, const gmx_ddbox_t
*ddbox
,
478 matrix box
, rvec izones_x0
, rvec izones_x1
,
479 real rlistlong
, real grid_density
)
483 set_grid_sizes(box
, izones_x0
, izones_x1
, rlistlong
, dd
, ddbox
, grid
, grid_density
);
485 grid
->ncells
= grid
->n
[XX
]*grid
->n
[YY
]*grid
->n
[ZZ
];
487 if (grid
->ncells
+1 > grid
->cells_nalloc
)
489 /* Allocate double the size so we have some headroom */
490 grid
->cells_nalloc
= 2*grid
->ncells
;
491 srenew(grid
->nra
, grid
->cells_nalloc
+1);
492 srenew(grid
->index
, grid
->cells_nalloc
+1);
496 fprintf(fplog
, "Grid: %d x %d x %d cells\n",
497 grid
->n
[XX
], grid
->n
[YY
], grid
->n
[ZZ
]);
501 m
= std::max(grid
->n
[XX
], std::max(grid
->n
[YY
], grid
->n
[ZZ
]));
502 if (m
> grid
->dc_nalloc
)
504 /* Allocate with double the initial size for box scaling */
505 grid
->dc_nalloc
= 2*m
;
506 srenew(grid
->dcx2
, grid
->dc_nalloc
);
507 srenew(grid
->dcy2
, grid
->dc_nalloc
);
508 srenew(grid
->dcz2
, grid
->dc_nalloc
);
512 for (i
= 0; (i
< grid
->ncells
); i
++)
518 static void calc_bor(int cg0
, int cg1
, int ncg
, int CG0
[2], int CG1
[2])
538 fprintf(debug
, "calc_bor: cg0=%d, cg1=%d, ncg=%d\n", cg0
, cg1
, ncg
);
539 for (m
= 0; (m
< 2); m
++)
541 fprintf(debug
, "CG0[%d]=%d, CG1[%d]=%d\n", m
, CG0
[m
], m
, CG1
[m
]);
547 void calc_elemnr(t_grid
*grid
, int cg0
, int cg1
, int ncg
)
550 int *cell_index
= grid
->cell_index
;
551 int *nra
= grid
->nra
;
555 ncells
= grid
->ncells
;
558 gmx_fatal(FARGS
, "Number of grid cells is zero. Probably the system and box collapsed.\n");
561 not_used
= ci_not_used(grid
->n
);
563 calc_bor(cg0
, cg1
, ncg
, CG0
, CG1
);
564 for (m
= 0; (m
< 2); m
++)
566 for (i
= CG0
[m
]; (i
< CG1
[m
]); i
++)
571 range_check_mesg(ci
, 0, ncells
, range_warn
);
578 void calc_ptrs(t_grid
*grid
)
580 int *index
= grid
->index
;
581 int *nra
= grid
->nra
;
582 int ix
, iy
, iz
, ci
, nr
;
585 ncells
= grid
->ncells
;
588 gmx_fatal(FARGS
, "Number of grid cells is zero. Probably the system and box collapsed.\n");
592 for (ix
= 0; (ix
< grid
->n
[XX
]); ix
++)
594 for (iy
= 0; (iy
< grid
->n
[YY
]); iy
++)
596 for (iz
= 0; (iz
< grid
->n
[ZZ
]); iz
++, ci
++)
598 range_check_mesg(ci
, 0, ncells
, range_warn
);
608 void grid_last(t_grid
*grid
, int cg0
, int cg1
, int ncg
)
612 int ci
, not_used
, ind
, ncells
;
613 int *cell_index
= grid
->cell_index
;
614 int *nra
= grid
->nra
;
615 int *index
= grid
->index
;
618 ncells
= grid
->ncells
;
621 gmx_fatal(FARGS
, "Number of grid cells is zero. Probably the system and box collapsed.\n");
624 not_used
= ci_not_used(grid
->n
);
626 calc_bor(cg0
, cg1
, ncg
, CG0
, CG1
);
627 for (m
= 0; (m
< 2); m
++)
629 for (i
= CG0
[m
]; (i
< CG1
[m
]); i
++)
634 range_check_mesg(ci
, 0, ncells
, range_warn
);
635 ind
= index
[ci
]+nra
[ci
]++;
636 range_check_mesg(ind
, 0, grid
->nr
, range_warn
);
643 void fill_grid(gmx_domdec_zones_t
*dd_zones
,
644 t_grid
*grid
, int ncg_tot
,
645 int cg0
, int cg1
, rvec cg_cm
[])
650 int zone
, ccg0
, ccg1
, cg
, d
, not_used
;
651 ivec shift0
, useall
, b0
, b1
, ind
;
656 /* We have already filled the grid up to grid->ncg,
657 * continue from there.
662 set_grid_ncg(grid
, ncg_tot
);
664 cell_index
= grid
->cell_index
;
666 /* Initiate cell borders */
669 for (d
= 0; d
< DIM
; d
++)
671 if (grid
->cell_size
[d
] > 0)
673 n_box
[d
] = 1/grid
->cell_size
[d
];
680 copy_rvec(grid
->cell_offset
, offset
);
684 fprintf(debug
, "Filling grid from %d to %d\n", cg0
, cg1
);
688 if (dd_zones
== NULL
)
690 for (cg
= cg0
; cg
< cg1
; cg
++)
692 for (d
= 0; d
< DIM
; d
++)
694 ind
[d
] = static_cast<int>((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.
703 else if (ind
[d
] >= grid
->n
[d
])
705 ind
[d
] = grid
->n
[d
] - 1;
708 cell_index
[cg
] = xyz2ci(nry
, nrz
, ind
[XX
], ind
[YY
], ind
[ZZ
]);
713 for (zone
= 0; zone
< dd_zones
->n
; zone
++)
715 ccg0
= dd_zones
->cg_range
[zone
];
716 ccg1
= dd_zones
->cg_range
[zone
+1];
717 if (ccg1
<= cg0
|| ccg0
>= cg1
)
722 /* Determine the ns grid cell limits for this DD zone */
723 for (d
= 0; d
< DIM
; d
++)
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.
733 if (grid
->ncpddc
[d
] == 0)
743 b1
[d
] = grid
->ncpddc
[d
];
748 b0
[d
] = grid
->ncpddc
[d
];
754 not_used
= ci_not_used(grid
->n
);
756 /* Put all the charge groups of this DD zone on the grid */
757 for (cg
= ccg0
; cg
< ccg1
; cg
++)
759 if (cell_index
[cg
] == -1)
761 /* This cg has moved to another node */
762 cell_index
[cg
] = NSGRID_SIGNAL_MOVED_FAC
*grid
->ncells
;
767 for (d
= 0; d
< DIM
; d
++)
769 ind
[d
] = static_cast<int>((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.
782 else if (ind
[d
] >= b1
[d
])
790 /* Charge groups in this DD zone further away than the cut-off
791 * in direction do not participate in non-bonded interactions.
797 if (cg
> grid
->nr_alloc
)
799 fprintf(stderr
, "WARNING: nra_alloc %d cg0 %d cg1 %d cg %d\n",
800 grid
->nr_alloc
, cg0
, cg1
, cg
);
804 cell_index
[cg
] = xyz2ci(nry
, nrz
, ind
[XX
], ind
[YY
], ind
[ZZ
]);
808 cell_index
[cg
] = not_used
;
817 void check_grid(t_grid
*grid
)
819 int ix
, iy
, iz
, ci
, cci
, nra
;
821 if (grid
->ncells
<= 0)
823 gmx_fatal(FARGS
, "Number of grid cells is zero. Probably the system and box collapsed.\n");
828 for (ix
= 0; (ix
< grid
->n
[XX
]); ix
++)
830 for (iy
= 0; (iy
< grid
->n
[YY
]); iy
++)
832 for (iz
= 0; (iz
< grid
->n
[ZZ
]); iz
++, ci
++)
836 nra
= grid
->index
[ci
]-grid
->index
[cci
];
837 if (nra
!= grid
->nra
[cci
])
839 gmx_fatal(FARGS
, "nra=%d, grid->nra=%d, cci=%d",
840 nra
, grid
->nra
[cci
], cci
);
843 cci
= xyz2ci(grid
->n
[YY
], grid
->n
[ZZ
], ix
, iy
, iz
);
844 range_check_mesg(cci
, 0, grid
->ncells
, range_warn
);
848 gmx_fatal(FARGS
, "ci = %d, cci = %d", ci
, cci
);
855 void print_grid(FILE *log
, t_grid
*grid
)
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
++)
868 fprintf(log
, "%5d %5d\n", i
, grid
->cell_index
[i
]);
870 fprintf(log
, "cells\n");
871 fprintf(log
, " ix iy iz nr index cgs...\n");
873 for (ix
= 0; (ix
< grid
->n
[XX
]); ix
++)
875 for (iy
= 0; (iy
< grid
->n
[YY
]); iy
++)
877 for (iz
= 0; (iz
< grid
->n
[ZZ
]); iz
++, ci
++)
879 index
= grid
->index
[ci
];
881 fprintf(log
, "%3d%3d%3d%5d%5d", ix
, iy
, iz
, nra
, index
);
882 for (i
= 0; (i
< nra
); i
++)
884 fprintf(log
, "%5d", grid
->a
[index
+i
]);