Comparative CHARMM and AMOEBA Simulations of Lanthanide

Subscriber access provided by UNIV OF DURHAM

Article

Comparative CHARMM and AMOEBA Simulations of Lanthanide Hydration Energetics and Experimental Aqueous-Solution Structures Baofu Qiao, S. Skanthakumar, and Lynda Soderholm J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b01018 • Publication Date (Web): 13 Feb 2018 Downloaded from http://pubs.acs.org on February 16, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Comparative CHARMM and AMOEBA Simulations of Lanthanide Hydration Energetics and Experimental Aqueous-Solution Structures Baofu Qiao,* S. Skanthakumar, L.Soderholm* Chemical Sciences and Engineering Division, Argonne National Laboratory, Argonne, Illinois 60439, USA. The submitted manuscript has been created by UChicago Argonne, LLC, operator of Argonne National Laboratory (“Argonne”). Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan. http://energy.gov/downloads/doe-public-accessplan

1 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Comparative CHARMM and AMOEBA Simulations of Lanthanide Hydration Energetics and Experimental Aqueous-Solution Structures Baofu Qiao,* S. Skanthakumar, L.Soderholm* Chemical Sciences and Engineering Division, Argonne National Laboratory, Argonne, Illinois 60439, USA.

ABSTRACT: The accurate understanding of metal ion hydration in solutions is a prerequisite for predicting its stability, reactivity, and solubility. Herein additive CHARMM force field parameters were developed to enable molecular dynamics (MD) simulations of lanthanide (Ln) speciation in water. Quantitatively similar to the much more resource intensive polarizable AMOEBA potential, the CHARMM simulations reproduce the experimental hydration free energies and correlations in the first shell (Ln-oxygen distance and hydration number). Comparisons of difference pair-distribution functions (dPDF) obtained from the two simulation approaches with those from high-energy Xray scattering (HEXS) experiments reveals good agreement of first-coordination sphere correlations for the Lu3+ ion (CHARMM only), but further improvement to both approaches is required to reproduce the broad, non-Gaussian distribution seen from the La3+ experiment. Secondcoordination sphere comparisons demonstrate the importance of explicitly including an anion in the simulation. This work describes the usefulness of less resource-intensive additive potentials in some complex chemical systems such as solution environments, where multiple interactions have similar energetics. In addition, 3-dimensional descriptions of the La3+ and Lu3+ coordination geometries are extracted from the CHARMM simulations and generally discussed in terms of potential improvements to solute-structure modeling within solution environments. KEYWORDS: metal ion; atomistic simulation; CHARMM potential; AMOEBA potential; highenergy X-ray scattering


Page 2 of 27

Page 3 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


INTRODUCTION The speciation of a metal ion in water, defined by its interactions with other solute and solvent molecules, determines the thermodynamic and kinetic properties underlying the observable aqueous chemistry.1-2 Predicting metal-ion solution behavior, notably complexes formed, their relative stabilities, reactivities and solubilities, requires a metrical knowledge of the complex free energies to accuracies of a few kcal/mol. Such requirements necessitate structural representations in not only the first coordination sphere but also the longer-range correlations, a particularly important inclusion for higher-valent metal ions.3 The complexity of treating such a large-scale system is best handled through molecular dynamics (MD) simulations and is often very resource intensive.4 It is thus important to understand the accuracy to be gained by the degree of model complexity and whether simpler, less expensive approaches will suffice for the problem under consideration. The accuracy of an MD calculation depends on the potentials employed. As part of a continuing effort to improve computational capabilities several groups are developing additive5-9 or polarizable10-17 potentials for trivalent metal ions (see review in the Supporting Information; the efforts on lower-valent metal ions have been reviewed recently. 4). Of the two, polarizable potentials are considered the more promising for reproducing experimental hydration free energies13 (or enthalpies12) and structures10-11, 13-17 due to their explicit inclusion of the atomic polarizability (and the multipole in the AMOEBA potential18).4 For instance, the AMOEBA parameters for La3+, Eu3+, Gd3+, Ac3+, Am3+ and Cm3+ ions developed by Marjolin et al.13 were able to reproduce both the experimental hydration free energies (except for Ac 3+) and the metal ion-water distance r(M3+-O) and the hydration number Nwater in the first two hydration shells. Unfortunately, such polarizable potentials are time and resource demanding and thus computationally expensive (Table S1), prohibiting their broad application for large or complex systems such as heavy-metal chemical separations.19 Even in relatively simple systems there exist only a limited number of studies quantitatively comparing polarizable with additive potentials,17, 20-21

none of which focus on Ln aqua speciation.

With this in mind, we report herein the development of additive CHARMM force-field parameters for the Ln3+ ions. We are aiming to simultaneously reproduce the experimental hydration free energies and atomic correlations in the first- and more-distant-hydration shells. We use two 3 ACS Paragon Plus Environment


Page 4 of 27

different approaches for simulating trivalent-lanthanide aqua speciation, employing either the polarizable AMEOBA13 or the additive CHARMM potentials.22-23 The simulation results are compared with available hydration free energies24 and structures.25-27 In addition, they are graphed together and compared with a high-energy X-ray scattering (HEXS) experiment.28-31 MATERIALS AND METHODS Additive CHARMM Simulations 𝑞𝑖 𝑞𝑗

𝐸𝑖𝑗 = 4𝜋𝜖

0 𝑟𝑖𝑗

+ 𝜀𝑖𝑗 [(

𝑅𝑚𝑖𝑛,𝑖𝑗 𝑟𝑖𝑗

12

)

− 2(

𝑅𝑚𝑖𝑛,𝑖𝑗 𝑟𝑖𝑗

6

) ]

(1)

In the present work, the additive CHARMM force field 22-23 was employed using the package NAMD 2.11.32 In the CHARMM potential, each monoatomic ion carries three (non-bonded) parameters, the atomic charge q and the Lennard-Jones (LJ) 12-6 parameters 𝑅𝑚𝑖𝑛 and 𝜀 (Eq. 1). For the trivalent metals, the atomic charge is fixed as q = +3e with e being the elementary charge. The LJ parameter 𝜀 represents the potential well depth and 𝑅𝑚𝑖𝑛 the distance at which the potential is minimum. The CHARMM potential uses the Lorentz-Berthelot combination rule between different types of atoms 𝜀𝑖𝑗 = (𝜀𝑖𝑖 𝜀𝑗𝑗 )

1⁄ 2

and 𝑅𝑚𝑖𝑛,𝑖𝑗 =

𝑅𝑚𝑖𝑛,𝑖𝑖 2

+

𝑅𝑚𝑖𝑛,𝑗𝑗 2

. The recommended

modified TIP3P water model33-35 was used with the structures constrained using the SETTLE algorithm.36 The CHARMM TIP3P water model was modified based on the original TIP3P water model37 with the LJ parameters of hydrogen included. Their LJ interactions with other atoms also follow the Lorentz-Berthelot combination rule. The non-bonded pair lists were searched up to the distance of 13.5 Å. The LJ potential was switched off from 10 to 12 Å. The short-range Coulomb potential was cut off at 12 Å, with the long-range interactions calculated by means of the particle mesh Ewald (PME) algorithm 38-39 under 3-dimensional periodic boundary (PBC) conditions. In the PME algorithm, a uniform background charge density (neutralizing plasma40) is used to neutralize the net charge when it is not zero. The combination of PME for the long-range electrostatic interactions and the PBC simulation box has been shown to not affect the accuracy of the calculated hydration free energy of Na+.41 The six-order interpolation was used along with the tolerance of 10-6 and the grid spacing of 1 Å for the Ewald algorithm. The NTP ensemble (constant number of particles, temperature and pressure) was employed. The Langevin dynamics was used with the temperature of 298.15 K and 4 ACS Paragon Plus Environment

