What Is the Structure of Kaolinite? Reconciling Theory and Experiment

Apr 23, 2009 - Institute, Australian Nuclear Science and Technology Organisation, PMB1, ... ReceiVed: NoVember 28, 2008; ReVised Manuscript ReceiVed: ...
0 downloads 0 Views 655KB Size
Download PDF
6756

J. Phys. Chem. B 2009, 113, 6756–6765

What Is the Structure of Kaolinite? Reconciling Theory and Experiment Claire E. White,† John L. Provis,*,† Daniel P. Riley,‡ Gordon J. Kearley,§ and Jannie S. J. van Deventer† Department of Chemical & Biomolecular Engineering, UniVersity of Melbourne, Victoria 3010, Australia; Department of Mechanical Engineering, UniVersity of Melbourne, Victoria 3010, Australia; and Bragg Institute, Australian Nuclear Science and Technology Organisation, PMB1, Menai, NSW 2234, Australia ReceiVed: NoVember 28, 2008; ReVised Manuscript ReceiVed: March 19, 2009

Density functional modeling of the crystalline layered aluminosilicate mineral kaolinite is conducted, first to reconcile discrepancies in the literature regarding the exact geometry of the inner and inner surface hydroxyl groups, and second to investigate the performance of selected exchange-correlation functionals in providing accurate structural information. A detailed evaluation of published experimental and computational structures is given, highlighting disagreements in space groups, hydroxyl bond lengths, and bond angles. A major aim of this paper is to resolve these discrepancies through computations. Computed structures are compared via total energy calculations and validated against experimental structures by comparing computed neutron diffractograms, and a final assessment is performed using vibrational spectra from inelastic neutron scattering. The density functional modeling is carried out at a sufficiently high level of theory to provide accurate structure predictions while keeping computational requirements low enough to enable the use of the structures in largescale calculations. It is found that the best functional to use for efficient density functional modeling of kaolinite using the DMol3 software package is the BLYP functional. The computed structure for kaolinite at 0 K has C1 symmetry, with the inner hydroxyl group angled slightly above the a,b plane and the inner surface hydroxyls aligned close to perpendicular to that plane. Introduction Kaolinite is an important mineral with numerous applications across many industry sectors. Its uses include the production of ceramics, cosmetics, paints, and paper, as well as being used as a sorbent for pollutants and comprising a major part of many soils worldwide.1,2 Metakaolin (dehydroxylated kaolinite) is also used as a supplementary cementitious material in concretes3 as well as in the production of geopolymer.4 The total annual industrial usage of kaolinite clays amounts to approximately 40 × 106 tonnes.2 Here, kaolinite is used as a case study to show how experimental and computational techniques can be employed simultaneously to obtain the total structure solution of a crystalline material. Kaolinite [Al4Si4O10(OH)8] is a 1:1 dioctahedral phyllosilicate (layered aluminosilicate). Each layer contains a sheet of SiO4 tetrahedra forming six-membered silicate rings connected to the sheet of AlO6 octahedra through the apical oxygens.5 Two of the hydrogens are bonded to the nonapical upper oxygens in the AlO6 octahedra. The other six hydrogens are bonded to the lower oxygens in the AlO6 octahedra, also forming hydrogen bonds with the basal oxygens in the silicate sheet of the next layer. The structure of kaolinite, depicting the tetrahedral and octahedral sheets, is given in Figure 1. Experimental investigations into the structure of kaolinite have reported conflicting results, including space group and hydroxyl bond lengths and orientations.6-12 Reasons for such discrepan* To whom correspondence should be addressed. Tel: +61 3 8344 8755. Fax: +61 3 8344 4153. Email: [email protected]. † Department of Chemical & Biomolecular Engineering, University of Melbourne. ‡ Department of Mechanical Engineering, University of Melbourne. § Australian Nuclear Science and Technology Organisation.

Figure 1. 34-atom unit cell of kaolinite showing the silicate tetrahedral sheet, alumina octahedral sheet, hydroxyl positions, and atom labels used throughout this paper.

cies have been discussed in those investigations and may include impurities in samples, temperature, preferential orientation of the crystallites, and difficulties in structure refinement. Over the past 20-30 years, numerous studies of the structure of kaolinite have been performed using ab initio computational methods.5,13-20 As is the case in the various published experimental investigations, these theoretically based studies have also reported similarly conflicting results. This has become particularly noticeable due to improvements in the accuracy of ab initio methods as computational power increased. Therefore, to investigate the impact of different functionals on the ability to conduct efficient yet accurate computations of the kaolinite structure, a range of functionals have been assessed in this paper. This type of investigation has become possible due to improvement in computing facilities and algorithms in recent years,

10.1021/jp810448t CCC: $40.75  2009 American Chemical Society Published on Web 04/23/2009

Structure of Kaolinite

J. Phys. Chem. B, Vol. 113, No. 19, 2009 6757

TABLE 1: Different Functionals Used in the Analysis of Kaolinitea combination

qualityb

functional

basis set

rel computation time (%)c

1 2 3 4 5 6 7 8

fine fine fine fine fine fine fine fine

DNP DNP DNP DNP DNP DNP DNP DNP

95 55 49 85 88 100 84 88

9 10

fine fine

GGA-BLYP LDA-PWC27 LDA-VWN28 GGA-PW9127 GGA-BP27,29 GGA-PBE30 GGA-BOP29,31 GGA-VWN-BP (GGA-BP, combination 5, with the local correlation replaced by VWN functional) GGA-RPBE32 GGA-HCTH33

DNP DNP

93 89

a Two classes of functionals have been investigated: generalized gradient-corrected (GGA) functionals and local spin-density approximation (LDA) functionals. b Quality of computation controls the following parameters: SCF convergence tolerance, resolution of the integration grid, how the k-point grid is derived, the combined geometry optimization convergence thresholds for energy change, maximum force, and maximum displacement. Fine quality uses 10-5 hartree convergence tolerance, 0.002 hartree/Å max. force, 0.005 Å max. displacement, 5 × 3 × 4 k-points, and 4.8 Å orbital cutoff. c Computational time expressed as percentage compared with computation that required the longest computational time (Combination 6, 806.30 min using a quad-core 2.0 GHz Intel Core 2 processor).

