ii. System setup for simulation of all molecules: All the simulations were performed using the GROMACS 2019.6 molecular dynamics software patched with plumed 2.6 . Each protein-ligand complex system was solvated by ~76000 TIP3P water molecules in a box with dimensions 16 X 14 X 11 nm . The physiological concentration (150 mM) of Na + and Cl -ions along with extra Cl -ions, were used to neutralize the system. iii. Equilibration and Simulation: Initially, each system was energy minimized using steepest descent for 10000 steps, followed by heating to 300K in 200 ps using Berendsen thermostat and barostat with a coupling constant of 0.6 ps. Restraints of 25 kcal/mol/Å 2 were applied on heavy atoms during the heating process. Thereafter, equilibration was carried out for 2 ns at constant temperature (300 K) and pressure (1 bar) without any restraints using the same thermostat and barostat with coupling constants of 0.2 ps for each. The last 100 ps of NPT simulation was used to calculate the average volume, which was used in the final 5 ns unrestrained NVT equilibration using the Nosé-Hoover thermostat with a coupling constant of 0.2 ps. During the simulation, the LINCS was used to constrain all the bonds, and the Particle Mesh Ewald (PME) method was used for electrostatics. The distance cut-offs for the van der Waals (vdW) and electrostatic long-range interaction were both kept at 10 Å. The timestep for all simulations was taken to be 2 fs. The equilibrated ligand-bound protein structures (systems) were initially simulated for 5 ns. For a particular system, if the ligand was found to be bound after 5 ns simulation, we proceeded with the free energy calculation. Otherwise, the ligand was rejected as being too weak. iv. Free energy calculations. To calculate the binding free energy of the molecules, well-tempered metadynamics simulation was performed after equilibration. We used two collective variables (CVs) distvec (𝑋) and angvec (𝜃) for the well-tempered metadynamics simulations. 𝑋 = 𝑏 ̂. 𝑑 ⃗ and 𝜃 = cos -1 (𝑏 ̂. 𝑑 ⃗ |𝑑|) ⁄ , where 𝑏 ̂ is the unit vector from the Centre of Mass (COM) of residues 127, 130, 131, 134, 139, 140 -147, 157, 162, 175, and 787 to the COM of residues 545, 553, 555, 556, 617-619, 622, 623, 680, 682, 687, 691, 759-761, 811, 813, and 814 that lie more toward the hotspot region (yellow region in the figure ). 𝑑 ⃗ is the vector from the COM of 127, 130, 131, 134, 139, 140 -147, 157, 162, 175, and 787 to the COM of the drug molecule. A hill height of 0.2 kJ/mol and a bias factor of 10 was chosen with a hills deposition rate of 2 ps for the metadynamics simulations. Gaussian widths for 𝑋 and 𝜃 were taken to be 0.6 Å and 0.002 rad, respectively. These CVs have been previously used in DNA intercalation and protein-small molecule binding studies that produced free energy values in good agreement with experiments. An intuitive representation of the CVs is shown in Fig. . An increase in 𝑋 ensures that the molecule moves away from the equilibrated bound structure towards the solution medium. The 𝜃 is designed to capture the directionality of the unbinding of the molecule to complement 𝑋. The minimum free energy path from the ligand-bound minimum to the unbound state in the free-energy landscape was obtained using the program MULE . The free energy convergence for unbinding of the drug in our enhanced sampling simulation takes anywhere between 50 and 200 ns. The metadynamics simulation length for each of the molecules is shown in Table . The 2D representation of small molecules interacting with the protein surface is constructed using ligplot .