Page 5 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the damping coefficient of 1 ps-1. The integration time step was 1 fs. The impulse multiple time step algorithm42 was used, where the non-bonded potentials were calculated every one-time step along with every two-time steps for the long-range electrostatic interactions. The Nosé-Hoover isotropic Langevin piston pressure coupling 43-44 was employed to reach the optimal system density at 1 bar. The initial structures were built using the package Packmol,45 where one M3+ ion with the unit mass was distributed at the center of the box with 1014 waters. The initial box edge length was 32 Å, which was converged to the equilibrium length of around 31.1 Å in the equilibration simulations. The LJ parameter set (𝑅𝑚𝑖𝑛 /2, 𝜀) was varied with 1.35 ≤ 𝑅𝑚𝑖𝑛 ⁄2 ≤ 2.0 Å with the interval of 0.025 Å and 0.01 ≤ 𝜀 ≤ 0.12 kcal/mol with the interval of 0.01 kcal/mol. For some 𝑅𝑚𝑖𝑛 /2, the value of 𝜀 was calculated up to 0.16 kcal/mol. For any given LJ parameter set (𝑅𝑚𝑖𝑛 /2, 𝜀), the energy minimization was first performed for 5000 time steps, followed by the equilibration of 20 ps. The final structure was used for the subsequent free energy perturbation (FEP) calculations 46 and the simulations for the hydration structures. a) FEP calculations for the M3+ hydration free energies In each FEP calculation, 43 intermediate states were employed: λ LJ = 0, 0.05, 0.15 …… 0.85, 0.95, 1 with λCoulomb = 0, followed by λCoulomb = 0, 1/60, 3/60, 5/60 …… 57/60, 59/60, 1 with λLJ = 1, where the corresponding interactions were turned off and on at λ = 0 and 1, respectively. The free energy difference between two adjacent intermediate states A and B is thus expressed as: 46 ∆𝐺 (𝐴 → 𝐵) = −𝑘𝐵 𝑇𝑙𝑛 〈𝑒𝑥𝑝 (−

𝐸𝐵 −𝐸𝐴 𝑘𝐵 𝑇

)〉𝐴 ,

(2)

where 𝑘𝐵 is the Boltzmann’s constant, T is the system temperature, E stands for the potential energy. The FEP method has already been implemented in NAMD. For each intermediate state, the system was first equilibrated for 20 ps, followed by another 40 ps for the FEP data collection. Here we are reporting the hydration free energies, ∆𝐺𝐹𝐸𝑃 , from the forward FEP calculations. They last approximately 400 ns in total. The convergence of the FEP calculations is discussed in the Supporting Information. The hydration free energy obtained in the conditions of PBC and PME for long-range coulombic interactions is sometimes termed as “intrinsic” free energy due to the lack of the water-vacuum interface. 47-49 To obtain the absolute hydration free energy, some corrections need to be included. 5 ACS Paragon Plus Environment


Page 6 of 27

The first is made to represent the change in the standard state concentrations from the gas phase (1 bar) to solution (1 mol/L): 48, 50-51 0.001

∆𝐺𝑠𝑠 = −𝑅𝑇𝑙𝑛 (0.024788) = 7.958 kJ/mol,

(3)

where R is the gas constant 8.314 J/mol/K, 0.024788 m3 is the volume of one mole ideal gas at the pressure of 1 bar and the temperature of 298.15 K based on the ideal gas law (𝑃𝑉 = 𝑅𝑇), and 0.001 m3 is the volume of one mole ideal gas at the concentration of 1 mol/L. Under the PBC simulation box with PME for the long-range electrostatic interactions, the following corrections have been proposed by Hünenberger and co-workers to obtain the hydration free energy52-53 ∆𝐺𝑃𝐵𝐶−𝑃𝑀𝐸 = ∆𝐺𝐴𝑁𝐴 + ∆𝐺𝐷𝑆𝐶 + ∆𝐺𝐷 .

(4)

The first contribution is the correction for the finite size effect, which is approximated as ∆𝐺𝐴𝑁𝐴 = −

𝜉𝐿𝑆 𝑞2 8𝜋𝜖0 𝜖𝑠 𝐿

≈ 6.2 kJ/mol,

(5)

where q = +3e for the Ln3+ ions, L is the cubic simulation box length (~ 31.1 Å here), 𝜖𝑠 is the solvent relative permittivity (92 for the CHARMM TIP3P water 54), 𝜉𝐿𝑆 ≈ −2.837297 is the cubic lattice-sum integration constant, (4𝜋𝜖0 )−1 = 138.93545585 kJ·nm·e-2·mol-1. The other contributions related to the finite size effect are negligible (a few kJ/mol

54

), and are thus not

included in the present work. 55 The discrete solvent contribution (∆𝐺𝐷𝑆𝐶 ) is necessary to correct for the artifacts originating from the lattice-sum Ewald electrostatics (close to the particle based P-summation 41, 52-53, 56) under the PBC conditions, and is written as ∆𝐺𝐷𝑆𝐶 = −

𝑞𝛾𝑠 𝑁𝑠 6𝜖0 𝐿3

≈ −224.8 kJ/mol,

(6)

where the quadrupole-moment trace of the TIP3P water 𝛾𝑠 = 0.00764 e·nm2, 54 and the number of water solvent 𝑁𝑠 = 1014 here. The correction to the inaccurate relative permittivity of the CHARMM TIP3P water model (𝜖𝑠 = 92), ∆𝐺𝐷 = 0.6 kJ/mol, 54 is ignored in the present work. Consequently, the hydration free energy is:


Page 7 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


∆𝐺 ℎ𝑦𝑑 = ∆𝐺𝐹𝐸𝑃 + ∆𝐺𝑆𝑆 + ∆𝐺𝑃𝐵𝐶−𝑃𝑀𝐸 ,

(7)

which can be compared with the experimental data. Different methods to correct the calculated hydration free energy of monoatomic ions have been proposed by Roux and co-workers47-49,

56

and Lu and Cui.41 Readers are referred to these

publications for details. Note that experimentally different choices of the proton hydration free energy as the reference will lead to some differences in the resulting metal-ion hydration free energy. Two decades ago, Marcus used the value ∆𝐺 ∗ (𝐻 +) = -1056 ± 6 kJ/mol.

24, 57

Later, values of -1098 kJ/mole 58 and -1112

kJ/mol59 were suggested, the latter of which was further justified by Camaioni and Schwerdtfeger60 and Truhlar and co-workers.61 In this present work we employ the value of ∆𝐺 ∗ (𝐻 +) = -1112 kJ/mol. Consequently the experimental absolute hydration free energies of the Ln 3+ ions are recalculated: 24 ∆𝐺𝑎𝑏𝑠 = ∆𝐺𝑐𝑜𝑛. − 1112 × 𝑞,

(8)

where q = +3 for the Ln3+ ions. ∆𝐺𝑐𝑜𝑛. is the conventional free energy, which could be obtained using the reported hydration free energy by Marcus (Table 2 in Ref.[ 24]) and the reference proton hydration free energy of -1056 kJ/mol. The re-calibrated hydration free energies of Ln 3+ ions are provided in Table 2 in the Discussion section. The calculations based on the proton hydration free energy of -1056 kJ/mol are provided in the Supporting Information. b) Simulation structures for hydrated M3+ions in aqueous solution Simulations were performed to calculate the hydration structures at each LJ parameter set (𝑅𝑚𝑖𝑛 , 𝜀). Each simulation ran 2 ns, which is 600 ns in total. The metal-centric PDF was calculated based on 2000 snapshots obtained after the system came to equilibrium. The optimal distances r(M3+-O) are obtained from the primary PDF peak positions. The hydration numbers, Nwater, are obtained by the summations of water oxygens within the hydration shell. Polarizable AMOEBA Simulation In the present work, the force field parameters from Dognon and co-workers were employed for the La3+ ions 13 together with the original water model,62 both in the framework of the AMOEBA force field.18,

63

The GPU version of the package OpenMM 7.0.1 7 ACS Paragon Plus Environment

64

was employed. In the


