In this work, we present a series of benchmark variational calculations for the ground and 19 lowest bound excited singlet S and P states of the beryllium atom. The nonrelativistic wave functions of the states that represent the motion of the nucleus and the four electrons around the center of mass of the atom are expanded in terms of up to 17 000 all-particle explicitly correlated Gaussians. The Gaussians are optimized independently for each state. The leading relativistic corrections to the energy levels are computed in the framework of the perturbation theory and they explicitly include the nuclear recoil effects. We also calculate the leading quantum electrodynamics (QED) corrections for each considered state. Using the obtained energy levels and the corresponding wave functions, we compute the transition frequencies, transition dipole moments, and oscillator strengths. A comparison with the available experimental data shows very good agreement. The results of this most comprehensive set of calculations of spectroscopic accuracy for Be to date may open up new applications pertinent to the precision tests of QED, determination of the nuclear charge radius, and modeling matter-radiation equilibria of the beryllium gas that has relevance to the physics of interstellar media.