Theoretical Prediction of Hydrogen-Bond Basicity pKBHX Using Quantum Chemical Topology Descriptors

Hydrogen bonding plays an important role in the interaction of biological molecules and their local environment. Hydrogen-bond strengths have been described in terms of basicities by several different scales. The pKBHX scale has been developed with the interests of medicinal chemists in mind. The scale uses equilibrium constants of acid···base complexes to describe basicity and is therefore linked to Gibbs free energy. Site specific data for polyfunctional bases are also available. The pKBHX scale applies to all hydrogen-bond donors (HBDs) where the HBD functional group is either OH, NH, or NH+. It has been found that pKBHX can be described in terms of a descriptor defined by quantum chemical topology, ΔE(H), which is the change in atomic energy of the hydrogen atom upon complexation. Essentially the computed energy of the HBD hydrogen atom correlates with a set of 41 HBAs for five common HBDs, water (r2 = 0.96), methanol (r2 = 0.95), 4-fluorophenol (r2 = 0.91), serine (r2 = 0.93), and methylamine (r2 = 0.97). The connection between experiment and computation was strengthened with the finding that there is no relationship between ΔE(H) and pKBHX when hydrogen fluoride was used as the HBD. Using the methanol model, pKBHX predictions were made for an external set of bases yielding r2 = 0.90. Furthermore, the basicities of polyfunctional bases correlate with ΔE(H), giving r2 = 0.93. This model is promising for the future of computation in fragment-based drug design. Not only has a model been established that links computation to experiment, but the model may also be extrapolated to predict external experimental pKBHX values.


INTRODUCTION
Recently, fragment-based drug design has emerged as a popular method for screening small chemical structures against a known drug target. 1−4 Modifications to the lead compound serve to improve the target selectivity and administration/absorption, distribution, metabolism, excretion/elimination, toxicology (ADMET) profile. The resulting modifications aim to increase the chance of discovering a novel drug candidate. The probability of discovering a novel drug candidate is increased if the modifications to the lead compound are informed.
One concept broadly applied to help medicinal chemists make informed modifications to lead compounds is that of bioisosterism. 5−9 A bioisostere is a substituent or a fragment of an active compound with similar physical or chemical properties. 5 The concept of bioisosterism can be applied during the lead optimization process. Bioisosteric modifications can be made to lead compounds with the aim being to improve the ADMET profile of the lead. Only the properties important to the biological activity in each specific lead optimization need to be enhanced. However, the properties that are enhanced to improve the ADMET profile of one lead might unfortunately deliver contrasting results in the case of a different lead. Therefore, successful modifications of a lead are unique and specific to the lead itself. The unique combination of modifications that are specific to each successful lead optimization is often unknown. As the specific modifications are unknown for each lead, a number of broad properties have been outlined 7 that can be modified practically in the lead optimization process. These properties relate to both the drug molecule and its environment, and serve as a starting point when considering bioisosteric substitutions. The properties are size, shape, electronic distribution, lipid solubility, water solubility, pK a , chemical reactivity, and hydrogen-bonding capacity.
Hydrogen bonding is an important noncovalent interaction between a biological molecule and its local environment. 10−13 Hydrogen bonding may influence binding affinity 14−16 and solvation, 17,18 membrane permeability, 19 and adsorption onto surfaces. 20 Hydrogen bonding influences both the distribution of a drug to its target and the binding of a drug to its target. After the drug is administered and absorbed into the body, it must be delivered to the target. As the drug is absorbed into the body, it interacts with the biological fluids that transport the drug around the body. Any exposed polar residues of the drug are capable of forming hydrogen bonds with the biological fluids as the drug becomes dissolved. The drug must then shed the solvent to bind the target. The hydrogen-bonding network involved in the displacement of the solvent with the target influences the thermodynamics of a drug interaction. The Gibbs energy of a drug interaction can be improved by focusing on hydrogen bonding. There is a favorable enthalpy change as the drug forms hydrogen bonds with the target. However, there is also an unfavorable enthalpy change as hydrogen bonds are broken as the drug sheds the layer of solvent upon binding. There is a favorable entropy change as the layer of solvent is shed upon binding. During binding, there is, however, an entropic penalty due to the loss of translational and rotational freedom. Favorable enthalpy and entropy changes must be maximized while minimizing unfavorable changes to improve the binding profile of a lead compound. One way in which this may be achieved is to modify the lead in such a way that the strongest hydrogen bonds are formed between the drug and its respective target. Knowledge of relative hydrogen-bond strengths is therefore important to medicinal chemists in lead optimization.

