Black phosphorus serves as an exemplary stacked bidimensional semiconductor, exhibiting anisotropic features in electronic and optical properties that demand special attention in theoretical investigations. Herein, we employed a series of computational protocols, starting with first-principles approaches (particularly density functional theory�DFT), combined with the solution of the Bethe−Salpeter equation within the tight-binding method to explore the structural stability and optoelectronic properties (bandgap, exciton binding energies, and optical absorption) of black phosphorus (P n ) across layers ranging from n = 1 to 6 and bulk. In our DFT investigations, we observed that empirical and semiempirical van der Waals models, contributing a dispersion energy component, revealed a myriad of differences and similarities in properties, such as interlayer nonbonded interactions. Notably, the many-body dispersion correction exhibited superior performance in connecting layered systems with the bulk. The magnitude of dispersion energies correlated with the stability during the aggregation process P (n−1) + P 1 → P n . Additionally, the bandgap, properly corrected through relativistic quasi-particle calculations, narrowed due to enhanced interlayer wave function overlap, a result of the dispersion energies promoting the shortening of interlayer distances. Subsequently, we utilized the band structure relativistically corrected as a starting point to obtain the Hamiltonian, achieved through the generation of maximally localized Wannier functions. This facilitated a screening of the electron−hole (e−h) pairwise interaction Coulomb potential, specifically the exciton binding energy. We identified an indirect impact of the dispersion energies on excitonic properties, which were effectively described by the Rytova−Keldysh model for the e−h Coulomb potential, aligning well with photoluminescence experiments.