Atomistic Simulations of Poly(ethylene oxide) in ... - ACS Publications

Dec 23, 2013 - Department of Chemistry, Columbia University, 3000 Broadway, New York, ..... Gustavo A. Chapela , Orlando Guzmán , Enrique Díaz-Herre...
0 downloads 0 Views 1MB Size
Download PDF
Article pubs.acs.org/Macromolecules

Atomistic Simulations of Poly(ethylene oxide) in Water and an Ionic Liquid at Room Temperature Jagannath Mondal,† Eunsong Choi,‡ and Arun Yethiraj*,§ †

Department of Chemistry, Columbia University, 3000 Broadway, New York, New York 10027, United States Department of Physics, University of Wisconsin, Madison, Wisconsin 53706, United States § Department of Chemistry, University of Wisconsin, Madison, Wisconsin 53706, United States ‡

ABSTRACT: The behavior of polymers in ionic liquids is of technological and scientific interest. In this work we present atomistic simulations for the properties of isolated poly(ethylene oxide) (PEO) in the ionic liquid 1-butyl 3-methylimidazolium tetrafluoroborate ([BMIM][BF4]) and compare to the properties of the same polymer in water, at room temperature, for degrees of polymerization N ranging from 9 to 40. PEO chains are much more expanded in [BMIM][BF4] than in water. The root-mean-square radius of gyration, Rg, scales as Rg ∼ Nν with ν ≈ 0.9, and the distribution of end-to-end distance is bimodal, with coexisting extended and hairpin-like conformations. The simulations are consistent with polyelectrolyte behavior, i.e., Rg ∼ N, but the chains might be too short to be in the true scaling regime. (For comparison, Rg ∼ N0.5 in water.) [BMIM][BF4] is a much better solvent than water: In [BMIM][BF4] the solvation free energy of the monomer is 50% more negative, and the potential of mean force between two PEO 9-mers is significantly more repulsive than in water; the repulsion comes from energetic polymer−solvent interactions. The simulations suggest that the conformational behavior of PEO in ionic liquids is different from that in other common solvents, and computational studies of long chains will be interesting. sensitive to the identity of the anions32 and displays an unexpected dependence on molecular weight and concentration. Lee and Lodge32 showed that PEO is soluble in 1-ethyl3-methylimidazolium bis(trifluoromethylsulfonyl)imide ([EMIM][TFSA]) and 1-butyl-3-methylimidazolium hexafluorophosphate ([BMIM][PF6]) over a wide range of temperatures but displays a lower critical solution temperature (LCST) in 1butyl-3-methylimidazolium tetrafluoroborate [BMIM][BF4] and in 1-ethyl-3-methylimidazolium tetrafluoroborate [EMIM][BF4]. The phase diagram is essentially independent of molecular weight, and the critical concentration occurs at approximately 80 wt % of PEO. This is very different from phase diagrams of aqueous polymer solutions where the critical temperature is at very low polymer concentrations and displays a strong molecular weight dependence. A systematic study of polymers in ionic liquids will provide physical insight into their phase behavior and other interesting properties. The starting point would be the properties of single chains, or dilute solutions, but experiments on such systems are difficult. Very few single molecule experiments on polymers are available,33,34 and we are aware of only two experiments addressing the conformations of PEO in ionic liquids.27,35 Triolo et al.27 studied the conformational properties of deuterated PEO (Mw = 27.3 kDa) in [BMIM][BF4] as a

I. INTRODUCTION Ionic liquids have generated considerable excitement for their varied potential applications,1,2 as solvents for synthesis and catalysis,3,4 as electrolytes,5−7 as media for separations,8,9 as sorption media for gas absorption, and as lubricants.10 Ionic liquids (ILs) are usually composed of a large organic cation and a small organic or inorganic anion and are liquids at room temperature. They possess a number of interesting and important physical properties such as low volatility, nonflammability, high conductivity, and thermal and chemical stability.11 The fundamental interest and technological importance of ionic liquids have resulted in an exponential growth in studies devoted to understanding their physical chemistry.12−14 There has been recent interest in composites of ionic liquids with polymers because polymeric materials can provide mechanical integrity that ionic liquids otherwise lack.15 Possible applications of these composites include membranes for fuel cells, separations, and batteries, gels for artificial muscles, and dielectrics for energy storage.15−19 There have been a number of recent experiments that focus on the fundamental physical chemistry16,20−27 but little theoretical work.28−31 In this paper we report a simulation study of an atomistic model of poly(ethylene oxide), PEO, in water and ionic liquid and compare to the properties of the polymer in aqueous solution. There are many interesting open questions regarding the behavior of PEO in ionic liquids. Experiments have shown, for example, that the phase behavior of PEO in ionic liquids is very © 2013 American Chemical Society

Received: August 9, 2013 Revised: November 19, 2013 Published: December 23, 2013 438

dx.doi.org/10.1021/ma4016714 | Macromolecules 2014, 47, 438−446

Macromolecules

Article

Figure 1. Chemical formula of (a) PEO, (b) ionic liquid, i.e., [BMIM][BF4], and (c) BMIM with atoms are labeled differently.