BACKGROUND
2.1. Quantum Chemical Topology (QCT). Quantum chemical topology is a method that extracts chemical information from wave functions by using 21,22 the language of dynamical systems. QCT generalizes the original quantum theory of atoms in molecules (QTAIM) 23,24 and builds on the idea of partitioning quantum mechanical property densities by means of a (gradient) vector field. This approach defines finite volumes in 3D space that express the partitioning of the quantum system. In this Article, we only use topological properties retrieved from the electron density (for details and definitions, see section 3). Bond critical points (BCPs), also denoted (3,−1), are featured in this study, together with the quantum mechanical functions evaluated at them. These local properties are complemented by global properties, which are obtained by 3D integration over the volume of a topological atom. These concepts have been reviewed many times before.
2.2. Hydrogen-Bond Basicity Scales. Despite the importance of hydrogen bonding in biological systems, there is a general lack of understanding of the relative strengths of hydrogen-bond basicity. 25 The first scales of hydrogen-bond basicity were set up by Taft and co-workers, 26,27 who investigated linear free energy relationships (LFER) between the reference acid 4-fluorophenol and a series of OH hydrogenbond donors (HBDs). The HBDs were hydrogen bonded to a series of bases in 1:1 complexes. The equilibrium constants of these complexes were used to define the relative basicities of the hydrogen-bond acceptors (HBAs). The corresponding scale was called pK HB . However, the pK HB scale contains only 117 basicity values. Second, due to the nature of the LFERs, pK HB values are limited to (−OH) HBDs. When the LFERs were repeated using 5-fluoroindole instead of 4-fluorophenol as the reference acid, family-dependent behavior was observed. 28 Therefore, the pK HB scale covers only a limited proportion of all possible hydrogen bonds that may be formed, and has been made little use of by medicinal chemists.
The log K B H scale, established by Abraham and co-workers, 29 was constructed using 34 reference acids complexed to a series of bases in a 1:1 ratio. The log K B H scale was intended to establish an overall generality of hydrogen-bond basicity. The log K B H values were obtained from 34 LFER equations in the form of eq 1: where the index i refers to the base, A represents the acid, L A and D A are the gradients, and log K i is the intercept, respectively, of the 34 straight line equations.
In this equation, SP relates to a set of solute water−solvent partition coefficients in a given system, R 2 is the excess molar refraction, π 2 H is the dipolarity or polarizability, and α 2 H and β 2 H are the hydrogen-bond acidity and hydrogen-bond basicity. The variable V x is McGowan's characteristic molecular volume (cm 3 mol −1 /100). 31 The parameters of eq 2 can be determined experimentally with the exception of the acidity or basicity terms. The solutes were restricted either to give an acidity term of 0 or to monoacidic solutes for which the acidity term is known. 18 Hence, the only unknown term is the overall basicity. The ∑β 2 H values were estimated applying β 2 H values to eq 2 and back-calculating to return successive ∑β 2 H values in an iterative process until a self-consistent set of equations and therefore ∑β 2 H were obtained. It was therefore claimed that the ∑β 2 H scale of overall hydrogen-bond basicity was more useful than the β 2 H scale when considering biological properties. 30 Despite being the most quoted scale of hydrogen-bond basicity, 25,32 and despite the authors' claim that it is more useful than 1:1 scales in the interpretation of biochemical and physiochemical properties, 30 there are problems in the interpretation of ∑β 2 H values. The ∑β 2 H scale assigns only one ∑β 2 H value to the whole molecule, and is based on partition coefficients rather than binding. Although medicinal chemists may make use of ∑β 2 H values to predict permeability, the ∑β 2 H scale offers no thermodynamic binding information.
It is clear that there is a need for a measure of hydrogen-bond basicity that targets the needs of medicinal chemists. The requirement for this scale is that it should combine the simplicity of the pK HB scale and its thermodynamic relevance with the analysis of polyfunctional bases as in the ∑β 2 H scale. The pK BHX database 25 seems to meet these requirements. The pK BHX scale is essentially an extension of the pK BH scale in that it is based on equilibrium constants of 1:1 complexes in the solvent CCl 4 where 4-fluorophenol is the reference acid. The pK BHX values are therefore related to Gibbs energy and have thermodynamic relevance. However, the pK BHX scale also accounts for polyfunctional bases.
An important difference between the pK BHX scale, on one hand, and the pK BH and ∑β 2 H scales on the other hand, is the experimental method used to calculate the equilibrium constants. The difference in experimental methods allows polyfunctional bases to be included in the analysis. Fourier transform infrared (FTIR) spectrometry is used in the determination of the pK BHX scale rather than 19 F NMR, UV, and dispersive IR techniques, 25 which are mostly used for pK BH and ∑β 2 H . The advantage of using FTIR is that multiple significant hydrogen-bond acceptor sites of polyfunctional bases can be analyzed. 25 During the formation of a hydrogen bond, the stretching vibration of the XH bond is shifted to a lower wavenumber, giving rise to peaks in a FTIR spectrum. Each peak represents a hydrogen-bond site. The change in the vibrational frequency of the XH bond can be translated into a pK BHX value through family-dependent pK BHX −IR shift relationships. 33−36 The family-dependent pK BHX −IR relationships correlate the change in IR vibrational frequency of the XH bond to pK BHX values for a family of HBAs and a reference HBD in 1:1 complexes. The FTIR spectroscopy method allows equilibrium concentrations of these 1:1 complexes to be measured in a certain frequency range, therefore allowing a pK BHX value to be obtained. This allows a pK BHX value to be assigned to each HBA site separately. It is clear that the pK BHX scale is able to analyze polyfunctional bases as well as the ∑β 2 H scale. An important difference, however, is that pK BHX values are thermodynamically significant because each basic site is considered separately. On the other hand, the ∑β 2 H scale does not consider each hydrogen-bond acceptor site individually. This feature is illustrated by the fact that the pK BHX database contains 1338 pK BHX values taken from 1164 bases. This large data set is advantageous to medicinal chemists as most drugs contain several hydrogen-bond acceptor sites consisting of mainly O, N, S, F, Cl, and sp 2 hydridized carbons. It is therefore a possibility that the pK BHX scale could be used as a descriptor characterizing hydrogen bonding from which bioisosteric replacements could be found.
The ∑β 2 H scale considers the entire molecule to be thought of as the hydrogen-bond acceptor, whereas the pK BHX scale considers each individual hydrogen-bond acceptor site independently. This has important implications for fragmentbased drug design. The ∑β 2 H scale depends on the entire molecule being intact to return a ∑β 2 H value. The pK BHX scale, however, allows for each hydrogen-bond acceptor site to be assigned a pK BHX value. Therefore, a fragmented drug-like molecule may be assigned pK BHX values but not ∑β 2 H values.

