We present a new divide-and-conquer algorithm to efficiently evaluate the Coulomb interaction (which can be extended to other convolution terms) in large systems, which is an essential part of self-consistent first-principle calculations.The total Coulomb potential V(r) =l/|r| is divided into a short range part Vs(r) and a smooth long range part VL(r).The system is divided into many cuboids, with a small box defined for each cuboid plus a buffer region.For the short range part, the interaction convolution integral is calculated directly using a Fast Fourier Transformation (FFT) in the small box.For the smooth long range part, the convolution integral is evaluated by a global FFT but on a coarse grid.The conversion between the dense grid and coarse grid values is done using small box FFTs with proper mask functions.Using this small box FFT method, the total Coulomb potentials of test quantum dot systems on 4803 grid (Fig.1) and 24003 grid are calculated.For the 24003 grid case, the calculation is carried out by tens of thousands of processors with a computational speed up close to 10 times when compared with direct global FFT calculations using the FFTW package with the maximally allowed number of processors.The maximum relative error is 4 × 10-5 while the absolute error is less than 0.1 meV (Fig.1).