model was employed, the simulations were performed at higher temperature (400 K), and the dependence of the properties on the degree of polymerization was not investigated. In this work we use an atomistic model of polymer and solvent, perform simulations at room temperature, and investigate the conformational properties as a function of the molecular weight for N ranging from 9 to 40. We focus on three aspects: (1) the solubility of the monomeric subunit of PEO,42 i.e., dimethoxyethane (DME) (chemical formula: CH3−O−CH2−CH2−O−CH3), (2) the conformational properties of PEO, and (3) the free energy of association of two PEO molecules, all in [BMIM][BF4]. For the purposes of comparison we also perform simulations in water at the same temperature. We find that DME is more soluble in [BMIM][BF4] than in water with a solvation free energy of −22 kJ/mol in [BMIM][BF4] compared to −14 kJ/mol in water. Water is close to a θ solvent for PEO; the root-meansquare radius of gyration, Rg, scales as Rg ∼ N0.48 over the range of N investigated. The PEO chains are swollen in [BMIM][BF4] with Rg ∼ N0.9, which reflects the increased stiffness of the chains in this solvent. The number of statistical segments, even for N = 40, is therefore quite low, and the scaling regime is probably not reached in our simulations. In fact, the distribution of the mean-square end-to-end distance is bimodal for PEO in [BMIM][BF4] with both extended and hairpin conformations. The free energy of association of two PEO molecules is more unfavorable in [BMIM][BF4] than in water due to high positive energies of association. The rest of this paper is organized as follows: the simulation model and methods are described in section II, results are presented in section III, and some conclusions are presented in section IV.

function of concentration using small-angle neutron scattering (SANS) and found that the radius of gyration, Rg, scaled with concentration, c, as Rg ∼ c−0.25 in the semidilute regime. This is not consistent with the expected scaling for chains in a good solvent, which is Rg ∼ c−0.125. In their own words, they had “no explanations for such a behavior”. Werzer et al.35 studied the conformational properties of PEO (Mw = 38 kDa) in ethylammonium nitrate using SANS and found that Rg ∼ c−0.25 in the dilute regime and Rg ∼ c−1 in the semidilute regime, behavior that cannot be explained using standard scaling theory of polymer solutions. The scaling of the size of an isolated chain, i.e., a chain in dilute solution, with degree of polymerization, N, can be inferred from the scaling of chain size with concentration in the semidilute regime. If R0 is the size of an isolated chain and R0 ∼ Nν, then the overlap threshold concentration, c*, scales as c* ∼ N/R03 ∼ N1−3ν. According to the Flory ideality hypothesis, chains are ideal (Rg ∼ √N) in semidilute solutions. Using the scaling hypothesis ⎛ c ⎞γ R g ∼ R 0⎜ ⎟ ∼ ⎝ c* ⎠

N

(1)

gives a relationship between the scaling exponents ν and γ, i.e.

ν + (3ν − 1)γ =

1 2

(2)

For chains in a good solvent, ν ≈ 3/5 and γ ≈ −1/8. For polyelectrolytes, ν = 1 and γ = −1/4. The Triolo et al.27 experiments (γ = −0.25) are therefore consistent with polyelectrolyte behavior, i.e., R0 ∼ N, which is surprising because PEO chains are not charged. The conformational properties of PEO in an ionic liquid is therefore an open question. In this work we study the conformational properties of isolated PEO molecules in an ionic liquid, [BMIM][BF4] (chemical formula shown in Figure 1). Although there are many simulation studies of neat ionic liquids36−39 and PEO in water,40−43 we are aware of only one simulation of the properties of PEO in an ionic liquid,30 where a united atom

II. SIMULATION MODEL AND METHODS The OPLSAA force field44−49 is used for PEO and [BMIM][BF4], and the TIP4P model is used for water.50 We use the OPLSAA force field because parameters for both PEO and [BMIM][BF4] are available. It is difficult to test a force field for polymers because experiments in dilute solution are generally 439

dx.doi.org/10.1021/ma4016714 | Macromolecules 2014, 47, 438−446

Macromolecules

Article

not available. The OPLS force field has been extensively used for 64 different ionic liquids. Finally, the force fields for PEO and the ionic liquid have been systematically developed under the OPLS framework and are therefore compatible. All simulations are carried out using the GROMACS 4.5.451 software package. There are three separate sets of simulations for (1) solvation free energy of a monomer, (2) conformational properties of the polymer, and (3) free energy of association of a pair of polymers, and the simulation protocol for these is described separately below. In all simulations, the Lennard-Jones interaction potential is cut off at 1.4 nm and shifted so that both the energy and force are zero. The particle-mesh Ewald (PME) summation52,53 is used for the long-range electrostatic interactions. The PME parameters are as follows: real space cutoff distance of 1.4 nm and interpolation order of 6, with a maximum fast Fourier transform grid spacing of 0.12 nm. All bonds involving hydrogen are fixed at their equilibrium distance using LINCS constraint algorithm. 54 For simulation with water, the SETTLE55 algorithm is used to keep the water molecules rigid. The system is evolved using a leapfrog algorithm with a time step of 2 fs. Solvation Free Energy of Dimethoxyethane. The solvation Helmholtz free energy, ΔF, is calculated via thermodynamic integration56,57 ΔF =

∫0

1

∂Uss(λ) ∂λ

