The modeling and simulation of macromolecules in solution often benefits from fast analytical approximations for the electrostatic interactions. In our previous work [G. Sigalov et al., J. Chem. Phys. 122, 094511 (2005)], we proposed a method based on an approximate analytical solution of the linearized Poisson-Boltzmann equation for a sphere. In the current work, we extend the method to biomolecules of arbitrary shape and provide computationally efficient algorithms for estimation of the parameters of the model. This approach, which we tentatively call ALPB here, is tested against the standard numerical Poisson-Boltzmann (NPB) treatment on a set of 579 representative proteins, nucleic acids, and small peptides. The tests are performed across a wide range of solvent/solute dielectrics and at biologically relevant salt concentrations. Over the range of the solvent and solute parameters tested, the systematic deviation (from the NPB reference) of solvation energies computed by ALPB is 0.5-3.5 kcal/mol, which is 5-50 times smaller than that of the conventional generalized Born approximation widely used in this context. At the same time, ALPB is equally computationally efficient. The new model is incorporated into the AMBER molecular modeling package and tested on small proteins.