In this work we report progress in the development and implementation of quantum-mechanical methods for calculating bound ground and excited states of small atomic systems. The work concerns singlet states with the L = 1 total orbital angular momentum (P states). The method is based on the finite-nuclear-mass (non-Born-Oppenheimer; non-BO) approach and the use of all-particle explicitly correlated Gaussian functions for expanding the nonrelativistic wave function of the system. The development presented here includes derivation and implementation of algorithms for calculating the leading relativistic corrections for singlet states. The corrections are determined in the framework of the perturbation theory as expectation values of the corresponding effective operators using the non-BO wave functions. The method is tested in the calculations of the ten lowest The determination of atomic energy levels and transition frequencies with spectroscopic precision (i.e., well below 1 cm −1 ) as well as the corresponding wave functions remains one of the most technically and computationally challenging tasks for the quantum theory of atoms and molecules. The amount of computations required grows very rapidly with each additional electron. As a result, retaining accuracy at manageable computational cost becomes a difficult task even for fewelectron systems. Many accurate methods of various natures have been developed in the electronic structure theory in the past several decades. However, in practical calculations most of them can reach only chemical accuracy (1 kcal/mol) and often cannot be applied to excited states. Yet some emerging physical problems related to precision measurements, atomic clocks, and high-resolution spectroscopy require very accurate quantum-mechanical calculations that exceed chemical accuracy by several orders of magnitude.Precise determination of the atomic energies, transition frequencies, and other basic properties requires precise calculation of the effects resulting not only from the strong Coulombic interactions between the particles but also from more subtle effects due to relativism, quantum electrodynamics (QED), and finite nuclear mass and size. There has been a steady progress Kramida and Martin [6]. Recently, ECGs were used to calculate the lowest S → P transitions of beryllium [7]. The calculations were performed with the infinite-nuclear-mass (INM) approach and the finite-mass effects were obtained using the perturbation theory. In that work, the leading relativistic corrections for the 2 1 S, 3 1 S, and 2 1 P states of beryllium were calculated using the INM wave functions. We should also mentioned our recent high-accuracy calculations of the lowest four 2 S states of the boron atom [8] where similar accuracy as achieved before for four-electron atoms was reached. In those calculations the finite-nuclearmass (FNM) approach was used. Thus, the finite-mass effects were directly incorporated in the nonrelativistic total energy, as well as in the relativistic corrections which were calcul...