Molecular Modeling of the Structure and Dynamics of the Interlayer

Jun 10, 2008 - Laure Aimoz , Christine Taviot-Guého , Sergey V. Churakov , Marina Chukalina , Rainer Dähn , Enzo Curti , Pierre Bordet , and Marika ...
0 downloads 0 Views 1MB Size
Download PDF
7856

J. Phys. Chem. B 2008, 112, 7856–7864

Molecular Modeling of the Structure and Dynamics of the Interlayer Species of ZnAlCl Layered Double Hydroxide J. Pisson,† J. P. Morel,† N. Morel-Desrosiers,† C. Taviot-Gue´ho,‡ and P. Malfreyt*,† Laboratoire de Thermodynamique des Solutions et des Polymeres, UMR CNRS 6003, and Laboratoire des Mate´riaux Inorganiques, UMR CNRS 6002, UniVersite´ Blaise Pascal (Clermont-Ferrand II), 24 aVenue des Landais, 63177 Aubiere Cedex, France ReceiVed: January 21, 2008; ReVised Manuscript ReceiVed: March 17, 2008

Molecular dynamics simulations of the ZnAl layered double hydroxide containing interlayer chloride anions have been performed in the NpT and NpzzT statistical ensembles for metal Zn/Al ratios of 2 and 3. We have monitored the interlayer spacing as a function of the number of intercalated water molecules for each statistical ensemble. We have studied how these profiles are affected by the method of calculation of the charges of the hydroxide layer atoms. Diffusion coefficients of the interlayer water molecules have been calculated for different Zn/Al ratios. The calculation of the chemical potential of the interlayer water molecules has been carried out for three amounts of interlayer water molecules. The calculation showed a qualitative agreement with the bulk water chemical potential within a range of interlayer water molecule contents. 1. Introduction Layered double hydroxides (LDHs), also known as anionic clays, correspond to a class of intercalation compounds with positively charged brucite-like hydroxide layers containing divalent and trivalent metal cations. The charge neutrality is ensured by interlayer anions that are easily exchangeable. The general formula of the LDHs is [M1-xIIMxIII(OH) 2]x+(Ax/ nn-) · mH2O, abbreviated hereafter as M RIIMIIIA, in which MII ) Mg, Zn, Co, Ni, Mn, etc., MIII ) Cr, Fe, V, Co, etc., An- is an inorganic or organic anion of charge n, and R is the MII/MIII molar ratio. The wide choice of composition for both the metal cations and the interlayer anions gives access to numerous LDH compositions, purely inorganic or hybrid phases. Consequently, LDHs are receiving increasing attention for a wide variety of applications1 such as biosensors,2,3 flame retardant materials,4 drug support,5,6 catalysts,7 and optical devices.8 Understanding the properties and exploiting the chemistry of LDHs require a detailed knowledge of the structure. However, because of the inherent structural disorder and generally low crystallinity of LDHs, complete structure determination using X-ray diffraction is rarely possible. As a consequence, molecular modeling methods have been increasingly used in the past decade for understanding the LDH structure9–19 and the dynamical behavior of interlayer water molecules and anions. The application of any computational molecular modeling technique requires the use of interatomic potentials that accurately account for the interactions of all atoms in LDHs and leads to addressing fundamental methodological questions concerning the choice of the statistical ensemble (NpT, NpzzT, and µVT) as well as the method of calculation for the atomic charges of the hydroxide layers. The ability of reproducing experimental properties depends on the care taken to resolve these methodological issues. It is of main interest to establish * Corresponding author. E-mail: [email protected]. † Laboratoire de Thermodynamique des Solutions et des Polymeres, UMR CNRS 6003. ‡ Laboratoire des Mate ´ riaux Inorganiques, UMR CNRS 6002.

an appropriate methodology before a more developed study of the LDHs with organic intercalated anions. The most direct way of predicting the interlayer spacing is to perform Monte Carlo simulations in the open µVTgrand canonical ensemble, where µ is the chemical potential of the water molecules in the interlayer regions. In this ensemble, the confined water molecules are in equilibrium with those in the reservoir. The fact that the chemical potential is imposed allows the number of water molecules to fluctuate during the simulations. This type of simulation does predict directly the layer spacing and the number of intercalated water molecules at the end of the simulation. However, this method is very different from the standard simulations where the composition of the layer spacing is fixed during the run. Simulations in the grand canonical ensemble require the use of configurational-bias methods to efficiently sample the number of intercalated water molecules because of the fact that the water molecules are confined in a small region and that the predominant interactions come from the electrostatic contributions. However, this methodology is not designed to study the dynamic properties of the intercalated molecules such as the diffusion coefficients. For this reason, this alternative was not retained here and we preferred to carry out molecular dynamics (MD) simulations in the NpT ensemble by changing the number of water molecules at each run. As a result, we monitored the interlayer distance as a function of the hydration state. This method has been successfully applied in recent works.9–19 A brief scan of the literature has shown that calculation of the charges of the metal atoms of LDH layers has been performed using two routes. One9,10,14,15,17,18 consisted of calculating the charges from an equilibration method (QEq)20 and taking into account the redistribution of the charges on metal ions among the surrounding ions. An alternative route11 used the formal full charge for the metal atoms of the hydroxide layers. From molecular simulations of clay minerals such as montmorillonite, it has been observed that the NpzzT ensemble is more often applied and used with full formal charges for the metal cations of the layers.21–23 On the other hand, the NpT ensemble is often used in conjunction with the QEq. Addition-

10.1021/jp800574d CCC: $40.75  2008 American Chemical Society Published on Web 06/10/2008

Interlayer Species of ZnAlCl Layered Double Hydroxide

