A second-order perturbation (2PT) approach to the spin-orbit interaction (SOI) is implemented within a density-functional theory framework. Its performance is examined by applying it to the calculation of the magnetocrystalline anisotropy energies (MAE) of benchmark systems, and its efficiency and accuracy are compared with the popular force theorem method. The case studies are tetragonal FeMe alloys (Me=Co, Cu, Pd, Pt, Au), as well as FeMe (Me=Co, Pt) bilayers with (111) and (100) symmetry, which cover a wide range of SOI strength and electronic band structures. The 2PT approach is found to provide a very accurate description for 3d and 4d metals and, moreover, this methodology is robust enough to predict easy axis switching under doping conditions. In all cases, the details of the bandstructure, including states far from the Fermi level, are responsible for the finally observed MAE value, sometimes overruling the effect of the SOI strength. From a technical point of view, it is confirmed that accuracy in the MAE calculations is subject to the accuracy of the Fermi level determination.The interplay of MCA and dimensionality is considered as a promising route for the development of spintronics. Furthermore, MCA plays a central role in the magnetic properties of two-dimensional materials, since it can prevent thermal fluctuations from destroying long-range magnetic order [8]. In this context, density functional theory (DFT) calculations that include SOI constitute a powerful theoretical tool to characterize the MCA. Nonetheless, the small magnitude of the associated magnetocrystalline anisotropy energy (MAE) imposes stringent convergence to the DFT calculations at the expense of high computational demands. In practice, fully relativistic DFT calculations that include SOI are carried out by first computing self-consistently the Hamiltonian of the system including the scalar relativistic term and next adding the spin-orbit contribution. The latter may be treated self-consistently (FSC) or non-self-consistently [9, 10] within the so-called force theorem (FT) [11] to obtain the MAE. FT is often considered to produce a good estimate of the FSC MAE value for every system, at least for bulk crystals, although it has been shown to be less accurate in the case of low dimensional systems [12]. In any case, this comparison between FT and SCF calculations can only be done for relatively simple systems, because the FSC tends to be computationally too demanding. Alternatively, Green's function methods that treat SOI at the level of second-order perturbation theory (2PT) have been formulated in the literature under different flavours [13][14][15][16][17]. These approaches rely on the knowledge of the spin-orbit effects on atomic orbitals (AOs) [18][19][20], which can be quasi-analytically accounted. However, to our knowledge, the versatility of those perturbative methods has not been exploited within a DFT environment. Considering the high computational expense involved in DFT calculations that include SOI under the aforemention...