QM/MM calculations of kinetic isotope effects in the chorismate mutase active site

Kinetic isotope effects have been computed for the Claisen rearrangement of chorismate to prephenate in aqueous solution and in the active site of chorismate mutase from B. subtilus. These included primary 13C and 18O and secondary 3H effects for substitutions at the bond-making and bond-breaking positions. The initial structures of the putative stationary points on the potential energy surface, required for the calculations of isotope effects using the CAMVIB/CAMISO programs, have been selected from hybrid QM/MM molecular dynamical simulations using the DYNAMO program. Refinement of the reactant complex and transition-state structures has been carried out by means of AM1/CHARMM24/TIP3P calculations using the GRACE program, with full gradient relaxation of the position of > 5200 atoms for the enzymic simulations, and with a box containing 711 water molecules for the corresponding reaction in aqueous solution. Comparison of these results, and of gas phase calculations, with experimental data has shown that the chemical rearrangement is largely rate-determining for the enzyme mechanism. Inclusion of the chorismate conformational pre-equilibrium step in the modelled kinetic scheme leads to better agreement between recent experimental data and theoretical predictions. These results provide new information on an important enzymatic transformation, and the key factors responsible for the kinetics of its molecular mechanism are clarified. Treatment of the enzyme and/or solvent environment by means of a large and flexible model is absolutely essential for prediction of kinetic isotope effects.     

The rate-determining step in the rhodium-xantphos-catalysed hydroformylation of 1-octene

The rate-determining step in the hydroformylation of 1-octene, catalysed by the rhodium-Xantphos catalyst system, was determined by using a combination of experimentally determined (1)H/(2)H and (12)C/(13)C kinetic isotope effects and a theoretical approach. From the rates of hydroformylation and deuterioformylation, a small (1)H/(2)H isotope effect of 1.2 was determined for the hydride moiety of the rhodium catalyst. (12)C/(13)C isotope effects of 1.012(1) and 1.012(3) for the alpha-carbon and beta-carbon atoms of 1-octene were determined, respectively. Both quantum mechanics/molecular mechanics (QM/MM) and full quantum mechanics calculations were carried out on the key catalytic steps, for "real-world" ligand systems, to clarify whether alkene coordination or hydride migration is the rate-determining step. Our calculations (21.4 kcal mol(-1)) quantitatively reproduce the experimental energy barrier for CO dissociation (20.1 kcal mol(-1)) starting at the (bisphosphane)RhH(CO)(2) resting state. The barrier for hydride migration lies 3.8 kcal mol(-1) higher than the barrier for CO dissociation (experimentally determined trend approximately 3 kcal mol(-1)). The computed (1)H/(2)H and (12)C/(13)C kinetic isotope effects corroborate the results of the energy analysis.     

Dihydrogen to dihydride isomerization mechanism in [(C5Me5)FeH2(Ph2PCH2CH2PPh2)]+ through the experimental and theoretical analysis of kinetic isotope effects

The isomerization of complex [Cp*Fe(dppe)(eta2-H2)]+, generated in situ by low-temperature protonation of Cp*Fe(dppe)H with either HBF4 or CF3COOH, to the dihydride tautomer trans-[Cp*Fe(dppe)(H)2]+ is irreversible and follows first-order kinetics in the -10 to +15 degrees C range with Delta H double dagger = 21.6 +/- 0.8 kcal mol(-1) and DeltaS double dagger = 5 +/- 3 eu. The isomerization rate constant is essentially independent of the nature and quantity of a strong acid. Density functional theory (DFT) calculations on various models, including the complete system at both the quantum mechanics/molecular mechanics (QM/MM) and full QM levels, probe the relative importance of steric and electronic effects for the relative stability of the nonclassical and classical isomers and identify two likely isomerization mechanisms: a "direct" pathway involving simultaneous H-H bond breaking and cis-trans isomerization and a "via Cp" pathway involving agostic C5Me5H intermediates. Both pathways are characterized by activation energies in close correspondence with the experimental value (21.3 and 22.2 kcal mol(-1), respectively). Further kinetic studies were carried out for the Cp*Fe(dppe)H + CF3COOD and Cp*Fe(dppe)D + CF3COOD systems at 273 K. The [Cp*Fe(dppe)(eta2-HD)]+ complex establishes a very rapid isotope redistribution equilibrium with the eta2-H2 and eta2-D2 analogues. The equilibrium constant value (K = 3.3 +/- 0.3) indicates a significant equilibrium isotope effect. Simulation of the rate data provides access to the individual isomerization rate constants kHH, kHD, and kDD for the three isotopomers, yielding kinetic isotope effects: kHH/kHD = 1.24 +/- 0.01 and kHD/kDD = 1.58 +/- 0.01 (and, consequently, kHH/kDD = 1.96 +/- 0.02). The analysis of the DFT-calculated frequencies, using the [Cp*Fe(dhpe)H2]+ model system, for the [Cp*Fe(dhpe)(eta2-XY)]+ isotopomers as well as transition states for the "direct" (TSdir) and "via Cp" (TSrot) pathways (X = H, D) allowed the computation of the expected isotope effects. A comparison with the experiment strongly suggests that the mechanism occurs via the "direct" pathway for the present system, although the small difference in the calculated energy barriers suggests that the "via Cp" pathway may be preferred in other cases.     

