Ab initio QM/MM molecular dynamics simulation of preferential K+ solvation in aqueous ammonia solution

A combined ab initio quantum mechanical/molecular mechanical (QM/MM) molecular dynamics simulation has been performed to investigate the solvation structure of K+ in 18.4% aqueous ammonia solution.

The chemically most relevant region, the first solvation sphere of K+, was treated by Born–Oppenheimer ab initio quantum mechanics using LANL2DZ basis sets, while the rest of the system was described by classical pair potentials.

Within the first solvation shell of K+, the QM/MM simulation reveals a polyhedral structure with an average coordination number of 7.6, consisting of 6.7 water and 0.9 ammonia molecules, compared to the corresponding value of 8.7 composed of 5.3 water and 3.4 ammonia molecules obtained by classical pair potential simulation.

The QM/MM results, in contrast to the classical simulation, clearly indicate a preference for water ligands and a higher flexibility of ligand arrangements in the first solvation shell of the ion.

The preference for ligands is discussed on the basis of detailed simulation results.

In addition, a “structure-breaking” behavior of the ion is well recognized by the detailed analysis on ligand exchange processes and the mean residence times of the ligands surrounding the ion.


Since ion solvation processes are of great importance in a wide variety of systems, hydrated ions have become an increasingly important subject of both experimental and theoretical studies.1–9

In general, comparison between experiments and theories is not always straightforward because most of the experimental methods for structural analysis have to be performed with solutions of relatively high concentrations, while the theoretical approaches mostly refer to very dilute solutions.

Theoretical methods, in particular Monte Carlo (MC) and molecular dynamics (MD) simulations, can provide valuable complementary information not accessible to experimental approaches both in the characterization of ion–solvent complexes and in the specific mechanism of the involved interactions.

During the past decades, numerous MC and MD simulations have been carried out to evaluate hydration structures, hydration energies, mobilities, librational and vibrational spectra and other properties.

However, most of the simulation works had relied on classical molecular mechanical potentials, and hence, deviation of the results from experiments was often found, strongly depending on the type and quality of the potential models employed in the simulations.10,11

To provide a realistic view of the solvated ion properties, it has been shown that the model interaction potential must include polarizability and many-body nonadditive contributions.12

With regard to this point, various schemes have been proposed, for example, the polarizable continuum model (PCM),13 which incorporates the many-body interactions in an average way, or a direct approach by calculating the energy hypersurface of many-body interactions and then fitting them to an analytical function.14,15

Both exemplary models can provide significant improvement of the results, however, for the PCM model, a major weakness is that it can not reproduce specific interaction with the surrounding solvent, such as hydrogen bonds.16

For the second model, construction of many-body potentials is rather difficult, and is hardly feasible for large molecular systems because of their complicated orientation dependence.15,17

To describe molecular interactions correctly, more sophisticated simulation technique incorporating quantum mechanical algorithms, namely the combined quantum mechanical and molecular mechanical (QM/MM) method has been introduced.18–20

Several hybrid QM/MM models combine either semiempirical,18,20–22 density functional,23 valence bond,24,25 or even ab initio Hartree–Fock19,26 methodology with commonly used force fields.

In recent years, a “Born–Oppenheimer ab initio QM/MM dynamics” has been proposed and successfully applied for elucidating structural and dynamical properties of various ions in solutions.27–41

This technique treats the active-site region, i.e. the solvation shell around the ion, quantum mechanically, while the environment consisting of further solvent molecules is described by classical molecular mechanical potentials.

By this approach, the complicated many-body contributions as well as the polarization effects within the solvation sphere of the ion can be reliably included.

The QM/MM results have clearly shown the important role of non-additive contributions and that inclusion of higher-order interaction terms is essential for the correct description of the solvated ions, even for monovalent ions in which many-body interactions could be expected to be weaker than in the case of di- and trivalent ions.29,30,32,33,38

For ion solvation in mixed solvent systems, Born–Oppenheimer ab initio QM/MM molecular dynamics simulations have been successfully applied to Li+,32 Na+,33 Mg2+34 and Ca2+35 in aqueous ammonia solution.

