Molecular Dynamics Simulation of Layerd Double Hydroxide (LDH) Mineral: Hydrotalcite and Water Systems

Abstract

Some clay minerals are used and expected as barrier materials for engineering and in nature. Hydrotalcite, Mg₆Al₂[(OH)₁₆|CO₃]4H₂O and a kind of LDH (layered double hydroxides), is one of the most effective candidates for the anion adsorbents and the barrier. In this study, the behavior of hydrotalcite was investigated by means of the molecular dynamics method. Cl⁻ and I-hydrotalcite − water systems were simulated for various mineral/water ratios. The structure and dynamic properties are predicted. Water at the surface of hydrotalcite shows the electric double layer composed of Stern layer of one H₂O molecular layer thickness and large self-diffusion coefficient of H₂O and diffusion layer of 2.5 nm thickness at the interface.

1 Introduction

Some clay minerals have been used and are expected to be used as barrier materials against cations, anions, water and aqueous solutions, gasses, etc. for engineering systems and in nature. Smectite is one of the most effective materials as cation adsorbents, water sealing, etc. in the potential engineering barriers for industrial and radioactive waste disposals, etc. On the other hand, anion mobility seems to be fairly large and the absorbability of anions is poor in smectite clays. Hydrotalcite, Mg₆Al₂[(OH)₁₆|CO₃]4H₂O, is one of the most effective candidates for the anion barriers.

Hydrotalcite is a mineral relating to simple hydroxide minerals, brucite (Mg(OH)₂) and gibbsite (Al(OH)₃) having a sheet structure composed of edge-sharing octahedra. Hydrotalcite and related minerals with only octahedral cations of Mg²⁺ and/or Al³⁺ are listed in Table 1 according to Struntz Mineralogical Tables [1]. The octahedral sites are occupied mainly by Mg²⁺ ions and some by Al³⁺ ions, and the substitution of divalent cations by trivalent cations causes the positive layer charges. Because of the solid solution of Al³⁺ and Mg²⁺ ions in the octahedral sites, and uncertainty of cation distribution in the octahedral sites, and of species of anions or anionic molecules and their amounts, and of amount of H₂O molecules in the interlayer region; the crystal chemistry is not known completely yet.

Table 1.  Hydrotalcite and related minerals. ([1] Strunts Mineralogical Tables, 2001)

Mineral	Chemical formula		System		S.G.		(a,b,c (Angstrom),β(°))		Z
Brucite	Mg (OH)2		Trigonal		P-3m1		(3.15 3.15 4.77)		1
Gibbsite	γ-Al (OH)3		Monoclinic		P21/n		(8.66 5.07 9.72 94.5)		8
Bayerite	α-Al (OH)3		Monoclinic		P21/n		(5.05 8.67 9.42 90.3)		8
Meixnerite	Mg6Al2(OH)18.4H2O		Trigonal		R-3m		(3.05 3.05 22.93)		3/8
Hydrotalcite	Mg6Al2[(OH)16|CO3]4H2O		Trigonal		R-3m		(3.05 3.05 22.81)		3/8
Manasseite	Mg6Al2[(OH)16|CO3]4H2O		Hexagonal		P63/mmc		(6.12 6.12 15.34) 1
Quintinite-2H	Mg4Al2[(OH)12|CO3]3H2O		Hexagonal		P6322		(10.57 10.57 15.14)		4
Quintinite-3T	Mg4Al2[(OH)12|CO3]3H2O		Trigonal		P3112		(10.80 10.80 22.71)		6
MD calculation results at 298.15 K, 0.1 MPa
	a/A	b/A (b/√3)		c/A		alpha/°		beta/°		gamma/°
Brucite	2.99836	5.18323 (2.99254)		4.88296
Gibbsite	8.82789	5.11201		9.97155		90.0076		98.5055		89.9978
Hydrotalcite	3.06379	5.27381 (3.04484)		15.25512

