A procedure has been proposed to construct numerically the exchange-correlation epsilon(XC)(r) and correlation epsilon(C)(r) energy densities of density functional theory using the correlated first- and second-order density matrices from ab initio calculations, epsilon(C)(r) as well as its kinetic and potential components have been obtained for the two-electron He atom acid H-2 molecule. The way various correlation effects manifest themselves in the form of epsilon(C)(r) has been studied. The epsilon(C)(r) have been compared with some density functional local and gradient-corrected models epsilon(C)(mod)(r). The investigation of the shape of the model energy densities epsilon(C)(mod)(r) has been extended to the Be-2 and F-2 molecules and the corresponding correlation energies E(C) have been calculated and discussed for a number of atomic and molecular systems. The results show the importance of a proper modeling of E(C)(r) in the molecular bond midpoint region.