Hydration Shell Structure and Dynamics of Curium(III) in Aqueous

ARTICLE pubs.acs.org/JPCA

Hydration Shell Structure and Dynamics of Curium(III) in Aqueous Solution: First Principles and Empirical Studies Raymond Atta-Fynn,*,† Eric J. Bylaska,† Gregory K. Schenter,‡ and Wibe A. de Jong† †

W. R. Wiley Environmental Molecular Sciences Laboratory and ‡Chemical and Material Sciences Division, Pacific Northwest National Laboratory, P.O. Box 999, Richland, Washington 99352, United States

bS Supporting Information ABSTRACT: Results of ab initio molecular dynamics (AIMD), quantum mechanics/molecular mechanics (QM/MM), and classical molecular dynamics (CMD) simulations of Cm3þ in liquid water at a temperature of 300 K are reported. The AIMD simulation was based on the CarParrinello MD scheme and GGA-PBE formulation of density functional theory. Two QM/ MM simulations were performed by treating Cm3þ and the water molecules in the first shell quantum mechanically using the PBE (QM/MM-PBE) and the hybrid PBE0 density functionals (QM/MM-PBE0). Two CMD simulations were carried out using ab initio derived pair plus three-body potentials (CMD-3B) and empirical Lennard-Jones pair potential (CMD-LJ). The AIMD and QM/MM-PBE simulations predict average first shell hydration numbers of 8, both of which disagree with recent experimental EXAFS and TRLFS value of 9. On the other hand, the average first shell hydration numbers obtained in the QM/MM-PBE0 and CMD simulations was 9, which agrees with experiment. All the simulations predicted an average first shell and second shell CmO bond distance of 2.492.53 Å and 4.674.75 Å respectively, both of which are in fair agreement with corresponding experimental values of 2.452.48 and 4.65 Å. The geometric arrangement of the 8-fold and 9-fold coordinated first shell structures corresponded to the square antiprism and tricapped trigonal prisms, respectively. The second shell hydration number for AIMD QM/MM-PBE, QM/MM-PBE0, CMD-3B, and CMD-LJ, were 15.8, 17.2, 17.7, 17.4, and 16.4 respectively, which indicates second hydration shell overcoordination compared to a recent EXAFS experimental value of 13. Save the EXAFS spectra CMD-LJ simulation, all the computed EXAFS spectra agree fairly well with experiment and a clear distinction could not be made between configurations with 8-fold and 9-fold coordinated first shells. The mechanisms responsible for the first shell associative and dissociative ligand exchange in the classical simulations have been analyzed. The first shell mean residence time was predicted to be on the nanosecond time scale. The computed diffusion constants of Cm3þ and water are in good agreement with experimental data.

1. INTRODUCTION Curium (Cm) is a strong R-particle emitter found in spent nuclear fuel and high-level radioactive waste.1 As such, it is responsible for many of the long-term radiological and thermal problems associated with nuclear waste storage and disposal.1 Like many actinides, Cm exhibits a rich chemistry leading to additional complications in containment and remediation. Its ground state electronic configuration is [Rn]5f76d17s2 and its most stable oxidation state is Cm(III) (i.e., [Rn]5f7), which corresponds to the 8S7/2 state.2 Unfortunately, the chemical behavior and long-term fate of Cm(III) in solution is not fully understood. A thorough understanding of the physical and chemical behavior of Cm(III) and its complexes in solution is crucial to the assessment of the economical and environmental aspects of nuclear waste treatment. In this work, we provide a detailed theoretical description of the hydration shell structure and ligand exchange dynamics of the Cm(III) ion in aqueous solution using ab initio methods and parameter-based force fields. r 2011 American Chemical Society

Experiments have investigated the speciation of Cm3þ in various chemical environments (see Edelstein et al.2 for a comprehensive review), yet there is some debate over the hydration shell structure of Cm3þ. Laser-induced fluorescence (LIF) experiments by Kimura and co-workers3 using different aqueous solutions yielded primary hydration numbers of 8 or 9. From EXAFS measurements, Allen et al.4 obtained a coordination number of 10 and a corresponding CmO distance of 2.45 Å using 0.25 M HCl solution. Timeresolved laser fluorescence spectroscopy (TRLFS) by Stumpf et al.5 found a Cm3þ coordination number of 9 in glycolic acid. Recent TRLFS experiments on Cm3þ in 0.1 M perchloric acid by Lindqvist-Reis et al.6 clearly indicated that near room temperature (293 K), the first hydration shell coordination was 90% 9-fold and 10% 8-fold. They further observed that at 473 K, the 9-fold Received: January 31, 2011 Revised: March 30, 2011 Published: April 18, 2011 4665

dx.doi.org/10.1021/jp201043f | J. Phys. Chem. A 2011, 115, 4665–4677

