Briefly, we used two stages of minimization. Firstly, each complex was optimized by 500 steps of steepest descent followed by another 500 steps of conjugate gradient minimization, keeping all atoms of the complex restrained to their initial coordinate with a weak harmonic potential (force constant 2.0 kcal mol−1 Å−2). The second stage of minimization was carried out without any restraints by performing 100 steps of steepest descent, followed by another 900 steps of conjugate gradient minimization to remove any residual steric clashes. Subsequently, all systems were gradually heated from 0 K to 300 K at the NVT ensemble with a force constant of 2.0 kcal mol−1Å−2 acting on all solute atoms. Next, 1.0 ns equilibration MD phase was carried out without any restraint. Finally, we performed 150 ns production simulations for all four systems at the NPT ensemble. Overall, we accumulated 15000 conformations for each simulation, and we used 500 snapshots from the last 50 ns trajectories for binding affinity estimation using the MM-PBSA scheme. The trajectory analysis was done by the AmberTools19 CPPTRAJ (Roe & Cheatham III, 2013) module of Amber18.