We fully agree with Spackman’s tribute to the considerable contribution of R. F. Stewart to the field of charge-density analysis and will restrict our comments to the second half of his letter which discusses a number of issues relevant to the results presented in our paper (Volkov, King et al., 2006). Mark Spackman (MS) states that ‘it is worth noting that the map presented by VKCF is based on a direct-space summation of contributions from only the eight nearest neighbours in the crystal’, implying that the figure presented will change significantly if more neighbouring molecules are included in the calculation. In Fig. 1 below, we present analogous maps obtained by taking into account the contributions of pseudoatoms within several unit cells (this option was available in the latest version of XD2006 at the time our original paper was published). Only small changes are observed on the periphery of the map, which do not in any way change any of the conclusions presented in our paper. For the calculation of the electrostatic potential (ESP), electric field (EF) and electric field gradient (EFG), we use the same EPMM approach (Volkov et al., 2004) developed for calculation of the electrostatic interaction energy. Exact formulae are used for atoms located within a certain distance (typically ~6 A for the ESP, EF and EFG in organic systems) from the point at which the properties are calculated. Beyond that distance, approximate formulae, based on the expansion of the ESP, EF and EFG in atomic multipole moments, are used. Table 1 includes a comparison between exact and EPMM methods, demonstrating the satisfactory accuracy of the much faster (Table 2) EPMM method. Table 3 shows the convergence of the x and y components of the EF vector in the plane of the map versus the number of unit cells included in the calculation. An excellent convergence is achieved when including atoms with fractional coordinates of 1 < x, y, z < 2, while the inclusion of eight neighbouring cells is also quite satisfactory, i.e. the r.m.s. deviation is only 0.02–0.03 e A , which is less than 1% of the r.m.s. values of the Ex (~5.4 e A ) and Ey (~6.5 e A ) components of the EF vector in the plane of the map. Note that a similar direct-space approach is used for the calculation of the electrostatic binding energy in the crystal (Volkov et al., 2007), which has also been implemented in XD2006 (Volkov, Macchi et al., 2006). The efficient combination of the exact formulae and the multipole moment approximation eliminates the convergence problems inferred by MS. MS quotes Cummins et al. (1976) that ‘computation of the EF in the crystal due to point dipoles is well known to require lattice summation techniques to achieve convergence’ and adds that ‘these combined approaches are entirely analogous to the Ewald approach used to achieve rapid convergence of lattice sums.’ We are of course well aware of the existence of the Ewald summation techniques which we used previously in the calculation of the lattice energy of ionic (Su & Coppens, 1995) and molecular (Abramov, Volkov, Wu & Coppens, 2000a,b; Abramov, Volkov & Coppens, 2000) crystals. However, the direct-space summation was found to be advantageous in terms of both speed and ease of implementation when combined with the new EPMM method (Volkov et al., 2007). MS states that we misunderstand the purpose of the methods described by Brown & Spackman (1994) as several algorithms published in this paper were used only for debugging of the VALRAY code and should not be considered of practical importance. Table 1 Root-mean-squared differences in electric field components Ex and Ey (e A ) in the plane of the map between exact and EPMM calculations .