Molecular Dynamics Modeling of CO2 and Poly(ethylene glycol) in

Sep 10, 2013 - Center for Computational Natural Sciences and Bioinformatics, International Institute of Information Technology, Gachibowli, Hyderabad ...
1 downloads 17 Views 3MB Size
Download PDF
Article pubs.acs.org/JPCC

Molecular Dynamics Modeling of CO2 and Poly(ethylene glycol) in Montmorillonite: The Structure of Clay−Polymer Composites and the Incorporation of CO2 Marimuthu Krishnan,*,†,‡ Moumita Saharay,‡ and R. James Kirkpatrick§ †

Center for Computational Natural Sciences and Bioinformatics, International Institute of Information Technology, Gachibowli, Hyderabad 500 032, India ‡ Department of Chemistry and §College of Natural Science, Michigan State University, East Lansing, Michigan 48824, United States S Supporting Information *

ABSTRACT: Composite materials composed of aluminosilicate clays with organic molecules or biomolecules in the interlayer galleries are readily synthesized and have many applications in agriculture, medicine, environmental science, engineering, and geochemistry. Detailed characterization of the molecular-scale structure and dynamics of the interlayer galleries is difficult experimentally because of static and dynamic disorder but can be obtained by molecular dynamics (MD) simulation using classical force fields. This Article presents an MD study of smectite clay montmorillonite (MMT) intercalated with poly(ethylene glycol) (PEG) with and without CO2 also present in the interlayer. Bulk PEG can incorporate large amounts of CO2 under supercritical conditions, and MMT−PEG composites have potential as gas-separation materials. The MD simulations of anhydrous MMT containing interlayer CO2 (MMT−CO2), PEG (MMT−PEG), and PEG + CO2 (MMT−PEG−CO2) provide new atomic-level insight into the molecular ordering of CO2 near the basal surface and CO2-induced changes in the structure, conformation, and energetics of PEG in the interlayer. The results show that the structural arrangement among the CO2 molecules in the MMT− CO2 system is similar to that in supercritical CO2 and is analogous to that of crystalline CO2. All Na ions in this system remain coordinated by the basal oxygens (Ob) but are also coordinated by the oxygens (OCO2) of the CO2 molecules, and a few are displaced ∼2 Å above their surface sites. The cooperative motion of the Na ions and CO2 molecules increases the Na diffusion and the CO2 reorientation in the interlayers. The critical interactions among cations, CO2, PEG, and the basal surface in the MMT−PEG−CO2 system result in a layer of CO2 trapped between the basal surface and cation-permeable PEG film formed on the basal surface. In the MMT−PEG and MMT−PEG−CO2 systems, PEG is highly disordered and exhibits conformational and orientational heterogeneity in the interlayer confinement with many of the molecules lying in a layer parallel and close to the basal surface of the clay. Some of the PEG molecules span from one basal surface to the other across the interlayer. Many of the Na ions are displaced from the basal surface and are coordinated by both basal oxygens and oxygens of the PEG. Some are completely displaced from the surface and are dissolved in the PEG. In the MMT−PEG−CO2 system, the CO2 molecules occur in well-defined layers parallel to the basal surface and are most commonly dissolved in the PEG. In the full MMT−PEG−CO2 system, Na ions have few OCO2 nearest neighbors, and the presence of CO2 causes the basal spacing to expand but does not change the conformation of the PEG. The MD results provide detailed and otherwise unobtainable information about the coordination environments and distributions of interatomic distances and angles in the interlayer.

■

INTRODUCTION Clay minerals are abundant synthetic and naturally occurring (alumino)silicate phases that readily form complexes with organic compounds, polymers, and biomolecules.1−7 Understanding the interactions between the organic and inorganic components of these complexes is important for the development of applications in fields as diverse as agriculture, geochemistry, environmental science and engineering, and medicine.8,9 Clay minerals have layered structures, and swelling clays such as montmorillonite (MMT) have expandable interlayers that can accommodate a wide range of species © 2013 American Chemical Society

including cations, water molecules, CO2, and many polymers, organic molecules, and biomolecules. Many stable clay−organic complexes are known.1−9 Obtaining detailed molecular-scale insight into the structure and dynamics in the interlayer galleries is difficult experimentally because of static and dynamic disorder, and computational approaches using molecular dynamics (MD) can be highly Received: May 30, 2013 Revised: September 2, 2013 Published: September 10, 2013 20592

