## Abstract

Combined quantum mechanics/molecular mechanics (QM/MM) methods are increasingly important for the study of chemical reactions and systems in condensed phases. Here, we have tested the accuracy of a density functional theory-based QM/MM implementation (B3LYP/6-311+G(d,p)/CHARMM27) on a set of biologically relevant interactions by comparison with full QM calculations. Intermolecular charge transfer due to hydrogen bond formation is studied to assess the severity of spurious polarization of QM atoms by MM point charges close to the QM/MM boundary. The changes in total electron density and natural bond orbital atomic charges due to hydrogen bond formation in selected complexes obtained at the QM/MM level are compared with full QM results. It is found that charge leakage from the QM atoms to MM atomic point charges close to the QM/MM boundary is not a serious problem, at least with limited basis sets. The results are encouraging in showing that important properties of key biomolecular interactions can be treated well at the QM/MM level employing good-quality levels of QM theory.

## 1. Introduction

Theoretical studies of chemical reactions in complex systems require a potential function that can describe electronic changes in the region of interest. To approach chemical accuracy (e.g. errors in energy differences of approx. 1 kcal mol^{−1}), correlated electronic structure methods are required, but the high cost of these methods limits the size of system studied (Claeyssens *et al*. 2006; Mulholland 2007). Combined quantum mechanical (QM) and molecular mechanical (MM) models have been widely used to overcome this limitation (Warshel & Levitt 1976; Aqvist & Warshel 1993; Gao & Thompson 1998; Monard & Merz 1999; Murphy *et al*. 2000). In these hybrid QM/MM approaches, the reactive part of the system is treated quantum mechanically and a classical description is applied to the surrounding environment (Warshel & Levitt 1976; Aqvist & Warshel 1993; Field 1997; Gao & Thompson 1998). QM/MM methods are well suited for studying enzyme-catalysed reactions that typically take place in a solvent and protein environment involving several thousand atoms (Warshel & Levitt 1976; Aqvist & Warshel 1993; Monard & Merz 1999; Murphy *et al*. 2000; Harvey 2004; Mulholland 2005).

The QM/MM Hamiltonian can be written as(1.1)where *H*_{QM} and *H*_{MM} are the QM and MM Hamiltonians that correspond to the atoms in the QM and MM regions, respectively. In equation (1.1), the QM/MM coupling term, *H*_{QM/MM}, describes the interaction between atoms in the QM and MM regions, and consists of electrostatic, van der Waals and bonded interactions (Warshel & Levitt 1976; Field *et al*. 1990; Aqvist & Warshel 1993; Thery *et al*. 1994; Humbel *et al*. 1996; Gao & Thompson 1998; Antes & Thiel 1999; Monard & Merz 1999; Murphy *et al*. 2000; Reuter *et al*. 2000). The bonded interactions are included only where the partitioning of the QM and MM regions intersects covalent bonds. Different techniques have been used to treat the broken bond within the QM part of the calculation, with the most common being the link atom (Singh & Kollman 1986; Field *et al*. 1990; Humbel *et al*. 1996; Antes & Thiel 1999; Swart 2003; Bathelt *et al*. 2005*b*) and frozen orbital (Thery *et al*. 1994; Murphy *et al*. 2000) approaches. Gao and co-workers have introduced a generalized hybrid orbital method, in which hybrid orbitals on boundary atoms are constructed. One of these orbitals is an active orbital included in the QM part of the QM/MM calculation, with the remaining three orbitals being used as frozen auxiliary orbitals (Gao *et al*. 1998; Pu *et al*. 2004, 2005). By comparing these techniques in QM/MM calculations, Reuter *et al*. (2000) concluded that the results obtained using link atoms with proper constraints are comparable to results from the local self-consistent field frozen orbital method (Assfeld & Rivail 1996). Here, we do not examine the effects of link atoms or schemes used to account for QM/MM partitioning along covalent bonds.

