Large-Scale Molecular Dynamics Study of Montmorillonite Clay

8248

J. Phys. Chem. C 2007, 111, 8248-8259

Large-Scale Molecular Dynamics Study of Montmorillonite Clay: Emergence of Undulatory Fluctuations and Determination of Material Properties James L. Suter, Peter V. Coveney,* H. Chris Greenwell,‡ and Mary-Ann Thyveetil Centre for Computational Science, Department of Chemistry, UniVersity College London, 20 Gordon Street, London, WC1H 0AJ, United Kingdom ReceiVed: January 13, 2007; In Final Form: March 16, 2007

Large-scale atomistic simulations, which we define as containing more than 100 000 atoms, are becoming more commonplace as computational resources increase and efficient classical molecular dynamics algorithms are developed. With the advent of grid-computing methods, it is now possible to simulate even larger systems efficiently. Using this new technology, we have simulated montmorillonite clay systems containing up to approximately ten million atoms whose dimensions approach those of a realistic clay platelet. This considerably extends the spatial dimensions of microscopic simulation into a domain normally encountered in mesoscopic simulation. The simulations exhibit emergent behavior with increasing size, manifesting collective thermal motion of clay sheet atoms over lengths greater than 150 Å. This motion produces low-amplitude, longwavelength undulations of the clay sheets, implicitly inhibited by the small system sizes normally encountered in atomistic simulation. The thermal bending fluctuations allow us to calculate material properties, which are hard to obtain experimentally due to the small size of clay platelets. Montmorillonite is commonly used as a filler in clay-polymer nanocomposites, and estimation of the elastic properties of the composite requires accurate knowledge of the elastic moduli of the components. We estimate the bending modulus to be 1.6 × 10-17 J, corresponding to an in-plane Young’s modulus of 230 GPa. We encounter a clay sheet persistence length of approximately 1400 Å, which dampens the undulations at long wavelengths for the largest system in our study.

1. Introduction Smectite clay minerals are layered aluminosilicates which normally occur as small platelets with a thickness of the order of 1 nm and lateral dimensions of ca. 100-1000 nm. They consist of negatively charged quasi two-dimensional sheets kept together by charge-balancing interlayer cations giving a stacked crystalline structure.1,2 Montmorillonite is among the most common clay minerals and has received much attention recently due to the increased thermochemical and mechanical properties of composite materials comprised of polymers containing dispersed clay platelets.3-7 Montmorillonite is a 2:1 clay mineral, consisting of an octahedral alumina or magnesia layer sandwiched between two silica layers. The clay material is defined by the amount and location of isomorphic substitution in the clay sheets, which confers the overall negative charge. In humid conditions, the intersheet cations can become hydrated, forcing the clay sheets apart. The charge-balancing ions in naturally occurring clays tend to be small inorganic cations, such as Na+ or K+, but these can be exchanged for more complex species. This is facilitated by a characteristic property of clays, their variable intersheet spacing. Owing to their small platelet size and stacking behavior, clay minerals lack long-range order, and as a result, they are difficult to characterize accurately by experimental methods, such as powder X-ray diffraction. Some indication of the structure of the material can be readily obtained experimentally (e.g., the * To whom correspondence should be addressed. E-mail: [email protected].: +44 (0)20 7679 4560/4850.Fax: +44 (0)20 7679 4603. ‡ Present address: Centre for Applied Marine Science, School of Ocean Science, University of Wales Bangor, Mean Bridge, Anglesey, LL59 5AB, United Kingdom.

spacing between sheets), but ascertaining the intersheet molecular arrangement is much more difficult. Molecular simulation is, therefore, an ideal way of gaining insight at an atomistic level into the molecular structure and dynamics of clay minerals.7-10 However, the complexity of clay minerals requires large models to simulate realistic systems.11 This is commonly achieved by simulating a small system derived by replication of the original unit cell to form a supercell (typically on the order of 10 000 atoms for classical molecular dynamics) subjected to periodic boundary conditions in all three spatial directions during simulation. It is assumed that such an approach will replicate the properties of the extended condensed phase system. For clay minerals, the thickness is microscopic (nanometers), whereas the lateral dimensions of the sheets reside in the macroscopic domain (micrometers). A question seldom asked is whether the supercell is large enough for this assumption to be correct, that is, to sample the motions of the many different length scales inherent to clay minerals. In the motivating study of Greenwell et al.,12 large-scale simulations of intercalated clay-polymer nanocomposites revealed the existence of longwavelength, low-amplitude undulations in the clay sheets. The undulations were not observed in smaller size models; indeed, it is often assumed that clay sheets are rigid bodies.7 Analogous finite size effects have been reported previously in atomistic and mesoscopic simulations of biological13-16 and nonbiological membranes.17 The observation of clay sheet flexibility enabling thermally excited undulations is of significance as they allow the calculation of materials properties, such as the elastic modulus,18

10.1021/jp070294b CCC: $37.00 © 2007 American Chemical Society Published on Web 05/17/2007

Molecular Dynamics Study of Montmorillonite Clay required for theoretical understanding of clay-polymer nanocomposites, that are otherwise hard to obtain due to the small grain size of the clay minerals and the ease of reaction.19 Despite the current interest in the material properties of smectite clays, there remains a lack of data about their elastic constants. Understanding the materials properties of composites requires detailed knowledge about each component.20-22 The calculation of these properties from the clay sheet undulations is the main focus of the present paper. Through increasing computational power and more efficient algorithms, we can now study, in atomistic detail, clay systems large enough to be classified as being in the “mesoscopic realm”. In this article, we utilize the computational resources involved in high-performance grid computing to investigate a series of increasingly large molecular dynamics simulations of hydrated sodium montmorillonite, subject to periodic boundary conditions, ranging from 7000 to 10 million atoms. Through largescale atomistic-level simulations, we can calculate emerging macroscale material properties, thereby creating a bridge between atomistic- and continuum-level models. In this paper, we define “large-scale” systems as ones whose simulation cells contain more than 105 atoms. The largest system size studied here, containing approximately 107 atoms, approaches the dimensions of clay sheets and tactoids typically observed in electron micrographs for montmorillonite (∼200 nm).7 This considerably expands the spatial scales of atomistic aqueous mineral simulation. The structure of the paper is as follows. In the next section, we describe the components necessary to perform large-scale molecular dynamics, such as the scalability of our molecular dynamics code, the algorithms we employ, and the gridcomputing infrastructure. In section 3, we briefly summarize the theoretical background, relating fluctuations in clay sheet undulations and thickness to macroscopic quantities such as the bending modulus. In section 4, we investigate a number of properties of the clay system for each system size in order to determine the extent of finite size effects. The structure of the intersheet region can be determined from molecular simulation by calculating radial distribution functions (RDF) of intersheet species and atomic density profiles perpendicular to the clay surface. These reveal the coordination environment of atoms and ions, for example, the hydration state of the cations and hydrogen bond interactions between water molecules and the clay surface. These properties have previously been used to illustrate the balance between cation hydration and sheet-sheet repulsion in swelling clays.23 We describe how these properties change with system size. We also calculate the stress-strain relationship for the four smallest models in our study by increasing the tension in the system and measuring the change in internal pressure. Then, in section 4, we analyze the collective motions of the clay sheet atoms leading to undulatory motions and calculate continuum properties such as the bending modulus and Young’s modulus from these mesoscopic phenomena. Finally, we discuss the implications of our findings and present our conclusions. 2. Large-Scale Simulation Methods 2.1. Scalable Grid-Based Molecular Dynamics. All simulations were carried out with the LAMMPS2006 package, written in C++.24,25 The LAMMPS code, developed at Sandia National Laboratories, is a classical molecular dynamics code designed for simulation of atomic and molecular systems on parallel supercomputers. We have used specifically designed algorithms available in LAMMPS that reduce the wall clock time required