dx.doi.org/10.1021/jp405321t | J. Phys. Chem. C 2013, 117, 20592−20609

The Journal of Physical Chemistry C

Article

electronegative ether oxygen that is readily polarized. The methylene groups are hydrophobic, whereas the ether oxygen is hydrophilic. This amphiphilic character leads to the flexibility and viscoeleastic properties of bulk PEG, to its ability to coordinate a wide range of ions and molecules, and thus to the range of potential applications of clay−PEG composites.17,56,60−62 PEG has a net dipole moment and a helical conformation in the bulk crystalline phase due to the intermolecular and intramolecular dipole−dipole interactions between the polar ether oxygens and van der Waals interactions between the CH2 groups. The C−C− O−C and C−O−C−C dihedral angles have a trans conformation, and the O−C−C−O dihedral angles have a gauche conformation.63 Upon hydration, PEG retains its helical conformation, and the ether oxygens can accept hydrogen bonds (H bonds) from the water molecules.64,65 In the molten state, PEG adopts a random coil conformation with a significant number of C−O dihedrals populating the gauche states.66−68 In MMT interlayers, we expect the ether oxygens of PEG to coordinate to the charge-balancing cations through ion−dipole interactions and also to the sites on the basal oxygen surface of the T−O−T layers to which the oxygens of water molecules coordinate. Our MD results presented here show both to be true. The interplay between MMT−PEG and PEG−PEG interactions is also expected to influence the PEG conformation and packing. The MMT−PEG interactions dominate near the clay surface, whereas the intermolecular interactions between PEG chains dominate near the center of the interlayer. Although various authors have proposed that PEG intercalated into Na−MMT adopts a helical conformation coiling around the interlayer cations, a distorted helical conformation, or a planar conformation,56−59 a complete consensus has not yet been reached on the precise conformation of PEG confined in clay interlayers. Our results show that the PEG conformations are highly disordered with short-range helical order near Na ions. CO2 is a linear, nondipolar molecule with a nonzero quadrupole moment (4.1 × 10−26 esu69) under ambient conditions. The oxygens of CO2 are partially electronegative relative to the carbon atom and thus interact favorably with electron-deficient atoms of neighboring molecules. In the supercritical phase (scCO2), CO2 interacts strongly with polar polymers because of the enhancement of its multipole moments at elevated pressures.70 This enhancement is due to the inhomogeneous distribution of near-neighbor CO2 molecules around the central CO2 molecule. The relative orientations of nearest-neighbor pairs of CO2 molecules change with increasing density from the gas phase to the condensed phases. For instance, the “slipped” parallel geometry is relatively more stable than the T-shaped geometry for an isolated CO2 dimer, whereas the molecules prefer T-shaped configurations in the higher-order clusters and condensed phases.71−73 In the supercritical state, the nearest neighbors (NN) and next-nearest neighbors (NNN) lie preferentially in the equatorial plane and polar regions of the central CO2, respectively, to optimize eletron donor−acceptor interactions between the central and neighboring molecules.74,75 At very high pressures, the nearest-neighbor coordination is reminiscent of the structure in crystalline CO2.74,75 Experimental phase equilibrium studies of binary PEG−CO2 mixtures have shown an increase in PEG solubility at elevated pressures in scCO2, whereas at temperatures and pressures below the critical point of CO2, 305 K and 7.4 MPa, PEG interacts weakly with CO2.76,77 Supercritical CO2 decreases the melting point of PEG78 and causes an 8−10-fold decrease in its viscosity.79 FTIR experiments show changes in the stretching