The results from QM/MM simulations have provided more reliable geometrical arrangements of the solvated ion structures as well as ligand preferences in good agreement with the qualitative expectation according to Lewis acid–base interactions.42–44

In comparison to the results obtained by classical pair potential simulations, the QM/MM simulations have indicated that the pairwise additive potentials, which lack suitable parameterization, are obviously insufficient for the description of the structural properties of solvated ions in such a solvent mixture.

In the present work, a Born–Oppenheimer ab initio QM/MM dynamics simulation was performed to characterize the preferential solvation of K+ in aqueous ammonia solution.

The biological importance of K+ is well known, e.g. in neuron signaling and osmotic stability of cells.45,46

In ion channels of the nervous system, the binding of potassium to ligands is an essential part of the mechanism for message transport.47

In order to understand such biological processes, detailed knowledge of ion structures in solution is required.

The solvation structure of K+ in water and in ammonia has been studied by X-ray and neutron diffraction experiments4,48–52 as well as by theoretical investigations.8,9,30,53–55

Experiments reported hydration numbers of K+ in the range 4–7,48–51 while a coordination number of 6 was found in ammonia.52

In addition, the experiments have shown that K+ acts as “structure-breaking” ion, i.e., the presence of this ion in the solvent can be regarded rather as a perturbation of the solvent network structures.

In theoretical investigations, classical and QM/MM simulations predicted hydration numbers of K+ in the range of 5–108,9,53–55 and 8.3,30 respectively.

Our previous QM/MM simulation30 has shown that the structure-breaking behavior of this ion is only reproduced, if the first hydration shell is treated quantum mechanically.

For the coordination number of K+ in ammonia, a classical MD simulation reported an average value of 8..756

Since no structural data for K+ in aqueous ammonia solution have been reported so far, structural information for this particular ion solvation by means of ab initio QM/MM dynamics could be expected to provide useful information regarding the functionality of this metal ion in biological systems with different coordination sites.


Based on Born–Oppenheimer ab initio dynamics, the system is partitioned into a part described by quantum mechanics (QM) and another part treated by means of molecular mechanical potentials (MM).

The total interaction energy of the system can be defined as Etotal = 〈ΨQM|ĤQM〉 + EMM + EQM–MM, where 〈ΨQM|ĤQM〉 refers to interactions within the solvation sphere of K+, treated by Born–Oppenheimer ab initio Hartree–Fock quantum mechanical calculations, while EMM and EQM–MM represent interactions within the MM region and between QM and MM regions, respectively, and are described by classical pair potentials.

Within the QM region, the accurate description of molecular interactions at Hartee–Fock level obviously depends only on the basis set quality.

The LANL2DZ basis sets57,58 were selected as a reasonable choice compromising between the quality of simulation results and the requirement of CPU time.28–30,32–35

Effects of electron correlation may play some role in the interaction between the ion and solvent especially when the number of solvent molecules in the QM region increases.

However, treatment at the correlated level is very time-consuming, and thus beyond the current computational feasibility.

In our previous QM/MM simulation,38 these effects were found to be rather marginal even for anion–water interactions, when the stabilization energies obtained by HF single point calculations of the anions plus their first hydration shells for several selected configurations were compared to those of MP2 calculations.

In the QM/MM procedure, an ab initio calculation was performed in each simulation step, providing quantum mechanical forces to be incorporated into the total force of the system.

Since the exchange of solvent molecules between the QM and MM regions can occur frequently during the simulation, the forces acting on each particle in the system are switched according to the defined region upon entering or leaving of the ligand and can be defined as: Fi = Sm(r)FQM + (1 – Sm(r))FMM, where FQM and FMM are quantum mechanical and molecular mechanical forces, respectively.

Sm(r) is a smoothing function,59 where r1 and r0 are the distances characterizing the beginning and the end of the QM region.

Sm(r) ensures a continuous change of forces in the transition between QM and MM regions.

For interactions inside the MM and between QM and MM regions, flexible models, which describe inter- and intramolecular interactions, were employed for water60,61 and ammonia.62