AMOEBA simulations, the mutual induced polarization was used with the convergence criterion of 10-6 Debye. The LJ interactions were calculated up to the distance of 9 Å with the long-range dispersion correction applied. The short-range electrostatic interactions were calculated up to 7 Å with the PME algorithm for the long-range electrostatic interactions. The Langevin dynamics was applied under the temperature of 298.15 K with the integration time step of 1 fs and the damping coefficient of 1 ps-1. Initially, one metal ion was embedded in the simulation box of 1014 waters with the edge length of 32 Å in each dimension. The energy minimization was firstly performed with the energy convergence criterion of 1 kJ/mol. The simulation under the NTV ensemble were done for a duration of 10 ps, followed by the NTP ensemble simulation with the Monte Carlo barostat coupling65-66 (the reference pressure 1 atm) for 40 ps. In the production simulation, 5 ns simulation was performed to collect the structural properties. High-Energy X-ray Scattering (HEXS) Experiments HEXS experiments were performed on aqueous solutions containing 0.5 m Ln(ClO4)3 in 1 m HClO4. Similar experiments on background aqueous solutions containing 2.5 molal HClO 4. Experiments were conducted in transmission geometry at the 11ID-beamline of the Advanced Photon Source (Argonne National Laboratory) with an incident photon energy of 91 keV, corresponding to a wavelength of 0.13702 Å, well above the lanthanide K-absorption edges. Scattering intensities were measured with an amorphous silicon flat panel X-ray detector mounted in a static position (2 = 0°), providing detection in momentum-transfer space Q up to 32 Å-1. The wide Q range afforded by the high-energy photons is important for two reasons, (i) normalizing intensities of background and sample scattering functions, S(Q), for accurate S(Q) determination; and (ii) adequate correlation peak resolutions in g(r) after Fourier transform.67-69 After scattering intensities are integrated and corrected for polarization, Compton scattering, and incomplete detector absorption, the resulting background S(Q)s are subtracted from the corresponding sample data to yield an SQ), which contains only Ln correlation information. These and other details of sample preparation and the data reduction method used to obtain the g (r) are described elsewhere. 28, 68

Calculation of HEXS dPDFs from CHARMM and AMOEBA Simulations To calculate a PDF from the CHARMM simulations for comparison with the HEXS dPDF, we simulated two additional systems with 0.496 m Ln (ClO4)3 (Ln = La, Lu) and 0.99 m HClO4, 8 ACS Paragon Plus Environment

Page 8 of 27

Page 9 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


compared to the 0.5 m and 2.5 m concentrations of Ln and ClO4- in order to mimic concentrations used experimentally. In these simulations, 20 La3+ (or Lu3+), 40 H3O+ and 100 ClO4- ions are dissolved in 2240 water molecules. CHARMM force field parameters were chosen for H3O+ 70 and the ClO4- ion71, the latter using the LJ parameters of chloride ion (atom type CLA) and nitrobenzene oxygen (atom type OG2N1) from the CHARMM 36 potential. The production simulations were performed for a duration of 40 ns with a saving frequency of every 10 ps for the data analysis. The pair distribution functions (g(r)) are thus calculated between the Ln3+ and the other elements. The obtained results for the La 3+ system are presented in Figure S2.a. For comparison with the calculated correlations shown in Figure S2a it is necessary to convert the experimental scattering data to difference pair-distribution functions (dPDFs), which only include correlations involving the Ln ion (Ln-O, Ln-H, Ln-Cl and Ln-Ln).68 The pair distribution functions obtained from simulations were treated in the same way as the experiments data in order to compare them to dPDFs, as shown in Figure S2b. Partial scattering functions S(Q) were obtained by Fourier transforming calculated individual g(r) from r space to Q space using 𝑆 (𝑄 ) = 1 +

4𝜋𝜌0 𝑄

∫ 𝑟(𝑔(𝑟) − 1)sin(𝑄𝑟)d𝑟

(9)

where ρ0 is atomic density. The CHARMM simulations include Ln-ClO4 and Ln-Ln correlations whereas the AMOEBA simulations do not. The total scattering function S(Q) for Ln correlations in aqueous solution was calculated by using weighing factors, which depend on relative atomic concentrations and X-ray form factors, as in the experiment. Finally, S(Q)s were back transformed to real space using Qmax = 20 Å-1 for both CHARMM and AMOEBA models for direct comparison with experimental data.

RESULTS AND DISCUSSION Force Field Definition for the Ln3+ Ions.



Figure 1. Calculated properties as a function of LJ parameters (Rmin/2, ε) of a model trivalent metal ion (M3+): (a) the absolute hydration free energy ∆𝐺 ℎ𝑦𝑑; (b) the optimal distance between ion and water oxygen r(M3+-O) and (c) the hydration number (Nwater) of metal ion, both in the first hydration shell. The raw data are provided in the Supporting Information

The calculated hydration free energies are collected in Figure 1.a, with the hydration structures (metal-water oxygen distance, and hydration number) in panels (b) and (c), respectively. The influences of Rmin/2 and 𝜀 on the hydration free energy and hydration structures are qualitatively similar as the previous works. 6-7


Page 10 of 27

Page 11 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


𝜎 (𝑅𝑚𝑖𝑛 /2, 𝜀) = ∑𝑖 (

𝑅 𝐸𝑥𝑝.−𝑆𝑖𝑚𝑢𝑙.( 𝑚𝑖𝑛,𝜀) 2

𝐸𝑥𝑝.

2

)

(𝑖 = ∆𝐺 ℎ𝑦𝑑 , 𝑟(M 3+ − O), 𝑁𝑤𝑎𝑡𝑒𝑟 )

(10)

The unique definition of the LJ parameters is a well-known problem.4 In the present work, we are employing the values of Rmin/2 = 1.90 Å, 1.65 Å and 1.40 Å for La3+, Eu3+ and Lu3+ reported by Durand et al.6 Subsequently, the least square fit (Eq. 10) was employed to obtain the optimal LJ parameter ε at given Rmin/2 for certain metal ion, which is referred by the smallest value of 𝜎 (𝑅𝑚𝑖𝑛 /2, 𝜀). For the other Ln3+ ions, the Rmin/2 were obtained by the interpolation method based on the ionic radius,11 and rounded to the nearest Rmin/2 in Figure 1. In these calculations, the hydration free energy data compiled by Marcus24 were employed. The hydration structure from the X-ray scattering data of Habenschuss and Spedding were used for Tb3+, Dy3+, Er3+, Tm3+ Lu3+, 25

and La3+, Pr3+, Nd3+, 26 and Sm3+, Eu3+. 27 The obtained LJ parameters are listed in Table 1.

Table 1. Optimal LJ Parameters for the Trivalent Lanthanide Ions. Rmin/2

ε

Rmin/2

ε

(Å)

(kcal/mol)

(Å)

(kcal/mol)

La3+

1.90

0.05

Tb3+

1.60

0.05

Pr3+

1.80

0.05

Dy3+

1.55

0.05

Nd3+

1.75

0.05

Er3+

1.50

0.10

Sm3+

1.675

0.08

Tm3+

1.45

0.07

Eu3+

1.65

0.03

Lu3+

1.40

0.12

Ln3+

Ln3+

Accuracy of the New CHARMM Force Field Parameters.



Page 12 of 27

Figure 2. The new CHARMM force field parameters are compared with AMOEBA13 results in reproducing (a) the experimental hydration free energy,24 (b) optimal Ln3+-water oxygen distance and (c) hydration number in the first hydration shells.25-27 The ratio of the simulation results to the corresponding experimental data are reported, with the value of 1 referring to the exact agreement. The original data are provided in Table 2. Table 2. Comparison between the Simulation and Previously Published Experimental or Estimated Results 3+

Ln La3+

Pr3+ Nd3+ Sm3+ Eu3+

CHARMM a AMOEBA b Published CHARMM a Published. CHARMM a Published. CHARMM a Published CHARMM a AMOEBA b

ΔGhyd kJ/mol -3175 -3163 -3303c -3265 -3413 c -3308 -3448 c -3328 -3493 c -3445 -3343

1st hydration layer r(Ln -O), Å Nwater 2.50 9.2 2.5 9 2.580d 9.13d 2.44 9.0 d 2.539 9.22d 2.41 9.0 2.513d 8.9d 2.40 9.0 d 2.474 8.8d 2.28 8.2 2.4 8.8 3+


2nd hydration layer r(Ln3+-O), Å Nwater 4.71 21 4.69 18 4.65e 18e 4.63 20.0 e 4.63 18e 4.59 19.9 4.60e 18e 4.59 19.4 e 4.56 18e 4.49 19.3 4.64 15

Page 13 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Tb

3+

Dy3+ Er3+ Tm3+ Lu3+