which has resulted in the ability to calculate the differences introduced into computations by the various approximations inherent in density functional modeling codes. The experimental techniques used to characterize kaolinite have provided useful results which may be used to validate computations. For instance, vibrational spectroscopy21-23 has provided very accurate bond angle and bond length data. However, this technique does not explicitly give atom positions, which are more commonly obtained from scattering methods. With respect to kaolinite and the determination of hydrogen atoms, X-ray scattering has limitations and neutron scattering has its difficulties, as will be discussed below. Hence, due to the various different experimental and computed structures available in the literature, density functional modeling has been utilized here to provide conclusive evidence for the locations of these hydrogen atoms, and to reconcile uncertainties in the structure of kaolinite, particularly with regard to the positions and orientations of the hydroxyl groups. This has been achieved by using selected experimental results presented in the literature, in conjunction with new computational and experimental data. Reasons for discrepancies between experiment and theory have been made clear, including tendencies for theory to overestimate hydrogen bond lengths, and the limitations inherent in experimental data acquisition and analysis. Methods Computational Procedures. The density functional modeling software DMol3 version 4.2 (Accelrys, San Diego, CA) was used to perform geometry optimizations starting from the experimental P1 crystal structure of kaolinite given by Young and Hewat (cell parameters a ) 5.1490 Å, b ) 8.9335 Å, c ) 7.3844 Å, R ) 91.930°, β ) 105.042°, γ ) 89.791°).12 The Young and Hewat kaolinite atomic structure (space group P1) was chosen as the starting point for density functional calculations to avoid potentially introducing spurious symmetry from a less-primitive unit cell (Bish9 conducted a more detailed refinement, but used space group C1). Investigations were performed to assess the performance of the different exchange-correlation functionals. From these investigations the most accurate computed structure of kaolinite was determined. Cell parameters were held constant at the values determined by Young and Hewat.12 It is well-

known that DFT had difficulty in reproducing long-range dispersive interactions like the long-range tails of van der Waals interactions (see, e.g., Grimme24 for a discussion). As a result, the unit cell will change size unrealistically (normally expands) and it is best to constrain the unit cell to the experimental values that are usually known accurately. This leads to a pressure in the cell (∼10 kbar) which conveniently accounts for the dispersion interactions. This has been tested directly for the calculation of weak rotational potentials for methyl groups. These depend significantly on van der Waals interactions and can be obtained with a precision of ∼90%; see, for example, Johnson et al.25 In strongly hydrogenbonded systems this is less important since the H-bonds tend to prevent the cell expansion. An initial investigation determined the influence of the various functionals. The DNP basis set (double numerical plus polarization),26 which consists of two atomic orbitals for each occupied atomic orbital plus a polarization p-function on all hydrogen atoms, was used in all calculations. Ten different functionals were used in this assessment (Table 1). The primary methodology used here to evaluate functional selection involves comparing computed neutron diffractograms. In order to compare neutron diffractograms, an experimentally determined kaolinite unit cell was required. The atom positions within the unit cell were taken from Bish,9 since this Rietveld refined structure of kaolinite was obtained from high-resolution neutron data collected at 1.5 K (reducing thermal effects compared to room-temperature data) and is based on careful refinement. The software package RIETICA34 was used to simulate neutron diffractograms from both experimental and DFT unit cells. The neutron wavelength employed was 1.9102 Å (as used by Young and Hewat12 and Bish9). All neutron diffractograms were produced using unit cell parameters determined by Young and Hewat12 to ensure consistency between the experimental and computed kaolinite structures. Leastsquares error analysis (sum of squared differences, χ2, Rp, and Rwp; eqs 1-4) was used to quantify the differences between the neutron diffractograms of computed and experimental structures and hence identify the most accurate computed structure.

6758

J. Phys. Chem. B, Vol. 113, No. 19, 2009

White et al.

N

sum of squared differences )

∑ (yi,comp - yi,Bish)2

(1)

i)1

N

χ2 )

∑ wi(yi,comp - yi,Bish)2 i)1

(2)

N N

Rp )

Rwp )

[

∑ |yi,comp - yi,Bish| i)1

∑ yi,comp i)1

]

1/2

N

∑

(3)

N

wi(yi,comp - yi,Bish)2

i)1

N

2 ∑ wiyi,comp i)1

diffractograms were collected at 20 K intervals from 300 to 80 K, using a Cryostream blower. Each data set was collected from 5 to 85° 2θ as the sum of two 15 min scans, with the detector position offset by 5° between scans to account for gaps in the coverage of the strip detector. Data were corrected for the capillary background contribution, and lattice parameters were obtained by Rietveld refinement using the software package RIETICA34 and a NIST LaB6 standard for wavelength calibration. The results were then linearly extrapolated to 15 K to obtain the lattice parameters used in the theoretical calculations of INS vibrational spectra (see Supplementary Data 1 in the Supporting Information for extrapolation results). Results and Discussion

(4)

Here, yi,comp is the intensity of the computational neutron pattern at 2θ ) i, yi,Bish is the intensity of the computed neutron diffractogram of the refined structure of Bish9 at 2θ ) i, wi is the weighting factor (wi ) 1/yi,comp), and N is the total number of yi,comp’s. Once the most accurate computed structure was determined, an additional check of accuracy was performed by comparing the inelastic neutron scattering (INS) vibrational spectra from experiment and theory. Johnston et al.35 performed INS on kaolinite at 15 K, providing the experimental vibrational spectrum. This spectrum was compared with the theoretical spectrum obtained from the most accurate computed structure calculated using DMol3, with the comparison providing information pertaining to the accuracy of the computed structure. For the vibrational calculations the structures (using DMol3 and VASP) were first minimized with the cell itself being constrained to the experimental values (see Experimental Procedures for determination of these cell values). A series of single-point energy calculations were made for a series of structures in which the crystallographically distinct atoms were displaced by 0.03 Å in positive and negative directions along the x, y, and z directions. These calculations gave the Hellmann-Feynman (HF) forces from which DMol3 or VASP36-38 calculate the γ-point vibrational frequencies and mean-square displacements within the harmonic approximation. A modified version of CLIMAX39 was used to calculate the INS spectrum from the vibrational frequencies and atomic displacements, taking account of the experimental resolution function, multiphonon modes, and the effects of phonon wings. The molecular dynamics simulation was made using VASP at low precision for a single unit cell with a single k-point. In order to populate the vibrational modes a target temperature of 150 K was chosen and an initial thermalization using fixed NVT was made for 8200 fs with a time step of 1 fs. An NVE run was then made using the final trajectories of the NVT simulation as starting values. The NVE trajectories were taken for calculation of the INS spectrum that was made using NMoldyn.40 Experimental Procedures. In order to obtain lattice parameters at 15 K to enable theoretical calculations of the INS vibrational spectra, a synchrotron X-ray powder diffraction study of a standard kaolinite sample (KGa-1b, Source Clays Repository, Columbia, MO) was conducted at the Powder Diffraction beamline at the Australian Synchrotron.41 The sample was loaded into a 0.5 mm glass capillary on a spinner head, and