The use of flexible models is to be favored over any of the popular rigid water and ammonia models, in order to ensure compatibility and a smooth transition, when solvent molecules move from the QM region with their full flexibility of ligand molecules to the MM region.

The pair potential function for water–ammonia interactions was adopted from Tanabe and Rode.63

The pair potential function for K+–H2O interactions was obtained from our previous work,30 and the pair potential function for K+–NH3 interactions was newly constructed in the present work using DZV+(d,p) basis set for NH357 and Los Alamos ECP plus DZ basis set for K+.58

These basis sets are consistent with the ones used in the construction of the K+–H2O potential.

The 1500 Hartree–Fock interaction energy points for various K+–NH3 configurations, obtained from Gaussian9864 calculations, were fitted to the analytical form of where A, B, C and D are the fitting parameters (see Table 1), ric denotes the distances between K+ and the i-th atom of ammonia and q are the atomic net charges.

The charges on K+, N and H of ammonia were set to 1.0, –0.8022 and 0.2674, respectively.

A classical molecular dynamics simulation using pair potentials was performed first, then a combined QM/MM molecular dynamics simulation starting from the equilibrium configuration obtained by the classical simulation was carried out.

For the QM/MM simulation, undoubtedly, the QM region is the most expensive computational part and hence the size of it has to be a compromise between the reliability of results and the available computational resources.

In general, the suitable QM size is the sphere including all ligands within the complete first solvation sphere of the ions.

A too small size of the QM region can lead to artifacts, which, however, disappear, if the choice of the diameter guarantees that all molecules of the first solvation shell plus eventually exchanging ligands are included in this region.27,67

This is always achieved, when the QM limit is set to the beginning or even slightly inside the second solvation shell.

In this case, also an increase of width of the transition region does not make any difference for the results.

As can be seen from the K–(N + O) radial distribution function (RDF) obtained by the classical simulation (Fig. 1), the first minimum of the RDF peak is located around 3.8 Å.

Within this region, there are about 6–8 water and 3–4 ammonia molecules, which seem rather convenient to perform the QM/MM simulation.

However, as our previous QM/MM simulation for K+ in water30 had predicted a somewhat larger hydration shell compared to the one obtained by classical pair potential simulation, a larger QM region with diameter of 8.4 Å was employed.

This QM size was assumed to be large enough to ensure a smooth convergence of the quantum mechanical forces to the pair potential forces beyond the QM region.

There are about 8–10 water and 3–4 ammonia molecules located within this region, leading to 8–12 min for computing the forces in each QM/MM step on a DEC Alpha XP1000 workstation.

The classical and combined QM/MM simulations were performed in a canonical ensemble at 293 K. This canonical ensemble has been realized by coupling to an external temperature bath.

The time step was set to 0.2 fs, which allows for explicit movement of hydrogen atoms of water and ammonia.

The basic cube, with a length of 18.56 Å, contained one K+, 37 ammonia and 163 water molecules, assuming the experimental density of 18.4% aqueous ammonia solution at the given temperature (0.9398 g cm–3).

The reaction-field method was employed for the treatment of long-range interactions.65

This method also implies a compensation of the electrical non-neutrality of the basic box.

The classical pair potential simulation started from an equilibrium configuration obtained from our previous QM/MM simulation for Ca2+ in aqueous ammonia solution.35

The system was equilibrated for 80 000 steps, and the simulation was continued for another 100 000 steps for collecting configurations every 10th step.

As the equilibrium configuration of a QM/MM simulation is usually not much different from that of the classically treated system, and differences refer mostly to the small subsystem of the ion and its environment, a time-span of 2–3 ps is sufficient to achieve re-equilibration, even for systems with stronger ion–solvent interactions than the one investigated here, as can be seen from the behavior of velocities and RDFs after that period.34,37,41

The QM/MM simulation started with a re-equilibration for 20 000 simulation steps, followed by another 30 000 steps to collect configurations every 5th step.

The 6 ps interval of our QM/MM simulation can be considered long enough for structural analysis since it has been shown that even a 2 ps simulation can deliver most of the structural and dynamical properties of similar solutions.28–30