Exploring SCC-DFTB paths for mapping QM/MM reaction mechanisms

A new first-order procedure for locating transition structures (TS) that employs hybrid quantum mechanical/molecular mechanical (QM/MM) potentials has been developed. This new technique (RPATh+RESD) combines the replica path method (RPATh) and standard reaction coordinate driving (RCD) techniques in an approach that both efficiently determines reaction barriers and successfully eliminates two key weaknesses of RCD calculations (i.e., hysteresis/discontinuities in the path and the sequential nature of the RCD procedure). In addition, we have extended CHARMM's QM/MM reaction pathway methods, the RPATh and nudged elastic band (NEB) methods, to incorporate SCC-DFTB wave functions. This newly added functionality has been applied to the chorismate mutase-catalyzed interconversion of chorismate to prephenate, which is a key step in the shikimate pathway of bacteria, fungi, and other higher plants. The RPATh+RESD barrier height (DeltaE=5.7 kcal/mol) is in good agreement with previous results from full-energy surface mapping studies (Zhang, X.; Zhang, X.; Bruice, T. C. Biochemistry 2005, 44, 10443-10448). Full reaction paths were independently mapped with RPATh and NEB methods and showed good agreement with the final transition state from the RPATh+RESD "gold standard" and previous high-level QM/MM transition states (Woodcock, H. L.; Hodoscek, M.; Gilbert, T. B.; Gill, P. M. W.; Schaefer, H. F.; Brooks, B. R. J. Comput. Chem. 2007, 28, 1485-1502). The SCC-DFTB TS geometry most closely approximates the MP2/6-31+G(d) QM/MM result. However, the barrier height is underestimated and possibly points to an area for improvement in SCC-DFTB parametrization. In addition, the steepest descents (SD) minimizer for the NEB method was modified to uncouple the in-path and off-path degrees of freedom during the minimization, which significantly improved performance. The convergence behavior of the RPATh and NEB was examined for SCC-DFTB wave functions, and it was determined that, in general, both methods converge at about the same rate, although the techniques used for convergence may be different. For instance, RPATh can effectively use the adopted basis Newton-Raphson (ABNR) minimizer, where NEB seems to require a combination of SD and ABNR.     

Interfacing Q-Chem and CHARMM to perform QM/MM reaction path calculations

A hybrid quantum mechanical/molecular mechanical (QM/MM) potential energy function with Hartree-Fock, density functional theory (DFT), and post-HF (RIMP2, MP2, CCSD) capability has been implemented in the CHARMM and Q-Chem software packages. In addition, we have modified CHARMM and Q-Chem to take advantage of the newly introduced replica path and the nudged elastic band methods, which are powerful techniques for studying reaction pathways in a highly parallel (i.e., parallel/parallel) fashion, with each pathway point being distributed to a different node of a large cluster. To test our implementation, a series of systems were studied and comparisons were made to both full QM calculations and previous QM/MM studies and experiments. For instance, the differences between HF, DFT, MP2, and CCSD QM/MM calculations of H2O...H2O, H2O...Na+, and H2O...Cl- complexes have been explored. Furthermore, the recently implemented polarizable Drude water model was used to make comparisons to the popular TIP3P and TIP4P water models for doing QM/MM calculations. We have also computed the energetic profile of the chorismate mutase catalyzed Claisen rearrangement at various QM/MM levels of theory and have compared the results with previous studies. Our best estimate for the activation energy is 8.20 kcal/mol and for the reaction energy is -23.1 kcal/mol, both calculated at the MP2/6-31+G(d)//MP2/6-31+G(d)/C22 level of theory.     