Published CHARMM a Published CHARMM a Published CHARMM a Published CHARMM a Published CHARMM a Published

-3528 c -3435 -3568 c -3482 -3593 c -3452 -3663 c -3533 -3683 c -3523 -3683 c

2.450d 2.29 2.409d 2.27 2.396d 2.30 2.369d 2.24 2.358d 2.25 2.338d

8.34d 8.3 8.18d 8.0 7.93d 8.2 8.19d 8.0 8.12d 8.0 7.97d

4.54e 4.53 4.52e 4.49 4.51e 4.51 4.49e 4.44 4.48e 4.45 4.45e

18e 19.7 16-18e 19.0 14-16e 18.9 14-16e 18.5 14-16e 18.7 14-16e

a) CHARMM simulation results in this work b) Reported values by Marjolin et al. 13 c) hydration free energy by Marcus24 with a further re-calibration using the proton ion hydration free energy of ∆𝐺 ∗ (𝐻 +) = -1112 kJ/mol. d) X-ray scattering experimental data.25-27 e) Values estimated by Smirnov and Trostin. 72

The accuracies of the new CHARMM force field parameters are depicted in Figure 2 and Table 2. The current CHARMM parameters reproduce the experimental hydration free energies and metal correlations for all the Ln3+ ions investigated with deviations (|(simulation/experiment)-1|) generally less than 5%. Also presented in Figure 2 are the AMOEBA results using the reported values by Marjolin et al.13 The current CHARMM simulations are observed to provide similar accuracy in reproducing these experimental data as the AMOEBA results, whose deviations vary from 1.5% to 5.5%. For instance, for the La3+ ion ∆𝐺 ℎ𝑦𝑑 = -3303 kJ/mol,24 r(La3+-O) = 2.58 Å

26

and Nwater = 9.13

26

experimentally, which correspond to -3175 kJ/mol (3.9%), 2.50 Å (-3.1%) and 9.2 (0.8%) from the CHARMM simulation, and -3163 kJ/mol (4.2%), 2.50 Å (-3.1%) and 9.0 (-1.4%) reported for the AMOEBA simulation.13 The CHARMM results compare favorably with existing additive potential calculations,5-9 where either the absolute (or relative6) hydration free energy,5, 7-8 or the metal-oxygen distance5-6,

8-9

, or the hydration number6,

9

were separately reported to be in

agreement with the corresponding experimental data. The improved accuracy of the CHARMM parameters proposed herein is ascribed to the inclusion of the corrections in calculating the absolute free energy for (1) the change in standard state concentrations from the gas phase (1 bar) to solution (1 mol/L), (2) the finite size effect of the



simulation box and (3) the discrete solvent contribution under the 3-dimensional PBC conditions and with PME for the long-range electrostatic interactions. Though the first correction has been broadly employed, the last two were not included in the earlier publications using the PBC simulation box.5, 8 Note that the third correction is not necessary in the condition of spherical water droplets.7, 41, 56 Our test calculations, which ignored the last two corrections, lead to lower accuracy (Supporting Information), demonstrating the critical role of the hydration free energy for a precise understanding of metal ion speciation.73 Coordination Geometry of Ln3+ Ions. In addition to the hydration free energy and metal correlations, the geometric arrangement of waters within the Ln first hydration shell can also play an important role in the Ln-complex energetics.10-11, 74-75 Unfortunately, there are no direct experimental techniques currently available to provide information about coordination geometries in aqueous solution, other than indirect evidence such as optical selection rules relevant centric symmetry or NMR spectral or temporal changes, consequently geometries are often inferred from the single-crystal structures of homoleptic Ln-water structures.76 Complicating the situation is the very fast water-exchange rate, known to be on the order of > 5 ns-1,

77-78

which makes them very difficult to quantify

experimentally. Moreover, unlike crystalline solid-state materials, solution dynamics and available equilibria can allow for the simultaneous presence of multiple coordination numbers and geometries, as broadly reflected in various compilations of experimental stability constants representing thermodynamic equilibria.27,

79-83

These issues complicate the experimental

determination of coordination geometries and point to the advantage of employing complementary simulation approaches. Toward this goal, angular distribution functions (ADF) of water oxygen-Ln3+-water oxygen angles are determined from the trajectory snapshots for the La3+ and Lu3+ ions, and presented in Figure 3(a, d). For the La3+ ion, our calculation from the CHARMM simulation finds approximately nine coordinating waters, and supports two optimal O-La3+-O angles of 70° and 137°. These peak angles agree with our AMOEBA simulations using the reported force field parameters by Marjolin et al.,13 and with the polarizable simulation results by Duvail et al. 11 and Morales et al.75 In Figure 3a (3d), the total number of O-La3+-O (O-Lu3+-O) angles has been rescaled to 36 (28)75 for the


Page 14 of 27

Page 15 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


comparison; and in our calculations only the coordination geometry of dominant nona-hydrated La3+ (octa-hydrated Lu3+) ions are considered.

Figure 3. Angular distribution functions of (a) La3+ and (d) Lu3+ systems using the CHARMM parameters proposed here and the AMOEBA parameters by Marjolin et al. 13 Also provided are the discrete distributions of TTP and GySA for La3+ and SA data for Lu3+. The polarizable simulation results reported by Duvail et al. 11 and Morales et al.75 are included for comparison. Examples of the coordination geometries of (b) TTP and (c) GySA around La3+ (pink bead), and (e) SA around Lu3+ (green bead) are provided. The semi-transparent red and blue lines provide a general summary of the coordination-geometry trends.

Duvail et al.11 argued the presence of the tricapped trigonal prism (TTP, Figure 3b) geometry for La3+ ion by the observation of two ADF peaks at 70° and 137° (Figure 3a). Nevertheless the similar ADF curve was assigned to the gyro-elongated square antiprism (GySA, Figure 3c) by Morales et al.75 As demonstrated in Figure 3a, neither TTP nor GySA geometry alone account for the La3+ speciation in solution. Instead the simulation trajectories are consistent with the simultaneous presence of both nine-coordinate geometries, an uncommon situation in simulation trajectories. More common are deformed coordination geometries that conform to neither TTP nor



GySA. This explains the great deviation of the ADF for La3+ ions from the TTP and GySA discrete distributions. The lack of discrete O-Ln3+-O-angles is ascribed to the absence of lattice energies for molecules in solution11, 75 and is consistent with known solution kinetics.77 In contrast, the Lu3+ solution structures in the CHARMM simulation are in good agreement with the expected eight-coordinate, square-antiprism (SA) geometry (Figure 3d,e),11, 75 similar to that seen in the solid state.76 An inspection of the simulation trajectory supports a stable SA geometry for around one hundred picoseconds. No AMOEBA result is provided for Lu3+ due to the current lack of available force field parameters. A comparison with results from Duvail et al.11 and Morales et al.75 reveal that even though the peaks at around 75° and 142° are observed in all the simulations, the peak at around 119° is better supported by our CHARMM simulation. Moreover, the ratios of the occurrences at the three angles for the SA discrete distribution (16 : 4 : 8 for the O-M-O angles of 74.86° : 118.53° : 141.59°)75 are quantitatively reproduced, with a slight underestimation of the probability at 118.53°. Comparisons of Difference Pair-Distance Functions (dPDF) from HEXS with those derived from Simulations.

Figure 4. dPDFs from HEXS experiments are compared with metal-centric PDFs computed from the CHARMM (La3+, Lu3+) and AMOEBA (La3+) simulations.67-68 The CHARMM patterns include a contribution from the anion, which was absent in the AMOEBA pattern.


Page 16 of 27

