We present results of theoretical investigation of the stability of small beryllium clusters, Be 2−6 . Due to the fact that the clusters are homonuclear, their stability is sufficiently represented by the binding energy D e calculated per atom. The knowledge of this quantity offers a possibility of calculating more complex energetic effects, such as destruction or fusion of the beryllium clusters and adding or detaching a beryllium atom. The most stable structures among those of the singlet multiplicity were studied using the highly correlated coupled clusters method with iterative single-and double-excitations and the perturbative triple excitations (CCSD(T)). The effect of electron correlation, basis set saturation, basis set superposition error, relativistic effects and the core correlation were carefully investigated and discussed. The optimized virtual orbital space (OVOS) with the controlled accuracy technique is extensively applied to alleviate the computational cost of calculations. In this specific case, the use of the OVOS technique led to speed-ups in computational time by more than an order of magnitude. Finally, we present the basis set extrapolated binding energies with all of the aforementioned corrections properly included, which are the most reliable data for this key property for Be 3−6 clusters to date. Furthermore, these values can serve as a valuable benchmark for testing the accuracy of less time consuming computational methods applicable for larger beryllium clusters.This paper is dedicated to Professor Alexander I Boldyrev on the occasion of his 60th birthday.

No abstract