J. Phys. Chem. C, Vol. 111, No. 23, 2007 8249

Figure 1. Scaling properties of LAMMPS on the HPCx IBM eServer pSeries 690 (1.7 GHz Power4+) calculated using a hydrated montmorillonite clay system of 1 033 900 (diamonds) and 9 455 000 (squares) atoms. Speedup ) (time(N)/time(N1)) where N is the number of processors used in the simulation and N1 is the number of processors for the reference calculation (here, N1 ) 64).

for large-scale simulations. The largest system in our study utilized the reversible reference system propagator algorithm (rRESPA) multi-time scale integrator to determine atomic forces and velocities. The computationally cheaper intramolecular and short-range intermolecular interactions were computed in the inner loop, and the more expensive long-range electrostatics were computed in the outer loop. This allowed the long-range electrostatics to be calculated up to eight times less frequently without losing accuracy. In our study, we have used an inner time-step loop of 0.5 fs and a corresponding outer-loop time step of 4 fs. Periodic boundary conditions were imposed in all spatial directions. The expensive long-range electrostatics were calculated using Ewald summation combined with the efficient particle-particle particle-mesh (PPPM) solver, which scales as N log N, where N is the number of atoms. Both rRESPA and PPPM have been recently implemented in simulations of clays.11 LAMMPS exhibited a near-linear scaling relationship between the speed-up in simulation time and both the number of processors and the system size.11,26 The parallelization was achieved using MPI and spatial domain decomposition. The simulation box was decomposed into subregions, which were then distributed across many different processors. In Figure 1, we report the scaling performance for large-scale hydrated clay montmorillonite systems run on an increasing number of processors on the IBM eServer pSeries 690 (1.7 GHz Power4+) at HPCx, the U.K.’s national high-performance computing service. The 1 033 900 atom system scales almost linearly until 512 processors, after which the scaling decreases. For the 9 455 000 atom system, the system exhibits near-linear behavior even at the largest number of processors (1024). Therefore, to optimally exploit the scaling of LAMMPS, we should use 1024 processors for simulations containing ten million atoms and 512 processors for the million atom simulations. The large number of processors required was made possible by the capability computing service provided by HPCx, which allows the simulations to be run over the full production partition of 1024 processors for codes with near-linear scaling. In our study, approximately 60 000 and 18 000 CPU hours were used to simulate our 9 455 000 and 1 055 000 atom systems, respectively. The simulations of the smaller systems in our study were performed using the computing resources available on the lower

8250 J. Phys. Chem. C, Vol. 111, No. 23, 2007

Suter et al.

Figure 3. AtomEye visualization of the 263 750 atom sodium montmorillonite system III after approximately 1 ns of simulation time. The z direction has been expanded 15 times to assist viewing. The size of the simulation cell is 450 Å × 260 Å × 24.1 Å. The bending of the clay sheets and also the roughness of the clay surface are clearly visible. The black lines indicate the periodic boundaries. The arrows indicate the surface of the clay layer near some of the adsorbed sodium ions, where surface oxygen atoms protrude further into the interlayer spacing than the surrounding clay surface atoms of the simulation cell. The color coding is the same as in Figure 2. Figure 2. AtomEye-rendered snapshot of sodium montmorillonite in (a) the yz plane and (b) the xz plane. The system shown contains 6752 atoms, system I, composed of two aluminosilicate sheets. Each sheet is made up of two tetrahedral (T) and one octahedral layer (O). The interlayer region contains charge-balancing sodium ions and water molecules. The distances along the bottom of the snapshot illustrate the average dimensions of system I. The distances on the top of each snapshot are the unit cell dimensions, replicated to form system I. For (a), the z distance shown corresponds to the system dimension, and for (b), the z distance shown is the unit cell dimension. Atoms are colored as follows: Mg ) green, Al ) blue, O ) red, Si ) gray, Na ) yellow, H ) white.

end nodes of the U.K. National Grid Service, the US TeraGrid, and the EU Distributed European Infrastructure for Supercomputing Applications. The recently developed Application Hosting Environment (AHE) was used for job submission.27 The AHE is an interoperable lightweight hosting environment for running unmodified applications such as LAMMPS on Globus- and UNICORE-based grid resources. For the smallest system in our study, consisting of 6752 atoms, the simulations used 8 processors of a 70 processor Linux cluster within the Centre for Computational Science at University College London (UCL). Job submission and subsequent analysis for the smallest system were performed using WEDS (WSRF Environment for Distributed Simulation).28 WEDS allows deployment of pre-existing application codes across multiple resources within a single administrative domain. Subsequent analysis of structural and dynamic properties was distributed over local workstations using the work-flow capabilities of WEDS. WEDS and AHE are freely available and can be downloaded from the RealityGrid website (www.realitygrid.org and www.omii.ac.uk). Visualization is an important tool for examining and analyzing the structure and behavior of molecular systems. Our systems were visualized using the AtomEye atomistic configuration viewer.29 AtomEye is a freely available atomistic visualization software package for UNIX platforms, available in serial and parallel formats. AtomEye treats atoms as spheres and bonds as cylinders. These objects are rendered as primitive objects rather than composites of polygons. The largest system in our study was visualized on an SGI Prism at UCL using six processors equipped with 8 GB of memory. The smaller systems were visualized on a PC with 4 GB of memory. Example AtomEye-rendered images are shown in Figure 2, and animations are available in the Supporting Information. For systems greater than 250 000 atoms, we were more concerned with the visualization of the entire system rather than

Figure 4. AtomEye visualization of sodium montmorillonite system V, containing approximately 10 million atoms, after 0.5 ns of simulation. The z direction has been expanded 20 times to allow visualization of thermal sheet fluctuations. Color coding is the same as in Figure 2.

the motion of individual atoms. To emphasize this, we visualized the system with the z coordinates scaled by a factor of 10 or 15. Figures 3 and 4 are images of the 263 750 and 9 495 000 atom systems, respectively. Visualization of the atomic configurations of large-scale simulations allowed us to observe effects not present in smaller models, namely, the collective motion of the clay atoms to produce long-wavelength undulatory modes. 2.2. Force-Field and Simulation Details. The hydrated montmorillonite system simulated in this study uses the recently developed CLAYFF force field,30 which has been shown to be highly effective in modeling the structure of many oxide, hydroxide, and clay phases, as well as the interaction of aqueous solutions and their dissolved species with mineral surfaces.31-37 The interatomic potentials in CLAYFF are derived from parametrization incorporating structural and spectroscopic data for a variety of simple hydrated compounds. All atoms are represented as point charges and are allowed complete translational freedom, permitting simulation of complex disordered systems. Metal-oxygen interactions are based on a simple Lennard-Jones (6-12) potential combined with electrostatics. The empirical parameters are optimized on the basis of known mineral structures, and the partial atomic charges are derived from periodic density functional theory (DFT) calculations of simple oxide, hydroxide, and oxyhydroxide model compounds.30 The charges on the O and H atoms vary with nearest-neighbor cation substitution. Flexibility of the OH groups in the clay lattice is accomplished with harmonic bond stretch and angle bend terms. For water, the flexible version of the SPC (single point charge) model is used.38 The systems were simulated for at least 1ns at 300 K and 1 atm using an isobaric-isothermal

Molecular Dynamics Study of Montmorillonite Clay


TABLE 1: Simulation Cell Composition and Dimensions for the Simulated Clay-Water Systems system

no. of atoms

unit cell replication

no. of water

no. of Na+

simulation length ns

starting supercell dimensions/Å3

I II III IV V

6752 51 062 263 750 1 055 000 9 495 000

