In practice, we calculate this integral numerically over an interval centered on CN A with width 4/ √ β 2 . We determined that 64 quadrature points are sufficient for a smooth behavior of C 6 for all elements, as a lower number of quadrature points create artificial noise patterns with high frequency in elements with strongly oscillating values for the reference polarizabilites. This was especially noted in carbon and nitrogen. We additionally note that a lower value of β 2 leads to smoother changes in the C 6 coefficients with respect to coordination number alterations. While it does introduce a deviation of the calculated C 6 coefficients from the reference values, we note that this behavior is also present in standard D4, and it was not found to be detrimental to energy quality. We do not explicitly discuss the calculation and implementation of dispersion gradients and higher derivatives here, but their analytical calculation from the proposed equations is straightforward and was done analogously to the derivation of gradients in standard D4. In the case of D4SL, we discourage the use of numerical gradients, as the use of numerical integration introduces tiny ripples into the PES, which are strongly amplified by numerical differentiation.