It should be noted that PBE, based on GGA, offers a balanced exchange-correlation description, making it more reliable for adsorption studies. It better captures long-range dispersion interactions than PBESOL, which overbinds in bulk systems. PBE ensures accurate adsorption energies, charge transfer, and electronic modifications, making it ideal for studying gas and water adsorption on 2D materials. A plane wave basis set was employed with a kinetic energy cutoff of 50 Ry for electron density and 500 Ry for charge density, ensuring precision in structural optimizations. Brillouin zone (BZ) sampling used an 8X8X1 Monkhorst-Pack k-point grid for geometry optimization and a denser 10X10X1 grid for electronic structure calculations. The system's atomic positions and cell parameters were fully relaxed, achieving energy convergence within 10 -10 eV. To accurately model van der Waals interactions, DFT-D3 dispersion corrections were applied. Furthermore, to address the PBE functional's known bandgap underestimation, the HSE functional was employed with a reduced k-point grid of 4X4X1. The study's unit cell included configurations of one, two and three BN atom pairs doped within the qTP-C 60 sheet. A 25 Å vacuum region was introduced along the z-direction to eliminate interactions between periodic images, ensuring accurate adsorption energy calculations and electronic property evaluations. The QE input files of all four systems are provided in SI. The formation energy (E f or ) was determined using Equation .