Protein ligand binding prediction typically relies on docking methodologies and associated scoring functions to propose the binding mode of a ligand in a biological target. Significant challenges are associated with this approach, including the flexibility of the protein-ligand system, solvent-mediated interactions and associated entropy changes. In addition, scoring functions are only weakly accurate due to the short time required for calculating enthalpic and entropic binding interactions. The workflow described here attempts to address these limitations by combining Supervised Molecular Dynamics (SuMD) with Dynamical Averaging Quantum Mechanics Fragment Molecular Orbital (DA-QM-FMO). This is illustrated using a set of five ligands targeting the SARS-CoV-2 Papain-like protease protein. This combination significantly increased the ability to predict the experimental binding structure of proteinligand complexes independent from the starting position of the ligands or the binding site conformation. We found that the predictive power could be enhanced by combining the residence time (SuMD) and interaction energies (DA-QM-FMO) as descriptors in a novel scoring function named the P-score.