10236
J. Phys. Chem. B 2001, 105, 10236-10242
An ab Initio Study of the Structure and Properties of Aluminum Hydroxide: Gibbsite and Bayerite Julian D. Gale,*,† Andrew L. Rohl,‡ Victor Milman,§ and Michele C. Warren| Department of Chemistry, Imperial College of Science, Technology and Medicine, South Kensington, SW7 2AY, U.K., School of Applied Chemistry, Curtin UniVersity of Technology, P.O. Box U1987, Perth, Western Australia 6845, Molecular Simulations, Inc., The Quorum, Barnwell Road, Cambridge, CB5 8RE, U.K., and Department of Earth Sciences, UniVersity of Manchester, Oxford Road, Manchester, M13 9PL, U.K. ReceiVed: May 10, 2001; In Final Form: June 28, 2001
The two most important polymorphs of aluminum hydroxide, namely gibbsite and bayerite, have been studied for the first time using solid state ab initio quantum mechanical methods, both using plane wave and localized basis set methodologies, within the framework of nonlocal density functional theory. The fully optimized structures have been determined for both phases, yielding improved hydrogen positions in the case of gibbsite for which the only previous information is from X-ray data. Mechanical properties have been calculated for gibbsite, including the full elastic constants tensor and the bulk modulus. The latter is found to be 55 GPa, which is significantly lower than a recent experimental estimate. Vibrational spectra have been calculated for both phases and assignments of the hydroxyl stretching modes are proposed.
Introduction The study of gibbsite has significance from both a mineralogical and commercial point of view. Gibbsite is an important ore of aluminum and is one of three minerals that make up bauxite, the raw material for the Bayer process for the extraction of aluminum. Bauxite is composed of aluminum oxide and hydroxide minerals such as gibbsite, boehmite, AlO(OH) and diaspore, AlO(OH), as well as clays, silt, and iron oxides, and hydroxides. Because the aluminum minerals are soluble in strong caustic solutions, they are separated from the rest of the ore by addition of sodium hydroxide. In the rate-limiting step of the Bayer process, they are precipitated out from the caustic solution as gibbsite. Aluminum hydroxide (Al(OH)3) has three other polymorphs; bayerite, nordstrandite, and doyleite,1-3 although bayerite is the only one with any industrial significance, since it is an undesirable product in various stages of the Bayer process. In addition, bayerite incorporates less sodium during its precipitation than gibbsite and thus is used commercially to prepare aluminum oxide catalyst supports.4 All of the polymorphs of aluminum hydroxide are composed of layers of aluminum octahedra with hydroxyl groups on either side that hydrogen bond the layers together. The difference between the polymorphs is the stacking sequence of the layers. Gibbsite and bayerite represent the two extremes of the stacking sequence, with nordstrandite and doyleite being intermediate structures. The structures of gibbsite and bayerite, both of which adopt the space group P21/n, are illustrated in Figure 1. If we label the stacking of two oxygen layers that contain the aluminum ions as AB, the figure shows that the bayerite structure consists of these double layers stacked directly on top * To whom correspondence should be addressed. Email:
[email protected]. FAX: +44 207 594 5804. † Department of Chemistry. ‡ School of Applied Chemistry. § Molecular Simulations, Inc.. | Department of Earth Sciences.
Figure 1. The structure of the two main polymorphs of aluminum hydroxide, (a) bayerite and (b) gibbsite, illustrating the difference in the stacking of the layers.
of each other resulting in an ABABABAB... structure. In the other extreme, gibbsite, each double layer is a reflection of the previous one, resulting in an ABBAABBA... stacking motif. Also note that the figure clearly shows that only two-thirds of the octahedral sites are occupied by aluminum ions in each double layer. Although the basic structures of all the polymorphs are known, there is still much to be characterized. From a structural perspective, the most widely used crystal structure for gibbsite is based on X-ray diffraction,5 thus leading to poor characterization of the hydroxyl groups. Beyond the structure, there is still much to be determined concerning the physical properties and relative stabilities of the polymorphs where there are either gaps in the experimental information or uncertainties concerning some of the values. Experimental determinations are, in general, hindered by difficulties in obtaining pure samples with large crystal dimensions. There have been a small number of previous theoretical attempts to study gibbsite using both interatomic potentials6,7 and quantum mechanics.8 However, the quantum mechanical work to date has been based on the cluster approach which limits
10.1021/jp011795e CCC: $20.00 © 2001 American Chemical Society Published on Web 09/28/2001
Properties of Aluminum Hydroxide its ability to distinguish between different polymorphs of aluminum hydroxide and to determine the structure in a manner that properly acccounts for the intrinsic periodicity. The aim of the present work is to perform the first study, to the best of our knowledge, of the energetics, structure, and properties of gibbsite and bayerite, as the most significant polymorphs of aluminum hydroxide, using ab initio methods within periodic boundary conditions. Methods In the present study, two different approaches have been used for the calculation of the properties of gibbsite and bayerite. Both are based on the use of density functional theory,9 which represents, at the present time and for systems of this complexity, the only practical quantum mechanical method for studying materials properties with the inclusion of electron correlation. All calculations have been performed using the same gradientcorrected functional, namely, that of Perdew, Burke, and Ernzerhof (PBE).10 The two approaches taken differ only in the basis sets and algorithms used in order to evaluate the electronic structure. On one hand, we have used plane waves to expand the wave functions in the solid state,11 whereas in the second case a localized real space basis was employed. Plane waves have the distinct advantage of being able to systematically converge the basis set by purely increasing the kinetic energy cutoff, while the real space basis set has the advantage of being able to tackle larger system sizes. Hence one of the aims of this study is to validate the real space approach for future studies of more complex systems, including surfaces and defects. Plane wave calculations were performed using the program CASTEP12 with a kinetic energy cutoff of 340 eV. The nuclei and core electrons were represented by ultrasoft pseudopotentials of the form described by Vanderbilt13 on the basis of valence electron configurations of 3s23p1, 2s22p6, and 1s1, for Al, O, and H, respectively. Pseudopotentials were generated based on all electron calculations utilizing both the LDA and PBE functionals for the exchange-correlation potential. The Brillouin zone was sampled using a Monkhorst-Pack grid14 with shrinking factors of 2 × 4 × 2 for gibbsite and 4 × 2 × 2 for bayerite. As a complement to the plane wave calculations, we have also utilized a methodology which is suitable for the study of much larger systems, as embodied in the program SIESTA.15-17 In particular, the approach taken is designed to yield improved scaling with system size which can be made to be linear, via a suitable choice of algorithm, when the prefactor makes this favorable.18 For the size of the unit cell being considered here, the use of linear scaling techniques for solving the SCF procedure is not beneficial, relative to conventional diagonalization. However, for future studies on more complex environments, this situation will be reversed. In the approach taken within SIESTA, calculations are performed using a localized basis set which consists of numerical tabulations of the exact solutions to the pseudo-atomic problem,19 again with a pseudopotential representation of the core electrons and nucleus. Furthermore, the orbitals can be radially confined in a balanced manner across all species by specifying the energy increase that is to be allowed. By adjusting this energy shift, an appropriate compromise between convergence to the unconfined results and computational speed can be achieved. Here an energy shift of 0.005 Ry was found to be a suitable choice. This spatial confinement also leads to the improved scaling with respect to system size. All atoms are represented by the Troullier-Martins norm-conserving pseudo-
J. Phys. Chem. B, Vol. 105, No. 42, 2001 10237 potentials.20 The basis set was extended for the calculations to be of the double-ζ form, with a norm of 0.15 for the outer exponent, and polarization functions were utilized for all atoms, including hydrogens. For the evaluation of the Hartree and exchange-correlation potentials a real space mesh is used with an associated kinetic energy cutoff. For the calculations in the present work, a mesh cutoff of 250 Ry was found to offer a high level of convergence. Identical Monkhorst-Pack k-point grids to those employed in the plane wave calculation were used in the real space basis case to approximate the integration across the Brillouin zone. For the polymorphs of aluminum hydroxide considered here, the structures were fully relaxed, including the cell parameters, in the athermal limit. The gamma point phonons have been determined using SIESTA by central finite differencing of the analytical first derivatives with a displacement of 0.04 Bohr about the optimized configuration. Given the size of the unit cell for these materials, the phonons at the gamma point are expected to be representative of the full spectra. Furthermore, the principal interest is in comparison of results with infrared and Raman spectra, which primarily explore low Q values. In the future, the application of linear response techniques21 will make it more facile to access phonons in other regions of the Brillouin zone. In addition, the mechanical properties of gibbsite have been predicted with the plane wave methodology. Finite differences are again used, in this case in the strain, and from calculation of the resulting stress, the full elastic constants tensor may be determined for the relaxed structure.22,23 Since this tensor describes the full set of second derivatives of the energy with respect to the strain, other properties such as the Young’s moduli and Poisson ratios can then easily be determined. Importantly, it is necessary to relax the internal degrees of freedom for each point of the finite differences with respect to applied strain. This approach has already been successfully applied to low symmetry systems: the elastic constants of two orthorhombic polymorphs were accurately reproduced and the values for a third, triclinic, polymorph predicted ahead of experiment.24 For gibbsite, we have used four different strain patterns, each with four amplitudes between -0.3 and +0.3%. For each strain pattern, the symmetry and k-points were adjusted to reflect the reduced symmetry of the distorted structure. For each of the stress components generated as a result of each strain pattern, a leastsquares fit of stress against strain yielded one of the elastic constants. The bulk modulus of gibbsite is also available directly from the elastic constants tensor. However, for comparison we have also determined the value by fitting the pressure-volume curve to a third-order Birch-Murnaghan equation of state over the pressure range 0-10 GPa. Prior to the present work, much of the data known concerning the physical properties of gibbsite had only been determined theoretically, using interatomic potential models. However, given many of these were parametrized largely on the basis of structural data, rather than with the benefit of curvature related information as would normally be desirable, it is interesting to evaluate how such models fare against the ab initio data. Hence we present comparisons for the mechanical properties as determined according to shell model interatomic potentials within the framework of an ionic description. Here all calculations were performed using the program GULP25 which determines such properties analytically. The particular parameters used are those fitted by Fleming et al7 specifically for aluminum hydroxide.
10238 J. Phys. Chem. B, Vol. 105, No. 42, 2001
Gale et al.
TABLE 1: Comparison of Experimental and Calculated Cell Parameters for Gibbsite and Bayerite According to Plane Wave Calculationsa gibbsite
bayerite
parameter exptl5 LDA/LDA LDA/PBE PBE/PBE exptl26 PBE/PBE a(Å) b(Å) c(Å) β(°) vol (Å3)
8.684 5.078 9.736 94.54 428.0
8.504 4.867 9.224 93.06 381.2
8.623 5.017 9.598 92.76 414.7
8.798 5.063 5.107 5.110 8.672 8.852 9.674 9.425 9.517 92.54 90.26 90.47 434.4 413.8 430.2
a In the case of gibbsite, results are included for the cases LDA/ LDA, LDA/PBE, and PBE/PBE, where the first functional was used in the generation of the pseudopotentials and the second for the actual optimization of the polymorph.
Results Bulk Structure. First of all, we consider the performance of the two density functional methodologies for the bulk structures, commencing with the plane wave technique, which will act as a reference due to the greater ease of systematic convergence with respect to the basis set. For both gibbsite and bayerite, their structures have been determined by diffraction methods, including the positions of the hydrogen atoms. However, in the case of gibbsite, the main structural determination was performed with X-rays5 and thus the O-H bond lengths are systematically underestimated, while for bayerite a deutrated sample was used in neutron diffraction.26 Using the pseudopotential plane wave method we have fully optimized the structures of gibbsite and bayerite, the final unit cell dimensions being compared against experiment in Table 1. Although we are principally concerned with the values determined using the GGA of Perdew, Burke and Ernzerhof, a number of other results are also given for gibbsite. Included are results for the cases where ultrasoft pseudopotentials generated based on all electron LDA reference calculations are used in a bulk calculation with the LDA (LDA/LDA) and PBE (LDA/PBE) functionals. These results serve to illustrate the consequences of utilizing pseudopotentials which are inconsistent with the final functional used for the solid-state calculation. Considering first of all the results for gibbsite, we see the usual trend for the LDA/LDA calculation to underestimate the unit cell volume, in this case dramatically by 10%. Because the structure consists of layers in the ab plane, the largest single error arises through an overcompression of the c axis by approximately 0.5 Å. This is due to the tendency of LDA to yield hydrogen bonds, which determine the interlayer spacing, that are too short,27 though the intralayer cell parameters are also contracted with respect to experiment. Performing a calculation using the PBE functional with LDA generated pseudopotentials leads to a cell volume which is too small by 3%, while the correct consistent PBE calculation gives an overestimation by only 1%. This is broadly consistent with what might be expected based on previous calculations, in that GGA functionals tend to overestimate cell volumes by a few percent, though the error here is a little lower than typical. Often the LDA/GGA combination is found to give fortuitously good results, since the LDA pseudopotentials counteract the overestimation of volume due to the GGA bulk calculation. However, for gibbsite this is not found to be the case. Bayerite, in contrast, demonstrates just under a 4% overestimation of the cell volume despite showing much better reproduction of the monoclinic angle than for gibbsite. This larger error can be attributed to the fact that in gibbsite the interlayer c axis is underestimated, which largely negates the fact that the unit cell is too large in
Figure 2. The labeling of the hydrogen/deuterium atoms within the unit cell of (a) bayerite and (b) gibbsite.
TABLE 2: Optimized Fractional Coordinates for the Asymmetric Unit Hydrogens in Gibbsite and Their Corresponding O-H Bond Lengths (in Angstroms) site
x
y
z
calcd r(O-H)
exptl r(O-H)5
H1 H2 H3 H4 H5 H6
0.0757 0.5729 0.4951 0.9497 0.2941 0.8058
0.1393 0.5548 0.1105 0.8156 0.7181 0.1624
0.8741 0.8977 0.7949 0.8866 0.7933 0.7970
0.974 0.981 0.984 0.981 0.986 0.988
0.749 0.773 0.838 0.875 0.883 0.861
the ab plane for both polymorphs, while in bayerite all dimensions are overestimated. Given that the hydrogen positions are poorly characterized for gibbsite, a further calculation has been performed in which the unit cell and heavy atom positions were held fixed at those from Saalfeld and Wedde5 and just the hydrogen positions were optimized. The numbering of the hydrogen/deuterium atoms within the structures of bayerite and gibbsite are illustrated in Figure 2. The calculated hydrogen fractional coordinates for gibbsite are given in Table 2 along with the hydroxyl bond lengths. As expected, the optimized hydrogen positions are much more consistent with typical hydroxyl bond lengths. Hence we believe that these fractional coordinates represent a more reliable set than those from the original experimental data and would be suitable for future structural refinement. In comparison, Kubicki and Apitz8 obtain O-H distances in the range 0.940.95 Å from their cluster calculation at the HF/3-21G(d,p) level. This systematic difference is to be expected given the tendency of Hartree-Fock calculations to underestimate bond lengths, while nonlocal density functional theory tends to lead to overestimation. Having determined precise reference data for gibbsite and bayerite with converged plane wave calculations, the localized basis set approach was used to calculate structures for gibbsite and bayerite. The results are given in Table 3 again compared to experimental data. As can be seen, all individual cell lengths are reproduced to better than 3%, while the volumes are overestimated by 4%. Generally the real space approach shows similar trends to the PBE/PBE plane wave calculation except that the overestimation of the volume is more uniform and slightly larger. Perhaps the most significant discrepancy is the monoclinic angle for gibbsite which is too small, although the agreement for bayerite is much better. In terms of the hydroxyl groups, the localized basis approach yields O-H bond lengths
Properties of Aluminum Hydroxide
J. Phys. Chem. B, Vol. 105, No. 42, 2001 10239
TABLE 3: Comparison of Experimental and Calculated Cell Parameters for Gibbsite and Bayerite Using the Localized Basis Set Methodology Embodied within SIESTA gibbsite
bayerite
parameter
exptl5
calcd
exptl26
calcd
a (Å) b (Å) c (Å) R (deg) β (deg) γ (deg)
8.684 5.078 9.736 90.00 94.54 90.00
8.887 5.171 9.696 90.00 91.47 90.00
5.063 8.672 9.425 90.00 90.26 90.00
5.149 8.915 9.422 90.00 90.49 90.00
that are systematically of the order of 0.01 Å longer than those from the plane wave calculations, though this discrepancy could be reduced by adding further basis functions. For calibration, the present computational conditions in SIESTA yield an O-H bond length of 0.972 Å (0.957 Å) and an angle of 104.3° (104.5°) for an isolated water molecule, where the experimental values are given in parentheses. Energetic Stability Perhaps the principal question in this area concerns the order of stability of the polymorphs. The general predominance of gibbsite in nature suggests that this may be the more stable form. Furthermore, bayerite slowly transforms into gibbsite on standing in contact with the appropriate solution at 300 K.28 Experimental calorimetric measurement of the difference in the heats of formation at 298 K also suggests that gibbsite is more stable by approximately 5 kJ/mol.4 In our calculations gibbsite is found to be lower in energy by 7.7 kJ/mol than bayerite at 0 K based on the plane wave calculations. A similar difference of 6.3 kJ/mol is obtained from the localized basis set approach, again confirming that both approaches are consistent with each other. The magnitude of the energy difference is small and near the limit of the likely accuracy of present density functionals, making it hard to be certain that the ordering is correct based on this alone. However, both measurements agree well with the experimental estimate thus giving added confidence. Furthermore, the similarity of both structures should lead to cancellation of systematic errors in the density functional used. It has been argued that bayerite should be the more stable phase due to the fact that it is more dense than gibbsite at ambient pressure. However, this logic is not conclusive, because there are counter examples of hydrogen bonded molecular crystalline systems in which the most dense polymorph is less stable.29 In such cases, the key criterion is optimization of the directionality of the interactions, rather than the straightforward number density. Beyond the difference in the internal energies of the materials, the influence of the thermodynamics of vibration should also be allowed for. As will be discussed later, the gamma point phonons have been determined for both gibbsite and bayerite using the real space approach. Hence we can readily calculate the contribution of the zero point energy to the stability difference. This term is found to be quite small at 1.2 kJ/mol and again favors gibbsite over bayerite. This additional preference for gibbsite can be readily understood, because the higher density of bayerite leads to a blue shift in the hydroxyl stretching modes that dominate the zero point energy term. Mechanical Properties There have been previous theoretical studies of aluminum hydroxide polymorphs that have attempted to empirically parametrize interatomic potentials to reproduce the structures
TABLE 4: Elastic Constants of Gibbsite Calculated Using Both ab Initio and Ionic force field Methodsa elastic constant
ab initio
force field
C11 C22 C33 C44 C55 C66 C12 C13 C15 C23 C25 C35 C46
130.9 ( 2.3 144.7 ( 6.6 120.0 ( 3.8 23.0 ( 2.2 28.7 ( 2.4 55.1 ( 0.7 56.0 ( 1.7 4.0 ( 2.5 -4.6 ( 1.1 6.3 ( 2.6 -1.2 ( 0.7 -4.6 ( 0.5 -0.6 ( 0.2
122.1 141.8 76.5 10.4 13.7 34.0 71.8 3.4 -8.8 8.5 -0.6 -1.5 -3.0
a For the ab initio values, the numerical uncertainties are also quoted. All values are in GPa.
of such phases.7 Potential models are particularly important for the prediction of morphology and surface-related properties, given the interest in the growth of gibbsite, where calculations on many surfaces have to be performed.30 This would represent a computationally challenging task for ab initio methods and therefore reliable force fields are crucial in this area. Inclusion of information concerning the curvature of the energy surface about the experimental structure is of paramount importance in the derivation of good empirical potentials.31 This requires data such as the elastic constants tensor, the bulk modulus, dielectric constants, and other physical properties which are related to the second derivatives of the energy with respect to either internal or external strain. In the case of aluminum hydroxide there is a paucity of data of this nature. Recently an X-ray diffraction study of the structure of gibbsite under applied pressure has been published32 which suggests that it ultimately undergoes a pressure induced phase transition to nordstrandite. From the data the authors also managed to extract the bulk modulus as being 85 GPa, though there was considerable noise in the curve of volume vs pressure. This represents the only experimental information concerning the mechanical stiffness of aluminum hydroxide of which we are aware. To better characterize pure gibbsite as a material, we have used ab initio plane wave calculations to calculate both the bulk modulus and to predict the elastic constants tensor. The latter information will prove particularly valuable in the future refinement of empirical potentials for this class of material, as explained previously. Although it is unlikely that the elastic constants of gibbsite will be measured experimentally for the foreseeable future, due to the small size of crystals that can be obtained, there are new methods for their determination that can use much smaller samples than other approaches,33 but which rely on refining initial estimates. Hence, the present values would prove useful for providing starting values for such experiments. The symmetry-unique elements of the elastic constants tensor are reported in Table 4 along with limits on the numerical uncertainties arising from the finite difference procedure. Also given for comparison are the same values as calculated using a recent shell model interatomic potential derived for the study of gibbsite.7 Although there are some significant differences between the ab initio and potential-based values, most particularly for the value of C33, the general trends in values and the absolute magnitude of the on-diagonal terms for strain within the layers are generally quite reasonable, especially considering there was no information included in the fit concerning the mechanical properties. However, there is clearly scope for refinement of the potential in the light of these results.
10240 J. Phys. Chem. B, Vol. 105, No. 42, 2001 Perhaps the most significant aspect of the elastic constants tensor is that it reveals information concerning the anisotropy in the mechanical stiffness of gibbsite that is not discernible from the bulk modulus. As expected, the crystal is most readily deformed along the direction corresponding to C33 (excluding shearing). However, what is more surprising is the low degree of anisotropy, which suggests that the force constants for distorting the interlayer hydrogen bonding are only marginally weaker than those for perturbing the octahedra within the layers. Conversely, from considering C44, C55, and C66, the distortion of the octahedra is almost twice as stiff as the interaction between the layers with respect to shearing. The bulk modulus of gibbsite has been determined in two ways, both with the plane wave methodology. Both the fitting of the pressure vs volume curve to an equation of state and the direct calculation from the elastic constants tensor yielded the same value, within numerical precision, of 55 GPa. This compares to the experimental estimate of 85 GPa and the value from the interatomic potentials of 45 GPa. While the overestimation of the cell volume is partly responsible for the discrepancy between theory and experiment, this is insufficient to explain such a large error. Indeed, the accumulated evidence from previous calculations of bulk moduli on a variety of materials is that GGA estimated values are rarely wrong by more than a few percent, except in pathological cases, such as van der Waals solids. Certainly an underestimation by 35% would be very surprising. Given the difficulties of accurately measuring the equation of state under pressure experimentally for gibbsite, primarily due to issues of sample quality, it would be wrong to suggest that the density functional approach had failed in this case. Recently, the bulk modulus of diaspore, R-AlOOH, has been calculated using the same approach as in the present work.34 The calculated value of 148 GPa agreed well with experimental values obtained from the determination of the elastic constants, but differs from values extracted from highpressure X-ray measurements, just as in the case of gibbsite. There is also a further possible explanation for the disagreement. Because of the conditions under which gibbsite is usually prepared, it is difficult to obtain a pure sample of the material. Very commonly there will be sodium, or other cation, impurities present in the experimental sample. If such cations act to bind the layers together more strongly than the hydrogen bonds then this would raise the bulk modulus. Further calculations are in progress to determine the influence of sodium doping on the properties of gibbsite which may resolve this issue. Vibrational Spectra One of the best techniques for in situ identification of different polymorphs of aluminum hydroxide is via characterization of the vibrational modes,35 particularly using Raman spectroscopy.36 Here there are distinct differences in the hydroxyl stretching bands between gibbsite and bayerite, with the peaks for the latter being systematically shifted to higher wavenumber. As discussed in the methodology section, the vibrations of both minerals have been determined using finite differences. Currently this yields only the mode frequencies and not their intensities, thus making full comparison against experiment difficult. Furthermore, it is well-known that there are small systematic errors in vibrational frequencies associated with the particular quantum mechanical Hamiltonian and computational conditions, such as the basis set used, partially as a consequence of deviations from experimental geometries, as well as differences due to anharmonicity. Hence, it is common to introduce scaling factors to correct, to a first approximation, for some of
Gale et al. TABLE 5: Hydroxyl Stretching Frequencies (in Wavenumbers) and Their Assignments for Gibbsite and Bayeritea atom
gibbsite
bayerite
H1 H2 H3 H4 H5 H6
3569, 3570, 3575, 3578 3423, 3425, 3426, 3426 3302, 3311, 3322, 3416 3437, 3438, 3440, 3440 3272, 3274, 3285, 3297 3189, 3192, 3207, 3212
3306, 3308, 3309, 3317 3324, 3326, 3335, 3355 3387, 3390, 3397, 3441 3464, 3465, 3468, 3469 3471, 3473, 3480, 3481 3610, 3612, 3615, 3618
a Note for the assignments, the specified atom represents the dominant hydrogen contributing to the eigenvector, but not necessarily the only one.
these differences, which are typically of the order of a few percent.37 In the present work, the reported frequencies will be unscaled and consequently direct quantitative agreement is not expected, nor should it be, since we would be comparing calculated harmonic frequencies with anharmonic experimental values. The focus will therefore not be on the absolute positions, but the trends within the modes and the assignment of the spectra. There is one further complication when calculating the vibrational frequencies at the gamma point with this method, in that the LO-TO splitting that should occur will be absent. Consequently, some of the calculated frequencies will be slightly in error. However, we can readily assess the extent of the problem using the interatomic potential methodology where it is trivial to run calculations exactly at the center of the Brillouin zone and fractionally away from it. Such calculations indicate that the largest shift of any mode is less than 30 cm-1, and this only occurs for a single hydroxyl shift. Most frequencies are only perturbed by a few wavenumbers, which is of the order of the numerical uncertainty and therefore negligible. Hence we will neglect the influence of the LO-TO splitting in the following results and discussion. The most readily identified peaks in the vibrational spectrum are the hydroxyl stretching modes because they are completely separated from all other frequencies. Calculated mode frequencies for gibbsite and bayerite are given in Table 5 along with assignments to particular hydrogens based on the dominant contributors to the eigenvectors. The nomenclature for the hydrogens refers to the labeling given in Figure 2. It should be noted that while all the stretching modes have been apportioned to particular hydrogens, in many cases there is mixing between the modes of different hydrogens. For example, in the case of gibbsite H3, H4 and H5 tend to contaminate each other’s eigenvectors significantly, while for bayerite there are combinations involving H1 with H2, and H4 with H5. However, some hydrogens give rise to quite distinct modes. For gibbsite, the peaks due to H1 stretches are separated from other modes by over a 100 cm-1 at the high-frequency end of the spectrum, while those due to H6 are confined to the low-frequency extreme of the stretching region. In the case of bayerite, there is no distinct gap at the low-frequency end of the hydroxyl stretches, but there is again a clear splitting at the other extreme due to H6. When the structure of both gibbsite and bayerite are considered, the hydrogens are oriented primarily either in the plane of the layers (xy plane) or perpendicular to this (z axis). Thus, the vibrational modes can be seen to be polarized accordingly based on the dominant hydrogen contributors. What is found is that of the 24 hydroxyl stretching modes, the 12 with the highest frequencies are those that lie mainly in the xy layers. From this it is possible to infer that the hydrogen bonding between layers
Properties of Aluminum Hydroxide is greater than that within the layers, since stronger interactions of this kind lead to an increased red shift of the frequencies. This is what would be expected based on the fact that there is more freedom to achieve the optimum hydrogen bonding distance between the layers since it is only a matter of finding the best compromise between the three symmetry-distinct hydrogens oriented in this direction. Recently Wang and Johnston38 have made an experimental assignment of the hydroxyl stretching modes of gibbsite based on the polarization direction of each band observed when compared to the orientation of the hydroxyl groups in the crystal structure. They identified six bands at 3623, 3526, 3519, 3433, 3370, and 3363 cm-1 which were assigned to H1, H2, H4, H6, H3, and H5, respectively. Their results are broadly consistent with the present findings. In particular, both studies find that the three highest modes are associated with the stretching modes perpendicular to the 001 direction and that H1 has the greatest frequency with a significant separation. The order of H2 and H4 differs, though given the closeness of the modes it is probably beyond the accuracy of the calculation to make an unambiguous assignment here. Similarly the ordering within the triad of vibrations parallel to 001 differs. Once prediction of the intensities for the modes becomes routine, for both infrared and Raman spectra, then it may be feasible to determine the validity of the present assignments. The quantitative description of the gibbsite vibrational spectrum, with respect to the mode positions, is somewhat limited in comparison to the experimental data. The highest frequency band is predicted to be approximately 50 cm-1 too low. This tendency to underestimate the stretching frequencies, in contrast to most ab initio calculations that typically tend to overestimate them, is largely related to the fact that the hydroxyl bond lengths are too long. Hence a determination of the vibrational frequencies using the plane wave approach would be more accurate. More discouraging is that the range of frequencies spanned by the hydroxyl stretching modes is too great in the calculation, with the lowest mode disagreeing by the order of 170 cm-1. Again some allowance must be made for the fact that not all modes in the calculated phonon spectrum will have nonzero intensity and for the width of the experimental bands. Although the quantitative positions of the hydroxyl stretching modes are not neccessarily accurate, it is still possible to consider the shift in modes that occurs between the two phases. Experimentally a shift of approximately +36 cm-1 is found in the highest peak when going from gibbsite to bayerite, which is reproduced by the calculations with a corresponding displacement of +40 cm-1. This shift can be understood on the basis of the smaller unit cell volume of bayerite in comparison to gibbsite. Similarly the hydroxyl bending modes demonstrate the same tendency with highest frequencies of 1138 and 1152 cm-1 for gibbsite and bayerite, respectively. Conclusions The energetic and structural properties of the two dominant polymorphs of aluminum hydroxide have been determined for the first time using ab initio techniques. Both plane wave and localized basis set approaches have been employed in order to evaluate their relative performance, as well as characterizing the intrinsic uncertainties associated with the choice of computational implementation. The present results are in accord with the available experimental calorimetric data by finding that gibbsite is the thermodynamically favored material at low temperatures and ambient pressure.
J. Phys. Chem. B, Vol. 105, No. 42, 2001 10241 Structurally, both the ab initio techniques employed demonstrated the systematic deviations that would be expected: the GGA calculations tended to overestimate the unit cell volume, while using the LDA leads to a more dramatic overcompression of the cell due to shrinking of the interlayer spacing. When contrasting the plane wave and localized basis set results, the former are found to be consistently superior for the structure due to the greater degree of convergence of the basis set. This can of course be rectified by increasing the localized basis from double-ζ plus polarization to include more flexibility. However, the present basis set represents an acceptable compromise between precision and computational speed that would be tractable for calculations on larger systems. The full elastic constants tensor has been determined in the case of gibbsite by finite differences. In turn, this has been used to yield the bulk modulus for this phase that agrees well with the value determined from fitting an equation of state to the calculated variation of the unit cell volume with applied isotropic external pressure. The calculated value of 55 GPa is much lower than a recent experiment determination where it was found to be 85 GPa. It is proposed that one possible cause of the difference in values, beyond experimental and calculation uncertainties, is the presence of impurities within the sample of gibbsite. Calculations are currently under way to study the location of sodium within gibbsite and the influence it would have on the bulk modulus. Acknowledgment. J.D.G. would like to thank the Royal Society for a University Research Fellowship and the EPSRC for provision of computing facilities. A.L.R. and J.D.G. also thank the ARC for the award of an IREX grant. The financial support of the Australian Federal Government through its Cooperative Research Centres Program is also gratefully acknowledged. References and Notes (1) Schoen, R.; Roberson, C. E. Am. Miner. 1970, 55, 43. (2) Saalfeld, H.; Jarchow, O. Neues Jahrb. Mineral., Abh. 1968, 109, 185. (3) Chao, G. Y.; Baker, J.; Sabina, A. P.; Roberts, A. C. Can. Mineral. 1985, 23, 21. (4) Wefers, K.; Misra, C. Oxides and Hydroxides of Aluminium. Alcoa Technical Paper No. 19, Revised; Alcoa Laboratories: 1987. (5) Saalfeld, H.; Wedde, M. Z. Kristallogr. 1974, 139, 129. (6) Teppen, B. J.; Rasmussen, K.; Bertsch, P. M.; Miller, D. M.; Schafer, L. J. Phys. Chem. B 1998, 101, 1579. (7) Fleming, S. D.; Rohl, A. L.; Lee, M. Y.; Gale, J. D.; Parkinson, G. M. J. Cryst. Growth 2000, 209, 159. (8) Kubicki, J. D.; Apitz, S. E. Am. Miner. 1998, 83, 1054. (9) Kohn, W.; Sham, L. J. Phys. ReV. 1965, 140, 1133A. (10) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1996, 77, 3865. (11) Payne, M. C.; Teter, M. P.; Allan, D. C.; Arias, T. A.; Joannopoulos, J. D. ReV. Mod. Phys. 1992, 64, 1045. (12) Milman, V.; Winkler, B.; White, J. A.; Pickard, C. J.; Payne, M. C.; Akhmatskaya, E. V.; Nobes, R. H. Int. J. Quantum Chem. 2000, 77, 895. (13) Vanderbilt, D. Phys. ReV. B 1990, 41, 7892. (14) Monkhorst, H. J.; Pack, J. D. Phys. ReV. B 1976, 13, 5188. (15) Sanchez-Portal, D.; Ordejon, P.; Artacho, E.; Soler, J. M. Int. J. Quantum Chem., 1997, 65, 453. (16) Artacho, E.; Sanchez-Portal, D.; Ordejon, P.; Garcia, A.; Soler, J. M. Phys. Status Solidi B 1999, 215, 809. (17) Artacho, E.; Gale, J. D.; Garcia, A.; Ordejon, P.; Sanchez-Portal, D.; Soler, J. M. SIESTA, Version 1.0.30; 2000. (18) Ordejon, P. Phys. Status Solidi B 2000, 217, 335. (19) Sankey, O. F.; Niklevski, D. J. Phys. ReV. B 1989, 40, 3979. (20) Troullier, N.; Martins, J. L. Phys. ReV. B 1991, 43, 1993. (21) Giannozzi, P.; de Gironcoli, S.; Pavone, P.; Baroni, S. Phys. ReV. B 1991, 43, 7231. (22) Wentzcovitch, R. M.; Ross, N. L.; Price, G. D. Phys. Earth Planet. Inter. 1995, 90, 101.
10242 J. Phys. Chem. B, Vol. 105, No. 42, 2001 (23) Karki, B. B.; Stixrude, L.; Clark, S. J.; Warren, M. C.; Ackland, G. J.; Crain, J. Am. Mineral. 1997, 82, 51. (24) Winkler, B.; Hytha, M.; Warren, M. C.; Milman, V.; Gale, J. D.; Schreuer, J. Z. Krist., 2001, 216, 67. (25) Gale, J. D. J. Chem. Soc., Faraday Trans., 1997, 93, 629. (26) Zigan, F.; Joswig, W.; Burger, N. Z. Krist., 1978, 148, 255. (27) Winkler, B.; Milman, V.; Hennion, B.; Payne, M. C.; Lee, M. H.; Lin, J. S. Phys. Chem. Miner. 1995, 22, 461. (28) Sato, T. J. Appl. Chem. 1961, 11, 207. (29) Yoshino, M.; Takahashi, K.; Okuda, Y.; Yoshizawa, T.; Fukushima, N.; Naoki, M. J. Phys. Chem. A 1999, 103, 2775. (30) Fleming, S. D.; Parkinson, G. M.; Rohl, A. L. J. Cryst. Growth 1997, 178, 402.
Gale et al. (31) Gale, J. D. Philos. Mag. B 1996, 73, 3. (32) Huang, E.; Lin, J. F.; Xu, J.; Huang, T.; Jean, Y. C.; Sheu, H. S. Phys. Chem. Miner. 1999, 26, 576. (33) Migliori, A.; Sarrao, J. L.; Visscher, W. M.; Bell, T. M.; Lei, M.; Fisk, Z.; Leisure, R. G. Phys. B 1993, 183, 1. (34) Winkler, B.; Hytha, M.; Pickard, C. J.; Milman, V.; Warren, M. C.; Segall, M. 2001 Eur. J. Mineral. 2001, 13, 343. (35) Frost, R. L.; Kloprogge, J. T.; Russell, S. C.; Szetu, J. L. Appl. Spectrosc. 1999, 53, 423. (36) Huneke, J. T.; Cramer, R. E.; Alvarez, R.; El-Swaify, S. A. Soil Sci. Am. J. 1980, 44, 131. (37) Scott, A. P.; Radom, L. J. Phys. Chem. 1996, 100, 16502. (38) Wang, S. L.; Johnston, C. T. Am. Miner. 2000, 85, 739.