equilibration run in the NPT ensemble at 300 K and 1 bar (temperature and pressure are maintained by using a Nosé− Hoover thermostat59,60 and Parrinello−Rahman barostat61,62). The box length after equilibration is L = 4.18, 5.01, and 5.6 nm for N = 9, 18, and 27, respectively, while for N = 40, the equilibrated box dimension is 6.78 nm × 4.19 nm × 4.19 nm. The number of solvent water molecules used to solvate the cubic box is 2401, 4071, 5761, and 3819 for N = 9, 18, 7, and 40, respectively. In each case properties are averaged over three independent 10 ns trajectories in the NVT ensemble at 300 K. Statistical uncertainties are obtained by block averaging over 5 ns blocks. For PEO in [BMIM][BF4], standard molecular dynamics simulations do not sample configurational space adequately. A visualization of the trajectories shows that the polymer gets trapped in some conformations. (We had encountered similar issues for polystyrenesulfonate in water.) Therefore, with the ionic liquid solvent, we perform replica exchange simulations.63−66 Initial configurations for the replica exchange simulations are generated as follows. A PEO chain in a randomly chosen conformation is first solvated by inserting first the BMIM cations and then an equal number of anions at randomly chosen positions in a box that is larger than the chain dimensions. The number of resulting solvent ion pairs, Ns, is 179, 251, 306, and 356 for N = 9, 18, 27, and 40, respectively. A cubic simulation cell, with linear dimension L, is used in all cases except N = 40. For N = 40 we use a rectangular parallelepiped box. In all cases the system is equilibrated for 3 ns at 450 K in the NPT ensemble. The final configuration from this simulation is used as the initial configuration for the replica exchange simulations in the NVT ensemble. The equilibrium box length is L = 3.92, 4.38, and 4.7 nm for N = 9, 18, and 27, respectively, while for N = 40, the equilibrated box dimension is 6.71 nm × 4.19 nm × 4.19 nm. For the replica exchange simulations, 64 replicas are used which are exponentially distributed between 300 and 550 K, with the temperatures chosen using the REMD calculator67 implemented on the Gromacs Web server. For each chain length, each of 64 replicas is run for 10 ns in NVT ensemble using a Nosé−Hoover thermostat,59,60 thereby spanning a total simulation length of 640 ns, with a replica exchange frequency of 1 ps. Each simulation is repeated for three independent velocity seeds . The acceptance ratio of replica exchanges ranges from 45 to 50% for N = 9 to 32−40% for N = 40. We find significant overlap in the energy distributions of adjacent replicas, which indicates successful sampling using replica exchange simulations. Conformational properties are obtained by dividing the trajectories into multiple 5 ns blocks and then averaging properties over those blocks. Umbrella Sampling Simulations for Polymer Conformations. The free energy of a polymer chain (N = 18) as a function of the end-to-end distance is studied using umbrella sampling simulations at 300 K. The reaction coordinate, R2, is the distance between the two ends of the polymer molecule, which are constrained to a given distance using an umbrella potential. We employ 58 windows for R2 ranging from 0.65 to 3 nm in steps of 0.05 nm. The force constant in the umbrella potential is set to 10 000 kJ/(mol nm2). A 6 ns simulation is performed for each window. The weighted histogram analysis method (WHAM)68,69 is used to unbias the histogram resulting from the umbrella sampling, and the free energy landscape is obtained as a function of R2.

dλ λ

(3)

where Uss(λ) is the interaction potential and interpolates between λ = 0 (vacuum) and λ = 1 (with solvent), and the brackets denote an ensemble average for a particular value of λ. Bonded interactions are interpolated linearly, and a soft-core approach58 is used for nonbonded interactions in order to remove singularities that arise from inserting sites. We calculate the solvation free energy of dimethoxyethane (DME), which has the chemical formula CH3−O−CH2−CH2− O−CH3. A single DME molecule is first solvated with solvent in a cubic box of linear dimension 3.2 nm. A steepest descent energy minimization and a short 500 ps equilibration run in the NPT ensemble is followed by 14 ns stochastic dynamics simulation at 300 K, also in the NVT ensemble. We use stochastic dynamics in this case because the solute−solvent interaction are scaled and stochastic dynamics can better represent the reduced interaction in the form of viscous drag and random Gaussian force. The simulations are performed for values of λ ranging from 0 to 1 in intervals of 0.1. For each value of λ, the free-energy derivative is computed at an interval of 10 ps and then block-averaged over the entire 14 ns trajectory. The integral in eq 3 is evaluated numerically. The difference between the solvation free energy in [BMIM][BF4] (ΔFIL) and water (ΔFW) gives the transfer free energy. The procedure is repeated at 285 K in order to compute temperature derivatives to obtain the entropic and energetic contributions to the free energy. Single Chain Simulations. Equilibrium Simulations and Replica Exchange Simulations. Standard molecular dynamics simulations are used for PEO in water at 300 K. For all chain lengths except N = 40 a cubic simulation cell is used with dimensions larger than the chain dimensions. For N = 40 the simulation cell is a rectangular parallelepiped. A randomly generated conformation is inserted into the box and solvated with water molecules. The simulations are subjected to a 1 ns 440

dx.doi.org/10.1021/ma4016714 | Macromolecules 2014, 47, 438−446

Macromolecules

Article

Potential of Mean Force of Association of Two PEO Chains. We calculate the potential of mean force between two PEO chains, for N = 9, using an umbrella sampling procedure described elsewhere.70 The initial configuration of the chain is obtained from the simulations of a single chain, which is then duplicated and solvated in a box of solvent with linear dimension L = 6 nm. The reaction coordinate (ξ) is the distance between the oxygen atoms of fifth subunit of the two molecules. We employ 21 windows ranging from 0.4 to 2.4 nm at a separation of 0.1 nm. The force constant for the umbrella sampling is 1000 kJ/(mol nm2). For ξ > 1.4 nm the two chains are kept in a side-by-side parallel configuration initially. For shorter distances the final configuration from the simulation at a larger separation is used, and the chains are pulled to the desired separation. The initial configuration is subjected to energy minimization using the steepest-descent technique followed by short equilibration in NPT ensemble at 300 K temperature and 1 bar pressure in the presence of the umbrella potential. Averages are calculated over a production run of 15 ns. The weighted histogram analysis method (WHAM)68,69 is used to unbias the histogram. The entropy is calculated from the finite difference temperature derivative of the PMF or ΔF(ξ) at each intersolute separation ξ, i.e. −ΔS(ξ) =

ΔF(ξ , T + ΔT ) − ΔF(ξ , T ) ΔT