4×4×2 11 × 11 × 2 25 × 25 × 2 50 × 50 × 2 150 × 150 × 2

512 3872 20 000 80 000 720 000

96 726 3750 15 000 135 000

2.0 2.5 1.8 4.3 0.6

41.37 × 71.64 × 24.62 113.76 × 197.01 × 24.62 258.55 × 447.76 × 24.62 517.10 × 895.10 × 24.62 1551.29 × 2686.53 × 24.62

(NPT) ensemble in which the temperature and pressure were maintained via a Nose´-Hoover thermobarostat. The pressure was set separately in all three directions. This should lead to a tensionless state. The thermostat and barostat relaxation times were 100 and 1000 fs, respectively. The simulation box remained rectilinear in shape throughout the simulations. The cutoff distance for van der Waals interactions is 10 Å. We used the Verlet integration scheme with a time step of 1.0 fs. As mentioned above, the largest simulation and the stress-strain calculations use the rRESPA algorithm to economize CPU time, with an outer time step for electrostatics of 4 fs and an inner time step for bond terms of 0.5 fs. Angle and dihedral terms were calculated with time steps of 1 and 2 fs, respectively. 2.3. Model Structures. The model clay structure used was a Wyoming-like montmorillonite with unit cell formula Na3[Al14Mg2][Si31Al]O80(OH)16.30 Isotopic substitutions confer a negative charge on the aluminosilicate sheet framework; the octahedral layer contains 2 Mg2+ ions for every 14 Al3+ ions, while the tetrahedral layers contain 1 Al substitution for every 32 Si atoms. The net negative charge in the clay sheets is compensated by sodium counterions in the interlayer spacing. Rendered images of clay systems are shown in Figure 2. The hydrated sodium montmorillonite unit cell contains 16 water molecules, sufficient to form a single monolayer of water.39 In this study, we restricted our simulations to models with two intersheet regions in order to suppress artificial features due to periodic effects which arise in a single intersheet space. For example, two independent sheets can allow clay sheets to translate relative to each other, giving a more realistic system. Increasing the size of the z direction through the use of two intersheet spacings also reduces any self-interaction through the periodic boundary. The systems sizes were (I) 6752 atoms (i.e., a 4 × 4 × 2 array of clay unit cells), (II) 51 062 atoms (11 × 11 × 2), (III) 263 750 atoms (25 × 25 × 2), (IV) 1 055 000 atoms (50 × 50 × 2), and (V) 9 495 000 (150 × 150 × 2). The initial model supercells were constructed by replicating the unit cell by the required amount (Table 1). The initial unit cell was energy-minimized before being replicated to create the starting structures of the systems listed in Table 1. Both the clay and the interlayer species were replicated in this way. Therefore, the starting structures for each system in this study vary only in their size. We chose to output the atomic positions every 2.5 ps during production runs as an optimum trade-off between the drop in performance with I/O (input/output) for such large systems and the need to collect data for good sampling and statistics. From this data, radial distribution functions, atomic density profiles perpendicular to the clay surface, diffusion coefficients, and undulatory motion modes were calculated. System wide data, such as the potential energy and lattice constants, were output every 0.01 ps. For stress-strain behavior, we outputted the system-wide pressure every 0.01 ps. Outputting the stress on each atom, used in eq 18 below, can be very expensive due to

looping through the neighbor list and interprocessor communication; the stress-per-atom was therefore output every 5 ps. Molecular dynamics simulations can provide time-dependent quantities directly from the trajectories of the particles. In this paper, we calculate the self-diffusion coefficient of the interlayer water molecules for systems I-IV. Diffusion properties are calculated via the Einstein equation for a fluid confined to two dimensions, which relates the self-diffusion constant D to the mean square displacement of each particle as a function of time

1d 〈|r(t) - r(0)|2〉 4 dt

D)

(1)

where 〈|r(t) - r(0)|2〉 is the mean square displacement (MSD) of the diffusing particle and t is the time elapsed in the simulation. For diffusive motion, the plot of MSD versus t should, therefore, be linear with a slope of 4D. 2.4. Stress-Strain Behavior. Our clay system is anisotropic so that the elastic properties depend on direction. In three dimensions, Hooke’s law reads 3

σij )

3

∑∑Cijklkl k)1 l)1

(2)

where σij is the stress, kl is the strain, and Cijkl is the elastic modulus tensor. There are 81 components of the elastic modulus tensor. We can assume from the unit cell symmetry that our clay system will have orthotropic symmetry. This reduces the number of independent elastic constants required to define the clay sheet to nine, comprising three Young’s moduli Ex, Ey, Ez, three Poisson ratios νyz, νzx, νxy, and three shear moduli Gyz, Gzx, Gxy

xx )

νyx νzx 1 σxx - σyy - σzz Ex Ey Ez

yy ) -

νxy νzy 1 σxx + σyy - σzz Ex Ey Ez

zz ) -

νxz νyz 1 σxz - σyy + σzz Ex Ey Ez yz )

1 σ 2Gyz yz

zx )

1 σ 2Gzx zx

xy )

1 σ 2Gxy xy

(3)

To calculate the elastic properties, we applied uniaxial extension and compression to our simulation cell by changing the length of the periodic simulation cell in the direction of the deformation at every time step. This “slow growth” is valid as long as the


Suter et al.

deformation rate is low. The coordinates of the atoms were rescaled to fit the new geometry. The response stress due to the external tension was computed from the stress tensor, whose components are40,41

σkl )

1

N

∑ V i)1

(

)

mi 1N Firij vivi + 2 2 j)1

∑

(4)

kl

where V is the volume, mi is the mass, vi is the velocity, rij ) (ri - rj), ri is the position, and Fi is the force on the ith particle. The starting structure was comprised of the initial flat clay sheets with the dimensions listed in Table 1. We started the molecular dynamics simulation with the cell length in one axial direction fixed. The system was allowed to relax with this constraint imposed for 0.5 ns, the other directions being at zero pressure. We then deformed the cell length at a constant strain rate using the relation xi ) (1 + i)x, where x is the cell length. The other directions were allowed to relax with the pressure set to 0 atm. This ensured the orthogonal directions did not contribute to the stress-strain behavior. For example, for extension and compression in the x direction, σxx . σyy, σzz and, therefore, xx ) Ex-1σxx. The Young’s modulus for the ith direction Ei was computed from the initial slope of the σii versus i curve. The stress-strain calculation was performed for 2 ns with the strain increasing at each time step. The strain rate is 0.05 ns-1. Although this strain rate is high, it is inherent to the small time scale available in molecular dynamics simulations. However, as long as the strain rate is less than the speed of sound in the material, it should not affect any results calculated. The Poisson’s ratio (νkl) of a material is given by

νkl ) -

k

lateral strain )(tensile strain)

(5)

Figure 5. The height function h(x,y,n) defines the height of the clay sheet above a flat reference plane. The black lines represent the clay sheets; d is the distance between the sheets.

or displacement function, h(x,y,n), where h(x,y,n) ) z(x,y,n) z0,n and where z0,n is the mean z coordinate of the nth clay sheet (z0,n ) ∑i,jz(x,y,n)/NiNj, where Ni and Nj are the number of grid points in the x and y directions, respectively). This is known as the Monge representation (Figure 5). The two-dimensional Fourier transform of h(x,y,n) was computed to yield the q-space mode amplitude, the square of which is the spectral intensity. As we required the thermal averages, the fluctuation spectra calculated for each sheet were averaged over all snapshots. 3. Theory of Clay Sheet Motion 3.1. Undulatory Modes. On sufficiently large length scales, we can forget about the microscopic properties of the clay sheet. We can instead describe the properties of a clay sheet as those of a continuous surface since the ratio of thickness to the inplane dimensions is much less than unity.42,43 Moreover, the amplitude of the fluctuations will be small relative to the lateral dimensions, and there is no spontaneous curvature as the clay is approximately symmetric on either side of the mid-plane. We can, therefore, express the free energy of the bending motion of a single sheet in the form of a curvature expansion, using the Landau-Helfrich bending free energy44

