Determination of New Cu+, Cu2+, and Zn2+ Lennard-Jones Ion

Aug 14, 2013 - Department of Chemical Engineering, Igualada School of Engineering, Universitat Politècnica de Catalunya, Pça Rei 15, Igualada. 08700...
0 downloads 0 Views 925KB Size
Article pubs.acs.org/JPCB

Determination of New Cu+, Cu2+, and Zn2+ Lennard-Jones Ion Parameters in Acetonitrile Juan Torras*,† and Carlos Alemán‡,§ †

Department of Chemical Engineering, Igualada School of Engineering, Universitat Politècnica de Catalunya, Pça Rei 15, Igualada 08700, Spain ‡ Department of Chemical Engineering, Barcelona School of Industrial Engineering, Universitat Politècnica de Catalunya, Diagonal 647, Barcelona E-08028, Spain § Centre for Research in Nano-Engineering, Universitat Politècnica de Catalunya, Campus Sud, Edifici C′, C/Pasqual i Vila s/n, Barcelona E-08028, Spain ABSTRACT: We present new Lennard-Jones (LJ) parameters for Cu+, Cu2+, and Zn2+ ionacetonitrile interactions. The adjustment of ion parameters is made to reproduce simultaneously experimental solvation free energy and structural properties, namely ion−N distance and coordination numbers. Initially, the methodology has been validated deriving parameters for wellstudied Na+ and Cl− ions in acetonitrile being compared with experimental and theoretical data. The transferability of parameters is checked by the calculation of thermodynamic and structural properties with three different acetonitrile models. The results obtained for transition metal ions show an overall agreement with reference values. The solvation free energy calculated with new LJ trained parameters using a six-site acetonitrile model, and two older three- and six-site acetonitrile models presents, respectively, percent differences of 0.4, 4.8, and 7.3% when compared with experimental values.



INTRODUCTION Divalent ionic metals in aqueous media have an important role in medicine and biology, being essential Zn2+, Mg2+, and Mn2+ ions in many enzymatic processes and Cu2+ and Fe2+ ions mainly due to their ease of oxidation−reduction. This has led to the publication of a considerable amount of work on the parametrization of the long-range nonbonded interactions among different water models of both monovalent1 and divalent2 metals. Moreover, with regard to the abundance of experimental data for nonaqueous solutions leading to important applications in complex formation,3 such as chemical separation and ion recognition, there are not many studies and parametrizations of divalent ions in nonaqueous environments. As a result, an extension of theoretical insight into phenomena of aprotic solvation may be of interest to the condensed phase study. Among all aprotic solvents, acetonitrile (AN) is being used because of their good solvating properties for many hydrophobic and hydrophilic species.4 This aprotic solvent presents an important solvation factor with respect to monovalent metal ions with d10 configuration due to its nature of nitrogen donor solvent, but it also solvates the alkali, alkaline-earth and zinc group metal ions, although more weakly than water.5 When this solvent is compared with water, AN is much less structured; however, there is a strong short-range order with an almost full occupancy of nearest-neighbor sites reported experimentally. They have been found weakly associated dimers and trimers but with little structure beyond the first few neighbors.6 © 2013 American Chemical Society

The reported values of thermodynamic functions for the transfer reaction of free ions from water to AN show that many ions are better solvated in water than in AN. However, while water is preferred by alkali and halide ions, Cu+, Ag+, and Au+ metal ions appear to prefer AN. In fact, the specific solvation of univalent d10 acceptors such as Cu+ is especially pronounced in this solvent.5,7 Solvated structures of divalent metals, such as Cu2+ and Zn2+ in AN, were reported with a metal coordination of 6 (rCu−N(eq) = 1.99 Å, rCu−N(ax) = 2.18−2.23 Å)8 and 6 (rZn−N = 2.11 Å),9 respectively. Moreover, it is shown that Cu+ has a tetrahedral solvation structure (rCu−N = 1.99 Å),10 whereas Cu2+ presents a distorted octahedral configuration.8 Some authors have studied the solvation of free metal ions such as Cu+,11 Cu2+,2 and Zn2+ in water2,12 and being proposed some Lennard-Jones (LJ) parameters. However, there is a lack of LJ parameters for metallic ions in nonaqueous solvent; many of the theoretical studies involving metal ions in aprotic solvent such as AN are using a set of LJ parameters being obtained in order to reproduce ion−water van der Waals (vdW) interactions.13 However, still a few parametrizations can be found related to monovalent ions, especially alkali metals and halogens in AN. Specifically, Richardi et al.14 studied the solvation of some monovalent ions by using two previous parametrizations of Received: March 13, 2013 Revised: July 17, 2013 Published: August 14, 2013 10513

dx.doi.org/10.1021/jp402545g | J. Phys. Chem. B 2013, 117, 10513−10522

The Journal of Physical Chemistry B

Article

Wiese et al.15 and Pálinkás et al.16 that were adjusted to reproduce ion-AN interactions. Subsequently, a work of Guardia et al.17 was devoted to study the solvation structure of Na+ and Cl− ions in AN by applying the Pálinkás parameters. Other later settings were being conducted by Cabaleiro-Lago et al.18 by obtaining new ion−solvent pair potentials derived from computational calculation of metal-solvent clusters, M(AN)n in gas phase with M = Li+, Na+, or Mg2+, using DFT and MP2 methodologies. Also, Fischer et al.19 studied the influence on the solvation structure of ions using two different interactions models for solvent. Thus, a polarizable model is compared with a nonpolarizable model (OPLS) of AN using LJ parameters derived by Pálinkás. Later, Spangberg et al.20 reported a new parametrization of Li+ and Na+ ions in AN using ab initio calculations to obtain the ion−solvent potential. More recently, Hsiao et al.21 conducted an experimental and theoretical study of Cu(AN)n2+ complexes using density functional theory, where both EXAFS and computational methods were combined to interpret Cu2+ ion structures in AN. More recently Seo et al. have conducted studies beyond infinite dilution of Li+ ion by studying dependence of salt aggregation in AN on salt concentration.22 Following the soft-core thermodynamic integration (TI) methodology proposed by Steinbrecher et al.23 we have developed new ion-AN LJ parameters for Cu2+, Zn2+, and Cu+, adjusting them in order to simultaneously reproduce both thermodynamic and structural properties (i.e., distance metalAN and coordination number). To our knowledge, this is the first time that LJ parameters are proposed modeling ion-AN solvation for these three transition metal ions. On the methodology checking, new LJ parameters for Na+ and Cl− were also derived and compared with experimental and theoretical data. The paper is organized as follows. First, theoretical and methodological aspects are presented, such as the solvation energy derivation by TI of the solvation process and its molecular dynamics simulation. Next, structural aspects and molecular dynamic results are presented and discussed. Finally, discussion on the modeling of ion-AN interaction using the new proposed set of LJ parameters is analyzed as well as their transferability to different model solvents.

Table 1 lists the experimental thermodynamic data for some monovalent and divalent ions used in this work whose Table 1. Experimental Data for the Gibbs Free Energy of Hydration ΔGhyd * , Transfer from Water to AN, ΔGtr*, and Solvation in AN, ΔG*solv for the Set of Ions (all values in kcal/ mol) ions

ΔG*hyd Marcus25

ΔG*hyd Kelly27

Cl− Na+ Cu+ Cu2+ Zn2+ H+

−81.3 −87.2 −125.5 −480.4 −467.2 −254.3

−74.5 −103.2 −141.3 −507.0a −490.4 −265.9

ΔG*tr

ΔG*solv Marcus