Page 17 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In addition to the first hydration shell, atomistic simulations are able to provide longer-distance correlation information, which can be directly compared with HEXS experiments. HEXS provides structural correlations for more distance interactions that are inaccessible to other experimental techniques.84 Very recently, reports have become available of HEXS data coupled with density functional theory85 and all-atom MD simulations86 in studying solution systems. The metal-centric PDF patterns obtained from the simulations are compared with each other and the dPDFs from experiments (Figure 4). In these CHARMM simulations the experimental conditions were simulated whereas in the AMOEBA simulation, only one system with one La 3+ ion dissolved in 1014 waters was investigated due to the lack of available AMOEBA parameters for the Lu3+, ClO4- and H3O+ molecules. Further detail is available in the Supporting Information. The PDF patterns for La3+ ion derived from AMOEBA and CHARMM simulations (Figure 4a) are similar, particularly for the first hydration shell. In both cases an analysis determines approximately nine coordinating waters at a distance of 2.5 Å and similar La-water distance distributions. Nevertheless their comparison with experiment shows a mismatch; the experiment has a significantly non-Gaussian peak shape, with correlations extending out to longer distances. This non-Gaussian feature is consistent with a mixed TTP/GySA coordination geometry (Figure 3a-c), though no conformational evidence is available to date. In contrast the experimental dPDF in the first hydration shell of Lu3+ ions is in good agreement with the CHARMM simulation, where eight coordinate waters are found,25 coordinating with a SA coordination geometry (Figure 3d). Additionally, second-sphere peaks are observed in Figure 4 at 4.7/4.5 Å for La3+/Lu3+, with much weaker but still observable third peaks at about 6.9/6.8 Å for La3+/Lu3+ and fourth peaks at about 9 Å (Figure S2.b). These diffuse peaks explicitly provide evidence for longer-range correlations for higher-valent ions in aqueous solution.84, 87-88 In terms of the second- and more-distant metal correlations, a comparison of the CHARMM simulations with and without the inclusion of 0.496 m Ln(ClO4)3 and 0.99 m HClO4 (inset Figure S2) favor the former to better reproduce the experimental HEXS dPDF. For instance, the experimental structural feature at the distance of around 5.9 Å (5.7 Å for Lu 3+) could only be reproduced by the simulations with the counterions and acid included. The ClO4- counterion is generally assumed to be a weak, non-complexing ligand89 yet here we find that the second 17 ACS Paragon Plus Environment


coordination shell is influenced by the ClO4-, stabilized by the formation of hydrogen bonds with the waters directly coordinating Ln3+ ions. This finding suggests that an accurate understanding of a metal ion’s speciation must take into account the effects of even weakly coordinating counterions such as ClO4-. CONCLUSIONS Within the context of modeling metal-ion speciation in water accurately and computationally efficiently, we developed additive CHARMM parameters for the series of trivalent lanthanide ions. Comparisons of the CHARMM simulations with polarizable AMOEBA calculations demonstrates similar accuracies in reproducing experimental hydration free energies, metal-ligand (water) distances and hydration numbers for the first hydration shells. Further comparison with the HEXS experiments show that the CHARMM simulations (AMOEBA calculations not performed) are capable of capturing the Lu3+-water distance distribution on the basis of the dPDF patterns, but further improvements for simulations of the La3+-water correlations would provide critical insights into the current understanding of the complex coordination environment. Similarly, the square antiprism coordination geometry of Lu3+ ion is well captured by the CHARMM simulation, but details of the La3+ geometry remain unclear. Correlations for the second- and more-distant coordination spheres are better reproduced by the CHARMM simulation in which the anions were included. Whereas these findings demonstrate the need for further increasing modeling capabilities they also highlight the importance of including the anion contribution when modeling metal-ion speciation, an expansion that favors the use of less resource intensive additive potentials. Opportunities for improvement still exist. For instance, an accurate understanding of the hydration behavior beyond a dissolved metal-ion’s first coordination shell is a topic under development. 3, 9093

The further development of experimental techniques and accurate quantum mechanics 94 are

highly desired in this regard. The current CHARMM additive potential might be further improved by, for example, the NBFIX method95 and the Lennard-Jones 12-6-4 potential8 though these approaches will increase the complexity and limit the transferability of the parameters. Further improvement of the AMOEBA polarizable potential might be affected through the new AMOEBA14 water model96 and the many-body solvation effects.97 Following the current status in the field of atomistic MD simulations, the force-field parameters of both additive-potentials (CHARMM, OPLS, AMBER) and polarizable potentials (AMOEBA, 18 ACS Paragon Plus Environment

Page 18 of 27

Page 19 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


DRUDE) are generally assumed to be transferable to other molecules and to more complex conditions. That is, the expectation is that parameters of simple molecules or fragments developed in vacuum or bulk solutions can simply be transferred directly to more complex systems. Contrary to this assumption, some studies have reported that force-field parameters are system-dependent. For instance, it has been reported that a force field proposed for UO22+ in aqueous solution needed to be modified from those optimized in vacuum.97 Notwithstanding such reports, assuming forcefield transference within similar systems continues to provide insights into observable chemical properties of metal ions in aqueous solution. From a broader perspective, 3-dimensional modeling of metal-ion coordination symmetries as described herein provides an opportunity to advance the general understanding of their solubility and reactivity in solution. Current perspectives do not holistically address solute coordination in these terms but instead describe the same fixed coordination environment used to model solidstate structures. Missing are the thermodynamic equilibria and kinetic exchange information that fundamentally define metal-ion solution speciation. For example, the inner-sphere coordination numbers for lanthanide ions in aqueous solutions are described as changing from nine for the lighter (larger) Ln to eight for the heavier (smaller) Ln, even though equilibria are known to exist 27, 69

and NMR experiments clearly demonstrate very fast exchange of inner/outer-sphere waters.

77-78

Our simulations do reveal water exchange in the first and second hydration shells. Their averaged, static description remain relatively invariant over the simulation time and the one-dimensional structures suggested in RDFs reflect only Ln-O distances. However the three-dimensional coordination geometries, determined from the calculations (angular distribution function, Figure 3), provide a more holistic description that is not currently captured by experiment. The geometric analysis presents coordination structures in the first hydration shell revealing interesting differences in the light and heavy extremes of the lanthanide series. The coordination around La 3+ presents a deformed TTP or GySA geometry (or a combination of both) (Figure 3a-c), whereas a simple SA geometry is sufficient to account for the Lu 3+ environment. (Figure 3d-e). Moreover, it is noteworthy that the agreement between the dPDFs of Lu 3+ from the HEXS experiment and the CHARMM simulation is consistent with a single SA coordination structure around Lu3+ ion. Thus, insights are gained from the geometrical description derived from the simulation that helps 19 ACS Paragon Plus Environment


to understand the discrepancy between the experimentally and computationally derived La 3+ dPDFs. This result points to an avenue for a more insightful description of the La3+ coordination geometry, with a more statistical structural description as a first step. Such an approach would acknowledge the potential for thermodynamic equilibria and capture the dynamics accessible in the calculation. However, the model as presented does not formally recognize the coordination kinetics. What is needed moving into the future is a new description of metal-ion speciation that incorporates the statistical aspect of solution behavior within a less rigid formalism.

ASSOCIATED CONTENT Supporting Information This material is available free of charge via the Internet at http://pubs.acs.org. Review of the existing force fields for Ln3+ ions; benchmark of CHARMM and AMOEBA simulation speeds. (PDF) Raw simulation data on FEP calculations (XLSX) AUTHOR INFORMATION Corresponding Authors * Email: (B.Q.) [email protected]; (L.S.) [email protected] ORCID B.Q.: 0000-0001-8870-5985; L.S.: 0000-0003-4435-2721 Notes The authors declare no competing financial interests.

ACKNOWLEDGMENT This work, supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Biosciences and Geosciences, under Contract DEAC02-06CH11357, used the Advanced Photon Source, a U.S. Department of Energy (DOE) 20 ACS Paragon Plus Environment

Page 20 of 27

Page 21 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Office of Science User Facility operated for the DOE Office of Science by Argonne National Laboratory under the same contract number. The computing resources, provided on Blues, a highperformance computing cluster operated by the Laboratory Computing Resource Center at Argonne National Laboratory, are gratefully acknowledged.