Substrate hydroxylation in methane monooxygenase: quantitative modeling via mixed quantum mechanics/molecular mechanics techniques

Using broken-symmetry unrestricted density functional theory quantum mechanical (QM) methods in concert with mixed quantum mechanics/molecular mechanics (QM/MM) methods, the hydroxylation of methane and substituted methanes by intermediate Q in methane monooxygenase hydroxylase (MMOH) has been quantitatively modeled. This protocol allows the protein environment to be included throughout the calculations and its effects (electrostatic, van der Waals, strain) upon the reaction to be accurately evaluated. With the current results, recent kinetic data for CH3X (X = H, CH3, OH, CN, NO2) substrate hydroxylation in MMOH (Ambundo, E. A.; Friesner, R. A.; Lippard, S. J. J. Am. Chem. Soc. 2002, 124, 8770-8771) can be rationalized. Results for methane, which provide a quantitative test of the protocol, including a substantial kinetic isotope effect (KIE), are in reasonable agreement with experiment. Specific features of the interaction of each of the substrates with MMO are illuminated by the QM/MM modeling, and the resulting effects upon substrate binding are quantitatively incorporated into the calculations. The results as a whole point to the success of the QM/MM methodology and enhance our understanding of MMOH catalytic chemistry. We also identify systematic errors in the evaluation of the free energy of binding of the Michaelis complexes of the substrates, which most likely arise from inadequate sampling and/or the use of harmonic approximations to evaluate the entropy of the complex. More sophisticated sampling methods will be required to achieve greater accuracy in this aspect of the calculation.     

Combined QM/MM and path integral simulations of kinetic isotope effects in the proton transfer reaction between nitroethane and acetate ion in water

An integrated Feynman path integral-free energy perturbation and umbrella sampling (PI-FEP/UM) method has been used to investigate the kinetic isotope effects (KIEs) in the proton transfer reaction between nitroethane and acetate ion in water. In the present study, both nuclear and electronic quantum effects are explicitly treated for the reacting system. The nuclear quantum effects are represented by bisection sampling centroid path integral simulations, while the potential energy surface is described by a combined quantum mechanical and molecular mechanical (QM/MM) potential. The accuracy essential for computing KIEs is achieved by a FEP technique that transforms the mass of a light isotope into a heavy one, which is equivalent to the perturbation of the coordinates for the path integral quasiparticle in the bisection sampling scheme. The PI-FEP/UM method is applied to the proton abstraction of nitroethane by acetate ion in water through molecular dynamics simulations. The rule of the geometric mean and the Swain-Schaad exponents for various isotopic substitutions at the primary and secondary sites have been examined. The computed total deuterium KIEs are in accord with experiments. It is found that the mixed isotopic Swain-Schaad exponents are very close to the semiclassical limits, suggesting that tunneling effects do not significantly affect this property for the reaction between nitroethane and acetate ion in aqueous solution.     

Isotope effects on the enzymatic and nonenzymatic reactions of chorismate

The important biosynthetic intermediate chorismate reacts thermally by two competitive pathways, one leading to 4-hydroxybenzoate via elimination of the enolpyruvyl side chain, and the other to prephenate by a facile Claisen rearrangement. Measurements with isotopically labeled chorismate derivatives indicate that both are concerted sigmatropic processes, controlled by the orientation of the enolpyruvyl group. In the elimination reaction of [4-2H]chorismate, roughly 60% of the label was found in pyruvate after 3 h at 60 degrees C. Moreover, a 1.846 +/- 0.057 2H isotope effect for the transferred hydrogen atom and a 1.0374 +/- 0.0005 18O isotope effect for the ether oxygen show that the transition state for this process is highly asymmetric, with hydrogen atom transfer from C4 to C9 significantly less advanced than C-O bond cleavage. In the competing Claisen rearrangement, a very large 18O isotope effect at the bond-breaking position (1.0482 +/- 0.0005) and a smaller 13C isotope effect at the bond-making position (1.0118 +/- 0.0004) were determined. Isotope effects of similar magnitude characterized the transformations catalyzed by evolutionarily unrelated chorismate mutases from Escherichia coli and Bacillus subtilis. The enzymatic reactions, like their solution counterpart, are thus concerted [3,3]-sigmatropic processes in which C-C bond formation lags behind C-O bond cleavage. However, as substantially larger 18O and smaller 13C isotope effects were observed for a mutant enzyme in which chemistry is fully rate determining, the ionic active site may favor a somewhat more polarized transition state than that seen in solution.     

