Computational study of the cerium(III) ion in aqueous environment.

This work comprises the first quantum chemical simulation study of the Ce3+ ion in aqueous environment. The structural and dynamical properties have been investigated by means of the quantum mechanical charge field (QMCF) molecular dynamics (MD) approach and the results, where applicable, have been compared to experimental data. Besides conventional analytical tools, angular radial distribution functions have been employed to gain deeper insight into the structure of the hydrate. The ion-oxygen stretching motion's wavenumber, further characterising the Ce-O bond, is in excellent agreement with experimental results, same as the structural values obtained from the simulation.

Introduction

The highly coordinating lanthanoid ions [1] are of particular interest nowadays with many topical fields of application. Especially the lighter ions, such as Ce3+, exert diverse biological effects mainly by their resemblance to calcium ions, thus enabling for in vivo substitution in biomolecules [2]. Long before the discoverage of Ce3+’s potential biological effects, the antiseptic characteristic of cerium ions was known [3,4], leading to the investigation of virus removal qualities of biogenic cerium [5]. Most recently, Ce3+’s anti-cancer potential [6], the character of cerium based complexes [7] and the influence of cerium salts on fibroblast and osteoblast proliferation [8] have been investigated. These fields of interest, all requiring in detail understanding of structural and dynamical features of the ion in solution, were a main reason for conducting a quantum mechanical simulation study of the hydration of Ce3+.

To date, characterisation of hydrated trivalent cerium has been carried out by various EXAFS experiments [9-11] and via mid-infrared/Raman spectroscopy [12]. Due to the nature of such experiments, however, ultrafast dynamical inter- and intramolecular processes and an in-detail structural description are not available. As an observation of such processes is still not feasible with modern experimental techniques, the application of theoretical simulations for this task appeared promising. To date, only one in vacuo geometry study of cerium–water clusters at the Becke3-LYP level [13] has been performed. Thus, it looked to be a fruitful task to conduct a simulation study enabling structure and dynamics based conclusions for a system simulated at ambient conditions.

Methods

Simulation method and setup

Similar as the conventional QM/MM-MD approach [14,15], the QMCF-MD ansatz [16-19] is also based on partitioning the simulation box into quantum mechanically and molecular mechanically treated regions. In this framework, the QM region is split into two sub-regions, resembling the chemically relevant first and second hydration spheres. The construction of non-Coulombic solute–solvent potentials, being a difficult and tedious task, can be renounced in this approach, as the radius of the quantum mechanically treated region is on the scale as the non-Coulombic cutoff distance, which is possible due to the high performance computational facilities available to date. Detailed information on the methodology of the QMCF-MD ansatz is given in the corresponding works [16-19], describing, besides the employment of calculation formalisms, the setup of interaction forces, the realisation of smooth particle migration between QM and MM regions, also the versatile applicability of this approach.

Prior to the simulation, geometry optimisations of cerium–water clusters similar as in the previously mentioned work [13] have been conducted with the Gaussian09 software package [20] in order to confirm the quality of the chosen effective core potential (ECP) basis-set for cerium [21]. The quantum mechanical part of the simulation was based on the Hartree–Fock level and the quality of this formalism was also proven along the mentioned cluster calculations, yielding only minimal electron correlation contributions when compared to optimisations based on multideterminantal calculation methods. The resulting equilibrium geometry (tri-capped trigonal prismatic hydration structure) of the [Ce(H2O)9]3+ cluster is in agreement with the reported experimental [10,11] and theoretical results [13], further underlining the quality of the chosen method and the basis set. The popular and well established Dunning double-ζ basis set augmented by polarising functions was employed for the description of hydrogen and oxygen in the QM region [22-24]. For a proper description of the solvent in the MM region, the flexible BJH-CF2 water model [25,26], was used.

Simulation protocol

The cubic simulation box, containing 1000 water molecules and one Ce3+ ion, was defined with a side length of 31.15 Å and a density of 0.997 g/cm3. The radii of the QM core and layer zones were set to 3.2 and 6.2 Å, respectively, already including a 0.2 Å thick smoothing region. The constant temperature of 298 K was maintained by the Berendsen weak-coupling algorithm with a 0.1 ps relaxation time [27] and the canonical NVT ensemble was chosen for the simulation. The long-range Coulombic cut-off (15 Å) was supplemented by the reaction field approach (ε = 78.36) [28], and the integration of the equations of motion was realised with an Adams–Bashforth predictor–corrector algorithm. A bulk-like environment, avoiding surface conditions, was established by utilising periodic boundary conditions and the minimum image convention. The initial equilibration period took 2 ps, followed by heating the hydrate to 600 K and a subsequent re-equilibration at room temperature for another 2 ps. The actual sampling was conducted for 26 ps at a simulation time step length of 0.2 fs. The rather high loop time of roughly 250 s per simulation step on the utilised 16 CPU core platform can mainly be attributed to the required open-shell treatment of the solvate, resulting in a total of 15 months of calculation time. The quantum mechanical part of the simulation and the cluster optimisations were both treated with the Gaussian09 software package [20].