informative. There have been few MD studies of clay−organic complexes.5,6,10−18 This Article presents an MD study of MMT intercalated with poly(ethylene glycol) (PEG) with and without CO2 also present in the interlayer. Bulk PEG can incorporate large amounts of CO2 under supercritical conditions, and MMT−PEG composites have potential as gas separation materials as well as in applications in the storage and delivery of drugs and in oil field drilling operations.19−27 The results provide a detailed picture of the structural interactions among the interlayer species and the clay layers. The well-known T−O−T structures of the (alumino)silicate layers of swelling clays consist of a layer of oxygen octahedra (O) sandwiched between two layers of oxygen tetrahedra (T). The octahedra typically contain Al3+, Mg2+, and other cations with charges of +1, +2, and +3, and the tetrahedra typically contain Si4+, Al3+, and occasionally other cations. The total structural charge on the layer is negative and is controlled by the charge and site occupancy of the octahedral cations and the extent of Al3+ for Si4+ substitution in the tetrahedral layers.2,5 This negative structural charge is balanced by the positive charge of interlayer and surface species. There are many methods for intercalating polymers, organic species, and biomolecules into clay interlayer galleries from solution or melts.1,28−32 These complexes often have mechanical and thermal properties significantly different from those of their individual components. The presence of the organic species can also greatly affect the affinity, interlayer structural environments and transport properties of small neutral species such as H2O and CO2.10,14,33−36 The interlayer spacings of clay−organic complexes are controlled by the complex interactions among the T− O−T layers and the interlayer species.12,37−40 These spacings can be readily determined by X-ray and neutron diffraction techniques, but these methods rarely provide detailed information about the interlayer gallery structure.7,29,36,41−44 Spectroscopic methods such as NMR and IR can provide insight into local structural environments in the galleries,7,25,45−49 but MD methods are often needed to provide a more detailed understanding of many structural, energy, and dynamic features of the interactions among the T−O−T layers and the interlayer species.5,6,11,14,23 The characterization of clay−polymer composites is especially difficult because of the large number of conformational degrees of freedom of the polymers. Clay−polymer composites have a range of potential industrial applications in nuclear waste disposal and gas storage and separation . There is, thus, significant interest in understanding the structure and energetics of clay composites containing polymers that interact positively with CO2. PEG is an environmentally benign polymer that shows high selectivity for CO2.19,50−53 It is used to make CO2-selective membranes and is a major component of the solvents that remove CO2 from flue gases.19,50−53 CO2 is well known to cause PEG to swell.54,55 Clay−PEG composites thus have attracted significant interest in recent years, but there is considerable uncertainty about the molecular-scale details of the conformation, orientation, and packing of PEG molecules and their interaction with the T−O− T layers and the structural and dynamical behavior of the charge balancing cations and CO2 molecules.7,44,47,56−59 PEG is also a useful model species for understanding the role of ether oxygens in the interaction of clays with organic species in the natural environment. PEG is a relatively simple polymer with a structural formula of HO−(CH2CH2O)n−OH. The monomeric unit, CH2CH2O, consists of two concomitant methylene groups attached to an 20593

dx.doi.org/10.1021/jp405321t | J. Phys. Chem. C 2013, 117, 20592−20609

The Journal of Physical Chemistry C

Article