Molecular dynamics simulations (MD) of halide bearing hydrotalcite have been carried out by some researchers. Wang et al. [2] reported the MD simulations of hydrotalcite with interlayer Cl⁻ and H₂O using CVFF_aug force field and SPC model. The water contents were 0 to 7 of H₂O/Cl. The swelling behavior as c-axis expansion and the hydration energy and H₂O potential energy were derived. Their simulations showed the stability of the hydration state at H₂O/Cl ratio of two. Wang et al. [3] simulated hydrotalcite, [Mg₃Al(OH)₃]Cl.3H₂O, to obtaine the power spectra using molecular dynamics methods with their modified CRAYFF force field and flexible SPC model. They calculated the power spectra of 0 to 600 cm⁻¹ to compare with experimental far-infrared spectra, and found good agreement between them. They mentioned the motion of Cl⁻ion in the interlayer region as remarkably similar to those of bulk aqueous chloride solutions.

In this paper, the preparations of initial structures of hydrotalcite for molecular dynamics (MD) simulations of hydrotalcite related systems were developed. The structure and properties of hydrotalcite with some different anions, mainly Cl⁻ and I⁻, and water in the vicinity of the surface were investigated by means of MD method.

2 Molecular dynamics simulations and interatomic potential model

The Ewald method was used for the summations of Coulomb interactions. The cut off distance of summations of long and middle range interactions was set at the smaller one of 1.5 nm or the inscribed sphere radius of the periodic basic cell. The short range interactions were summed within the radius of 0.75 nm. Integration of equation of atom motions was performed by the velocity Verlet algorithm with a time increment of 0.4 fs because of hydrogen atoms. The NVT (in some of the cases of including vacuum space in the systems) and NPT ensembles were employed, where N is the number of atoms in the simulation (basic) cell, V the basic cell volume, T temperature, and P pressure. Pressure and temperature with the number of atoms were kept constant during the simulations with NPT emsembles. Temperature and pressure were controlled by scaling of atom velocities and basic cell edges and the shape. Molecular dynamics (MD) simulations were carried out at 293 K and 0.1 MPa. The MD simulation codes, MXDORTO and MXDTRICL [4] were used. At least, a 100,000 steps calculation was performed for the initial relaxation for each system. Subsequent 100,000 to 2,000,000 steps simulations for each system were carried out to obtain ensemble averaged properties.

Before the MD calculations of surface-water drop systems, wide areas of hydrotalcite surfaces were prepared, and water drops in vacuum were prepared separately and equilibrated for 3000 H₂O molecule systems. MD calculations were started from the system containing the sheet of clay mineral and a water drop. The clay surface-external water systems were also simulated to obtain the interface properties.

The interatomic potential model used in this study is the pair atom central force model with full freedom of atom motions which was described by Kumagai et al. [5], Nakano et al. [6], and Kawamura [7]. The three body force term was added for H₂O molecule. The two body terms represent Coulomb, van der Waarls, non-bonding short-range repulsive, and radial covalent terms (three terms) in the sequence of the following formula:   

$\begin{array}{c}{u}_{ij}\left(r_{ij}\right)=\frac{z_i{z}_j{e}^{\text{2}}}{\text{4}\pi {\epsilon }_{\text{0}}{r}_{ij}}\\ -\frac{c_i{c}_j}{r_{ij}^{\text{6}}}\\ +f_{\text{0}}\left(b_i+b_j\right)\text{exp}\left(\frac{a_i+a_j-r_{ij}}{b_i+b_j}\right)\\ +D_{1ij}\text{exp}\left(-{\beta }_{\text{1}ij}{r}_{ij}\right)+D_{2ij}\text{exp}\left(-{\beta }_{\text{2}ij}{r}_{ij}\right)+D_{3ij}\text{exp}\left[-{\beta }_{\text{3}ij}{\left(r_{ij}-r_{\text{3}ij}\right)}^{\text{2}}\right]\end{array}$

The three body term represents the angular part of covalent interaction;   

