We used the PBE functional with the all-electron def2-SVP basis sets for all the DFT calculations in the ORCA4.1 . The PBE functional has been used to describe the interactions of phosphorene with a wide range of adsorbates. The DFT-D3(BJ) procedure included dispersion corrections into the PBE functional for energies and gradients . Molecular structures were optimized without geometrical constraints and verified through frequency analyses, where positive vibrations were associated with all vibrational modes. Convergence tolerance values of 1•10 -8 and 1•10 -6 Ha were used for SCF and geometry optimizations; the geometries converged with tolerance values in gradients and coordinate displacements of 3•10 -4 Ha/Bohr and 4•10 -3 Bohr, respectively. Continuum solvent effects were included by the universal continuum solvation model (SMD) based on the quantum mechanical charge density of a solute molecule interacting with a continuum description of the solvent . Excited states were obtained with the simplified timedependent density functional theory (sTD-DFT) approach combined with the meta-hybrid TPSSh functional (10% Hartree-Fock exchange) ; configuration state functions were included up to an energy threshold of 6.2 eV. Phosphorene nanoflakes (P126H30) were used for adsorption studies and with a surface area of at least 1585 Å 2 considering its electron density, which is relatively larger to obtain well-converged adsorption energies concerning the fullerenes surface area (up to 456 Å 2 ). Adsorption energies (Eads) were computed as: