Structure and Dynamics of AlPO4-5 and Other Aluminophosphates

Jul 3, 2007 - ... Hélioparc Pau−Pyrénées, 2, Avenue Pierre Angot, 64053 Pau Cedex 9, France, and Department of Chemistry and Biochemistry and Cen...
0 downloads 0 Views 592KB Size
10972

J. Phys. Chem. C 2007, 111, 10972-10981

Structure and Dynamics of AlPO4-5 and Other Aluminophosphates: Classical Molecular Dynamics and ab initio Calculations Patrice Bordat,† Johanna Kirstein,‡ Pierre Labe´ guerie,† Mohammadou Merawa,† and Ross Brown*,† Institut pluridisciplinaire de recherche sur l’enVironnement et les mate´ riaux, umr 5254 du C.N.R.S., l’UniVersite´ de Pau et des pays de l’Adour, He´ lioparc Pau-Pyre´ ne´ es, 2, AVenue Pierre Angot, 64053 Pau Cedex 9, France, and Department of Chemistry and Biochemistry and Center for NanoScience (CeNS), UniVersity of Munich, Butenandtstr. 5-13 (E), 81377 Munich, Germany. ReceiVed: February 1, 2007; In Final Form: May 7, 2007

While the topology of the molecular sieve AlPO4-5 (structure type AFI) is well established, details of the structure determined experimentally, such as the precise space group, have varied significantly since its first synthesis 25 years ago. Classical molecular dynamics simulations with several force fields, periodic ab initio energy minimization, and ab initio magnetic-shielding calculations are therefore applied here to models of calcined AlPO4-5 and other aluminophosphates and compared to experimental data to help to resolve experimental discrepancies. Empirical models with partial atomic charges are overall in better agreement with a range of diffraction and NMR data for AlPO4-5 and other aluminophosphate sieves, than is a model with formal cation charges and polarizable oxygens. An order-disorder transition is found in the model, which may explain disagreements in the literature. It is predicted to occur close to room temperature from a triclinic to an on average hexagonal space group. Periodic ab initio and density functional theory calculations confirm that the space group of the most stable structure is not P6cc.

1. Introduction Since their synthesis in the early 1980s, aluminum phosphate molecular sieves have been investigated with many experimental and theoretical methods. Beside their uses as molecular sieves and in heterogeneous catalysis, they are fascinating materials in their own right. A specially useful feature of porous materials synthesized by sol-gel processes, such as silicas and zeolites, is the possibility to encapsulate large organic molecules such as dyes, combining advantages of the guest compound with those of the rigid framework. Such hybrid materials have demonstrated or have potential applications as laser media, nonlinear optical components, or sensors for medical and environmental applications. In particular, AlPO4-5 has been used as a host for a number of dyes, including coumarins and DCM,1 azo dyes,2 p-nitroaniline,3,4 rhodamines,5,6 and azo-benzene.7 AlPO4-5 was recently used as a host for a series of oxazine dyes for which differing degrees of alignment with the host channels were reported, depending subtley on the molecular structure of the guest. Rotational and spectral jumps of single guest molecules were also observed by confocal microscopy.8 A detailed understanding of the structure of the composite materials and of deformations that may be induced by the guest molecule would clearly be helpful to discuss the interesting dynamics of these systems. Yet, although AlPO4-5 was the first of the aluminophosphate molecular sieves to be synthesized, its structure continues to attract attention. The material is synthesized by a sol-gel process from orthophosphoric acid and a source of alumina in the presence of a templating * To whom correspondence should be addressed. E-mail: ross.brown@ univ-pau.fr. † Institut pluridisciplinaire de recherche sur l’environnement et les mate´riaux. ‡ Univeristy of Munich.

molecule. The as-synthesized material contains both water and template molecules, which may be removed by calcination. Figure 1 shows a molecular model of AlPO4-5 (AFI structure type). The most obvious features are puckered 12-membered rings of alternating, corner-sharing AlO4 and PO4 tetrahedra. These rings are connected by oxygen bridges into sheets or gaskets in the a, b plane, forming in the process smaller, fourand six-membered rings. The sheets are connected by oxygen bridges along the c crystal axis, to produce one-dimensional channels of which the largest, formed by the 12 rings, have a van der Waals diameter of about 7.3 Å (see Figure 1). The topology of the AFI framework (cf. International Zeolite Association structure commission) has the space group P6/mcc. There are thus 6-, 3-, and 2-fold axes perpendicular to the sheets at the centers of the 12, 6, and 4 rings, and considering the bond connectivity, there is only one topologically distinct kind of each T-atom, accompanied by four distinct oxygen sites. The original structural determination on AlPO4-5 by powder X-ray diffraction on material containing tetrapropylammonium template ions, led to the space group P6cc.9 X-ray diffraction on single crystals (same template) by Qiu et al.10 confirmed this group. On the other hand, some authors have found orthorhombic, pseudohexagonal (b/a ) x3) subgroups of P6cc, attested by splitting of some of the diffraction peaks of the parent group. Thus, Ohnishi et al.11 interpreted their powder X-ray diffraction data in terms of the subgroup Ccc2. Their material underwent a reversible phase transition to the hexagonal form (P6cc) around 400 K. The template in their material was tropine with fluoride anions. Mora et al.12 and Ikeda et al.,13 applying neutron powder diffraction, found the space group P6cc for material containing triethylammonium template ions. A lower symmetry, P6, was found by Klap et al.,14 using X-ray diffraction on single crystals containing the same template.

10.1021/jp070888e CCC: $37.00 © 2007 American Chemical Society Published on Web 07/03/2007

Dynamics of AlPO4-5 and Other Aluminophosphates

Figure 1. Connectivity of AlPO4-5. (a) Orthogonal view down the c axis of a sheet of 12 rings, showing alternation of Al and P and (b) view down a, showing intersheet oxygen bridges connecting sheets of 12 rings. Al and P, large dark and light gray disks; O, small dark disks. Coordinates from Ikeda et al.13 (Pnn2 space group).