The Journal of Physical Chemistry A coordination reduces to 60% and the 8-fold coordination increases to 40%. Finally, recent EXAFS experiments for 0.523 M Cm3þ solution prepared in 1 M perchloric acid by Skanthakumar and coworkers7 showed that the primary hydration shell comprised nine water molecules arranged in a trigonal-tricapped prism (TTP) around Cm3þ. In the TTP structure, six water molecules are located at the vertices of the trigonal prism with a CmO bond length of 2.45 Å, and three water molecules are located at the capping faces and lying in the equatorial plane with a CmO bond length of 2.55 Å. Thus, experimental data on Cm3þ in the aqueous environment suggests that the primary shell hydration number of Cm3þ is either 8, 9, or 10, with the a hydration number of 9 occurring frequently. The only experimental information on the secondary coordination shell of Cm3þ is due to Skanthakumar et al.7 They measured the second shell coordination number to be 13 ( 4 and a corresponding CmO distance of 4.65 Å using high-energy X-ray scattering (HEXS). The actinides and their fellow f-electron brethren, the lanthanides, show a similar trend in the size of the ionic radii in their þ3 oxidation states, where their ionic radii decrease smoothly across each series as the 5f (4f) shell is gradually filled.8 It is well-known in experimental literature that the first shell hydration number of the trivalent lanthanides in aqueous solution is 9 at the beginning of the series (La3þ) and eventually becomes 8 (Lu3þ) at the end of the series. Specifically, the coordination number changes discontinuously from 9 to 8 at Gd3þ, the so-called “gadolinium break”.9 Generally, the first shell hydration number of an ion solution is correlated with the ionic radius, which makes this property a good predictor of solvation structure. Although the Gd3þ valence is isoelectronic to the Cm3þ valence, the ionic radius of Cm3þ (0.97 Å) is closer to that of Sm3þ (0.964 Å) or Eu3þ (0.95 Å)—both lanthanides have first shell hydration numbers close to 9—than to that of Gd3þ (0.938 Å). Thus, one would expect the hydration number of Cm3þ to be 9 or close to 9. Only a few theoretical calculations on hydrated Cm3þ have been reported in the literature.1013 Mochizuki and Tatewaki10 performed gas-phase all-electron DiracHatreeFock cluster calculations on [Cm(OH2)n]3þ (n = 1, 2, 4, 6). They observed that with increasing n, the CmO bond distance increases and the binding energy per water molecule decreases. Wiebke et al.11 used gradient-corrected density function theory and the COSMO implicit solvent model to calculate the geometry, hydration energies, and binding energies of [Cm(OH2)n]3þ (n = 79) and [Cm(OH2)n-1 3 OH2]3þ (n = 8, 9) clusters. While the CmO bond distances showed an overall agreement with experiment, they predicted a primary hydration number of 8 using the DFT-BP86 exchange-correlation functional, which is slightly undercoordinated compared to recent experiments.6,7 Their prediction was based on the difference ΔEn between the hydration energies of [Cm(OH2)n]3þ(aq) and [Cm(OH2)n-1 3 OH2]3þ(aq). Specifically, ΔE8 = 3.45 kcal/mol and ΔE9 = þ3 kcal/mol, implying that [Cm(OH2)8]3þ(aq) is more stable. In fact, Wiebke et al.11 predicted the first shell coordination for all the þ3 actinide ions, Ac3þLr3þ, to be 8 at the DFT-BP86 level of theory. On the other hand, their second-order MøllerPlesset perturbation (MP2) calculations predicted preferred first shell coordination numbers of 9 and 8 for Ac3þMd3þ and No3þLr3þ respectively and, hence, a break in the coordination number. Strictly speaking, however, the break in the coordination number is no longer analogous to the Gd break in lanthanides, since the ionic radius of Gd3þ (0.938 Å) corresponds to that of

ARTICLE

Cf3þ (0.934 Å).8 In other words, on the basis of the ionic radii, the Gd break in lanthanides is analogous to the Cf3þ break in the actinides. Nevertheless, the results demonstrate the sensitivity of the first shell coordination number to the applied level of theory. Yang and Bursten12 studied the hydration number and geometry of Cm3þ using hybrid B3LYP density functional theory, the conductor-like continuum polarizable model, and [Cm(OH2)n]3þ (n = 810) and [Cm(OH2)n 3 (OH2)m]3þ clusters. They predicted the primary hydration number of 9 with a TTP geometry, which is in agreement with TRLFS6 and EXAFS7 experiments. In the same work, they performed a rigid hydrated ion model classical molecular dynamics (CMD) simulation at 298 K. In the rigid hydrated ion model, the entire first shell geometry, that is, [Cm(OH2)9]3þ, was held rigid and treated as a single solute interacting with water molecules in the second and bulk solvation shells. The rigid hydrated ion model is unphysical, since crucial first shell geometrical parameters such the mean CmO bond distances and OCmO bond angles and dynamical parameters such ligand exchange rates and the selfdiffusion constant of Cm3þ cannot be measured. However, their simulation predicted the average second shell CmO distance to be 4.65 Å, which matches with experiment,7 but the predicted second shell hydration number of 21 deviated significantly from the experimental value of 13.7 Hagberg et al.13 performed CMD simulations at 300 and 473 K using potentials constructed from ab initio data. They reported average first shell and second shell CmO distances of 2.55 and 4.9 Å, respectively, which are both longer than the EXAFS value of 2.48 and 4.65 Å. Using a second shell radial cutoff of 3.0 Å (3.4 Å), their simulation predicted the first hydration shell to be 95% 9-fold and 5% 8-fold (90% 9-fold and 10% 8-fold) at 300 K and 80% 9-fold and 20% 8-fold (60% 9-fold and 50% 8-fold) at 473 K. The results at both temperatures are in good agreement with the aforementioned TRLFS experiments.6 Furthermore, they reported the most probable second shell hydration to be 16. While the theoretical gas phase and molecular dynamics simulations studies summarized above give some information about Cm3þ hydration, the solvent shell geometry and dynamics as well as the electronic structure require a more detailed investigation. In this work, we employ parameter-free density functional theory based CarParrinello ab initio molecular dynamics (AIMD) simulations, partially parametrized quantum mechanics/molecular mechanics (QM/MM) MD simulations, and fully parametrized classical molecular dynamics (CMD) simulations to examine the solvent shell structures and dynamical processes of Cm3þ immersed in bulk water. Specifically, we investigated the first shell geometry, hydrogen-bonding network, electronic structure of the solvating waters, mean ligand residence times, and exchange mechanisms. The remainder of this paper is organized as follows. In section 2, we discuss the various computational methods and simulation details. In section 3 we present and discuss the results. The concluding remarks are presented in section 4.

2. COMPUTATIONAL METHODOLOGY Five liquid-phase simulations were performed: one AIMD, two QM/MM MD, and two CMD simulations. A cubic box containing one Cm3þ ion and 64 water molecules was employed for the AIMD and QM/MM MD simulations, while a box containing one Cm3þ ion and 512 water molecules was employed for the CMD simulations. All simulations employed periodic boundary conditions to mimic bulk conditions, and the density of water in each box was approximately 1 g/cm3. The 4666

dx.doi.org/10.1021/jp201043f |J. Phys. Chem. A 2011, 115, 4665–4677

The Journal of Physical Chemistry A

ARTICLE

Table 1. Summary of the MD Simulation Parameters simulationa

system size

Ld

Δte

equilibration time

collection time

AIMD

Cm3þ þ 64H2O

12.42

5 au

2.5 ps

18.2 ps

QM/MM-PBE

Cm3þ þ 9H2O/55H2O b

12.42

5 au

8.5 ps

45.1 ps

QM/MM-PBE0

Cm3þ þ 9H2O/55H2O c

12.42

5 au

9.0 ps

43.9 ps

CMD-3B

Cm3þ þ 512H2O

24.83

0.5 fs

1 ns

4 ns

CMD-LJ

Cm3þ þ 512H2O

24.83

0.5 fs

1 ns

4 ns

All simulations were carried out in the canonical ensemble at a temperature of 300 K. b Cm3þ and nine first shell waters were treated at the QM level using the PBE functional. c Cm3þ and nine first shell waters were treated at the QM level using the PBE0 hybrid functional. d L is the box length in Å which was chosen to yield to water density of approximately1 g/cm3. e Δt is the simulation time step. a

þ3 charge on the metal cation was neutralized by a uniform background charge of opposite sign as opposed to explicitly using counteranions. The initial configurations were generated by placing the Cm3þ at the center of the box and successively placing water molecules in the periodic box at random positions subject to the constraint that RCmO g 2.5 Å, ROO g 2.7 Å, and the dipole vectors of the water molecules inside or close to the first shell should point away from Cm3þ. With such an initial structure, equilibration is achieved very quickly because the structure is geometrically robust. The FORTRAN code used to generate the initial configurations is provided in the Supporting Information. Well-equilibrated configurations from the CMD simulations were used as the initial coordinates for AIMD and QM/MM simulations. The AIMD simulation was carried out within the DFT pseudopotential plane-waves formalism using the CarParrinello algorithm14 as implemented in the NWChem computational package.15 The PBE16 formulation for the GGA exchange-correlation functional was employed. The kinetic energy and charge density cutoffs were 120 and 240 Ry, respectively. The valence electron interactions with the atomic O and H cores were approximated using the generalized norm-conserving Hamann pseudopotentials17,18 modified into a separable form, as suggested by Kleinman and Bylander.19 The s and p cutoff radii of H were rcs = 0.8 au and rcp = 0.8 au, respectively, and the s, p, d cutoff radii of O were rcs = 0.7 au, rcp = 0.7 au, rcd = 0.7 au. A TroullierMartins pseudopotential20 was generated for Cm3þ using [Rn]5f76d0 as the reference state. 78 electrons were treated in the core region and 6s26p65f76d0 was treated as the valence region. Cutoff radii of 1.3, 1.55, 2.35, and 1.2 au were used for the s, p, d, and f angular momentum channels, and the local part of the pseudopotential was chosen to be the s state. The CarParrinello equations of motion were integrated in the presence of NoseHoover thermostats,21,22 coupled separately to the electronic and ionic degrees of freedom. A fictitious electronic mass of 600 au was used. To facilitate the integration of the equations of motion, all the hydrogen atoms were replaced by deuterium. A summary of the AIMD simulation parameters can be found in Table 1. In the QM/MM MD simulations, the QM region comprised [Cm(OH2)9]3þ (first hydration shell) while the MM region comprised the remaining 55 water molecules in the second hydration shell and the bulk region. The QM region was described using DFT with pseudopotentials and plane wave basis. The pseudopotential parameters and plane wave energy cutoffs were the same as in the AIMD calculations. The PBE16 (QM/MM-PBE) and hybrid PBE0 density functionals23 (QM/MM-PBE0) were employed. The molecular dynamics simulation within the region was carried out using the CarParrinello algorithm. A detailed explanation on the

implementation of the QM/MM methodology in the NWchem code has been described elsewhere.24 The PBE and PBE0 functionals were employed to probe the possible dependence of the first hydration shell properties on density functionals. The MM region was described with classical molecular dynamics using the flexible SPCToukanRahman potential for water25,26 (see the Supporting Information). The nonelectrostatic interaction between the QM and MM regions at the interface was described by a Lennard-Jones interaction, and the corresponding parameters for the atoms in the QM region are given by εCm3þ = 0.013 kcal/mol, εO = 0.1554 kcal/mol, εH = 0.044 kcal/mol, σCm3þ = 3.7042 Å, σO = 3.1655 Å, and σH = 0.70 Å. The summary of the simulation parameters for the two QM/MM MD simulations can be found in Table 1. The QM/ MM simulation scheme (temperature control, integration of equation of motion, etc.) was the same as the AIMD simulation. We must mention that there are some artifacts in the QM/MM method that can degrade the results. For example, the interaction that couples the QM and MM boundaries (the coupling is described by electrostatic and Lennard-Jones potentials) sometimes promotes an unphysical exchange of water molecules between the QM and MM regions. As a result, particles that make transitions from the MM region into the QM region will still be described with an MM potential (because the QM and MM atoms are fixed throughout the simulations), and this is problematic. In this work, any data collected after an MM water molecule makes a permanent transition into the QM region is discarded. Classical molecular dynamics (CMD) simulations were carried out to probe long time scale phenomena, such as particle exchange mechanisms, mean residence time, and translational diffusion. The intermolecular water interaction was modeled by the rigid body SPC/E potential.27 We developed two potentials to describe the Cm(III)water interaction (see Supporting Information for the details). The potentials were employed to elucidate the importance of including three-body terms in the classical description of the system. In one simulation, the Cm(III)water interaction was described using the two- and threebody analytic potentials (CMD-3B) fitted from ab initio restricted open shell HartreeFock (ROHF) calculations.28 The total potential energy for the CMD-3B simulation, UCMD-3B, is expressed as UCMD-3B ¼ Uele þ ULJ ðO OÞ þ U2B ðCm3þ H2 OÞ þ U3B ðO Cm3þ OÞ

ð1Þ

where Uele is the electrostatic interaction energy term composed of the solventsolvent and solventsolute interactions, ULJ(OO) is the Lennard-Jones potential describing the nonelectrostatic short-range interaction between the oxygen 4667



ARTICLE

Table 2. Two and Three-Body Interaction Parameters Used for the Classical MD Simulations Two-Body Interaction Parameters CMD-3B 1

ACmO (kcal mol )

90 139.305 454 3

ACmH (kcal mol1)

19 533.436 139 1

BCmO (Å1)

CMD-LJ

3.411 177 7

BCmH (Å1)

3.239 302 2

CCmO (kcal mol1 Å4)

618.274 311 58

CCmH (kcal mol1 Å4)

948.021 166 78

DCmO (kcal mol1 Å6) DCmH(kcal mol1 Å6)

482.957 869 87 349.304 082 54

εCm3þ (kcal mol1)

0.013

σCm3þ (Å)

3.700 Three-Body Interaction Parameters CMD-3B

R (kcal mol1) β (Å1)

188.2527 0.925 965 73

σ (Å1)

0.056 691 77

Figure 1. CmO radial distribution function (upper panel) and running coordination number (bottom panel) plots per simulation.

centers of SPC/E water molecules, U2B(Cm3þH2O) is the nonelectrostatic short-range solutesolvent two-body interaction, and U3B(OCm3þO) is nonelectrostatic short-range solventsolutesolvent three-body interaction. The explicit analytic forms for the two-body U2B(Cm3þH2O) and threebody U3B(OCm3þO) solutesolvent interactions are given by X R¼O, H1, H2

U2B ðCm3þ H2 OÞ ¼ ACmR expð BCmR rCmR Þ þ

CCmR DCmR þ rCmR 4 rCmR 6

U3B ðO Cm3þ OÞ ¼ R expð βr1 βr2 γr 3 Þ

ð2Þ ð3Þ

The two-body parameters A, B, C, D, and the three-body parameters R, β, and γ are listed in Table 2. In the second classical simulation, the Cm(III)water interaction was described using a purely empirical pair Lennard-Jones potential (CMD-LJ). The CMD-LJ Lennard-Jones parameters were obtained by adjusting the Rappe et al. unified force field (UFF) parameters29 for Cm3þ so that the first and second peak positions in the simulated CmO radial distribution function approximately match the experimental first and second shell CmO distances, respectively. In this case, the total pair potential energy, UCMD-LJ, is expressed as UCMD-LJ ¼ Uele þ ULJ ðO OÞ þ ULJ ðCm3þ OÞ

ð4Þ

where Uele and ULJ(OO) have the same definition as in eq 1 and ULJ(Cm3þO) is the nonelectrostatic solutesolvent pair Lennard-Jones interaction. The Lennard-Jones parameters for Cm3þare listed in Table 2 and Lennard-Jones parameter for O is the same as in the SPC/E interaction. The effective LJ parameters εCm3þO and σCm3þO for ULJ(Cm3þO) were obtained using the LorentzBerthelot mixing rules: εCm3þO = (εCm3þεO)1/2 and σCm3þO = 1/2(σCm3þ þ σO) The NoseHoover thermostat chain21,30 was used to control the temperature in both CMD simulations. The length of

the thermostat chain was 5. A YoshidaSuziki integrator31,32 of order 9 was used to propagate thermostat positions and velocities. The long-range electrostatic interactions were treated using the classical Ewald summation technique.33,34 The Ewald convergence parameter was set to R = 5/L Å1. The real space Ewald summation was carried out in the central unit cell using the standard minimum image convention. In the reciprocal summation, truncation of the modulus of the reciprocal lattice vector was determined by |η| e 27, where the reciprocal lattice vector is given by k = 2πη/L = (2π/L)(η1,η2,η3), and ηi = 0, (1, (2... . The RATTLE algorithm35 was used to constrain the geometry and velocities of the water molecules. The relevant input parameters for the CMD simulations are listed in Table 1.

3. RESULTS AND DISCUSSIONS 3.1. Hydration Shell Properties. In the radial distribution function (RDF) plot in Figure 1, well-defined isolated peaks can be seen in the 23 and 45.3 Å ranges for the first and second solvation shells, respectively. The first shell parameters in Table 3 lists the average first shell Cm3þO separation, rCm-OI, to be 2.50 Å for the AIMD simulation and in the 2.522.53 Å range for the QM/MM and CMD simulations. The corresponding deviations in the Cm3þO, σCm-OI, are 0.11 and 0.080.11 Å. Compared to experimental data, all the bonds are slightly longer but shorter when compared to the past theoretical data of Hagberg et al.13 The average AIMD and QM/MM Cm3þO bond lengths are also slightly contracted compared to the gasphase values of [Cm(OH2)8]3þ and [Cm(OH2)9]3þ, respectively. Also OH bond lengths and HOH bond angles show very small deviations from the gas-phase values. A key feature of the gCmO(r) RDF in Figure 1 is the character of the peak height and width in the first shell. The AIMD peak height is the shortest, mainly due to the fact that it has the lowest first shell coordination number (the coordination number is proportional to the area under the curve), followed by the QM/MM-PBE curve with a primary coordination number of 8. The QM/MM-PBE0 and CMD-3B curves with coordination number 9 have nearly the same character, but the height difference between CMD-LJ first shell and second shell peaks is the largest, indicating that the first 4668



ARTICLE

Table 3. Average Geometric Parameters of the First and Second Solvation Shell of Cm3þ per Simulation rCm-Oa

σCm-Ob

rCm-H

Average First Shell Parameters rOH σOH — HOHc

σCm-H

σ — HOH

ψtiltd

CNe

[Cm(OH2)8]3þf

2.51

0.95

105.8

[Cm(OH2)9]3þf

2.56

0.95

105.4

AIMD

2.50

0.11

3.12

0.14

0.97

0.03

106.2

5.3

30

8

QM/MM-PBE

2.49

0.10

3.11

0.13

0.97

0.02

105.8

4.7

30

8

QM/MM-PBE0 CMD-3B

2.52 2.53

0.10 0.11

3.12 3.18

0.13 0.16

0.95 1.00

0.02 0.00

106.0 109.5

4.7 0.0

30 23

9 9

CMD-LJ

2.52

0.08

3.18

0.11

1.00

0.00

109.5

0.0

16

ref 13

2.55

8 9

expt3 expt4

9 9 8/9

2.45

10

expt5

9

expt6

8/9 g

expt7

2.48

3.05

9 Average Second Shell Parametersh rOIOII σOIOII — OIHIOII

rCm-OII

σCm-OII

rCm-HII

σCm-HII

AIMD

4.71

0.31

5.03

0.47

2.85

0.18

162

8.9

15.8

QM/MM-PBE

4.64

0.30

5.05

0.47

2.81

0.17

163

8.9

17.7

QM/MM-PBE0 CMD-3B

4.68 4.70

0.28 0.32

5.10 5.11

0.47 0.47

2.83 2.84

0.17 0.17

162 162

8.9 9.2

17.7 17.4

CMD-LJ

4.70

0.30

5.16

0.40

2.82

0.17

162

9.2

16.4

ref 12

4.65

ref 13

4.90

expt7

4.65

5.25

2.65

σ — OIHIOII

CN

21.0 16.0

2.742.81

13 ( 4

rA-B is the average AB bond distance in Å. b σ denotes the standard deviation. c — ABC denotes the average ABC bond angle (in deg). d ψtilt is the water tilt angle of the first shell water with respect to the Cm3þO bond (in deg) e CN is the hydration shell coordination number based on cutoff distances of 3.3 and 5.3 Å for the first and second coordination shells, respectively. f DFT-PBE plane wave pseudopotential gas-phase calculations. g Eight-fold coordination occurs with 10% probability, and 9-fold coordination occurs with 90% probability. h I and II denote the first and second solvation shells, respectively. a

shell water molecules in the CMD-LJ simulations are more tightly bound to Cm3þ. The plots of the running coordination number nCmO(r) in Figure 1, together with the data in the last column of Table 3, clearly indicate that the average first shell coordination number is 8 for the AIMD and QM/MM-PBE simulations and 9 for the QM/MM-PBE0, CMD-3B, and CMDLJ, respectively. We did not observe the experimental results by Lindqvist-Reis et al.,6 that is, 10% 8-fold and 90% 9-fold first coordination near room temperature. Rather, from the statistics of the long-time MD simulations, we obtained the occurrences of 8-fold, 9-fold, and 10-fold first shell coordination in CMD-3B simulation to be 0.2%, 99.3%, and 0.5%. For the CMD-LJ simulation, we observed no occurrences of 8-fold first shell coordination; the statistics corresponded to 99.7% 9-fold coordinated and 0.3% 10-fold coordinated first shell. Assuming that the energy barrier for a change in coordination number is low, then the CMD-3B simulation predicts the free energy, ΔG9f8, required to change from 9-fold first shell coordination to 8-fold to be ΔG9f8 = kBT ln(0.2%/99.3%) ≈ 3.7 kcal/mol. Similarly, the free energy change from 9-fold coordination to 10-fold coordination, ΔG9f10 is 3.2 and 3.5 kcal/mol for the CMD-3B and CMD-LJ simulations, respectively. By defining an appropriate reaction coordinate for the coordination number, we can estimate the change in free energy for a change in the coordination number from constrained first principles simulations. The

free energy analysis is currently not within the scope of this work and will be addressed in a separate communication. An average static measure that can be used to gauge the angular flexibility of water molecules around the metal ion in the first hydration shell is the average tilt angle ψtilt of the water molecules with respect to the Cm3þO bond. It is often measured in neutron diffraction experiments. First shell water molecules with large tilt angles tend to have a tetrahedral hydrogen-bonding structure whereas molecules with small tilt angles tend to have a trigonal hydrogen-bonding structure. Listed in Table 3 are the values of the average first shell water tilt angle, ψtilt, per simulation. As can be seen from Table 3, the AIMD and QM/MM simulations have similar ψtilt values of 30. Furthermore, all the AIMD and QM/MM MD tilt angles are larger than the corresponding values from the classical simulations (23 for CMD-3B and 16 for CMD-LJ). We note that the tilt angle for the CMD-LJ is nearly twice as small as that of the AIMD and QM/MM MD simulations, implying that CMD-LJ first shell is rigid. This is in agreement with the observation made in the RDF analysis above. Furthermore, the CMD-3B tilt angle is much closer to the QM tilt angles compared to that of CMD-LJ, indicating the importance of including three-body and higherbody terms in the interaction potential for highly charged metal ions in solution. The differences in the tilt angles imply that the first shell waters in the AIMD and QM/MM simulations 4669



ARTICLE

Figure 3. First shell OCmO bond angle distribution for 8-fold coordinated structure (top) and 9-fold coordinated structures. The delta functions correspond to the angular distribution of the ideal crystals.

Figure 2. (a) Two views of the first shell square antiprism (SAP) geometry of [Cm(OH2)8]3þ from the AIMD simulation. The face containing Oi (i = 1, 4) are rotated by approximately 45 with respect to the face containing Oi (i = 5, 8) (b) Two views of the first shell tricapped trigonal prism (TTP) geometry of [Cm(OH2)9]3þ from the CMD-3B simulation. Oi (i = 1, 6) are located at the prismatic positions and Oi (i = 7, 9) are located at the capping positions.

experience ionwater interactions and hydrogen-bonding interactions with the second shell and bulk waters, which differs from the CMD simulations. Furthermore, the large values of the tilt angles in the AIMD and QM/MM MD simulations suggest that the hydrogen-bonding features of the first shell water molecules exhibit some tetrahedral characteristics. On the other hand, the first shell hydrogen-bonding network in the CMD simulations is expected to have trigonal characteristics, particularly in the CMD-LJ simulation (due mainly to its rigid first shell geometry). Analysis of the hydrogen-bonding character in the networks is presented is section 3.4. The average CmO distances of 4.674.71 Å in the second solvation listed in Table 3 are fairly close to the corresponding experimental value by Skanthakumar et al.7 and the result of Yang and Bursten,12 both of which are 4.65 Å. The second shell CmO distance of 4.9 Å due to Hagberg et al.13 is quite large. Also from Table 3, we observe that the average secondary shell coordination numbers for all our simulations in the 15.817.6 range and that of previous calculations (16 and 21) are overestimated, compared to the experimental value of 13(4). Since counterions are employed in experiments and not in theoretical calculations, the discrepancy between the simulated and experimental secondary shell coordination numbers is possibly due to counterion presence in the second shell. The average separation between the first and second shell oxygen atoms, rOIOII, is slightly overestimated compared to the experimental data. 3.2. Geometry of the 8- and 9-Fold Coordinated First Shell Structures. As pointed earlier, the isoelectronic correspondence between the actinides and lanthanides, as well as the similarities in

their ionic radii, implies that there are similarities in the geometries of their solvent shell structures. For triply charged lanthonoids, it is well-known that the typical 8-fold coordinated stoichiometries are the square antiprism (SAP), bicapped trigonal prism (BTP), and triangulated dodecahedron (TD), while the typical 9-fold stoichiometries are the tricapped trigonal prism (TTP) and the capped square antiprism (CSAP).9,3638 We observed that SAP is the dominant 8-fold stoichiometry and TTP is the dominant 9-fold stoichiometry in our simulations and only a few BTP and CSAP structures were found. Due to strong thermal effects on the atomic motions, however, ideal crystal units could hardly be identified, and the structures we identified were either mildly or severely distorted structures. Examples of a mildly distorted SAP structure from the AIMD simulations and TTP structure from the CMD simulation are shown in Figure 2. The figures were generated using the VMD software.39 An ideal SAP is simply obtained from a cubic prism by rotating one flat square face rotated through 45 with respect to the other opposite parallel face. Assuming that the SAP in Figure 2a is ideal, the two such faces will correspond to O1O2O3O4 and O5O6O7O8. An ideal TTP is characterized by six oxygen atoms located at the prismatic positions at equal distances from the ion and three oxygen atoms at equal distances from the ion emanating out of the capping faces and lying in the equatorial plane. Assuming that the TTP structure in Figure 2b is ideal, then O1, O2, O3, O4, O5, O6 are the prismatic positions (O1O2O3 and O4O5O6 both equilateral triangles) and O7, O8, O9 are the capping positions ( — O7CmO8 = — O7 CmO9 = — O8CmO9 = 120). The first shell OCmO angular distribution function (ADF) can be used to identify the geometrical characteristics of the first solvation shell structure. The ADF provides a basis for determining the dominant structural units that make up the first solvation shell, albeit mostly distorted relative to the ideal geometry. Specifically, we can predict the dominant first hydration shell stoichiometries (8-fold [Cm(OH2)8]3þ and 9-fold [Cm(OH2)9]3þ) by comparing the ADF from the simulations to ideal crystal stoichiometries. Our claim that the dominant 8- and 9-fold stoichiometries are SAP and TTP, respectively, is evidenced by the ADF plots in Figure 3. In the top panel, the ADF 4670



ARTICLE

Figure 4. Calculated and experimental plots of the k3-weighted EXAFS spectra. The experimental result is due to Skanthakumar et al.7

Figure 5. Calculated and experimental plots of the modulus of the k3weighted Fourier transform EXAFS spectra.

for AIMD simulation (with first shell coordination number CN = 8), the QM/MM-PBE0MD simulations (with CN = 9), and the ideal SAP (CN = 8) are shown. A similar plot is shown in the lower panel, except that the ADF of the ideal TTP (CN = 9) is shown. In the top panel, it can clearly be seen that the SAP crystal peaks match the AIMD curve with the same CN = 8 but show a marked deviation from the QM/MM-PBE0 curve with CN = 9. Conversely, in the lower panel, the TTP crystal peaks match the QM/MM-PBE0 curve with CN = 9 and deviates from the AIMD curve with CN = 8. In both plots, the ideal crystal peaks in the 90120 range, which are the least frequent, are not immediately visible in the simulated curves, due to thermally induced distortions. The ADF plots for the other simulations yielded similar results. 3.3. Comparison of the Simulated EXAFS Spectra to Experimental Data. Using the MD-EXAFS method4042 we generated the EXAFS data of our simulated structures and compared them to the experimental data due to Skanthakumar et al.7 In the MD-EXAFS method, the Cm and O atomic coordinates of each of the snapshots collected over the production times per each simulation were used as input to calculate the EXAFS spectra using the FEFF8 ab initio scattering code.43,44 The resultant spectra were combined to create the ensemble average spectra. In Figure 4, the simulated k3-weighted EXAFS spectra k3χ(k) is shown for each simulation, with the experimental data superimposed on each plot for purposes of comparison. We observe that the AIMD, both QM/MM MD, and CMD-3B simulations clearly match the experimental oscillating frequencies and amplitudes over a large k range, while the amplitudes of the CMD-LJ simulation show some deviation from experiment in the k = 512 Å1 range. The fact that the CMD-3B simulations, which employ simple nonpolarizable pair and three-body potentials, yield EXAFS spectra that match experimental ones fairly well makes the potential model a fairly good one for aqueous phase simulation. Figure 5 depicts the modulus |χ h(R)|3 and the imaginary part Im[χ (R)] of the Fourier transform of k χ(k). The AIMD h

results show a good match with experiment, followed closely by the QM/MM results. The CMD-3B results show a fairly good match, but the CMD-LJ results show a marked deviation from experiment. The differences between the CMD-3B and CMD-LJ results emphasize the inadequacy of pure pair potentials in accurately describing the hydration shell characteristics of highly charged heavy ions. There is no apparent difference between the QM/MM-PBE results and the QM/MM-PBE0 results, despite the fact that they have different first shell coordination numbers. The AIMD simulation shows the best overall agreement with the EXAFS data, but the first shell coordination number is 8, while the coordination number measured from EXAFS experiment is 9. From this observation we cautiously conclude that either (i) both 8- and 9-fold coordinated first shell of Cm3þ ions coexist in aqueous solution (more specifically, the coordination number fluctuates between 8 and 9) or (ii) the EXAFS measurement is not too sensitive the small changes in the coordination number. This is plausible since the CMD simulations predicted the change in the free energy required to change the first shell coordination number from 9 to 8 to be ∼3.7 kcal/mol. 3.4. Hydrogen-Bond Structure in the Hydration Shells. A hydrogen-bond structural analysis in each network was performed to further evaluate the structural arrangement of the water molecules near and away from the Cm3þ ion. The criteria for a hydrogen-bond formation were (i) the maximum Odonor Oacceptor separation is 3.2 Å and (ii) the minimum OdonorH Oacceptor angle is 140, where Odonor and Oacceptor denote the donor and acceptor oxygens, respectively; H is the hydrogen atom bonded to Odonor. The analyses were performed by computing the acceptor bond distribution per water, nacceptor, for the first, second, and bulk solvent shells. Specifically, the procedure for computing nacceptor at a molecular site at a distance r from Cm is as follows: (1) At each postequilibration time step i, the total number of water molecules Ni whose O atoms are at a distance r from Cm3þ are located. (ii) Then for each of the Ni molecules calculate the total number of the acceptor hydrogen 4671



ARTICLE

Figure 6. Distribution of the acceptor hydrogen bonds per water located at a distance r from Cm with first shell, second shell, and bulk water for (a) AIMD, (b) QM/MM-PBE, (c) QM/MM-PBE0, (d) CMD-3B, and (e) CMD-LJ.

bonds (hij) with all the water molecules located in the shell of interest, where Oacceptor belongs to one of the Ni water molecules and Odonor belongs to a water molecule in the shell of interest. The total number of acceptor hydrogen bonds at time step i is i Hi = ∑N j=1hij. (iii) Finally nacceptor is computed as T P

nacceptor ðrÞ ¼

i¼1 T P i¼1

Hi ð5Þ Ni

where T is the total number of time steps. Figure 6 depicts a shell-by-shell acceptor hydrogen bond per water distribution for all the simulations. For the purposes of discussion, the corresponding CmO RDF for each simulation is superimposed on each plot (curve with shaded region). Except for subtle differences, all simulations display quite similar results for the acceptor hydrogen bond distribution with the first shell waters (red blocks). The value of nacceptor for water molecules forming acceptor bonds with the first shell water molecules (red dots) has an average near 1 in the 3.5 Å< r < 4.2 Å range, which implies that, on the average, a water molecule in the second hydration shell is bound to a proton in the first shell. The AIMD and QM/MM MD first shell water molecules form acceptor bonds with the second shell water (green dots in the 2 Å < r < 3 Å range in Figure 6ac), while no such first shell bonds can be seen in the plots for the CMD simulations. This implies that QM treatment of the first shell results in a mix of a few tetrahedral and dominant trigonal second shell structures, while the CMD simulation results in a mainly trigonal structure. The observed tetrahedral structure is in accordance with the large tilt angles in the AIMD and QM/MM simulations, since large tilt angles provide the first shell oxygen atoms the freedom to form acceptor hydrogen bond with a second shell donor. An example of the second

Figure 7. (a) Tetrahedral hydrogen-bonding network structure of a snapshot from the AIMD simulation (a limited number of second shell water molecules are depicted). Black bonds denote hydrogen bonds between a first shell acceptor and second shell donors; purple bonds denote a second shell water molecule forming two acceptor bonds with first shell water molecules. (b) Trigonal hydrogen-bonding network structure of the second shell water molecular for a snapshot from the QM/MM-PBE0 simulation. First shell, second shell, and bulk oxygen are colored red, orange, and gold, respectively.

shell tetrahedral hydrogen bond structure from the AIMD simulation is shown in Figure 7a, where the first shell is partially clothed by a constant density contour. In Figure 7a, the black bonds denote hydrogen bonds formed between first shell acceptor water molecules and second shell donors. The purple bonds denote hydrogen bonds formed between a second shell acceptor and two first shell donors. The purple bonds explain why the acceptor bond distribution with first shell donors (red dots in Figure 6) is ∼1.5, whereas it 1.21.3 in the other simulations. Because of limited cell size, the complete distribution of the number of acceptor bonds with donors from the bulk (blue dots in the r > 5 Å region) cannot 4672



ARTICLE

Table 4. Computed Average Dipole Moments of Water Molecules in the First, Second, and Bulk Solvation Shells Using the Maximally Localized WannierBoys Methoda dipole moment (D) system size

method

CN1

M O 3þ

first shell H2O

b

bulk H2O

3.6 ( 0.0

[Cm(OH2)8]3þ(g)

PBE

8

2.51

[Cm(OH2)9]3þ(g)

PBE

8

2.56

3.4 ( 0.0

Cm3þ þ 64H2O Cm3þ þ 64H2O

AIMD QM/MM-PBE

8 8.5

2.50 2.52

3.6 ( 0.1 3.4 ( 0.2

Cm3þ þ 64H2O

QM/MM-PBE0

9

2.52

3.2 ( 0.2

La3þ þ 64H2Oc

AIMD

8.5, 9

2.58

3.5

Y3þ þ 64H2Od

AIMD

8

2.38

3.6 ( 0.3

64H2O

AIMD

a

second shell H2O

2.9 ( 0.3

2.9 ( 0.3

2.9 ( 0.3 2.9 ( 0.6e, 3.0f

3þ

3þ

Except for the dipole moments of the 0 K gas-phase structures [Cm(OH2)8] (g) and [Cm(OH2)9] (g), all the dipole moments were computed for the aqueous phase at 300 K. CN1 denotes the first shell hydration number of the metal ion. b M3þO denotes the average first shell trivalent metal ionoxygen distance in (Å). c Reference 38. d Reference 51. e Reference 52. f Reference 53.

be fully shown for the AIMD and QM/MM simulations. However, the picture for the CMD simulations clearly indicate that the acceptor bond distribution increases with distance and saturates at around 1.6, which is quite close to the value for bulk water. It is surprising that the experimentally measured average second shell coordination number is 13 while our simulations yield ∼1617. If the first solvation shell accommodates nine water molecules, then the maximum number of second shell waters forming acceptor bonds with the first shell waters is 18 (perfect acceptor structure) and the minimum number is nine (in this case, each second shell water forms two acceptor bonds with a first shell water). As pointed out earlier, counteranions were employed in experiments while the present simulations do not. The fact that the experimental second shell coordination is 13 implies that on the average, four or five s shell waters form two acceptor bonds with the first shell water. It is likely that such a hydrogen-bond network is due to the presence of counterions in the second shell. That is, the presence of counterions in the shell alters the hydrogen-bonding network between the first and second shells and consequently reduces the number of second shell water molecules. The effect of perchlorate counterions on the hydration shell properties is currently being investigated, and preliminary results indicate that the second shell coordination number is reduced when the perchlorate ions reside in the shell. 3.5. Electronic Structure of the Solvating Waters. Using the maximally localized WannierBoys orbital technique,4550 we estimated the average dipole moment of the first shell water molecules in the AIMD and QM/MM simulations as well as the average dipole moments of the water molecules in the second and bulk solvation shells in AIMD simulations. For each simulation, 20 snapshots (equally spaced over the entire simulation time) were used to estimate the dipole moment. The results are presented in Table 4. Also listed in Table 4 are the theoretical estimates of dipole moments for the first shell water molecules of þ3 lanthanide metal ions La3þ,38 Y3þ,51 and also that of bulk water.52,53 The data suggest that the electric field of the þ3 metal ions is strong enough to induce dipole moments of 0.20.6 D in the first shell water relative to bulk water. Regarding the results obtained from our simulation, we observe that the AIMD first shell water dipole moments is the largest, followed by QM/MMPBE and then QM/MM-PBE0. These results are consistent with [Cm(OH2)8]3þ and [Cm(OH2)9]3þ gas-phase values. No

Figure 8. Time evolution of the instantaneous first shell coordination number (CN).

dipole moments are induced beyond the first shell in the AIMD simulation. Examination of the difference charge density distribution revealed a similarity between the gas-phase first shell and the solution-phase first shell. Furthermore, the ionwater interaction is dominated by electrostatic (chargecharge, charge dipole, and dipoledipole) effects in the first shell. This is evidenced by the similarities between the gas-phase and liquid dipole moments. We infer that the presence of the second shell and bulk water does very little to alter the dipole moments of the first shell. 3.6. Dynamical Properties. Water exchanges between the first and second shell were observed in the classical MD simulations. No exchange was observed in the AIMD and QM/MM MD simulations, since the time scale on which exchange processes occur is much longer than the simulation time. To gain initial insights into the nature of the exchange processes, we plotted the instantaneous first shell coordination number (CN) of each CMD simulation with respect to time. This is depicted in Figure 8. An associative exchange event between the first and second shell corresponds to a temporary increase in the first shell 4673



Figure 9. (Top) Trajectories of the molecules involved in the dissociative exchange mechanism in the CMD simulation. OIN is the incoming molecule and OOUT is the outgoing molecule. (Bottom) Visual representation for the exchange process: (A) pre-exchange geometry, (B) intermediate exchange geometry, and (C) postexchange geometry.

coordination number followed by a restoration to the original coordination number. A dissociative exchange event corresponds to a temporary decrease in the first shell coordination number followed by a restoration to the original coordination number. An intermediate exchange event corresponds to simultaneous particle interchange, leading to coordination number preservation. Clearly, Figure 8 can provide information on only associative and dissociative exchanges; it does not provide any information on intermediate exchanges. We verified separately that intermediate first exchanges occur in both CMD simulations. The top panel of Figure 8 suggests that all three types of exchange occurred in the CMD-3B simulations, while no dissociative exchanges occurred in the CMD-LJ simulation. If we define t* to be the minimum period for which an exchange event continuously persists, then the number of first shell events occurring in the CMD-3B and CMD-LJ simulations were 196 and 54, respectively, for t* = 0 ps (t* = 0 ps corresponds to all instantaneous crossing across the boundary separating the first and second shell). For a much longer persistence time of t* = 2 ps, the number of first shell events in the CMD-3B and CMD-LJ simulations reduce to 38 and 14, respectively. The second shell events in all five simulations consist of all the different exchange mechanisms. Figure 9 shows a typical dissociative water molecule exchange mechanism between the first and second shell in the CMD-3B simulation. The graph in the top panel shows the variation of the CmO distances of the incoming and outgoing water molecule over a specified period of time. We have used configurations at selected times A, B, and C to illustrate the first shell geometrical

ARTICLE

Figure 10. Same as Figure 9 but an associative exchange process.

changes involved in the exchange mechanism. Initially in scenario A, the incoming water molecule (OIN) is in the second shell and the outgoing water molecule (OOUT) is located at the equatorial capping position in a first shell structure in a distorted TTP stoichiometry. The two molecules are nearly diametrically opposite relative to the ion (angle OINCmOOUT is 170). The distance between Cm and OOUT is the longest among first shell CmO distances. After 1.4 ps (scenario B), OOUT leaves the first shell into the second shell, with OIN still in the second shell, and the resulting 8-fold first shell structure temporarily reorganizes into a mildly distorted SAP stoichiometry. In the next picosecond (scenario C), OIN relocates into the first shell and the stoichiometry changes back into a distorted TTP. We thus predict from the classical simulations that the first shell 9 f 8f9 dissociative ligand exchange proceeds via a TTP f SAPfTTP first shell geometry transformation. Figure 10 shows a typical first shell associative ligand exchange mechanism in the CMD simulations. Initially in A, the incoming water molecule (OIN) is in the second shell and the outgoing water molecule (OOUT) is located in the capping position in a first shell distorted TTP structure. As was the case in the dissociative mechanism, the two molecules were initially nearly diametrically opposite relative to the ion, with the OINCmOOUT angle being 180. After 1.8 ps (scenario B), OIN relocates into the first shell, resulting in a coordination number of 10, with OOUT still within the first shell. The intermediate 10-fold first shell structure corresponds to a nearly ideal bicapped SAP (SAP with two waters capping the 4674



ARTICLE

Table 5. First and Second Shell Direct (τt*D) and Probabilistic (τt*I ) Mean Residence Times for t* = 0.4, 0.8, 1.2, 1.6, 2.0 ps t* = 0.4 ps

t* = 0.8 ps τ0.8 D

τ0.8 I

t* = 1.2 ps τ1.2 D

τ1.2 I

t* = 1.6 ps τ1.6 D

t* = 2 ps

τ0.4 D

τ0.4 I

τ1.6 I

τ2.0 D

τ2.0 I

CMD-3B

0.60

1.52

0.90

1.56

0.95

1.56

CMD-LJ

2.25

3.97

2.25

3.97

2.25

3.97

0.95

1.56

0.95

1.56

2.57

3.97

2.57

3.97

AIMD

3.9

8.4

9.6

10.3

12.5

QM/MM-PBE

4.9

8.8

10.6

12.5

13.1

QM/MM-PBE0

5.4

CMD-3B

3.9

8.9

5.8

10.6

7.1

11.7

8.2

12.4

9.1

13.1

CMD-LJ ref 12

5.1

11.9

7.2

13.6

8.4

14.6

9.4

15.3

10.2

15.8 161.0

Mean Residence Times in the First Solvation Shell (in nanoseconds)

Mean Residence Times in the Second Solvation Shell (in picoseconds)

8.6

10.9

faces). Finally, after 2 ps (scenario C), OOUT departs the first shell and the first shell rearranges into a distorted TTP structure. Hence, the first shell associative ligand 9f10f9 exchange mechanism corresponds to a TTPfbicapped SAPfTTP structural transformation. Other associative and dissociative mechanisms that we examined proceeded via similar routes. Clearly, the intermediate structures (SAP for the dissociative exchange and bicapped SAP for the associative mechanisms) are short-lived, with a typical lifetime of ∼37 ps. Furthermore, associative or dissociative exchange mechanisms involved a capping water in the TTP structure with the longest CmO bond distance. The mean residence time (MRT) is the average time a ligand resides in a coordination shell of the solute. For all simulations, we computed the MRT τt*D using the so-called “direct approach” by Hofer and co-workers.54 For the long-time CMD simulations, the probabilistic approach by Impey55 and co-workers was also used to compute the MRT τt*I . τt*I was not computed for the AIMD and QM/MM simulations due to the lack of sufficient statistics. The description of the procedure used to compute the MRT is provided in the Supporting Information. As mentioned earlier, the time parameter t* in both methods is used to prevent short-lived visitations to and departures from a shell from being counted as real events, that is, only visitations persisting continuously for a period longer than t* are used to compute the MRT. The question of which value of t* should be used to calculate the MRT is not welldefined in the literature. Impey and co-workers55 suggested a value of t* = 2.0 ps based on their calculation which yielded a mean first passage time of 1.8 ps in bulk water. On the basis of systematic QM/ MM-HF studies of several cations in solution, Hofer et al.54 proposed a value of t* = 0.5 ps. We computed τt*D and τt*I for t* = 0.4, 0.8, 1.2, 1.6, and 2.0 ps. The results are reported in Table 5. From the first shell MRT in CMD simulations we see, for all configurations, that τt*I converges to 1.56 ns (3.97 ns) for the CMD3B (CMD-LJ) simulation when t* is at least 0.8 ps, whereas τt*D converges to 0.95 ns (2.57 ns) when t* is at least 1.6 ps. We also observe that τt*D is consistently less than τt*I , but their individual trends are similar. The observed differences between the two methods have been observed by Hofer et al.54 Over all, the first shell MRT for Cm3þ is on the nanosecond time scale. Because of the overall better performance of CMD-3B potential over CMD-LJ (as gauged by the first shell static properties and comparisons to EXAFS data), we deem the MRT for CMD-3B to be more reliable. There is no experimental value for the mean residence time for Cm3þ. However, Gd3þ, which is isoelectronic to Cm3þ, has an

13.1

14.4

experimentally measured first shell MRT of 0.8356 and 0.943 ns,57 which appears to be on the same time scale predicted by our simulations. We must also mention that the experimentally measured MRT for other 9-fold coordinated þ3 lanthanoids are also on the nanosecond time scale: 2.02 ns for Tb,58 2 ns for Pr,59 and 2.5 for Nd.59 As was the case for the first shell, the mean residence times for second shell indicate that τt*D is consistently less than τt*I . More astonishingly, our largest reported value of τt*I (15.8 ps) is an order of magnitude less than the simulated value of 161 ps by Yang and Bursten12 using the rigid hydrated ion model and TIP3P water potential. We believe that the use of the rigid hydrated ion model in ref 12 introduces serious artifacts in the simulation and therefore cast doubts on the reliability of the rather large residence time. The diffusion coefficient, which provides information on the mobility of ions in solution, allows us to check the validity of the classical MD potentials. Intuitively, we expect the Cm(III) ion to diffuse more slower in CMD-LJ simulations compared to CMD3B, since the CMD-LJ first shell is quite rigid. The self-diffusion coefficient Di of species i was computed from the slope of the mean square displacement in the long time limit * + Ni 1 d 1X 2 Di ¼ lim jri ðt þ t0 Þ ri ðt0 Þj 6 t f ¥ dt Ni i ¼ 1

ð6Þ t0

where Ni is the number of species i present in the solution, Æ æt0 denotes an average over all time origins t0, and ri(t) is the center of mass position of molecule/atom i at time t. The diffusion constant was computed by using the algorithm by Frenkel and Smit.60 In Table 6, the self-diffusion coefficients of Cm3þ and water are reported along with the corresponding experimental values. The experimental diffusion constant of Cm3þ at 298 K was determined in an aqueous solution of Gd(ClO4)3HClO4 (pH 2.5) using the open end-capillary method (OECM).61 Clearly the CMD-3B value of 0.52 105 cm2/s matches fairly well with the experimental value of 0.60 105 cm2/s, while the values from the CMD-LJ simulations are smaller than experiment. As expected, the diffusion of Cm3þ obtained from then CMD-LJ simulation has a smaller value of 0.38 105 cm2/s. Also all the computed diffusion constants of water match well with the experimental value of 2.3 105 cm2/s.62,63 The diffusion constants for the AIMD and AIMD/CMD were not reported because we not do believe that the total simulation time was long enough to reach the diffusive regimeparticularly for 4675



ARTICLE

Table 6. Calculated and Experimental Diffusion Coefficients (in units of 105 cm2/s) for Cm3þ (DCm3þ) and Water (DH2O) method

a

DCm3þ

DH2O

CMD-3B

0.52

2.47

CMD-LJ expt

0.38 0.60a

2.46 2.30b,c

Reference 61. b Reference 62. c Reference 63.

Cm3þ, which diffuses at a much slower rate compared to a water molecule.

4. SUMMARY AND CONCLUSIONS In this paper, AIMD, QM/MM MD, and CMD simulations of the curium(III) ion in liquid water have been reported. The AIMD and QM/MM MD simulations employed cubic periodic cells containing 64 water molecules, while the CMD simulation employed periodic cells containing 512 water molecules. The simulations were carried out at in the NVT ensemble at T = 300 K, and lengths of the cubic cells were chosen to correspond to a water density close to 1 g/cm3. One CMD potential comprising of pair potential with three-body corrections was developed from ab initio ROHF, while the other potential was based on purely empirical pair Lennard-Jones potential. The AIMD and QM/MM-PBE simulations yielded a first shell hydration number of 8, while the QM/MM-PBE0 and CMD simulations yielded a first shell coordination number of 9. While both coordination numbers have been observed in experiments, the most recent experiments indicate a first shell hydration number of 9. However, the simulated second shell hydration number was overestimated by ∼24 more than the experimental value of 13, with the low experimental value attributed to the presence of counterions in the second shell. All simulations reproduced the first and second hydration shell distances fairly well. Specifically, the simulated average first shell CmO bond distances of 2.502.53 Å compared well to the recent EXAFS value of 2.48 Å, while the simulated average second shell CmO bond distances of 4.674.75 Å compared well to the recent HEXS value of 4.65 Å. The simulations predicted that the 8- and 9-fold first coordinated geometries were square antiprism and tricapped trigonal pyramid, respectively. The AIMD and QM/MM simulations yielded an average water dipole moment tilt angle of 30, the CMD-3B yielded a tilt angle of 23, and the CMD-LJ simulations yielded a rigid first shell with a tilt angle of 16. The structure of the second shell water molecules is mainly trigonally coordinated via hydrogen bonds to the first shell water molecules, with QM treatment of the first shell exhibiting some tetrahedral characteristics. The computed EXAFS spectra from the CMD-LJ simulation showed a significantly poor match with recent EXAFS experimental data, while all other computed EXAFS spectra showed a decent match with experiment. In comparison to the experimental EXAFS, it thus appears that the simulated EXAFS spectra are much more sensitive to the tilt angle rather than the first shell hydration number. Specifically, there is no clear distinction between the computed EXAFS spectra for simulations with average first shell coordination numbers of 8 or 9 in comparison to the experimental EXAFS data. The insensitivity of the simulated EXAFS to the first shell coordination suggests the

likely statistical coexistence of 8- and 9-fold coordinated structures in aqueous solution. First shell exchange events were identified and characterized in CMD simulations. Associative exchange mechanisms (coordination 9fcoordination 10fcoordination 9) occurs via TTP fbicapped SAPfTTP geometry transformation, while dissociative exchange processes (coordination 9fcoordination 8fcoordination 9) occurs via TTPfSAPfTTP geometry transformation. The intermediate structures in exchange process have a lifetime less than 7 ps. CMD-3B simulations yielded a diffusion constant of Cm3þ, which is in good agreement with experiment. On the basis of a persistence time of t* = 2 ps, the mean ligand lifetime in the first shell was predicted to be in the 0.951.56 ns range, while the mean second shell residence time was predicted to be in 9.115.8 ps range. The electronic structure of analysis of the AIMD and QM/ MM simulations indicated that the Curium ion strongly polarizes the first shell water molecules, causing the dipole moments to increase by ∼0.30.5 D relative to bulk water. Similarities were found between gas-phase electronic structures [Cm(OH2)8]3þand [Cm(OH2)9]3þ and the first shells of the AIMD and QM/MM simulations. respectively, implying that the Cmwater interaction is confined mainly to the first shell.

’ ASSOCIATED CONTENT

bS

Supporting Information. Strategy for obtaining model potential parameters for the CMD simulations, description of the MRT computation, and the FORTRAN code for the generating initial structures. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*Phone: 1-509 371 6472. Fax: 1-509 376 3650. E-mail: raymond. [email protected].

’ ACKNOWLEDGMENT R.A.-F. would like to thank D. F. Johnson and C. Mundy for stimulating discussions and N. Govind for technical help. The authors would like to acknowledge Patrick Nichols for providing the initial pseudopotential for Cm3þ and L. Soderholm and her group for graciously making the experimental EXAFS data available to us. This research was performed using the Molecular Science Computing Capability in the William R. Wiley Environmental Molecular Science Laboratory, a national scientific user facility sponsored by the U.S. Department of Energy’s Office of Biological and Environmental Research and located at the Pacific Northwest National Laboratory, operated for the Department of Energy by Battelle. This work was supported by the BES Heavy Element Chemistry program and the Division of Chemical Sciences, Geosciences, and Biosciences, Office of Basic Energy Sciences, U.S. Department of Energy. ’ REFERENCES (1) Katz, J. J, Seaborg, G. T., Morss, L. R., Eds. The Chemistry of Actinide Elements, 2nd ed.; Chapman and Hall: New York, 1986; Vol. 2. (2) Edelstein, N. M.; Klenze, R.; Fangh€anel, T.; Hubert, S. Coord. Chem. Rev. 2006, 250, 948. 4676


The Journal of Physical Chemistry A (3) Kimura, T.; Choppin, G. R. J. Alloys Compd. 1994, 213-214, 313. Kimura, T.; Choppin, G. R.; Kato, Y.; Yoshida, Z. Radiochim. Acta 1996, 72, 61. (4) Allen, P. G.; Bucher, J. J.; Shuh, D. K.; Edelstein, N. M.; Craig, I. Inorg. Chem. 2000, 39, 595. (5) Stumpf, T.; Fangh€anel, T.; Grenthe, I. J. Chem. Soc., Dalton Trans. 2002, No. 20, 3799. (6) Lindqvist-Reis, P.; Klenze, R.; Schubert, G.; Fanghanel, T. J. Phys. Chem. B 2005, 109, 3077. (7) Skanthakumar, S.; Antonio, M. R.; Wilson, R. E.; Soderholm, L. Inorg. Chem. 2007, 46, 3485. (8) Vertes, A., Nagy, S., Klencsar, Z., Eds. Handbook of Nuclear Chemistry; Kluwer Academic Publishers: Dordrecht, Netherlands, 2003, Vol. 2. (9) Richens, D. T. The Chemistry of Aqua Ions, 1st ed.; Wiley: Chichester, UK, 1997. (10) Mochizuki, Y.; Tatewaki, H. J. Chem. Phys. 2002, 116, 8838. (11) Wiebke, J.; Moritz, A.; Cao, X.; Dolg, M. Phys. Chem. Chem. Phys. 2007, 9, 459. (12) Yang, T.; Bursten, B. E. Inorg. Chem. 2006, 45, 5291. (13) Hagberg, D.; Bednarz, E.; Edelstein, N. M.; Gagliardi, L. J. Am. Chem. Soc. 2007, 129, 14136. (14) Car, R.; Parrinello, M. Phys. Rev. Lett. 1985, 55, 2471. (15) Valiev, M.; Bylaska, E. J; Govind, N.; Kowalski, K.; Straatsma, T. P.; van Dam, H. J. J.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T. L.; De Jong, W. A. Comput. Phys. Commun. 2010, 181, 1477 (http://www. nwchem-sw.org). (16) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865. (17) Hamann, D. R.; Schl€uter, M.; Chiang, C. Phys. Rev. Lett. 1979, 43, 1494. (18) Hamann, D. R. Phys. Rev. B 1989, 40, 2980. (19) Kleinman, L.; Bylander, D. M. Phys. Rev. Lett. 1982, 48, 1425. (20) Troullier, N.; Martins, J. L. Phys. Rev. B 1991, 43, 1993. (21) Nose, S. Mol. Phys. 1984, 52, 255. Hoover, W. G. Phys. Rev. A 1985, 31, 1695. (22) Blochl, P. E.; Parrinello, M. Phys. Rev. B 1992, 45, 9413. (23) Ernzerhof, M; Scuseria, G. E. J. Chem. Phys. 1999, 110, 5029. Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 6158. (24) Cau€et, E.; Bogatko, S.; Weare, J. H.; Fulton, J. L.; Schenter, G. K.; Bylaska, E. J. J. Chem. Phys. 2010, 132, 194502 (see also the planewave DFT section of NWchem manual at http://www.nwchem-sw. org). (25) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, Netherlands, 1981; p 331. (26) Toukan, K.; Rahman, A. Phys. Rev. B 1985, 31, 2643. (27) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (28) Inada, Y.; Mohammed, A. M.; Loeffler, H. H.; Rode, B. M. J. Phys. Chem. A 2002, 106, 6783. (29) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A., III; Skiff, W. M. J. Am. Chem. Soc. 1992, 114, 10024. (30) Martyna, G. J.; Klein, M. L.; Tuckerman, M. E. J. Chem. Phys. 1992, 97, 2635. Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Mol. Phys. 1996, 87, 1117. (31) Yoshida, Y. Phys. Lett. A 1990, 150, 262. (32) Suzuki, M. J. Math. Phys. 1991, 32, 400. (33) Ewald, P. Ann. Phys. 1921, 369, 253. (34) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1989. (35) Andersen, H. C. J. Comput. Chem. 1983, 52, 24. (36) Kowal, Th.; Foglia, F.; Helm, L.; Merbach, A. E. J. Phys. Chem. 1995, 99, 13078. (37) Duvail, M.; Souaille, M.; Spezia, R.; Cartailler, T. J. Chem. Phys. 2007, 127, 034503. Duvail, M.; Vitorge, P.; Spezia, R. J. Chem. Phys. 2009, 130, 104501.

ARTICLE

(38) Terrier, C.; Vitorge, P.; Gaigeot, M.-P.; Spezia, R.; Vuilleumier, R. J. Chem. Phys. 2010, 133, 044509. (39) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33 (http://www.ks.uiuc.edu/Research/vmd/). (40) Palmer, B. J.; Pfund, D. M.; Fulton, J. L. J. Phys. Chem. 1996, 100, 13393. (41) Dang, L. X.; Schenter, G. K.; Glezakou, V. A.; Fulton, J. L. J. Phys. Chem. B 2006, 110, 23644. (42) McCarthy, M. I.; Schenter, G. K.; ChaconTaylor, M. R.; Rehr, J. J.; Brown, G. E. Phys. Rev. B 1997, 56, 9925. (43) Zabinsky, S. I.; Rehr, J. J.; Ankudinov, A.; Albers, R. C.; Eller, M. J. Phys. Rev. B 1995, 52, 2995. (44) Rehr, J. J.; Albers, R. C. Rev. Mod. Phys. 2000, 72, 621. (45) Berghold, G.; Mundy, C. J.; Romero, A. H.; Hutter, J.; Parrinello, M. Phys. Rev. B 2000, 61, 10040. (46) Foster, J. M.; Boys, S. F. Rev. Mod. Phys. 1960, 32, 300. (47) Silvestrelli, P. L.; Marzari, N.; Vanderbilt, D.; Parrinello, M. Solid State Commun. 1998, 107, 7. (48) Wannier, G. H. Phys. Rev. 1937, 52, 0191. (49) Marzari, N.; Vanderbilt, D. Phys. Rev. B 1997, 56, 12847. (50) Vanderbilt, D.; Kingsmith, R. D. Phys. Rev. B 1993, 48, 4442. (51) Ikeda, T.; Hirata, M.; Kimura, T. J. Chem. Phys. 2005, 122, 024510. (52) Badyal, Y. S.; Saboungi, M.-L.; Price, D. L.; Shastri, S. D.; Haeffner, D. R.; Soper, A. K. J. Chem. Phys. 2000, 112, 9206. (53) Silvestrelli, P. L.; Parrinello, M. Phys. Rev. Lett. 1999, 82, 3308. (54) Hofer, T. S.; Tran, H. T.; Schwenk, C. F.; Rode, B. M. J. Comput. Chem. 2004, 25, 211. (55) Impey, R. W.; Madden, P. A.; McDonald, I. R. J. Phys. Chem. 1983, 87, 5071. (56) Helm, L.; Merbach, A. E. J. Chem. Soc., Dalton Trans. 2002, 633. (57) Lide, D., R., Ed. CRC Hand book of Chemistry and Physics, 77th ed.; CRC Press: Boca Raton, FL, 1996. (58) Laurenczy, G.; Merbach, A. E. Helv. Chim. Acta 1988, 71, 1971. Cossy, C.; Helm, L.; Merbach, A. E. Inorg. Chem. 1988, 27, 1973. Ohtaki, H.; Radnai, T. Chem. Rev. 1993, 93, 1157. (59) Helm, L.; Merbach, A. E. Chem. Rev. 2005, 105, 1923. (60) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed.; Academic Press: San Diego, CA, 2001. (61) Latrous, H.; Oliver, J. J. Radioanal. Nucl. Chem. 1992, 156, 291. (62) Lide, D., R., Ed. CRC Hand book of Chemistry and Physics, 84th ed.; CRC Press: Boca Raton, FL, 2003. (63) Eisenberg, D. S.; Kauzmann, W. The Structure and Properties of Water; Clarendon Press: Oxford, UK, 1969.

4677


Hydration Shell Structure and Dynamics of Curium(III) in Aqueous

Recommend Documents