A method is developed to compute the electrostatic interaction energy of charged colloidal particles in a concentrated aqueous dispersion. The linearized Poisson-Boltzmann equation is used for the potential distribution with the assumption of low particle potential. The boundary element method is used to solve the integral representation of the linearized Poisson-Boltzmann equation under the condition of constant surface potential or constant surface charge density. The method is successfully tested for the system of two spherical particles for which the analytical solution is available. We apply the boundary element method to the systems of three spherical particles, FCC lattice, and BCC lattices. The results show that the pairwise additivity assumption is not good for strongly interacting particles. In addition, by comparing the electrostatic interaction energies of the FCC and BCC lattices at given concentrations, we may explain the possibility of the phase transition between the two phases.