The 2D binding energy landscapes were used to fit the potential. To validate the parameterization method, we calculated the binding energy profile along a path perpendicular to the molecular plane using the optimized parameters. This path is aligned along the z-direction and centered at (x min , y min ), where x min and y min represent the position of the minimum binding energy in the x-y plane. The point (x min , y min ) was determined from the 2D binding energy landscape. The distance z was systematically varied from z = 0 to z = 6 Å in increments of 0.1 Å. At each point along this path, the binding energy was calculated using DFT and MS. For the calculations in the MS framework, both corrected and uncorrected force fields were considered. For all the systems, we estimated the value of the binding energy, ∆E(z min ), by identifying the global minimum of the binding energy profile at z min .