Figure 2. Distribution of two different average dihedral angles obtained from simulation of PEO18 chain in [BMIM][BF4]. This distribution is similar irrespective of chain length.

The distribution of C−C−O−C dihedral angle is also periodic in nature with two most prominent peaks at −180° and 180° accompanied by two small peaks at −60° and 60°. PEO is more expanded in [BMIM][BF4] than in water. Simulation results for the radius of gyration, Rg, defined as Rg =

(4)

1 N

∑ (ri − rCM)2 (5)

i

where ri is the position of site i on a chain and rCM is the position of the chain center-of-mass, are presented in Table 1 and Figure 3.

where F is the Helmholtz free energy, ξ is the reaction coordinate, T is the temperature, and ΔT is the temperature difference between two simulations. We perform simulations for T = 285 and 300 K, and ΔT = 15 K. The energetic contribution is obtained from ΔU(ξ) = ΔF(ξ) + TΔS(ξ).

Table 1. Simulations Results for the Radius of Gyration, Rg Rg (nm)

III. RESULTS AND DISCUSSION The simulation result for the solvation free energy, ΔFIL, of DME in [BMIM][BF4] is −22.86 ± 0.57 kJ/mol, and that in water, ΔFW, is −14.35 ± 0.10 kJ/mol. The experimental value71,72 for ΔFW is −18 kJ/mol. The simulation results therefore overestimate the solvation free energy of DME in water. We are not aware of experimental results for ΔFIL. The simulations suggest that DME is more soluble in [BMIM][BF4] than in water, with a transfer free energy (from water to [BMIM][BF4]) of −8 kJ/mol. In both [BMIM][BF4] and water, the energy of solvation and the entropy of solvation are both negative; i.e., the process is energetically favorable but entropically unfavorable: ΔUIL = −102 kJ/mol, TΔSIL = −76 kJ/mol, ΔUW = −80 kJ/mol, and TΔSW = −62 kJ/mol. DME therefore has a very favorable energetic interaction with the solvent for both [BMIM][BF4] and water, but more so for the ionic liquid. This carries an entropic penalty because of the ordering of the solvent molecules, which is also more significant for the ionic liquid when compared to water. The scattering from the [BMIM][BF4] shows strong peaks, which are signatures of local ordering, as has been previously reported earlier by Siqueira and Ribeiro73 and Borodin and Smith.74 The symmetric nature of the distribution of dihedral angles gives us confidence in the conformational equilibration of the polymer molecules. The distribution (not shown) of C−O and C−C bond length shows a single sharp peak at about 0.15 nm. The distributions of C−O−C and O−C−C angles (not shown) also show a single sharp peak at about 112°. Figure 2 depicts the distribution dihedral angle O−C−C−O, which is periodic and has two symmetric prominent peak at about −60° and 60°.

a

N

Mwa

9 18 27 40

398 794 1190 1762

in water 0.524 0.746 0.933 1.089

± ± ± ±

0.007 0.021 0.046 0.056

in [BMIM][BF4] 0.570 0.864 1.210 1.728

± ± ± ±

0.009 0.021 0.033 0.022

Mw is the molecular weight.

Ideal chain statistics are observed in water, but PEO is significantly more expanded in the [BMIM][BF4]. In water, we obtain a scaling exponent ν, where Rg ∼ Nν of ν = 0.48 ± 0.07, which is consistent with the values of ν = 0.515 obtained previously using a different (CHARMM) force field and longer

Figure 3. Logarithmic plot of simulation results for the variation of Rg with N in water (filled circles) and [BMIM][BF4] (filled squares). Lines are fits to the three highest values of N and give scaling exponents of 0.48 ± 0.07 and 0.87 ± 0.1 in water and [BMIM][BF4], respectively. 441

dx.doi.org/10.1021/ma4016714 | Macromolecules 2014, 47, 438−446

Macromolecules

Article

chains.42 In [BMIM][BF4] we obtain a scaling exponent of ν = 0.87 ± 0.01, which is significantly higher. (Similar results are obtained if we include the data for N = 9 or if we analyze the RMS end-to-end distance.) The scaling exponent ν is universal and one expects only four values: 1/3, 1/2, 0.588, and 1 for collapsed, ideal, self-avoiding walk, and rod-like chain statistics. An exponent of 0.87 suggests implies the chains are not long enough in the [BMIM][BF4] to be in the true scaling regime. The simulations are not, however, definitive on the value of the asymptotic N → ∞ exponent. The behavior of the single chain structure factor, S(q), is consistent with the scaling of chain size with molecular weight. In the scaling regime, S(r) ∼ m/r3 where r is the size of a correlation volume with a monomer at the origin and m is the number of monomers within this volume. If r ∼ mν, then S(r) ∼ r1/ν−3 and S(q) ∼ q−1/ν. For an ideal chain, ν = 1/2, S(q) ∼ q−2 in the scaling regime, and a Kratky plot of q2S(q) vs q is flat in the scaling regime. If ν > 1/2, then q2S(q) ∼ q2−1/ν is an increasing function of q in the scaling regime. Figure 4 depicts a

extended as well as hairpin-like conformations in [BMIM][BF4].

Figure 6. Representative snapshots of hairpin (left) and extended conformations (right) for N = 18 in [BMIM][BF4] as obtained from umbrella sampling simulations. These snapshots closely resemble those obtained from REMD simulations.

A bimodal distribution of chain conformations is also reflected in the free energy (or potential of mean force) as a function of the end-to-end distance, depicted in Figure 7. Two

Figure 4. Kratky plot of PEO in water and [BMIM][BF4] for N = 18.