Dependence of transition state structure on substrate: the intrinsic C-13 kinetic isotope effect is different for physiological and slow substrates of the ornithine decarboxylase reaction because of different hydrogen bonding structures

Ornithine decarboxylase is the first and the rate-controlling enzyme in polyamine biosynthesis; it decarboxylates l-ornithine to form the diamine putrescine. We present calculations performed using a combined quantum mechanical and molecular mechanical (QM/MM) method with the AM1 semiempirical Hamiltonian for the wild-type ornithine decarboxylase reaction with ornithine (the physiological substrate) and lysine (a "slow" substrate) and for mutant E274A with ornithine substrate. The dynamical method is variational transition state theory with quantized vibrations. We employ a single reaction coordinate equal to the carbon-carbon distance of the dissociating bond, and we find a large difference between the intrinsic kinetic isotope effect for the physiological substrate, which equals 1.04, and that for the slow substrate, which equals 1.06. This shows that, contrary to a commonly accepted assumption, kinetic isotope effects on slow substrates are not always good models of intrinsic kinetic isotope effects on physiological substrates. Furthermore, analysis of free-energy-based samples of transition state structures shows that the differences in kinetic isotope effects may be traced to different numbers of hydrogen bonds at the different transition states of the different reactions.     

Ensemble-averaged QM/MM kinetic isotope effects for the S(N)2 reaction of cyanide anions with chloroethane in DMSO solution

The existence of solvent fluctuations leads to populations of reactant-state (RS) and transition-state (TS) configurations and implies that property calculations must include appropriate averaging over distributions of values for individual configurations. Average kinetic isotope effects 〈KIE〉 for NC(-) + EtCl → NCEt + Cl(-) in DMSO solution at 30?°C are best obtained as the ratio 〈f(RS)〉/〈f(TS)〉 of isotopic partition function ratios separately averaged over all RS and TS configurations. In this way the hybrid AM1/OPLS-AA potential yields 〈KIE〉 values for all six isotopic substitutions (2° α-(2)H(2), 2° β-(2)H(3), α-(11)C/(14)C, leaving group (37)Cl, and nucleophile (13)C and (15)N) for this reaction in the correct direction as measured experimentally. These thermally-averaged calculated KIEs may be compared meaningfully with experiment, and only one of them differs in magnitude from the experimental value by more than one standard deviation from the mean. This success contrasts with previous KIE calculations based upon traditional methods without averaging. The isotopic partition function ratios are best evaluated using all (internal) vibrational and (external) librational frequencies obtained from Hessians determined for subsets of atoms, relaxed to local minima or saddle points, within frozen solvent environments of structures sampled along molecular dynamics trajectories for RS and TS. The current method may perfectly well be implemented with other QM or QM/MM methods, and thus provides a useful tool for investigating KIEs in relation to studies of chemical reaction mechanisms in solution or catalyzed by enzymes.     

Geometry optimization based on linear response free energy with quantum mechanical/molecular mechanical method: applications to Menshutkin-type and Claisen rearrangement reactions in aqueous solution

The authors present a method based on a linear response theory that allows one to optimize the geometries of quantum mechanical/molecular mechanical (QM/MM) systems on the free energy surfaces. Two different forms of linear response free energy functionals are introduced, and electronic wave functions of the QM region, as well as the responses of electrostatic and Lennard-Jones potentials between QM and MM regions, are self-consistently determined. The covariant matrix relating the QM charge distribution to the MM response is evaluated by molecular dynamics (MD) simulation of the MM system. The free energy gradients with respect to the QM atomic coordinates are also calculated using the MD trajectory results. They apply the present method to calculate the free energy profiles of Menshutkin-type reaction of NH3 with CH3Cl and Claisen rearrangement of allyl vinyl ether in aqueous solution. For the Menshutkin reaction, the free energy profile calculated with the modified linear response free energy functional is in good agreement with that by the free energy perturbation calculations. They examine the nonequilibrium solvation effect on the transmission coefficient and the kinetic isotope effect for the Claisen rearrangement.     