l

where l is the strain in the direction of the uniaxially applied deformation and k is the resulting strain in the transverse direction. In conventional materials, k is negative (contracts) when l is positive (expands) and conversely, k is positive (expands) when l is negative (contracts). This gives a positive νkl. In this study, imposing an extension or compression in the y direction allows us to calculate the in-plane νxy and out-of-plane νzy Poisson ratios from the strains that develop in the relevant directions. Correspondingly, extension and compression in the x direction allow us to calculate νyx and νzx. 2.5. Spectral Analysis of Undulatory Modes. Spectral analysis was performed by mapping the system onto a twodimensional grid for each clay sheet using snapshots from the molecular dynamics simulation generated at 2.5 ps intervals. The atoms in the clay sheet are defined as those comprising the clay framework, that is, not water atoms or sodium ions. We chose a grid spacing of 4.7 Å in both the x and y directions. For the initial (flat) structure, the xy plane contained the clay sheet, while z was perpendicular to the clay surface. We determined the lowest and highest clay species for each species within the grid spacing (x,y); the position of the clay sheet was then calculated as the mean z value of the uppermost and lowermost clay sheet atom. The grid size needs to be large enough to contain a clay sheet atom on both the upper and lower surface of the sheet at each grid point, but it must be sufficiently small to sample high-frequency modes. The center of the nth clay sheet at each grid point z(x,y,n) was replaced by a height

F)

∫ dA[γ + 21k(H - C0)2 + kGK]

(6)

where γ is the external tension, H ) (C1 + C2)/2 is the mean curvature, C0 is the spontaneous curvature, K ) C1C2 is the Gaussian curvature, while k and kG are the bending and saddlesplay moduli, respectively, with C1 and C2 as the local principal curvatures of the surface; C0 ) 0 for a symmetric clay sheet. Fluctuations will not change the topology of the clay sheet so the last term in the integrand in eq 6 is a constant.44 In this study, we have used an NPT ensemble with the lateral and normal directions scaled independently to atmospheric pressure, approximating a state of zero external surface tension. In the Monge representation, the free energy per unit area for an anisotropic system (with different kx and ky)45 can be expressed in terms of the height function

F)

[

]

2 ∂2h(x,y) 1 ∂ h(x,y) kx + k y 2 ∂x2 ∂y2

2

(7)

In Fourier space, this becomes

F)

1 4A

[kxq4x |h(qx)|2 + kyq4y |h(qy)|2] ∑ ∑ q q x

(8)

y

where A is the area in the xy plane. From the equipartition theorem, we know that the energy of each mode q is equal to



(1/2)kBT, and therefore, one can calculate the average mean squared amplitudes for undulating systems

〈|h(qx)|2〉 )

kBT Akxq4x

〈|h(qy)|2〉 )

kBT

(9)

Akyq4y

Here, and throughout, angular brackets 〈‚〉 refer to thermal averages. If a macroscopic surface tension was present, an additional term γq2 would need to be added to the denominator of eq 9, leading to a reduction in the amplitude of the bending undulations. Although the surface with a finite bending stiffness is locally smooth, the orientational correlation decays over long length scales. The spatial decorrelation of the normals to the surface is characterized by the persistence length ξip, such that at length scales larger than ξp in direction i, the surface becomes crumpled

ξip ) ae(2πki/kBT)

(10)

where a is a microscopic distance. The persistence length ξip is thus extremely sensitive to the rigidity constant ki. The persistence length is essentially a correlation length associated with order normal to the surface. Wavelengths close to the persistence length will cause a softening of the bending rigidity. The bending rigidity of a thermotropic liquid crystal is k ∼ 10-20 J, corresponding to 2πk/(kBT) ∼ 12. As ξp depends exponentially on k, the persistence length is 103a ∼ 100 Å.46 A clay sheet is much stiffer than a liquid crystal; therefore, we expect the persistence length to be much larger than any of the dimensions in the systems in our study because of the exponential dependence of eq 10. For multisheet systems, there will also be a potential corresponding to the interaction with other sheets at long wavelengths. This interaction will tend to restore order in the sheet system. This energy term is assumed to depend quadratically on (hn - hn+1), where n is the sheet number of our multisheet system. For N sheets, the free energy per unit area then becomes N-1

F)

[( 1

∑ kx n)1 2

∂2hn ∂x2

+ ky

)

∂2hn ∂y2

2

B + (hn - hn+1)2 2

]

(11)

where k is the bending modulus of a single sheet and B the compressibility modulus. The fluctuations are most conveniently studied in Fourier space in both the z direction and in the x and y directions

〈|h(qx, qz)|2〉 )

nkBT 1 (12) × N‚A 2B[1 - cos(q d)] + k q4 z

x x

nkBT 1 (13) 〈|h(qy, qz)|2〉 ) × N‚A 2B[1 - cos(q d)] + k q4 z

y y

where d is the average interlayer spacing. Integrating eq 13 with respect to qz for i ) x, y gives us

si(qi) )

1 N2

∑ q z

〈|h(qi, qz|2〉 )

1 N-1

∑ 〈hj(qi)‚hn+j(qi)〉

N j)0

(14)

s0(qi) describes correlations within sheets, whereas sn(qi) describes correlations between the sheets. For the two sheet system in our study, s0(qi) is the average of the single sheet

fluctuation spectrum for both sheets. This is the quantity we shall calculate here. In an infinite slab, we can approximate the sum in eq 14 by an integral; combining this with eq 11, we obtain47

s0(qi) )

kB T Akiq4i

[

1+

4 (ξiBqi)4

]

-1/2

(15)

where ξiB ) (ki/B)1/4 is the in-plane correlation length in the i direction. We therefore expect two types of behavior to be manifested in the long-wavelength regime. At smaller wavelengths within this regime, fluctuations of the sheets are incoherent, and the sheets behave like free unconstrained sheets with an intensity proportional to q-4. At longer wavelengths, greater than the in-plane correlation length, the fluctuation spectrum is proportional to q-2. The longer wavelengths correspond to collective motions of the clay sheets whose amplitudes grow more slowly with wavelength than those of free sheets. Below a crossover wavelength, we expect to observe single or collective motion of clay surface atoms, and the continuum picture described above for the long-wavelength regime will no longer be valid. In this short-wavelength regime, the restoring force will be related to the microscopic surface free energy (γp) because the modes tend to increase the local surface area, and we expect their fluctuations to follow a (γpq)-2 behavior. The crossover between the short-wavelength and long-wavelength regimes is denoted qp. We can now divide the undulations into three regions

{

(ξiB)2/2(kcq2i ) qi < (ξiB)-1 kBT (ξiB)-1 < qij < qpi (16) × (kcq4i )-1 s0(qi) ) A (γpq2i )-1 qi > qpi Because clay systems are solid-like, the in-plane diffusion coefficients are zero for the clay sheet atoms, and there is an in-plane shear modulus, which will resist bending. Out-of-plane fluctuations may become coupled to in-plane phonon degrees of freedom at very long wavelengths due to the interaction between shear and compression elastic moduli of the clay, leading to a reduction in out-of-plane amplitude. At wavelengths below this flattening transition (ξf), eq 16 holds. Above the flattening transition, the bending modulus becomes scale dependent, diverging with increasing wavelength. This effect competes with the softening associated with the persistence length ξp, the net result being an ultraviolet stable fixed point associated with the crumpling transition.48 4. Analysis and Results 4.1. Potential Energy Profile. In Figure 6, we show the potential energy profiles as a function of system size for 0-1 ns of simulation time scaled by the standard deviation to allow comparison. The average potential energy and standard deviation are listed in Table 2. We see the potential energy initially decreasing for each system size. The larger systems show the presence of metastable states which relax to more stable states at 550 ps for system III and at 200 ps for system IV. This is due to the development of thermal undulatory fluctuations of the clay sheets, which are discussed in section 4.5. 4.2. Layer Spacing. The layer spacing d is an important characteristic of clay minerals, which is directly amenable from powder X-ray diffraction experiments. The variation of the


