We use the linearized Poisson-Boltzmann equation to construct an analytical estimate of the electrostatic interaction energy for two rigid, globular proteins with arbitrary charge distributions represented by spherical multipole moments. Our analysis avoids invoking the superposition approximation and therefore is not limited to weakly interacting double layers. The interaction energy is a function of the protein properties (radii, multipole moments, and dielectric constant); the geometry (protein separation and orientations); and the solution properties (Debye screening length and dielectric constant). The proteins may be different. This is an extension of previous treatments for constant surface charge or for constant surface potential particles that are more applicable to isotropic colloids. We give example calculations for the interactions of Ribonuclease A molecules interacting in an electrolyte. These results suggest that electrostatic effects may be largely responsible for attractions between proteins, even if they are like charged. A set of routines using MATHEMATICA is available for performing the interaction energy and electrostatic potential calculations.