  In this work we propose a novel composite method for accurate calculation of
the energies of many-electron atoms. The dominant contribution to the energy
(pair energies) are calculated by using explicitly correlated factorisable
coupled cluster theory. Instead of the usual Gaussian-type geminals for the
expansion of the pair functions, we employ two-electron Hylleraas basis set.
This eliminates the need for massive optimisation of nonlinear parameters and
the required three-electron integrals can now be calculated relatively easily.
The remaining contributions to the energy are calculated within the algebraic
approximation by using large one-electron basis sets composed of Slater-type
orbitals. The method is tested for the beryllium atom where the accuracy better
than $1\,$cm$^{-1}$ is obtained. We discuss in details possible sources of the
error and estimate the uncertainty in each energy component. Finally, we
consider possible strategies to improve the accuracy of the method by one to
two orders of magnitude. The most important advantage of the method is that is
does not suffer from an exponential growth of the computational costs with
increasing number of electrons in the system and thus can be applied to heavier
atoms preserving a similar level of accuracy.