J. Phys. Chem. B, Vol. 112, No. 26, 2008 7857

ally, the recent molecular simulations9,10,14–19 of the LDHs have been carried out only in the NpT ensemble. It becomes interesting to understand why the NpzzT ensemble and the QEq have not been used together for the modeling of LDHs so far. In the present work, we aimed to clarify this point by making a comparison between the two different routes. Our present work is different from the previous studies9–19 in a number of respects. First, we compare two different statistical ensembles (NpT and NpzzT) and two different ways for calculating the charges of the atoms of the hydroxide layers. We also pay attention to how the ordered and disordered distribution of the metal cations within the LDH hydroxide layers may affect the structural and dynamical properties of the intercalated species. The results can be compared to X-ray diffraction data. Having established the potential model, the statistical ensemble, and the method for the charge calculation, future studies of intercalated organic anions will be possible. The paper is arranged as follows. Section 2 presents the potential models, the computational procedures, and the statistical ensembles used in this paper. Section 3 establishes the thermal equilibrium of the simulation cell and discusses the value of the layer spacing as a function of the amount of interlayer water molecules for different methodologies. We complete this section with the calculation of the diffusion coefficient of the water molecules in the interlayer region, and we present a preliminary attempt for the calculation of the chemical potential of the intercalated water molecules. Section 4 contains our conclusions and presents future work. 2. Methodology of Simulation 2.1. Potential Models. The LDH layers, chloride anions, and water molecules were modeled using the all-atom (AA) version of the Dreiding force field.26,27 The general potential function is of the form

U)

∑ 21 kb(r-r°)2

(1)

∑

(2)

TABLE 1: Parameters Used with the Dreiding Force Field bond

r° (Å)

kb (kJ mol-1Å-2)

Al-O Zn-O O-H

2.05 2.05 0.98

2930 2930 2930

angle

θ° (deg)

kθ (kJ mol-1 rad-2)

Zn-O-H Al-O-H O-Zn-O(1) O-Zn-O(2) O-Zn-O(3) O-Al-O(1) O-Al-O(2) O-Al-O(3)

120 120 97.1 82.9 180 97.1 82.9 180

419 419 419 419 419 19 19 419

atoms

Rii° (Å)

Dii° (kJ mol-1)

Al Zn O H Cl

4.39 4.54 3.41 3.19 3.95

1.30 0.23 0.40 0.064 1.186

represents the position of the point charge qi. The prime just after the summation term indicates that the sum is performed over all periodic images, n ) (nxLx, nyLy, nzLz) with nx, ny, and nz integers, and over atoms j, except j ) i if n ) 0. The electrostatics interactions were calculated using the Ewald sum method28,29 with the different contributions given by the formula

1 Uelec ) 8πε0 + -

bonds

+

1 kθ(cos θ-cos θ°)2 2 angles

∑

+

1 kφ{1 - cos[l(φ - φ°)]} 2 dihedrals

+

∑∑

N-1

N

i)1 j)i+1

(3)

{ [( ) ( ) ] Dij°

Rij° rij

12

- 2.0

Rij° rij

6

∞

+

∑ ′ |r l

qiqj ij + nL|

} (4)

where kb, kθ, and kφ are the force constants for deformation of bonds, angles, and dihedrals, respectively. The equilibrium values of bond distances and valence angles correspond to r° and θ°, respectively. These parameters are given in Table 1 for the atoms of the hydroxide layers. The expression in eq 2 is replaced by kθ(1 + cos θ) for linear equilibrium geometries. In the dihedral angle term, l is the periodicity and δ is the phase factor. For dihedral angles involving one or two metal atoms, this dihedral contribution is zero. It results that there is no torsional energy contribution in the different hydroxide layers. The inter- and intramolecular interactions consisted of a van der Waals repulsion-dispersion term calculated using the first term of eq 4. The parameters for the interactions between unlike atoms were calculated by using the Lorentz-Berthelot mixing rules (quadratic and arithmetic rules for Dij° and Rij° parameters, respectively). The last term in eq 4 corresponds to the electrostatic energy of the system where rij ) ri - rj and ri

-

(

∞

∑ ∑ ∑ ∑ ∑ ′ |r i

a

j*i

b

∑ 12 exp -

1 8πε0

∑∑∑

k)0 k

i

R 4π3⁄2ε0

a b)a

k2 4R2

i

|

qia exp(ikria)

a

2

qiaqib erf(Rriaib) riaib

∑ ∑ qia2 i

|n|)0

( ) |∑ ∑

1 2Vε0

)

qiaqjb erfc(R|riajb + n|) iajb + n|

(5)

a

where n is the lattice vector of a periodic array of MD cell images. a and b denote atoms of molecules i and j, respectively. In the first term of eq 5, when i ) j, the summation discounts any excluded atoms b of atom a. Let us remind one that the excluded atoms b of atom a are atoms that are linked through a bond, angle, or torsion to atom a. The third term in eq 5 indicates that the summations run over only the excluded atoms b of atom a in the molecule i. riajb is the distance between the atoms a and b belonging to the two different molecules i and j. qia and qib represent the charges on atoms a and b, respectively. 2.2. Computational Procedures. The hydroxide layers of the simulation cell were constructed using the atomic positions from the crystal structure provided by Rietveld refinement30 of the powder X-ray diffraction pattern data (Figure 1) of a Zn 2AlCl LDH having the composition Zn0.677Al0.323(OH)2Cl0.297(CO3)0.013 · 0.591H2O.31 The parameters are given in Table 2 and Figure 1. The simulation supercell consisted of 10 × 9 × 1 rhombohedral (R3jm) unit cells, containing three metal hydroxide layers. When the Zn/Al ratio is 2, each hydroxide layer contains 60 Zn2+ and 30 Al3+. The composition of each hydroxide layer was 68 Zn2+ and 22 Al3+ in the case of a modeling of ratio Zn/Al close to 3. The compositions of the ZnAl LDH supercells were Zn180Al90(OH)540 and Zn204Al66(OH)540 for Zn/Al ratios of 2