Computational Predictions of Hydrogen-Bond
Basicity Values. There have been a number of computational models used to predict hydrogen-bond basicity. The minimum value of the electrostatic potential at the hydrogen-bond acceptor atom 37,38 has been used to obtain a family-dependent model for the hydrogen-bond basicity parameter β used in linear solvation energy relationships. 39 A more general model displaying a reasonable family-independent correlation, r 2 = 0.81, between log K HB and the minimum value of the electrostatic potential has also been established. 40 Basicity is typically more difficult than acidity to model due to the variety of different HBA sites. By considering the nitrogen HBAs separately, the correlation can be improved when correlating pK HB to the minimum of the electrostatic potential, 41 and refined further by also including the interaction energy, giving r 2 = 0.99. Platts 42 developed a model to predict hydrogen-bond basicity based on DFT calculations on a series of bases taken from Abraham's list of ∑β 2 H values. 43 Theoretical properties were calculated for both the isolated base and the base complexed with HF. Once again, strong family-dependent relationships were observed. The atomic multipole moments of nitrogens occurring in the isolated bases gave the best correlations, r 2 = 0.786. The combination of the minimum electrostatic potential at the nucleus of the HBA and an atomic charge (according to QCT) of the isolated base gave a better correlation, yielding r 2 = 0.887 for oxygen HBAs. However, with the use of a HF probe, improvements were made with the oxygen model yielding r 2 = 0.939. The binding energy calculated from the supramolecular complex correlated reasonably well, r 2 = 0.901, to give a good family-independent model. It could be argued that the ∑β 2 H scale might not be the best measure of basicity for the medicinal chemist to predict. This is due to the lack of HBA site-specific information and thermodynamic information offered by the scale. A recent study repeated the supramolecular approach of Platts, this time using phenol as the HBD, and found 43 good general correlations (r 2 = 0.94) between log K β values and the hydrogen-bond binding energy for nitrogen and oxygen HBAs.
Platts later refined his initial model correlating free energy of formation with pK HB values. 44 Fits of r 2 = 0.969 were observed for a family-independent model using HF as the HBD. The final refinement to this model was made by combining the minimum of the electrostatic potential at the HBA and properties of the hydrogen-bond critical point (BCP), 45 leading to r 2 = 0.97. Devereux et al. suggested 46 this model to be too expensive computationally for use in QCT-aided drug design. 47 These authors constructed a model 46 by combining the minimum, median, and mean of the electrostatic potential at the HBA with the minimum of the local ionization energy on the HBA surface and correlating with β 2 H values. They observed family-dependent models with non-nitrogen bases producing fits of r 2 = 0.931, while nitrogen bases gave r 2 = 0.871. However, the results were found to be less accurate than Platts' model, 42 which correlated properties of hydrogen-bonded complexes with ∑β 2 H values. It would therefore seem that the most accurate models used to predict hydrogen-bond basicity must be constructed from properties of the hydrogen-bonded complex. Also, by using β 2 H values, the model is limited to the 1:1 hydrogen-bonded complexes only. However, in drug··· target interactions, 1(HBA):n(HBD) complexes are likely to occur, and an improvement to the model could be made by correlating to a measure of hydrogen-bond basicity that accounts for all major hydrogen-bond acceptor sites of a particular molecule.
The previous models are limited by 1:1 ratios. However, drug design warrants the inclusion of polyfunctional bases and therefore needs a model that goes beyond the 1:1 ratio. Such a refined model must keep the accuracy of Platts' model by considering properties of the hydrogen-bonded complex, while taking into account polyfunctional hydrogen-bond acceptors. A good candidate to fulfill the above requirement is the pK BHX scale. The advantage of using pK BHX values over ∑β 2 H values is that a pK BHX value is listed for each hydrogen-bond acceptor site. This important feature allows for the correlation between a