$\begin{array}{c}{u}_{jij}\left({\theta }_{jij},r_{ij}\right)=-f_k\left\{\text{cos}\left[\text{2}\left({\theta }_{jij}-{\theta }_{\text{0}}\right)\right]-\text{1}\right\}\sqrt{k_{\text{1}}{k}_{\text{2}}}\\ k_{\text{1}}=\frac{\text{1}}{\text{exp}\left[g_r\left(r_{ij\text{1}}-r_m\right)\right]+\text{1}}\end{array}$

The parameters appearing in this formula, z, c, a, and b, for atoms, D₁, D₂, D₃, β₁, β₂, β₃, and r₃ for atom pairs, f_k, θ₀, g_r, and r_m for three atoms, were presented for the molecular simulations of systems relating to H₂O, CO₃^2-, clay minerals, and related crystalline materials (Table 2).

Table 2.  Parameters of interatomic potential model for MXDORTO/MXDTRICL

Hydrotalcite sheet
Atom	w/103kg/mol−1	z/e	a/nm	b/nm	c/(kJ/mol)0.5nm−3
O	16.00	-1.19375	0.1804	0.0151	0.05606
Al	26.98	1.950	0.0958	0.0056	0.0
Mg	24.305	1.640	0.1049	0.0057	0.00409
H	1.01	0.460	0.0123	0.0042	0.0
Cl	35.45	−1.000	0.2200	0.0160	0.05729
I	126.90	−1.000	0.2350	0.0160	0.06138
Atom-atom	D1/kJ.mol−1	β1/nm−1	D2/kJ.mol−1	β2/nm−1	D3/kJ.mol−1	β3/nm−1	r3/nm
O-Al	151515.0	50.0	−8145.0	22.6
O-H	57394.9	74.0	−2189.3	31.30	34.74	128.0	0.1283
Atom-atom-atom	fk/10−19J	θ0/°	rm/nm	gr/nm−1
H-O-H	1.15	99.5	0.143	92.0
H2O molecule
Atom	w/103kg/mol−1	z/e	a/nm	b/nm	c/(kJ/mol)0.5nm−3
O	16.00	−0.920	0.1728	0.01275	0.05606
H	1.01	0.460	0.0035	0.00440	0.0
Atom-atom	D1/kJ.mol−1	β1/nm−1	D2/kJ.mol−1	β2/nm−1	D3/kJ.mol−1	β3/nm−1	r3/nm
O-H	57394.9	74.0	−2189.3	31.3	34.74	128.0	0.1283
Atom-atom-atom	fk/10−19J	θ0/°	rm/nm	gr/nm−1
H-O-H	1.15	99.5	0.143	92.0
CO32- ion molecule
Atom	w/103kg/mol−1	z/e	a/nm	b/nm	c/(kJ/mol)0.5nm−3
C	12.01	1.240	0.0485	0.0067	0.000
O	16.00	−1.080	0.1879	0.0169	0.04501
Atom-atom	D1/kJ.mol−1	β1/nm−1	D2/kJ.mol−1	β2/nm
C-O	33484.0	60.0	−12138.0	24.0

Using therse parameters, the lattice parameters and crystal structures of brucite, gibbsite, and hydrotalcite were reasonably reproduced as shown in Table 1.

3 Crystal structure of hydrotalcite

Although with the rather simple chemical formula of hydrotalcite, Mg₆Al₂[(OH)₁₆|CO₃]4H₂O, the relation between the crystal structure and the chemistry shows complex features. The crystal data of hydrotalcite by Allmann and Jepsen [8] and by Bellotto et al. [9] are shown in Table 3. The chemical formulas of these two reports are different from each other. Both of the crystal structure data have site occupancy factors (SOF) for Al, Mg, CO₃^2-, and H₂O, and there is no description of their distributions. These structure data cannot be used for MD calculations as they stand, because the "true" atom positions in the unit cell should be given properly to perform MD simulation.

Table 3.  The crystal data of hydrotalcite reported by Allmann and Jepsen (1969) [6] and Belloto et al. (1996) [7]. SOF: site occupancy factor, ITF: isometric temperature factor.