ΔG*solv Kelly

10.0b 3.5b −12.5b 12.0c 13.2c

−71.3 −83.7 −138.0 −468.4 −454.2

−64.5 −99.7 −153.8 −495.0 −477.2

Taken from ref 25 and adjusted according to the ΔGsolv * (H+) value of −265.9 kcal/mol. bReference 7b. cReference 5.

a

* ) have been tabulated by hydration free energies (ΔGhyd Marcus25 and Kelly et al.27 The former is commonly adopted by the molecular dynamics simulation community and the latter based on the most recent and recommended benchmark value of ΔGhyd * (H+). To obtain ion solvation free energies in AN (ΔGsolv * ) a previously methodology described elsewhere13b,28 was followed. Thus, by combining the data for the aqueous free energies of solvation (ΔG*hyd) for the ions with the data for the free energy of transfer from water to AN (ΔGtr*), we obtained the free energies of solvation for our set of ions in AN, (ΔG*solv).5,7b All experimental data collected by the described procedure are displayed in Table 1. Ion-Acetonitrile Lennard-Jones Parameterization. All MD simulations were performed with AMBER 12 software29 where the electrostatic and van der Waals (vdW) interactions between nonbonded atoms are described within an additive pairwise scheme by a 1-12-6 potential, Uij(rij) =



qiqj rij

⎛ R * ⎞12 ⎛ R * ⎞6 ij ⎟ − 2εij⎜ ij ⎟ + εij⎜⎜ ⎟ ⎜ r ⎟ r ⎝ ij ⎠ ⎝ ij ⎠

(1)

where R*ij and εij are the Lennard-Jones parameters, rij is the distance between atoms i and j, qi and qj are the partial charges on atoms i and j, respectively. Also, The Lorentz−Berthelot mixing rules were used for cross interactions between different ions and between ions and solvents sites: εij = (εii * εjj)1/2 and Rij* = (Rii* + Rjj*)/2 = Ri* + Rj*. Three AN solvent models were used here, two all-atom models and one united atom model, being all of them compatible with the AMBER force field of Cornell et al.30 The most recent six-site AN model was the Nikitin et al. model.31 In that model the point charges were obtained under the RESP approach at MP2 6-311++G(3df,3p) level, with values of −0.5503, 0.4917, −0.5126, and 0.1904 e, for CMe, CN, N, and H atoms, respectively. This solvent model was specifically tuned to fit the experimental density and evaporation heat of AN and is used in this work to derive the new ion-AN LJ parameters. The second all-atom model was the one derived by Grabuleda et al.32 The point charges were obtained under the RESP approach at HF/6-31G* level (−0.237, 0.382, 0.577, −0.490, and 0.115 e, for CMe, CN, N, and H atoms, respectively). Finally, the united atom AN model was derived by Blas et al.33 and Guardia et al.34 using point charges of 0.206, 0.247, and −0.453 e, for CMe, CN and N, respectively. The last two AN

METHODS Solvation Free Energies of Ions. Solvation free energies are thermodynamically defined for neutral salts only. Thus, direct experimental measurements can only provide relative values of ion solvation free energies. Additional extrathermodynamic assumptions are made to separate the solvation free energy of a salt into contributions from the cation and anion in water, which often leads to tabulated values derived with respect to the hydration free energy of the proton ΔG*hyd(H+). However, there is no a single generally accepted * (H+). In fact, current estimations vary by more value for ΔGhyd * (H+). The than 12 kcal/mol among reported values of ΔGhyd estimates of Tissandier et al.24 −265.9 kcal/mol, and Marcus25 −254.3 kcal/mol are widely used. That presents a problem for both ions tabulated and thermodynamic data accuracy which is * (H+). given by those scattered estimated values of ΔGhyd Consequently, there is a significant scatter in the data tabulated according to these values as evolved ΔG*hyd(H+). Most recently Camaioni and Schwerdtfeger26 have recommended using * (H+) = −265.9 kcal/mol, a benchmark value derived by ΔGhyd Tissandier et al.24 and recently reproduced by Kelly et al.27 10514

dx.doi.org/10.1021/jp402545g | J. Phys. Chem. B 2013, 117, 10513−10522

The Journal of Physical Chemistry B

Article

(3)

was equilibrated and was followed by a MD simulation for each TI window. This procedure was repeated at each one of scanned LJ parameters of a given ion to obtain expected solvation free energies. Molecular Dynamics Simulation Protocol. Two different MD simulation protocols were conducted depending on it was either performing a scan of the free energy surface or it was being checked its final ion-AN LJ parameters. Thus, the first one and the quickest protocol involves an ion-AN system made of a solvent box containing 512 AN molecules with the reference ion placed on the center. All simulations were conducted using a 2.0 fs time step while atomic bond distance involving hydrogen atoms were constrained via the SHAKE algorithm.36 The system temperature was constrained to 298 K using the Langevin thermostat, whereas a mean pressure of 1 atm was controlled by the Berendsen barostat.37 The nonbonded cutoff distance was 15 Å and the long-range electrostatics treatment was computed by means of Ewald summations. The default parameters for the tolerance, 10−5, were used on the direct part of the Ewald sum and the Ewald coefficient was automatically determined as a function of tolerance and cutoff values.38 The energy of the system was initially minimized for 2500 steps using steepest descent, followed by 0.2 ns of equilibration using an NPT ensemble to reach constant density, and 0.5 ns of production dynamics at a mean temperature of 298 K. The coordinates of production trajectories were saved every 5000 steps (10 ps intervals) for subsequent analysis. The equilibrated ion-AN system was then used in the solvation free energy derivation. Following with the previous protocol, a MD trajectory of 100 ps for each TI window were obtained. The average derivative of Hamiltonian (⟨dV/dλ⟩) was calculated by using the data collected from each window and integrated numerically by Gaussian quadrature to obtain the free energy difference between the two states of both TI steps; charge neutralization and disappearing. The α parameter, used in this simulation for the modified version of the vdW equation that is used to smoothly switch off nonbonded interactions of the ion in the second step of this TI, is set to 0.5. Once the free energy surface was scanned and final LJ parameters for each ion were obtained, a more extended protocol to derive thermodynamic and structural data was conducted. Thus, following a similar scheme to the one shown above, a bigger ion-AN system containing 1024 AN molecules was used to check new parameters. Also, a longer MD trajectory on the production of the ion-AN system for 2 ns was obtained to derive structural data. Similarly, a longer MD trajectory for each TI window of 0.5 ns long was performed to improve thermodynamic data.

All solvation free energies of ions involved in this work have been calculated in a two-stage TI approach, following the methodology proposed by Steinbrecher et al.23 Thus, in the first stage a charge neutralization is performed by means of a linear perturbation of the ion to neutrality in AN. On the second stage, the neutralized-ion disappears from the AN, being removed the van der Waals potentials using a soft-core vdW potential function. The TI simulations were performed by using AMBER 12 program.29 Then, the integral of eq 2 is evaluated by numerical integration using a twelve-point Gaussian quadrature scheme. The solvation free energy of ions is obtained computing the expected value ⟨dV/dλ⟩λ at each one of the twelve λ-windows being applied on eq 2. Thus, given a reference ion and its tested LJ parameters, the ion-AN system

RESULTS AND DISCUSSION This section shows the results obtained deriving a new set of LJ parameters of ions in AN, and their structural analysis. Initially, one halide and one alkali ion (Cl− and Na+) are chosen to test the general methodology as weakly solvated ions. The results for a short set of transition metal ions (Cu+, Cu2+, and Zn2+) with a stronger interaction with the cyano moiety of solvent molecule are given in a subsequent section and, finally, some transferability studies on the LJ parameters are commented on. Cl− and Na+ Ions in Acetonitrile. Because of the large series of TI runs that are necessary to scan the (R*, ε) free energy surface for a given ion, an initial study in order to reduce and optimize the computational time was conducted. Thus,

models were used to study the transferability of the new LJ parameters to a different AN solvent model. Beside of the different final density of those models (0.774, 0.726, and 0.747 g/cm3, respectively), we want to notice the main difference on the point charge distribution. In order to derive a new set of LJ parameter (R*, ε), the following iterative process was conducted: 1. Initially, for each ion a restricted exploratory space scanning is carried out with the two LJ parameters. Thus, an initial grid of values from 0.3 up to 1.2 Å for R*/2 parameter and from 0.0001 up to 1 kcal/mol for ε parameter are tested. 2. For each one of previous scanned points, a new ion-AN system is build, equilibrated, and a MD production trajectory is obtained to derive the following three physical variables: the solvation free energy, the most probable distance between ion and nitrogen atom from AN, and the ionic coordination number. 3. Those observables are used to train a specific neural network for each studied ion trying to reproduce expected experimental values. Then, a new candidate of LJ parameters is derived. 4. The point 2 and 3 are repeated until matching physical variables derived in step 2 with their experimental values. Cl− and Na+ ions have been trained to reproduce solvation free energy only whereas transition metal ions were trained to reproduce solvation structure and thermodynamics. The solvation free energy was obtained by TI. The ion−nitrogen distance was derived from the first maximum on the radial distribution functions (rdf) for nitrogen atom of AN around the ions. Finally, coordination number is obtained by counting the number of AN molecules within the distance given by the first minimum on the rdf of AN center of mass with respect to the ions. Thermodynamic Integration. According to the thermodynamic integration formalism, the free energy difference between two Lambda-coupled states can be computed from * = ΔGsolv

∫0

1

∂V(λ) ∂λ

dλ λ

(2)

where λ is a nonspatial coordinate to couple the potential functions of the initial and final states unperturbated systems, V0 and V1, respectively. There are several ways to couple the two end-point potential functions into the mixed potential, V(λ). Thus, choosing a simple function f(λ) that goes from zero to one, could be used to obtain V(λ), despite having shown to be a difficult task35 V (λ) = f (λ)V1 + [1 − f (λ)]V0



10515

dx.doi.org/10.1021/jp402545g | J. Phys. Chem. B 2013, 117, 10513−10522

The Journal of Physical Chemistry B

Article

Table 2. Lennard-Jones Parameters (R* in Å, ε in kcal/mol) Derived in this Work, and the Solvation Free Energy using 512 and 1024 AN in the Simulation Box (ΔG*solv in kcal/mol), Coordination Number (CN), First Peaks of rdf’s Between Ion−Nitrogen (dion−N in Å), Ion-Cyano Carbon (dion‑CN in Å), Ion-Methyl Carbon (dion‑Me in Å), and Shell Radius Rshell (in Å) in AN ΔG*solv Cl‑ set1a MDb MOZc MDd MCe Na+ set1a MDg MOZc MDd QMh MCe MDi

R*/2

ε

512AN

1024AN

exp.

2.24

0.004

−64.5

−64.2

−64.5

1.40

0.050

−99.7

−99.2

−99.7

CN 8.9 11 9.3 10.7 9.5f 6.1 6.7, 7.0 8.2 6.4 5.5f 6.7 5.6f 5.2−6.5

dion−N 5.32

dion−CN 4.37

dion−Me

Rshell

3.55

5.6

6.1 3.8 2.29 2.7

3.37 3.7

4.79 5.1

4.2

4.4 2.3−2.5 2.6 2.4−2.6

3.5−3.6

a

LJ parameters adjusted based on Kelly hydration free energy of the proton, set1. bMarkovich et al. ref 39. cRichardi et al. ref 14 using the molecular Ornstein−Zernike theory. dGuardia et al. ref 17. eFischer et al. ref 19. fData derived from polarizable solvent model. gTroxler et al. ref 3. hCabaleiro et al. ref 18. iSpångberg et al. ref 20

different sizes of periodic boxes and different cutoffs on the long-range interactions were used to check the solvation free energy. From TI calculations performed on the Na+ ion we * of about 0.2 kcal/mol observed little differences on ΔGsolv between experiences using 13, 15, and 20 Å of cuttof distances. Also, deriving ΔG*solv by using periodic box sizes containing 512 and 1024 AN molecules did not lead to bigger differences. Hereafter, all scanned points (R*, ε) from the surface of free energy, involving TI to derive new LJ parameters, were conducted using a cuttof distance of 15 Å in a periodic box containing 512 AN molecules. Table 2 shows new LJ parameters derived in this work for the Cl− and Na+ ions in AN that specifically have been fitted to reproduce experimental ΔG*solv obtained from the tabulated ΔG*hyd of Kelly (Table 1), hereafter “set1”. From proposed LJ parameters and using the smaller solvent box, a new equilibration and its subsequent production run is performed during 2 ns in a larger box containing 1024 AN molecules in order to derive ion-AN structural parameters (Table 2). The resulting solvation free energies obtained using smaller solvent boxes are very close to the respective experimental values and do not change significantly when a TI is performed on solvent boxes with 1024 AN molecules. Solvent structure around Cl− and Na+ ions are studied and compared with bibliographical data. The ion-AN center of mass rdf for both vdW fittings are shown in Figure 1. For both ions, there is a well-defined first solvation shell. However, the solvation shell of Na+ displays a more distinct structure than that of Cl−. Specifically, the set1 fitting parameters, which lead a * difference between ions, presents a first maximum bigger ΔGsolv on the Na+ g(r) function plot that is higher and shifted toward smaller ion−solvent distances than that of halide ion. We used the positions of the first rdf minima to determine the radius of solvation shells (Rshell). Since the running coordination numbers (n(r)) gives the mean number of molecules within a sphere of radius r centered on the ion, the coordination number (CN) were obtained as the plateau value of (n(r)) at distances of Rshell (see Figure 1). These data are summarized in Table 2. From simulations on Cl− and Na+ ions we obtained a CN of 8.9

Figure 1. Ion-solvent center of mass radial distribution functions and running coordination. Results derived with the new proposed LJ parameters. Na+ ion in solid brown line and Cl− ion in dashed black line.

and 6.1 using the set1 fitting parameters, respectively. Comparing CN values with the available bibliography, these are closer to theoretical data but still away from experimental values of 4 and 3.2 for the Na+ from FTIR spectroscopy40 and Cl− from thermodynamic data,14,41 respectively. Table 2 lists solvent structure differences between the Cl− and Na+ ions. On the last two decades few bibliographic data were published on the study of structural data of ions in AN, being alkali and halide ions the most referred. As a reference point, a comparison between new derived LJ parameters with bibliographical data has been summarized in Table 2. In general, rdf’s using new parametrized ions are similar to previously published distribution functions for AN. However, available data are very mixed by using different methodologies and solvent models, for example, polarizable AN14,19,20 and nonpolarizable AN.3,17,19,39 In fact, in most cases used LJ parameters were fitted based on the Marcus tabulation.3,14,17,19,39 More recently, since the ΔG*hyd(H+) issue has been solved,26 additional studies using the Kelly tabulation can be found.13b,20,28a,42 However, in some cases extrapolation of 10516

dx.doi.org/10.1021/jp402545g | J. Phys. Chem. B 2013, 117, 10513−10522

The Journal of Physical Chemistry B

Article

Table 3. Lennard-Jones Parameters (R* in Å, ε in kcal/mol) Derived in this Work, Ionic charge, Solvation Free Energy (ΔGsolv * in kcal/mol), Coordination Number (CN), and the Average Ion−Nitrogen Distance (dion−N in Å) using a 512 and 1024 AN Simulation Box (Standard Error of ΔG*solv Is Also Shown) ΔG*solv

a

ions

q

R*/2

ε

Cl− Na+ Cu+ Cu2+

−1 +1 +1 +2

2.240 1.400 0.230 0.413

0.004 0.050 63.0 410.0

−64.5 −99.7 −153.7 −495.1

Zn2+

+2

0.4565

285.0

−476.7 ± 2.9

b

c

512AN ± ± ± ±

1024AN 0.5 0.8 0.8 2.3

dion−N

CN exp.

MD

exp.

MD

0.2 0.3 0.5 1.4

−64.5 −99.7 −153.8 −495.0

8.9 6.1 4.2 6.0

3.2a 4b 4c 6d

5.32 2.29 1.78 1.95

−476.9 ± 1.4

−477.2

6.3

6e

1.99

−64.2 −99.3 −153.1 −493.1

± ± ± ±

d

exp.

1.99c 1.99(eq)d 2.23,2.40(ax)d 2.11e

e

From ref 14. From ref 40. From refs 8b and 10 From refs 8, 10, and 45 From ref 9.

occupation of d orbitals on metal ions also affects the solvation structure because of the difference in the ligand field stabilization contribution. Thus, Cu2+ (d9) presents an axially elongated octahedral coordination, Cu+ (d10) a tetrahedral coordination in AN, but Zn2+ (d10) has been shown an octahedral structure in AN.8b,10 The additional difficulty that presents these transition metal ions has been previously addressed on aqueous environment.2,11,12 Also, beside usual LJ parametrization, polarization of model solvent has been considered11 and more recently a water−copper distance-dependent correction term was added on the force field to address metal−solvent interaction.44 However, to our knowledge, there is no previous LJ parametrization for transition metals ions of Cu and Zn in acetonitrile. The ion-AN LJ parameters in Table 3 lists the best fit achieved in order to reproduce not only experimental solvation free energy but coordination number and ion−N distances. From Table 3 one can be easily derive an averaged cumulative error of about 10% between the three transition metal ions when is compared with experimental data of solvent structure and thermodynamical data. Despite the simplicity of the used LJ model, a larger error of about 23% was derived when existent LJ parameters, specifically trained for water, were used to simulate ion-solvation in a box of 512 AN molecules (not shown). There is a clear contrast between the complex behavior exhibiting these ionic species and the model simplicity to parametrize the vdW interactions on usual molecular dynamics simulations. On the two steps of thermodynamical integration, the first stage of charge neutralization, where the Coulomb interactions between metallic ion and the AN model are taken into account, becomes the most important contribution to the solvation free energy with about 98% of contribution for the alkali and halide ions. On the other hand, for transition metal ions with a more complex solvent interaction, the contribution of charge neutralization to the solvation free energy is only about 88−90%. This is reflected on the higher ε and shorter R* obtained to fit a much higher solvation free energy than those of previous considered ions and still keeping target structure. However, in some cases a compromise between thermodynamics and structure has been conducted. Cu+. The ion-solvent center of mass rdf (Figure 2a) using the new proposed LJ parameters showed that the first solvation shell ends at 3.45 Å. Thus, the corresponding first shell CN is 4.2, close to experimental data. From the Cu+−N rdf (Figure 2b) can be derived the most probable Cu−N distance (dCu−N) being the first peak located at 1.78 Å, which is 10.5% lower from experimental value of 1.99 Å.8b,10 This distance penalty of 0.21 Å on the new proposed LJ parameters allows to fit

ion−water LJ parameters were used in AN solvation studies also.3 The new proposed parametrization fits with the experimental ΔGsolv * and presents a closer approximation to experimental CN than those of previous parametrizations using similar solvent models. However, the reported CN values are higher than those studies involving polarizable AN models (about 10%) and even more to experimental values. In fact, this observable was not considered on the LJ parametrization of Na+ and Cl− ions in present work because of uncertainty on its experimental values. Several points come together to explain the differences between experimental and simulated values of CN. On one hand, it is well-known that observed CN increases with decreasing concentration. Seo et al.41 point out that the difficult to determine the number of primary solvation sites around the Na+ ion in acetonitrile is due to lack of experimental data at low concentration. In fact, the highest experimental value (CN = 2.37) derived by them on their work at concentrations of about 0.6 m (molality) leads them to propose an even higher value on the infinite solution, perhaps around 4, similar than those expected for Li+ in AN.41 On the other hand, theoretical CN values are usually derived at infinite dilution whereas in experimental measures, which are derived at finite concentration, the presence of the counterion into the first solvation shell leads to a decrement of CN on those weakly solvated ions.19,22 However, used force field also takes its role at this point when observing the lower values of CN obtained by using polarizable force fields on theoretical simulations (about 5.5, see Table 2). Similarly, the use of a more elaborated nonbonded potential including three-body interactions leads to a lower value of coordination (CN = 5.220). Consequently, those differences have been attributed to the cumulative reasons of concentration, solvent model, and the use of different nonbonded potentials modeling vdW interactions. Cu+, Cu2+, and Zn2+ Ions in Acetonitrile. The solvation study of transition metal ions is particularly difficult to be modeled because of the ligand field stabilization effects and of the charge transfer from the ligands to the metal.7a,43 In AN, the solvation scheme is more complicated. Transfer thermodynamic studies show that, whereas univalent d10 ions (e.g., Cu+) are strongly solvated in AN with negatives values of free energies, enthalpies, and entropies of transfer to acetonitrile, both Cu2+ and Zn2+ metal ions present positive values of Gibbs free energy of transfer and entalphy of transfer. Thus, Cu2+ is, as expected, more weakly solvated by AN than by the other solvents such as water, methanol, dimethyl sulfoxide, and pyridine. Cu+ ion, not stabilized through complex formation, is not so stable in solvents like water and methanol and is more stable in soft donor solvents such as AN.5 In addition, the 10517

dx.doi.org/10.1021/jp402545g | J. Phys. Chem. B 2013, 117, 10513−10522

The Journal of Physical Chemistry B

Article

derived from EXAFS spectroscopy of 3.12 and 4.61 Å, respectively.10 Also, thermodynamic integration perfectly fits experimental ΔG*solv from the Kelly tabulation using 512 AN molecules in the solvation box. A little difference of 0.4% (2 kcal/mol) is appreciated when a bigger box with 1024 AN molecules is used to perform thermodynamic integration, but still within the accepted uncertainties associated with experimental data of 2−3 kcal/mol.46 Zn2+. The results of a large trajectory of Zn2+-AN system to study the ion−solvent structure using the new set of proposed LJ parameters for this ion are shown in Figure 2. In the same way as for the above transition metal, Zn2+ ion also shows a stable octahedral structure, in concordance with the observed CN at the first solvation shell of 6.3 (Figure 2a). Thus, the first solvation shell ends at 3.88 Å. Moreover, the first peak from Zn2+−N rdf (Figure 2b) is located at 1.99 Å, close to experimental data derived from XAFS spectra measurements by Inada et al.8b that point to dZn−N = 2.11 Å. In addition, those authors noted that, although it has not been probed a stable octahedral structure in AN, the Zn−N observed bond distance is quite reasonable for the 6-coordinate geometry. Thermodynamic integrations using both size of simulation box perfectly fit experimental values with errors of 0.1% (0.3−0.5 kcal/mol) on the ΔGsolv * . We want to notice that any modification on the LJ parameters involving an increase on the dZn−N, keeping experimental ΔG*solv as fitting target, leads to a rapid increase on the coordination number. Thus, it is necessary a compromise between equilibrium distance and CN when selecting new LJ parameters. This translates to an observed absolute error of 0.1 Å on the dZn−N regarding experimental data. Self-Diffusion Coefficients. To investigate the mobility of the ions, the self-diffusion coefficients (Dion and DAN) were calculated from the corresponding mean square displacement functions (MSD) by following the Einstein relation

Figure 2. Ion-solvent center of mass radial distribution functions and running coordination (a), and ion−nitrogen radial distribution functions (b). Results derived with the new proposed LJ parameters of Cu+ (brown line), Cu2+ (blue line), and Zn2+ (green line).

experimental ΔGsolv * within an error of 0.4% and a proper CN. However, when other set of optimized LJ parameters were tested that reproduced more precisely experimental ΔG*solv as well as experimental Cu−N distance, the CN rose up to 6.0, being structural data very sensible to the Cu−N distance. As seen from previous section, the Cu+ ion is the strongest solvated in AN from those studied transition metal ions and showing the largest deviation from structural experimental data. From additional Cu+−CN and Cu+−CMe rdfs (not shown), we obtained the most probable Cu−CN and Cu−CMe distances, which are located at 2.90 and 4.30 Å, respectively. Those distances are in agreement with experimental values derived from EXAFS spectroscopy of 3.01 and 4.41 Å, respectively.10 Cu2+. Figure 2a, using the new proposed LJ parameters, plots the Cu2+-AN center of mass rdf with a well-defined first solvation shell, that ends at 3.87 Å, and a CN of 6.0 in concordance with experimental values. Also, a more structured bulk-solvent can be observed if it is compared with that of Cu+ ion. Analyzing Cu2+−N rdf (Figure 2b) we can observe a first peak at dCu−N = 1.95 Å which is in agreement with experimental results. However, it cannot be reproduced in simulation with the current parameters the distorted octahedral structure being shown from experimental results derived by fluorescent EXAFS measurements (dCu−N(eq)=1.99 and dCu−N(ax)=2.40 Å) and with those obtained by EXAFS spectrum measured in the transmission mode (dCu−N(eq)=1.99 and dCu−N(ax)=2.23 Å).8b From additional Cu2+-CN and Cu2+-CMe rdf’s (not shown) we obtained the most probable Cu−CN and Cu−CMe distances, which are located at 3.08 and 4.50 Å, respectively. Those distances are in very good agreement with experimental values

D = lim

t →∞

1 2 ⟨| r (⃗ t ) − r (0) ⃗ |⟩ 6t

(4)

The system-size dependence of the diffusion coefficient is well-known and the effect was studied theoretically by Dünweg et al.,47 who found to first order the diffusion coefficient to be inversely proportional to the linear box dimensions. Later, Spångberg et al.48 proposed a simple method to correct the ion diffusion coefficients derived from the previous Einstein relation but using the following equation corr sim Dion = Dion

exp DAN sim DAN

(5)

This simple procedure allows to reduce the system-size dependence of the final ionic values, as well as corrects for errors in the solvent model itself.48 However, long simulation times are required to obtain reliable estimates of the selfdiffusion coefficient of a solute at infinite dilution. Wang et al.49 proposed a practical sampling strategy in order to obtain a reliably calculation of Dion but using several independent short trajectories. Following this approach, 20 independent MD sampling runs on a periodical box containing 512 AN molecules were performed with trajectories of 0.5 ns each (time step 0.5 fs) in a NVE ensemble after equilibration at 298 K. Then, a mean MSD are calculated by running average of MSD on the 20 trajectories after dividing the data into ten independent blocks. Finally, Dion is estimated by a least-squares fitting of mean MSD simulation time. Here the standard error 10518

dx.doi.org/10.1021/jp402545g | J. Phys. Chem. B 2013, 117, 10513−10522

The Journal of Physical Chemistry B

Article

of the mean calculated from averages is provided in order to give an estimate of the error contained in these values. The self-diffusion coefficient of both the AN bulk and simulated ionic solutions in AN are collected in Table 4. The Table 4. Self-Diffusion Coefficients for AN (Dsim AN ) and the a Ions (Dsim ) in the AN Bulk and Ionic Solutions ion Dsim ion

model AN Na+

Cl‑ Cu+ Cu2+ Zn2+

this work this work P2c E2c MDd this work MDd this work this work this work

Dsim AN

Dcorr ion

3.09 ± 0.00 1.18 1.42 1.23 1.2 1.82 1.8 1.44 0.74 0.76

± 0.02 ± 0.29 ± 0.14

4.29 ± 0.03 4.39 ± 0.05 2.6

± 0.02 2.6 ± 0.02 ± 0.01 ± 0.01

1.67 1.44 1.22 2.0e 2.57 3.0e 2.04 1.05 1.08

± 0.03 ± 0.29 ± 0.14

Exp. 4.37b 2.05f

± 0.03

2.50f

± 0.03 ± 0.02 ± 0.02

1.72f 1.32f 1.30f

Figure 3. Solvation free energy evolution against metal−nitrogen distance of Cu2+ in AN on the scanned points (R*, ε) of the free energy surface. All points have ε values lower than 1.0 kcal/mol. Two different AN models are shown; brown squares as the six-site model of Grabuleda et al.32 and blue diamonds as the six-site model of Nikitin et al.31

a Also, ionic self-diffusion coefficients extrapolated to infinite system size, using the experimental bulk value, are also shown (Dcorr ion ). All units are in 10−9 m2/s and standard error is shown. bHolz et al. ref 50. cP2: pure pair (two-body). E2: effective pair potential derived from ab initio calculations of many distorted ion-AN clusters, Spångberg et al. ref 20. d Guardia et al. ref 17. eValues corrected by us. fExperimental values calculated from conductivity measurements ref 51 and the relation Di = kBT Λi/|zi|.

parametrization (ε < 1 kcal/mol). Also, taking into account that Coulomb interaction between Cu2+ and the AN model contributes in more than 98% to ΔG*solv, the main differences between the two AN models are located on the point-charge distribution. In spite of Nikitin AN model allows to recover major part of experimental ΔGsolv * using only the electrostatic interaction, it is not possible to reach experimental thermodynamic values unless a high penalty is being paid on the solvation structure (i.e., ∼1 and ∼1.4 Å of Cu−N equilibrium distance for both models, respectively, instead of expected experimental ∼2 Å, as well as a very strange solvation structure around the solvent). This point has been compounded with the new widely accepted values of ΔGhyd * of Kelly (see Table 1) which are almost about 20 kcal/mol higher for transition metals that perhaps those reported by Marcus, the latter being the most used in previous ion parametrization.2 Consequently, a new set of LJ parameters are proposed in this work where the vdW contribution to ΔGsolv * is increased up to 10−12%. This election leads to higher ε parameters than usual, for example, 410 kcal/mol for Cu2+. In order to illustrate this point, we have computed the interaction energy (BE) between Cu2+ and one AN molecule in gas phase by using three different computational approaches. Figure 4 shows the ab initio BE, taken as a reference point, which was computed at the MP2/6-311+G** level. The BE is compared with that obtained using the proposed LJ parameters for Cu2+ and the Nikitin AN model. Thus, ab initio BE presents a minimum at 1.817 Å with a BE = −136.9 kcal/mol after correcting the basis superposition set (BSSE). On the other hand, classical interaction curve shows a similar behavior comparing with ab initio calculation but with a minimum located at 1.9 Å and a BE = −83 kcal/mol. As can be seen from Figure 4, the vdW interaction (green triangles) behaves as expected, avoiding an approaching between particles and contributing about −9.6 kcal/mol to the BE at the equilibrium distance, and vanishing quickly as Cu−N distance increases. Despite of the high value of ε parameter, the physical behavior on the metal-AN interaction still remains mainly Coulombic (about 88%) with the solvent, but in turn its solvation structure, dynamics, and thermodynamics remain close to experimental values. Transferability of Lennard-Jones Sets. Two additional AN models have been used in order to assess if the ion-AN LJ

standard error obtained from statistical data derived from of the mean calculated from averages after dividing the data into 10 independent blocks is low. Also, this methodology allows obtaining an acceptable linearity between the MSD of ions in solution versus time to determine self-diffusivity coefficients (R2 + > 0.993). As can be seen, Dsim ion (Na ) presents close values than those obtained with similar vdW parametrization (i.e., P2 and E2, and MD). However, the corrected self-diffusion coefficients are higher than those previously computed values20 but lower than obtained from Guardia et al.,17 which match experimental data after self-diffusion coefficient correction. On the other − hand, the Dcorr ion (Cl ) value obtained in this work matches experimental results after size correction, whereas transition metal ions present similar deviation to experimental data. Metal−Acetonitrile Interaction. It is well-known the strong interaction of AN with univalent d10 acceptors such as transition metal ions being especially pronounced for Cu+.7a That makes difficult to treat those interactions with a simple LJ potential plus the coulomb interaction between point-charges assigned to particles on a classical model. Usually, on the MD force field the vdW parametrization is mainly oriented to avoid excessive approaching between unbonded particles. Consequently, they are not required too large values for the ε parameter that in turn do not contribute that much to the total energy of the interacting particles around of the equilibrium distances. In fact, this behavior it is more pronounced for the case of free ions in solution where only the coulomb interaction between particles seems to be enough to reproduce usual structure and solvation energy. However this assumption strongly depends on the point-charges assignation for a given system. In the present case, a proper election on the charge distribution of the solvent can makes a difference. In order to illustrated this behavior, Figure 3 shows the solvation free energy of Cu2+ in AN against the equilibrium distance of Cu−N, with two different models of AN; Grabuleda et al. and Nikitin et al. but using usual ε values on the LJ 10519

dx.doi.org/10.1021/jp402545g | J. Phys. Chem. B 2013, 117, 10513−10522

The Journal of Physical Chemistry B

Article

interactions in water.1,2,52 Regarding solvation structure around the ion, all three AN models show results in very good agreement with available experimental data. Then, taking into account that solvation structure and thermodynamics are similar between all three models, it is suggested the transferability of the new ion-AN parameters proposed in this work to the others two AN models.



CONCLUSIONS We have presented a new set of optimized LJ parameters based on TI MD simulations to reproduce currently available structural and thermodynamic data of halide ion Cl−, alkali metal ion Na+, and some transition metal ions such as Cu+, Cu2+, and Zn2+. All derived parameters have been obtained using a recent all-atom AN model which reproduce the observed solvation free energy in excellent agreement with experimental ΔG*solv values based on tabulated ΔG*hyd by Kelly within an error of 0.4%. Structural data of ion-AN solvation obtained are in agreement with experimental data. However, there is dependence between CN and ion−N distance on transition metal ions that lead to a compromise between those structural data to optimize LJ parameters. Concretely, this behavior is more pronounced in Cu+ ion which is strongly solvated in AN, though such behavior is also observed on Zn2+ ion with a moderate increment on CN when Zn−N distance rise up to expected values. Transferability of LJ parameters has been studied using two different models of AN: a three- and six-site AN models. Both models show ion solvation structure in agreement with experimental values, and the new LJ set can be transferred with an average error on ΔG*solv of 4.8 and 7.3% in the new three- and six-site AN models, respectively.

Figure 4. Evolution of the interaction energy of a [Cu-AN]2+ cluster versus the Cu−N distance by using three different interaction energy calculations; the brown squares show a MP2/6-311+G** level of calculation, blue diamonds is using the proposed parameterization of Cu2+ (Coulomb+vdW), and green triangles show only the vdW interaction energy.

parameters developed here could be directly transferred to other six-site (Grabuleda et al.,32 all atom) and three-site (Blas et al.,33 and Guardia et al.,34 united atom) AN models. Tested models are usually used to simulate AN systems being compatible as well with AMBER force field. Thus, additional 2 ns simulation and its thermodynamic integration on all involved halide, alkali, and transition metal ions were conducted. Table 5 lists all structural and thermodynamic results comparing the three used models using a big simulation box containing 1024 AN molecules. The main feature that differentiates the three AN models among them is the different distribution of point charges of each one. Thus, both different charge distribution and the important weight that represents Coulomb interaction on the solvation free energy lead to a different error on the recovering of experimental solvation free energy. In fact, whereas the new proposed LJ parameters, using the Nikitin et al.31 model of AN, reproduce experimental ΔG*solv of Kelly (Table 1) with a mean error of 0.4%, the united-atom model present an average error of 4.8% and the all-atom model from Grabuleda et al.32 shows an average error of about 7.3%. Moreover, transition metal dications present higher errors, mainly due to a more complicated modeling of the ligand field stabilization effects and charge transfer. It is interesting to notice that united-atom model reproduce ΔGsolv * close to experimental data based on ΔG*hyd tabulated by Marcus (Table 1), which have been previously used as reference values to parametrize vdW



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work has been supported by MICINN and FEDER funds (MAT2009-09138 and MAT2009-11513) and by the DIUE of the Generalitat de Catalunya (2009SGR925, 2009SGR1208, and XRQTC). The authors are indebted to the Barcelona Supercomputer Center (BSC) for the computational resources provided. Support for the research of C.A. was received through the prize “ICREA Academia” for excellence in research funded by the Generalitat de Catalunya.

Table 5. Thermodynamics and Structure Comparison between Three Different AN Solvent Models Using the New Proposed Lennard-Jones Parametersa six-site AN model 1b

ions −

Cl Na+ Cu+ Cu2+ Zn2+

six-site AN model 2c

three-site AN model 3d

ΔG*solv

dion−N

CN

ΔG*solv

dion−N

CN

ΔG*solv

dion−N

CN

−64.2(0.4) −99.2(0.5) −153.1(0.4) −493.1(0.4) −476.9(0.1)

5.32 2.29 1.78 1.95 1.99

8.9 6.1 4.2 6.0 6.3

−61.6(4.5) −94.0(5.8) −139.6(9.2) −452.0(8.7) −436.8(8.5)

5.46 2.46 1.91 2.11 2.13

8.5 5.9 4.1 6.0 6.0

−59.9(7.2) −98.2(1.5) −146.3(4.9) −469.3(5.2) −453.3(5.0)

5.23 2.45 1.87 2.06 2.10

8.0 6.1 4.1 6.0 6.0

Solvation free energy (ΔGsolv * in kcal/mol), coordination number (CN), and the average ion−nitrogen distance (dion−N in Å) using a 1024 AN simulation box. ΔG*solv errors are shown in parentheses (in %). bNikitin et al. ref 31. cGrabuleda et al. ref 32. dBlas et al. ref 33 and Guardia et al. ref 34. a

10520

dx.doi.org/10.1021/jp402545g | J. Phys. Chem. B 2013, 117, 10513−10522

The Journal of Physical Chemistry B



Article

Test of Different Combination Rules. Ber. Bunsen-Ges. Phys. Chem. 1989, 93 (12), 1464−1467. (16) Pálinkás, G.; Riede, W. O.; Heinzinger, K. Z. Naturforsch. Teil A 1977, 32, 1137. (17) Guàrdia, E.; Pinzón, R. On the solvation shell of Na+ and Cl− ions in acetonitrile: A computer simulation study. J. Mol. Liq. 2000, 85 (1−2), 33−44. (18) Cabaleiro-Lago, E. M.; Ríos, M. A. Ab initio study of M(CH3CN)n clusters (M = Li+, Na+, Mg2+) in the gas phase. Chem. Phys. 2000, 254 (1), 11−23. (19) Fischer, R.; Richardi, J.; Fries, P. H.; Krienke, H. The solvation of ions in acetonitrile and acetone. II. Monte Carlo simulations using polarizable solvent models. J. Chem. Phys. 2002, 117 (18), 8467−8478. (20) Spångberg, D.; Hermansson, K. The solvation of Li+ and Na+ in acetonitrile from ab initio-derived many-body ion−solvent potentials. Chem. Phys. 2004, 300 (1−3), 165−176. (21) Hsiao, Y.-W.; Ryde, U. Interpretation of EXAFS spectra for sitting-atop complexes with the help of computational methods. Inorg. Chim. Acta 2006, 359 (4), 1081−1092. (22) (a) Seo, D. M.; Borodin, O.; Han, S.-D.; Ly, Q.; Boyle, P. D.; Henderson, W. A. Electrolyte Solvation and Ionic Association. J. Electrochem. Soc. 2012, 159 (5), A553−A565. (b) Seo, D. M.; Borodin, O.; Han, S.-D.; Boyle, P. D.; Henderson, W. A. Electrolyte Solvation and Ionic Association II. Acetonitrile-Lithium Salt Mixtures: Highly Dissociated Salts. J. Electrochem. Soc. 2012, 159 (9), A1489−A1500. (23) (a) Steinbrecher, T.; Mobley, D. L.; Case, D. A. Nonlinear scaling schemes for Lennard-Jones interactions in free energy calculations. J. Chem. Phys. 2007, 127 (21), 214108−13. (b) Steinbrecher, T.; Joung, I.; Case, D. A. Soft-core potentials in thermodynamic integration: Comparing one- and two-step transformations. J. Comput. Chem. 2011, 32 (15), 3253−3263. (24) 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 Cluster-Ion Solvation Data. J. Phys. Chem. A 1998, 102 (40), 7787− 7794. (25) 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. (26) 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. (27) 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. (28) (a) Bryantsev, V. S.; Diallo, M. S.; Goddard, W. A., III Calculation of Solvation Free Energies of Charged Solutes Using Mixed Cluster/Continuum Models. J. Phys. Chem. B 2008, 112 (32), 9709−9719. (b) Ben-Naim, A. Standard Thermodynamics of Transfer. Uses and Misuses. J. Phys. Chem. 1978, 82 (7), 792−803. (29) Case, D. A.; Darden, T. A.; T.E. Cheatham, I.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M.; Roberts, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Swails, J.; Goetz, A. W.; Kolossváry, I.; Wong, K. F.; Paesani, F.; Vanicek, J.; Wolf, R. M.; Liu, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Cai, Q.; Ye, X.; Wang, J.; Hsieh, M.-J.; Cui, G.; Roe, D. R.; Mathews, D. H.; Seetin, M. G.; Salomon-Ferrer, R.; Sagui, C.; Babin, V.; Luchko, T.; Gusarov, S.; Kovalenko, A.; Kollman, P. A. AMBER 12; University of California: San Francisco, 2012. (30) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117 (19), 5179−5197. (31) Nikitin, A. M.; Lyubartsev, A. P. New six-site acetonitrile model for simulations of liquid acetonitrile and its aqueous mixtures. J. Comput. Chem. 2007, 28 (12), 2020−2026.

REFERENCES

(1) Horinek, D.; Mamatkulov, S. I.; Netz, R. R. Rational design of ion force fields based on thermodynamic solvation properties. J. Chem. Phys. 2009, 130 (12), 124507−21. (2) Babu, C. S.; Lim, C. Empirical Force Fields for Biologically Active Divalent Metal Cations in Water. J. Phys. Chem. A 2006, 110 (2), 691− 699. (3) Troxler, L.; Wipff, G. Conformation and Dynamics of 18-Crown6, Cryptand 222, and Their Cation Complexes in Acetonitrile Studied by Molecular Dynamics Simulations. J. Am. Chem. Soc. 1994, 116 (4), 1468−1480. (4) Sazonov, V. P.; Shaw, D. G.; Sazonov, N. V.; Skrzecz, A.; Lisov, N. I. IUPAC-NIST Solubility Data Series. 78. Acetonitrile Binary Systems. J. Phys. Chem. Ref. Data 2002, 31 (4), 989−1133. (5) Chaudhry, M.; Persson, I. Transfer thermodynamic study on the copper(II) ion from water to methanol, acetonitrile, dimethyl sulfoxide and pyridine. J. Chem. Soc., Faraday Trans. 1994, 90 (15), 2243−2248. (6) Reimers, J. R.; Hall, L. E. The Solvation of Acetonitrile. J. Am. Chem. Soc. 1999, 121 (15), 3730−3744. (7) (a) Persson, I. Solvation and complex formation in strongly solvating solvents. Pure Appl. Chem. 1986, 58 (8), 1153−1161. (b) Johnsson, M.; Persson, I. Determination of Gibbs free energy of transfer for some univalent ions from water to methanol, acetonitrile, dimethylsulfoxide, pyridine, tetrahydrothiophene and liquid ammonia; standard electrode potentials of some couples in these solvents. Inorg. Chim. Acta 1987, 127 (1), 15−24. (8) (a) Inada, Y.; Sugimoto, Y.; Nakano, Y.; Itoh, Y.; Funahashi, S. Formation and Deprotonation Kinetics of the Sitting-Atop Complex of Copper(II) Ion with 5,10,15,20-Tetraphenylporphyrin Relevant to the Porphyrin Metalation Mechanism. Structure of Copper(II)−Pyridine Complexes in Acetonitrile As Determined by EXAFS Spectroscopy. Inorg. Chem. 1998, 37 (21), 5519−5526. (b) Inada, Y.; Nakano, Y.; Inamo, M.; Nomura, M.; Funahashi, S. Structural Characterization and Formation Mechanism of Sitting-Atop (SAT) Complexes of 5,10,15,20-Tetraphenylporphyrin with Divalent Metal Ions. Structure of the Cu(II)−SAT Complex As Determined by Fluorescent Extended X-ray Absorption Fine Structure. Inorg. Chem. 2000, 39 (21), 4793− 4801. (9) Inada, Y.; Niwa, Y.; Iwata, K.; Funahashi, S.; Ohtaki, H.; Nomura, M. Solvation structure of metal ions in nitrogen-donating solvents. J. Mol. Liq. 2006, 129 (1−2), 18−24. (10) Persson, I.; Penner-Hahn, J. E.; Hodgson, K. O. An EXAFS spectroscopic study of solvates of copper(I) and copper(II) in acetonitrile, dimethyl sulfoxide, pyridine, and tetrahydrothiophene solutions and a large-angle x-ray scattering study of the copper(II) acetonitrile solvate in solution. Inorg. Chem. 1993, 32 (11), 2497− 2501. (11) Ponomarev, S. Y.; Click, T. H.; Kaminski, G. A. Electrostatic Polarization Is Crucial in Reproducing Cu(I) Interaction Energies and Hydration. J. Phys. Chem. B 2011, 115 (33), 10079−10085. (12) Li, X.; Tu, Y.; Tian, H.; Agren, H. Computer simulations of aqua metal ions for accurate reproduction of hydration free energies and structures. J. Chem. Phys. 2010, 132 (10), 104505. (13) (a) McDonald, N. A.; Duffy, E. M.; Jorgensen, W. L. Monte Carlo Investigations of Selective Anion Complexation by a Bis(phenylurea) p-tert-Butylcalix[4]arene. J. Am. Chem. Soc. 1998, 120 (20), 5104−5111. (b) Böes, E. S.; Livotto, P. R.; Stassen, H. Solvation of monovalent anions in acetonitrile and N,N-dimethylformamide: Parameterization of the IEF-PCM model. Chem. Phys. 2006, 331 (1), 142−158. (c) de Araujo, A. S.; Piro, O. E.; Castellano, E. E.; Danil de Namor, A. F. Combined Crystallographic and Solution Molecular Dynamics Study of Allosteric Effects in Ester and Ketone p-tertButylcalix[4]arene Derivatives and Their Complexes with Acetonitrile, Cd(II), and Pb(II). J. Phys. Chem. A 2008, 112 (46), 11885−11894. (14) Richardi, J.; Fries, P. H.; Krienke, H. The solvation of ions in acetonitrile and acetone: A molecular Ornstein–Zernike study. J. Chem. Phys. 1998, 108 (10), 4079−4089. (15) Wiese, H.; Brickmann, J. Lennard-Jones (12,6) Potentials for the Non-Ionic Interaction of Alkali Cations and Halide Anions  A 10521

dx.doi.org/10.1021/jp402545g | J. Phys. Chem. B 2013, 117, 10513−10522

The Journal of Physical Chemistry B

Article

J. Am. Chem. Soc. 1965, 87 (12), 2529−2534. (b) Harkness, A. C.; Daggett, H. M., Jr. The Electrical Conductivities of Some Tetra-nalkylammonium Salts in Acetonitrile. Can. J. Chem. 1965, 43 (5), 1215−1221. (c) Yeager, H. L.; Kratochvil, B. Conductance study of copper(I) and silver(I) salts of symmetrical anions in acetonitrile. J. Phys. Chem. 1969, 73 (6), 1963−1968. (d) Diamond, A.; Fanelli, A.; Petrucci, S. Kinetics of association of nickel(II), copper(II), and zinc(II) perchlorates in acetonitrile. Inorg. Chem. 1973, 12 (3), 611− 619. (52) Joung, I. S.; Cheatham, T. E. Determination of Alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations. J. Phys. Chem. B 2008, 112 (30), 9020− 9041.

(32) Grabuleda, X.; Jaime, C.; Kollman, P. A. Molecular dynamics simulation studies of liquid acetonitrile: New six-site model. J. Comput. Chem. 2000, 21 (10), 901−908. (33) Blas, J. R.; Márquez, M.; Sessler, J. L.; Luque, F. J.; Orozco, M. Theoretical Study of Anion Binding to Calix[4]pyrrole: the Effects of Solvent, Fluorine Substitution, Cosolute, and Water Traces. J. Am. Chem. Soc. 2002, 124 (43), 12796−12805. (34) Guàrdia, E.; Pinzón, R.; Casulleras, J.; Orozco, M.; Luque, F. J. Comparison of Different Three-site Interaction Potentials for Liquid Acetonitrile. Mol. Simul. 2001, 26 (4), 287−306. (35) Blondel, A. Ensemble variance in free energy calculations by thermodynamic integration: Theory, optimal “Alchemical” path, and practical solutions. J. Comput. Chem. 2004, 25 (7), 985−993. (36) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23 (3), 327−341. (37) Berendsen, H. J. C.; Postma, J. P. M.; Gunsteren, W. F. v.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81 (8), 3684−3690. (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 (12), 10089−10092. (39) Markovich, G.; Perera, L.; Berkowitz, M. L.; Cheshnovsky, O. The solvation of Cl−, Br−, and I− in acetonitrile clusters: Photoelectron spectroscopy and molecular dynamics simulations. J. Chem. Phys. 1996, 105 (7), 2675−2685. (40) Barthel, J.; Deser, R. FTIR study of ion solvation and ion-pair formation in alkaline and alkaline earth metal salt solutions in acetonitrile. J. Solution Chem. 1994, 23 (10), 1133−1146. (41) Seo, J.-S.; Cheong, B.-S.; Cho, H.-G. Solvation of LiClO4 and NaClO4 in deuterated acetonitrile studied by means of infrared and Raman spectroscopy. Spectrochim. Acta, Part A 2002, 58 (8), 1747− 1756. (42) Riccardi, D.; Guo, H.-B.; Parks, J. M.; Gu, B.; Liang, L.; Smith, J. C. Cluster-Continuum Calculations of Hydration Free Energies of Anions and Group 12 Divalent Cations. J. Chem. Theory Comput. 2012, 9 (1), 555−569. (43) Rode, B. M.; Schwenk, C. F.; Hofer, T. S.; Randolf, B. R. Coordination and ligand exchange dynamics of solvated metal ions. Coord. Chem. Rev. 2005, 249 (24), 2993−3006. (44) Zang, J.; Nair, S.; Sholl, D. S. Prediction of Water Adsorption in Copper-Based Metal−Organic Frameworks Using Force Fields Derived from Dispersion-Corrected DFT Calculations. J. Phys. Chem. C 2013, 117 (15), 7519−7525. (45) Inamo, M.; Kamiya, N.; Inada, Y.; Nomura, M.; Funahashi, S. Structural Characterization and Formation Kinetics of Sitting-Atop (SAT) Complexes of Some Porphyrins with Copper(II) Ion in Aqueous Acetonitrile Relevant to Porphyrin Metalation Mechanism. Structures of Aquacopper(II) and Cu(II)−SAT Complexes As Determined by XAFS Spectroscopy. Inorg. Chem. 2001, 40 (22), 5636−5644. (46) 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. (47) Dünweg, B.; Kremer, K. Molecular dynamics simulation of a polymer chain in solution. J. Chem. Phys. 1993, 99 (9), 6983−6997. (48) Spångberg, D.; Hermansson, K. Effective three-body potentials for Li+(aq) and Mg+(aq). J. Chem. Phys. 2003, 119 (14), 7263−7281. (49) Wang, J.; Hou, T. Application of molecular dynamics simulations in molecular property prediction II: Diffusion coefficient. J. Comput. Chem. 2011, 32 (16), 3505−3519. (50) Holz, M.; Mao, X.-a.; Seiferling, D.; Sacco, A. Experimental study of dynamic isotope effects in molecular liquids: Detection of translation-rotation coupling. J. Chem. Phys. 1996, 104 (2), 669−679. (51) (a) Coetzee, J. F.; Cunningham, G. P. Evaluation of Single Ion Conductivities in Acetonitrile, Nitromethane, and Nitrobenzene Using Tetraisoamylammonium Tetraisoamylboride as Reference Electrolyte. 10522

dx.doi.org/10.1021/jp402545g | J. Phys. Chem. B 2013, 117, 10513−10522