7858 J. Phys. Chem. B, Vol. 112, No. 26, 2008

Pisson et al.

TABLE 2: Rietveld Refinement Results Zn2Al(OH)6Cl · 2H2O R3jm

space group lattice parameters a (Å) c (Å) Z temperature (K) λ(Cu KR1) (Å) Rwp (%) RBragg (%) χ2

3.07558(6) 23.248(1) 3 293 1.540 598 9.9/14.6 (BC) 6.3 1.7

a

atom

site

x

y

z

Al Zn OH Cl Ow

3a 3a 6c 18g 36i

0 0 0 0.189(4) 0.29(4)

0 0 0 0 0.503(8)

0 0 0.3760(1) 1 /2 1 /2

a

Biso (Å2) 1.6 1.6 1.6 5.0 5.0

BC ) background corrected.

Figure 2. Snapshots of the simulation supercell of ZnAlCl LDH, consisting of 10 × 9 × 1 rhombohedral (R3jm) unit cells. There are three hydroxide layers and three interlayers. Blue octahedra are Al atoms, brown octahedra are Zn atoms, and the white sticks are OH groups of the hydroxide layer with the O atom represented in red. In the interlayer, green balls are Cl - and water molecules are represented by gray sticks. The structure a represents a Zn2AlCl equilibrium configuration with a complete ordering of Al and 60 water molecules per interlayer, whereas the structure b shows a Zn3AlCl equilibrium configuration with 22 water molecules per interlayer.

Figure 1. Experimental X-ray diffraction pattern (red cross), calculated pattern (solid black line), Bragg reflections (green ticks), and difference profiles (solid blue line) for Zn2Al(OH)6Cl · 2H2O.

and 3, respectively. A total of 90 chloride anions (30 per each interlayer) for Zn/Al ) 2 or 66 Cl- (22 per each interlayer) for Zn/Al ) 3 and a variable number of water molecules are added to complete the simulation cell. Snapshots of the simulation supercells of Zn2AlCl and Zn3AlCl are shown in Figure 2. The water molecules were modeled using the TIP3P model.32 For water, atomic charges (+0.417 e for H and -0.834 e for O) were taken for the TIP3P model and a partial charge of -e was allocated to the chloride anion. These charges are constant along the different simulations. For a Zn/Al ratio equal to 2, the use of the periodic boundary conditions did not allow one to completely satisfy both the disordering of the trivalent cations and the cation avoidance rule.33 Accordingly, we chose to simulate two pseudosystems to study the influence of the trivalent cation avoidance rule. In the system abbreviated as (Zn2AlCl)ord, the cation avoidance rule is satisfied with all of the Al-Al distances greater than the intermetallic distance. Nevertheless, this simulation cell shows an ordered distribution of metal cations with six nearest-neighbor Zn atoms around Al atoms, while each Zn atom is surrounded by three nearest-neighbor Zn atoms and three nearest-neighbor Al atoms. In the system abbreviated as (Zn2AlCl)dis, the distribution of Al atoms is not ordered and some Al-Al distances may not respect the distance criterion; i.e., some Al3+ may occupy adjacent octahedra. When the Zn/Al ratio is equal

to 3, the smaller number of Al atoms allowed one to respect the cation avoidance rule and the disordered distribution of Al in the simulation box. The periodic boundary conditions were applied in the three dimensions. The long-range electrostatic interactions were calculated by the Ewald summation technique. The parameters of this method were R ) 0.3208 Å-1 (convergence parameter) within a relative error of 10-6 and kmax ) {9 × 9 × 8} (the reciprocal space vectors). The equations of motions were integrated using the Verlet Leapfrog algorithm with 1 fs as the time step. The cutoff radius was 10 Å, and the Verlet list sphere radius was fixed to 12 Å. A typical run consisted of 50 ps of equilibration followed by a production phase of an additional 70 ps. The structural and thermodynamic properties were calculated over 7000 configurations saved during the acquisition phase. The configurations were generated using the parallel version of the DL_POLY_MD package34 by using up to 16 processors at a time. 2.3. Statistical Ensemble. The MD simulations were performed at constant composition in the isothermal-isobaric ensemble. The temperature was fixed at T ) 298 K. Within this ensemble, two alternatives could be used. One consisted of imposing the uniaxial pressure component pzz perpendicular to the LDH layers. In this case, only the interlayer spacing was allowed to vary and the x and y dimensions of the simulation cell were kept constant. Within this methodology, the MD simulations were performed more strictly at the NpzzT ensemble. In the standard NpT ensemble, the pressure p was applied in the three directions and the three dimensions of the simulation box were allowed to vary. In all of the cases, p and pzz were

Interlayer Species of ZnAlCl Layered Double Hydroxide

J. Phys. Chem. B, Vol. 112, No. 26, 2008 7859

Figure 3. Temperature profiles of different components across the slab along different directions, (a) Txx(x), (b) Txy(x), (c) Tyy(y) and (d) Tzz(z), calculated in the Zn2AlCl with a complete ordered arrangement of Al and 60 water molecules in each interlayer region.