Kratky plot for PEO in [BMIM][BF4] and water for N = 18. In water, the behavior of the structure factor is similar to that of an ideal chain, but in [BMIM][BF4] the chain is more expanded. Of course, for these short chains there is not a well-defined scaling regime. The wiggles in S(q) for PEO in [BMIM][BF4] are due to insufficient sampling and are more pronounced for longer chains. The simulations provide a picture of the conformations of PEO in [BMIM][BF4]. The distribution of chain sizes is bimodal in the [BMIM][BF4] solvent. Figure 5 depicts the distribution P(R), of the end-to-end distance R2, for N = 18 and 40. In both cases, there are two peaks in P(R), suggesting that there are coexisting ensembles of configurations. An examination of snapshots (see Figure 6) confirms that there are

Figure 7. Free energy landscape as a function of end-to-end distance for N = 18 in [BMIM][BF4] obtained from umbrella sampling simulations at 300 K.

minima are observed for R2 ≈ 1.2 and 2.3 nm, which correspond to the peaks in the distribution of R2 in the REMD simulations (see Figure 5a). An examination of representative snapshots at R2 corresponding to the two free energy minima confirms the presence of hairpin and extended conformations. The agreement in the results obtained from REMD and umbrella sampling simulations strongly supports the existence of two subpopulations (hairpin and extended) of conformations for PEO in [BMIM][BF4].

Figure 5. Distribution of end-to-end distance for (a) N = 18 and (b) N = 40 in [BMIM][BF4]. 442

dx.doi.org/10.1021/ma4016714 | Macromolecules 2014, 47, 438−446

Macromolecules

Article

and various atoms on the ionic liquid. The atoms are labeled differently, as depicted in Figure 1c. Similar results are seen for other chain lengths. The positions of the peaks in the pair correlation functions are sensitive to the size of the various sites, and the heights of the peaks tell us about the correlation (or local density) between the sites. In all cases the correlations are not very strong; i.e., g(r) is never much greater than 1, with g(r) ≈ 1.6 for the highest peak. There are also no long-ranged correlations between the polymer and the ionic liquid cation. The polymer chain is correlated with the methyl group and side of the imidazole ring that has the carbon atom bonded to two nitrogen atoms. This signifies a strong energetic preference for the solvent cations to orient so that the butyl tails are away from the polymer molecules. This can be seen from Figure 10, where for both the oxygen and carbon atoms of the PEO, the peak in g(r) is most prominent with the methyl carbon, the C1 carbon on the imidazole ring, and the H1 atom (attached to the C1 atom). This depicts that the PEO chains are most strongly correlated with the CM and C1 region of the cation. This observation of a prominent correlation between the C1 region and PEO are consistent with results of a recent simulation study75 of PEO in a different ionic liquid. A pictorial view of a representative snapshot of a PEO chain (Figure 11) supports the conclusions inferred from g(r). In the first solvation shell of the polymer the cations reside with the methyl groups in close proximity to the polymer surface and the butyl groups pointed away. Figure 12 depicts the pair correlation functions between the O and C atoms of the PEO and the atoms of the anion. There are two peaks in the g(r) involving F atoms because of the geometry of the BF4 molecule. Other than that, the pair correlation functions reflect the size of the anion, and there is no preferential solvation of cation or anion at the polymer surface.

The distribution of end-to-end distance for PEO in water shows only one peak. Figure 8 depicts the end-to-end

Figure 8. Distribution of end-to-end distance for N = 18 in pure water.

distribution function of PEO in water for N = 18 and shows a single (albeit broad) peak centered at approximately 1.5 nm. We do not see any evidence of hairpin conformations in snapshots taken from these simulations. These two features suggest that the conformations of PEO in [BMIM][BF4] are quite different from those in water. The conformational properties suggest that [BMIM][BF4] is a very good solvent for PEO while water is a θ or marginal solvent. An examination of the potential of mean force between two PEO molecules confirms this conclusion; i.e., the PMF between PEO molecules is repulsive in [BMIM][BF4] and of very small magnitude (although repulsive) in water. Figure 9a depicts the PMF between two PEO chains (with N = 9) in [BMIM][BF4 ] and water. The entropic and energetic contributions to the PMF are shown in Figure 9b. In [BMIM][BF4] the PMF between PEO molecules is strongly repulsive and fairly long-ranged. In water, the PMF has a low magnitude, less than 0.1 kcal/mol at separations greater than a few angstroms. In both [BMIM][BF4] and water, the energetic contribution is unfavorable (positive) and the entropic contribution is favorable. One would expect PEO to be completely miscible in water and [BMIM][BF4] at this temperature, which is consistent with recent high-energy Xray diffraction experiment and molecular dynamics simulation.75 The pair correlation functions provide insight into which part of the solvent is correlated with different atoms on the PEO chain for N = 18. Figure 10 depicts the pair correlation functions between the oxygen and carbon atoms of the PEO

IV. SUMMARY AND CONCLUSIONS We present atomistic simulations for the behavior of isolated PEO chains in [BMIM][BF4] and compare to the properties of these oligomers in water. The monomer is significantly more soluble in [BMIM][BF4] than water with a larger negative free energy of solvation. The solvation of the monomer is energetically favorable with an entropic penalty. PEO is significantly more expanded in [BMIM][BF4] than in water. For the range of chain lengths studied, Rg ∼ N0.87 in [BMIM][BF4] compared to Rg ∼ N0.5 in water. The

Figure 9. (a) Potential of mean force between two PEO chains (N = 9) in [BMIM][BF4] and water and (b) energetic (ΔU) and entropic (−TΔS) contributions to the PMF. The solid lines correspond to [BMIM][BF4], and the dashed lines correspond to water. 443

dx.doi.org/10.1021/ma4016714 | Macromolecules 2014, 47, 438−446

Macromolecules

Article