Several methods have been used to describe the electrostatic interaction between QM and MM regions (Gao & Thompson 1998), of which an ‘electrostatic embedding’ scheme is the most common. In this model, interactions with the point charges on the MM atoms are included in the one-electron Hamiltonian of the QM region (Bakowies & Thiel 1996; Antes & Thiel 1999)(1.2)where the Hamiltonian is expressed in atomic units, such that the usual 1/4π*ε*_{0} term equals 1. *Q*_{M} represents the magnitude of the point charges on the MM atoms and *R*_{M} and *r*_{N} represent the positions of the MM atoms and the QM electrons, respectively. This model directly allows for the electronic polarization of the QM region by the MM environment. One concern about this approach is that the use of point charges (as opposed to a more extensive multipolar expansion), parametrized to give accurate results in MM force-field calculations, may not yield an accurate enough description of the polarization effects within the QM region. Also, there has been debate in the literature (Singh & Kollman 1986; Field *et al*. 1990; Waszkowycz *et al*. 1991; Das *et al*. 2002; Laio *et al*. 2002) on possible unphysical electrostatic interactions between QM and MM atoms near to the partition region, and on charge leakage from the QM region to MM atoms. Clearly, using point charges on MM atoms far from the QM/MM boundary to construct the potential felt by the electrons in the QM region is reasonable. However, since Pauli repulsion corresponding to electrons in the MM atoms is not accounted for in the above scheme, it does not describe the short-range interactions between QM electrons and MM atoms correctly. One particular concern is that positive charges on MM atoms close to the QM/MM boundary may act as a trap for (QM) electrons. In the literature, this excess polarization nearer to the QM/MM boundary has been suggested as a possible source of error in QM/MM calculations (Das *et al*. 2002; Laio *et al*. 2002; Guallar *et al*. 2003; Biswas & Gogonea 2005; Dulak & Wesolowski 2006).

A few models have been introduced in QM/MM methods to damp the interaction between the QM electrons and the MM point charges close to the boundary, e.g. by using a delocalized charge distribution instead of a point charge (Eichinger *et al*. 1999; Das *et al*. 2002; Biswas & Gogonea 2005). In a related approach, Laio *et al*. (2002) introduced a modified electrostatic potential in the interaction Hamiltonian of their hybrid Car–Parrinello molecular dynamics model. In both techniques, the electron density of the MM atoms is represented as a delocalized distribution such that the potential due to the point charges representing the MM atoms is finite at the QM/MM boundary, and asymptotically converges to the potential of a point charge far from the boundary. Excluding the point charge of some of the MM atoms close to the QM/MM boundary from QM/MM interactions is also suggested by several research groups (Singh & Kollman 1986; Field *et al*. 1990; Waszkowycz *et al*. 1991). However, several published results reveal that the best results are obtained when all MM point charges interact with the QM region, except those situated on link atoms (Vasilyev 1994; Reuter *et al*. 2000). Modification of MM atomic charges on neighbouring (bonded) groups is also often recommended (Field *et al*. 1990; Mulholland & Richards 1997). However, no systematic investigation has been performed to find the extent of polarization effects at the QM/MM boundary due to the point charges of MM atoms, except one by Curutchet *et al*. (2003) employing Kitaura–Morokuma (KM) energy decomposition analysis (Kitaura & Morokuma 1976) at the HF/6-31G(d) level for a set of hydrogen-bonded systems.

Electron density maps have been widely used to analyse chemical properties (Geerlings *et al*. 2003; Leyssens *et al*. 2005; Merino *et al*. 2005). Unlike local quantities, such as bond length, total electron density maps can yield information about the global rearrangement of electrons due to complex formation or chemical reaction. Hence, electron density plots can be expected to be a sensitive tool for investigating charge polarization of the QM region and charge transfer between the QM and MM regions. In the present study, the change in total density due to hydrogen bond formation in selected systems calculated at the QM/MM level is compared with the full QM result. Changes in electron density upon dimer formation calculated at the QM level have been used by a number of groups to examine the polarization of one system by the other (Yamabe & Morokuma 1975; Krijn & Feil 1988; Hannachi *et al*. 1991; Lin *et al*. 1994; Bartha *et al*. 2003). In this work, we compare the density change calculated at the QM/MM and full QM levels in models of typical polar groups present at the QM/MM boundary in calculations on biomolecular systems. The QM treatment in each case uses density functional theory, with the popular B3LYP functional. The calculations enable us to assess whether the QM/MM method treats polarization correctly, and whether the QM/MM treatment introduces spurious charge leakage effects. We have also computed density changes integrated over individual atoms, by using changes in atomic charges calculated using the natural bond orbital (NBO) method (Reed *et al*. 1988). In contrast to the density plots, this gives quantitative information about electronic charge rearrangement on hydrogen bond formation.