((Mg4Al2)(OH)12(CO3)(H2O)3)0.5	Allmann and Jepsen (1969)
Trigonal	R-3m	a = 3.054, b = 3.054, c = 22.81, = 120°			Z = 3
Atom	SITE	x	y	z	SOF
Mg	3a	0	0	0	0.6667
Al	3a	0	0	0	0.3333
O	6c	0	0	0.3771 (5)	1.
H	6c	0	0	0.4144 (86)	1.
O	18h	0.092 (25)	−0.092 (25)	0.5	0.167
C	6c	0	0	0.167	0.083
(Mg0.833Al0.167)(OH) 2 (CO3)0.083(H2O)0.75 Bellotto et al. (1996)
Trigonal	R-3m	a = 3.0808 (3), b = 3.0808 (3), 23.784 (4), γ = 120°				Z = 3
Atom	SITE	x	y	z	SOF	ITF (U)
Mg	3a	0	0	0	0.8333	0.0319 (8)
Al	3a	0	0	0	0.1667	0.0319 (8)
O	6c	0	0	0.3754 (2)	1.	0.036 (1)
H	6c	0	0	0.432 (1)	1.	0.059 (7)
O	18h	0.1260 (9)	0.8740 (9)	0.5	0.1667	0.021 (3)
C	6c	0.3333	0.6667	0.5	0.0416	0.059 (7)
H	6c	0.3333	0.6667	0.5	0.75	0.059 (7)
MD simulation of Mg6Al2 (OH) 16CO3.4H2O: a = 3.0643A b = 3.0446A C = 20.3532A

When Mg:Al = 3:1, the octahedral sheet may have ordered structure. Figure 1 in which the a-edge length was doubled, and the b-edge length was changed because of converting from the trigonal cell into orthorhombic cell, shows one of the ordered models (Al atoms at sites 1 and 2), and another model is that of Al atoms at sites 1 and 3 under the condition that Al atoms do not adjoin each other. When Mg:Al = 2:1, the structure should include adjoining of Al octahedra.

Figure 1.

 Unit structure of hydrotalcite as octahedral sheet.

In our model, the unit cell (R-3, Z = 3 for Mg₆Al₂[(OH)₁₆|CO₃]4H₂O) [8], (a,b,c) = 0.61 nm × 1.06 nm × 2.28 nm was multiplied by (4,2,1) as the MD basic cell. We calculated for three types of hydrotalcite, CO₃^2-, Cl⁻, and I-hydrotalcites.

The (NPT)-MD calculations of model CO₃^2--hydrotalcite at 293 K and 0.1MPa, show the average unit cell parameters of a = 3.069, b = 3.042, and c = 23.218Å(Angstrom) giving the density of 2.00362 × 10³ kg/m³ which is well reproducing the reference parameters. The structure snapshot is shown in Figure 2. The octahedral sheet structure (3 layers parallel to the a-b plane) is stably maintained, and interlayer sheet also formed with CO₃^2- and H₂O molecules. The single sheet of interlayer structure projected onto the a-b plane is shown in Figure 3. While the CO₃^2- units line regularly in the Figure, the arrangement of H₂O molecules is irregular.

Figure 2.

 Snapshot of structure of hydrotalcite, Mg₆Al₂[(OH)₁₆|CO₃]4H₂O by MD simulation. A portion surrounded by dashed line means an octahedral sheet.

Figure 3.

 The structure of interlayer CO₃^2- and H₂O from middle layer of Figure 2.

We investigated the cases that relatively small amounts of H₂O were inserted in interlayer regions of hydrotalcite. Mg₃Al (OH)₈[Cl or I].nH₂O where n = 1 to 20 are investigated by MD calculations (the structures of n = 1 to 4 are shown in Figure 4). The structures of these two series are slightly different because of the difference of anion size. In the structure of n = 2 and 3 in I-hydrotalcite, the rows of H₂O molecules are clearly separated in the single interlayer region. Two dimensional diffusion coefficients D_xy of anions and H₂O molecules were calculated (Table 4). The diffusion coefficients of Cl⁻hydrotalcite are smaller than those of I-hydrotalcite. The diffusion coefficient of H₂O and halide ions in hydrotalcite interlayer is extremly small for 1H₂O of Cl⁻hydrotalcite. The diffusion coefficients increase with n for both anions and H₂O molecules. In the systems of 10H₂O and 20H₂O, the diffusion coefficients of H₂O molecules are almost the same and the same with bulk water.