An interesting feature of all the refinements in the group P6cc is the large, anisotropic displacement parameter of the oxygens in the intersheet Al-O-P bridges. It has long been interpreted as evidence of disorder of the oxygen,9,13 but could also result from an intimate mixture of three domains of symmetry P6.14 A further difficulty of structural determinations assuming the space group P6cc is that the same oxygen atoms form nearly linear Al-O-P bond angles, far outside the usual range of 140160°. These angles were closer to normal when the structural refinement was performed in subgroups of lower symmetry. The calcined material was examined by Mora et al.12 using high-resolution X-ray and neutron powder diffraction. They reported splitting of some of the parent, P6cc diffraction peaks and found that their samples were best described by an orthorhombic subgroup, Pcc2. Richardson et al.,15 applying neutron powder diffraction, found it necessary to assume disorder of Al and P (symmetryP6/mcc) to achieve convincing bond lengths in material calcined after synthesis with tetrapropylammonium templates. However, subsequent NMR studies by Peeters et al. showed that Al and P alternate strictly in the framework.16 These authors applied powder X-ray diffraction (no space group was reported) and 27Al and 31P magic-angle spinning and double-rotation NMR. Their conclusion was that there are (at least) three independent T-sites, revealed to varying degrees in different samples, by shoulders or splitting of the NMR lines. This observation is incompatible with the space group P6cc, which requires there be but one T-site. The assignments Pnn2, suggested by Ikeda et al.13 for calcined samples, or Pcc2,12 both with 6 T-sites, might be compatible with the NMR data, depending on the sensitivity of the NMR spectra to deformations

J. Phys. Chem. C, Vol. 111, No. 29, 2007 10973 of the lattice with respect to P6cc. The presence of three T-sites is also consistent with the assignment Ccc2 by Ohnishi et al.11 for the as-synthesized material and with the description of the as-synthesized material as a mixture of domains with P6 symmetry.14 On the other hand, Liu et al.17 examined both the assynthesized material (tetraethylammonium template ions) and carefully annealed, calcined samples, by powder X-ray diffraction and electron diffraction. They reported no evidence of splitting of the parent X-ray diffraction peaks, but did find highly structured diffuse electron diffraction. Applying lattice dynamical calculations, they interpreted the diffuse diffraction as due to the thermal excitation of many rigid unit modes of the material. The local, instantaneous symmetry may then differ from P6cc, allowing for normal angles in the intersheet oxygen bridges, while maintaining an average P6cc symmetry. The discrepancies about the space group, the phase transition reported by Ohnishi et al., and the apparent disorder of one of the oxygens suggest a degree of flexibility of the AFI framework. Indeed, condensation of rigid unit modes was suggested to explain the incommensurate modulation along the c axis of the all-silica analogue, SSZ-24.18 AlPO4-5 is thus an appealing candidate for theoretical work, which might provide some new insights into the experimental discrepancies. Molecular dynamics simulations could be particularly useful to capture features of the complicated dynamics of AlPO4-5 inaccessible to simple energy minimizations. Several papers have reported ab initio studies of different aluminophosphates. For example, the structure and dynamics of water in VPI-5 were studied in cluster models19 and CarParrinello molecular dynamics simulations.20 Hydrated AlPO434 was studied recently with ab initio molecular dynamics methods. Unfortunately, AlPO4-5 has a rather large unit cell, so it has not yet been widely studied with quantum chemical methods. It was one of 12 aluminophosphate structures in a study with periodic Hartree-Fock methods, of the correlation between structural parameters and the charge distribution of the oxygen atoms, but determined only at the assumed experimental geometry.21 Surprisingly, although there are several molecular dynamics studies of the diffusion of guest molecules in AlPO4-5,22-24 particularly with respect to single-file diffusion at high loadings in one-dimensional channels,25 force field models of the host itself, and indeed of aluminophosphate sieves in general, are sparse and largely restricted to energy minimizations.26-28 As we shall see below, the flexibility of the framework may cause structures at finite temperature to differ from those determined by energy minimizations. AlPO4-5 was one of several aluminophosphates examined by de Vos Burchart et al.,26 who performed energy minimizations with a consistent molecular mechanics force field. Starting from AlPO4-5 structures with space group P6cc, they obtained two minima, one with space group P6 and the other (slightly higher in energy) with the trigonal group P31c. Ruiz-Salvador et al.27 performed energy minimizations with different force fields. They found that the force field of Gale and Henson29 (with polarizable oxygens) yielded a minimum with space group P6 with the plane-bridging oxygen bond angles in the range 143-150°. They also found a very closely related orthorhombic structure (Pcc2) only 0.12 kcal/mol higher in energy. On the other hand, optimization with the force field of van Beest et al.30 resulted in a structure with successive 12-ring sheets laterally shifted with respect to one another, as was already observed by de Man et al.,31 and is inconsistent with

10974 J. Phys. Chem. C, Vol. 111, No. 29, 2007 all the experimental results at room temperature or above. These results once again point to a flexible structure, possibly with many close energy minima on the potential energy surface. The molecular dynamics simulations presented below confirm this point. The only dynamical simulation study of aluminophosphate sieve frameworks appears to be that of Demontis et al.,32 who described AlPO4-5 with a harmonic force field model. They could not decide between the groups P6cc and Pcc2, but did find that the intersheet oxygens are dynamically disordered. While the average atomic positions, such as would result from an X-ray diffraction experiment, defined linear oxygen bridges, the instantaneous bond angles at the oxygens were of the order of 160°. Disorder is thus dynamical in their models. This observation recurs in studies of natural zeolites, both for the framework33 and for structural water.34 The possibility of disorder of one of the oxygens and the report of a reversible phase transition for AlPO4-511 led us to undertake molecular dynamics simulations of this material. As a first step, we discuss here molecular dynamics simulations on the empty framework for comparison with results on calcined samples. After describing the molecular dynamics techniques in Section 2, we present the results in Section 3. With some of the results being in good agreement with experiment while others are at variance with experimental findings, we attempted to check the plausibility of the simulations in several ways. First, we performed periodic ab initio calculations on the bare framework of AlPO4-5, assuming different space groups. Second, we calculated the 27Al NMR spectrum of our model and compared it with the data of Peeters et al.16 Finally, to judge the general reliability of the force fields used, we also performed successful dynamical simulations on the sieves AlPO4-34 and AlPO4-53 (B and C phases) and on the quartz analogue, berlinite. Section 4 summarizes the conclusions, which despite some surprising features of the model structures, agree with a variety of experimental data and point to an order-disorder transition near room temperature, from a triclinic deformation of the P6cc group as the low-temperature phase, to an on average hexagonal high-temperature phase. 2. Methods 2.1. Molecular Dynamics. Molecular dynamics calculations were applied with three empirical force fields, one with formal charges on all atoms, polarizable oxygen atoms, and 6-exp functions to represent the covalent part of the interaction, derived by Gale and Henson (GH),29 the second with partial charges and a Morse potential energy function, due to Kitao et al. (MSQ),35 and the third (BKS) with partial charges and 6-exp interactions, by van Beest, Kramer, and van Santen.30 Force fields were used as specified in the literature, except for the MS-Q model which specifies charge equilibration in response to the changes in the instantaneous bonding of the centers,36 a feature not taken into account in the DL_POLY37 molecular dynamics suite used for all the simulations below. However, when this equilibration procedure was applied by Kitao et al. to another aluminophosphate sieve, VPI-5, the range of charges for crystallographically different atoms of a given type was only 0.03 electrons. We therefore used fixed average values from ref 35, qAl ) +1.4078, qP ) +0.5678, and qO ) -0.4939. The cation charges in MS-Q are thus in reverse order to the formal charges and to the Mulliken charge analysis of the DFT results on AlPO4-5 (P6cc, see below), which yields Al ) +1.903, P ) +2.224, and O )-1.032. The charges in the BKS model are qAl ) +1.4, qP ) + 3.4, qO ) -1.2. These data underline

