Numerical integration of the exchange-correlation (XC) contribution to the Fock matrix can be efficiently carried out on GPUs. In this work, we follow the approach proposed by Ochsenfeld and co-workers, where the evaluation of densities on the grid points and XC contribution to the Fock matrix are done by matrix multiplication which can make efficient use of the GPUs. Moreover, the O(N 3 ) computational cost in XC integration, where N is the system size, can be greatly reduced by integrating over even-sized grid batches, such that the computational cost becomes asymptotically linear due to efficient screening of basis functions. The key to such an efficient implementation is to make sure that the DFT grid points of the whole system are divided as evenly as possible into batches of spatially adjacent grid points. In this work, this is realized by repeatedly bisecting the box of grid points into smaller and smaller boxes, until the number of grid points in each box is below a given threshold. The bisection of boxes is done as evenly as possible (in terms of number of grid points), such that the final boxes of grid points have almost constant computational cost on average. Through matrix multiplication the majority of the computation steps of XC numerical integration is done on GPUs, while the evaluation of functional derivatives, for which we use the Libxc library, is done on the node CPU.