Results and discussion

The hydrate proved stable throughout the simulation trajectory and even during the heating period, no hydrolysis reactions have been observed. Very few ligand displacement reactions occurred between the first and the second hydration sphere, yielding a mean first shell coordination of ∼9.1 water molecules around the solute.

The RDF in Figure 1 indicates well defined first and second hydration spheres with their maxima being located at 2.61 Å and 4.95 Å, respectively. A weak third hydration sphere has been identified as well and its existence is proven by the plots [29-31] in Figure 2. Table 1 compares the obtained data with experimental results and the values from the theoretical approach by Dinescu and Clark [13], attesting the quality of this QMCF-MD based treatment of the hydrate. The first hydration sphere is mainly (90% relative occurence) characterised by a ninefold coordination and the rarely observed tenfold coordination is attributed to the few ligand exchange reaction between the innermost two hydration spheres. The mean first shell coordination number of ∼9.1 water molecules surrounding the ion is in good agreement with most of the experiments [10,11], as shown in Table 2. The twelvefold coordination from the 1995 EXAFS experiment reported by Solera et al. [9] has been corrected by the works of Persson et al. [10], Allen et al. [11] and the comprising works of Helm and Merbach [32-34]. It is generally accepted today, that lighter lanthanoid ions are mostly ninefold coordinated while with higher element number, the first shell coordination number is lowered down to ∼8 in the case of lutetium [10]. For the much more dynamic second hydration sphere, many ligand exchange reactions were observed. The corresponding shell boundaries (Table 2), indicating a broader layer of hydration, prove the weaker influence of the cation at a longer distance.

Figure 1

Ce3+—O and Ce3+—H RDFs and their running integrations.

Figure 2

Local density corrected three-body distribution functions for the first, second and third shell of hydration. Overlay of the O–O pair distribution function for pure water (dashed line) is given for comparison [41].

Table 1

Hydration sphere radii () and maximum shell peaks () of the Ce3+–hydrate in Å – all values refer to oxygen.

Study	r_min,1	r_max,1	r_min,2	r_max,2	r¯_1	r¯_2
QMCF-MD (this work)	2.3	3.4	3.9	5.8	2.61	4.95
EXAFS [9]	n.a.	n.a.	n.a.	n.a.	2.55	n.a.
EXAFS [10]	n.a.	n.a.	n.a.	n.a.	2.538	n.a.
EXAFS [11]	n.a.	n.a.	n.a.	n.a.	2.52	n.a.
B3LYP-Opt [13]	n.a.	n.a.	n.a.	n.a.	2.60	n.a.

Table 2

Minimum () and maximum () coordination numbers and mean coordination numbers () of the Ce3+–hydrate.

Study	CN_min,1	CN_max,1	CN_1	CN_2
QMCF-MD (this work)	9	10	9.1	21.4
EXAFS [9]	n.a.	n.a.	12	n.a.
EXAFS [10]	n.a.	n.a.	9	n.a.
EXAFS [11]	n.a.	n.a.	9.3	n.a.
B3LYP-Opt [13]	8	9	8–9	n.a.

The O—Ce—O angular distribution function (ADF) in Figure 3 shows the well defined maxima of the first shell hydrate structure in the regions of and . The narrowness of the peaks reflects a strong solute–solvent interaction and the values of the peak maxima resemble the two main structural motifs observed: tri-capped trigonal prism and capped square antiprism (Figure 4a). The valley between the two peaks, ranging from roughly to , indicates numerous structural interconversions between the two structures. These interconversions result from either minor ultrafast geometric modifications within the hydrate or from pseudorotations, requiring substantially longer times in the scale of 1 ps for establishment of an altered geometry. The ARD plot [35] in Figure 4a, which was drawn over the full trajectory length, shows, that the hydrate is indeed characterised by frequent structural interconversions. The capped square antiprismatic structure, however, is slightly more dominant, as it can be noticed from the higher localisation densities in the respective spots (yellow colour for the prismatic positions and red colour for the capping position).1 The lighter densities within the white circles reflect the slightly less frequently established tri-capped trigonal prismatic structure, resulting from ligand reorientation. It has to be mentioned, however, that the respective densities in the ARD plot have to be treated with care, as both reported structures share almost identical properties in terms of bond lengths and angles, ultimately leading to a somewhat diffuse ligand localisation. This can be explained by the dynamic character of the liquid state and the frequent formation and break-up of hydrogen bonds. Computation of [Ce(H2O)9]3+ clusters in tri-capped trigonal prismatic and capped square antiprismatic configurations at Hartree–Fock level yielded an energy difference of 3.05 kcal/mol (see supporting information). At first sight this appears large compared to the energy stored within one degree of freedom being 0.296 kcal/mol at room temperature (1/2 kBT). However, the first shell water molecules plus the ion store an equivalent of 24.9 kcal/mol (84 degrees of freedom) which is more than sufficient to overcome the barrier to transform the structure from a tri-capped trigonal prismatic to a capped square antiprismatic coordination, whenever concerted motion of some of the ligands and/or the ion occurs.