Evaluation of Published Experimental Structures. Experimental techniques used to determine the structure of kaolinite have been successful in resolving the positions of silicon, aluminum, and oxygen atoms; however, determination of the exact positions of hydrogen is limited by the technique used. X-ray diffraction techniques, both laboratory and synchrotron based, are unable to provide an accurate hydrogen position information due to the weak interaction of X-rays with hydrogen atoms. Neutron diffraction techniques are more successful due to the fact that (i) hydrogen atoms interact more strongly with neutrons than with X-rays and (ii) the hydrogen atoms can be isostructurally exchanged for deuterium which, due to its greater coherent scattering length, results in a higher useful diffraction signal from these atoms. However, even though theoretically the position of the deuterium atoms can be accurately defined, in reality this has not proven to be the case. As the lightest element in kaolinite, deuterium (and the OD bond) is strongly affected by temperature, which contributes to discrepancies in experimentally refined kaolinite structures. However, even structure refinements from data collected at the same temperature report conflicting results. There have been numerous experimental investigations into the structure of kaolinite, reporting varying results for the positions of hydrogen atoms and the bond lengths of the hydroxyls.6-12 The results are displayed in Tables 2 and 3, which show a range of reported orientations and bond lengths for both the inner hydroxyl and the inner surface hydroxyl groups (see Figure 1 for hydroxyl classification). There is general consensus that the inner hydroxyls lie close to parallel to the a,b plane, whereas the inner surface hydroxyls are almost perpendicular to the a,b plane. However, the exact orientations of the hydroxyls are not consistent between the various investigations (Table 2). This is also true for the hydroxyl bond lengths reported in the literature for kaolinite (Table 3). Reasons for the discrepancies in determination of the unit cell structure of kaolinite may also include different sample environments, sources of samples (giving differences in crystallinity, purity and stacking faults) and methods of structure refinement. Another discrepancy in the literature on the kaolinite structure regards the space group of the triclinic unit cell, reported as being either P1 or C1. Giese and Datta6 used C1 symmetry and minimization of the electrostatic energy to find the positions of the hydrogens. Suitch and Young11 carried out Rietveld refinement43,44 of X-ray and neutron data on a P1 unit cell of kaolinite, using two starting models (Adams and Hewat,42 denoted AHR, and Giese and Datta,6 GDR). The sample used for neutron diffraction was not deuterated. These results showed that the space group was P1 due to the positions of the inner hydrogens.

Structure of Kaolinite

J. Phys. Chem. B, Vol. 113, No. 19, 2009 6759

TABLE 2: Experimentally Determined Hydroxyl Group Orientations in Kaolinitea experimental structure Switch and Young (1983) atom pairs (Figure 1)

Giese and Datta (1973)b

AHRc

GDRc

Young and Hewat (1988)

temperature space group O(11)-H(1) O(12)-H(2) O(13)-H(3) O(14)-H(4) O(15)-H(5) O(16)-H(6) O(17)-H(7) O(18)-H(8)

C1 14 77 14 77 14 77 14 77

ambient P1 -23 56 61 54 10 65 52 85

ambient P1 -31 61 55 77 11 51 53 62

1.5 K P1 12 63 50 66 -22 75 87 50

Bish (1993) 1.5 K C1 0 73 68 60 0 73 68 60

Akiba et al. (1997)

Neder et al. (1999)

ambient C1 -21 68 61 62 -21 68 61 62

ambient C1 12 64 73 47 12 64 73 47

a All angles are given in degrees relative to the a,b plane. b Giese and Datta6 hydroxyl group orientation results were determined by minimizing the electrostatic energy. These results were not experimentally determined, but nor are they from ab initio computations, hence they are discussed with the experimental results. c Starting models AHR and GDR refer to structures determined by Adams and Hewat42 and Giese and Datta.6

TABLE 3: Experimentally Determined Hydroxyl Bond Lengths (in Å) in Kaolinite experimental structure Suitch and Young (1983)

atom pairs (Figure 1)

Giese and Datta (1973)

AHRa

GDRa

Young and Hewat (1988)

Bish (1993)

Akiba et al. (1997)

Neder et al. (1999)

temperature O(11)-H(1) O(12)-H(2) O(13)-H(3) O(14)-H(4) O(15)-H(5) O(16)-H(6) O(17)-H(7) O(18)-H(8)

b -

ambient 0.85 1.48 1.11 0.93 1.15 0.86 0.83 1.26

ambient 0.88 1.21 1.1 1.17 1.27 1.52 1.02 1.17

1.5 K 1.01 1.02 1.09 0.91 1.01 0.96 1.05 1.05

1.5 K 0.975 0.982 0.976 0.975 0.975 0.982 0.976 0.975

ambient 0.93 0.94 0.90 0.83 0.93 0.94 0.90 0.83

ambient 0.75 0.76 0.77 0.88 0.75 0.76 0.77 0.88

a Starting models AHR and GDR refer to structures determined by Adams and Hewat42 and Giese and Datta.6 b Oxygen positions not given by Giese and Datta,6 hence no bond lengths can be compared.

Young and Hewat12 carried out refinement of neutron diffraction data collected from a nondeuterated kaolinite sample, concluding that the unit cell is P1. Bish9 claimed that if kaolinite was P1 then there would be extra peaks in the X-ray and neutron diffraction patterns, which are not visible in the data of Young and Hewat,12 and hence the unit cell is C1. Bish9 also reported that some of the bond lengths in the Young and Hewat12 structure were unreasonable, falling outside accepted ranges. Hence, Bish9 refined the C1 structure from the Young and Hewat12 data set using soft constraints, with the resulting structure of the non-hydrogen atoms being of space group C1. However, as noted by Bish,9 the hydroxyl bond distances were undoubtedly affected by the use of soft distance constraints. Akiba et al.7 used a synthetic sample of deuterated kaolinite to obtain a time-of-flight neutron powder diffraction pattern. The pattern was refined using the Bish and von Dreele8 model of kaolinite (space group C1) as the starting structure for the non-hydrogen atoms and the Young and Hewat12 structure for initial positions of the hydrogen atoms. During the refinement, constraints were used on Al-O and Si-O distances as well as on one of the hydrogen atoms. Their conclusion was that kaolinite was of space group C1. Neder et al.10 obtained synchrotron X-ray diffraction patterns from single crystals of natural kaolinite, and the results indicated C-centering and space group C1. However, Neder et al.10 stated that the refinement was not straightforward, and even with constraints the position of the inner hydrogen atom could not be refined.