Figure 4.

 Snapshots of Mg₃Al (OH)₈.Cl.nH₂O (left) and Mg₃Al (OH)₈.I.nH₂O (right).

Table 4.  Two dimension diffusion coefficients D_xy of anion and H₂O molecule in interlayer region of Cl- and I-hydrotalcite at 293 K, 0.1 MPa. The nH₂O means Mg₃Al (OH)₈[Cl or I].nH₂O

nH2O	Diffusion coefficient Dxy/cm2/s
	halide ion		H2O
1H2O	Cl	2.7E-7	9.1E-7
	I	17.7E-7	45.3E-7
2H2O	Cl	42.9E-7	69.1E-7
	I	48.6E-7	74.7E-7
3H2O	Cl	26.3E-7	47.0E-7
	I	52.9E-7	134.1E-7
10H2O	Cl	94.2E-7	165.2E-7
	I	88.1E-7	179.6E-7
20H2O	Cl	108.3E-7	162.8E-7
	I	91.4E-7	180.6E-7

4 Wetting of water on hydrotalcitesurface

Hydrotalcite surface − water systems were investigated by MD simulations to obtain the wetting behavior, the local properties, structure, density, diffusion coefficient, and viscosity of water as functions of distance from hydrotalcite surface.

Wetting of water on the Cl⁻ and I-hydrotalcite surface was investigated using MD calculations. A water drop was positioned initially just in the vicinity of the surface with the whole drop velocity at zero (Figure 5), then MD calculations were started. The water drop is composed of 3000 H₂O molecules and well relaxed previously. The structural snap shots of steady state related structures are displayed in Figure 6. Both Cl⁻ and I-hydrotalcite show very good wetting behavior where wet angles are well less than 30 degrees. In both structures, most of the anions attach to the hydrotalcite sheet surfaces directly and only a small number of anions are fully hydrated unlike the case of Na-smectite surfaces where the surface has very good wetting properly and the all cations are fully hydrated [10].

Figure 5.

 Initial structure of Cl-hydrotalcite surface with A water drop (3000 H₂O molecules) system for MD calculations.

Figure 6.

 Snapshots of a sheet of Cl-(upper) and I-hydrotalcite (lower) surface − water drop systems, where the number of H₂O molecules are 3000 for both systems.

5 Local structure and properties ofwater between hydrotalcite surfaces

As the hydrotalcite − intergranular (external or outside interlayer region) water systems, we simulated Mg₆Al₂(OH)₁₆[Cl or I]-nH₂O, n = 40, 200 and 400 systems. The results for Mg₆Al₂(OH)₁₆[Cl or I]₂-200H₂O system are shown in Figures 7 and 8. Local structure and properties were calculated every 0.25 nm thick layer of the water parallel to the hydrotalcite surface whose thickness is nearly equal to the diameter of an H₂O molecule. The atom distributions perpendicular to hydrotalcite surface are displayed in Figure 7. Cl⁻ ions distribute only near the mineral layer surfaces. The orientations of H₂O molecular dipoles are strongly biased at the interfaces and the directions are clearly from the mineral layer surface to inside water (Figure 8, middle). The electrical double layer was formed as the Stern layer of one molecular layer of surface water where the molecular orientations are strongly biased and the diffusion layer of about 2.5 nm thickness where Cl⁻ ion distribute there (Figure 7). The I-hydrotalcite shows the same structure in our calculations. The density profile of water shows the decrease at the surface two molecular layers of water unlike the case of Na-smectite [8]. Diffusion coefficients of water show the maxima at the surfaces, and decrease with distance from the surfaces to three H₂O molecule distance, and increase slightly again with distance in the region where Cl⁻ ions exist. The viscosity shows the reverse behavior to that of diffusion coefficient. This behavior was reported for brucite surface (Mg (OH)₂) by Sakuma et al. [11]. It is, however, more prominent for hydrotalcite surface than for brucite surface. The diffusion coefficient of H₂O is larger in Cl⁻hydrotalcite than in I-hydrotalcite and the viscosity of external (pore) water of Cl⁻hydrotalcite is smaller than that of I-talcite as shown in Figure 9. This can be explained by the difference of ionic radii of halide ions. Namely, the large anion I⁻ has weaker interaction between halide ion and H₂O molecule than Cl⁻ ion.

