The binding energies obtained by each method were benchmarked against DLPNO-CCSD(T 0 )/augcc-pVTZ calculations with the NormalPNO criterion, carried out on top of the reference ωB97X-D/6-31++G(d,p) geometries. The mean absolute error of these results is presented in Figure , where "reference geometry" refers to geometries optimized at ωB97X-D/6-31++G(d,p), and "optimized geometry" refers to geometries optimized with the same methods utilized for calculating the binding energies. Hence, the "optimized geometry" illustrates both the error in the binding energies as well as changes in the geometry. For the DFT methods with the small 6-31++G(d,p) basis, PW91 performs significantly better than M06-2X and ωB97X-D. However, this is reversed for the larger 6-311++G(3df,3pd) basis set. This could indicate that substantial cancellation of errors is present in the PW91/6-31++G(d,p) calculations. Interestingly, the PBEh-3c and r 2 SCAN-3c methods, which did well in reproducing the geometries, present large errors in the binding energies. As expected, the semi-empirical methods GFN1-xTB and GFN2-xTB exhibited poor performance in calculating binding energies, especially when calculating energy on top of geometries optimized by the same method. Our results indicate that B97-3c represents a good compromise between accuracy and efficiency, and is expected to be suitable for larger systems. If better accuracy is required, ωB97X-D/6-311++G(3df,3pd) or coupled cluster binding energies are needed. However, these methods are extremely expensive for large clusters. These results are consistent with our recent studies and illustrate that the benchmark findings are most likely transferable to other SA-base systems as well.