without CO2 present. The initial structures of these systems were generated by removing all PEGs from the MMT−PEG−CO2 system except one per interlayer or removing all PEGs except one per interlayer and all CO2 molecules. The MMT−CO2 and MMT−PEG−CO2 systems contain 4.17 CO2/Na, comparable to the expected CO2 coordination around Na in clay confinement.14 The MMT−PEG and MMT− PEG−CO2 systems contain 1.04 11-monomer PEG molecules/ Na, comparable to experimentally observed composites.7,44 The MD simulations use the CLAYFF force field13 for the clay substrate, the CHARMM force field91,92 for the PEG, and the new force field parameters of Cygan et al.14 for CO2. The partial charges on the atoms of MMT and PEG are given in Table S1 (Supporting Information). The simulations were carried out using NAMD 2.9.93 They were performed in the NPT ensemble at a temperature of 318 K and a pressure of 8 MPa, just above the critical point of bulk CO2 (305 K, 7.4 MPa), using a Langevin thermostat and barostat with a damping coefficient of 5 ps−1. The nonbonded pair interaction potentials were truncated at 10 Å and smoothed between 8 and 10 Å using a cubic switching function. Periodic boundary conditions were applied. Electrostatic interactions were computed using the particle mesh Ewald (PME) method94 with a real space cutoff of 10 Å. The reciprocal space interactions were computed on 42 × 72 × 52 Å3 (MMT/ PEG and MMT−PEG−CO2 systems) and 42 × 72 × 27 Å3 (MMT−CO2 system) grids using sixth-degree B splines. The equations of motion were integrated with a time step of 1 fs. After initial energy minimization, the systems were equilibrated for 3 ns followed by a 6 ns production run. The interlayer structures of the modeled systems were visualized using the VMD software95 and analyzed with atomic probability density plots perpendicular to the clay layers (Zdensity plots), projections of the positions of the atoms within defined mean Z levels onto the clay basal surface (atomic projection maps), radial distribution functions between atoms, the distribution of C−C and C−O dihedral angles of PEG, the distribution of the angle between the basal clay surface and the molecular axis of the CO2 molecules, and the probability densities of the distributions of the O atoms of CO2 (OCO2) around each other and Na ions and CCO2 around the ether oxygens of PEG (OPEG). To characterize the lengths of individual PEG molecules, we define two conformation lengths: the Z length and the XY length. Let xmin, xmax, ymin, ymax, zmin, and zmax be the minimum and maximum values of the positions of the C and O atoms of an individual polymer chain in the x, y, and z directions. These values define a bounding box around the molecule with lengths xmax − xmin, ymax − ymin, and zmax − zmin. We define Z length = zmax − zmin, quantifying the length of the molecule perpendicular to the surface, and XY length = ((xmax − xmin)2 + (ymax − ymin)2)1/2, quantifying the length of a chain parallel to the surface. To calculate the 3D atomic density distributions around PEG, a fixed frame of reference was chosen such that the C−O−C unit of PEG lies on the XY plane, one of the O−C bond vectors is oriented along the X-axis, and the ether oxygen is at the origin. At each time step, this coordinate transformation was performed for both the C−O−C unit and its environment. The density maps were averaged over all snapshots in the MD trajectory and all C− O−C units of all PEG chains. To determine the distribution of CO2 around a central CO2 molecule, a fixed frame of reference was chosen such that the carbon of the central CO2 is at the

and bending modes of scCO2 resulting from the interactions between PEG and scCO2.80 There is also a significant decrease in the interfacial tension between CO2 and PEG at pressures above the critical pressure of CO2.81 The solubility of PEG in scCO2 decreases with increasing temperature and increasing PEG molecular weight.76,77 The molecular-scale theoretical understanding of these experimental observations is incomplete. It is known that CO2 can be incorporated into clay interlayers at elevated pressure,14,33−36,82,83 but the molecular-scale understanding of its behavior is just emerging and much remains to be learned about its structure, dynamics, and energetics, the interactions between CO2 and the T−O−T layers and chargebalancing cations, and the effects of pressure and temperature. Although previous experimental studies have used scCO2 as a medium for the organomodification of MMT to produce polymer−clay composites,84−89 to our knowledge there have been no reported studies of the structure, dynamics, and energetics of clay-PEG−CO2 composites. This Article presents a computational MD study of MMT− containing interlayer CO2, PEG, and PEG + CO2 that provides significant new insight into the structure and energetics of the interactions among these species, charge-balancing Na+, and the T−O−T layers.

■

SIMULATION DETAILS Classical MD simulations were carried out using a model system consisting of 128 crystallographic unit cells of MMT (8a × 8b × 2c) with the structural formula Na0.75(Al3.5Mg0.5) [Si7.75Al0.25]O20(OH)4. In this model, structural charge development occurs because of both octahedral and tetrahedral substitution. The sites of octahedral Mg for Al substitution and tetrahedral Al for Si substitution were chosen randomly but were evaluated individually to ensure that the substitution sites are separated from each other by at least one nonsubstituted site. This model system contains two interlayer galleries each containing 48 Na ions. In the initial structure, the Na ions (24/ basal surface) were placed above the centers of the hexagonal rings of (Si,Al)-O tetrahedra 2 Å from the basal surface. The models of the base MMT, MMT−CO2, MMT−PEG, and MMT−PEG−CO2 composites were generated by suitably expanding the interlayers of the dry Na−MMT (up to 13 Å) and randomly packing the desired molecules in the interlayer regions.90 Each PEG molecule consisted of 11 CH2−O−CH2 units with the two ends capped by CH2−OH groups. Their molecular weight was 546.6 g/mol. The initial structure of the PEG molecule was generated by extracting an 11-mer segment from a protein−PEG complex (PDB ID 3V4P) and capping the two ends with CH2−OH units. All PEG chains in the initial structures of MMT−PEG and MMT−PEG−CO2 systems had identical starting conformations but were translated and rotated randomly to pack the interlayer galleries optimally with few or no hard contacts between neighboring chains.90 In most of the models presented here, the initial positions of the Na ions were close to the basal surfaces. To evaluate whether the systems reached an equilibrium configuration during the simulations, we ran an additional MMT−PEG−CO2 system in which the Na ions were initially located randomly throughout the interlayer. The results for the two types of systems are essentially identical. To examine the interactions among MMT, PEG, and CO2 in the absence of PEG−PEG intermolecular interactions, additional simulations were carried out for systems containing a single PEG molecule intercalated in each of the two interlayers with and 20594