The switching function of eqn. (2) was applied within an interval of 0.2 Å (i.e., in eqn. (3), r1 = 4.0 and r0 = 4.2 Å).

Results and discussion

Fig. 1 displays and compares the over-all K–(N + O) and K–H (from both water and ammonia molecules) RDFs and their corresponding integration numbers obtained from both classical and combined QM/MM simulations.

The classical simulation gives a pronounced first K–(N + O) peak at 2.93 Å, corresponding to an average coordination number of 8.7.

The first solvation shell is not clearly separated from the bulk, indicating an interchange of ligand molecules between the first solvation shell of K+ and bulk.

In the QM/MM simulation, a less pronounced unsymmetrical first K–(N + O) peak, with a maximum at shorter distance of 2.79 Å, is observed.

The shape of the K–(N + O) peak indicates an asymmetry of the K+ solvation.

The first solvation shell of K+ is also not distinctly separated from the bulk, showing even more frequent ligand interchange at the boundary between the solvation shell and bulk.

An integration, calculated up to the first K–(N + O) minimum, yields an average coordination number of 7.6.

In both pair potential and QM/MM simulations, a second solvation shell of K+ is not found, as can be seen from further curves of K–(N + O) RDFs.

The ligand compositions of the solvation shell of K+ can be analyzed by plotting the K+–H2O and K+–NH3 RDFs and their corresponding integration numbers separately, as shown in Figs. 2 and 3, respectively.

For the K+–H2O RDFs (Fig. 2), the QM/MM simulation reveals a broader unsymmetrical first K–O peak with maximum at 2.81 Å, compared to the corresponding peak at 2.75 Å of the pair potential simulation.

The unsymmetrical K–O peak observed in the QM/MM simulation corresponds to weak K+–water interactions and a high mobility of water molecules in the solvation shell of K+.

In addition to the broader K–O peak, a larger and broader K–H(W) peak obtained from the QM/MM simulation (Fig. 2b) supports well a high flexibility and rather freely-oriented arrangements of the water molecules.

The average numbers of water molecules in the solvation shell of K+ are about 5.3 ± 0.2 and 6.7 ± 0.1 for pair potential and QM/MM simulation, respectively.

For the K+–NH3 RDFs (Fig. 3), the pair potential simulation produces a sharp first K–N peak with maximum at 2.96 Å, giving an average coordination number of 3.4 ± 0.2.

In the QM/MM simulation, a less pronounced unsymmetrical first K–N peak is exhibited at a much shorter K–N distance of 2.77 Å, corresponding to approximately 0.9 ammonia molecules.

As in the case of K–O peak, the asymmetrical K–N peak in the QM/MM simulation also indicates a high ammonia mobility, associated with an over-all distorted structure of the solvation shell.

Fig. 4 shows the probability distributions of the solvation numbers of K+, calculated up to the first minimum of the K–(N + O) RDFs.

In the pair potential simulation, the first solvation shell of K+ prefers a coordination number of 8 (followed by 9 and 7 in decreasing amounts), consisting of 5 (and 6 in decreasing amount) water and 3 (followed by 4 in smaller amount) ammonia molecules.

In the QM/MM simulation, the first solvation shell favors a coordination number of 7 (in addition to 8 and 6 in smaller amounts), consisting of 6 (followed by 7, 8 and 5 in smaller amounts) water and one ammonia molecule.

As can be seen from Fig. 4, the 8-, 9- and 7-coordinated complexes are mostly dominant in the pair potential simulation, while numerous possible species, varying from 5- to 10-fold coordinated complexes, exist in the QM/MM simulation.

In comparison to the statistical average distribution of ligands in the solution of approximately 4.4∶1 for water and ammonia molecules, respectively, the QM/MM simulation obviously indicates a preference for H2O with a water-to-ammonia ratio of 7.4∶1.

In contrast to the QM/MM results, the corresponding ratio of 1.6∶1 is observed in the pair potential simulation, which reflects a preference for NH3.

In our previous works, we have shown that the assumption of pairwise additive approximations for the K+(H2O)n30 and K+(NH3)n56 complexes had led to energetic errors of approximately 17% and 13%, respectively, for n = 6, and these errors increased rapidly to 32% and 20% for n = 8.