Thus, it is clear that the space group of kaolinite has been debated for some time, with most recent investigations concluding a space group of C1. Hence, one of the aims of this investigation is to independently confirm the space group of kaolinite by computational techniques. Evaluation of Published Computational Structures. Numerous computational studies using ab initio methods or density functional theory (DFT) have been performed on the atomic structure of kaolinite.5,13-20 However, to date no study has assessed an extensive range of functionals with respect to the geometry and total energy of the kaolinite unit cell and by comparison with experimental data, with a view toward finding the most efficient method for conducting further computations using the structure. Tables 4 and 5 show the computed bond angles and bond lengths of the hydroxyls from previous ab initio and DFT studies of kaolinite. The results from the various investigations exhibit numerous differences; however, they do follow a trend similar to the experimental work, where the inner hydroxyls are close to parallel to the a,b plane while the inner surface hydroxyls are close to perpendicular to the a,b plane. Such variation in computed structural properties indicates the importance of evaluating the parameters available in ab initio or DFT modeling before drawing conclusions about the structure being investigated. Hu and Michaelides45 investigated various properties of the kaolinite crystal structure and performed DFT calculations using the PBE functional and a plane-wave basis set (CASTEP code).

6760

J. Phys. Chem. B, Vol. 113, No. 19, 2009

White et al.

TABLE 4: Theoretically Determined Hydroxyl Group Orientations in Kaolinitea computational structure Benco et al. (2001)

computational method space group used space group determined inner hydroxyls inner surface hydroxylse

Hess and Saunders (1992)

Hobbs et al. (1997)d

Balan et al. (2001)c

horizontal structureb

vertical structureb

Sato et al. (2004)

Hu and Michaelides (2008)c,d

HF/STO-3G, HF/6-21G C1

CA-PZ/DZP P1 C1 3.76 91.16 78.24 64.58

PBE/PW C1

PW91/PW C1

PW91/PW C1

PBE/PW C1

25.2 5.0 84.5 78.0

0.6 71.8 67.5 75.4

-8.9 70.8 67.8 52.2

PBE/PW P1 C1

3.1 fixed at 57.7

a All angles are given in degrees relative to the a,b plane. b Benco et al.16 reported two possible configurations of kaolinite. The “vertical structure” has all three inner surface hydroxyls pointing toward the adjacent silicate layer. The “horizontal” structure has one of the three inner surface hydroxyls almost parallel to the silicate layer, with the remaining two hydroxyls perpendicular. c Balan et al.15 and Hu and Michaelides45 did not report hydroxyl group orientations. d Hobbs et al.14 and Hu and Michaelides45 were the only studies to report investigating the space group of kaolinite. e Values given for each of the three inner surface hydroxyls in each study should be compared with all other inner surface hydroxyls from other studies as it is not reported which value corresponds to which specific inner surface hydroxyl group.

TABLE 5: Theoretically Determined Hydroxyl Bond Lengths (in Å) in Kaolinite computational structure

computational method space group used space group determined inner hydroxyl bond lengths inner surface hydroxyl bond lengthsa

Hess and Saunders (1992)

Hobbs et al. (1997)

Balan et al. (2001)

Castro and Martins (2005)

Tosoni et al. (2006)

Hu and Michaelides (2008)

HF/STO-3G, HF/6-21G C1

CA-PZ/DZP P1 C1 0.9905

PBE/PW C1

RHF/various, B3LYP/various cluster

B3LYP/various C1

0.980

0.962-0.993

0.969

PBE/PW P1 C1 0.974

0.9873

0.975

0.963-0.987

0.964

0.969

0.965 0.962

0.970 0.969

0.99

0.976 0.975

a Values given for each of the three inner surface hydroxyls in each study should be compared with all other inner surface hydroxyls from other studies as it is not reported which value corresponds to which specific inner surface hydroxyl group.

The PBE functional was selected because the authors believed it to be particularly accurate for the calculation of hydrogen bonding. The study of Tosoni et al.19 presents geometry optimization of kaolinite using the B3LYP hybrid method, stating that B3LYP gives more reliable results than other DFT methods when dealing with OH group vibrational properties, especially when hydrogen bonds are present. However, it was shown by Paier et al.46 that B3LYP performs clearly worse than other hybrid functionals for periodic systems. Balan et al.15 calculated the theoretical infrared spectrum of kaolinite using PBE functional and plane-wave basis set, which is very similar to the method used by Hu and Michaelides.45 However, during the full geometry optimization, cell parameters were fixed to experimental values, and for the final optimization only hydrogen atoms were allowed to relax. Similar to the technique used with neutron diffraction patterns in the Results and Discussion section of this paper, Balan et al.15 used simulated X-ray powder diffraction patterns to validate their computations. However, the studies of Hu and Michaelides45 and Balan et al.15 did not assess a range of functionals in order to find the most appropriate choice. Hess and Saunders13 used a periodic ab initio quantum chemical method to predict the position of the inner hydrogen atoms in the crystal structure of kaolinite. They used a Hartree-Fock program for periodic systems, CRYSTAL, with the STO-3G basis set to compute the optimal position of the inner hydrogen atoms. The STO-3G and 6-21G basis sets gave the same results qualitatively, although the localized nature of the minimal STO-3G basis set may not properly account for long-range forces. This study displays the issues associated with

computational methods 15 years ago, primarily due to the restriction in computing power. Hobbs et al.14 evaluated the kaolinite structure using periodic DFT calculations, with a local approximation for the functional (with Perdew-Zunger parametrization of the Ceperley-Alder electron gas results47). Two experimentally determined unit cells of kaolinite were geometrically optimized (Table 5): the structures of Young and Hewat12 (P1) and Bish9 (C1). As part of this process, the c dimensions of the cell parameters were adjusted manually in order to assess the energy changes; however, there was no significant change in the c dimension for the geometrically optimized kaolinite unit cell compared with the experimental values. The resulting geometrically optimized structures derived from the structures of Young and Hewat12 and Bish9 both possessed overall C1 symmetry. For comparison with the Hess and Saunders13 Hartree-Fock calculation, Hobbs et al.14 performed a H-atom geometry optimization of the Bish and von Dreele8 kaolinite structure while keeping the heavy atoms fixed. The energies associated with the inner OH positions were compared to results from Hess and Saunders,13 with Hobbs et al.14 obtaining different results for the inner OH angle, attributed mainly to improvements in the basis sets available at the time of each investigation. In other investigations, Castro and Martins20 investigated the surface of kaolinite by employing a cluster model and the ONIOM method. Benco et al.16 used the GGA functional PW91 to investigate hydroxyl group orientations and hydrogen bond formation. Sato et al.5 used the PBE functional to investigate kaolinite and its bulk modulus, providing further insight into

