  Average atom models are widely used to make equation of state tables and for
calculating other properties of materials over a wide range of conditions, from
zero temperature isolated atom to fully ionized free electron gases. The
numerical challenge of making these density functional theory based models work
for any temperature, density or nuclear species is formidable. Here we present
in detail a hybrid Green's function/orbital based approach that has proved to
be stable and accurate for wide ranging conditions. Algorithmic strategies are
discussed. In particular the decomposition of the electron density into
numerically advantageous parts is presented and a robust and rapid self
consistent field method based on a quasi-Newton algorithm is given. Example
application to the equation of state of lutetium (Z=71) is explored in detail,
including the effect of relativity, finite temperature exchange and
correlation, and a comparison to a less approximate method. The hybrid scheme
is found to be numerically stable and accurate for lutetium over at least 6
orders of magnitude in density and 5 orders of magnitude in temperature.
