The fragment molecular orbital (FMO) method is combined with the multipole method (MM) to reduce the quadratic scaling of the electrostatic terms. MM is applied to the energy and gradient of the well-separated interfragment interactions, and the gradient of the one electron electrostatic potential. The accuracy is controlled using the rigorous criteria based on the multipole expansion. Test calculations on ice surfaces using the 6-31G** and 6-31++G** basis sets show that MM is accurate and significantly accelerates FMO, although some quadratically scaling terms remain to be improved. (C) 2012 Elsevier B.V. All rights reserved.