dx.doi.org/10.1021/jp405321t | J. Phys. Chem. C 2013, 117, 20592−20609

The Journal of Physical Chemistry C

Article

Figure 1. Z-density plots. (a) Base MMT system, (b) MMT−CO2 system, (c) MMT−PEG system, (d) MMT−PEG−CO2 system, (e) MMT−PEG− CO2 system with an initially random distribution of Na, and (f) MMT−PEG−CO2 system with only 1 PEG molecule/interlayer.

Figure 2. XY atomic density maps. Blue lines represent the Al,Si−O atoms on the basal surface of the clay. Atomic densities shown are for those atoms that occur between the basal surface and the first minimum in the Z-density profile away from the basal surface for that species. (a) Base MMT (red = Na), (b) MMT−CO2 (red = Na, yellow = Si, blue = Al, white = Ob, green = OCO2), (c) MMT−PEG (red = Na, yellow = OPEG), (d) MMT−PEG−CO2 system (red = Na, green = OCO2, yellow = OPEG, white = OPEG of representative PEG molecules with all OPEG on the basal surface), (e) MMT−PEG− CO2 system with one PEG molecule/interlayer (red = Na, green = OCO2, yellow = OPEG next to the basal surface, and purple = CPEG next to the basal surface). The same XY-density maps viewed through the basal surface are shown in Figure S1 (Supporting Information).

Na occurs in two layers, one near each basal surface (Z = −0.72 and 0.6 Å). In all of the Z-density plots presented here, Z = 0 denotes the center of the interlayer. The presence of the Na probability density at Z = 0 indicates that there is a nonzero probability for the Na ions to hop between the two basal surfaces, as also observed for smectite Na-hectorite.15 The distributions of the atoms of the T−O−T layers in this simulation and all others

origin and one of the C−O bond vectors is oriented along the z axis.

■

RESULTS AND DISCUSSION

Base MMT. The Z-density plot for the base Na−MMT system with no PEG or CO2 (Figure 1a) shows that the basal d spacing is 9.4 ± 0.02 Å, which is consistent with earlier simulations (9.4 Å)13 and experimental (9.5−9.8 Å)96,97 results. 20595

dx.doi.org/10.1021/jp405321t | J. Phys. Chem. C 2013, 117, 20592−20609

The Journal of Physical Chemistry C

Article

Figure 3. Radial distribution functions (solid lines) and running coordination numbers (dashed lines) in the base MMT and MMT−CO2 systems. (a) Na−O of base MMT, (b) Na−O of MMT−CO2, all Na, (c) Na−O of MMT−CO2, Na in layer closest to the basal surface, (d) Na−O of MMT−CO2, Na displaced from basal surface, (e) OCO2−Ob, MMT−CO2, all OCO2 (black), only OCO2 of CO2 in the shoulder at Z = ±1.34 Å (red), (f) CCO2−Ob, MMT−CO2, all CCO2 (black), only CCO2 of CO2 with at least one OCO2 in the shoulder at Z = ±1.34 Å.