Bordat et al. that it is the overall performance of the combination of the partial charges and the Morse or 6-exp potential energy well that is to be judged, rather than the component interactions individually. Energy minimizations unconstrained with respect to symmetry were performed as a preliminary, after which constant temperature and stress simulations were performed with the Berendsen thermostats and barostats.38 Ewald summation was used to evaluate the long range forces. Short-range terms were cut off at 12 Å (10 Å for the GH model). Periodic boundary conditions were applied. A time step of 0.5 fs was used throughout. Unless otherwise stated, results for AlPO4-5 refer to models with 2 × 2 × 3 unit cells in the a, b, and c directions for both the hexagonal and the orthorhombic starting configurations. No significant differences were observed in larger simulation boxes, including extension along the c axis to six unit cells. Symmetryconstrained energy minimizations and calculations of the phonon spectra were performed with gulp 1.3.39 2.2. Ab initio Calculations. Symmetry-constrained, periodic ab initio and density fucntional theory (DFT) energy minimizations were performed with a development version of the crystal03 code,40 in which orbitals are represented as linear combinations of Bloch functions. Each function is built from contracted Gaussian-type local atomic orbitals (AO). Here, allelectron basis sets were used for Al, P, and O atoms:41,42 Al 8-511(1), P 8-51(1), and O 6-31(1), where the first figure refers to an s shell, the others to sp shells, and d shell contractions are given in parentheses. The exponents of the outer sp and d shells were reoptimized at the Hartree-Fock (HF) level to the following values (in bohr-2 units) and used as such for the B3LYP Hamiltonian: (Al, sp) ) 0.576, 0.280, (Al, d) ) 0.470; (P, sp) ) 0.135, (P, d) ) 0.746; (O, sp) ) 0.274, (O, d) ) 0.600. Standard values for the computational tolerances were adopted for all steps of the calculation, except for the vibrational frequencies, that require more severe conditions. In the geometry optimization, a structural relaxation procedure consisting of two independent steps was performed iteratively. In the first step, the cell parameters (a, c) were optimized with the atoms at fixed fractional positions. Cell optimization was carried out by means of a modified Polak-Ribiere algorithm in which the energy gradients were evaluated numerically by means of the centraldifference formula.43 In the second step, atomic positions (x, y, z fractional coordinates) were fully relaxed at fixed cell parameters. Forces on atoms were obtained by using the analytical HF44 and DFT45 energy gradients and were used to relax the atoms to equilibrium by using a modified conjugate gradient algorithm proposed by Schlegel.46 Convergence was tested on the root mean square and the absolute value of the largest component of the gradients and the estimated displacements. The thresholds for the maximum force, the RMS force, the maximum atomic displacement, and the RMS atomic displacement on all atoms were set to 0.00045, 0.00030, 0.00180, and 0.00120 a.u. in that order. Optimization of the atomic positions was considered complete when these four conditions were satisfied. The crystal symmetry was maintained during the optimization process. The two-step structure optimization process was repeated until both the cell parameters and the atomic positions convergence criteria were satisfied.47 DFT calculations were performed with the B3LYP hybrid functional.48,49 We refer to a recent paper50 for a complete discussion of the computational conditions and other numerical aspects concerning the calculation of the phonon frequencies at the Γ ) (0,0,0) point. In the case of ionic compounds, long-range Coulomb

Dynamics of AlPO4-5 and Other Aluminophosphates

Figure 2. Skeleton view of AlPO4-5 projected onto the a, b plane: (a) the Pcc2 data of Mora et al.;12 (b) the GH model minimized from the Pnn2 data of Ikeda et al.; (c) same with the MS-Q model.

effects due to coherent displacement of the crystal nuclei are neglected in the above calculation, as a consequence of imposing the periodic boundary conditions.51 The mass-weighted Hessian needs to be corrected to obtain the longitudinal optical modes. 52 The electronic dielectric and Born charge tensors have therefore been calculated with the two Hamiltonians. The grid used for the numerical integration of the exchange-correlation term contained 75 radial points (the default value is 55) and 974 angular points (default is 434) in a Lebedev scheme. The condition for the self-consistent field (SCF) convergence was set to 10-10 hartree. The shrinking factor of the reciprocal space net was set to 8, corresponding to 75 reciprocal space points at which the Hamiltonian matrix was diagonalized. The total energies obtained with this mesh were fully converged. NMR chemical shifts were calculated with the gauge independent atomic orbital method in Gaussian 9853 and were applied to clusters cut out of the bulk materials with dangling bonds saturated with optimized hydrogens (mode redundant algorithm to preserve the geometry of the cut-out clusters). The SCF convergence was set to 10-10 hartree with numerical integration on an ultrafine grid. 3. Results and Discussion 3.1. AlPO4-5. 3.1.1. Energy Minimizations. A model of 2 × 2 × 3 unit cells of the calcined structure refined in the Pnn2 space group by Ikeda et al.13 was set up and first optimized

