In this work, we conducted three independent MFES for each molecule. Each MFES was performed with a ladder of 11 equidistant distributions (λ = 0.0, 0.1 . . . 0.9, 1.0) connecting the reference and target potential (two end states and nine bridging distributions). 5,000 samples were generated from each distribution by performing 5 ns simulations, and writing coordinates every 1 ps. Simulations were performed using Langevin dynamics, with a collision rate of 10 ps - and a time step of 1 fs. All simulations were performed in vacuum, and initial velocities were drawn from a Maxwell-Boltzmann distribution at the simulation temperature of 300 K. To ensure that our samples were drawn from the equilibrium Boltzmann distribution, we discarded the first 20% of each simulation (resulting in 4,000 samples per simulation) and selected only every 5th sample for the equilibrium free energy estimates. The resulting dataset contained 11×800 equidistant samples for the 11 alchemical states. The free energy difference between the MM and NNP potential was then calculated using the MBAR (multistate Bennett acceptance ratio) method as implemented in the pymbar package. The free energy estimate is calculated as the average of three independent simulations.