REFERENCES: 1. Ohtaki, H.; Radnai, T., Structure and Dynamics of Hydrated Ions. Chem. Rev. 1993, 93 (3), 11571204. 2. Marcus, Y., Effect of Ions on the Structure of Water: Structure Making and Breaking. Chem. Rev. 2009, 109 (3), 1346-1370. 3. Gutowski, K. E.; Dixon, D. A., Predicting the Energy of the Water Exchange Reaction and Free Energy of Solvation for the Uranyl Ion in Aqueous Solution. J. Phys. Chem. A 2006, 110 (28), 8840-8856. 4. Li, P.; Merz, K. M., Metal Ion Modeling Using Classical Mechanics. Chem. Rev. 2017, 117 (3), 1564-1686. 5. van Veggel, F. C. J. M.; Reinhoudt, D. N., New, Accurate Lennard-jones Parameters for Trivalent Lanthanide Ions, Tested on [18]crown-6. Chem. Eur. J. 1999, 5 (1), 90-95. 6. Durand, S.; Dognon, J.-P.; Guilbaud, P.; Rabbe, C.; Wipff, G., Lanthanide and Alkaline-earth Complexes of EDTA in Water: A Molecular Dynamics Study of Structures and Binding Selectivities. J. Chem. Soc. Perkin Trans. 2 2000, (4), 705-714. 7. Won, Y., Force Field for Monovalent, Divalent, and Trivalent Cations Developed under the Solvent Boundary Potential. J. Phys. Chem. A 2012, 116 (47), 11763-11767. 8. Li, P.; Song, L. F.; Merz, K. M., Parameterization of Highly Charged Metal Ions Using the 12-6-4 LJ-Type Nonbonded Model in Explicit Water. J. Phys. Chem. B 2015, 119 (3), 883-895. 9. Migliorati, V.; Serva, A.; Terenzio, F. M.; D’Angelo, P., Development of Lennard-Jones and Buckingham Potentials for Lanthanoid Ions in Water. Inorg. Chem. 2017, 56 (11), 6214-6224. 10. Duvail, M.; Souaille, M.; Spezia, R.; Cartailler, T.; Vitorge, P., Pair Interaction Potentials with Explicit Polarization for Molecular Dynamics Simulations of La3+ in Bulk Water. J. Chem. Phys. 2007, 127 (3), 034503. 11. Duvail, M.; Vitorge, P.; Spezia, R., Building a Polarizable Pair Interaction Potential for Lanthanoids(III) in Liquid Water: A molecular Dynamics Study of Structure and Dynamics of the Whole Series. J. Chem. Phys. 2009, 130 (10), 104501. 12. Martelli, F.; Abadie, S.; Simonin, J. P.; Vuilleumier, R.; Spezia, R., Lanthanoids(III) and Actinoids(III) in Water: Diffusion Coefficients and Hydration Enthalpies From Polarizable Molecular Dynamics Simulations. Pure Appl. Chem. 2013, 85 (1), 237-246. 13. Marjolin, A.; Gourlaouen, C.; Clavaguéra, C.; Ren, P. Y.; Piquemal, J.-P.; Dognon, J.-P., Hydration Gibbs Free Energies of Open and Closed Shell Trivalent Lanthanide and Actinide Cations From Polarizable Molecular Dynamics. J. Mol. Model. 2014, 20 (10), 2471. 14. Duvail, M.; Ruas, A.; Venault, L.; Moisy, P.; Guilbaud, P., Molecular Dynamics Studies of Concentrated Binary Aqueous Solutions of Lanthanide Salts: Structures and Exchange Dynamics. Inorg. Chem. 2010, 49 (2), 519-530. 15. Duvail, M.; Guilbaud, P., Understanding the Nitrate Coordination to Eu3+ Ions in Solution by Potential of Mean Force Calculations. Phys. Chem. Chem. Phys. 2011, 13 (13), 5840-5847.



16. Ferru, G.; Reinhart, B.; Bera, M. K.; Olvera de la Cruz, M.; Qiao, B.; Ellis, R. J., The Lanthanide Contraction beyond Coordination Chemistry. Chem. Eur. J. 2016, 22 (20), 6899-6904. 17. Spezia, R.; Migliorati, V.; D’Angelo, P., On the development of polarizable and Lennard-Jones force fields to study hydration structure and dynamics of actinide(III) ions based on effective ionic radii. J. Chem. Phys. 2017, 147 (16), 161707. 18. Ponder, J. W.; Wu, C.; Ren, P.; Pande, V. S.; Chodera, J. D.; Schnieders, M. J.; Haque, I.; Mobley, D. L.; Lambrecht, D. S.; DiStasio, R. A.; Head-Gordon, M.; Clark, G. N. I.; Johnson, M. E.; Head-Gordon, T., Current Status of the AMOEBA Polarizable Force Field. J. Phys. Chem. B 2010, 114 (8), 2549-2564. 19. Wilson, A. M.; Bailey, P. J.; Tasker, P. A.; Turkington, J. R.; Grant, R. A.; Love, J. B., Solvent extraction: the coordination chemistry behind extractive metallurgy. Chem. Soc. Rev. 2014, 43, 123-134. 20. Wang, L.-P.; Head-Gordon, T.; Ponder, J. W.; Ren, P.; Chodera, J. D.; Eastman, P. K.; Martinez, T. J.; Pande, V. S., Systematic Improvement of a Classical Molecular Model of Water. J. Phys. Chem. B 2013, 117 (34), 9956-9972. 21. Wu, J. C.; Chattree, G.; Ren, P., Automation of AMOEBA Polarizable Force Field Parameterization for Small Molecules. Theor. Chem. Acc. 2012, 131 (3), 1138. 22. Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; Mackerell, A. D., CHARMM General Force Field: A Force Field for Drug-Like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields. J. Comput. Chem. 2010, 31 (4), 671-690. 23. Best, R. B.; Zhu, X.; Shim, J.; Lopes, P. E. M.; Mittal, J.; Feig, M.; MacKerell, A. D., Optimization of the Additive CHARMM All-Atom Protein Force Field Targeting Improved Sampling of the Backbone ϕ, ψ and Side-Chain χ1 and χ2 Dihedral Angles. J. Chem. Theory Comput. 2012, 8 (9), 3257-3273. 24. Marcus, Y., A Simple Empirical Model Describing the Thermodynamics of Hydration of Ions of Widely Varying Charges, Sizes, and Shapes. Biophys. Chem. 1994, 51 (2), 111-127. 25. Habenschuss, A.; Spedding, F. H., The Coordination (Hydration) of Rare Earth ions in Aqueous Chloride Solutions from X-Ray Diffraction. I.TbCl3,DyCl3, ErCl3,TmCl3,and LuCl3. J. Chem. Phys. 1979, 70 (6), 2797-2806. 26. Habenschuss, A.; Spedding, F. H., The Coordination (Hydration) of Rare Earth ions in Aqueous Chloride Solutions from X-Ray Diffraction. II. LaCl3, PrCl3, and NdCl3. J. Chem. Phys. 1979, 70 (8), 37583763. 27. Habenschuss, A.; Spedding, F. H., The Coordination (Hydration) of Rare Earth ions in Aqueous Chloride Solutions from X-Ray Diffraction. III. Samarium(III) chloride, europium(III) chloride, and series behavior. J. Chem. Phys. 1980, 73 (1), 442-50. 28. Soderholm, L.; Skanthakumar, S.; Neuefeind, J., Determination of Actinide Speciation in Solution Using High-Energy X-Ray Scattering. Anal. Bioanal. Chem. 2005, 383 (1), 48-55. 29. Skanthakumar, S.; Antonio, M. R.; Wilson, R. E.; Soderholm, L., The Curium Aqua Ion. Inorg. Chem. 2007, 46 (9), 3485-3491. 30. Wilson, R. E.; Skanthakumar, S.; Burns, P. C.; Soderholm, L., Structure of the Homoleptic Thorium(IV) Aqua Ion [Th(H2O)10]Br4. Angew. Chem. Int. Ed. 2007, 46 (42), 8043-8045. 31. Skanthakumar, S.; Antonio, M. R.; Soderholm, L., A Comparison of Neptunyl(V) and Neptunyl(VI) Solution Coordination: The Stability of Cation-Cation Interactions. Inorg. Chem. 2008, 47 (11), 45914595. 32. Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K., Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26 (16), 17811802. 33. Neria, E.; Fischer, S.; Karplus, M., Simulation of activation free energies in molecular systems. J. Chem. Phys. 1996, 105 (5), 1902-1921.


Page 22 of 27

Page 23 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


34. MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M., All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102 (18), 3586-3616. 35. Mark, P.; Nilsson, L., Structure and Dynamics of the TIP3P, SPC, and SPC/E Water Models at 298 K. J. Phys. Chem. A 2001, 105 (43), 9954-9960. 36. Miyamoto, S.; Kollman, P. A., SETTLE: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. J. Comput. Chem. 1992, 13, 952-962. 37. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L., Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79 (2), 926-935. 38. Darden, T.; York, D.; Pedersen, L., Particle Mesh Ewald: An N·log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089-92. 39. Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L., A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577-93. 40. Cisneros, G. A.; Karttunen, M.; Ren, P.; Sagui, C., Classical Electrostatics for Biomolecular Simulations. Chem. Rev. 2014, 114 (1), 779-814. 41. Lu, X.; Cui, Q., Charging Free Energy Calculations Using the Generalized Solvent Boundary Potential (GSBP) and Periodic Boundary Condition: A Comparative Analysis Using Ion Solvation and Oxidation Free Energy in Proteins. J. Phys. Chem. B 2013, 117 (7), 2005-2018. 42. Bishop, T. C.; Skeel, R. D.; Schulten, K., Difficulties with Multiple Time Stepping and Fast Multipole Algorithm in Molecular Dynamics. J. Comput. Chem. 1997, 18, 1786-1790. 43. Martyna, G. J.; Tobias, D. J.; Klein, M. L., Constant Pressure Molecular Dynamics Algorithms. J. Chem. Phys. 1994, 101 (5), 4177-4189. 44. Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R., Constant Pressure Molecular Dynamics Simulation: The Langevin Piston Method. J. Chem. Phys. 1995, 103 (11), 4613-4621. 45. Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M., PACKMOL: A Package for Building Initial Configurations for Molecular Dynamics Simulations. J. Comput. Chem. 2009, 30, 2157-2164. 46. Zwanzig, R. W., High‐Temperature Equation of State by a Perturbation Method. I. Nonpolar Gases. J. Chem. Phys. 1954, 22 (8), 1420-1426. 47. Lamoureux, G.; Roux, B., Absolute Hydration Free Energy Scale for Alkali and Halide Ions Established from Simulations with a Polarizable Force Field. J. Phys. Chem. B 2006, 110 (7), 3308-3322. 48. Yu, H.; Whitfield, T. W.; Harder, E.; Lamoureux, G.; Vorobyov, I.; Anisimov, V. M.; MacKerell, A. D.; Roux, B., Simulating Monovalent and Divalent Ions in Aqueous Solution Using a Drude Polarizable Force Field. J. Chem. Theory Comput. 2010, 6 (3), 774-786. 49. Lin, Y.-L.; Aleksandrov, A.; Simonson, T.; Roux, B., An Overview of Electrostatic Free Energy Computations for Solutions and Proteins. J. Chem. Theory Comput. 2014, 10 (7), 2690-2709. 50. Hummer, G.; Pratt, L. R.; García, A. E., Free Energy of Ionic Hydration. J. Phys. Chem. 1996, 100 (4), 1206-1215. 51. Kelly, C. P.; Cramer, C. J.; Truhlar, D. G., Single-Ion Solvation Free Energies and the Normal Hydrogen Electrode Potential in Methanol, Acetonitrile, and Dimethyl Sulfoxide. J. Phys. Chem. B 2007, 111 (2), 408-422. 52. Rocklin, G. J.; Mobley, D. L.; Dill, K. A.; Hünenberger, P. H., Calculating the binding free energies of charged species based on explicit-solvent simulations employing lattice-sum methods: An accurate correction scheme for electrostatic finite-size effects. J. Chem. Phys. 2013, 139 (18), 184103.



53. Kastenholz, M. A.; Hünenberger, P. H., Computation of methodology-independent ionic solvation free energies from molecular simulations. II. The hydration free energy of the sodium cation. J. Chem. Phys. 2006, 124 (22), 224501. 54. Zhang, H.; Jiang, Y.; Yan, H.; Yin, C.; Tan, T.; van der Spoel, D., Free-Energy Calculations of Ionic Hydration Consistent with the Experimental Hydration Free Energy of the Proton. J. Phys. Chem. Lett. 2017, 8 (12), 2705-2712. 55. Misin, M.; Fedorov, M. V.; Palmer, D. S., Hydration Free Energies of Molecular Ions from Theory and Simulation. J. Phys. Chem. B 2016, 120 (5), 975-983. 56. Simonson, T.; Hummer, G.; Roux, B., Equivalence of M- and P-Summation in Calculations of Ionic Solvation Free Energies. J. Phys. Chem. A 2017, 121 (7), 1525-1530. 57. Marcus, Y., Thermodynamics of solvation of ions. Part 5.-Gibbs free energy of hydration at 298.15 K. J. Chem. Soc. Faraday Trans 1991, 87 (18), 2995-2999. 58. Zhan, C.-G.; Dixon, D. A., Absolute Hydration Free Energy of the Proton from First-Principles Electronic Structure Calculations. J. Phys. Chem. A 2001, 105, 11534-11540. 59. Tissandier, M. D.; Cowen, K. A.; Feng, W. Y.; Gundlach, E.; Cohen, M. H.; Earhart, A. D.; Coe, J. V.; Tuttle, T. R., The Proton's Absolute Aqueous Enthalpy and Gibbs Free Energy of Solvation from ClusterIon Solvation Data. J. Phys. Chem. A 1998, 102 (40), 7787-7794. 60. Camaioni, D. M.; Schwerdtfeger, C. A., Comment on “Accurate Experimental Values for the Free Energies of Hydration of H+, OH-, and H3O+”. J. Phys. Chem. A 2005, 109 (47), 10795-10797. 61. Kelly, C. P.; Cramer, C. J.; Truhlar, D. G., Aqueous Solvation Free Energies of Ions and Ion−Water Clusters Based on an Accurate Value for the Absolute Aqueous Solvation Free Energy of the Proton. J. Phys. Chem. B 2006, 110 (32), 16066-16081. 62. Ren, P.; Ponder, J. W., Polarizable Atomic Multipole Water Model for Molecular Mechanics Simulation. J. Phys. Chem. B 2003, 107 (24), 5933-5947. 63. Shi, Y.; Xia, Z.; Zhang, J.; Best, R.; Wu, C.; Ponder, J. W.; Ren, P., Polarizable Atomic MultipoleBased AMOEBA Force Field for Proteins. J. Chem. Theory Comput. 2013, 9 (9), 4046-4063. 64. Eastman, P.; Friedrichs, M. S.; Chodera, J. D.; Radmer, R. J.; Bruns, C. M.; Ku, J. P.; Beauchamp, K. A.; Lane, T. J.; Wang, L.-P.; Shukla, D.; Tye, T.; Houston, M.; Stich, T.; Klein, C.; Shirts, M. R.; Pande, V. S., OpenMM 4: A Reusable, Extensible, Hardware Independent Library for High Performance Molecular Simulation. J. Chem. Theory Comput. 2013, 9 (1), 461-469. 65. Chow, K.-H.; Ferguson, D. M., Isothermal-Isobaric Molecular Dynamics Simulations with Monte Carlo Volume Sampling. Comput. Phys. Commun. 1995, 91 (1), 283-289. 66. Åqvist, J.; Wennerström, P.; Nervall, M.; Bjelic, S.; Brandsdal, B. O., Molecular Dynamics Simulations of Water and Biomolecules with a Monte Carlo Constant Pressure Algorithm. Chem. Phys. Lett. 2004, 384 (4–6), 288-294. 67. Egami, T.; Billinge, S. I. L., Underneath the Bragg Peaks : Structural Analysis of Complex Materials, 2nd Edition. Pergamon: 2012; p 422 pp. 68. Skanthakumar, S.; Soderholm, L., Studying actinide correlations in solution using high energy xray scattering. MRS Symp. Proc. 2006, 893, 411-416. 69. Soderholm, L.; Skanthakumar, S.; Wilson, R. E., Structures and energetics of erbium chloride complexes in aqueous solution. J. Phys. Chem. A 2009, 113, 6391-6397. 70. van Keulen, S. C.; Gianti, E.; Carnevale, V.; Klein, M. L.; Rothlisberger, U.; Delemotte, L., Does Proton Conduction in the Voltage-Gated H+ Channel hHv1 Involve Grotthuss-Like Hopping via Acidic Residues? J. Phys. Chem. B 2017, 121 (15), 3340-3351. 71. Cadena, C.; Maginn, E. J., Molecular Simulation Study of Some Thermophysical and Transport Properties of Triazolium-Based Ionic Liquids. J. Phys. Chem. B 2006, 110 (36), 18026-18039.