Figure 10. Comparison of radial distribution function of (a) heavy atoms of [BMIM] with O atoms of PEO, (b) heavy atoms of [BMIM]with C atoms of PEO, and (c) hydrogen atoms with O atoms of PEO, for N = 18.

this temperature is essentially a theta solvent, with ideal chain scaling observed over the chain lengths studied. The potential of mean force between two PEO molecules is repulsive and much stronger than that between two PEO molecules in water, where the PMF is essentially zero over all but the shortest separations. The repulsion comes from energetic interactions between the polymers and the solvent which overcomes the favorable entropy of association. The scaling of chain size with N we observe is unexpected. Since there are no long-ranged intramolecular interactions, one expects self-avoiding walk statistics for long chains. It is possible that these chains are too short to be in the asymptotic long chain limit. On the other hand, similar scaling behavior is inferred in experiments for the same system with much larger degree of polymerization (N ≈ 600),27 although in the experiments the variation of chain size with concentration (not degree of polymerization) was measured. It is tempting to think that the polymer might be associated with one of the ions of the solvent, thus effectively a polyelectrolyte in nature, but we see no evidence of strong association. An examination of the pair correlation functions indicates that the polymer is associated with the methyl part and the C1 carbon atom of the imidazole ring. If one wishes to alter the solvent quality, then this is the region of the cation that one would change. There has been a suggestion that [BMIM][BF4] is a heterogeneous solvent with hydrophobic and polar regions. If this is true, then PEO would preferentially segregate to the polar regions, i.e., away from the alkyl tails. If such a segregation did occur, it might explain the surprising conformational properties because the chain conformations might be slave to the underlying structure of the (heterogeneous) solvent. The simulations of PEO in [BMIM][BF4] presented here are restricted to fairly short chains because the dynamics are slow and replica exchange simulations are necessary. It would be

Figure 11. A representative snapshot of spatial distribution of BMIM cations (colored blue) within 0.4 nm of polymer surface of N = 18. Most of the cations are oriented with butyl groups away from the polymer surface.

Figure 12. Comparison of radial distribution function of O and C atoms of PEO, with anions for N = 18.

distribution of end-to-end distance is bimodal, with coexisting extended and hairpin-like conformations. In contrast, water at 444

dx.doi.org/10.1021/ma4016714 | Macromolecules 2014, 47, 438−446

Macromolecules

Article

(18) Ueki, T.; Watanabe, M. Macromolecules in ionic liquids: Progress, challenges, and opportunities. Macromolecules 2008, 41, 3739−3749. (19) Gray, F. Polymer Electrolytes; RSC Materials Monographs: Cambridge, 1997. (20) Rodriguez, H.; Rogers, R. Liquid mixtures of ionic liquids and polymers as solvent systems. Fluid Phase Equilib. 2010, 294, 7−14. (21) Harner, J.; Hoagland, D. Thermoreversible gelation of an ionic liquid by crystallization of a dissolved polymer. J. Phys. Chem. B 2010, 114, 3411−3418. (22) Rodriguez, H.; Francisco, M.; Rahman, M.; Sun, N.; Rogers, R. Biphasic liquid mixtures of ionic liquids and polyethylene glycols. Phys. Chem. Chem. Phys. 2009, 11, 10916−10922. (23) Sarkar, A.; Trivedi, S.; Pandey, S. Polymer molecular weightdependent unusual fluorescence probe behavior within 1-butyl-3methylimidazolium hexafluorophosphate + poly(ethylene glycol). J. Phys. Chem. B 2009, 113, 7606−7614. (24) Sarkar, A.; Trivedi, S.; Baker, G.; Pandey, S. J. Phys. Chem. B 2008, 112, 14927−14936. (25) Tsuda, R.; Kodama, K.; Ueki, T.; Kokubo, H.; Imabayashi, S.; Watanabe, M. Chem. Commun. 2008, 40, 4939−4941. (26) Ueki, T.; Karino, T.; Kobayashi, Y.; Shibayama, M.; Watanabe, M. Difference in lower critical solution temperature behavior between random copolymers and a homopolymer having solvatophilic and solvatophobic structures in an ionic liquid. J. Phys. Chem. B 2007, 111, 4750−4754. (27) Triolo, A.; Russina, O.; Keiderling, U.; Kohlbrecher, J. Morphology of poly(ethylene oxide) dissolved in a room temperature ionic liquid: A small angle neutron scattering study. J. Phys. Chem. B 2006, 110, 1513−1515. (28) Kirchner, B. Ionic Liq. 2009, 290, 213−262. (29) Costa, L.; Ribeiro, M. Molecular dynamics simulation of polymer electrolytes based on poly(ethylene oxide) and ionic liquids. II. Dynamical properties. J. Chem. Phys. 2007, 127, 164901. (30) Costa, L.; Ribeiro, M. Molecular dynamics simulation of polymer electrolytes based on poly(ethylene oxide) and ionic liquids. I. Structural properties. J. Chem. Phys. 2006, 124, 184902. (31) White, R. P.; Lipson, J. E. G. Macromolecules 2013, 46, 5714− 5723. (32) Lee, H.; Lodge, T. Lower critical solution temperature (LCST) phase behavior of poly(ethylene oxide) in ionic liquids. J. Phys. Chem. Lett. 2010, 1, 1962−1966. (33) Wang, D.; Yuan, Y.; Mardiyati, Y.; Bubeck, C.; Koynov, K. From single chains to aggregates, how conjugated polymers behave in dilute solutions. Macromolecules 2013, 46, 6217−6224. (34) M, V.; Habuchi, S. Conformation and physics of polymer chains: a single-molecule perspective. NPG Asia Mater. 2010, 2, 131− 142. (35) Werzer, O.; Warr, G. G.; Atkin, R. Conformation of poly(ethylene oxide) dissolved in ethylammonium nitrate. J. Phys. Chem. B 2011, 115, 648−652. (36) Kowsari, M. H.; Alavi, S.; Ashrafizaadeh, M.; Najafi, B. Molecular dynamics simulation of imidazolium-based ionic liquids. I. Dynamics and diffusion coefficient. J. Chem. Phys. 2008, 129, 224508. (37) Ludwig, R.; Maginn, E.; Balasubramanian, S. Editorial: Ionic liquids: The fundamentals and forces driving their rise. ChemPhysChem 2012, 13, 1603−1603. (38) Rai, N.; Maginn, E. J. Vapor-liquid coexistence and critical behavior of ionic liquids via molecular simulations. J. Phys. Chem. Lett. 2011, 2, 1439−1443. (39) Castner, E.; Margulis, C. J.; Mark, M.; Wishart, J. F. Ionic liquids: Structure and photochemical reactions. Annu. Rev. Phys. Chem. 2011, 62, 85−105. (40) Smith, G. D.; Bedrov, D.; Borodin, O. Conformations and chain dimensions of poly(ethylene oxide) in aqueous solution: A molecular dynamics simulation study. J. Am. Chem. Soc. 2000, 122, 9548−9549. (41) Smith, G. D.; Bedrov, D.; Borodin, O. Molecular dynamics simulation study of hydrogen bonding in aqueous poly(ethylene oxide) solutions. Phys. Rev. Lett. 2000, 85, 5583−5586.