maintained at 1 atm using the Berendsen algorithm35 with coupling constants 0.1 ps (temperature) and 0.5 ps (pressure). 3. Results and Discussion At the start of the simulation, the water molecules were placed randomly in the interlayer region within the 36i position with a partial occupancy of 2/36. We performed four series of MD simulations. Each series consisted of four simulations with a water content varying from 30 to 75 water molecules per each interlayer region. A series is characterized by the statistical ensemble (NpT or NpzzT) and the charge calculation method (QEq or full charges for the metal cations). In a first step, we focused on the thermal equilibrium of the MD configurations. The evolution of the temperature was monitored by calculating an average slab temperature defined in eq 6 along the γ direction

kBTRβ(γ) )

〈

N-1

∑ H(γi) mi(vi)R · (vi)β i)1

N-1

∑ H(γi) i)1

〉

(6)

where kB is Boltzmann’s constant and mi and (vi)R are the molecular mass and R component of the velocity of the centerof-mass of water molecules or chloride anions, respectively. H(γi) is a top-hat function defined such that

H(γi) )

{

∆γ ∆γ < γi < γ + 2 2 0 otherwise 1 for γ -

(7)

where ∆z is the slab thickness and γ represents x, y, or z directions. The analysis of temperature tensor component Txx(x) and Tyy(y) profiles establishes in parts a and b of Figure 3 that

the temperature can be considered as homogeneous at each x and y, respectively. The average values of Txx and Tyy calculated over the slabs are 302 ( 11 and 301 ( 11 K along the x and y directions, respectively. We also checked that the profile of one off-diagonal temperature tensor component Txy(y) is constant and fluctuates around zero as expected for a system at thermal equilibrium. These plots show that the LDHs are well thermalized in the plane parallel to the hydroxide layers. In the direction normal to the layers, the system presents a heterogeneity with chloride anions and water molecules located in neighboring slabs along z. This means that the profile is not continuous along z, and we decided to increase the thickness of the slab to obtain an accurate estimate of Tzz. Nevertheless, the average of Tzz in slabs representing the hydroxide layers and the interlayer region matches very well with the input temperature. We observed that the diagonal elements of the temperature tensor are equivalent, showing an identical equipartition of the kinetic energy in the three directions. Concerning the mechanical equilibrium of the system, we checked that the off-diagonal elements of the pressure tensor are close to zero and that the on-diagonal elements of the pressure tensor are identical and equal to the trace of the pressure tensor as required for a system at mechanical equilibrium. We checked that the way of modeling the hydroxide layer with atoms completely movable allows a good reproduction of the thermodynamic equilibria. This contrasts with the way of treating the atoms of the hydroxide layers as fixed on the rigid lattice; this approach does not allow the exchange of energy and momentum among the atoms of the layers and those of the interlayer space. Figure 4 shows the evolution of the interlayer spacing as a function of time during the production phase for each interlayer region of the simulation cell in the case of the Zn2AlClord. This figure confirms that the configurations are well-equilibrated with a good convergence of the interlayer distance. The time average

7860 J. Phys. Chem. B, Vol. 112, No. 26, 2008

Pisson et al.

Figure 4. Trajectories of the interlayer spacing for each interlayer region over a simulation time of 50 ps.

TABLE 3: Partial Charges of the Atoms of the Hydroxide Layers Calculated from QEq a QEq (e) atoms

Zn2AlClord

Zn2AlCldis

Zn3AlCl

formal ZnAlCl (e)

Al Zn O H

1.54 0.45 -0.56 0.32

1.44 0.50 -0.56 0.32

1.53 0.44 -0.53 0.3

3 2 variable variable

a When the charges on Zn and Al are 2 and 3, respectively, the charge of the O atoms of the hydroxide layers are variable.

value of this distance is 7.75 ( 0.02 Å for each interlayer region, in agreement with the X-ray diffraction data. The fact that the trajectories are identical for the three layers establishes that the way of modeling this LDH with a periodically repeated simulation cell is valid. We now compare the simulations carried out in the NpT and NpzzT ensembles on Zn2AlClord, Zn2AlCldis, and Zn3AlCl with partial charges on the atoms of the hydroxide layers calculated using the charge equilibration approach (QEq20). The atomic charges are listed in Table 3. This approach allows the charges of the hydroxide layers to respond to changes in their environment and explains why the calculated charges are slightly different between the ordered and disordered arrangements of metal cations in LDHs. We also give for comparison the atomic charges in Zn3AlCl. We checked that the partial charges on the hydroxide layers are consistent with those calculated on the atoms of Mg2AlCl9,17,18 LDHs. Figure 5 presents the evolution of the interlayer distance as a function of the interlayer water content for the NpT and NpzzT simulations with hydroxide layers atom charges calculated from the QEq method. We report for comparison the experimental interlayer distance obtained from X-ray diffraction analysis in the case of ZnAlCl with Zn/Al ratios of 2 and 3.36 Interestingly, we observe that the NpzzT simulations cannot reproduce the experimental interlayer spacing with QEq charges. This contrasts with the NpT simulations that show a good agreement between the computed and experimental interlayer distance. We also observe that the number of water molecules in the interlayer region required for a good prediction of the interlayer distance ranges from 45 to 60. This hydration state is consistent with the chemical formula, leading to 60 interlayer water molecules if we assume 90 cells. In the case of ZnAlCl with a Zn/Al ratio of 3, we observe in Figure 5 that the experimental layer distance is not reached by the simulation. This number of water molecules should be close to 80. When increasing Zn2+/Al3+, a higher water content is expected as a result of the decrease of the layer charge density and of the

Figure 5. Interlayer distances for the ZnAlCl layer hydroxide with molecular ratios of Zn/Al of 2 and 3 as a function of the interlayer water content. The simulations are performed in the NpT and NpzzT ensembles, with the atomic charges of the layers calculated from QEq. The charges on the metal cations are given for completeness in Table 3.

amount of interlayer anions, leaving more space to accommodate additional water molecules. A significant difference between the NpT and NpzzT simulations is the behavior of the interlayer spacing as a function of the hydration state. In other words, the NpzzT simulations show a marked evolution of the interlayer spacing with the water content, whereas the interlayer spacing appears to be less dependent on the number of water molecules in the NpT simulations. This may be explained by the fact that imposing the pressure in the NpT simulation results in changes of the dimension of the simulation box in the three dimensions, whereas only the z dimension is allowed to vary in the NpzzT simulation. Yet, the amplitude of variation is larger in the case of the NpzzT ensemble, whereas it is weakened in the three directions for the NpT ensemble. The NpT simulations allow one to reproduce the evolution of the interlayer space with respect to the water content for the ZnAlCl LDHs with only the use of a single set of charges calculated from the QEq approach. When the atomic charges on the atoms of the hydroxide layers are calculated from QEq, the simulations in the NpzzT overestimate significantly the interlayer spacing whatever the interlayer water content. To further investigate the NpzzT simulations, we applied full formal charges on metal atoms of the hydroxide layers and changed the charge of the O atom of the hydroxide groups in such a way as to predict interlayer spacing. Interestingly, we observe in Figure 6 that it is not able to reproduce the expected interlayer separation distances for a range of charges between -2.0 and -1.0 e. Additionally, parts a and b of Figure 6 suggest that the NpzzT simulations can give interlayer spacings in line with the experiments for Zn/Al ratios of 2 (Figure 6a) and 3 (Figure 6b). Depending on the hydration state and on the Zn/ Al ratio, we conclude that the optimum charge for the O atom of the hydroxide groups is between -1.4 and -1.6 e. This value of the oxygen atomic charge is consistent with those calculated in montmorillonite (-1.42 e24,25,37 and -1.72 e22). We also observe that only the range of charges between -2.0 and -1.6 e allows one to reproduce the increase of the interlayer space as a function of the water content. We demonstrated here that the methodology of simulation is very sensitive to the method of calculating the effective charges on the atoms of the layers.

Interlayer Species of ZnAlCl Layered Double Hydroxide

Figure 6. Interlayer spacing as a function of the charge of the O atoms of the hydroxide layers for different water amounts in the interlayer region calculated in the NpT and NpzzT ensembles: (a) Zn2AlCl with an ordered arrangement of metal cations with a Zn/Al ratio equal to 2; (b) Zn2AlCl with a Zn/Al ratio of 3.

It is then possible to make the interlayer distance calculated from NpzzT simulations in line with the experimental data as long as full formal charges are assigned to the metal atoms of the layers. It also explains why, in the field of molecular simulation of clay minerals, NpT is often associated with the QEq method9,10,14–19 and that the use of full formal charges on the metal cations of layers is done with the NpzzT statistical ensemble.21–25 Nevertheless, the two methodologies lead to similar values for both the interlayer spacing and the interlayer water content, in total agreement with the experimental data, i.e., X-ray diffraction analysis and chemical analysis. Figure 7a shows the total atomic density profiles along the direction normal to the hydroxide surfaces for the Zn2AlClord with 60 interlayer water molecules. We report for comparison the experimental distances d003 and d006 determined from the positions of the 003 and 006 X-ray diffraction reflections. The d003 interlayer distance is 7.76 Å and corresponds to twice the value of d006. The d006 spacing is then the distance between the hydroxide layer and the center of the interlayer spacing. For comparison, the one-dimensional electron density map F(z) of the Zn2AlCl structure along the c stacking axis was constructed. F(z) values were obtained by a Fourier transform of the 00l X-ray diffraction reflection intensities.39 We observe that the simulation supercell conserves its bidimensionality during the course of the simulation and that the simulation is capable of reproducing the lamellar features of the LDH. We also confirm that the different distances d003 and d006 agree very well with those deduced from MD simulations. We checked that the total atomic density profile was symmetrical, leading to three identical local compositions of hydroxide layers and of interlayer spacings. The highest peaks correspond to the atomic densities of the metal atoms; the two peaks close to the

J. Phys. Chem. B, Vol. 112, No. 26, 2008 7861

Figure 7. (a) Total density profiles along the z direction of the supercell simulation (black solid line) and electronic density profile (red solid line). The distances corresponding to the 003 and 006 reflections are given for comparison. (b) Atomic density profiles of the intercalated species as indicated in the legend.

highest peaks represent the O atoms of the hydroxy groups, whereas the H atoms are represented by the smallest peaks. The peak in the midplane of the interlayer corresponds to the local density of water molecules and Cl atoms. We observe that the experimental electronic density profile matches that of the MD simulation. Nevertheless, the electronic density profile is not able to distinguish between the contributions of the hydroxy groups of the layers because fo the poor crystallinity of the Zn2AlCl sample. Part b of Figure 7 zooms on the local composition of the interlayer spacing as a function of z. As expected, the cloride anions and water molecules are located in the center of the interlayer region. Besides, the distance between the layer of chloride anions and the center of the hydroxide layers agrees with that determined from the 006 reflection. Each chloride anion interacts preferentially with 3 × 2 hydroxyl groups at a distance of 2.94 Å of adjacent layers, forming a prismatic environment, and with two water molecules at 2.93 Å. The position of the peak of the atomic density profile of the H atoms of water molecules is located at the same z as that of the density profile of the water O atom. This indicates that the water molecules are parallel to the hydroxide surfaces. The fact that the peak of density of the O atom of the water molecules is wider than that of the chloride anions may indicate that the mobility of the water molecules is higher than that of the chloride anions. To investigate the structure of the interlayer region, we calculated in Figure 8 the contour map of normalized atomic probability densities in the plane of the interlayer (x, y). On the same figure are given the x and y coordinates of the 36i crystallographic sites in the R3jm space group. One can see that

7862 J. Phys. Chem. B, Vol. 112, No. 26, 2008

Pisson et al.

Figure 8. Contour maps of normalized atomic probability densities in the plane of the interlayer (x, y) of Zn2AlClord with 60 intercalated water molecules. Red contours are for chloride atoms, and yellow contours are for O atoms of water molecules. The small white points correspond to the positions of the 36i crystallographic sites in the R3jm space group.

the positions of the maxima of the two-dimensional density profiles of the water molecules coincide with the statistically equivalent positions of the crystallographic sites. This result is consistent with, first, the structural refinement based on XRD data and, second, the MD simulations of Zn-containing LDH performed by Pappalardo et al.16 The water molecules are located on the 36i positions with a partial occupancy of 2/36, and chloride anions are located on the 18g positions with an occupancy of 1/18. These partial occupancies on high-multiplicity positions are related to the structural disorder of the interlayer species. The dynamic properties of the water molecules of the interlayer space were also studied. Interestingly, we calculated in Figure 9a the trajectory of a water molecule in a (x, y) plane in a bulklike system and in the interlayer space of Zn2AlCl. This curve is very instructive by clearly showing that the water molecules are confined in a small space. The translational selfdiffusion coefficient can be determined by two methods: from mean-squared displacement of the molecular centers of mass using the Einstein relation (eq 8)

1 Dmsd ) lim tf∞ 6NN0t

〈∑ ∑ N0

N

t0)1 i)1

[rcom,i(t + t0) - rcom,i(t0)]2

〉

(8)

Figure 9. (a) Trajectory of a water molecule in the (x, y) plane of the interlayer space in an isotropic bulklike system and in the Zn2AlCl ord, as indicated in the legend. (b) Mean-squared displacements calculated in each of the three directions for the interlayer water molecules. (c) Velocity autocorrelation functions of the intercalated water molecules calculated in the x, y, and z directions. The velocity autocorrelation function of isotropic bulk water molecules is given for comparison.

or from the integration of the velocity autocorrelation function38

1 Dvacf ) 3NN0

