The PPP method was implemented in the Python package coulson which is freely available on GitHub with a permissive open-source MIT license. Further details on coulson will be reported elsewhere. The PPP wavefunction was converged with the self-consistent field method using the naive variational principle. CIS calculations used the converged SCF wavefunction as reference and modelled vertical excitations based on the pre-computed geometries. To make the calculations more robust, we derived a spline interpolation of the overlap integrals, which will be reported separately. Evaluated on the azaphenalene dataset (vide infra), it shows good accuracy with 𝑅 2 = 1.000 and RMSE = 0.004 eV for Δ𝐸 DSP ST,CIS (Figure ). All calculations have been performed on a MacBook Air laptop computer with an M2 processor. We further integrated coulson with PySCF 55 to provide alternative algorithms for the SCF convergence and CIS and found the results consistent within numerical accuracy. To give an indication of the computational cost, the SCF calculation for azaphenalene (13 heavy atoms, 14 electrons), takes 3.63±0.37 ms (average and standard deviation over 100 runs, respectively). On top of this, DSP takes 0.20±0.04 ms, and CIS + DSP takes 11.36±0.73 ms. Oscillator strengths were calculated with the dipole length approximation. Some compounds exhibit negative triplet excitation energies and in some case even negative singlet excitation energies, which indicates a restricted-unrestricted instability of the ground state wavefunction. We have here taken the pragmatic approach to ignore these issues, while we will highlight below some examples where it occurs and if it has any effect on the overall results. This approach can be partly justified as the reference data also does not include any stability analysis. The results should anyway be indicative of the gap between the lowest singlet and triplet states with open-shell character, regardless of if they are lower in energy than the closed-shell singlet.