interesting to study this system for longer chains, a variety of ionic liquid solvents, and for a range of temperatures. It is probably necessary to develop a coarse-grained force field to achieve that goal. Coarse-grained models for PEO in water and for pure [BMIM][BF4] are available, and developing such a force field for PEO in [BMIM][BF4] is a promising direction. It is important also to validate the force fields by comparison to experiment of, for example, the transfer free energy of the monomer from water to the ionic liquid. These are interesting future directions that we hope are motivated by the interesting conformational properties observed in this work.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail [email protected] (A.Y.). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This material is based upon work supported by the National Science Foundation under Grant CHE-1111835. We are grateful for generous computational support from trestles and gordon machine in San Diego Supercomputer Center (SDSC) under Grant TG-CHE090065 and the UW Madison Chemistry Department Cluster Phoenix under Grant CHE-0840494.

■

REFERENCES

(1) Plechkova, N.; Seddon, K. Applications of ionic liquids in the chemical industry. Chem. Soc. Rev. 2008, 37, 123−150. (2) Wassercheid, P. Volatile times for ionic liquids. Nature 2006, 439, 797−797. (3) Zhang, Z. Catalysis in Ionic Liquids. Adv. Catal. 2006, 49, 153− 237. (4) Welton, T. Ionic liquids in catalysis. Coord. Chem. Rev. 2004, 248, 2459−2477. (5) Armand, M.; Endres, F.; MacFarlane, D.; Ohno, H.; Scrosati, B. Ionic-liquid materials for the electrochemical challenges of the future. Nat. Mater. 2009, 8, 621−629. (6) Galinski, M.; Lewandowski, A.; Stepniak, I. Ionic liquids as electrolytes. Electrochim. Acta 2006, 51, 5567−5580. (7) Buzzeo, M.; Evans, R.; Compton, R. ChemPhysChem 2004, 5, 1106−1120. (8) Han, X.; Armstrong, D. Ionic liquids in separations. Acc. Chem. Res. 2007, 40, 1079. (9) Smiglak, M.; Metlen, A.; Rogers, R. The second evolution of ionic liquids: from solvents and separations to advanced materials-energetic examples from the ionic liquid cookbook. Acc. Chem. Res. 2007, 40, 1182−1192. (10) Zhu, F.; Liang, Y.; Liu, W. Ionic liquid lubricants: designed chemistry for engineering applications. Chem. Soc. Rev. 2009, 38, 2590−2599. (11) Huddleston, J.; Visser, A.; Reichert, W.; Willauer, H.; Broker, G.; Rogers, R. Green Chem. 2001, 3, 156−164. (12) Clare, B.; Sirwardana, A.; MacFarlane, D. R. Ionic Liq. 2009, 290, 1−40. (13) Gomes, M.; Lopes, J.; Padua, A. Ionic Liq. 2009, 290, 161−183. (14) Zhao, H.; Xia, S.; Ma, P. J. Chem. Technol. Biotechnol. 2005, 80, 1089−1096. (15) Lodge, T. A. Unique platform for materials design. Science 2008, 321, 50−51. (16) Winterton, N. Solubilization of polymers by ionic liquids. J. Mater. Chem. 2006, 16, 4281−4293. (17) Kubisa, P. Ionic liquids as solvents for polymerization processes - Progress and challenges. Prog. Polym. Sci. 2009, 34, 1333−1347. 445

dx.doi.org/10.1021/ma4016714 | Macromolecules 2014, 47, 438−446

Macromolecules

Article