their surface positions in the absence of CO2 and are coordinated to OCO2. A visual inspection of the atomic trajectories reveals that all Na ions shuttle between positions on the basal surface (large peak) and closer to the CO2 layer (small peak) during the course of the simulation, demonstrating a dynamic origin of the computed bimodal Z distribution. The residence time of the Na ions is longer on the basal surface than in the layer closer to CO2, showing that the Z = ±1.35 Å state is energetically less stable than the Z = ±2.55 Å state and that the interaction of Na+ with the basal surface is stronger than with CO2. The potential of mean force, F(Z), associated with a given Z state of Na can be determined using the Boltzmann inversion formula: F(Z) = −kBT ln P(Z), where P(Z) denotes a normalized Z-density profile, kB is the Boltzmann constant, and T is the temperature. Using this formulation, we find that the Z = ±1.35 Å state is 0.65 kcal/mol less stable than the Z = ±2.55 Å state and that the activation energy for Na to hop from the surface site to the CO2 layer is ∼1.5 kcal/mol. The lower-bound estimate of the freeenergy barrier for Na to move across the CO2 layer is ∼7 kcal/ mol. Because the affinity of water molecules for Na ions is greater than that for CO2, it is most probable that the water molecules will displace the CO2 molecules from the coordination shells of Na ions in the hydrated MMT-CO2 systems. Consequently, the shoulders in the OCO2 distribution are expected to disappear as observed in a recent study on a one-layer hydrate of MMT− CO2.98 Similarly, the shoulders in the OCO2 distribution at Z = ±1.35 Å are also due to dynamic effects. No CO2 molecule is observed to be oriented permanently at a high angle with respect to the basal surface. The trajectories also reveal a dynamic correlation between the hopping of Na ions between the two positions and the reorientation of CO2 molecules from parallel to tilted orientations. This correlation causes Na ions that have moved

are narrow and well-defined, demonstrating the stability of the aluminosilicate MMT structure during the simulations. The Na ions in each layer of this system lie above the center of the hexagonal rings of the Si,Al−O tetrahedral of the nearest basal surface and execute small-amplitude vibrations about the center of these rings (Figure 2a). Because of the composition of the clay, the Na ions do not occupy every ring. The XY maps for all systems shown in Figure 2 are for one basal surface and show the positions of the atoms of a given type located between the plane of basal oxygens (Ob) and the first minimum in the corresponding Z distribution. In the base MMT system, Na has a total oxygen coordination of about 10 within an NN distance of 3.5 Å (Figure 3a). Six of these are Ob from the basal surface with which the Na is associated, 3 are Ob from the opposite basal surface, and 1.0 is the hydroxyl oxygen (OOH) below the hexagonal ring with which it is associated. MMT−CO2 System. For the MMT−CO2 system, the interlayer spacing is 12.9 ± 0.06 Å, an expansion sufficient to accommodate one layer of CO2 molecules (Figure 1b). A recent simulation study of hydrated MMT−CO2 has shown d spacings of ∼12.7 and 12.2 Å for a one-layer hydrate with two and three CO2 molecules per unit cell, respectively.14,98 The Z-density distribution for CCO2 has a relatively broad peak centered at Z = 0, showing that the centers of the CO2 molecules are at the midplane of the interlayer. The OCO2 Z densities show a broad peak centered at Z = 0 that is split into two subpeaks at Z = ±0.25 Å and shoulders near Z = ±1.3 Å. The Na ions in this system occur primarily in two layers adjacent to the two basal surfaces (Z = ±2.55 Å), comparable to their location in the base MMT system. A few Na ions are displaced from the surface (peaks at Z = ±1.35 Å) with Z levels that overlap with the shoulders in the OCO2 distribution. These Na ions are displaced about 1.3 Å from 20596

dx.doi.org/10.1021/jp405321t | J. Phys. Chem. C 2013, 117, 20592−20609

The Journal of Physical Chemistry C

Article

Figure 4. Distribution of angles in the NN coordination environments of Na displaced from the basal surface in the MMT−CO2 system. (a) The angle θ between the surface normal and the Na−OCO2 vector, (b) the angle ϕ between two Na−OCO2 vectors of the same first-neighbor coordination shell, (c) the angle η between the surface normal and the OCO2−CCO2 vector. All NN OCO2 (black), only NN OCO2 of CO2 in the shoulder at Z = ±1.34 Å (red), and NN OCO2 not in the shoulder at Z = ±1.34 Å (blue).