Figure 7.

 Distribution of atoms vertical to Cl-hydrotalcite surface in the hydrotalcite-200H₂O system.

Figure 8.

 Results of MD calculations of hydrotalcite Cl-Mg₆Al₂(OH)₁₆-200 H₂O system at 293 K and 0.1 MPa. Structure snapshot (upper), H₂O molecular orientation statistics and distribution of H₂O molecules and Cl⁻ ions (middle) where theta means angles between dipoles of H₂O molecules and c-axis and phi means angles between the normal vector of H₂O molecular plane and c-axis, and diffusion coefficient of H₂O (), viscosity of water (), and distribution of H₂O molecules () (bottom). The vertical axis is shown in the left side of the Figure with the units of 1 × 10⁻⁵ cm/s for diffusion coefficient and 1 × 10⁻³ for viscosity respectively.

Figure 9.

 Diffusion coefficient of H₂O (), viscosity (), and density of H₂O () of hydorotalcite, X-Mg₆Al₂(OH)₁₆-40H₂O (X=Cl (left) and I (right)) systems.

6 Conclusions

The crystal structure of hydrotalcite was reasonably modeled and reproduced using molecular dynamics simulations with the appropriate interatomic potential model of hydrotalcite sheet, H₂O molecule, halide anions, and CO₃^2- ion. Wetting of water drop on hydrotalcite surfaces show fairly good wetting behavior with wetting angle well less than 30 degrees. We performed hydrotalcite − intergranular (pore) water systems, and calculated local density, diffusion coefficients, and viscosity of water. The behavior of hydrotalcite surface − external water resembles that of brucite surface. But wetting property of hydrotalcite is much better than that of brucite. The decrease of the viscosity of water at the interface between mineral suface and pore water is more prominent in hydrotalcite than in brucite.

References

• 1   H. Struntz, E. H. Nickel, "Struntz Mineralogical Tables. Chemical-Structural Mineral Classification System", 9th ed. Stuttgart, Schweizerbart (2001)
• 2   J. Wang, A. G. Kalinchev, R. J. Kirkpatrick, X. Hou, Chem. Mater., 13, 145 (2001).
• 3   J. Wang, A. G. Kalinichev, J. E. Amonette, R. J. Kirkpatrick, Am. Mineral., 88, 398 (2003).
• 4   K. Hirao, K. Kawamura, "Material Design using Personal Computers" (Japanese), Shokabou, 217p. (1994)
• 5   Kumagai, N., Kawamura, K. and Yokokawa, T., 12, pp. 177- 186 (1994).
• 6   M. Nakano, K. Kawamura, Y. Ichikawa, Appl. Clay Sci., 23, 15 (2003).
• 7   K. Kawamura, Low Temp. Sci., 65, 3 (2005).
• 8   Allmann, R., Jepsen, H.P, . Neues Jahrbuch fuer Mineralogie. Monatshefte (Band=Jahr), 544–551 (1969)
• 9   M. Bellotto, B. Rebours, O. Clause, J. Lynch, D. Bazin, E. Elkaim, J. Phys. Chem., 100, 8527 (1996).
• 10   K. Kawamura, Chikyuukagaku, 42, 115 (2008).
• 11   H. Sakuma, T. Tsuchiya, K. Kawamura, K. Otsuki, Surf. Sci., 536, L396 (2003).