J. Phys. Chem. C, Vol. 111, No. 29, 2007 10975

Figure 3. Orientation of the interplane Al-O-P bridges of AlPO4-5 in the energy minima of the GH model (2 × 2 × 3 cells) viewed down the c axis. (a) Matchsticks are drawn on O’s representing the sagitta from the line of centers of Al and P to the O atom in each bridge; (b) projections of all the sagittae in the lowest energy structure found; (c) one of the many defect structures (this one less than 1 meV per AlPO4 unit above the lowest minimum), showing partial disorder of the interplane bridges.

with the force field of GH,29 by coupling to a cold heat bath in the isothermal, isostress (1 bar) ensemble with the Berendsen thermostat and barostat. No symmetry constraints were imposed. Figure 2 compares the overall relaxed structure to the experimental data and to the structure obtained with the MS-Q model. The GH model spontanteously yields a structure in good overall agreement with the experimental data. The unit cell is orthorhombic, with 24 distinct AlO4 tetrahedra in place of the 12 in the structure from ref 13. The final energy of this model was -268.070 eV per AlPO4 unit. The model tetrahedra are much less deformed than the experimental ones. The intersheet oxygen bridge angles range from 152 to 158°. Unconstrained minimization from the X-ray data of Mora et al. converged to a structure at first sight hexagonal (b/a exceeded x3 by only 0.005). But closer examination revealed an orthorhombic unit cell with a doubling of the a and b unit cell parameters (a ) 2 × 13.80 Å, b ) 2 × 23.97 Å, c ) 8.41 Å, R ) β ) γ ) 90 ( 0.001°, -268.0642 eV per AlPO4 unit) and 24 distinct tetrahedra. The intersheet oxygen angles were in the range 156 ( 2°. In fact, simulated annealing between room temperature and 0 K with the GH model quickly revealed the presence of a number of secondary minima. Over 10 were found lying within 0.005 eV per AlPO4 unit of the lowest structure (Pnn2 above). They are best characterized by the orientation of the “elbow” formed by the intersheet Al-O-P bend. In all the

10976 J. Phys. Chem. C, Vol. 111, No. 29, 2007

Bordat et al.

Figure 5. Distributions of instantaneous interplane oxygen bridge angles in molecular dynamics simulations of AlPO4-5 at 300 K, with the GH (solid line) and the MS-Q model (dotted). The thermally averaged bridge angles are nonlinear but not their instantaneous values.

Figure 4. Orthographic projections of all the AlO4 tetrahedra viewed in a standard orientation, down the intersheet Al-O bond (close to the c axis) with Al at the origin, the 12 ring on the left (broad gray arc) and the Al-O pointing out of the ring, Oout, on the right. (a) Superposition of the minimized GH model (circles) and the BKS model (crosses); (b) the same at 300K; gray cloud of points shows instantaneous positions of the oxygens relative to Al in the BKS model; circles and crosses are average positions in the GH and the BKS models.

minima, all the Al-O-P bends are close to 156° and the bridge elbows point tangentially round the periphery of the 12 rings. The lowest energy structure was found with all the elbows in alternate rows of 12 rings parallel to a pointing alternately clockwise and anticlockwise (Figure 3). In other minima, elbows point clockwise or anticlockwise in different rings in the same row or in the upper and lower halves of the same ring. Clearly, in this model, the bridges are rather free to swing about the c axis. Their precise orientation in the real material might be determined by the kind and amount of occluded template molecules and the amount of adsorbed water, which would help to explain the disagreement about the precise structure of the material. Contrary to the GH model, simulated annealing with the MS-Q and the BKS models, both with partial charges, yielded in each case a single minimum, which was obtained in three settings, related by rotations of 120° about the c axis. The minima in the two models were similar to each other and to that reported in the literature27,31 (see Figure 2c). The minima differ qualitatively from those yielded by the GH model, by the staggered arrangement of adjacent sheets in the unit cell. In this case, it is the intersheet-bridging oxygens that lie above each other along the c direction, rather than the aluminums (or

Figure 6. Sample of the simulated synchronized swinging of the AlO-P intersheet bridges at 300 K in the MS-Q model of AlPO4-5. Points: azimuths of the intersheet Al(-O-)P vectors projected onto the a, b crystal plane; only half shown for clarity, the other half being on average offset by 180°. Continuous line: fluctuations of the ratio b/a of the cell sides about x3.

phosphorus) in the GH model. Figure 2 illustrates the minimum obtained for the MS-Q model. The minimized structures, though close to hexagonal, are definitely triclinic (e.g., a ) 14.097, b ) 14.138, c ) 8.389 Å; R ) 89.58, β ) 90.40, γ ) 121.65 ( .01° for the BKS model) starting from the P6cc structure of Ikeda et al.13 There are 12 distinct AlO4 tetrahedra. This result was independent of the system size up to 4 × 4 × 6 unit cells in the a, b, and c directions (as was the dynamics, see Figure 7 below). 3.1.2. Molecular Dynamics. We also performed molecular dynamics runs at finite temperatures between 50 and 500 K. Simulations with the GH model, which were started with an orthorhombic box with 2 × 2 × 3 Pnn2 cells, converged quickly at all temperatures to average P6cc symmetry, as shown by b/a ) x3 and the presence of only one T-site and a 6-fold axis, illustrated in Figure 4 by the projected AlO4 tetrahedra at room temperature. Figure 5 shows the distribution of the bridge angles, average 159.5°, standard deviation 8°. Note that although the instantaneous intersheet Al-O-P bridges are not linear, their swinging around the Al-P (c) axis leads to an average structure with linear bridges. The time-averaged structure has symmetry P6cc. Furthermore, although the intersheet bridges swing around the c axis, they have no particular instantaneous orientation. The barriers to rotation of the Al-O-P elbows about c must be very low in the GH model, because this picture held for all simulations down to the lowest attempted here, 50 K.

Dynamics of AlPO4-5 and Other Aluminophosphates

J. Phys. Chem. C, Vol. 111, No. 29, 2007 10977