Page 24 of 27

Page 25 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


72. Smirnov, P. R.; Trostin, V. N., Sructural Parameters of the Nearest Surrounding of Tri- and Tetravalent Actinide Ions in Aqueous Solutions of Actinide Salts. Russ. J. Gen. Chem. 2012, 82 (7), 12041213. 73. Skanthakumar, S.; Jin, G. B.; Lin, J.; Vallet, V.; Soderholm, L., Linking Solution Structures and Energetics: Thorium Nitrate Complexes. J. Phys. Chem. B 2017, 121 (36), 8577-8584. 74. Lindqvist-Reis, P.; Apostolidis, C.; Rebizant, J.; Morgenstern, A.; Klenze, R.; Walter, O.; Fanghänel, T.; Haire, R. G., The Structures and Optical Spectra of Hydrated Transplutonium Ions in the Solid State and in Solution. Angew. Chem. Int. Ed. 2007, 46 (6), 919-922. 75. Morales, N.; Galbis, E.; Martínez, J. M.; Pappalardo, R. R.; Sánchez Marcos, E., Identifying Coordination Geometries of Metal Aquaions in Water: Application to the Case of Lanthanoid and Actinoid Hydrates. J. Phys. Chem. Lett. 2016, 7 (21), 4275-4280. 76. Harrowfield, J. M.; Kepert, D. L.; Patrick, J. M.; White, A. H., Structure and Steriochemistry in 'fBlock' Complexes of High Coordination Number. VIII The [M(unidentate) 9] System. Crystal Structures of [M(OH2)9][CF3SO3]3, M = La, Gd, Lu, Y. Aust. J. Chem. 1983, 36, 483-92. 77. Helm, L.; Merbach, A. E., Water exchange on metal ions: experiments and simulations. Coord. Chem. Rev. 1999, 187, 151-181. 78. Helm, L.; Berbach, A. E., Inorganic and Bioinorganic Solvent Exchange Mechanisms. Chem. Rev. 2005, 105, 1923-1959. 79. Waite, T. D., Thermodynamics of the Iron System in Seawater. IUPAC Ser. Anal. Phys. Chem. Environ. Syst. 2001, 7 (Biogeochemistry of Iron in Seawater), 291-342. 80. Grenthe, I.; Fuger, J.; Konings, R. J. M.; Lemire, R. J.; Muller, A. B.; Nguyen-trung, C.; Wanner, H., Chemical Thermodynamics of Uranium. Elsevier: Amsterdam, 1992; Vol. 1. 81. Neuefeind, J.; Soderholm, L.; Skanthakumar, S., Experimental Coordination Environment of Uranyl(VI) in Aqueous Solution. J. Phys. Chem. A 2004, 108 (14), 2733-2739. 82. Kremer, C.; Torres, J.; Dominguez, S., Lanthanide Complexes with Oda, Ida, and Nta: From Discrete Coordination Compounds to Supramolecular Assemblies. J. Mol. Struct. 2008, 879 (1-3), 130149. 83. Gianguzza, A.; Giuffre, O.; Piazzese, D.; Sammartano, S., Aqueous Solution Chemistry of Alkyltin(IV) Compounds for Speciation Studies in Biological Fluids and Natural Waters. Coord. Chem. Rev. 2012, 256 (1-2), 222-239. 84. Kalaji, A.; Skanthakumar, S.; Kanatzidis, M. G.; Mitchell, J. F.; Soderholm, L., Changing Hafnium Speciation in Aqueous Sulfate Solutions: A High-Energy X-ray Scattering Study. Inorg. Chem. 2014, 53 (12), 6321-6328. 85. Yang, K. R.; Matula, A. J.; Kwon, G.; Hong, J.; Sheehan, S. W.; Thomsen, J. M.; Brudvig, G. W.; Crabtree, R. H.; Tiede, D. M.; Chen, L. X.; Batista, V. S., Solution Structures of Highly Active Molecular Ir Water-Oxidation Catalysts from Density Functional Theory Combined with High-Energy X-ray Scattering and EXAFS Spectroscopy. J. Am. Chem. Soc. 2016, 138 (17), 5511-5514. 86. Fujii, K.; Ueki, T.; Hashimoto, K.; Kobayashi, Y.; Kitazawa, Y.; Hirosawa, K.; Matsugami, M.; Ohara, K.; Watanabe, M.; Shibayama, M., Microscopic Structure of Solvated Poly(benzyl methacrylate) in an Imidazolium-Based Ionic Liquid: High-Energy X-ray Total Scattering and All-Atom MD Simulation Study. Macromolecules 2017, 50 (12), 4780-4786. 87. Soderholm, L.; Skanthakumar, S.; Gorman-Lewis, D.; Jensen, M. P.; Nagy, K. L., Characterizing solution and solid-phase amorphous uranyl silicates. Geochim. Cosmochim. Acta 2008, 72 (1), 140-150. 88. Hu, Y.-J.; Knope, K. E.; Skanthakumar, S.; Kanatzidis, M. G.; Mitchell, J. F.; Soderholm, L., Understanding the Role of Aqueous Solution Speciation and Its Application to the Directed Syntheses of Complex Oxidic Zr Chlorides and Sulfates. J. Am. Chem. Soc. 2013, 135 (38), 14240-14248. 89. Rosenthal, M. R., The Myth of the Non-Coordinating Anion. J. Chem. Educ. 1973, 50 (5), 331.



90. Schwenk, C. F.; Rode, B. M., Extended ab initio Quantum Mechanical/Molecular Mechanical Molecular Dynamics Simulations of Hydrated Cu2+. J. Chem. Phys. 2003, 119 (18), 9523-9531. 91. Bryantsev, V. S.; Diallo, M. S.; Goddard, W. A., Computational Study of Copper(II) Complexation and Hydrolysis in Aqueous Solutions Using Mixed Cluster/Continuum Models. J. Phys. Chem. A 2009, 113 (34), 9559-9567. 92. Kuta, J.; Clark, A. E., Trends in Aqueous Hydration Across the 4f Period Assessed by Reliable Computational Methods. Inorg. Chem. 2010, 49 (17), 7808-7817. 93. Parmar, P.; Samuels, A.; Clark, A. E., Applications of Polarizable Continuum Models To Determine Accurate Solution-Phase Thermochemical Values Across a Broad Range of Cation Charge - The Case of U(III-VI). J. Chem. Theory Comput. 2015, 11 (1), 55-63. 94. Dolg, M., Computational Methods in Lanthanide and Actinide Chemistry. John Wiley & Sons Ltd: Chichester, UK, 2015. 95. Luo, Y.; Roux, B., Simulation of Osmotic Pressure in Concentrated Aqueous Salt Solutions. J. Phys. Chem. Lett. 2010, 1 (1), 183-189. 96. Laury, M. L.; Wang, L.-P.; Pande, V. S.; Head-Gordon, T.; Ponder, J. W., Revised Parameters for the AMOEBA Polarizable Atomic Multipole Water Model. J. Phys. Chem. B 2015, 119 (29), 9423-9437. 97. Rai, N.; Tiwari, S. P.; Maginn, E. J., Force Field Development for Actinyl Ions via Quantum Mechanical Calculations: An Approach to Account for Many Body Solvation Effects. J. Phys. Chem. B 2012, 116 (35), 10885-10897.


Page 26 of 27

Page 27 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For Table of Contents Only

Comparative CHARMM and AMOEBA Simulations of Lanthanide Hydration Energetics and Experimental Aqueous-Solution Structures Baofu Qiao,* S. Skanthakumar, L.Soderholm* Chemical Sciences and Engineering Division, Argonne National Laboratory, Argonne, Illinois 60439, USA.


Comparative CHARMM and AMOEBA Simulations of Lanthanide

Recommend Documents