Figure 3

O—Ce3+—O ADF evaluated for the first hydration sphere.

Figure 4

Ce3+—O ARD plot drawn over the full trajectory length.

To facilitate understanding of an ARD plot only depicting a diffuse mean structure if drawn over a whole simulation trajectory, Figure 4b shows how similar the two observed hydrate structures actually are: following the four screenshots on the right, each depicting the same configuration, the capped square antiprism (visualised via the red bonds) is tilted vertically. The fourth hydrate structure now resembles a tri-capped trigonal prism where the corners of the prism are formed by the blue oxygens. For improved visualisation, the bonds have been altered in this picture. In order to obtain a well defined trigonal prism, only very minor movement of two oxygens (see arrows in Figure 4b) is required, demonstrating the subtle difference between the tri-capped trigonal prismatic and the capped square antiprismatic configuration.

For the second hydration sphere, the ARD plot confirms the previously mentioned wider shell boundaries and the much more dynamic character within this layer and the third sphere of hydration.

Table 3 lists the computed mean ligand residence times and compares the values to the ones obtained for other highly charged cations. With only two ligand displacements observed between the first and the second layer of hydration in the 26 ps long trajectory, the mean first shell ligand residence time was computed as >118 ps. While the statistics based on such scarce exchange events may not be accepted as being accurate, they still do indicate a tendency of the hydrate having a faintly dynamic character. This hypothesis can be further underlined by examining the ligand exchange statistics between the second and the third hydration sphere (Table 3, values and ). denotes the mean ligand residence time computed with a ligand displacement lasting at least 0.5 ps and indicates the number of required exchange attempts until a lasting ligand displacement is achieved.

Table 3

First shell () and second shell () MRTs and corresponding values obtained for various polarising ions in aqueous solution.

Ion	τ10.5 (ps)	τ20.5 (ps)	R_EX,1	R_EX,2
CeQMCF3+ (this work)	>118	2.6	3.0	5.5
AlQMCF3+[36]	n.a.	17.7	n.a.	21.1
UQMCF4+[37]	n.a.	8.1	n.a.	5.9
ZrQMCF4+[38]	n.a.	5.5	n.a.	6.8
LaQMCF3+[39]	16.6	2.3	2.4	2.3
HfQMCF4+[40]	n.a.	15.4	n.a.	10.6
H2O_QM/MM[41]	1.7	n.a.	11.2	n.a.

From the well-equilibrated trajectory, the cerium–oxygen bond was characterised by calculating the corresponding stretching motion’s wavenumber and the bond force constant (Table 4). As a Raman experiment has been conducted in the past [12], the theoretically derived results can be directly compared to the experimental data, yielding excellent agreement with the literature value ( vs. ). Irrespective of differences in the electronic structure of the compared ions, the bond force constant serves as a useful tool for measuring and comparing the ion–oxygen bond strength on an objective basis. Furthermore, as the frequency depends on the second derivative of the energy with respect to the nuclear coordinates, it can be seen as a very sensitive probe of the accuracy of the simulation. The corresponding Ce—O force constant is 106 compared to 109 derived from the Raman data.

Table 4

Comparison of theoretically derived stretching frequencies Q and force constants k (top) with experimental results (bottom), ∗ value has been unscaled prior to comparison [44,45].

Ion	Q_IonO(cm^-1)	k_IonO(Nm^-1)
AlQMCF3+∗[36]	560	194
ZrQMCF4+[38]	484	188
LaQMCF3+[39]	354	106
HfQMCF4+[40]	n.a.	212
FeQMCF2+[42]	357	93
FeQMCF3+[42]	513	193
CeQMCF3+ (this work)	354	106
CeRaman3+[43]	359	109

Conclusion

This work presents the first ab initio simulation treatment of trivalent cerium in aqueous environment. The significant sampling period and the substantially large QM treated region both contribute to the quality of the reported structural and dynamical data. The quality of the employed QMCF-MD approach was proven with this study not only because of the fact that quantum mechanical simulations of highly charged ions are a rather delicate task, but also because of the excellent agreement of the obtained results with experimental data, in particular the vibrational mode of the ion–oxygen stretching motion. Analysis of the angular radial distribution of ligands yielded deeper insight into the actual hydration structure and its dynamical character. Besides two well defined layers of hydration, a weak third hydration sphere has been identified not only from reviewing the pair distribution function of ligands, but also from investigating the local density corrected three-body functions.