Figure 6. Potential energy profiles as a function of time for each system size. The profiles have been scaled by the standard deviation for each system to allow direct comparison. The mean values and standard deviations are listed in Table 2. The systems are (top to bottom) IV, III, II, and I, respectively.

average d spacing once a constant value has been reached is shown in Table 2 for each system size with more than 1 ns of simulation time. We see that the largest system IV has a slightly larger mean d spacing than those of systems II and III. This can be explained by examining the undulatory motion of the clay sheets, which is available as an MPEG movie in the Supporting Information. System III initially forms a standing wave in the y direction. By 0.6 ns, the clay sheet has evolved to an essentially flat state, possessing dynamic low-amplitude undulations in both the x and y directions. After this point, the d spacing is constant at 12.031 Å. System IV adopts a largely flat state, from 0.23 to 0.32 ns, where the d spacing varies between 12.038 and 12.030 Å (shown in Figure 7). After 0.35 ns and for the rest of the simulation, system IV supports a standing wave in the x direction and has a higher d spacing than the flat state (12.049 Å). However, we can see in the potential energy profile (Figure 6) that the flat state (0.230.32 ns) has a lower potential energy than other sheet conformations. We expect that, eventually, system IV will relax back to this state over a long enough time scale. Subsequent analysis of clay sheet conformations will use data from this flat state. 4.3. Radial Distribution Functions, Density Profiles, and Water Self-Diffusion. In order to characterize the variation of the structure of the clay interlayer with system size, we determined radial distribution functions (RDF) and density profiles perpendicular to the clay surface. The first peak of the sodium ion radial distribution function with water oxygen atoms, corresponding to the first coordination sphere, is found at 2.4 Å for all system sizes. We find the first peak in the RDF between the sodium ions and the clay surface oxygen atoms is at approximately the same distance as the sodium-water oxygen radial distribution function, indicating the coordination sphere of sodium is composed of both oxygen atoms from the interlayer water and the clay surface. There is no discernible difference with system size. The density profiles of the interlayer species are shown in Figure 8. With increasing system size, it is possible to observe a slight broadening of the peaks of each species, as evidenced by the overlap between the water hydrogen and tetrahedral layer oxygen atom distributions and the broadening of the surface oxygen peak for system IV. This broadening of the density profiles was observed by Greenwell et al.;49 it is due to the

Suter et al. thermal undulations in clay sheets present only at large system sizes, which cause the z position of each species to becomes less well localized. In Table 2, we report the diffusion coefficients of interlayer water molecules for systems that have been simulated for greater than 1 ns, I-IV. The MSD are calculated for 1 ns after equilibrium has been reached. We find, as expected, that diffusion is only in the x and y directions. Diffusion in the direction perpendicular to the clay surface is zero. We observe that the self-diffusion coefficients increase only slightly with increasing system size. The maximum MSD for the water molecules does not exceed the lateral dimensions of any system at 1 ns. It should be noted that the value of ca. 2.3 × 10-10 m2 s-1 calculated for the largest system model is similar to the self-diffusion coefficient reported from experiments using quasi-elastic neutron scattering.39,50 We observe a small increase in the diffusion coefficient with increasing system size; we can tentatively assign this to the thermal undulations of the clay sheet present only in the larger models enhancing the mobility of the water molecules. 4.4. Stress-Strain Behavior. In this section, we derive stress-strain curves from molecular simulations using the method described in section 2.4 to predict the elastic moduli of our clay systems.51,52 In Figure 9, we show the stress-strain curve for systems I-IV in the x and y directions. We have plotted a running average of the stress, with each value averaged over 10 ps. All system sizes show almost identical stress-strain behavior for positive strain values. Buckling occurs at smaller compressive strain values as the system size increases. We find an average Young’s modulus of 172 GPa in the x direction (Ex) and 182 GPa in the y direction (Ey). Table 3 shows the elastic constants generally decreasing as the size of the simulation increases. This is presumably due to flattening of the thermal undulations, which become more prevalent as the system size increases. The force required to flatten the undulations is lower than that required to stretch the clay sheet. Similarly, we can compute Poisson ratios from eq 5 and find that νzx ≈ νyz ) 0.17 and νxy ≈ νyx ) 0.36. It is also possible to estimate elastic constants from the thermal fluctuations in the shape and size of our simulation cell. We assume that the state of zero strain is defined by the average shape of the simulation cell and the instantaneous strain is given by the change in the cell length, that is, ij ) ∆l/l. The components of the elastic modulus tensor in eq 2 can be extracted using the Parrinello-Rahman fluctuation formula53

Cijkl )

kBT

1 〈V〉 〈ijkl〉

(17)

We have used eq 17 to estimate Ex and Ey in Table 3 from the strain fluctuations in the x and y directions, respectively. Although this method generally requires long simulation times to correctly sample all configurations of the simulation cell,54 we see good agreement with the Young’s moduli calculated using the deformation approach described above. The Young’s moduli calculated above are for the clay sheets and interlayer spacing. However, composite theory of clay polymer nanocomposites containing exfoliated clay sheets requires knowledge of the elastic constants of a clay platelet.20 We therefore calculate the Young’s modulus for an isolated clay sheet by modifying eq 4 such that the stress response to the



TABLE 2: Potential Energy, Dimensions, and Self-Diffusion Coefficients for Interlayer Water Averaged over the Simulation after Equilibrium Is Reached for Systems with Greater than 1 ns of Simulation Timea system

potential energy eV

average x (Lx) Å

average y (Ly) Å

average z (Ly) Å

self-diffusion coeff. Dw 10-10 m2 s-1

I II III IV

