the simplest COHSEX self-energy of Eq. (52) and COHSEX correlation energy of Eq. (51) contain the non-local density matrix ρ$_{0}$ that projects the action of the screened Coulomb interaction W onto the occupied states. If the occupied states are derived from localized states such as 3 d or 4 f orbitals of transition metals, then the projection is automatically onto this manifold. Namely, already at the COHSEX level, we expect to recover the benefits of an LDA+U type treatment. This is no surprise since a static and localized approximation to the GW self-energy yields the kernel of the LDA+U approach [7]. As an added benefit, a COHSEX type approach should automatically include a properly screened U parame- ter, as opposed to U being an externally chosen parameter in the usual LDA+U method. The dynamic formulae of Eqs. (44) and (45) extend these results to include the frequency dependence of the screening. What is missing from these approximate results are the con- tributions to the RPA correlation energy which are physically distinct from those stemming from screened interband transitions: these are the neglected residue contributions in Eq. (43) from the poles from the screening modes themselves. Unfortunately, at present we are un- able to provide simple physical examples of situations where these neglected contributions play the dominant role and our key approximation to fail for basic physical reasons. This question is a subject of present investigation. In this section, we describe numerical results on atomic systems. The aim of this section is not to present an exhaustive and comprehensive treatment. That is the subject of a future investigation. Rather, the main aim is to demonstrate the numerical efficacy of the plasmon formula Eq. (36) for the GW -RPA correlation energy and to test the main approximation of high frequency screening used in the approximate results of Eqs. (44) and (45). Our atomic code uses a standard non-relativistic approach. The atomic eigenfunctions are assumed to take a spherical form: for an eigenstate with spin index σ = ± 1, the spatial part is R$_{nlσ}$ ( r ) Y$_{lm}$ ( θ, φ ). The radial part R$_{nl}$ ( r ) is represented numerically on a radial grid of exponentially-spaced points [63, 64]. The spherical harmonics Y$_{lm}$ are included up to a some maximum angular momentum l$_{max}$ , typically l$_{max}$ = 3 below. When building up the one-particle density matrix ρ$_{σ}$ for spin σ , we allow for m -dependent state fillings f$_{nlmσ}$ ∈ { 0 , 1 } . The local spin-density approximation (LSDA) [62] or unrestricted Hartree-Fock (HF) [65] equations are solved by minimizing the appropriate total energy over arbitrary occupied orthonormal radial functions. We have tested our code with available high quality LSDA [64] and HF [66] data to ensure agreement to at least one part in 10 8 in total energies. When computing various correlation energies below, we make a spherical approximation which is consistent with assuming spherical eigen- functions. This means that when we compute the polarizability ( P or P$_{t}$ ) and the screened interaction ( W or W$_{t}$ ), we assume the fillings have no m -dependence ( i.e. , shell occupancies f$_{nlσ}$ instead of orbital occupancies f$_{nlmσ}$ ): different m -channels will not interact so that taking both l and m as good quantum numbers is a self-consistent assumption. However, ρ$_{σ}$ ( r, θ, φ, r $^{′}$, θ $^{′}$, φ$^{′}$) = ∑ n,l,m f$_{nlmσ}$R$_{nlσ}$ ( r ) R$_{nlσ}$ ( r $^{′}$) Y$_{lm}$ ( θ, φ ) Y$_{lm}$ ( θ $^{′}$, φ$^{′}$)$^{∗}$, 11. NUMERICAL TESTS: ATOMS 28