The development of quantitative models for understanding physical properties of alloys requires a proper treatment of magnetic interactions, which is of paramount importance for the microstructural stability, especially in steels and high-entropy alloys containing magnetic elements. These magnetic interactions also control the defects behavior which affects the mechanical properties and the response under irradiation. Current interatomic potentials for molecular dynamics (MD) simulations still lack an adequate formulation to include magnetism into the simulations. In this paper, the universal equation of states (UES) is revisited and generalized by including ferromagnetic (FM) and antiferromagnetic (AFM) configurations with the aim of proposing a new formulation to develop interatomic potentials with magnetic contribution. For the case of Fe, given a fixed magnetic configuration and magnitude of the magnetic moment, the energy of the system is calculated by means of three parameters, namely, the energy, volume, and corresponding scaling volume (directly related to the bulk modulus) at the local ground state of the corresponding lattice. These parameters depend on three terms: firstly, a distance-dependent function, which gathers the nonmagnetic influence of the surrounding atoms; secondly, a magnetically dependent function, contributing to the energy by means of the magnetic nature of the atom, irrespective of the magnetic moment magnitudes of its surrounding atoms; and finally, a term which is magnetically and distance dependent simultaneously, which describes the influence of the magnetic state of the surrounding atoms on the energy considering their interatomic distance. This latter term is built via two functions which cannot be disconnected: one dependent on the distance between two atoms (a decreasing function with the distance in absolute value) multiplied by another function which is dependent on the magnetic moment of these two atoms. In this way, the magnetic influence of a distant atom scales with its distance. The new formulation is tested for magnetic iron, where 18 240 spin polarized density functional theory (DFT) calculated energies for different lattices, volumes, and magnetic moments in FM and AFM configurations showed that the generalized UES (GUES) accurately describes the energy of the system. The root-mean-square error of the GUES is in the range of $5.9\ifmmode\times\else\texttimes\fi{}{10}^{\ensuremath{-}3}\phantom{\rule{0.16em}{0ex}}\mathrm{eV}$ over all DFT calculated energies, showing good accuracy and allowing us to propose a formulation for developing a magnetic interatomic potential (MIP) in Fe. The potential is developed for FM configuration in iron, aiming at studying the stability of the ferritic and austenitic phases but also defects and other configurations of special relevance as, for instance, in irradiation conditions for fusion or fission applications. The distance and magnetic functions of the GUES are tabulated to obtain a MIP, which describes the DFT calculated energies for different lattices, volumes, and magnetic moments in FM configuration. Further, the MIP has been validated in other crystal lattices (A15 and C15), elastic constants, stresses in the lattice, vacancies, interstitials, forces at different temperatures, transformation paths between body-centered cubic, face-centered cubic, hexagonal close-packed, and simple cubic structures as well as $\ensuremath{\gamma}$ surfaces. We conclude that the results in this paper pave the way to develop MIPs with accuracy and predictability beyond the state of the art in MD.