-1.26 × 106 ( 6.5 × 101 -9.50 × 106 ( 3.8 × 102 -4.91 × 107 ( 1.5 × 103 -1.96 × 108 ( 1.1 × 103

41.50 ( 0.024 114.12 ( 0.029 259.30 ( 0.025 518.61 ( 0.024

72.00 ( 0.040 198.02 ( 0.043 450.13 ( 0.044 900.31 ( 0.037

24.102 ( 0.05 24.078 ( 0.013 24.062 ( 0.006 24.097 ( 0.003

1.962 ( 0.0024 2.157 ( 0.0014 2.258 ( 0.0006 2.334 ( 0.0003

a

Equilibrium was designated to have occurred when no drift in the potential energy was observed, or at 0.5 ns, whichever was the greater.

Figure 7. Average d spacing of each clay model system. Each system contains two clay sheets; therefore, the d spacing is half the c dimension of the simulation cell. The solid line refers to the average d spacing for all systems after a constant value has been reached or after 0.5 ns, whichever is the greater, except for system IV, where the d spacing for the flat state (0.23-0.32 ns) is reported.

Figure 9. The stress in the x (upper panel) and y (lower panel) direction as a function of strain for systems I-IV during deformation in the x and y axial directions, respectively. Each point is a running average over 10 ps time intervals. The gradient of the stress-strain curve is the Young’s modulus, Ex, for strain in the x direction, and Ey, for strain in the y direction.

Figure 8. Atom density profiles perpendicular to the clay sheets of the interlayer species: Na (dashed), water oxygen (dot-dashed), and water hydrogen (dotted). The clay oxygen atoms are also indicated (solid). The systems are I, II, III, and IV from top to bottom. For system IV, the density distribution is noticeably broader.

stepwise strain is from the clay sheet alone, not clay plus interlayer species. The stress tensor elements are therefore

σkl )

1 Vclay

(

Nclay

∑ i)1

)

mi 1N Firij vivi + 2 2 j)1

∑

(18)

kl

where Nclay is the number of atoms in the clay sheets, Vclay ≡ n‚hclayA, hclay is the thickness of the clay sheet, n is the number of clay sheets in the system, and i labels the atoms in the clay sheets. The value of the Young’s modulus for a single clay sheet is therefore inversely proportional to the clay thickness hclay in eq 18. As a first estimate, we use the average distance between

the uppermost and lowermost clay layer atoms for each grid point fitted to the clay sheet, as described in section 2.5. Averaged over all system sizes, we find hclay to be 6.69 Å. However, this may underestimate the thickness as it does not include the size of the surface oxygen atoms. If we include the van der Waals radius of the oxygen surface atoms (1.34 Å), the clay sheet has a thickness of 9.37 Å. In Table 3, we estimate effective Young’s moduli of a single clay sheet using both values of hclay. We find an average Ex of 320 and 233 GPa for the lower and higher hclay values, respectively, with 344 and 247 GPa as the corresponding Ey moduli. 4.5. Undulatory Motions and Elastic Constants. We now proceed to undertake a quantitative analysis of the thermal undulations of the clay sheets observed through visualization. This requires us to define the position of the clay surface as a function of the xy plane (see Figure 5) by mapping the coordinates onto a 2D grid (x,y). Figure 10 shows the time development of the height fluctuations of a single sheet within system IV, containing 1 055 000


Suter et al.

TABLE 3: In-Plane Elastic Constants in GPa Calculated from (i) Stress-Strain Behavior and Average Strain Fluctuations for the Whole System Containing the Clay Framework Plus Interlayer, Output Every 0.01 ps, and (ii) Using the Stress-Strain Behavior of the Clay Sheet Elements (eq 18), Output Every 5 psa (i) clay and interlayer

(ii) clay sheet

system

Exb

Ey

Exc

Eyd

Exe hclay ) 9.37 Å

I II III IV

174 ( 0.6 173 ( 0.2 171 ( 0.1 182 ( 0.1

187 ( 0.8 190 ( 0.3 189 ( 0.1 192 ( 0.1

172 175 171 173

183 182 181 180

210 ( 32 236 ( 13 235 ( 5 231 ( 2

Ey hclay ) 9.37 Å

Ex hclay ) 6.69 Å

Ey hclay ) 6.69 Å

254 ( 38 246 ( 13 245 ( 3 245 ( 2

294 ( 48 331 ( 6 330 ( 7 322 ( 3

356 ( 51 344 ( 18 343 ( 4 343 ( 2

a This latter value is dependent on the thickness of the clay sheet, h clay. We have chosen two values, 6.69 and 9.37 Å, corresponding to the distance between surface oxygen atoms and this value plus twice the van der Waals radius of an oxygen atom. b Using eq 4. c Using eq 17 with 〈xxxx〉. d Using eq 17 with 〈yyyy〉. e Using eq 18.

Figure 10. Evolution of long-wavelength undulatory modes illustrated by the height function h(x,y,1) for a single sheet of the 1 055 000 atom system IV. We see after 0.3 ns that the clay sheet is still in a predominantly flat phase, but the longest wavelength undulations possible in the x direction (Lx) are visible after 500 ps. In the y direction, only smaller wavelength undulations can be seen, even after 2 ns.

Figure 11. Height functions h(x,y,1) for a single clay aluminosilicate sheet, as determined for clay system sizes after 1 ns of simulation for (a) I 6752 atoms, (b) II 51 062 atoms, (c) III 263 750 atoms, and after 0.5 ns for (d) V 9 495 000 atoms.

atoms. Figure 11 illustrates height functions after 1 ns of simulation for systems I-III and after 0.5 ns for system V. In Figure 10, we see the development of long-range undulatory modes within several hundred picoseconds. At the end of the simulation, an undulatory mode in the x direction with wavelength Lx can be observed, although the longest wavelength Ly cannot be clearly seen in the y direction. Animations of the height functions are available in the Supporting Information.

Figure 12. Spectral intensity per undulatory mode versus the wave vector in the x and y directions for systems I-IV (top panel) and in the y direction for system V (lower panel). The dashed line is a fit to the undulatory q-4 modes for q < 0.04 with k ) 1.6 × 1017 J. This indicates that the clay sheet is acting like a single elastic sheet. The spectral intensity for system IV is from the predominantly flat state (0.23-0.32 ns). The dash-dotted line shows q-2 behavior at higher q values. This is due to surface atom interactions with the interlayer species. The local peak centered at q ) 0.35 Å-1 is due to artificial periodicity resulting from the way isomorphic substitutions are included in all of the models. Isomorphic substitutions in real clays are randomly distributed; therefore, this peak can be disregarded. The crossover between q-4 and q-2 behavior is at approximately q ) 0.04. This is of such a long wavelength that only systems III-V possess a significant number of longer wavelength modes.

Figure 12 shows the spectral intensity in both the x and y directions for systems I-IV, that is, s0(qx) ) 〈|h(qx)|2A〉 and s0(qy) ) 〈|h(qy)|2A〉, where A is the averaged projected area (i.e, 〈Lx × Ly〉), versus q ) 2π/λ, where λ is the undulation wavelength (eq 14). Only the y direction is shown for system V in Figure 12 as only undulations in this direction were observed during the short simulation timespan.

Molecular Dynamics Study of Montmorillonite Clay Several features can be seen in the fluctuation spectrum. A local maximum is observed at q ) 0.35 Å-1, which corresponds to a wavelength of approximately 18 Å. This distance is the same as the unit cell dimensions in the y direction (Figure 2). Recall that the unit cell was replicated to produce the starting structure of the larger systems. As a result, this also corresponds to the distance between sites of isomorphic substitution in the tetrahedral layer and the associated sodium ions adsorbed on the clay surface in the starting structure. It is this feature that results in the observed periodicity. Cation adsorption leads to local raising of the position of the clay sheet atoms (Figure 3); the oxygen atoms on the clay surface adjacent to the adsorbing sodium ions protrude further into the interlayer region than other oxygen atoms in the clay surface. As the change in h(x,y,n) is purely a function of isomorphic substitution and sodium ion adsorption, the height of the fluctuations is expected to be independent of system size. Thus, when scaled by the projected area h‚A to give the spectral intensity, we see the peak increasing with system size. By contrast, the spectral intensity for all other modes is independent of system size. Increasing the system size introduces new modes with longer wavelengths without changing the shorter wavelength modes. We can observe two distinct regimes. At q values smaller than 0.04 Å-1, we have a spectral distribution proportional to q-4, corresponding to the mesoscopic undulations predicted by eq 16. For q values greater than 0.04 Å-1, we observe q-2 behavior, indicating surface clay atom protrusions into the interlayer spacing. The crossover between the two regimes, qp, is of such a long wavelength (∼150 Å) that it can only be seen in the largest size simulations, which fit our definition of large-scale simulations comprising more than 100 000 atoms, systems III-V. Fitting q values lower than qp in Figure 12 to q-4 behavior (dashed line) for systems I-IV, we estimate the bending modulus to be approximately 1.6 ( 0.02 × 10-17 J. The fitting of the q-2 behavior for q > 0.04 Å-1 gives an approximate surface free energy of 8 ( 0.1 N m-1 (dot-dashed line). This value is extracted from shorter wavelength undulations and is, hence, unaffected by system size. In passing, we note that the bending modulus is approximately 3 orders of magnitude higher than that found for the bending of lipid bilayers13-16 and nonbiological membranes.17 This is due to the much greater resistance to bending of the clay sheets, illustrated by the smaller bending amplitudes. The attraction between aluminosilicate atoms is much greater than that between lipids, and the energies of compression and stretching are far higher. For the largest system in our study, comprising approximately 10 million atoms, we see that the longest wavelength contributes much less than expected from the q-4 behavior. This is caused by the existence of a flattening transition point (ξf). At wavelengths greater than ξf, the bending modulus is enhanced, causing the surface to become flattened over long distances. This is due to the coupling of in-plane phonon modes with outof-plane displacements.55 We estimate ξyf to be approximately 1400 Å, after which the bending rigidity would become scaledependent. We can now attempt to translate our calculated bending modulus into the elastic constants of a macroscopic body and compare the results to bulk measurements. To estimate the inplane Young’s moduli for a clay sheet, we use the theory of bending of thin films43


ki )

Eih3clay 12(1 - ν2)

(19)

where Ei is Young’s modulus in the i direction, ν is the Poisson’s ratio, hclay is the thickness of the sheets, and ki is the bending modulus in the i direction. To calculate Young’s moduli using eq 19, it is therefore necessary to define the thickness of the clay sheet. In section 4.4, we calculated the in-plane Young’s moduli for a clay sheet using the stress-strain behavior of the clay sheet elements (eq 18). The Young’s moduli calculated using this method are inversely proportional to our choice of clay sheet thickness. Equation 19 gives us a new relationship between the Young’s modulus and the clay sheet thickness. We can, therefore, define an effective clay thickness such that eqs 19 and 18 give the same Young’s moduli. Using ν ) 0.36 in eq 19 for the Poisson’s ratio calculated from the stress-strain behavior in section 4.4 and comparing with the results of eq 18, we find an identical Young’s modulus with an effective thickness of approximately 9.0 Å. This is remarkably close to the thickness we calculated by including the van der Waals radius of surface oxygen atoms (9.37 Å). Experimentally, the thickness reported for collapsed montmorillonite clay is 9.8 Å,19 although this will include a small interlayer gallery where the sodium ions reside. Using 9.0 Å for hclay, we conclude that for clay sheets, Exx = Eyy ≈ 230 ( 4 GPa. It is hard to determine Young’s moduli from experiment for clay minerals due to the small size of clay platelets, but the values we calculated from the bending modulus are in line with previous simulations by Mantevitch and Rutledge of a single 288 atom sheet of montmorillonite under slow rates of deformation. These authors predicted a bulk modulus of 270 GPa for a bulk clay including intersheet gallery and 400 GPa for a single sheet.51 In their study, the calculated value was Eihclay, leading to uncertainty in the Young’s modulus due to the uncertain value of hclay, which they estimated to be 6.78 Å. They also estimated the bending modulus from the critical stress for nonlinear compression to be 1.25 × 10-17 J. Differences may be due to force-field effects and/or the excessively small system size used in their study. The system used in their study was an uncharged clay sheet in vacuo; the presence of the sodium ions and interlayer water molecules in this study could lead to a stiffer and thicker clay sheet. As experimental data have not been forthcoming, Chen and Evans surveyed estimates of smectite clay platelets19 by focusing on analogous materials and moduli-density relations. These authors concluded that most estimates of Young’s moduli lie in the range of 178-265 GPa. The calculated Young’s moduli from our simulations using the bending modulus are in agreement with this figure. 5. Conclusions We have demonstrated in this paper that the dynamics of clay sheets takes place on a wide variety of length scales, which can be observed only through large-scale atomistic simulation. We find that “local” properties, such as radial distribution functions and diffusion coefficients, are essentially unchanged with system size. Analyzing the dynamics of the clay sheet at length scales less than 150 Å, interactions of surface atoms with ions, atoms, and molecules in the interlayer lead to single and collective atomic protrusions into the interlayer spacing. This results in a

8258 J. Phys. Chem. C, Vol. 111, No. 23, 2007 “roughening” of the clay surface. These protrusions are described by an effective microscopic surface free energy of 8 N m-1. At larger length scales, we find emergent behavior of collective undulatory motions of the clay sheets, which is well described by approximating each clay sheet as a single, continuous surface. We parametrize this motion by a bending modulus, which we find to be 1.6 × 10-17 J. At length scales of approximately 0.14 µm, we see a sharp reduction in the amplitude of the longest wavelength undulation, which we assign to a wavelength dependence of the bending modulus beyond a flattening transition point. This is due to the coupling of the in-plane phonon modes and the out-of-plane undulations, which causes the surface to be stiffened over long distances. We have calculated the stress-strain behavior by increasing the strain in the system and measuring the resultant internal stress. From this, we find the in-plane Young’s modulus (Ex, Ey) of the clay plus interlayer to be 170-190 GPa. We have calculated the stress response from the atoms within the clay layer, allowing us to estimate Young’s moduli for a single clay sheet, if the sheet thickness is specified. We find Young’s moduli in the range of 230-360 GPa depending on our definition of clay thickness. It is also possible to calculate Young’s moduli from the bending modulus of the clay sheet calculated from the undulatory motion, using the Poisson ratios calculated in the stress-strain behavior simulations. This allows us to estimate an effective clay thickness, for which both the stress-strain behavior and bending modulus results produce the same in-plane Young’s moduli for the clay sheet. We find an effective clay sheet thickness hclay ≈ 9.0 Å and the resultant Young’s moduli Ex, Ey to be on the order of 230 GPa. In future simulation work, we plan to study the effect of clay edges explicitly. We plan to simulate montmorillonite clay without periodic boundaries, with the clay sheets terminated and surrounded by water molecules. Edge fluctuations for solid sheets are expected to be higher than those far from the boundaries.56 It may be possible to observe a crumpling transition,57,58 which has been previously reported for other flat surfaces such as tethered membranes.59 We also hope to perform very large-scale simulations of montmorillonite clay platelets embedded in a polymer matrix, a situation that arises in claypolymer nanocomposites. The high-performance grid-computing infrastructure described here should render such investigations feasible on a large enough scale to capture similar mesoscopic effects to those described in this paper. Acknowledgment. This work was partially supported by EPSRC (GR/T27488/01), which provided access to HPCx (www.hpcx.ac.uk). Our work was also supported by the National Science Foundation under NRAC Grant MCA04N014, utilizing resources on the US TeraGrid (www.teragrid.org), and the DEISA Consortium (co-funded by the EU, FP6 Project 508830) within the DEISA Extreme Computing Initiative (www. deisa.org). M.-A.T. was funded by an EPSRC Ph.D. studentship associated with RealityGrid (www.realitygrid.org). We are grateful to Fiona Reid at EPCC, University of Edinburgh, for benchmarking data for LAMMPS on HPCx and Randy Cygan of Sandia National Laboratory for providing the montmorillonite clay model used in this study. Supporting Information Available: Radial distribution functions of interlayer species, MPEG animations of height functions of systems I-IV and MD trajectories for 0-2 ns of systems II-IV with a z expansion factor of 10 and 15 for

Suter et al. systems III and IV, respectively. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Sposito, G. The Surface Chemistry of Solids; Oxford University Press: New York, 1984. (2) Brindley, G. W.; Brown, G. Crystal Structures of Clay Minerals and Their X-Ray Identifaction; Mineralogical Society: London, 1980. (3) Giannelis, E. P. AdV. Mater. 1996, 8, 29. (4) Alexandre, M.; Dubois, P. Mater. Sci. Eng. 2000, 28, 1. (5) LeBaron, P. C.; Wang, Z.; Pinnavaia, T. J. Appl. Clay Sci. 1999, 15, 11. (6) Greenwell, H. C.; Bowden, A. A.; Chen, B.; Boulet, P.; Evans, J. R. G.; Coveney, P. V.; Whiting, A. J. Mater. Chem. 2006, 16, 1082. (7) Greenwell, H. C.; Jones, W.; Coveney, P. V.; Stackhouse, S. J. Mater. Chem. 2006, 16, 708. (8) Skipper, N. T.; Chang, F. R. C.; Sposito, G. Clays Clay Miner. 1995, 43, 294. (9) Skipper, N. T.; Chang, F. R. C; Sposito, G. Clays Clay Miner. 1995, 43, 285. (10) Stackhouse, S.; Coveney, P. V.; Sandre´, E. J. Am. Chem. Soc. 2001, 123, 11764. (11) Boulet, P.; Coveney, P. V.; Stackhouse, S. Chem. Phys. Lett. 2004, 389, 261. (12) Greenwell, H. C.; Harvey, M. J.; Boulet, P.; Bowden, A. A.; Coveney, P. V.; Whiting, A. Macromolecules 2005, 38, 6189. (13) Marrink, S. J.; Mark, A. E. J. Phys. Chem. B 2001, 105, 6122. (14) Lindahl, E.; Edholm, O. Biophys. J. 2000, 79, 426. (15) Goetz, R.; Gompper, R.; Lipowsky, R. Phys. ReV. Lett. 1999, 82, 221. (16) den Otter, W. K.; Briels, W. K. J. Chem. Phys. 2003, 118, 4712. (17) Boek, E. S; Padding, J.; den Otter, W. K.; Briels, W. K. J. Phys. Chem B 2005, 109, 19851. (18) Lipowsky, R. Nature 1991, 349, 475. (19) Chen, B.; Evans, J. R. G. Scr. Mater. 2006, 54, 1581. (20) Sheng, N.; Boyce, M. C.; Parks, D. M.; Rutledge, G. C.; Abes, J. I.; Cohen, R. E. Polymer 2004, 45, 487. (21) Zhu, L.; Narh, K. A. J. Polym. Sci., Part B: Polym. Phys. 2004, 42, 2391. (22) Buryachenko, V.; Roy, A.; Lafdi, K.; Anderson, K.; Chellapilla, S. Compos. Sci. Technol. 2005, 65, 2435. (23) Boek, E. S.; Coveney, P. V.; Skipper, N. T. J. Am. Chem. Soc. 1995, 117, 12608. (24) Plimpton, S. J. J. Comput. Phys. 1995, 117, 1. (25) Plimpton, S. J. Large-Scale Atomic/Molecular Massively Parallel Simulator. http://lammps.sandia.gov (accessed 12/4/06), Sandia National Laboratories: Albuquerque, NM, 2005. (26) Hein, J. I.; Reid, F.; Smith, L.; Bush, I.; Guest, M.; Sherwood, P. Philos. Trans. R. Soc. London, Ser. A 2005, 363, 1833. (27) Coveney, P. V.; Saksena, R.; Zasada, S. J.; McKeown, M.; Pickles, S. Comput. Phys. Commun. 2007, 176, 406. (28) Coveney, P. V.; Vicary, J.; Chin, J.; Harvey, M. Philos. Trans. R. Soc. London, Ser. A. 2005, 363, 1807. (29) Li, J. Modell. Simul. Mater. Sci. Eng. 2003, 11, 173. (30) Cygan, R. T.; Liang, J. J.; Kalinichev, A. G. J. Phys. Chem. B 2004, 108, 1255-1266. (31) Greathouse, J. A.; Cygan, R. T. Phys. Chem. Chem. Phys. 2005, 7, 3580. (32) Kalinichev, A. G.; Kirkpatrick, R. J.; Cygan, R. T. Am. Mineral. 2000, 85, 1046. (33) Wang, J.; Kalinichev, A. G.; Kirkpatrick, R. J.; Hou, X. Chem. Mater. 2001, 13, 145. (34) Kalinichev, A. G.; Kirkpatrick, R. J. Chem. Mater. 2002, 14, 3539. (35) Kirkpatrick, R. J.; Yu, P.; Hou, X.; Kim, Y. Am. Mineral. 1999, 84, 1186. (36) Wang, J.; Kalinichev, A. G.; Amonette, J. E.; Kirkpatrick, R. J. Am. Mineral. 2003, 88, 398. (37) Wang, J.; Kalinichev, A. G.; Kirkpatrick, R. J. Geochim. Cosmochim. Acta 2004, 68, 3351. (38) Berendsen, H. J. C.; Postma, J. P. M.; von Gunsteren, W. F.; Hermans, J. Interaction Models for Water in Relation to Protein Hydration. In Intermolecular Forces; Pullman, B., Ed.; D. Reidel: Dordrecht, Holland, 1981. (39) Skipper, N. T.; Refson, K.; McConnell, J. D. C. J. Chem. Phys. 1991, 94, 7434-7445. (40) Qi, D.; Hinkley, J.; He, G. Modell. Simul. Mater. Sci. Eng. 2005, 13, 493. (41) Frankland, S. J. V.; Harik, V. M.; Odegard, G. M.; Brenner, D. W.; Gates, T. S. Compos. Sci. Technol. 2003, 63, 1655.

Molecular Dynamics Study of Montmorillonite Clay (42) McWhirter, J. L.; Ayton, G.; Voth, G. A. Biophys. J. 2004, 87, 3242. (43) Landau, L. D.; Lifshitz, E. M. Theory of Elasticity; Pergamon: New York, 1970. (44) Helfrich, W. Z. Naturforsch., C: Biosci. 1973, 28, 693. (45) Xing, X.; Mukhopadhyay, R.; Lubensky, T. C.; Radzihovsky, L. Phys. ReV. E 2003, 68, 021108. (46) de Gennes, P. G.; Taupin, C. J. Phys. Chem. 1982, 86, 2294. (47) Loison, C.; Mareschal, M.; Kremer, K.; Schmid, F. J. Chem. Phys. 2003, 24, 13138. (48) Bowick, M.; Catterall, S. M.; Falcioni, M.; Thorleifsson, G.; Anagnostopoulos, K. J. Phys. I 1996, 10, 1321. (49) Greenwell, H. C.; Harvey, M. J.; Boulet, P.; Bowden, A.; Coveney, P. V.; Whiting, A. Macromolecules 2005, 38, 6189.

J. Phys. Chem. C, Vol. 111, No. 23, 2007 8259 (50) Skipper, N.; Lock, P.; Titiloye, J.; Swenson, J.; Mirza, Z.; Howells, W.; Fernandez-Alonso, F. Chem. Geol. 2006, 230, 182. (51) Manevitch, O. L.; Rutledge, G. C. J. Phys. Chem. B 2004, 108, 1428. (52) Grigoras, S.; Gusev, A. A.; Santos, S.; Suter, U. W. Polymer 2002, 43, 489. (53) Parrinello, M.; Rahman, A. J. Chem. Phys. 1982, 76, 2662. (54) Sprik, M.; Impey, R. W.; Klein, M. L. Phys. ReV. B 1984, 29, 4368. (55) Aronovitz, J. A.; Lubensky, T. C. Phys. ReV. Lett. 1988, 60, 2364. (56) Abraham, F. Phys. ReV. Lett. 1991, 67, 1669. (57) Kantor, Y.; Nelson, D. R. Phys. ReV. Lett. 1987, 26, 2774. (58) Vliegenthart, G.; Gompper, G. Nat. Mater. 2006, 5, 216. (59) Zhang, Z.; Davis, H. T. Phys. ReV. E 1996, 53, 1422.

Large-Scale Molecular Dynamics Study of Montmorillonite Clay

Recommend Documents