〈

N0

N

∫ ∑ ∑ [Vcom,i(t + t0) Vcom,i(t0)] ∞

0

t0 )1 i)1

〉

dt (9)

where rcom,i and Vcom,i are the position and velocity of the center of mass of water molecule i, respectively. The mean-squared displacement was calculated over a trajectory of 500 ps using N0 ) 50 time origins. We observe very few differences between the diffusion coefficients calculated from the NpT and NpzzT simulations; these differences are within the standard deviations of the diffusion coefficients. Figure 9b presents the meansquared displacement of the water of the interlayer region in the x, y, and z directions. This figure shows that the meansquared displacements in the x and y directions are equivalent, whereas it is slightly increased in the z direction. The diffusion

coefficient of water calculated from the total mean-squared displacement is 5.1 × 10-12, 2.0 × 10-12, and 6.1 × 10-12 m2 s-1 for Zn2AlClord, Zn2AlCldis, and Zn3AlCl, respectively. The diffusion coefficients of the chloride anions are equal to 3.8 × 10-12, 1.8 × 10-12, and 3.6 × 10-12 m2 s-1 for Zn2AlClord, Zn2AlCldis, and Zn3AlCl, respectively. We also checked that the diffusion coefficients calculated from the velocity autocorrelation functions are equivalent to those calculated from eq 5 within the statistical fluctuations. These calculations suggest that the water molecules have a relatively higher mobility than the chloride anions. Furthermore, these diffusion coefficients appear to be 3 orders of magnitude smaller than those calculated in other mineral clays such as montmorillonites24,25 and 10 times smaller than those calculated in MgAl LDH containing organic guest anions.10