(42) Lee, H.; Venable, R.; MacKerrel, A., Jr.; Pastor, R. Molecular dynamics studies of polyethylene oxide and polyethylene glycol: hydrodynamic radius and shape anisotropy. Biophys. J. 2008, 95, 1590−1599. (43) Lee, H.; deVries, A.; Marrink, S.; Pastor, R. A coarse-grained model for polyethylene oxide and polyethylene glycol: Conformation and hydrodynamics. J. Phys. Chem. B 2009, 113, 13186−13194. (44) McDonald, N. A.; Jorgensen, W. L. Development of an all-atom force field for heterocycles. properties of liquid pyrrole, furan, diazoles, and oxazoles. J. Phys. Chem. B 1998, 102, 8049−8059. (45) Briggs, J. M.; Matsui, T.; Jorgensen, W. L. Monte Carlo simulations of liquid alkyl ethers with the OPLS potential functions. J. Comput. Chem. 1990, 11, 958−971. (46) Chandrasekhar, J.; Spellmeyer, D. C.; Jorgensen, W. L. Energy component analysis for dilute aqueous solutions of lithium(1+), sodium(1+), fluoride(1−), and chloride(1−) ions. J. Am. Chem. Soc. 1984, 106, 903−910. (47) Jorgensen, W. L.; Maxwell, D.; Tirado-Rives, J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (48) Sambasivarao, S. V.; Acevedo, O. Development of OPLS-AA force field parameters for 68 unique ionic liquids. J. Chem. Theor. Chem. 2009, 5, 1038−1050. (49) Ren, C.-L.; Tian, W.-D.; Szliefer, I.; Ma, Y.-q. Specific salt effects on poly(ethylene oxide) electrolyte solutions. Macromolecules 2011, 44, 1719−1727. (50) Jorgensen, W.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926. (51) Hess, B.; Kutzner, C.; Van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. J. Chem. Theory Comput. 2008, 4, 435− 447. (52) Darden, T.; York, D.; Pederson, L. G. Particle mesh Ewald: An Nlog(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089. (53) Essman, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pederson, L. G. A smooth particle mesh Ewald method. J. Chem. Phys. 1995, 103, 8577. (54) Hess, B.; Bekker, H.; Berendsen, H.; Fraaije, J. G. E. LINCS: A linear constraint solver for molecular simulations. J. Comput. Chem. 1997, 18, 1463−1472. (55) Miyamoto, S.; Kollman, P. Settle: An analytical version of the SHAKE and RATTLE algorithm for rigid water models. J. Comput. Chem. 1992, 13, 952−962. (56) van Gunsteren, W. F.; Berendsen, H. J. C. Thermodynamic cycle integration by computer simulation as a tool for obtaining free energy differences in molecular chemistry. J. Comput.-Aided Mol. Des. 1987, 1, 171−176. (57) Martini coarse-grained force field: Extension to carbohydrates. J. Chem. Theor. Chem. 2009, 5, 3195−3210. (58) Beutler, T. C.; Mark, A. E.; van Schaik, R. C.; Gerber, P. R.; van Gunsteren, W. F. Avoiding singularities and numerical instabilities in free energy calculations based on molecular simulations. Chem. Phys. Lett. 1994, 222, 529−539. (59) Nosé, S. A molecular dynamics method for simulations in the canonical ensemble. Mol. Phys. 1984, 52, 255−268. (60) Hoover, W. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A 1985, 31, 1695−1697. (61) Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 52, 7182. (62) Nosé, S.; Klein, M. L. Constant pressure molecular dynamics for molecular systems. Mol. Phys. 1983, 50, 1055−1076. (63) Hukushima, K.; Nemoto, K. Exchange Monte Carlo method and application to spin glass simulations. J. Phys. Soc. Jpn. 1996, 65, 1604− 1608.

(64) Sugita, Y.; Okamoto, Y. Replica-exchange molecular dynamics method for protein folding. Chem. Phys. Lett. 1999, 314, 141−151. (65) Seibert, M.; Patriksson, A.; Hess, B.; van der Spoel, D. Reproducible polypeptide folding and structure prediction using molecular dynamics simulations. J. Mol. Biol. 2005, 354, 173−183. (66) Seibert, M.; van der Spoel, D. Protein folding kinetics and thermodynamics from atomistic simulations. Phys. Rev. Lett. 2006, 96, 238102−238105. (67) Patriksson, A.; van der Spoel, D. A temperature predictor for parallel tempering simulations. Phys. Chem. Chem. Phys. 2008, 10, 2073−2077. (68) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. The weighted histogram analysis method for freeenergy calculations on biomolecules. I. The method. J. Comput. Chem. 1992, 13, 1011−1021. (69) Grossfield, A. WHAM: the weighted histogram analysis method, version 2.0; http://membrane.urmc.rochester.edu/content/wham. (70) Mondal, J.; Yethiraj, A. Driving force for the association of amphiphilic molecule. J. Phys. Chem. Lett. 2011, 2, 2391−2395. (71) Kelly, C. P.; Cremer, C. J.; Truhlar, D. G. SM6:A density functional theory continuum solvation model for calculating aqueous solvation free energies of neutrals, ions, and solute and water clusters. J. Chem. Theory Comput. 2005, 1, 1133−1152. (72) Vorobyov, I.; Anisimov, V. M.; Greene, S.; Venable, R.; Moser, A.; Pastor, R. W.; Mackerrel, A. D. Additive and classical Drude polarizable force fields for linear and cyclic ethers. J. Chem. Theory Comput. 2007, 3, 1120−1133. (73) Siqueira, L. J. A.; Ribeiro, M. C. C. Charge ordering and intermediate range order in ammonium ionic liquids. J. Chem. Phys. 2011, 135, 204506. (74) Borodin, O.; Smith, G. S. Structure and dynamics of N-methylN-propylpyrrolidinium bis(trifluoromethanesulfonyl)imide ionic liquid from molecular dynamics simulations. J. Phys. Chem. B 2006, 110, 11481−11490. (75) Asai, H.; Fujii, K.; Nishi, K.; Sakai, T.; Ohara, K.; Umebayashi, Y.; Shibayama, M. Solvation structure of poly(ethylene glycol) in ionic liquids studied by high-energy X-ray diffraction and molecular dynamics simulations. Macromolecules 2013, 46, 2369−2375.

446

dx.doi.org/10.1021/ma4016714 | Macromolecules 2014, 47, 438−446