peak is due to the coordination of Na ions by all other NN CO2 molecules. Most of the Na ions that are displaced from the basal surface are coordinated by 4 OCO2, and a few are coordinated by 3 or 5 (Figure S2a, Supporting Information). The Na ions on the basal surface are coordinated predominantly by either 2 or 3 OCO2’s (Figure S2b). The distribution of the angle, ϕ, between the Na−O vectors of two NN O atoms (OCO2 and/or OPEG) of a near-surface Na ion is a measure of the shape of the first coordination shell around Na. For Na coordination by OCO2 in the MMT−CO2 system, the P(ϕ) distribution for all Na atoms has two maxima at ϕ = 71 and 126° (Figure 4b). The magnitude of the maximum at ϕ = 71° is approximately 3 times greater than that at ϕ = 126°. The distribution for just those Na in the layer displaced from the basal surface has a maximum at about 126°, demonstrating that this maximum is due to CO2 molecules that are tilted toward the basal surface, and the maximum at ϕ = 71° is due to those CO2 molecules that have their molecular axes almost parallel to the basal plane. CO2 is a linear molecule, and its orientation in the interlayer is an important structural parameter. For the MMT−CO2 system, the OCO2 and CCO2 Z-density distributions suggest that at any instant most of the CO2 molecules are oriented parallel or nearly parallel to the T−O−T layers but some are at a high angle to it. The distribution of the angle, ψ, between the molecular axis of CO2 and the normal to the clay basal surface shows this to be true (Figure 5). The mean ψ value is 90° (e.g., the mean orientation of

from a basal surface position to return to the basal surface at a different hexagonal ring, thus facilitating Na diffusion parallel to the clay surfaces. No Na ion is observed to move across the CO2 layer during the simulation. This observation is in line with a recent MD simulation study of a hydrated MMT−CO2 system, which demonstrated that an increase in the concentration of CO2 prohibits the migration of the Na ions from the basal surface to the center of the interlayer.98 Many of the Na ions in the MMT−CO2 system occupy the same types of positions as in the base Na−MMT, but a few others show vibrational displacement, especially toward tetrahedral Al sites, as shown by the elongated probability distributions (Figure 2b). The OCO2’s are located above the centers of the hexagonal rings, and CO2 packs near the basal surfaces, thus reflecting the periodicity of the underlying clay structure. This structure is reflected in the OCO2−Ob RDF and running coordination numbers, which show that OCO2 of the CO2 molecules in the Zdensity shoulder at Z = ±1.34 Å have six NN Ob’s and welldefined periodicity at longer distances, reflecting the structure of the basal surface (Figure 3e). The RDFs and running coordination numbers for all OCO2’s and CCO2’s are much less structured and are similar to those of a homogeneous system (Figure 3e,f). In the MMT−CO2 system, Na on average has a coordination of about nine within a nearest-neighbor (NN) distance of 3.5 Å with 6 Ob and 3 OCO2 (Figure 3b). For those Na ions in the layer closest to the basal surface (Figure 1b), the total coordination is about 9.5 within a NN distance of 3.5 Å. Of these, 6 are Ob from the basal surface with which the Na is associated, 1.0 is the hydroxyl oxygen under the six-membered ring, and about 2.5 are OCO2 (Figure 3c). This coordination is consistent with the typical Na+ being located above the center of a 6-membered tetrahedral ring with approximately 3 nearby CO2 molecules at any instant. For Na ions in the layer displaced from the surface, the NN coordination is 7.5 with 3.4 Ob and 4.1 OCO2 (Figure 3d). This coordination is consistent with Na being displaced from the center of the six-membered rings, as shown in Figure 2b. The distribution of the angle, θ, between the normal to the clay basal surface (Z axis) and the vector between a displaced Na ion and the coordinating NN OCO2 atoms is broad with a maximum at θ = 60° and a shoulder at θ = 82° (Figure 4a). The shoulder in this distribution is due to the coordination of Na that is displaced from the basal surface (small peak in Figure 1b) and coordinated by OCO2 at the same Z level (shoulder in Figure 1b). The main

Figure 5. Distribution of the ψ angles between the surface normal and the molecular axis of CO2 for the MMT−CO2 (black) and MMT− PEG−CO2 (red) systems. 20597

dx.doi.org/10.1021/jp405321t | J. Phys. Chem. C 2013, 117, 20592−20609

The Journal of Physical Chemistry C

Article

the CO2 axis is parallel to the surface), but the ψ angle distribution is relatively broad with the population between 77 and 103° being nearly flat. There are very few CO2 molecules oriented at a high angle to the basal surface with ψ values of less than 20° or greater than 160°. This distribution is reflected in the overlap of the OCO2 and CCO2 Z-density distributions for this system (Figure 1b). A similar distribution that peaked at around 90° was observed for the hydrated MMT−CO2 systems.98 The distribution of the angle, η, between the surface normal and the OCO2−CCO2 bond vector of the CO2 molecules coordinated to Na ions that are displaced from the surface further quantifies the tilting of CO2 molecules in the vicinity of the basal surfaces (Figure 4c). The OCO2 at the same Z level as the Na ions has a maximum at 65°, and the others have a maximum at 83°. The peak at η = 65° indicates that a few CO2 molecules proximal to the basal surface are tilted such that one of their OCO2 atoms is pointed toward the nearby surface and the other OCO2 points toward the center of the interlayer. The coordination among the CO2 molecules in the MMT− CO2 system is very similar to that in supercritical CO2.70 As for supercritical CO2, the probability density isosurfaces for the location of OCO2 molecules around a central CO2 show a belt at the level of the equatorial plane (blue in Figure 6) and hemispherical caps on either end of the molecule (red in Figure 6). As described by Saharay et al.,70 this coordination is also analogous to that in crystalline CO2, in which there are 6 OCO2’s in the equatorial plane and the O of the central CO2 is coordinated to 3 CCO2’s of the coordinating molecules with their