Figure 9. NMR spectrum of 27Al in AlPO4-5. Points, DOR spectrum redrawn from ref 16; dotted line, present work, AlO4 tetrahedra obtained from empirical BKS model, shifts relative to Al(NO3)3; solid line, shifted 7 ppm for comparison of the width with the experimental data.

Figure 7. (a) Sample fluctuations of the ratio b/a about x3 (straight lines) in isothermal, isostress (1 atm) simulations of models of 2 × 2 × 3 cells of AlPO4-5 with the MS-Q forcefield. Curves at 300 and 500 K are vertically offset for clarity by 0.05 and 0.1; a simulation at 300 K with 4 × 2 × 6 cells is also shown. (b) Distributions of the ratio b/a at increasing temperatures, vertically offset for clarity, sampled over time in molecular dynamics runs of at least 100 ps (except 40 ps at 220 K). The arrow marks x3.

Figure 8. Comparison of the average atomic positions in the simulated AlO4 tetrahedra at room temperature (BKS model, crosses) and those in the P6cc structural refinement of Klap et al.14 (circles). Same projection as in Figure 4.

The behavior of the MS-Q model at finite temperatures was more complex. The initial ratio b/a of the minimized structure is 1.755, significantly over x3. On heating to temperatures up to 500 K, the ratio was observed to jump randomly bewteen two temperature-dependent values, x3 - (T), and x3 + (T). Figure 7a illustrates this behavior at 300 K. The lower value was much more probable at low temperatures, but as the temperature was raised the upper and lower values became equiprobable and eventually (T) vanished as the ratio fluctuated randomly around x3 (see Figure 7b). The jumps in b/a coincide with synchronous 120° swings of the cant of all the intersheet

Figure 10. Skeleton view down the 3-fold axis of the averaged structure of AlPO4-34 at room temperature (MS-Q force field, scaled to the same unit cell as the experimental data) superposed on the experimental structure of Tuel et al.60,61