For these calculations, we have evaluated the electron density in hydrogen-bonded dimers, and compared it with the density in the corresponding monomers, at the geometry of the dimer. For this purpose, we have used the geometry of the dimer optimized at the QM level for evaluating the density difference in both the QM and the QM/MM calculations. In principle, we could have used the QM/MM-optimized geometry for the latter comparison. We have optimized this geometry in every case and found that both the geometry and the bond energy derived from the QM/MM calculations are in reasonable agreement with the corresponding values obtained in all-QM calculations. However, there are some modest differences in geometry especially in the length of the hydrogen bond.1 It is preferable to avoid this source of bias when comparing the QM and QM/MM differences between the densities in the fragments and the dimer. This being said, we should point out that, for some of the complexes, the QM and QM/MM density differences were also computed using a QM/MM-optimized geometry, and the results are qualitatively very similar indeed to those reported here based on the use of a QM geometry.

## 2. Computational details

QM calculations were performed using the Jaguar v. 5.0 (Jaguar 2002) and Gaussian v. 03 (Frisch *et al*. 2004) programs. First, the structures of the hydrogen-bonded systems shown in scheme 1 were optimized at the B3LYP/6-311+G(d,p) level of QM theory. The performance of the B3LYP functional for hydrogen-bonded systems has been studied extensively (e.g. Rablen *et al*. 1998; Rappe & Bernstein 2000). Results with basis sets larger than 6-31+G(d) are typically comparable to results from accurate electron correlated molecular orbital methods (Sim *et al*. 1992; Rablen *et al*. 1998). The change in electron density due to hydrogen bonding was obtained by subtracting the sum of the total electron densities of the constituent monomers (using their geometries in the bonded complex) from the total electron density of the hydrogen-bonded system, Δ*ρ*_{QM}=*ρ*_{DA,QM}-(*ρ*_{D,QM}+*ρ*_{A,QM}), where D and A denote donor and acceptor monomers, respectively. The GaussView program was used to plot the difference in densities Δ*ρ*, with an isodensity value of 0.002 e bohr^{−3}. All density calculations and NBO calculations were performed at the geometry optimized at the full QM level (B3LYP/6-311+G(d,p)).

In the QM/MM modelling, QM calculations were performed using Jaguar v. 5.0 at the B3LYP/6-311+G(d,p) level. The Tinker (Ponder & Case 2003) molecular mechanics program with the CHARMM27 all-atom force field (MacKerell *et al*. 1998, 2000) was used to evaluate MM terms. The results from QM and MM calculations were combined using our own QM/MM program, QoMMMa (Harvey 2004). This program creates input for both programs and automatically extracts the required information from output. For all the hydrogen-bonded species studied, two separate QM/MM calculations were carried out. In the first, the hydrogen bond donor was treated using the QM method, with the acceptor treated using the MM method. In the second calculation, this was reversed (QM treatment of the acceptor and MM treatment of the donor). Note that, for these systems, the QM/MM boundary always separates species that are not covalently bonded, so there was no need to use link atoms. Full geometry optimization was carried out at the QM/MM level for each complex (see the electronic supplementary material).

