  A novel energy minimization formulation of electrostatics that allows
computation of the electrostatic energy and forces to any desired accuracy in a
system with arbitrary dielectric properties is presented. An integral equation
for the scalar charge density is derived from an energy functional of the
polarization vector field. This energy functional represents the true energy of
the system even in non-equilibrium states. Arbitrary accuracy is achieved by
solving the integral equation for the charge density via a series expansion in
terms of the equation's kernel, which depends only on the geometry of the
dielectrics. The streamlined formalism operates with volume charge
distributions only, not resorting to introducing surface charges by hand.
Therefore, it can be applied to any spatial variation of the dielectric
susceptibility, which is of particular importance in applications to
biomolecular systems. The simplicity of application of the formalism to real
problems is shown with analytical and numerical examples.
