Abstract
This paper presents results of a molecular dynamics simulation study of dehydrated 2:1 clay minerals using the Parrinello-Rahman constant-pressure molecular dynamics method. The method is capable of simulating a system under the most general applied stress conditions by considering the changes of MD cell size and shape. Given the advantage of the method, it is the major goal of the paper to investigate the influence of imposed cell boundary conditions on the molecular structural transformation of 2:1 clay minerals under different normal pressures. Simulation results show that the degrees of freedom of the simulation cell (i.e., whether the cell size or shape change is allowed) determines the final equilibrated crystal structure of clay minerals. Both the MD method and the static method have successfully revealed unforeseen structural transformations of clay minerals upon relaxation under different normal pressures. It is found that large shear distortions of clay minerals occur when full allowance is given to the cell size and shape change. A complete elimination of the interlayer spacing is observed in a static simulation. However, when only the cell size change is allowed, interlayer spacing is retained, but large internal shear stresses also exist.
1. Introduction
In the last two decades, computer simulation methods, including the Monte Carlo (MC) method and the Molecular Dynamics (MD) method, have been successfully used in the studies of dehydrates/hydrates of 2:1 clay minerals [1–4]. Computational molecular models based on well-calibrated intermolecular potentials, supported by ever-increasing computer capabilities, have offered unmatched advantages in probing the molecular mechanisms and surface properties of the clay-water system [2, 3, 5–13]. Of the two methods, however, the MC method has received a much wider application in previous studies, mainly because it is relatively easy for the method to calculate desired properties of the system under any specified pressure, temperature, and volume. Based on the repeated random sampling and probability theory, the MC method is a procedure for evaluating average properties of phase space equilibrium of constant temperature ensembles, such as the canonical (NVT) ensemble and the isothermal-isobaric (NPT) ensemble. Average macrostate properties of a clay-mineral system, particularly under the conditions of volume and/or pressure changes, can thus be readily predicted by the method. A major limitation of the MC method, however, is that it only deals with position variables of the system and provides no information about the time dependency of properties and their fluctuations.
In comparison, the MD method, which is based on the Newtonian equations of motion of particles, is a classical method to calculate average properties of a molecular system varying as a function of time. For the conventional MD method, the energy, volume, and number of particles are fixed for a particular system, and it is assumed that time averages of properties of the system are equal to the microcanonical (NVE) ensemble averages of the same properties. When such a system moves along its trajectory, the pressure and temperature will change. Such a limitation hinders the application of the MD method to the study of clay minerals to a large extent.
To remove such a limitation and render the MD method applicable to more general conditions, Andersen [15] proposed new MD methods which can handle the simulation at constant pressure and/or temperature. Following Andersen [15], Parrinello and Rahman [16] extended the method by taking into account the MD cell shape change. This is achieved by introducing a tensor defining fully the cell size and shape, which are incorporated as extra degrees of freedom into Hamiltonian equations of motion.
The P-R method has been frequently used in the studies of structural transformations of solids or liquids with volume and temperature fluctuations under constant external stress conditions. For example, in the same paper, Parrinello and Rahman [16] applied the method to the simulation of a lattice model of Nickel subjected to uniaxial compression and tensile load and uncovered unforeseen patterns of crystal transition. Yamakov et al. [17, 18] studied the nanocrystalline aluminum using the P-R method and showed that it was a powerful new tool for elucidating and quantifying the atomic level mechanisms controlling the complex dislocation and grain-boundary processes in heavily deformed nanocrystalline materials. Success was also achieved by Tsuneyuki et al. [19] in applying the method to the simulation of pressure-induced structural transitions of various polymorphs of silica, which experiences discontinuous volume reduction and novel structural phases with increasing hydrostatic pressures. Little study, however, has been published on the application of the method to clay minerals. In the development and employment of a general force field in the simulation of clay minerals, Cygan et al. [20] adopted the P-R technique, but no attempt was made in exploring the effects of external stress and cell boundary condition on the structural transition behavior of clay minerals.
In this paper, we apply for the first time the P-R method to the simulations of 2:1 clay minerals. The focus of the study is on the interlayer molecular structural transformations of dehydrated dioctahedral smectites bearing octahedral layer charges. Given the advantage of the method, it is our major goal to predict the behavior of dehydrated smectite sheets under a constant normal pressure but varying degrees of cell boundary constraints. In particular, the effects of cell volume and shape changes on the structural transformation of clay minerals will be fully investigated. The coupled changes of cell configuration and atomic positions driven by the external and internal stress imbalance, which cannot be predicted by the conventional MD method, provide new insights on the atomic-scale behavior of clay minerals.
To better study the equilibrium structures of clay minerals after a large number of time steps, over which cumulative computational errors may affect the simulation results, we also study the problem by using a static version of the method, in which time variations of system properties are disregarded, and only atom positional variables are involved. As a result, the static simulation provides results that represent the behavior of clay minerals at zero Kelvin.
2. MD Simulation Method
2.1. Molecular Model of Na-Saturated Smectite
The unit-cell formula of Na-saturated dioctahedral smectite is [21]: where is the overall layer charge neutralized by the Na+ counterions, with and equal to the number of tetrahedral substitutions of Si by Al and octahedral substitutions of Al by Mg, respectively. In this paper, we studied the behavior of true micas bearing an octahedral layer charge of −2.0 (i.e., , ) [22]. The crystallographic structure of mica is monoclinic (C2/c) with unit cell dimensions Å, Å, and Å. The initial atom positions and charges are listed in Table 1 (adopted from Skipper et al. [2]).
To appropriately represent the infinitely extended clay sheets using very limited number of repeating unit cells in the simulation cell, 3D periodic boundary conditions must be imposed to ensure that every atom interacts with all the other atoms in the simulation cell as well as their infinitely repeating images (including images of itself) outside of the simulation cell [2, 23, 24]. Practically, this representation scheme can be greatly simplified by applying a cut-off to the short-range contributions to the interatomic potentials so that every atom only interacts with its nearest images (called the “minimum image convention”). As a result, the minimum length of any side of the simulation cell is twice the value of the cut-off. In this study, we have selected 9 Å for the cut-off, resulting in a 21.12 Å × 18.3 Å rectangle in the basal plane including eight repeating unit cells. In the dimension normal to the basal plane (i.e., [] direction), we selected to include two single sheets of mica because we found that this was the minimum size that the simulation cell could have to be faithfully representative of the macroscopic clay system. The dimension of the simulation cell in the [] direction, therefore, is twice the sum of clay sheet thickness (i.e., c = 6.56 Å) and interlayer spacing (d). This selection was based on an extensive parametric study, which showed that other smaller simulation cells with less number of atoms had numerical problems and could not converge into a stable condition in most cases. Under this arrangement, a typical simulation cell for the dehydrated mica contains a total number of 672 atoms. Figure 1 shows the crystal structure of the dehydrated mica contained in a typical simulation cell.
(a)
(b)
(c)
2.2. Interatomic Potential Functions
In the current study, we have adopted a uniform potential functional form for all atom-atom pairs (i.e., two-body term) in the simulation system as follows [25, 26]:
Equation (2) consists of Coulomb energy (the first term), short-range repulsion with Born-Mayer-Higgins (the second term), van der Waals attraction (the third term), and covalent Morse terms (the last two terms). In (2), q is the charge for each atom, is the interatomic distant, C2N-1m-2 is the dielectric constant for vacuum, N is the constant for unit conversion, a, b, and c are material constants of different atom species, and and () denote the depths and shape of the Modified Morse type potential, respectively. This comprehensive potential function has been successfully used in the MC simulations [2, 9] and MD simulations [25, 26] of smectite-type clay minerals. To determine the parameters , we have followed the scheme by Skipper et al. [2], which is based on the optimal use of MCY water-water model [27] by representing the clay mineral solely by water molecules (under the consideration that the valence electrons in hydroxyl, tetrahedral, and octahedral groups are centered on the O atoms, which will dominate any short-range interaction with a clay mineral layer). The resulting parameter values are listed in Table 2. For the octahedral substitution of Al by Mg, the charge on the cation site is simply reduced from +3e to +2e whereas the tetrahedral substitution of Si by Al corresponds to the charge reduction from +1.2e to +0.2e. The parameters in (2), as listed in Table 3, are determined so as to reproduce the structural and physical properties of water and several oxide crystals [25, 28, 29].
In the handling of long-range coulomb interaction in the first term in (2), a real-space cut-off cannot be applied to the summation of this term over an infinite number of atom pairs without causing serious nonconvergence problems [24]. Therefore, the Ewald summation technique [23, 30, 31] was used to treat this problem. In the Ewald sum, the total Coulomb energy is divided into two rapidly convergent series, one in real space and the other in reciprocal space. The rate of convergence depends on a common parameter in both series, and a balance must be achieved between accuracy and the speed of convergence by selecting an appropriate value of [30, 32, 33]. We used the value of 0.1 Å-1 for in this study, with the resulting error of Ewald sum energy less than 0.5%.
3. Simulation Results
In order to carry out a physically sound MD simulation for the clay mineral system, a parametric study has been made to determine the time step t. It is found that the time step essentially controls the convergence of the system and a value of sec for t is determined to be appropriate. In the following sections, we present results of MD simulations of dehydrated Na-mica under structural relaxation under different normal pressures. It is our major interest to study the effects of cell boundary constraints on the model behavior, that is, whether the cell size or shape is allowed to change. Results from static simulations are also presented for a comparison study.
3.1. Relaxation Behavior with Cell Size Change
A first MD simulation of the dehydrated mica relaxed under zero external stress (i.e., ) is carried out, where only the cell size change is allowed (i.e., no shear distortion is allowed, , for ). The initial interlayer spacing () is 9.5 Å. The time variations of the relaxation behavior are shown in Figure 2. The system starts from the initial atom positions shown in Table 1 with zero velocities and evolves according to the dictates of the equations of motion of the system. Driven by the large imbalance between the initial internal stress () and zero external stress (), the MD cell starts to distort away from its initial shape. This is reflected in the changes of diagonal terms of the cell matrix (Figure 2(d)). The nondiagonal terms are equal to zero due to the prohibition of any shear distortion. Structural equilibrium is deemed to be achieved after about 20 000 time steps, where the average particle momentums in the three Cartesian directions start to fluctuate around a mean value of kg·m/sec (Figure 2(b)). Time evolutions of the average unbalanced interparticle forces are shown in Figure 2(a). Because of the effects of thermal motions, these average unbalanced interparticle forces do not diminish completely but fluctuate around certain values under the equilibrium condition. As can be readily understood and inferred, in a static simulation where all thermal motions are removed, these unbalanced interparticle forces should approach zero upon equilibrium. This will be verified later.
(a)
(b)
(c)
(d)
The evolutions of the internal stress tensor shown in Figure 2(c) indicate that all the diagonal terms reduce remarkably from high initial values and end up with fluctuations around low values near zero. These low-amplitude, nonzero fluctuations at equilibrium, are again attributed to the thermal motions of the system and will likewise disappear in a static simulation. The nondiagonal stress terms, however, retain almost constant values throughout the whole simulation. This is particularly the case for (and ), which suggests a large shear stress in the a-b plane. The existence of these deviatoric stress terms is mainly due to the restraint from shear distortions of the MD cell.
2D and 3D views of the MD cell structure with all the atom positions at equilibrium are illustrated in Figure 3. It is seen that relaxation under zero external stresses results in moderate rearrangement of atom positions under the equilibrium condition. The average values of , , and during the last 5 000 time steps are 2.86 nm, 1.94 nm and 3.55 nm, respectively, indicating that the MD cell undergoes expansions in all the three axial directions. The final interlayer spacing is estimated to be about 7.5 Å.
(a)
(b)
(c)
(d)
3.2. Relaxation Behavior with Both Cell Size and Shape Changes
Figure 4 shows the results of an MD run in which full freedom is given to the MD cell size and shape changes. All the other conditions remain the same as the last simulation. It is noted that two distinct features appear that are different from the last simulation: (2) shear distortions of the MD cell take place as indicated by the gradual increases of nondiagonal terms of the h matrix, with the largest shear distortion occurring in the sheet plane (Figure 4(d)); (2) all the terms of the internal stress tensor , including deviatoric stress terms, tend to approach zero (with mild fluctuations) at the equilibrium state. The above two features, which are predicted by the current MD method for the first time, suggest the inherent correlations between the boundary constraint and external stress conditions imposed on the MD cell, and resulted cell deformation and internal stress conditions. The time averages of the matrix h during the last 10 000 time steps are: based on which we calculated the angles between the otherwise orthogonal axes to be 70 (±3), 85 (±2), and 85 (±2). It is found that the time averages of , and are very close to those in the last simulation, indicating that there is little correlation between the axial and shear deformations of the MD cell under zero external pressures.
(a)
(b)
(c)
(d)
The equilibrated MD cell structure in the simulation is shown in Figure 5. Overall, it is seen that more atom positional rearrangements are obtained as the system moves along its trajectory towards equilibrium. Shear distortions, especially in the a-b plane, are distinctly visualized in the corresponding 2D views (Figures 5(a)–5(c)). The interlayer spacing, when examined carefully, is found to be slightly less than that in the last simulation (i.e., by comparing Figure 5(c) and Figure 3(c)).
(a)
(b)
(c)
(d)
3.3. Results of Static Simulations
The two simulations described in Sections 3.1 and 3.2 are recalculated using the static version of the method, which provides better pictures of the relationship between boundary constraints and resulted internal stresses, because all the thermal fluctuations are removed. However, the observed behavior is only valid at zero Kelvin. The results of the two static simulations are shown in Figures 6 and 7, respectively.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
It is clear that with the removal of thermal fluctuations, all the curves become much more smooth and well defined in terms of the trend of convergence. The average and maximum unbalanced interparticle forces, which are now used as major criteria for judgment of convergence, monotonically approach zero towards equilibrium (Figures 6(a) and 6(b) and Figures 7(a) and 7(b)). It is more interesting to examine the evolutions of the internal stress tensor and cell matrix h. Firstly, for the case where merely cell size change is allowed, it is clear that all the diagonal stress terms have been reduced to zero as a result of cell dimension changes while all the nondiagonal stress terms retain nearly constant, nonzero values at equilibrium state (Figure 6(c)) due to the constraint of cell shear distortion. Particularly, we note the high shear stress (and ) = 1.7 GPa in the sheet plane, which is about 55% higher than the value in the corresponding MD run. Final values of cell dimensions are = 2.46 nm, = 1.81 nm, and = 3.17 nm (Figure 6(d)), which are also smaller than those in the MD run. These differences suggest that the system becomes more “compacted” in an idealized static condition.
In the second case where full freedom is given to the cell size and shape change, it is interesting to observe that the entire tensor decays to zero after about 16 000 iterations (Figure 7(c)), and in order to allow that to happen, remarkable changes have occurred to the matrix h (Figure 7(d)): (2) and undergo a significant increase and decrease to reach a value of 3.36 nm and 2.24 nm, respectively, with the latter essentially eliminating any existing interlayer spacing (see Figure 8(a)); (2) and undergo significant increases to reach an approximate value of 9.8 Å, implying a much larger shear distortion in the sheet plane (see Figure 8(b)). Changes in and , and and are much smaller and close to each other in each symmetric pair. As compared to the MD runs, it is quite remarkable to observe these severe crystal transformations of the clay mineral that are required to take place to dissipate all the internal stress terms to reach an equilibrium status, as predicted by static simulations. However, it should be noted that these results only represent idealized behavior at zero temperature and should be cautiously interpreted with regard to their implications to more realistic material behavior at finite temperatures.
(a)
(b)
3.4. Effects of Normal Pressures
In addition to the relaxation behavior of clay minerals under zero normal pressure, we are more interested to find out how the interlayer spacing at equilibrium will vary with the normal pressure placed on the sheet plane. Figure 9 shows the interlayer spacing-versus-normal pressure relationships from two monotonic axial compression simulations, where the cell shear distortion is allowed and prohibited in the first and second simulations, respectively. The axial compression starts with the fully relaxed structure under zero normal pressure and applies an incremental axial load of 0.2 GPa per load step. The interlayer spacing is recorded after the structure achieves an equilibrium in each step. It is seen that the interlayer spacing decreases approximately in a linear fashion with the normal pressure within the whole pressure range when no cell shear distortion is allowed (Figure 9(b)). For the case with cell shear distortion, a nonlinear, slower decrease of interlayer spacing is found in the low pressure range of 0~4 GPa, which is followed by a higher rate, linear reduction in the high pressure range (Figure 9(a)). These results clearly indicate that the compressibility of dehydrated mice sheets in the [] direction is chiefly due to the existing interlayer spacing.
(a)
(b)
(c)
The above results agree well with the previous experimental observations on monoclinic dioctahedral micas by Gatta et al. [14] and Zanazzi and Pavese [34]. To further demonstrate this deformation mechanism and validate the current MD model, we compare our simulation results with the K-polyhedron volume-versus-hydrostatic pressure relationship observed recently by Gatta et al. [14] in an in situ single-crystal X-ray diffraction study on a 3T-phengite mica, as shown in Figure 9(c). In this experimental study, Gatta et al. [14] found that the compressibility of the phengite mica was mainly attributed to the reduction of K-polyhedron volume, which is equivalent to the reduction of interlayer spacing observed in the current study. Although the chemical formula of the 3T-phengite mica ((()()()) is much more complicated than that of the Na-saturated mica used in the current study, it is clear in Figure 9 that the same structural deformation mechanism is displayed for micas subject to hydrostatic or normal compressions.
Figure 10 shows the result of the same axial compression simulation without cell shear distortion using the static method. Interestingly, a drastically different behavior featured with an instantaneous drop of the interlayer spacing from 9.3 Å to 4.1 Å at a normal pressure of about 1.5 GPa is observed. Before and after this collapsing behavior, the rate of compression is very low and keeps a constant within the whole pressure range. At zero temperature, such a brittle collapse behavior can be readily understood because the energy barrier, in the absence of any thermal agitation, is solely a function of the positional variables of the system. The very low rate of compression other than the collapsing drop is mainly attributed to the volume reduction of the individual mica sheets, which possess a very high compressive modulus.
4. Concluding Remarks
In this paper, we have applied the Parrinello-Rahman constant-pressure molecular dynamics method to the study of dehydrated mica, in which unforeseen patterns of crystal structural transformation of mica relaxed under constant normal pressures have been uncovered. We are particularly impressed to see the large shear distortion of the MD cell and the accompanied disappearance of deviatoric stress terms under zero boundary constraint (or its opposite scenario). Under monotonic axial loading conditions, an approximate linear reduction of interlayer spacing with the applied normal pressure is observed from an MD simulation while a brittle-type collapse of the interlayer spacing at a certain normal pressure is yielded by a static simulation. These results reveal the inherent relationships between the external applied and internal generated stress tensors and imposed cell boundary conditions and attest exactly to the capability of the method to predict phase transitions of solids in relation to particle interactions.
We have noticed the differences of results from the static method and the full MD method. The differences are attributed to the thermal agitations of the system, without which a much more well-defined correlation of stress and structural deformation is observed. It is particularly noted that the interlayer spacing of mica is completely eliminated when the maximum degrees of freedom of the MD cell are given in a static simulation. This behavior, however, may not be valid in a finite temperature condition, as predicted otherwise by an MD simulation. The existence of interlayer spacing, as shown by both methods, has a close relationship with the compressive modulus behavior of clay mineral.
The successful application of the Parrinello-Rahman constant-pressure MD method to a model phyllosilicate system which is governed by multiple pair interactions between different species of atoms (instead of a single pair interaction among uniform atoms in metals) demonstrates that the method is a robust and a general one that is not limited by the type of the interatomic potential, but could be applied to many kinds of solids. Many other applications are conceived to be possible, such as the studies of monotonic/cyclic loading of clay minerals, the interactions and properties of the clay-water-cation system, or the nanoindentation of clay minerals. These proposed topics are currently under investigation.
Acknowledgments
The study presented in this paper was supported by the Strategic Research Grant of City University of Hong Kong under Project no. 7002557 and the Petroleum Research Fund, American Chemical Society under Grant no. ACS PRF# 43848-AC8. These supports are gratefully acknowledged.