Structure of Kaolinite

J. Phys. Chem. B, Vol. 113, No. 19, 2009 6761

TABLE 6: Comparison of Atom Positions between GGA-BLYP (Combination 1) and LDA-PWC (Combination 3)a GGA-BLYP (combination 1) (Å)

LDA-PWC (combination 3) (Å)

atom name (Figure 1)

a

b

c

a

b

c

∆r (Å)

Al1 Al4 Si1 Si4 O1 O2 O3 O4 O5 O11 O12 O13 O14 H1 H2 H3 H8

1.867 4.402 0.281 2.923 0.582 0.915 0.401 1.306 1.276 0.588 0.148 0.505 0.495 1.096 0.601 0.595 3.129

4.411 2.945 2.990 1.460 3.125 5.911 4.429 1.938 6.888 8.675 1.483 4.248 7.619 0.583 1.523 4.453 2.817

3.396 3.379 0.498 0.511 2.167 2.175 -0.212 -0.052 -0.212 2.259 4.336 4.337 4.378 2.346 5.320 5.313 5.290

1.860 4.390 0.298 2.942 0.573 0.932 0.358 1.358 1.322 0.572 0.101 0.502 0.504 1.082 0.609 0.557 3.154

4.397 2.934 3.009 1.476 3.137 5.907 4.448 2.014 6.856 8.649 1.473 4.208 7.632 0.565 1.525 4.455 2.824

3.388 3.368 0.515 0.527 2.171 2.179 -0.193 -0.029 -0.194 2.269 4.317 4.316 4.371 2.356 5.294 5.284 5.294

0.018 0.020 0.031 0.030 0.016 0.018 0.051 0.095 0.060 0.032 0.052 0.046 0.017 0.025 0.028 0.048 0.027

a

If the distance between a set of equivalent atoms is greater than 0.01 Å, then the two functionals result in different kaolinite structures (i.e., beyond convergence tolerance). If the distance is less than 0.01 Å for all sets of atoms, then the two structures must be considered indistinguishable at this level of computational accuracy.

the mechanical properties of clay minerals which can be difficult to obtain experimentally. Due to the large variation in application of computational procedures to the study of kaolinite, it is difficult to critically evaluate the success of each study via a simple summary such as Table 5. Nevertheless, some studies are more thorough than others, and advancements in basis sets and functionals have enabled more and more detailed calculations as time has passed. Outcomes of DFT Computations. Analysis of all the geometrically optimized structures of kaolinite provided conclusive evidence that the space group of this structure is C1 (see Supplementary Data 2 in the Supporting Information). Starting a refinement from a primitive P1 unit cell, the atoms move to their energy-minimizing positions, which can then be analyzed to see if they do in fact display some degree of symmetry. In each geometrically optimized structure, the majority of the atoms within the unit cell obey the translational operator for C1 to within the tolerance of the computations (0.005 Å). For those atoms which fall outside the tolerance (at most four atoms out of O13, O14, H1, H3, H8), the deviation is still very small (0.01 Å?

2 3 4

LDA-VWN LDA-PWC GGA-PW91

yes yes yes

5

GGA-BP

yes

6 7 8

GGA-PBE GGA-BOP GGA-VWN-BP

yes yes yes

9 10

GGA-RPBE GGA-HTCH

yes yes

which atoms? all all Si1, Si4, O3, O4, O5, O11, O12, O13, H1, H2, H3 Si1, Si4, O3, O4, O5, O11, H1, H2, H3 O3, O4, O5, O11, H3 O4 Si1, Si4, O3, O4, O5, O11, H1, H2, H3 O11, O12, O13, H8 all except Al4, O1, O11

that the computation using the LDA-PWC functional does give a significantly different kaolinite structure to the structure given by combination 1 (GGA-BLYP), since the displacement between equivalent atoms does exceed 0.01 Å for at least one pair of equivalent atoms (in fact, for this comparison all displacements exceed 0.01 Å). Table 7 gives a summary of the comparison of atom positions for combinations 1-10, which were achieved by undertaking the same process as shown in Table 6. From these results it is apparent that each functional gives a different structure for kaolinite. Hence, there is the need for experimental validation in order to select the most appropriate functional. Comparison of the simulated neutron diffractograms produced from the computationally derived kaolinite unit cells and the refined unit cell obtained by Bish9 has been carried out for all combinations given in Table 1. As with any refined structure for such a complex phase, the results obtained by Bish9 cannot be guaranteed to be a perfect representation of the true structure. However, as this is the best available experimentally derived structure, it will be used as the key point of comparison. This structure provides a very close fit to the high-quality experimental neutron diffraction data used to refine the structural parameters. Comparing the diffractograms of computed structures with the simulated diffractogram of a high-quality refined structure (rather than directly with experimental data) provides

6762

J. Phys. Chem. B, Vol. 113, No. 19, 2009

White et al.

Figure 2. Simulated neutron diffractograms for the computationally derived kaolinite unit cell (GGA-BLYP) and the unit cell obtained by Bish,9 using λ ) 1.9102 Å. Lattice parameters for the computed structure neutron diffractograms were fixed to the values refined by Young and Hewat12 to avoid differences due to Bragg peak offsets.

TABLE 8: Least-Squares Error Analyses of Neutron Diffractograms, Comparing Computed Kaolinite Unit Cells with Experimental Kaolinite Unit Cell Determined by Bish9 combination

functional

sum of squared differences

1 2 3 4 5 6 7 8 9 10

GGA-BLYP LDA-VWN LDA-PWC GGA-PW91 GGA-BP GGA-PBE GGA-BOP GGA-VWN-BP GGA-RPBE GGA-HTCH

1864 2179 2167 1928 1935 1921 1905 1938 1967 2211

χ2

Rp

Rwp

0.1623 0.1670 0.1664 0.1613 0.1651 0.1634 0.1675 0.1651 0.1718 0.1736

0.2334 0.2340 0.2337 0.2317 0.2340 0.2334 0.2364 0.2339 0.2385 0.2387

0.3538 0.3564 0.3558 0.3520 0.3565 0.3546 0.3596 0.3565 0.3643 0.3642