OCO2 located above the O of the central CO2. This coordination is possible in our MMT−CO2 system because not all of the OCO2’s lie on the same Z level and not all of the CO2 molecules have the same orientation, as discussed above. The probability density isosurfaces for the OCO2 that coordinate the Na ions in the MMT−CO2 system that are in the layer displaced from the basal surface show four lobes at about 90° to each other (Figure 7), consistent with the four OCO2 coordination in the running coordination numbers (Figures 3d and S2a). Although the nature of the Na−CO2 interactions and the locations and dynamics of the Na ions are sensitive to the degree of water content in the interlayer galleries,14,98 the anhydrous MMT−CO2 model used here characterizes the specific interactions between CO2 and MMT in the absence of water molecules. MMT−PEG System. The presence of interlayer PEG results in a substantial reorganization of the interlayer gallery relative to the base MMT system, including the displacement of many of the Na ions from the basal surface. The Z-density plot for the MMT−PEG system has a much richer structure that reflects the distribution and conformation of the PEG molecules and the coordination of the Na ions by OPEG as well as by Ob (Figure 1c). The basal d spacing is 23.7 ± 0.08 Å and is expanded substantially from the base system to accommodate the PEG. The Z-density peaks for OPEG are not symmetrical about Z = 0, reflecting the conformation of the molecules. There are OPEG peaks at Z = ±4.75, −1.35, and 1.0 Å (with a shoulder at Z = 1.9 Å). CPEG shows multiple peaks at Z = ±5.6, −1.4, and 1.8 Å with significant intensity at Z = 0 and between the main peaks. Na also shows multiple peaks at Z = ±8.35, ±6.9, and ±3.3 Å. These peaks represent, respectively, Na in positions similar to those in the base MMT system (called layer1 in the discussion below), those further from the surface and coordinated by both Ob atoms and OPEG (layer2), and those fully displaced from the surface but proximal to the PEG adlayers and coordinated by OPEG only (layer3). No Na+ was observed to migrate to the center of the interlayer during the course of the simulation. The affinity of OPEG for cations is well known and is the reason for the dissolution of alkali salts in the bulk polymer that makes it an excellent polymer electrolyte.17,30,31,38,99,100 The relative intensities of the OPEG and CPEG peaks show that many of the monomers lie close to the basal surface and that CPEG is typically closer to the surface than is OPEG, although the peaks overlap. Previous experimental studies have demonstrated that the intercalation of PEO into MMT increases the mobility of cations by weakening the cation−clay interactions.4,42,56,101 Past MD simulation studies have observed a dynamic exchange of the interlayer cations between the surface and the gallery center in the hydrated system whereas the Na ions in the dry MMT were observed to coordinate mainly by the PEO oxygens located proximal to the clay surfaces.16,57 However, the time scales of the previously published simulations are too short (