Characteristics of CO3(2-)-water hydrogen bonds in aqueous solution: insights from HF/MM and B3LYP/MM MD simulations

Two combined quantum mechanics/molecular mechanics (QM/MM) molecular dynamics simulations, namely HF/MM and B3LYP/MM, have been performed to investigate the local hydration structure and dynamics of carbonate (CO(3)(2-)) in dilute aqueous solution. With respect to the QM/MM scheme, the QM region, which contains the CO(3)(2-) and its surrounding water molecules, was treated at HF and B3LYP levels of accuracy, respectively, using the DZV+ basis set, while the rest of the system is described by classical MM potentials. For both the HF/MM and B3LYP/MM simulations, it is observed that the hydrogen bonds between CO(3)(2-) oxygens and their nearest-neighbor waters are relatively strong, i.e., compared to water-water hydrogen bonds in the bulk, and that the first shell of each CO(3)(2-) oxygen atom somewhat overlaps with the others, which allows migration of water molecules among the coordinating sites to exist. In addition, it is observed that first-shell waters are either "loosely" or "tightly" bound to the respective CO(3)(2-) oxygen atoms, leading to large fluctuations in the number of first-shell waters, ranging from 1 to 6 (HF/MM) and 2 to 7 (B3LYP/MM), with the prevalent value of 3. Upon comparing the HF and B3LYP methods in describing this hydrated ion, the latter is found to overestimate the hydrogen-bond strength in the CO(3)(2-)-water complexes, resulting in a slightly more compact hydration structure at each of the CO(3)(2-) oxygens.     

An application of the novel quantum mechanical/molecular mechanical method combined with the theory of energy representation: An ionic dissociation of a water molecule in the supercritical water

A novel quantum chemical approach recently developed has been applied to an ionic dissociation of a water molecule (2H(2)O-->H(3)O(+)+OH(-)) in ambient and supercritical water. The method is based on the quantum mechanical/molecular mechanical (QM/MM) simulations combined with the theory of energy representation (QM/MM-ER), where the energy distribution function of MM solvent molecules around a QM solute serves as a fundamental variable to determine the hydration free energy of the solute according to the rigorous framework of the theory of energy representation. The density dependence of the dissociation free energy in the supercritical water has been investigated for the density range from 0.1 to 0.6 g/cm(3) with the temperature fixed at a constant. It has been found that the product ionic species significantly stabilizes in the high density region as compared with the low density. Consequently, the dissociation free energy decreases monotonically as the density increases. The decomposition of the hydration free energy has revealed that the entropic term (-TDeltaS) strongly depends on the density of the solution and dominates the behavior of the dissociation free energy with respect to the variation of the density. The increase in the entropic term in the low density region can be attributed to the decrease in the translational degrees of freedom brought about by the aggregation of solvent water molecules around the ionic solute.     

Reaction path potential for complex systems derived from combined ab initio quantum mechanical and molecular mechanical calculations

Combined ab initio quantum mechanical and molecular mechanical calculations have been widely used for modeling chemical reactions in complex systems such as enzymes, with most applications being based on the determination of a minimum energy path connecting the reactant through the transition state to the product in the enzyme environment. However, statistical mechanics sampling and reaction dynamics calculations with a combined ab initio quantum mechanical (QM) and molecular mechanical (MM) potential are still not feasible because of the computational costs associated mainly with the ab initio quantum mechanical calculations for the QM subsystem. To address this issue, a reaction path potential energy surface is developed here for statistical mechanics and dynamics simulation of chemical reactions in enzymes and other complex systems. The reaction path potential follows the ideas from the reaction path Hamiltonian of Miller, Handy and Adams for gas phase chemical reactions but is designed specifically for large systems that are described with combined ab initio quantum mechanical and molecular mechanical methods. The reaction path potential is an analytical energy expression of the combined quantum mechanical and molecular mechanical potential energy along the minimum energy path. An expansion around the minimum energy path is made in both the nuclear and the electronic degrees of freedom for the QM subsystem internal energy, while the energy of the subsystem described with MM remains unchanged from that in the combined quantum mechanical and molecular mechanical expression and the electrostatic interaction between the QM and MM subsystems is described as the interaction of the MM charges with the QM charges. The QM charges are polarizable in response to the changes in both the MM and the QM degrees of freedom through a new response kernel developed in the present work. The input data for constructing the reaction path potential are energies, vibrational frequencies, and electron density response properties of the QM subsystem along the minimum energy path, all of which can be obtained from the combined quantum mechanical and molecular mechanical calculations. Once constructed, it costs much less for its evaluation. Thus, the reaction path potential provides a potential energy surface for rigorous statistical mechanics and reaction dynamics calculations of complex systems. As an example, the method is applied to the statistical mechanical calculations for the potential of mean force of the chemical reaction in triosephosphate isomerase.     