the significant advantage that the experimental background (which is significant in this case) need not be considered. Examples of the calculated neutron diffractograms are given in Figure 2. Least-squares error analysis results are given in Table 8. From these results, if χ2, Rp, and Rwp are used to assess the accuracy of the computed structure then GGA-PW91 (combination number 4) would be the most accurate structure. However, the sum of the squared differences for this combination is relatively large compared with other combinations (1928 cf. 1864 for GGA-BLYP); hence this computed structure cannot so straightforwardly be considered the most accurate. The reason for the sum of the squared differences not following the same trend as χ2, Rp, and Rwp for combination 4 is due to the absence of the weighting factor wi in this equation and the process of squaring the difference (eq 1). Without the weighting factor and with squaring, the high intensity peaks dominate. Hence, this term provides an analysis method which concentrates on the high-intensity peaks and therefore presents an important way to measure the contributions of these peaks. Combination 4 is not providing an accurate fit to these important high-intensity peaks. Of the remaining combinations, the most accurate computed kaolinite structure is given by GGA-BLYP, with GGA-PBE the next most accurate. Even though GGA-VWN-BP and GGABP report lower total energies than GGA-BLYP (Table 9), these

TABLE 9: Total Energies of Geometry Optimized Kaolinite Structures combination

functional

total energy (hartrees)

1 2 3 4 5 6 7 8 9 10

GGA-BLYP LDA-PWC LDA-VWN GGA-PW91 GGA-BP GGA-PBE GGA-BOP GGA-VWN-BP GGA-RPBE GGA-HCTH

-3489.05 -3471.18 -3471.26 -3488.75 -3489.06 -3486.72 -3488.65 -3489.16 -3488.76 -3488.80

functionals do not perform as well in replicating the neutron diffraction pattern. Hence, total energy should not be the only method of assessing functional performance, as it does not always result in a more accurate structure. The local spin-density approximation functionals (LDA-PWC and LDA-VWN) produce the highest total energies, which is in agreement with previous studies comparing total energies produced by LDA and GGA functionals.48 This shows that LDA functionals cannot accurately reproduce total energies, giving markedly higher energies than the more advanced functionals based on the GGA approximation. Hence, taking into account total energies, equivalent atom positions, and comparison of neutron diffractograms, the most accurate functional to use for density functional calculations of kaolinite among the candidates tested is the generalized gradient approximation (GGA) functional BLYP. Structure of Kaolinite. As discussed previously, the most appropriate (accurate and efficient) functional choice for kaolinite when using DMol3 is GGA-BLYP. The resulting geometrically optimized structure is given in Table 10, which displays fractional coordinates of the seventeen atoms defining the C1 unit cell of kaolinite. Atom names correspond to those given in Figure 1. This structure compares well with previous structures of kaolinite provided in the literature. Hydroxyl group orientations and bond lengths are given in Table 11.

Structure of Kaolinite

J. Phys. Chem. B, Vol. 113, No. 19, 2009 6763

TABLE 10: Fractional Atomic Coordinates of Kaolinite C1 Structure Geometrically Optimized Using DMol3 (GGA-BLYP) atom name (Figure 1)

fractional a

fractional b

fractional c

Al1 Al4 Si1 Si4 O1 O2 O3 O4 O5 O11 O12 O13 O14 H1 H2 H3 H8

0.36266 0.85492 0.05464 0.56769 0.11303 0.17770 0.07793 0.25365 0.24774 0.11411 0.02877 0.09810 0.09615 0.21276 0.11667 0.11558 0.60760

0.49376 0.32966 0.33468 0.16344 0.34982 0.66164 0.49580 0.21699 0.77104 0.97111 0.16604 0.47556 0.85284 0.06527 0.17043 0.49851 0.31528

0.45990 0.45758 0.06743 0.06914 0.29339 0.29455 -0.02869 -0.00707 -0.02877 0.30590 0.58718 0.58733 0.59286 0.31768 0.72050 0.71954 0.71641

TABLE 11: Hydroxyl Bond Lengths and Orientations of DMol3 Geometrically Optimized Kaolinite Structure (GGA-BLYP)a atom pair (Figure 1)

bond length (Å)

bond orientation

O(11)-H(1) O(12)-H(2) O(13)-H(3) O(18)-H(8)

0.974 0.972 0.970 0.970

5.057 78.272 84.585 68.582

a

All angles are given in degrees relative to the a,b plane.

When comparing the computed hydroxyl group orientations (Table 11) with experimentally determined structure in the literature (Table 2), the computed hydroxyl groups orientations match best with the data of Bish,9 with the largest difference in hydroxyl bond orientation being approximately 23°. However, this difference is far from insignificant, and reasons for such a large difference can be attributed to two possibilities. The first possibility is that the approximations inherent in DFT calculations may lead to imprecise results. These approximations and subsequent errors in bond orientations and bond lengths have been well documented in the literature.49-51 The second possibility is that errors were induced in treatment of the experimental data, either during the experiment or during the refinement of the data. Rietveld refinement is a complex procedure and authors have reported the need to impose various constraints when refining the structure of kaolinite, which can bring about errors. Without further investigation and knowledge of the exact method undertaken during Rietveld refinement to obtain the Bish structure,9 it is difficult to specify the exact cause of this difference in bond orientations. Previously reported theoretically determined hydroxyl bond orientations (Table 4) are similar to those computed in this study (Table 11), apart from the structure of Hess and Saunders13 who kept the bond orientations of the inner surface hydroxyls fixed at 57.7° relative to the a,b plane. The study that is closest to those reported in Table 11 is that of Hobbs et al.,14 which was based on the density functional computational method using a local approximation for the functional (with Perdew-Zunger parametrization of the Ceperley-Alder electron gas results47). Hence, the study of Hobbs et al.14 is similar to the study reported here in terms of the approach to refining the structure of kaolinite. However, the study presented here has shown that the GGA class of functionals provide a more realistic structure