Interlayer Species of ZnAlCl Layered Double Hydroxide

J. Phys. Chem. B, Vol. 112, No. 26, 2008 7863

Part c of Figure 9 shows the center of mass velocity autocorrelation functions of water molecules in the interlayer region in the three directions. We report for comparison the velocity correlation function of bulk water molecules. Interestingly, we observe very significant changes in these correlation functions between the bulk and interlayer water molecules. The velocity correlation functions of the water molecules in the interlayer region exhibit negative regions. The velocity autocorrelation functions calculated in the x and y directions are comparable, characterized by a first negative peak, and tend to a zero value with dumped oscillations. In the z direction, the amplitude of the first negative region is much more pronounced, and we observe a clear oscillatory behavior with greater amplitudes. The negative portion of the velocity autocorrelation function can be ascribed to a reversal of the velocities with respect to an initial value updated at each time origin after a time lapse that appears to be smaller in the x and y directions. The more marked negative lobe of the velocity autocorrelation functions in the direction normal to the hydroxide layers indicates that the cage effect is less pronounced in the z direction, leading to a more significant rattling of the water molecules along this direction. We have completed this study by the calculation of the chemical potential of water molecules of the interlayer region. This part constitutes a preliminary attempt to show the feasibility of such a calculation even in the case of water molecules in a highly confined space. The chemical potential in the NpT statistical ensemble can be calculated by µ ) kBT ln F

(

h2 2πmkBT

)µid + µex

)

3⁄2

- kBT ln

[

〈V exp(-∆U ⁄ kBT)〉NpT 〈V 〉

]

(10)

(11)