To calculate the change in electron density due to hydrogen bonding at the QM/MM level Δ*ρ*_{QM/MM}, it is necessary to define a ‘total’ density of the hydrogen-bonded system at this level. We take this density *ρ*_{DA,QM/MM} as the sum of the densities obtained in two QM/MM calculations in which each monomer in turn was treated as the QM region. Hence, Δ*ρ*_{QM/MM}=(*ρ*_{D*A,QM/MM}+*ρ*_{DA*,QM/MM})−(*ρ*_{D,QM}+*ρ*_{A,QM}), where the asterisk denotes the partner treated at the QM level in a QM/MM calculation. These density difference calculations were carried out using the QM-optimized geometries of the complexes. The resulting density difference Δ*ρ*_{QM/MM} was then plotted and compared with Δ*ρ*_{QM}. As discussed in the text above, making this comparison with a consistent set of geometries (we chose the QM geometries) means that the comparison is not affected by the structural differences that would arise when using different levels of theory. We note that, in the QM/MM calculations, it is possible to define separate density differences for each of the constituent monomers, i.e. Δ*ρ*_{D,QM/MM}=*ρ*_{D*A,QM/MM}−*ρ*_{D,QM} and Δ*ρ*_{A,QM/MM}=*ρ*_{DA*,QM/MM}−*ρ*_{A,QM}. As this is not straightforwardly possible at the QM level, and as the contributions from Δ*ρ*_{D,QM/MM} and Δ*ρ*_{A,QM/MM} are easily seen in the figures below, we have not plotted these separate density differences here.

## 3. Results and discussion

The structures of the hydrogen-bonded systems studied in the present investigation are shown in scheme 1. The changes in total electron density Δ*ρ*_{QM} and Δ*ρ*_{QM/MM} upon hydrogen bond formation were calculated using the QM and QM/MM methods at the B3LYP/6-311+G(d,p) level of theory, and are shown in figures 1 and 2. Geometries of the hydrogen-bonded complexes optimized at the full QM level were used in all density calculations. Similar results were obtained when QM/MM-optimized structures were used. The reference density for each monomer was computed at the corresponding geometry in the hydrogen-bonded complex. Red regions denote depletion of electron density and blue ones correspond to a gain of electron density on formation of the hydrogen-bonded complex.

The main changes observed upon complex formation are a polarization of electron density within the hydrogen bond acceptor towards the region where the acceptor lone pair is situated, and a polarization of the hydrogen bond donor so as to increase the positive charge on the hydrogen. In all cases, corresponding zones of increase and decrease in density, respectively, are found in the appropriate parts of the complex. In turn, this polarization of the boundary region between the two fragments leads to changes in density elsewhere within a given fragment. For example, in cases where a water molecule accepts a hydrogen bond (figure 1*b–d*,*g*,*i*), the increase in electron density in the region above the molecular plane and near the oxygen atom is accompanied by a decrease in electron density on the hydrogen atoms, as hydrogen bonding increases the polarity of the O–H bonds. For a given X–H bond within a donor monomer where the H forms a hydrogen bond to an acceptor, a common pattern of four alternating regions of decrease and increase in density is observed. The first of these regions is the already noted decrease in density (increase in positive charge) near the H atom, and is followed by a region of increase in density along the X–H bond, close to the electronegative X atom, as the bond become more polarized towards X. There is then a region of decrease in density around the X atom as its electron density polarizes away from the partial negative charge of the acceptor. Finally, there is an increase in density on the distal side of the X atom, also corresponding to this polarization of the atomic density. As expected, in the hydrogen-bonded systems that include an ionic monomer (figures 1*e–g* and 2), the polarization effects are somewhat larger than in the systems consisting of neutral monomers.

As can be seen in figures 1 and 2, the difference densities calculated at the full QM level and using the QM/MM formalism are remarkably similar in all cases. The position and size of the regions of decreased and increased density are very consistent from one method to the other for every dimer considered. As we use the same isovalue to plot all of the surfaces in figures 1 and 2, small changes in the density difference could lead to significant changes in the appearance of the plots. The fact that only slight changes occur shows that the QM/MM calculations provide a good description of the electronic polarization effects. At the same time, it should of course be noted that small differences can be detected between the QM and QM/MM plots, especially when one of the monomers is charged.

These results can be related to energy decomposition analyses of Curutchet *et al*. (2003). For a set of hydrogen-bonded systems similar to that studied here, they found that, for systems consisting of neutral monomers, the polarization energy amounts to only 10 per cent of the electrostatic energy, but that in cases with charged monomers this proportion rises to approximately 40 per cent. They also observed a close agreement between the polarization energies calculated for each monomer through the energy decomposition analysis and through a QM/MM approach similar to that used in the present investigation.

