for small systems but approaches linear scaling (i.e. (N N ) cut ∝ N ) in the limit of large system size. The integrals are made and contracted on-the-fly with the bonding orbitals' pseudo-densities to make bonding-specific exchange integrals K i µν with compute effort scaling as O(O(N N ) 2 cut ). The K i µν matrices are then transformed into the virtual space as K i ab in Eq. 7 with compute cost scaling as O(OV N 2 + OV 2 N ). Asymptotically this is the dominant step in this method unless more careful thresholding is considered. Then, we divide by the appropriate denominator to get T i ab in Eq. 8 (with O(OV 2 ) effort). Lastly, we diagonalize T i for each bonding orbital to get the AB2 antibonding orbitals as in Eq. 9 with O(OV 3 ) effort. Note that the last step can in principle be made O(OV 2 ) since we are only solving for the eigenvector with the largest amplitude in each matrix. We can contrast this procedure with the modified FNO approach which has a dominant computational step that scales as the 5 th power of molecule size: constructing P ab in Eq. 10 with complexity of O(O 2 V 3 ).