We describe a new algorithm for the generation of 3D grids for the numerical evaluation of multicenter molecular integrals in density functional theory. First, we use the nuclear weight functions method of Becke [A. D. Becke, J. Chem. Phys. 88, 2547 (1988)] to decompose a multicenter integral integral F(r) dr into a sum of atomic-like single-center integrals. Then, we apply automatic numerical integration techniques to evaluate each of these atomic-like integrals, so that the total integral is approximated as integral F(r) dr approximate to Sigma(i) omega(i)F(r(i)). The set of abscissas r(i) and weights wi constitutes the 3D grid. The 3D atomic-like integrals are arranged as three successive monodimensional integrals, each of which is computed according to a recently proposed monodimensional automatic numerical integration scheme which is able to determine how many points are needed to achieve a given accuracy. When this monodimensional algorithm is applied to 3D integration, the 3D grids obtained adapt themselves to the shape of the integrand F(r), and have more points in more difficult regions. The function F(r), which, upon numerical integration, yields the 3D grid, is called the generating function of the grid. We have used promolecule densities as generating functions, and have checked that grids generated from promolecule densities are also accurate for other integrands. Our scheme is very reliable in the sense that, given a relative tolerance epsilon, it generates 3D grids which are able to approximate multicenter integrals with relative errors smaller than epsilon for all the molecules tested in this work. Coarser or finer grids can be obtained using greater or small tolerances. For a series of 21 molecules, the average number of points per atom for epsilon=2.0.10(-3), epsilon=2.0.10(-4), epsilon=2.0.10(-5), epsilon=2.0.10(-6), and epsilon=2.0.10(-7) is respectively 3141 (2.9.10(-4)), 10271 (2.4.10(-5)), 27184 (3.1.10(-6)), 72266 (1.9.10(-7)), and 164944 (5.2.10(-9)) (in parentheses are the maximum errors obtained when integrating the density). It is possible to reduce the number of points in the grid by taking advantage of molecular symmetry. It seems that our method achieves a given accuracy with fewer points than other recently proposed methods.