In order to get quantitative information on intermolecular charge transfer upon hydrogen bonding, changes in atomic charge were calculated using NBO charges at the QM and QM/MM levels for selected hydrogen-bonded systems. The results are summarized in figure 3. While these values provide quantitative information that is not easily obtainable from the plots of figures 1 and 2, they also contain much less detail since some changes in electron density distribution occur within the region surrounding a given atom. The changes in atomic charges are similar for formamide with water and for the water dimer in both QM and QM/MM formalisms, except for slight differences for the atoms directly involved in hydrogen bonding. Even for the latter, the values from the QM and QM/MM calculations agree to within less than 0.01 arb. units.

In the case of the guanidinium–acetate ion pair, the computed changes in atomic charges obtained from full QM and QM/MM levels differ more, with a maximum difference of 0.1 arb. units. This is in part again due to the larger polarization effects in this charged system compared with the neutral–neutral cases mentioned above. However, part of the discrepancy is also due to the fact that, at the QM/MM level, the NBO charges on each fragment are calculated separately, and the sum of their changes for each fragment must be exactly zero. In other words, intermolecular charge transfer is not accounted for, at least in terms of the NBO charges. The QM-derived NBO charges involve significant charge transfer and this accounts for a large part of the discrepancy between the charge changes in the QM and QM/MM calculations in the ion pair case. It is interesting to note, however, that the density change maps shown in figure 2 are very similar for the QM and QM/MM calculations, so clearly some of the ‘charge transfer’ for the guanidinium–acetate ion pair is reproduced in the QM/MM calculation as a polarization of the hydrogen bond acceptor.

It is also important to study basis set effects on the change in electron density. The change in electron density due to hydrogen bond formation between guanidinium and acetate was calculated at the QM and QM/MM levels using three different basis sets, 6-31G(d), 6-311+G(d,p) and aug-cc-pvtz. The results are shown in figure 2. For all these calculations the geometry optimized using B3LYP/6-311+G(d,p) was used. Both in QM and QM/MM calculations the polarization at the carboxylate oxygens in acetate towards the guanidinium hydrogen is slightly increased with increasing basis set size. This extra polarization is expected in more flexible basis sets.

As mentioned above, there has been considerable concern about the possibility of charge leakage from the QM region to point charges on MM atoms when using the simple Hamiltonian term of equation (1.2). This ‘charge leakage’ phenomenon could also be described as being an unphysical polarization of the QM density. It has been found that the use of diffuse functions can sometimes lead to an SCF convergence problem in QM/MM calculations (Mulholland *et al*. 2000), which could be caused by such effects. Much effort has been put into developing screened point charges to avoid this pathological effect. By contrast, our calculations described above and other work with the QM/MM method (Harvey 2004; Bathelt *et al*. 2005*a*,*b*; Claeyssens *et al*. 2005; Mulholland 2005; Strickland *et al*. 2006) use simple, unscreened charges. To examine whether charge leakage occurs in these QM/MM calculations, we plot the density difference for the ion pair case of figure 2, but using only the QM/MM charge density in which the acetate is treated in the QM region (figure 4). Leakage is most likely to occur from the electron-rich acetate system towards the positive point charges used to represent the MM atoms of the guanidinium H-bond donor. The plot shows in both QM and QM/MM calculations the expected polarization of the electron density towards the guanidinium protons. This occurs to a somewhat greater extent using the larger basis. In any case, no significant build-up of density is observed on the guanidinium hydrogen atoms, suggesting that pathological charge leakage does not in fact occur. This is partly due to the localized nature of the basis set in these QM/MM calculations, as basis functions centred on the acetate oxygens do not describe charge density around the guanidinium protons effectively.2 Significant problems from leakage can be expected only when a very large (or plane wave) basis set is used, and a large positive MM charge is close to a QM region with significant electron density. However, this would be a very unusual arrangement in a standard QM/MM set-up for biomolecular systems. For more standard systems, the confining effect of the atom-centred basis set, and the low charges on typical acidic and basic protein side chains suggest that this problem is unlikely to be very serious. These results indicate that the hydrogen bond interactions are treated well at the QM/MM levels. Similarly, comparison of *ab initio* QM/MM and QM results (Ranaghan *et al*. 2004), and QM energy decomposition analysis for interactions along the reaction path of the enzyme chorismate mutase (Szefczyk *et al*. 2004) showed hydrogen bonds (including charged hydrogen bonds) to be treated well at the QM/MM levels.