than the LDA functional as used by Hobbs et al.,14 without excessive increases in computational requirements. This can be seen in Table 1, where GGA functionals took approximately double the time compared to the LDA functionals, which is not unreasonable given the accuracy achieved by the majority of GGA functionals compared with LDA functionals (Table 8). The experimental and theoretical reported values for hydroxyl bond lengths (Tables 3 and 5) vary widely from 0.75 to 1.45 Å, with the values calculated in this study falling between 0.970 and 0.974 Å (Table 11). As was the case for computed hydroxyl bond orientations when compared to experimental values in literature, the values computed in this study agree most closely with the experimental values of Bish9 in Table 3, which reaffirms the choice of this experimental structure as the structure to use for validation of computations. The theoretically determined hydroxyl bond lengths do not vary as widely as the experimentally determined lengths (0.962-0.993 Å, Table 5), which, due to the ability to provide realistic structural characteristics, shows the advantage of using DFT techniques to investigate the structural properties of such materials. The values calculated in this study are also in good agreement with the DFT computations of Hu and Michaelides,45 Balan et al.,15 and Tosoni et al.19 The hydroxyl bond lengths obtained in this study are within 0.003 Å of the lengths obtained by Hu and Michaelides,45 and within 0.006 Å of the results of Balan et al.15 Tosoni et al.19 carried out geometry optimization of kaolinite using a polarized double-ζ Gaussian basis set, B3LYP functional, and a very fine convergence tolerance of 10-8 hartree (cf. 10-5 hartree in this study). The bond lengths calculated here are systematically slightly longer than those of Tosoni et al.,19 which may be attributed to overestimation of bond lengths by the hybrid functional B3LYP used by those authors. Finally, the dynamical behavior of the computed kaolinite structure was used as additional method to assess the accuracy of density functional calculations. Balan et al.15 have previously used experimental results from infrared and Raman spectroscopy as the method of assessing the accuracy of their computations on the kaolinite structure, and achieved acceptable agreement between theory and experiment. However, these authors compared mainly vibrational frequencies and because of the low symmetry of the cell there are almost 100 fundamental normal modes. We compare our calculated frequencies (DMol3 and VASP) with those of Balan et al.15 in Supplementary Data 4 and it is clear that with so many frequencies plausible agreement can always be found. In an H-bonded system with anharmonicity, overtones and combinations will also arise and the calculation of IR intensities becomes intractable. The best solution is to use the inelastic neutron spectrum (INS) which is relatively easy to calculate, including overtones and combinations from the eigenvectors and eigenvalues of the dynamical matrix.39 In the present case this is important for the hydroxyl groups since INS is primarily sensitive to modes in which H-atoms are displaced. It is worth remembering that with overtones and combinations there are several hundred spectral peaks that merge into a spectral profile. Observed and calculated spectral profiles are compared in Figure 3. The region between 200 and 1100 cm-1 is important for the orientations of the hydroxyl groups and it is clear that DMol3 gives rather poor agreement in this region. Using essentially the same structure, VASP produces a much better spectrum (note nothing is fitted), which is rather surprising given that previous work has shown that DFT packages basically give very similar results when used at the same accuracy, with the same finite displacements. Calculated frequencies from VASP are included as tick marks

6764

J. Phys. Chem. B, Vol. 113, No. 19, 2009

Figure 3. Observed and calculated INS spectra. The experimental curve was taken from Johnston et al.,35 and the calculations were made from the crystal structure without adjustment or fitting of any parameters. The tick marks indicate the calculated frequencies of the fundamental vibrational modes using VASP (third curve from the top). It should be noted that the minimized structures using DMol3 and VASP obtained exactly the same structure to within convergence tolerance (see Supplementary Data 3 for comparison).

in Figure 3. The main difference between the observed spectrum and that calculated with VASP is the low-frequency doublet that is calculated about 100 cm-1 is too high. This discrepancy almost certainly arises from dispersion that is expected to arise in this region for hydrogen-bonded systems such as kaolinite. The model used here assumes that all vibrations in the unit cell are exactly in phase with all other unit cells, the so-called Γ-point. Interactions via the H-bonds cause the frequencies to vary depending on the phase relationship between neighboring cells. This can be calculated, see for example Kearley et al.,52 but this is a lot of work for such a low-symmetry system and would only be worthwhile if a high-resolution INS spectrum was available. Interestingly, the INS spectrum calculated from an ab initio MD simulation (VASP at low precision) does indeed show much better agreement for the low-frequency doublet (Figure 3) presumably because it can access much more conformational space than the finite displacement (typically 0.01 Å) used to construct the dynamical matrix. In summary, the vibrational spectrum calculated from the best structure (albeit with VASP, not DMol3) is not inconsistent with the measured INS spectrum in the spectral region where hydroxyl orientations are important. A high-resolution INS spectrum and a full lattice dynamics calculation are required to fully establish consistency and this work is currently being planned. Conclusions Density functional modeling of the unit cell structure of kaolinite has been performed in order to attain an accurate computationally derived structure of this important crystalline clay mineral and hence address uncertainties regarding the exact geometry of the hydroxyl groups in this structure. However, density functional modeling such as DMol3 requires user input with choice of exchange-correlation functional and other parameters. Hence, computations were carried out using various functionals, which were assessed in terms of total energies and equivalent atom positions, and compared with an experimentally determined kaolinite structure using simulated neutron diffractograms. The computed structure that agreed best with the experimental structure was generated by using the generalized gradient approximation BLYP functional. This computation confirmed that kaolinite is of space group C1. The inner

White et al. hydroxyl group is nearly parallel to the a,b plane (5.057°), whereas the inner surface hydroxyl groups are closer to perpendicular (78.272°, 84.585°, 68.582°). Hydroxyl bond lengths vary between 0.970 and 0.974 Å. This therefore represents a computationally inexpensive method whereby an accurate representation of the kaolinite structure can be generated for use in larger-scale simulations. Comparison of this computed structure of kaolinite with those given in the literature showed that the experimental structure that matched best was the one determined by Bish.9 However, there are still significant differences between the hydroxyl group orientations and bond lengths. The computed structure does match better with theoretical kaolinite structures presented in the literature, particularly those presented by Hu and Michaelides,45 Balan et al.,15 and Tosoni et al.19 Such differences between experimental and theoretical results can be attributed to approximations inherent in the theoretical calculations, and difficulties and systematic errors associated with experiments. Calculation of inelastic neutron spectra using DMol3 did not give satisfactory results, although the VASP package did generate reasonable agreement with experimental data for a structure that was identical to within convergence tolerance, and this is an area which requires further investigation. From the density functional calculations performed on the structure of kaolinite, it can be concluded that this method can generate accurate structural properties for reasonably complex and lowsymmetry layered aluminosilicate structures, although experimental validation is necessary to ensure that the results are meaningful. Acknowledgment. This work was funded in part by the Australian Research Council (ARC) (including some funding via the Particulate Fluids Processing Centre, a Special Research Centre of the ARC), in part by a University of MelbourneANSTO Collaborative grant, and in part by a studentship paid to C.E.W. by the Centre for Sustainable Resource Processing via the Geopolymer Alliance. The authors acknowledge Dr. Margaret Elcombe for generously sharing her unpublished neutron powder diffraction data for kaolinite, and Dr. Kia Wallwork (Australian Synchrotron) for assistance in collecting the powder diffraction data for kaolinite at low temperatures. Supporting Information Available: The linear extrapolations performed on the lattice parameters determined from synchrotron X-ray powder diffraction is available (Supplementary Data 1). A comparison of calculated atom positions showing the close adherence to the C1 space group when refining a structure from P1 symmetry is also available (Supplementary Data 2). The comparison of the atomic coordinates of the final DMol3 and VASP structures is given in Supplementary Data 3, along with a table comparing the observed and calculated vibrational frequencies (Supplementary Data 4). This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Murray, H. H. Appl. Clay Sci. 2000, 17, 207. (2) Murray, H. H. Clays. In Ullmann’s Encyclopedia of Industrial Chemistry; John Wiley & Sons, Inc.: Weinheim, Germany, 2006. (3) Richardson, I. G. Cem. Concr. Res. 1999, 29, 1131. (4) Duxson, P.; Ferna´ndez-Jime´nez, A.; Provis, J. L.; Lukey, G. C.; Palomo, A.; van Deventer, J. S. J. J. Mater. Sci. 2007, 42, 2917. (5) Sato, H.; Ono, K.; Johnston, C. T.; Yamagishi, A. Am. Mineral. 2004, 89, 1581. (6) Giese, R. F.; Datta, P. Am. Mineral. 1973, 58, 471. (7) Akiba, E.; Hayakawa, H.; Hayashi, S.; Miyawaki, R.; Tomura, S.; Shibasaki, Y.; Izumi, F.; Asano, H.; Kamiyama, T. Clays Clay Miner. 1997, 45, 781.