where V is the volume of the system, T is the absolute temperature, and kB and h are Boltzmann’s and Planck’s constants, respectively. p refers to the pressure of the system, and m is the molecular mass. The most common method for the calculation of the excess contribution of the chemical potential is the Widom test particle40 and consists of inserting a ghost particle randomly into a simulation box and calculating the energy of its interaction with the N particles. The operational expression for the calculation of the chemical potential is given in eq 10, where ∆U ) U(rN+1) - U(rN) represents the energy of interaction of the (N + 1)th particle with the other N particles. The angular brackets 〈...〉NpT denote a canonical ensemble average over the configuration space of the N particles, F is the mean number density expressed as the number of particles per cubic meter. The first term in eq 10 is the ideal gas chemical potential, and the second term defines the excess chemical potential of the system. When the simulation is performed in the NpzzT ensemble, the volume V in the second part of eq 10 is replaced by Lz. The difficulty here arises from the insertion of a charged water molecule. In other words, the calculation of the ∆U energy part implies the calculation of the dispersion-repulsion energy contribution as well as of the electrostatic part from the different contributions of the Ewald summation method. The insertion of a water molecule with its partial charges does not allow the convergence of the calculation. To avoid this problem, we have decomposed the calculation into two steps: the first step consists of inserting the uncharged water molecule and calculating the contributiontothechemicalpotentialfromthedispersion-repulsion energy contributions, and the second stage takes advantage of gradually increasing the partial charges on the atoms using a perturbation process with the FEP 41,42 and TI approaches.42–44

Figure 10. Total chemical potential of water molecules as a function of the number of water molecules of the interlayer space in Zn2AlClord.

The ideal contribution calculated from the first term of eq 10 at T ) 298 K, p ) 1 atm, and F ) 1 g cm-3 is -19.1 kJ mol-1. In a first step, we have calculated the excess part of the chemical potential for water bulk molecules. The dispersion-repulsion contribution of the excess chemical potential was calculated using the Lennard-Jones parameters because of the insertion of the O atom of the water TIP3P model. This contribution is then equal to + 9.2 kJ mol-1. The apparition of the partial charges on the H and O atoms of a water molecule was carried out in 11 consecutive windows. The perturbations were performed in both directions (“double-wide sampling”); the difference between the forward (+∆λ) and backward (-∆λ) simulations gives a lower-bound estimate of the error in the calculations. The electrostatic contribution to the chemical potential is then equal to -34.7 ( 0.5 kJ mol-1 for the FEP and TI approaches. The calculated total chemical potential of water at T ) 298 K is -44.6 kJ mol-1, which compares very well with the value of -45.5 kJ mol-1 given by Ben-Naim and Marcus.45 This was extended to the calculation of the chemical potential of the interlayer water molecules of Zn2AlClord. The ordering arrangement of the water molecules requires the chemical potential calculation for various water molecules of the interlayer in order to obtain a consistent average value of µ. Therefore, the calculation of µ becomes very time-consuming, and the average value of µ shows important fluctuations. Figure 10 shows that the calculated chemical potential is smaller than the corresponding bulk water chemical potential. For a total amount of water ranging from 50 to 60, it is then possible to obtain qualitatively within the statistical fluctuations a chemical potential close to the bulk value. The convergence of the chemical potential could be further improved by calculating the chemical potential for each intercalated water molecule. The purpose of the present study was not to do such calculations but to demonstrate that it is feasible even for confined water molecules. This calculation confirms what was reported in the introduction of this paper concerning the difficulty of equilibrating these types of systems in the µVT ensemble. Let us recall that the chemical potential has been already calculated in montmorillonites.46,47 However, in such mineral clays owing to their swelling properties, the interlayer water molecules sample much more of the interlayer space and consequently make the chemical potential calculation easier to equilibrate. 4. Conclusions The structural features and dynamic properties of the interlayer species of the ZnAlCl LDHs have been studied

7864 J. Phys. Chem. B, Vol. 112, No. 26, 2008 from MD simulations. We have shown that the MD configurations are well thermal equilibrated and that the mechanical equilibrium is respected. We have established that the MD method is capable of reproducing both the interlayer spacing and the interlayer water molecule content in the NpT and NpzzT statistical ensembles. The agreement with the X-ray diffractions data and chemical analysis depends on the method used for calculation of the charges of the metal cations of the hydroxide layers. More precisely, reliable values are obtained in the NpT ensemble when the charges of the metal cations are calculated from the equilibration charge method and in the NpzzT ensemble with full formal charges for the metal cations. We have checked that water molecules of the interlayer space sample the 36i crystallographic sites in agreement with the R3jm space group reported for this LDH. Furthermore, the diffusion coefficient of these water molecules is 3 orders of magnitude smaller than that of the intercalated water molecules in montmorillonite. We have also given a possible methodology for the calculation of the chemical potential of the interlayer water molecules. Despite the high confinement of the water molecules, we have shown that this type of calculation remains feasible from a methodological viewpoint. Nevertheless, the calculation of the chemical potential cannot be used to predict the amount of water molecules in the interlayer space of LDHs because of the large error bars relative to a poor sampling of the interlayer space. Given that the NpzzT ensemble used with full formal charges on the metal cations of the hydroxide layers is capable of reproducing the main structural and dynamic features of the intercalated species of ZnAlCl LDHs, we aim at using it in a forthcoming work for the study of the intercalation of dicarboxylate anions. Acknowledgment. The authors acknowledge the Institut du De´veloppement et des Ressources en Informatique Scientifique IDRIS (CNRS) for a generous allocation of CPU time on parallel computers. P.M. is indebted to F. Goujon, A. Ghoufi, and Y. Israeli for helpful discussions. References and Notes (1) Leroux, F.; Forano, C.; Prevot, V.; Taviot-Gue´ho, C. Layered Double Hydroxide Assemblies and Related Nanocomposites in Encyclopedia of Nanosciences and Nanotechnology; American Chemical Society: Washington, DC, 2008;in press. (2) Guliants, V. V.; Carreon, M. A.; Lin, Y. S. J. Membr. Sci. 2004, 235, 53. (3) Forano, C.; Vial, S.; Mousty, C. Curr. Nanosci. 2006, 2, 283. (4) Leroux, F.; Taviot-Gue´ho, C. J. Mater. Chem. 2005, 15, 3628. (5) Oh, J. M.; Choi, S. J.; Kim, S. T.; Choy, J. H. Bioconjugate Chem. 2006, 17, 1411. (6) Del Hoyo, C. Appl. Clay Sci. 2007, 36, 103. (7) Greenwell, H. C.; Holliman, P. J.; Jones, W.; Velasco, B. Catal. Today 2006, 114, 397. (8) Latterini, L.; Nochetti, M.; Aloisi, G. G.; Costantino, U.; Elisei, F. Inorg. Chim. Acta 2007, 360, 128. (9) Aicken, M.; Bell, I. S.; Coveney, P. V.; Jones, W. J. AdV. Mater. 1997, 9, 496.