First X-ray characterization and theoretical study of pi-alkyne,alkynyl-hydride,and vinylidene isomers for the same transition metal fragment [Cp*Ru(PEt3)2]+

The reaction of the chloro-complex [CpRuCl(PEt(3))(2)] with acetylene gas in methanol gave the pi-alkyne complex [CpRu(eta(2)-HCtbd1;CH)(PEt(3))(2)][BPh(4)] (1), which has been structurally characterized by X-ray analysis. The alkyne complex undergoes spontaneous isomerization even at low temperature, yielding the metastable alkynyl-hydride complex [CpRu(H)(Ctbd1;CH)(PEt(3))(2)][BPh(4)] (2), as the result of the oxidative addition of the alkyne C-H bond. This compound has also been structurally characterized despite it tautomerizes spontaneously into the stable primary vinylidene [CpRu(=C=CH(2))(PEt(3))(2)][BPh(4)] (3). This species has been alternatively prepared by a two-step deprotonation/protonation synthesis from the pi-alkyne complex. Moreover, the reaction of the initial chloro-complex with monosubstituted alkynes HCtbd1;CR (R = SiMe(3), Ph, COOMe, (t)Bu) has been studied without detection of pi-alkyne intermediates. Instead of this, alkynyl-hydride complexes were obtained in good yields. They also rearrange to the corresponding substituted vinylidenes. In the case of R = SiMe(3), the isomerization takes place followed by desilylation, yielding the primary vinylidene complex. X-ray crystal structures of the vinylidene complexes [CpRu(=C=CH(2))(PEt(3))(2)][BPh(4)] (3) and [CpRu(=C=CHCOOMe)(PEt(3))(2)][BPh(4)] (10) have also been determined. Both, full ab initio and quantum mechanics/molecular mechanics (QM/MM) calculations were carried out, respectively, on the model system [CpRu(C(2)H(2))(PH(3))(2)](+) (A) and the real complex [CpRu(C(2)H(2))(PEt(3))(2)](+) (B) to analyze the steric and electronic influence of ligands on the structures and relative energies of the three C(2)H(2) isomers. QM/MM calculations have been employed to evaluate the role of the steric effects of real ligands, whereas full ab initio energy calculations on the optimized QM/MM model have allowed recovering the electronic effects of ligands. Additional pure quantum mechanics calculations on [CpRu(C(2)H(2))(PH(3))(2)](+) (C) and [CpRu(C(2)H(2))(PMe(3))(2)](+) (D) model systems have been performed to analyze in more detail the effects of different ligands. Calculations have shown that the steric effects induced by the presence of bulky substituents in phosphine ligand are responsible for experimentally observed alkyne distortion and for relative destabilization of the alkyne isomer. Moreover, increasing the phosphine basicity and sigma donor capabilities of ligands causes a relative stabilization of an alkynyl-hydride isomer. The combination of both steric and electronic effects, makes alkyne and alkynyl-hydride isomers to be close in energy, leading to the isolation of both complexes.     

Quantum mechanics/molecular mechanics electrostatic embedding with continuous and discrete functions

A quantum mechanics/molecular mechanics (QM/MM) implementation that uses the Gaussian electrostatic model (GEM) as the MM force field is presented. GEM relies on the reproduction of electronic density by using auxiliary basis sets to calculate each component of the intermolecular interaction. This hybrid method has been used, along with a conventional QM/MM (point charges) method, to determine the polarization on the QM subsystem by the MM environment in QM/MM calculations on 10 individual H(2)O dimers and a Mg(2+)-H(2)O dimer. We observe that GEM gives the correct polarization response in cases when the MM fragment has a small charge, while the point charges produce significant over-polarization of the QM subsystem and in several cases present an opposite sign for the polarization contribution. In the case when a large charge is located in the MM subsystem, for example, the Mg(2+) ion, the opposite is observed at small distances. However, this is overcome by the use of a damped Hermite charge, which provides the correct polarization response.