Article pubs.acs.org/JPCB
Solvation Energy of Ions in Polymers: Effects of Chain Length and Connectivity on Saturated Dipoles near Ions Lijun Liu† and Issei Nakamura*,‡ †
State Key Laboratory of Polymer Physics and Chemistry, Changchun Institute of Applied Chemistry, Chinese Academy of Sciences, Changchun, Jilin 130022, China ‡ Department of Physics, Michigan Technological University, Houghton, Michigan 49931, United States S Supporting Information *
ABSTRACT: We illustrate the effects of chain connectivity on the solvation energy of ions immersed in polymer liquids by developing a new coarse-grained molecular dynamics simulation. Our theory accounts for the dielectric response of the polymers through the connection of dipolar, monomeric units with nonlinear springs. In stark contrast to the standard Born solvation energy of ions, our results depend substantially on the chain length of the polymers. We also demonstrate the marked difference in the solvation energies of the ions immersed in non-polymeric particle mixtures, single-component polymers, polymer blends, and block copolymers. Thus, we suggest that the chain architecture of polymers is a key factor in ion solvation, whereas this feature is often inadequately considered in main theory and simulation literature. Our results are consistent with those predicted by previous coarse-grained mean-field theories when the dipole moment of the polymer compositions is relatively small. However, we also demonstrate that the strong ion−dipole and dipole−dipole interactions cause the chain-like association of the monomeric units, resulting in a qualitative discrepancy between the mean-field theory and simulation. Such a strong electrostatic correlation may reverse the dependence of the chain length on the solvation energy of the ions in the polymers.
1. INTRODUCTION Ion solvation in liquids may give rise to a major driving force that causes significant changes in various thermodynamic properties. This mechanism can typically be rationalized by the delicate balance between the entropy and enthalpy of solvents and ions. Among other effects, electrostatic interactions, such as charge screening and the reorientation of solvent dipoles, are a key factor to consider when analyzing diverse phenomena driven by ion solvation. Nevertheless, the complexity of the electrostatic interactions in ion-containing liquids is often elusive in theoretical and computational modeling.1−5 In conventional molecular dynamics simulations based on implicit solvent models, the electrostatic interaction between two charges with q1e and q2e separated by the distance, r, is written in the form of the Coulomb potential, q1q2e2/(4πε0εr), with the elementary charge, e, and the vacuum permittivity, ε0. In this case, solvents are typically treated as a dielectric continuum. Note that the electrostatic interaction between the charges is modeled by introducing the bulk dielectric constant ε. However, this expression does not adequately account for the saturated dipoles near the ions caused by the strong electrostatic field around the ions, called the “dielectric saturation”.2,4−6 Accordingly, the three-body interaction between the solvent and ion is ignored.3 For example, the © XXXX American Chemical Society
dielectric functions near the like and unlike charges are substantially different because of the difference in the electrostatic fields in close proximity to the ion pairs. Generally, the dielectric function, ε(r), should be elaborated to precisely account for the true nature of the polarization in complex fluids. Indeed, our recent coarse-grained mean-field theory for ioncontaining polymers7 suggested that the correlation between the dielectric saturation and the chain connectivity can be critical in consideration of ion solvation in polymer liquids; however, the parameterization with the bulk dielectric constant, ε, cannot adequately capture this correlation. The solvation energy of ions with charge qe can be calculated q 2e 2 ⎡ 1 ⎤1 as ΔG = 8πε ∫ dr ⎣ ε(r) − 1⎦ r 2 , with the distance from the ion 0 as the coordinate, r. If we assume a spatially uniform dielectric response in a medium with ε(r) = ε, ΔG can be cast into the standard Born solvation energy form ΔG Born =
q 2e 2 8πε0a
( 1ε −1).
Here, a is the radius of the ion. Importantly, ΔGBorn decreases as the dielectric constant, ε, is increased. Thus, the ions tend to be solvated by a higher dielectric component to reduce the Received: January 20, 2017 Revised: March 14, 2017 Published: March 17, 2017 A
DOI: 10.1021/acs.jpcb.7b00671 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
between the particles is attractive when the intermolecular distance rij is larger than 21/6σ. However, in this work, we use only the repulsive part of the potential to highlight the attractive interactions arising purely from the dipolar nature and to compare the results more directly with our previous coarse-grained mean-field theory.7 Thus, to obtain the purely repulsive force, we set the cutoff radius rc = 21/6σ. We write the electrostatic interactions at the distance, r, between the ion and dipole (U iid ) and between
enthalpic cost, with its entropic cost counterbalanced. Despite its simple form, the Born solvation energy has been invoked over many decades to explain various experiments.8−10 To account for the saturated dipoles near ions, the Langevin− Debye theory may be employed to modify the solvation energy.2,3,11 However, the dielectric response of ion-containing polymers involves further complexity and is therefore inadequately represented. The fact that the correlation between the reorientation of monomer dipoles and the intramolecular connectivity of a polymer is ignored in main theory and the simulation literature is of particular concern.12 In this context, the effect of chain architecture may significantly influence the position-dependent dielectric function, ε(r). Indeed, our recent coarse-grained mean-field theories predicted this feature;7 the dielectric functions of the non-polymeric fluid, polymer blend, and block copolymer melt are qualitatively different and exhibit a notable dependence on the chain length. Nevertheless, some approximations employed in the theory remain unclear; for example, the validity of coarse-graining the solvation shells around ions on the atomistic scale is uncertain. If the dielectric function involves a substantial effect of chain architecture, the current theory and simulation for ion-containing polymers would require drastic modification to explore novel phase behaviors. In this article, we examine the hypothetical effect of chain architecture on the solvation energy of ions through the dielectric response using molecular dynamics simulations. To illustrate our points, we draw upon the Stockmayer-like fluid of polymerized dipolar particles but ignore the attractive longrange term in the dispersion force. However, in the case of real polymers, monomeric units are connected through covalent bonding that restricts the rotational orientation of the dipoles attached to the monomers. If the dipoles are located on a flexible side chain, classified as “type C” by Stockmayer,12 the dielectric response of a polymer may be qualitatively analyzed in terms of freely rotating dipoles.13 However, our rationale for invoking the Stockmayer-like fluid to model polymers is that this theory can succinctly account for the spatial inhomogeneity of the dielectric response, ranging from the atomistic to the nanometer scale. Moreover, polymerization and phase separation in the Stockmayer fluid were shown to be consistent with the Flory−Huggins theory developed by Dudowicz and Freed.14 In their model, dipolar particles can freely associate to form linear polymer chains. However, the microscopic nature of dipoles and their rotational orientation relative to the backbone chain can be coarse-grained by introducing the effective Flory χ parameter over a wide range of model parameters. Thus, our model resorting to the polymerized dipolar particles should also provide valuable information and theoretical insight into the solvation energy of ions dissolved in polymers.
qe μ⃗i · r ⃗
id the dipole and dipole (Udd ij ) as follows: Ui = − 4πε
1 ⎡ ⎢μi⃗ ·μj⃗ 4πε0r 3 ⎣
r3
0
3(μi⃗ · r )( ⃗ μj⃗ · r )⃗ ⎤
and
, where μ⃗ i and μ⃗j are the dipole ⎦⎥ moments of particles i and j, respectively. To model the polymers, we connect adjacent dipolar particles with nonlinear springs based on the finite extensible n o n l i n e a r el a s t i c ( F E N E ) p o t e n t i a l a s f o l l o w s : 2⎤ ⎡ r 1 UFENE(rb) = − 2 kR 02 ln⎢1 − Rb ⎥, where rb is the distance 0 ⎦ ⎣ between the covalently bonded monomers. To avoid bond crossing, we use spring constant k = 30ϵLJ/σ2 and cutoff radius R0 = 1.5σ.15 We write the Lagrangian for the system consisting of N particles as follows: Uijdd =
−
r2
( )
L=
∑ 1 mvi⃗ 2 + 2
i
− +
·2 1 ∑ 1 Iμi⃗ + 1 ∑ ∑ μi⃗ ·T·μj⃗ 2 2 i≠j j μ i 2
1 ∑ ∑ UijWCA − 2 i≠j j
∑
UFENE +
∑ λi(μi⃗ 2
− μ2 )
i
bonds
∑ μi⃗ ·E⃗
(1)
i
where I is the moment of inertia of the dipoles, m is the particle mass, μ is the permanent dipole moment, and vi⃗ is the velocity of the particle, i. The first two terms describe the translational and rotational kinetic energies of the particles, respectively. The third, fourth, and fifth terms represent the dipole−dipole interaction, and the WCA and FENE interparticle potentials, respectively. For computational simplicity, we introduced a Lagrange multiplier λi in the sixth term to treat the rotational part of the motion;16,17 with this constraint, the magnitude of the dipole moment remains constant at each time step while solving the equations of motion in the Cartesian coordinates. The last term represents the interaction between the dipole and the external field. In the present case, this term provides the ion−dipole interaction. T is given by the dipole interaction tensor T =
1 4πε0
(
3 r ⃗r ⃗ r5
−
I r3
) with the unit matrix, I. Note that
the Lagrange equation for coordinate x is given by ∂L d ∂L − ∂x = 0. We thus obtain the equations of motion for dt ∂x ̇ the translational degrees of freedom
2. THEORY AND SIMULATION We consider a single ion with charge qe immersed in dipolar particles. The excluded volume interaction between the particles is taken into account through a truncated and shifted purely repulsive Lennard-Jones (LJ) potential [so-called Weeks−Chandler−Andersen (WCA) potential], UijWCA = 4ϵLJ[(σ/rij)12 − (σ/rij)6] + ϵLJ at rij ≤ rc and UWCA = 0 at rij ij > rc. Here, σ symbolizes the diameter of the dipolar particle, ϵLJ is the depth of the potential well in the original LJ potential, and rij = |ri⃗ − rj⃗ | is the intermolecular distance between particles i and j. Note that without the condition given by rc, the force
··
mri ⃗ = μi⃗ ·∑ ∇T·μj⃗ + μi⃗ ·∇E ⃗ + j≠i
WCA
∑ Fij⃗ j≠i
+ Fi⃗
FENE
(2)
where F⃗WCA = −∇iUWCA and F⃗FENE = −∑bonds∇iUFENE are the ij ij i forces due to the excluded volume interactions and the bond connection between the particles, respectively. For nonpolymeric particles, F⃗FENE is set to 0. Similarly, the equations i of motion for the rotational degrees of freedom are given by B
DOI: 10.1021/acs.jpcb.7b00671 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B I ·· μ⃗ = μ2 i
∑ T·μj⃗
respectively. The time step used for solving the equation of motion is Δt = 2 fs. The simple truncation of the electrostatic interactions can be comparable to the Ewald method when the effective range of the potential is relatively shorter than the system size.20 For computational simplicity, therefore, we introduce a truncation radius set to half of the box size L. Note that the average value of the ion−dipole interaction per particle in the bulk phase can be approximately evaluated by ⟨Ũ idi ⟩ = ∫ dΩid Ũ idi exp(−Ũ idi / kBT)/∫ dΩid exp(−Ũ idi /kBT), where dΩid denotes dθid sin θid and Ũ idi = Uidi /ε accounts for the effect of the dielectric screening through the bulk dielectric constant ε. Here, ε may be estimated from the Debye equation for polar gases, ε = 1 + ρμ2/(3ε0kBT), where ρ is the number density.21 When the distance between two particles is L/2 = 22.5 Å, we obtain ⟨Ũ idi ⟩ = −4.19 × 10−3kBT for the dipole moment μ = 3 D and ⟨Ũ idi ⟩ = −3.36 × 10−3kBT) for μ = 0.5 D, respectively. Given that the number of the particles at r = L/2 from the ion is estimated by M = 4πr2Δrρ with the average intermolecular distance Δr = [3/ (4πρ)]1/3, we find M to be 125.1, and hence M⟨Ũ idi ⟩ = −0.52kBT and −0.42kBT for μ = 3 and 0.5 D, respectively. Incidentally, the corresponding values at r = L are −0.13kBT and −0.11kBT for μ = 3 and 0.5 D, respectively. Moreover, M⟨Ũ dd i ⟩ for the dipole−dipole interaction is smaller than M⟨Ũ idi ⟩ because the dipole−dipole interaction is significantly weaker than the ion−dipole interaction. However, we consider all of these energy values to be negligible relative to the solvation energies of the ions, −100kBT [see Figures 2 and 6a]
+ E ⃗ + 2λiμi⃗ (3)
j≠i
We use the Velocity-Verlet algorithm to integrate the equation ·
1
··
of motion x (⃗ t + Δt ) = x (⃗ t ) + x (⃗ t )Δt + 2 x (⃗ t )Δt 2 , in which variable x⃗ represents r ⃗ or μ⃗ . Δt is the time step. We use the periodic boundary condition. The Lagrange multiplier λi is given by λi = ( −C2 + C22 − C1C3 )/C1, in which C1 = 4α2μ2Δt3, C2 = 2αΔt[μ⃗i(t)·F⃗(t)Δt + μ2 + αμ⃗i(t)·R⃗ (t)Δt2], and C3 = α2R⃗ 2(t)Δt3 + 2αμ⃗ i(t)·R⃗ (t)Δt + 2[F⃗(t)·μ⃗i(t) + αF⃗(t)· R⃗ (t)Δt2] + F⃗2(t)Δt, α = μ2/I, R⃗ (t) = ∑j≠iT·μ⃗j(t) + E⃗ , where · ·· F (⃗ t ) = μi⃗ (t − Δt ) + 0.5μi⃗ (t − Δt )Δt .16,17 The system size is set to L3 with L = 45 Å. In general, unavoidable energy drift may arise from numerical integration artifacts that are caused by the finite time step, Δt. Note that according to the energy equipartition theorem, the average 3 kinetic energy per particle is given by 2 kBT , where kB and T are the Boltzmann constant and temperature, respectively. Thus, we can avoid numerical divergence by maintaining a constant temperature. In this study, we set the temperature to be 300 K using velocity rescaling schemes.18 We calculate the solvation energy of ions as a transfer energy from a pure liquid to an ion-containing liquid. Thus, the solvation energy is given by ΔG = Uid + ΔUdd; the total ion− dipole interaction Uid ≡ ∑iUidi and the relative dipole−dipole 1 1 interaction ΔU dd ≡ 2 ∑i ≠ j ∑j Uijdd − 2 ∑i ≠ j ∑j Uijdd(q = 0) 1
1
with 2 ∑i ≠ j ∑j Uijdd and 2 ∑i ≠ j ∑j Uijdd(q = 0) defined as the total dipole−dipole interactions in the systems with and without the ions, respectively. Here, the free-energy shift is caused primarily by the enthalpic change, but we use notation ΔG instead of “ΔU” to be consistent with the main literature on ion solvation in this field of study, specifically with the generalized Born model suggested by Still et al.19 Our system consists of a single monovalent ion and 512 solvent dipoles. The diameters of the ion and solvent are 3 and 4.5 Å, respectively. Here, we define the first solvation shell by the contact between the ion and solvent, that is, it is 2.25 Å from the surface of the ion, as illustrated in Figure 1a,b. Similarly, the second solvation shell is 3.15 Å from the first solvation shell. Taking liquid water as an example, the particle mass, parameter ϵLJ in the WCA potential, and the moment of inertia are set to m = 2.99 × 10−26 kg, ϵLJ = 5. 27 × 10−21 J, and I = 0.025mσ2 = 1.514 × 10−46 kg·m2 with σ = 4.5 Å,
Figure 2. Simulation results of the non-polymeric particles. Peak value of the ion−dipole radial distribution function (red squares), coordination number of the solvents in the first solvation shell Nc (red circles), and solvation energy of the ions (blue triangles) as a function of the dipole moment. To guide the eye, many of the particles are not visualized in the snapshots.
and −10kBT (see Figure 2). Thus, the electrostatic part of the solvation energy of the ions is insignificant when the particle distance is larger than the truncation radius L/2.22 In terms of the Ewald method, this case corresponds to the situation in which the contribution from periodic images due to the periodic boundary conditions is negligible. To further validate the truncation radius, we increased L from 45 to 54 Å, increasing the truncation radius from L/2 = 22.5 to 27 Å. We found that the difference in the solvation energy is less than 4%. We fit the autocorrelation functions of the end-to-end vector, ⟨R⃗ end(t)·R⃗ end(0)⟩/⟨R2end(0)⟩, with an exponential function A e(−t/τ). Here, re⃗ nd(t) and r1⃗ (t) designate the monomer positions of two chain ends, R⃗ end(t) ≡ re⃗ nd(t) − r1⃗ (t) is the end-to-end vector, and τ denotes the chain relaxation time. In the case of
Figure 1. (a) Schematic illustration of the first and second solvation shells. θid is the angle of the dipole, μ⃗ , relative to the radial axis from the ion in the center, and θdd is the angle between the neighboring dipoles. (b) Definitions of the first and second solvation shells. C
DOI: 10.1021/acs.jpcb.7b00671 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B the chain length N = 16, τ is smaller than 104 ps. We thus equilibrate the system for 400 ps for the chain length N ≤ 16. Then, the simulation is run for 2 ns. Similarly, in the case of N = 64 and 32, we equilibrate the system for 2 ns, on the basis of τ < 0.77 ns. Then, the simulation is run for 5 ns. We averaged on the order of 104 samples.
ion−dipole and dipole−dipole interactions. In this case, the dielectric response near the ions likely weakens, as deduced from the Debye equation for polar gases ε = 1 + ρμ2/ (3ε0kBT).21 This effect should also cause an increase in ΔG. We thus suggest that the solvation energy of the ions is qualitatively measured by the peak value of g(r), whereas this feature is often ignored in theoretical modeling of ion solvation. Indeed, a previous study based on the Stockmayer fluid showed that dipolar particles exhibit a chain-like association without ions at very low temperatures.24 However, our novel finding is that the chain-like structure is stabilized in the proximity of ions at room temperature, driven by both the ion−dipole and dipole−dipole interactions. We refer to recent coarse-grained mean-field theories of ioncontaining liquid mixtures that predicted the local enrichment of a higher dielectric component near ions.25,26 Indeed, our simulation shows this feature (Figure 4). Thus, our rationale for
3. RESULTS We start from a single ion in non-polymeric particles because such a system contains some key features for understanding ion solvation in polymers. This system corresponds to a dilute ionic solution with the periodic boundary condition imposed on the simulation box. In Figure 2, we show the peak value of ion− dipole radial distribution function g(r) and the solvation energy of an ion as a function of dipole moment μ. Note that the local structure of the solvents near the ions is substantially correlated with the peak value of g(r). The decline in the solvation energy of the ions, ΔG, with relatively small dipoles (