Undoubtedly, the many-body contributions can play a significant role for the preferential solvation of the K+, and thus the results predicted by classical pairwise additivity are by no means sufficient for a correct description of this ion's solvation in such binary solvent systems.

Fig. 5 displays the O–K–O, N–K–N and O–K–N angular distributions, calculated up to the first minimum of the K–(N + O) RDFs.

In both classical and combined QM/MM simulations, the solvation shell structure of K+ is predicted to be approximately polyhedral.

The rather broad angular distributions indicate the extreme lability and irregular structure of the K+ solvation.

In order to describe the geometrical arrangement of the K+ solvation, the angle θ, defined by the K⋯O and K⋯N axes and the dipole vectors of ligands, has been used to investigate the preferential orientations of solvated ligands.

The dipole-oriented arrangements of water and ammonia molecules within the solvation shell of the ion are given in Figs. 6.

Both classical and combined QM/MM simulations point out that ammonia ligands stick more rigidly to the dipole-oriented arrangement than water molecules.

In the QM/MM simulation, a more flexible arrangement is observed, in particular for water molecules.

As can be seen from Fig. 6a, the QM/MM simulation clearly indicates a higher flexibility of the solvated water's orientation, proving the weak structure-forming ability of K+ even for its nearest water molecules.

The higher mobility of the solvated ligands, as observed in the QM/MM simulation, would increase the number and rate of ligand exchange processes.

Several ligand exchange processes can be seen when the K–O and K–N distances were separately plotted during the QM/MM simulation, as shown in Figs. 7 and 8, respectively.

Within the simulation time of 6 ps, 13 water and 3 ammonia molecules were involved in the exchange processes, indicating the extremely fast dynamics of the K+ solvation.

For the exchange of water molecules, both associative and dissociative interchange (defined as Ia and Id) mechanisms are observed.

In the case of ammonia molecules, the first exchange, which took place around 2 ps, shows the Id mechanism, while the second exchange, which occurred beyond 5 ps, seemingly represents the Ia type.

In order to estimate the rate of ligand exchange processes at K+, the mean residence times (MRT) of the ligands were evaluated within the restrictions set by the available period of the QM/MM simulation.

Based on a direct accounting66 and setting the minimum time a ligand has to remain inside/outside a shell upon a migration process as t*, the MRT values of 2.1 and 1.9 ps for t* = 0.5 and of 15.3 and 5.9 ps for t* = 2.0 are obtained for water and ammonia molecules, respectively.

The corresponding values for water ligands of hydrated Cs+ and for pure water66,67 are 1.7 and 1.8 ps for t* = 0.5, and 10.6 and 13.7 ps for t* = 2.0.

These data show the particular lability of NH3 ligands in the first shell of K+, but also that the water ligands in this shell have a similar mobility as in pure water.

Overall, this represents a well-recognizable structure-breaking effect, although not as pronounced as in the case of Cs+.

The great lability of NH3 ligands in the neighborhood of K+ gives some valuable hints at the binding properties of this ion in biological systems involving N-binding sites.

In a previous paper,30 it has been shown that the main reason of K-specific ion channels to transport this ion preferentially to Na+ is the easy removal of the ion's hydration shell.

The results presented here indicate in addition that O-binding sites should be considerably more competitive in K+ binding than N sites.


The combined ab initio QM/MM molecular dynamics simulation has provided more detailed information on the solvation shell structure of K+ in 18.4% aqueous ammonia solution.

The QM/MM simulation has revealed a flexible polyhedral structure with an average coordination number of 7.6.

The solvation shell contains approximately 6.7 water and 0.9 ammonia molecules, which reflects a clear preference for water molecules.

In addition, the simulation has supplied information on the stability of the K+ solvate, in particular the “structure-breaking” behavior of the ion and its preference for O-binding sites.

The classical simulation based on pairwise additive approximations is apparently not adequate to reproduce well enough the flexible structure and a realistic solvation shell and, therefore, the relevance of n-body and polarization effects, even for weakly interacting ions as K+, could be demonstrated once more.