## 4. Conclusions

We have explored the key features of standard QM/MM implementations (such as that used in our group; Harvey 2004) to model biologically important hydrogen bond interactions. A set of hydrogen-bonded complexes containing typical amino acid side chains has been studied using both fully QM (DFT) and hybrid QM/MM methods. The change in total electron density upon hydrogen bond formation is a sensitive measure of the polarization effects between the two partners. For each hydrogen-bonded pair, we have found that the density changes computed using QM and QM/MM methods are very similar, showing that the QM/MM method provides a good description of polarization effects. This is confirmed by calculating differences in NBO charges upon hydrogen bond formation.

In one case that might *a priori* be expected to lead to non-physical charge leakage from the QM region to the point charges on the adjacent MM atoms, no evidence for such an effect is apparent in the electron density change. This is because the atom-centred nature of the basis set makes charge transfer (as opposed to mere polarization) unfavourable. Charge leakage would be expected to occur only in cases where very large (or plane wave) basis sets are used and large MM point charges are present at the QM/MM boundary, but this should not be the case in most biomolecular QM/MM studies.

For reliable modelling of biological reactions using QM/MM methods, it is vital that the interaction between the reacting system (QM) and its environment (MM) are treated carefully. Hydrogen bonds are among the most important interactions in biological reactions, and the present results indicate that these interactions can be modelled well by QM/MM methods.

## Supplementary Material

Distance between hydrogen bond donor atom and hydrogen atom, R(D–H) and hydrogen bond angle Å(D–H⋯A): calculated at QM/MM level with CHARMM27, and at full QM level (B3LYP/6-311+G(d,p).

## Acknowledgements

The authors thank the EPSRC for financial support (grant number EP/C532910/1). J.I.M. thanks Eusko Jaurlaritza and MEC for funding. A.J.M. also thanks the BBSRC and the IBM High Performance Computing, Life Sciences Outreach Programme. J.N.H. and F.R.M. thank EPSRC and the Royal Society, respectively, for research fellowships.

## Footnotes

One contribution of 9 to a Theme Supplement ‘Biomolecular simulation’.

↵One contribution to this difference in geometry is the fact that the QM/MM calculations can be described as accounting for polarization effects

*twice*at the QM–MM interface, e.g. where a hydrogen bond donor is in the QM region and the corresponding acceptor is in the MM region. Including the MM point charges in the QM Hamiltonian polarizes the X–H bond of the donor, strengthening and shortening the distance between the two partners. At the same time, the van der Waals non-bonded parameters for the H-bond donor will have been chosen in the MM force field so as to reproduce the H-bond distance. Hence in cases such as these, the QM/MM geometry often involves a too close contact, e.g. for the water dimer and the water–*N*–methylacetamide cases, the H-bonding H–O distances are of, respectively, 1.71 and 1.69 Å with QM/MM versus 1.93 and 1.87 Å with QM. This problem can be avoided by using van der Waals parameters for the QM atoms that are specifically optimized for QM/MM calculations (e.g. Luque*et al*. 2000; Martin*et al*. 2002; Freindorf*et al*. 2005; MacKerell 2005). However, it should be noted that tests indicate that this effect does not usually have a large impact on calculated relative energies (see Riccardi*et al*. 2004).↵We note that the energy of an electron interacting with a point charge of +0.46 (as for a guanidinium hydrogen in the CHARMM force field), and described by the 6-31G basis functions of an oxygen atom centred at 1.6 Å from the point charge, is positive by 130 kcal mol

^{−1}. Using a larger basis set such as 6-311+G(2d) on the adjacent centre does lead to a negative energy of −34 kcal mol^{−1}, but even this is well above the basis limit for a charge of 0.46 of −67 kcal mol^{−1}.

- Received June 5, 2008.
- Accepted August 4, 2008.

- © 2008 The Royal Society