This work develops a multiscale modeling framework for defects in crystals with general geometries and boundary conditions in which ionic interactions are important, with potential application to, e.g., ionic solids and electric field interactions with materials. The overall strategy is posed in the framework of the Quasicontinuum multiscale method; specifically, the use of a finite-element inspired kinematic description enables a significant reduction in the large number of degrees of freedom to describe the atomic positions. The key advance of this work is a method for the efficient and accurate treatment of nonlocal electrostatic charge-charge interactions without restrictions on the geometry or boundary conditions. Electrostatic interactions are long-range with slow decay, and hence require consideration of all pairs of charges making a brute-force approach computationally prohibitive. The method proposed here accounts for the exact charge-charge interactions in the near-field and uses a coarse-grained approximation in the far-field. The coarse-grained approximation and the associated errors are rigorously derived based on the limit of a finite body with a small periodic lengthscale, thereby enabling the errors in the approximation to be controlled to a desired tolerance. The method is applied to a simple model of Gallium Nitride and it is shown that electrostatic interactions can be approximated with a desired level of accuracy using the proposed methodology.