METHOD AND COMPUTATIONAL DETAILS
A training set of 41 bases was selected from the pK BHX database. 25 The bases were chosen to include a wide range of pK BHX values for nitrogen, oxygen, and sulfur acceptor atoms. The strength of the nitrogen bases ranges from weaker nitrile and aniline acceptors, to stronger pyridine and amines. The oxygen bases include weaker phenol and alcohols, and stronger carbonyl and amide acceptors. A small set of sulfide and thiol sulfur acceptors also formed part of the training set. The distribution of the pK BHX values within the training set was chosen to roughly match that of the entire pK BHX database. 25 Hydrogen-bond complexes (1:1) were formed between the 41 training bases and 6 possible HBDs. The latter consist of 4 OH donors (water, methanol, 4-fluorophenol, and serine), complemented by methylamine and hydrogen fluoride. The pK BHX scale was set up using 4-fluorophenol as the reference HBD, and is thought to be applicable to all OH, NH + , and NH donors. 25 Therefore, the chosen HBDs in this study, with the exception of hydrogen fluoride, should give complexes that the pK BHX scale can be applied to.
DFT calculations were performed using Gaussian 03. 48 A set of hydrogen-bonded complexes were set up, based on those investigated by Platts. 42 They were generated by placing a hydrogen atom of the HBD approximately 2.0 Å away from the primary HBA site. The essence of the hydrogen-bond system is then conveniently denoted by D−H···A, where H···A is the hydrogen bond and D−H is the covalent HBD bond. The hydrogen-bonded complexes were geometry optimized at the B3LYP/6-311++G(2d,p) level of theory. Properties of the complexes were obtained and correlated with their respective pK BHX values. The properties include the hydrogen-bond length r(A···H), the change in the HBD bond length Δr(H−D), and the hydrogen-bond energy ΔE, the details of which are explained below. Note that the symbol Δ denotes the difference of the complex and the constituent monomers forming it. Atomic properties are obtained by a volume integral over an atomic basin Ω. The kinetic energy densities G(r) and K(r), and the atomic energy E(Ω), are defined in eqs 3, 4, and 5: where G(r) is one type of kinetic energy density, while K(r) is another type. The symbol ∫ dτ′ refers to the integration over the spatial and spin coordinates of all electrons except one, dr denotes a three-dimensional integration over the basin of the topological atom Ω, while R represents the virial ratio correction factor. 23 It is an important property of a topological atom that its kinetic energy is well-defined (which is not true for an arbitrary subspace). Indeed, G(Ω) = K(Ω).
The atomic charges are defined in eqs 6 and 7: where Z is the charge of the nucleus inside the topological atom. All atomic integrations used to analyze atomic properties were calculated with the AIMAll suite. 49 The atomic properties used here are ΔE(H), or the change (going from the monomers to the complex) in the HBD energy, and Δq(H), or the change in charge of the hydrogen atom in the HBD. Bond properties were evaluated and analyzed only at BCPs associated with H···A hydrogen bonds and BCPs of the covalent HBD bonds (D−H), both in the hydrogen-bond system D−H···A. Properties analyzed at the hydrogen BCP are the electron density ρ(A···H), the Laplacian of the electron density ∇ 2 ρ(A···H), the kinetic energy density at the BCP G(A···H), and the kinetic energy density G(r) at the BCP, divided by the electron density G/ρ(A···H). Properties of the HBD critical point are the change in the electron density Δρ(H−D), and the change in the Laplacian of the electron density Δ∇ 2 ρ(H−D). Properties of the BCPs were obtained from the MORPHY98 program. 50,51 This program also generated Figure 1, which shows an atomic basin of a HBD and all BCPs of a hydrogen-bonded complex.
Basis set superposition error (BSSE) calculations were carried out to accurately estimate the hydrogen-bond binding energy for the water series only. As the BSSE of hydrogen-bond energy turned out to be very small, BSSE calculations were not repeated for complexes without water. The method used to perform the BSSE calculations is that of Handy and coworkers. 52  calculates the BSSE of the hydrogen-bond energy using the counterpoise correction method. 53 In eq 9, the asterisk refers to the geometry of the HBA and HBD in the optimized complex, and CP stands for counterpoise and refers to the energy of the HBA and HBD in the presence of ghost HBD and HBA Gaussian functions,   54−56 To investigate this criticism in the context of the current work, a further set of calculations using the MPWB1K/6-311++G(2d,p) level was performed for the training set of bases where water was used as a HBD.
To account for polyfunctional bases, calculations were performed at the B3LYP/6-311++G(2d,p) level on a set of 11 bases with 2 hydrogen-bond acceptor sites giving 22 acceptor sites in total. Initially, 22 1:1 complexes were generated. The properties calculated from the 1:1 complexes have been compared to properties calculated from the same set of bases in 2:1 complexes. The HBD chosen for these calculations was methanol. Table 1 list the names of the 41 bases (i.e., HBAs) and their corresponding pK BHX values.