Al-O-P bridges, passing through the azimuths (60° ( 120° when b/a ) x3 (see Figure 6). In between swings, the bridges oscillate in one of the three settings found above in the simulated annealing. We checked that these conclusions did not depend on the thermostat by comparing with simulations using the Hoover-Nose´ thermostat. Results were very similar. The areas under the peaks of the histograms in Figure 7b are not entirely significant, being subject to sampling errors dependent on the length of the simulation. What is significant is the coalescence of the peaks as the temperature is raised. In fact, the model undergoes an order-disorder transition to P6cc between 400 and 500 K, as judged by the loss of coherent orientation of the intersheet bridges. The same behavior was noted in a larger simulation of 4 × 2 × 6 Pnn2 unit cells at 300 K (see Figure 7a). Throughout all the simulations below the transition, there were several distinct T-sites (see Figure 4b). Disorder close to the transition has both dynamic and static characteristics: While the simulation is in any one of the three settings obtained by energy minimization, the instantaneous intersheet bridges are significantly nonlinear (see the spread of values in Figure 5), but the bridge angles from thermally averaged coordinates (available in the Supporting Information) are linear to within 5°. The bridges in the simulation thus display dynamic disorder while being on average linear. But if one accepts the model, one would expect the real material, which is much bigger than our simulations, to be composed of finite domains with the three settings, presumably with much longer switching times than those observed here, which would amount to static disorder. 3.1.3. AlO4 Tetrahedra. Although the force fields with partial atomic charges (BKS, MS-Q) cause shearing between successive sheets of 12 rings, the details of the structure, as judged by the

10978 J. Phys. Chem. C, Vol. 111, No. 29, 2007

Bordat et al.

TABLE 1: Experimental and Simulated Cell Parameters of Aluminophosphates material or model

temperature (K)

space group

298 0 300 0 300

P3221

rt 0 300 0 0 300

Pnn2 Pnn2 P6cc P1 P1

290 0 290 0 300

R3

a (Å)

b (Å)

c (Å)

R (°)

β (°)

γ (°)

energy/ AlPO4 unit (eV)

berlinite expa GH MS-Q

expb GH MS-Q BKS

4.9438 4.925 4.925 5.123 5.14 13.7711 13.89 13.825 13.986 14.097 14.3

10.950 10.93 10.96 11.311 11.35 AlPO4-5 23.846 23.77 23.955 24.505 14.138 14.3

8.3769 8.41 8.441 8.489 8.389 8.53

90 90 90 90 90 90.0 90.0 90.0 89.59 89.58 89.6

120. 120. 120. 120. 120.

-268.192 -29.672

-268.07 89.76 90.40 90.3

88.86 121.65 121.2

-29.322 -117.995

AlPO4-34 expc,d GH MS-Q

expe GH MS-Q

expe MS-Q MS-Q a

9.414 9.326 9.325 9.525 9.53

R3

rt 0 300 0 300

Pbca

rt 0 300

C121

94.61 94.84 94.85 94.45 94.4

18.024 17.965 17.99 18.355 18.36

AlPO4-53(B) 13.917 9.655 13.90 9.68 13.89 9.67 14.159 9.825 14.19 9.84

90 90 90 90 90

16.444 16.24 16.45

AlPO4-53(C) 5.108 13.485 5.27 13.95 5.30 14.06

90 90. 90.

-267.989 -29.1044

-268.065 -29.3274

88.259 88.0 88.2

90 90. 90.

-29.365

Ref 63. b Ref 13. c Ref 60. d Ref 61. e Ref 62.

shape and size of the TO4 tetrahedra, are actually quite close to those yielded by the model with formal charges (GH). Here, we discuss the data on the AlO4 tetrahedra, conclusions from the data on the PO4 units (not shown) being similar. Thus, the tetrahedra are very similar in all three models, both at the minimum and at 300 K (Figure 4a,b). The tetrahedra are quite regular at the minima of the models, but are skewed on average at 300 K. They compare well with the experimental data seen in Figure 8, where the AlO4 tetrahedra of the BKS model at 300 K are overlaid on data from the P6cc structural refinement of Klap et al.14 Tetrahedra from other experimental data, such as the Ccc2 refinement of Mora et al.,12 are also very similar (data not shown). Differences between tetrahedra within a given model are larger for BKS and MS-Q than for GH and remain larger at room temperature. Summarizing Subsection 3.1, simulations of AlPO4-5 with the MS-Q and the BKS force fields suggest a phase transition may exist close to room temperature in agreement with the observations of Ohnishi et al.11 Symmetry lower than P6cc allows for several T-sites in the low-temperature phase, which could explain the observation of several lines in the NMR data of Peeters et al.16 Yet, while both these models are locally very similar to the GH model, the predicted shearing between alternate a, b plane gaskets in the structure is hard to reconcile with available structural data. To put the results on AlPO4-5 with these force fields in perspective, we therefore (i) undertook periodic ab initio energy minimizations assuming different constrained space groups, as a (we hope temporary) replacement for an experimental structure determination at low temperature; (ii) calculated ab initio estimations of the chemical shifts of 27Al from the BKS optimized structure for comparison with the NMR spectra of Peeters et al.; and (iii) simulated two other

porous aluminophosphates and R-berlinite (analagous to R-quartz), using the MS-Q and the GH models. 3.2. Periodic ab initio Calculations on AlPO4-5. Energy optimizations of berlinite and of AlPO4-5 were performed using an optimized all electron basis (see Section 2.2 for details). Calculations on berlinite converged readily to a minimum in the experimental space group in good agreement with experimental data (see Table 1). Full details of the calculated elastic constants and phonon spectra will be published elsewhere.54 Optimizations of AlPO4-5 were also attempted, assuming different space groups reported in experimental papers. Calculations were difficult because of the large numbers of atoms in the unit cells. In fact, although optimization was attempted for example from the minima with distorted (b/a ) x3 () structures from MD minimizations, satisfactory convergence to a stationary point was achieved only for the P6cc space group. Agreement with experimental data was then fair, but calculation of the phonon frequencies at the Γ point revealed several imaginary values. This shows that the P6cc structure is a saddle, not a minimum. The result is consistent with the molecular dynamics simulations and experimental data (e.g., ref 17), which tend to show that the P6cc group at room temperature is a dynamical average structure. 3.3. NMR Spectra. Peeters et al.16 reported at least three lines or shoulders in the NMR spectra of 27Al in some of their AlPO4-5 samples, which are inconsistent with the P6cc space group. As a further check on the plausibility of the unexpected structure obtained with the MS-Q and BKS empirical models, chemical shifts of 27Al were evaluated by ab initio methods. Magnetic-shielding constants (especially that of aluminum) were calculated using the gauge including atomic orbitals approach55-57 as implemented in Gaussian 98.53 In preliminary

Dynamics of AlPO4-5 and Other Aluminophosphates

J. Phys. Chem. C, Vol. 111, No. 29, 2007 10979

Figure 11. Simulated phase transition (MS-Q model) of AlPO4-53B to AlPO4-53C: (a) experimental structure of AlPO4-53B; (b) end of simulation at 600 K; (c-e) after 50, 75, and 150 ps of simulation at 750 K; (f) experimental structure of AlPO4-53C. All views down the c crystal axis of AlPO4-53B.

calculations, AlO4 tetrahedra were extracted from the BKS model of AlPO4-5, accurately optimized with the gulp code.39 Hydrogens were added to saturate the valence of the oxygens and were optimized while retaining the structure of each tetrahedron. A comparison of Hamiltonians, HF, B3LYP, B3PW91, and of basis sets up to 6-31+G** showed that the relative chemical shifts were insensitive to details of the calculations. The three aluminums with smallest, largest, and intermediate chemical shifts in the preliminary calculations were then studied in more detail in calculations on clusters containing up to four phosphorus atoms in the second coordination shell around the central Al atom (sample clusters available in the Supporting Information). They were extracted from the BKS model and were saturated with hydrogens to investigate the influence of second neighbors. These calculations58 showed that (i) the chemical shifts increased from ca. 500 to 540 ppm on passing from clusters with 0 to 4 phosphorus neighbors; (ii) chemical shifts converged to stable values with 3-4 P atoms in the second coordination shell of Al; and (iii) chemical shifts calculated with four P neighbors were linearly related to those found in the preliminary calculations without phosphorus. These results were then applied to the estimation of the full chemical shifts of all the tetrahedra, including the influence of phosphorus second neighbors, from the calculations of the AlO4H4 clusters. Final chemical shifts are reported below relative to Al(NO3)3, a widely used reference solution,59 for which we find a computed magnetic shielding of 562.5471 ppm. The 27Al spectrum shown in Figure 9, obtained by summing Lorentzian lines at the calculated shifts (half-width at half-maximum 10 ppm), is in good agreement with the spectrum reported by Peeters et al.,16 tending to enhance our confidence in the structure of AlPO4-5 predicted with the empirical models with partial atomic charges. While computational cost prevented thermal averaging of the NMR spectra by using instantaneous tetrahedra from the molecular dynamics simulations, the agreement with experiment is at least qualitatively consistent with the results of Peeters et al.16

3.4. Other Aluminophosphates. 3.4.1. Berlinite. The GH force field for aluminophosphates, which was derived with reference to experimental data on berlinite, performs well on this material, and our calculations reproduce in all respects the full data in ref 29. The MS-Q force field performs adequately, maintaining the hexagonal symmetry, but overestimates cell sides by some 3.5%. Cell sides of all systems examined were overestimated with the MS-Q model and underestimated with the GH model (see Table 1). 3.4.2. AlPO4-34. The structures of calcined and rehydrated AlPO4-34 have been determined by Tuel et al.60,61 Simulations were performed with the GH and the MS-Q force fields on models of 3 × 3 × 3 rhombohedric unit cells (972 atoms). The rhombohedric form of the calcined material is well preserved by both the MS-Q and the GH force fields (see Table 1). The relative positions of contents of the unit cell are also well reproduced by both the GH and the MS-Q force fields. Figure 10 shows for example an overlay of the experimental and the averaged simulated structure (MS-Q, 300 K) after rescaling to the same cell side as the experimental data. 3.4.3. AlPO4-53B(C). Kirchner et al.62 reported the structures of as synthesized AlPO4-53(A) containing methylamine, of the dehydrated material, AlPO4-53(B), and of a new phase obtained by heating to 700 °C. The thermal transformation from form B to form C involves bond breaking and making in a first-order transition. Both the GH and the MS-Q models described form B adequately in simulations of 2 × 3 × 4 and 2 × 3 × 6 unit cells respectively (see Table 1). Atomic fractional coordinates were also well reproduced (data not shown). On performing a series of NσT simulations, a spontaneous transition to a new bonding topology equivalent to form C was observed at 750 K with the MS-Q model. Energy minimization from the experimental data on form C yielded a slightly lower energy than form B, consistent with the stability of form C on cooling to room temperature. On heating back to room temperature, this model accurately reproduced the experimental data. Figure 11 compares the simulated transition with the experimental results.

10980 J. Phys. Chem. C, Vol. 111, No. 29, 2007 On the other hand, no transition was observed in simulations up to 2000 K with the GH model. 4. Conclusion The present simulations were undertaken to better understand the discrepancies in the literature on the structure of AlPO4-5. Three empirical force field models were compared, two with partial atomic charges and one with formal cation charges and polarizable oxygen. Overall, models with partial charges agree better with a variety of data in the literature on AlPO4-5. They accurately reproduce other aluminophosphate materials. The main difference of the partial charge models of AlPO4-5 compared to the formal charge model is shearing of alternate sheets of 12 rings in one of three directions at angles of 120° in energy-minimized structures. The shear is averaged out in an order-disorder transition occurring around 400 K in the simulations. If such a transition really exists, it may explain the discrepancies between different reported space groups and the sensitivity of the material to the precise conditions of preparation (e.g., not all samples showed splitting of the 27Al NMR spectrum in the work of Peeters et al.16). Periodic ab initio calculations confirm that the P6cc space group is a dynamical average structure, because constrained minimization in that group leads to a saddle point, not a minimum. While most Al-O-P angles in AlPO4-5 are within the usual range for other aluminophosphates, all the models employed here lead to nearly linear thermally averaged intersheet bridges and two of them provide instantaneous values close to 180°. While surprising, this behavior may indicate a degree of frustration of interactions, subject to longer range constraints induced by the overall topology of the aluminophosphate network. The resulting strain may explain the sensitivity to occluded matter such as water or templates and the spectral and orientational dynamics of single probe dye molecules in AlPO4-5. Acknowledgment. We acknowledge helpful discussions with Pr. Henk van Koningsveld (Delft) and thank Dr. Florence Porcher (Nancy) for providing unpublished results and careful reading of the manuscript. Part of the calculations were performed on the machines of the M3PEC consortium, Universite´ de Bordeaux I. Supporting Information Available: Fractional and cartesian coordinates of models minimized with the BKS force field and sample clusters used for the NMR shielding calculations. This material is available free of charge via the Internet at http:// pubs.acs.org. References and Notes (1) O ¨ . Weiβ, J.; Loerke, U. W.; Marlow, F.; Schu¨th, F. J. Solid State Chem. 2002, 167, 302-309. (2) Braun, I.; Schulz-Ekloff, G.; Bosckstette, M.; Wo¨hrle, D. Zeolites 1997, 19, 128-132. (3) Klap, G. J.; van Klooster, S. M.; Wu¨bbenhorst, M.; Jansen, J. C.; van Bekkum, H.; van Turnhout, J. J. Phys. Chem. B 1998, 102, 95189524. (4) Wu¨bbenhorst, M.; Klap, G. J.; Jansen, J. C.; van Bekkum, H.; van Turnhout, J. J. Chem. Phys. 1999, 111, 5673-5640. (5) Bockstette, M.; Wo¨hrle, D.; Braun, I.; Shultz-Ekloff, G. Microporous Mesoporous Mater. 1998, 23, 83-96. (6) Ghanadzadeh, A.; Zanjanchi, M. A. Spectrochim. Acta., Part A 2001, 57, 18665-1871. (7) Hoffmann, K.; Resch-Genger, U.; Marlow, F. Microporous Mesoporous Mater. 2000, 41, 99-106. (8) Seebacher, C.; Hellriegel, C.; Bra¨uchle, C.; Ganschow, M.; Wo¨hrle, D. J. Phys. Chem. B. 2003, 107, 5445-5452.

Bordat et al. (9) Bennett, J. M.; Cohen, J. P.; Flanigen, E. M.; Pluth, J. J.; Smith, J. V. Vol. 218 of ACS Symp. Ser., page 109, American Chemical Society, 1155 16th Street, NW Washington DC 20036, 1983. (10) Qiu, S.; Pang, W.; Kessler, H.; Guth, J.-L. Zeolites 1989, 9, 440. (11) Ohnishi, N.; Qiu, S.; Terasaki, O.; Kajitani, T.; Hiraga, K. Microporous Mater. 1993, 2 (1), 73-74. (12) Mora, A. J.; Fitch, A. N.; Cole, M.; Goyal, R.; Jones, R. H.; Jobic, H.; Carr, S. W. J. Mater. Chem. 1996, 6 (11), 1831. (13) Ikeda, T.; Miyazawa, K.; Izumi, F.; Huang, Q.; Santoro, A. J. Phys. Chem. Solids 199, 60, 1531-1535. (14) Klap, G. J.; van Koningsveld, H.; Graafsma, H.; Schreurs, A. M. M. Microporous Mesoporous Mater. 2000, 38, 403-412. (15) Richardson, J. W.; Pluth, J. J.; Smith, J. V. Acta. Crystallogr., Sect. C 1987, 43, 1469-1472. (16) Peeters, M. P. J.; van de Ven, L. J. M.; de Haan, J. W.; van Hooff, J. H. C. J. Phys. Chem. 1993, 97, 8254-8260. (17) Liu, Y.; Withers, R. L.; Lore´n, L. Microporous Mesoporous Mater. 2001, 42, 103-111. (18) Liu, Z.; Fujita, N.; Terasaki, O.; Ohsuna, T.; Hiraga, K.; Camblor, M. A.; Dı´az-Caban˜as, M.-J.; Cheetham, A. K. Chem.sEur. J. 2002, 8, 4549-4556. (19) Prasad, S.; Vetrivel, R. J. Phys. Chem. 1994, 98, 1579-1583. (20) Fois, E.; Gamba, A.; Tilocca, A. J. Phys. Chem. B 2002, 106, 48064812. (21) Larin, A. V.; Vercauteren, D. P. J. Mol. Catal., A 2001, 166, 7385. (22) Scholl, D. S.; Fichthorn, K. A. J. Chem. Phys. 1997, 107, 43844389. (23) Kantola, J.-H.; Vaara, J.; Rantola, T. T.; Jokisaari, J. J. Chem. Phys. 1997, 107, 6470-6478. (24) Scholl, D. S. Chem. Phys. Lett. 1999, 305, 269-275. (25) Gupta, V.; Nivarthi, S. S.; McCormick, A. V.; Ted Davis, H. Chem. Phys. Lett. 1995, 247, 596-600. (26) de Vos Burkchart, E.; van Bekkum, H.; van de Graff, B.; Voigt, E. T. C. J. Chem. Soc., Faraday Trans. 1992, 88 (18), 2761-2769. (27) Ruiz-Salvador, A. R.; Sastre, G.; Lewis, D. W.; Catlow, R. A. J. Mater. Chem. 1996, 6 (11), 1837-1842. (28) Demontis, P.; Mellot-Draznieks, C.; Fe´rey, G. Curr. Opin. Solid State Mater. Sci. 2003, 7, 13-19. (29) Gale, J. D.; Henson, N. J. Chem. Soc., Faraday Trans. 1994, 90, 3175. (30) van Beest, B. W. H.; Kramer, G. J.; van Santen, R. A. Phys. ReV. Lett. 1990, 64, 1955. (31) de Man, A.; Jacobs, W.; Gilson, J.; van Santen, R. Zeolites 1992, 12, 826-836. (32) Demontis, P.; Gulı´n-Gonzalez, J.; Suffriti, G. B.; Tilocca, A.; de las Pozas, C. Microporous Mesoporous Mater. 2003, 5, 427-423. (33) Alberti, A. In New DeVelopments in Zeolite Science and Technology; Murakami, Y., Iijiima, A., Ward, J. W., Eds.; Kodansha Elsevier: Tokyo, 1986; pp 437-441. (34) Demontis, P.; Gulı´n-Gonzlez, J.; Suffritti, G. B. J. Phys. Chem. B. 2006, 110 (14), 7513-7518. (35) Kitao, O.; Demiralp, E.; Cagin, T.; Dasgupta, S.; Mikami, M.; Tanabe, K.; Goddard, III, W. A. Comp. Mater. Sci. 1999, 14, 135-137. (36) Rappe´, A. K.; Goddard, W. A., III. J. Phys. Chem. 1991, 95, 33583363. (37) Smith, W., Forester, T. R. The DL_POLY_2 User Manual; CCLRC, Daresbury Laboratory, Daresbury: Warrington, England, 2001. (38) Berendsen, H. J. C.; Potsma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (39) Gale, J.; Rohl, A. Mol. Simul. 2003, 29, 291. (40) Saunders, V. R., Dovesi, R., Roetti, C., Orlando, R., ZicovichWilson, C., Harrison, N., Doll, K., Civalleri, B., Bush, I., D’Arco, P., Llunell, M. Crystal 03 User Manual; Universita de Torino: Torino, Italy, 2003. (41) Catti, M.; Dovesi, R.; Pavese, A.; Saunders, R. J. Phys.: Condens. Matter 1991, 3, 4151. (42) Harrison, N. M.; Saunders, V. R. J. Phys.: Condens. Matter 1992, 4, 3873. (43) Zicovich-Wilson, C. M. Loptcg (shell procedure for numerical gradient optimisations); Universidad Autonoma del Estado de Morelos: Morelos, Mexico, 1998. (44) Doll, K.; Harrison, N. M.; Saunders, V. R. Int. J. Quantum Chem. 2001, 82, 1. (45) Doll, K. Comput. Phys. Commun. 2001, 137, 74. (46) Schlegel, H. B. J. Comput. Chem. 1982, 3, 214. (47) Civalleri, B.; D’Arco, P.; Orlando, R.; Saunders, V. R.; Dovesi, R. Chem. Phys. Lett. 2001, 348, 131. (48) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (49) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B. 1988, 37, 785. (50) Pascale, F.; Zicovich-Wilson, C. M.; Lopez Gejo, F.; Civalleri, B.; Orlando, R.; Dovesi, R. J. Comput. Chem. 2004, 25, 888. (51) Born, M.; Huang, K. Dynamical Theory of Crystal Lattices; OUP: Oxford, 1954.

Dynamics of AlPO4-5 and Other Aluminophosphates (52) Umari, P.; Pasquarello, A.; Corso, A. D. Phys. ReV. B. 2001, 63, 094305. (53) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Andres, J. L.; Head-Gordon, M.; Replogle, E. S.; Pople, J. A.; Gaussian 98, revision a.7; Gaussian, Inc.: Pittsburgh, PA, 1998. (54) Labe´guerie, P.; Merawa, M. To be published. (55) Dodds, J. L.; Mcweeny, R.; Sadlej, A. J. Mol. Phys. 1980, 41, 1419.

J. Phys. Chem. C, Vol. 111, No. 29, 2007 10981 (56) Wolinski, K.; Hilton, J. F.; Pulay, P. J. Am. Chem. Soc. 1990, 112, 8251. (57) Ditchfield, R. Mol. Phys. 1974, 27, 789. (58) Bordat, P.; Brown, R. 2006. (59) Brouwer, D. H.; Che´zeau, J.-M.; Fyfe, C. A. Microporous Mesoporous Mater. 2006, 88, 163. (60) Tuel, A.; Caldarelli, S.; Meden, A.; McCusker, L. B.; Baerlocher, C.; Ristic, A.; Rajic, N.; Mali, G.; Kaucic, V. J. Phys. Chem. B 2000, 104, 5697-5705. (61) Poulet, G.; Sautet, P.; Tuel, A. J. Phys. Chem. B 2002, 106, 85998608. (62) Kirchner, R. M.; Grosse-Kunstleve, R. W.; Pluth, J. L.; Wilson, S. T.; Broach, R. W.; Smith, J. V. Microporous Mesoporous Mater. 2000, 39, 319-332. (63) Muraoka, Y.; Kihara, K. Phys. Chem. Miner. 1997, 24, 243253.