Pisson et al. (10) Newman, S. P.; Williams, S. J.; Coveney, P. V.; Jones, W. J. Phys. Chem. B 1998, 102, 6710. (11) Fogg, A. M.; Rohl, A. L.; Parkinson, G. M.; O’Hare, D. Chem. Mater. 1999, 11, 1194. (12) Wang, J.; Kalinichev, A. G.; Kirkpatrick, R. J.; Hou, X. Chem. Mater. 2001, 13, 145. (13) Trave, A.; Selloni, A.; Goursot, A.; Tichit, D.; Weber, J. J. Phys. Chem. B 2002, 106, 12291. (14) Newman, S. P.; Cristina, T. D.; Coveney, P. V. Langmuir 2002, 18, 2933. (15) Greenwell, H. C.; Jones, W.; Newman, S. P.; Coveney, P. V. J. Mol. Struct. 2003, 647, 75. (16) Lombardo, G. M.; Pappalardo, G.; Punzo, F.; Costantino, F.; Costantino, U.; Sisani, M. Eur. J. Inorg. Chem. 2005, 5026. (17) Mohanambe, L.; Vasudevan, S. Langmuir 2005, 21, 10735. (18) Mohanambe, L.; Vasudevan, S. J. Phys. Chem. B 2005, 109, 15651. (19) Kumar, P. P.; Kalinichev, A. G.; Kirkpatrick, R. J. J. Phys. Chem. B 2006, 110, 3841. (20) Rappe´, A. K.; Goddard, W. A., III; Skiff, W. M. J. Phys. Chem. B 1991, 95, 3358. (21) Boek, E. S.; Coveney, P. V.; Skipper, N. T. J. Am. Chem. Soc. 1995, 117, 12608. (22) Skipper, N. T.; Chang, F. R. C.; Sposito, G. Clays Clay Miner. 1995, 43, 285. (23) Dufreche, J. F.; Marry, V.; Bernard, O.; Turq, P. Colloids Surf., A 2001, 195, 171. (24) Marry, V.; Turq, P.; Cartailler, T.; Levesque, D. J. Chem. Phys. 2002, 117, 3454. (25) Marry, V.; Turq, P. J. Phys. Chem. B 2003, 107, 1832. (26) Mayo, S. L.; Olafson, D.; Goddard, W. A., III J. Phys. Chem. B 1990, 94, 8897. (27) Rappe´, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A., III; Skiff, W. M. J. Am. Chem. Soc 1992, 114, 10024. (28) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1989. (29) Smith, E. R. Proc. R. Soc. London A 1981, 375, 475. (30) Rodirguez-Carjaval, J. Recent DeVelopments of the Program FULLPROF, in Commission on Powder Diffraction IUCr), in Newsletter 2001, 26, 12. (31) Isra¨eli, Y.; Taviot-Gue´ho, C.; Besse, J. P.; Morel, J. P.; MorelDesrosiers, N. J. Chem. Soc., Dalton Trans. 2000, 791. (32) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D. J. Chem. Phys. 1983, 79, 926. (33) Radoslovich, E. W. Am. Mineral. 1963, 48, 76. (34) DL_POLY is a parallel MD simulation package developed at the Daresbury Laboratory Project for Computer Simulation under the auspices of the EPSRC for the Collaborative Computational Project for Computer Simulation of Condensed Phases (CCP5) and the Advanced Research Computing Group (ARCG) at the Daresbury Laboratory. (35) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, A.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (36) Roussel, H. Thesis, Universite´ Paris-Sud XI, Paris, France, 1999. (37) Smit, D. E. Langmuir 1998, 14, 5959. (38) Quarrie, D. A. Mc. Statistical Mechanics; University Science Book: Herndon, VA, 2000. (39) Wittingham, M. S.; Jacobson, A. J., Intercalation Chemistry; Academic Press: New York, 1982. (40) Widom, B. J. Chem. Phys. 1963, 39, 2802. (41) Zwanzig, R. W. J. Chem. Phys. 1954, 22, 1420. (42) Mezei, M.; Beveridge, D. L. Ann. N.Y. Acad. Sci. U.S.A. 1986, 482, 1. (43) Pearlman, D. A. J. Comput. Chem. 1994, 15, 105. (44) Mezei, M.; Swaminathan, S.; Beveridge, D. L. J. Am. Chem. Soc. 1978, 100, 3255. (45) Ben-Naim, A.; Marcus, Y. J. Chem. Phys. 1984, 81, 2016. (46) Hensen, E. J. M.; Tambach, T. J.; Bliek, A.; Smit, B. J. Chem. Phys. 2001, 115, 3322. (47) Odriozola, G.; Guevara-Rodriguez, J. Langmuir 2004, 20, 2010.

JP800574D