Free energy perturbation (FEP) based on molecular mechanics (MM) simulations can be an effective guide during the hit-to-lead stages of a drug design campaign as it provides a formally rigorous means to compute protein-ligand binding free energies. In practice, the predictive ability of such simulations is effectively limited by two major factors, 1) the accuracy of the underlying MM force field that is used for the rapid calculation of the system energy, and 2) finite simulation times that can limit the conformational space explored. In the expectation of making such calculations routinely reliable, the development of enhanced sampling methods is an active area of research, yet virtually all FEP simulations employ transferable biological force fields, such as AMBER, OPLS, GROMOS and CHARMM, all with quite similar functional forms and parameter fitting strategies. These biological force fields, alongside their small molecule counterparts, have had wide success to-date thanks to meticulous fitting of parameters to reproduce quantum mechanical (QM) and experimental properties of sets of small organic molecules. However, there is room for improvement. It is widely acknowledged that atomic point charges are sensitive to their (local and long-ranged) environment, which is why small molecule force fields typically employ atomic charges that are fit to the molecular electrostatic properties (e.g. ESP or CMx charges ), on a case-by-case basis. Interestingly, this leads to a disconnect between protein and small molecule force fields, in which the former sets of atomic charges are read from a transferable library, and the latter are derived using methods that are not always consistent with the underlying biological force field. Also, standard libraries of parameters describing torsional rotation about flexible bonds are often blamed for observation of unphysical conformations in MM simulations, and these parameters are often re-derived specifically for the molecule under study. With regards to these issues, there has been recent interest in molecule-specific, or bespoke, force fields in which the parameters that govern the dynamics of the system are not assigned from a library based on predetermined atom types, but instead are inferred directly from QM calculations specifically for the molecule of interest. One such example is the QUantum mechanical BEspoke (QUBE) force field, which has a particular focus on scalability to large system sizes and applications in the condensed phase. The QUBE force field shares its functional form with OPLS, so that it retains the favorable computational efficiency of transferable force fields, but differs in that as many parameters as feasible are derived directly from routine, molecule-specific QM calculations. The ground state electron density of the molecule under study is first computed in a weak implicit solvent to simulate the effect of environmental polarization. The density is then partitioned into a set of approximately spherical atom-centered basins via the density derived electrostatic and chemical (DDEC) atoms-in-molecule (AIM) approach, from which we compute the environment-specific non-bonded parameters, including (atom-centered and off site) atomic charges and Lennard-Jones parameters. Since the DDEC method is implemented in the linear-scaling density functional theory code, ONETEP, we can derive these parameters consistently for both small molecules and also systems comprising thousands of atoms, such as proteins. QUBE bond and angle force field parameters are derived directly from the QM Hessian matrix of small molecules, as described previously, and flexible torsions may be parametrized by fitting to constrained one-dimensional QM dihedral scans. Parameter assignment is automated by the QUBEKit software package. To date, the first generation of the QUBE force field has undergone extensive benchmarking against established performance metrics, such as the prediction of the condensed phase thermodynamic properties (density, heat of vaporization and free energy of hydration) of over 100 small organic molecules. A custom library of bonded parameters for protein simulations has been developed and validated via the comparison of molecular dynamics trajectories with NMR observables. In all of these cases, QUBE performed to a similar standard as established and optimized transferable force fields. In the context of FEP calculations, QUBE has been applied to the study of the benchmark L99A mutant of T4 lysozyme, achieving a mean unsigned error (MUE) of 0.85 kcal/mol in the prediction of the absolute binding free energies of six benzene derivatives. However, typical hit-to-lead studies in drug discovery scenarios are significantly more complex than the above study in terms of the sizes of the ligands, the nature of their interactions, and their conformational flexibility. In this letter, we therefore retrospectively calculate the relative binding free energies of a series of 17 drug-like inhibitors of p38α MAP kinase (Figure ). This represents a typical optimization scenario involving both polar and non-polar substitutions around a benzene ring, with activities that span 2-3 orders of magnitude (Table ). As we shall discuss, the binding pose is determined to a large extent by two flexible dihedral angles (φ 1 and φ 2 , Figure ), which impose complex sampling requirements on the simulations. This set of transformations has been the target for a range of activity prediction methods including FEP calculations, which were used to demonstrate the importance of the initial water placement during Monte Carlo (MC) simulations using the OPLS force field. 29 Here, we begin by parametrizing molecule-specific force fields for the 17 p38 kinase inhibitors (plus compound 18, which does not have experimental data for comparison but is a useful FEP intermediate) using the QUBEKit software. Namely, non-bonded (charge and Lennard-Jones) parameters are derived using atoms-in-molecule partitioning of the ground state electron density as described previously. Parametrization of the protein non-bonded parameters is performed using the same protocols, while bonded parameters are read in from Bond and angle parameters of the small molecules (1-18) are derived using the modified Seminario method computed using the QM Hessian matrix at the optimized geometries. Finally, parameters describing rotation about the two flexible dihedral angles φ 1 and φ 2 are fit to constrained QM potential energy scans. Figure shows the results of the torsion parameter optimization for ligand 1. The fit to the underlying QM data is very good with an average root mean square deviation between sampled QM and QUBE torsional scans of 0.07 kcal/mol. For comparison, typical errors in excess of 1.5 kcal/mol are observed using small molecule transferable force fields. By deriving the QUBE force field directly from QM, our goal is to provide accurate and automated molecule-specific parameters that reproduce as closely as possible the full QM potential energy surface. Figure shows the correlation between QUBE and QM relative energies of structures 3 and 10 extracted from Monte Carlo simulations (see later). The correlation between QUBE and QM energetics is similar to that previously reported, and significantly QUBE does not predict any physically unreasonable structures (either bound to the protein or in water) whilst retaining the fixed MM functional form that provides us with a practical method for deployment in free energy predictions. Additional analysis of torsion scans and correlations between QM and QUBE energetics for the remaining ligands may be found in the Supporting Information. Having parametrized the 18 inhibitors, we turn now to the computation of their relative binding free energies to p38 kinase. Free energy calculations were performed using the MCPRO software. The ligand binding site is expected to be hydrated, and so the JAWS water placement algorithm 32 was used to optimize the initial solvent distribution. As reported previously, the majority of the ligands 1-17 are expected to bind in pose 1 (Figure ).