Numerical simulation of diffraction on thick holographic gratings in shift-multiplexed optical data storage application is presented. The grating is generated by the interference of a spherical reference wave and a plane signal wave corresponding to a single pixel of the input data page. To describe diffraction on this weak-index-modulated grating, we use the volume integral equation in the first Born approximation. This description yields a convolution integral that can be efficiently evaluated by a 3D fast Fourier transform (FFT) technique. For a $51.2\mathrm{\mu }\mathrm{m}$ recording layer thickness, a serial-divided single personal computer code was built based on parallel FFT coding principles. Diffracted electric field and Poynting-vector distributions are calculated for probe beams spatially shifted with respect to the reference beams. The shift selectivity curves show significant differences from previous analytical calculations based on paraxial propagation and infinite gratings, as they have monotonic decrease in all three directions instead of sinclike functions with Bragg nulls. With the chosen numerical aperture of 0.6 and linear polarization, both the scalar and vector calculations provided similar results within 5%.