Structure of Kaolinite (8) Bish, D. L.; von Dreele, R. B. Clays Clay Miner. 1989, 37, 289. (9) Bish, D. L. Clays Clay Miner. 1993, 41, 738. (10) Neder, R. B.; Burghammer, M.; Grasl, T.; Schulz, H.; Bram, A.; Fiedler, S. Clays Clay Miner. 1999, 47, 487. (11) Suitch, P. R.; Young, R. A. Clays Clay Miner. 1983, 31, 357. (12) Young, R. A.; Hewat, A. W. Clays Clay Miner. 1988, 36, 225. (13) Hess, A. C.; Saunders, V. R. J. Phys. Chem. 1992, 96, 4367. (14) Hobbs, J. D.; Cygan, R. T.; Nagy, K. L.; Schultz, P. A.; Sears, M. P. Am. Mineral. 1997, 82, 657. (15) Balan, E.; Saitta, A. M.; Mauri, F.; Calas, G. Am. Mineral. 2001, 86, 1321. (16) Benco, L.; Tunega, D.; Hafner, J.; Lischka, H. Am. Mineral. 2001, 86, 1057. (17) Benco, L.; Tunega, D.; Hafner, J.; Lischka, H. J. Phys. Chem. B 2001, 105, 10812. (18) Sato, H.; Ono, K.; Johnston, C. T.; Yamagishi, A. Am. Mineral. 2005, 90, 1824. (19) Tosoni, S.; Doll, K.; Ugliengo, P. Chem. Mater. 2006, 18, 2135. (20) Castro, E. A. S.; Martins, J. B. L. Int. J. Quantum Chem. 2005, 103, 550. (21) Frost, R. L. Clays Clay Miner. 1995, 43, 191. (22) Frost, R. L.; Fredericks, P. M.; Kloprogge, J. T.; Hope, G. A. J. Raman Spectrosc. 2001, 32, 657. (23) Frost, R. L.; Kloprogge, J. T. Spectrochim. Acta, Part A 2001, 57, 163. (24) Grimme, S. J. Comput. Chem. 2004, 25, 1463. (25) Johnson, M. R.; Prager, M.; Grimm, H.; Neumann, M. A.; Kearley, G. J.; Wilson, C. C. Chem. Phys. 1999, 244, 49. (26) Delley, B. J. Chem. Phys. 1990, 92, 508. (27) Perdew, J. P.; Wang, Y. Phys. ReV. B 1992, 45, 13244. (28) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200. (29) Becke, A. D. J. Chem. Phys. 1988, 88, 2547. (30) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1996, 77, 3865. (31) Tsuneda, T.; Suzumura, T.; Hirao, K. J. Chem. Phys. 1999, 110, 10664.

J. Phys. Chem. B, Vol. 113, No. 19, 2009 6765 (32) Hammer, B.; Hansen, L. B.; Norskov, J. K. Phys. ReV. B 1999, 59, 7413. (33) Boese, A. D.; Handy, N. C. J. Chem. Phys. 2001, 114, 5497. (34) Hunter, B. A.; Howard, C. J. A computer program for Rietveld analysis of X-ray and neutron powder diffraction patterns; Lucas Heights Laboratories, Australia, 1998. (35) Johnston, C. T.; Bish, D. L.; Eckert, J.; Brown, L. A. J. Phys. Chem. B 2000, 104, 8080. (36) Kresse, G.; Furthmu¨ller, J. Software VASP 1999. (37) Kresse, G.; Furthmuller, J. Phys. ReV. B 1996, 54, 11169. (38) Kresse, G.; Furthmuller, J. Comput. Mater. Sci. 1996, 6, 15. (39) Kearley, G. J. A review of the analysis of molecular vibrations using INS. Proceedings of the Third Workshop on Neutron Scattering Data Analysis, Chilton, Didcot, UK; 1995. (40) Kneller, G. R.; Keiner, V.; Kneller, M.; Schiller, M. Comput. Phys. Commun. 1995, 91, 191. (41) Wallwork, K. S.; Kennedy, B. J.; Wang, D. AIP Conf. Proc. 2007, 879. (42) Adams, J. M.; Hewat, A. W. Clays Clay Miner. 1981, 29, 316. (43) Rietveld, H. M. J. Appl. Crystallogr. 1969, 2, 65. (44) Rietveld, H. M. Acta Crystallogr. 1967, 22, 151. (45) Hu, X. L.; Michaelides, A. Surf. Sci. 2008, 602, 960. (46) Paier, J.; Marsman, M.; Kresse, G. J. Chem. Phys. 2007, 127. (47) Perdew, J. P.; Zunger, A. Phys. ReV. B 1981, 23, 5048. (48) Corso, A. D.; Pasquarello, A.; Baldereschi, A.; Car, R. Phys. ReV. B: Condens. Matter 1996, 53, 1180. (49) Ireta, J.; Neugebauer, J.; Scheffler, M. J. Phys. Chem. A 2004, 108, 5692. (50) Hamann, D. R. Phys. ReV. Lett. 1996, 76, 660. (51) Demuth, T.; Jeanvoine, Y.; Hafner, J.; Angyan, J., G. J. Phys.: Condens. Matter 1999, 11, 3833. (52) Kearley, G. J.; Johnson, M. R.; Tomkinson, J. J. Chem. Phys. 2006, 124, 044514.

JP810448T