RESULTS AND DISCUSSION
Each of the HBAs listed in Table 1 has been used to form a hydrogen-bonded complex with each of the six HBDs, which are water, methanol, 4-fluorophenol, methylamine, serine, and HF. Examples of these complexes are given in Figure 2, except for HF, illustrating each possible element (S, O, and N) as an acceptor. Table 2 summarizes the correlations of calculated properties (energetic, geometric, atomic, and BCP) with the experimental pK BHX values. With the exception of the hydrogen fluoride set, the quantity ΔE(H) consistently outperforms all other properties, giving r 2 values of 0.96 for water, 0.95 for methanol, 0.91 for 4-fluorophenol, 0.93 for serine, and 0.97 for methylamine. Particular attention was given to the electronic energy of the hydrogen bond, ΔE. However, poor correlations across the board of OH and NH donors were observed. The BSSE corrected energy offered no improvement to the correlations. However, the plot of ΔE against pK BHX where hydrogen fluoride was used as the probe gave a correlation of r 2 = 0.60. Although 0.60 is a weak correlation, it is significantly better than any of the values obtained for OH and NH donors (which are generally about 0.30 with the exception of 4fluorophenol with r 2 = 0.52). Lamarche and Platts 44 obtained a value of r 2 = 0.91 for the correlation between computed ΔE and Taft's pK HB values using hydrogen fluoride as the HBD. Taking into account the stronger correlations of ΔE with pK BHX and the observations of Lamarche and Platts when hydrogen fluoride is used as the HBD, it could be said that any relationship between ΔE and basicity is highly questionable due to the use of hydrogen fluoride as a probe. Indeed, the computed ΔE of a hydrogen fluoride complex correlates reasonably well to basicity values, but a link cannot be made between theory and experiment as the basicity scales are not applicable to hydrogen fluoride complexes.
The MPWB1K density functional has been shown to perform better than B3LYP when computing the interaction energies of hydrogen-bonded complexes. 54−56 Table 2 reveals that hydrogen-bond binding (interaction) energy correlates with pK BHX for MPWB1K calculations but not B3LYP calculations. However, as the binding energies of the complexes used here do not correlate with pK BHX values as strongly as ΔE(H), it is adequate to use B3LYP. Furthermore, Table 2 shows how ΔE(H) calculated from a MPWB1K wave function performs equally as well as ΔE(H) calculated from a B3LYP wave function when correlating ΔE(H) to pK BHX .   Figure 3 shows the plots of ΔE(H) against pK BHX for water, methanol, 4-fluorophenol, serine, and methylamine. However, pK BHX values did not correlate with ΔE(H) when HF was used as the HBD, as made from the very low correlation coefficient r 2 = 0.04. It is claimed 25 that the pK BHX scale is applicable to hydrogen-bonded complexes provided the HBD's functional group is either OH, NH, or NH + . Therefore, hydrogen fluoride was chosen as a control probe because hydrogen-bonded complexes, where hydrogen fluoride is the HBD, have no link to the experimental pK BHX scale. The relationship between ΔE(H) and pK BHX values gains credibility as it holds only for HBDs that the pK BHX scale may be applied to. Therefore, a strong link between theory and experiment has been established through the partitioned atomic energy of the HBD. Essentially, the energy assigned to a theoretically derived atomic basin is related to a molecular free energy quantity that is derived from experiment. It is also astonishing that ΔE(H) is able to predict the pK BHX of an external data set.  An external data set of 41 compounds has been chosen to demonstrate that the QCT property ΔE(H) can predict pK BHX values of HBAs chosen from the pK BHX database. The HBAs that form the external data set are listed in Table 3. Figure 4 plots the predicted versus experimental pK BHX values for the HBAs listed in Table 3. Predicted pK BHX values are calculated from eq 11: 180.94 (H) 3.7074 BHX (11) which corresponds to the least-squares line in panel B in Figure  3. This panel plots ΔE(H) against pK BHX , where methanol is the HBD, complexed to the HBAs of Table 2. Methanol was chosen as the probe to use in the model for three reasons: (i) its simplicity, (ii) its use in obtaining secondary LFER derived experimental values, and (iii) the strong correlation of r 2 = 0.95 between the ΔE(H) of these complexes and their pK BHX values.
Although the correlation between ΔE(H) and pK BHX for the methylamine set (r 2 = 0.97) is stronger than the methanol set, it is challenging to model specific hydrogen bonds with methylamine due to the presence of an additional hydrogen attached to the nitrogen as compared to OH donors that have a single hydrogen. This additional hydrogen may cause the molecule to rotate during optimization and form interactions other than the desired one, a characteristic often observed in this work due to no geometrical constraints placed upon the optimization. Therefore, it is easier to describe pK BHX in terms of ΔE(H) on water or methanol HBDs. An F-test was performed, and the variances between the water and methanol models were not found to be significantly different. Following, a paired t test assuming equal variances was performed. Again, the water and methanol models were not found to be significantly different. Therefore, there is no statistical advantage to choosing methanol as the probe over water.   Water would in fact be computationally less expensive than methanol. However, methanol has been chosen due to its usefulness for obtaining secondary pK BHX values from LFERs. The advantage that methanol can be used in such a way justifies the slightly lengthier calculations. Therefore, by using methanol as the probe, the computation is kept in close proximity to the experimental procedure. Figure 4 shows a relatively strong correlation, r 2 = 0.90, for the plot of predicted versus experimental pK BHX values. This demonstrates impressive extrapolation of the original model. The only constraint put on the compounds in Table 3 is that the HBA must be either an oxygen, a nitrogen, or a sulfur atom. Not only does the model based on ΔE(H) show impressive external predictability, but also good generality. This is because all HBAs can be grouped together in one data set. There seems to be no need to separate the HBAs into atomic groups (i.e., families) nor build individual models for each HBA type.
The results discussed so far are for 1:1 HBD:HBA complexes where the bases have only one major HBA site. The attractiveness of the pK BHX scale to the medicinal chemist lies in its inclusion of data for the more realistic cases of polyfunctional bases with two or more major HBA sites. Table 4 shows a list of 11 HBAs, each with two major HBA sites (so leading to 22 entries). Two approaches were used to test the validity of the ΔE(H) model for polyfunctional bases. The first approach simply involved generating a 1:1 complex for each HBA site. The second approach aimed to mimic experiment in which the bases are surrounded by an excess of acid allowing each HBA site to reach equilibrium with a HBD. The way in which the experiment was mimicked consisted of generating 2:1 HBD:HBA complexes where both HBA sites are occupied in the same calculation. Figure 5 shows the results of the two approaches. The first approach in which 1:1 complexes were generated marginally outperforms the second approach in which 2:1 complexes were computed. This may be surprising as 2:1 complexes are closer to the experimental approach. However, the fact that 1:1 complexes slightly outperform 2:1 complexes is beneficial to the computational chemist because 1:1 complexes are computationally less expensive than the 2:1 complexes.

CONCLUSIONS
It has been shown that the energy, E(H), assigned to a theoretically derived atomic basin of a molecule in the gas phase is related to a molecular free energy quantity that has been derived from experiment, pK BHX . Furthermore, the relationship breaks down when hydrogen fluoride is used as the hydrogen-bond donor. Therefore, the relationship between computation and experiment is only observed with computations that use hydrogen-bond donors that the experimental scale is applicable to.
It has also been shown that pK BHX values outside the initial data set are accurately predicted by the model set up with the initial data set. Therefore, ΔE(H) is able to predict the pK BHX values of an external data set.
The attraction of the pK BHX scale is that a pK BHX value is listed for each hydrogen-bond acceptor site on polyfunctional bases. It has also been found that ΔE(H) relates to pK BHX values of polyfunctional bases. This is an important observation when considering site-specific thermodynamic data and could help to improve computational drug design and development.

Notes
The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS
A.J.G. is grateful for a studentship funded by the BBSRC with sponsorship of Syngenta Ltd. ■ REFERENCES Figure 5. Correlation between the QCT ΔE(H) value and the pK BHX values for the 11 polyfunctional bases (each with two HBA sites) listed in Table 4