|Home | About | Journals | Submit | Contact Us | Français|
Variations in hydrogen-bond strengths are investigated for complexes of nine para-substituted phenols (XPhOH) with a water molecule and chloride ion. Results from ab initio HF/6-311+G(d, p) and MP2/6-311+G(d, p)//HF/6-311+G(d, p) calculations are compared with those from the OPLS-AA and OPLS/CM1A force fields. In the OPLS-AA model, the partial charges on the hydroxyl group of phenol are not affected by the choice of para substituent, while the use of CM1A charges in the OPLS/CM1A approach does provide charge redistribution. The ab initio calculations reveal a 2.0-kcal/mol range in hydrogen-bond strengths for the XPhOHOH2 complexes in the order X = NO2 > CN > CF3 > Cl > F > H >OH >CH3 > NH2. The pattern is not well-reproduced with OPLS-AA, which also compresses the variation to 0.7 kcal/mol. However, the OPLS/CM1A results are in good accord with the ab initio findings for both the ordering and range, 2.3 kcal/mol. The hydrogen bonding is, of course, weaker with XPhOH as acceptor, the order for X is largely inverted, and the range is reduced to ca. 1.0 kcal/mol. The substituent effects are found to be much greater for the chloride ion complexes with a range of 11 kcal/mol. For quantitative treatment of such strong ion-molecule interactions the need for fully polarizable force fields is demonstrated.
In our development of non-nucleoside inhibitors of HIV-1 reverse transcriptase (NNRTIs), high sensitivity to substitution at the 4-position in the phenyl ring has been found for the thiazole series 1 and the pyrimidines 2, as summarized in Table 1.1,2 Specifically, for the thiazoles, there is a 50-fold enhancement in activity as the substituent X is made more electronegative in going from X = H to CN, while a 1500-fold enhancement is obtained in the pyrimidine series. The structures of the complexes of such NNRTIs with HIV-RT have been well established through X-ray crystallography and computation.3-5 A key feature is a short hydrogen bond between the amino group of the NNRTI and the carbonyl oxygen of Lys101 of HIV-RT (Figure 1). The β-nitrogen in the heterocycle is also in a longer hydrogen bond with the backbone NH of Lys101.
Questions that then arise are (a) how sensitive are such hydrogen-bond strengths to substitution in the phenyl rings, and (b) are such effects adequately reflected in the force-field calculations that are often used to examine the energetics of protein-ligand binding?6-8 For example, successful guidance of lead-optimization by performing free-energy perturbation calculations to predict the effects of changes in substituents on rings and of choices of heterocycles on binding affinities is expected to require proper representation of such effects.1,9,10 Electronic polarization is a central issue here since change in X alters the charge distribution including for the key hydrogen-bond donating hydrogen of the ligand.7,11 Related effects on acidities of substituted benzoic acids and phenols led to the development of the Hammett equation and the σ and σ− substituent constants.12
As summarized below, these ideas were pursued by performing ab initio calculations on prototypical hydrogen-bonded systems and comparing the results to those obtained from the non-polarizable OPLS-AA force field13 and its OPLS/CM1A variant,8 which incorporates polarized partial atomic charges that are obtained from the quantum mechanical CM1A procedure.14 The CM1A method, which is based on AM1 wavefunctions, was derived to reproduce dipole moments for organic molecules in the gas phase.14 The CM1A charges when enhanced by 14% (1.14*CM1A) were also found to perform well for computing free energies of hydration of 25 diverse organic molecules15 in explicit TIP4P water16 with all other force-field parameters taken from OPLS-AA. The systems chosen for initial study of substituent effects on hydrogen bonding are the phenol-water and phenol-chloride ion complexes, 3 – 5.
The principal interest here is comparison of ab initio and force field predictions for the effects of the substituents X on the hydrogen-bond strengths. Ab initio and density functional theory calculations were carried out with Gaussian 03 and all geometrical degrees of freedom were optimized for the complexes and separated components.17 In Table 2, results for the PhOHOH2 (3) complex and water dimer are compared at the HF/6-31G(d), B3LYP/6-31G(d), HF/6-311+G(d, p), and MP2/6-311+G(d, p)// HF/6-311+G(d, p) levels. The interaction energies at the latter two levels bracket what is accepted as the true value for the water dimer, −5.4 ± 0.7 kcal/mol from experiment18 and −5.1 ± 0.2 kcal/mol from theory.19 The ΔE results for the phenol-water complex are also similar to those from the highest-level calculations in a prior study, i.e. −7 to −8 kcal/mol at the MP2/aug-cc-pVDZ level.20 Thus, the substituent effects were explored with the HF/6-311+G(d, p) and MP2/6-311+G(d, p) calculations. Counterpoise corrections have not been made since they are expected to show little variation with the choice of substituent X.
The corresponding force field calculations were carried out first with the substituted benzenes described with the OPLS-AA force field21 and with the water molecule represented by the TIP4P model.16 Complete energy minimizations were carried out with the BOSS program22 except that the internal geometry of the water molecule is fixed in the TIP4P model, r(OH) = 0.9572 Å and HOH = 104.52°. Notably, in the reported OPLS-AA model for para-substituted benzenes, the net charge on the substituent plus attached benzene carbon atom is zero.21 This permits transferability that simplifies the modeling of arbitrary substituted benzenes, but it ignores associated polarization effects. Thus, the partial atomic charges on the COH group of all para-substituted phenols are the same (Figure 2). Though testing for numerous mono- and di-substituted benzenes has revealed modest average errors for computed free energies of hydration (0.5 kcal/mol) and pure liquid heats of vaporization (1.0 kcal/mol) and densities (0.02 g/mL),21,23 differential polarization effects are expected to be more apparent upon examination of specific hydrogen-bond strengths as in protein-ligand binding.
For the OPLS/CM1A approach,8 OPLS-AA parameters are used except for the partial atomic charges, which are obtained from the CM1A method.14 A sequence of geometry optimizations and CM1A calculations is performed with a BOSS script until the charges are converged. For neutral molecules, it is noted again that the CM1A charges are scaled by a factor of 1.14 for use in the OPLS/CM1A force field.15 For ions, the CM1A charges are not scaled to avoid non-physical net charges.24 The program also symmetrizes the charges for equivalent atoms, e.g., the charges are averaged for equivalent methyl hydrogens or the ortho carbons and hydrogens in Figure 2. Without the symmetrization, artifacts arise for general molecular modeling such as introduction of spurious minima in conformational searching. Optimization of the complexes is then performed with the converged charges and with the internally rigid TIP4P water molecule. Further polarization of the charge distribution for the substituted benzenes upon complex formation with water is not carried out. For the much stronger phenolchloride ion interactions, the importance of a full treatment of polarization effects is considered below.
The computed interaction energies ΔE for the complexes 3 with the phenol as the hydrogen-bond donor are summarized in Table 3 and the optimized OO distances are in Table 4. The trend in the ab initio results is as expected with electron-withdrawing substituents acidifying the hydroxyl group, increasing the hydrogen-bond strengths (Table 3), and decreasing the hydrogen-bond lengths (Table 4). However, there are fine points. For example, the π-donating character of the amino substituent outweighs its σ-withdrawing character to yield the weakest hydrogen bond. For fluorine and chlorine, the opposite pattern seems to be operative as the hydrogen bonds are stronger than for phenol (X = H) in those cases. As shown in Figure 3, the ab initio hydrogen-bond strengths roughly follow the trend of Hammett σ constants, though this is not fully expected in view of the differences in the processes, i.e., substituent effects on phenol-water hydrogen-bond strengths and on acidities of substituted benzoic acids in aqueous solution.
The substituent effects on the hydrogen-bond strengths are substantial with a 2 kcal/mol range from both the HF and MP2 calculations for the complexes 3. In view of the constancy of the OPLS-AA partial charges for the phenolic hydroxyl group, it is not surprising that the range for ΔE is compressed to 0.7 kcal/mol and the ordering of the values is poor. There is also negligible variation in the OO distances in Table 4 with OPLS-AA, while the HF results show a reduction of the hydrogen-bond lengths by 0.06 Å in going from X = NH2 to NO2 for the complexes 3. In contrast, use of the 1.14*CM1A charges nicely corrects the problems with the interaction energies and yields absolute values roughly midway between the HF and MP2 results (Table 3 and Figure 3). The level of accord was not anticipated, but it suggests that the polarization of the charge distributions by the CM1A method is remarkably accurate within the context of the simple point charge model for the force fields (single atom-centered partial charges). The hydrogen-bond lengths with OPLS/CM1A also decrease with increasing strength, though the range is less than from the HF optimizations (Table 4). The absolute hydrogen-bond lengths are 0.15 – 0.20 Å shorter from the force fields than from the ab initio calculations, which is normal for fixed-charge force fields that are intended for condensed-phase simulations.7,8,13,16
Turning to the complexes 4 with the phenol as the hydrogen-bond acceptor, the trends for hydrogen-bond strengths and lengths from the ab initio calculations are now opposite with electron-withdrawing substituents weakening the basicity of the phenolic oxygen. Thus the MP2 results for 3 range from −7.8 to −9.7 kcal/mol in going from X = NH2 to NO2, while the corresponding values for 4 are −5.4 to −4.9 kcal/mol (Table 5), i.e., opposite in trend, much weaker, and in a narrower range. Qualitatively similar results are obtained from the HF calculations, though the range for the complexes 4 is 1.0 kcal/mol. As before, the OPLS-AA results are too invariant, while the OPLS/CM1A model is successful in showing the weakening of the hydrogen bond for 4 with increasing electron-withdrawing character for the substituent X.
A key point from Tables Tables33 and and55 is the increasing gap between the substituted phenol's ability to act as a hydrogen bond donor and acceptor with increasing electron-withdrawing character for X. E.g., for p-cyanophenol as donor and acceptor the difference in hydrogen-bond strengths is 4.6 kcal/mol from the MP2 results and 4.3 kcal/mol with OPLS/CM1A, while the difference is only 1.7 kcal/mol from the OPLS-AA calculations. It is clear that (a) such modulation of hydrogen-bonding ability is important for proper description of intermolecular interactions, and (b) its accurate description requires methodology that allows polarization of the charge distributions. It is also apparent from Figure 2 and Tables Tables33 and and55 that the hydrogen-bond strengths are very sensitive to small changes in the atomic charges. The variations for the hydroxyl oxygen and hydrogen are only ca. 0.01 e with OPLS/CM1A; the variation for the ipso carbon is actually much greater, 0.1 e. For phenol itself, if the OPLS-AA charges for the hydroxyl oxygen and hydrogen are changed by 0.01 e, the strength of the hydrogen-bond for the complex 3 changes by ca. ±0.3 kal/mol in the expected manner. This sensitivity is well known and has always been a challenge in the development of force fields.7,8
Naturally, the substituent effects are much enhanced for the complexes with chloride ion, 5 (Table 6). It is noted that the OPLS chloride ion parameters (q = −1.0 e, σ = 4.02 Å, ε = 0.71 kcal/mol) that were used here are from a recent, comprehensive study of the hydration of halide and alkali ions.24 For the complexes with the substituted phenols, the ranges for the interaction energies are −17 to −28 (HF), −23 to −35 (MP2), −14 to −19 (OPLS-AA), and −15 to −27 kcal/mol (OPLS/CM1A). Some experimental data from high-pressure mass spectrometry are also available for complexes 5, as listed in Table 6.25 For anion-molecule complexes like these, conversion of the computed electronic energy change ΔE to ΔH at 298 K involves a correction of ca. +0.9 kcal/mol.26 It is apparent that the MP2 results are in close accord with the experimental data, while the HF and force-field results significantly underestimate the hydrogen-bond strengths. However, the OPLS/CM1A approach again does much better than OPLS-AA for the magnitude, pattern, and range for the interaction energies. In this case, polarization of the substituted phenol by the chloride ion can be expected to be significant and the fixed-charge OPLS-AA and OPLS/CM1A models are both inadequate. This is the case in spite of the fact that the optimal interaction energy for Cl− with a TIP4P water molecule of −13.0 kcal/mol is in good agreement with the best available estimates.24 The larger phenols are more polarizable than a water molecule.
For proper treatment of such strong ion-molecule interactions, it is accepted that a fully polarizable force field is required.7,11 Thus, we have been exploring the addition of inducible dipoles to the OPLS models. A simple approach has been taken by which an inducible dipole can be added to non-hydrogen atoms. Furthermore, the electric field that determines the inducible dipoles is only computed from the permanent charges (eq 1) and the total polarization energy is given by the usual formula, eq 2. The key approximation
is that the induced dipoles do not contribute to the electric field, which simplifies the computations since iterative solution for the dipoles is not required. Addition of the induced dipoles to OPLS-AA and OPLS/CM1A yields OPLS-AAP and OPLS/CM1AP. The implementations are residue-based in that the electric field at an atom is determined by the charges on all other atoms not in the same residue. For the complex 5, the substituted phenol and chloride ion are treated as separate residues, so the induced dipoles for the phenol are only determined by the field from the chloride ion. The same polarization model has been used by others27,28 and it performed well in a previous study of ours for reproducing solvent effects for the gauche/anti equilibrium for 1,2-dichloroethane in multiple solvents and the free energy of solvation of water in cyclohexane.29
Modest parameter optimization has been carried out for the polarizabilities α to reproduce gas-phase complexation energies (MP2/6-311G(d, p)) for ca. 30 ion-molecule complexes focusing on cation-π interactions. This led to setting αi = 1.0 Å3 for carbon and αi = 1.5 Å3 for heteroatoms. With these choices, the phenol-chloride ion complexes were optimized yielding the results in Table 7. Inclusion of the induced dipoles is found to enhance the interaction energies by 5-6 kcal/mol. The hydrogen-bond lengths are also shortened by ca. 0.1 Å to yield the values that are listed in Table 7. The agreement between the ab initio and OPLS/CM1AP results is certainly respectable, while the OPLS-AAP results still suffer from the underlying problems with the invariant partial charges in the OPLS-AA model. For the phenol-water complexes, addition of the inducible dipoles to the force fields strengthens the hydrogen bonds uniformly by 0.4 - 0.5 kcal/mol and shortens them by ca. 0.02 Å.
Substituent effects on the interaction energies for complexes of phenols with water and chloride ion have been investigated. For the complexes with water, the OPLS/CM1A model was found to yield good reproduction of ab initio results, while the OPLS-AA force field with invariant partial charges for the hydroxyl group compresses the substituent effects. For the complexes with chloride ion, the interaction energies and substituent effects are much magnified, and the need for explicit treatment of the intermolecular polarization energy is apparent. Addition of inducible dipoles on non-hydrogen atoms was found to enhance the interaction energies by 5-6 kcal/mol. The resultant OPLS/CM1AP model performed well and warrants further investigation. It is emphasized that the ability to predict accurately substituent effects on intermolecular interactions is central to key applications of molecular modeling, for example, in the design of drugs, materials, and catalysts.
Gratitude is expressed to the National Science Foundation (CHE-0446920) and National Institutes of Health (GM32136) for support.
Supporting Information Available: An Excel file is provided with the ab initio energies for the substituted phenols and complexes 3 - 5. These materials are available free of charge via the Internet at (http://pubs.acs.org/instruct/jacsat.pdf).