Jorge
Arce-Molina
ab,
Ricardo
Grau-Crespo
c,
Dewi W.
Lewis
*d and
A. Rabdel
Ruiz-Salvador†
*a
aInstitute of Materials Research and Engineering (IMRE), University of Havana, Havana 10400, Cuba. E-mail: rruisal@upo.es
bDept. Metallurgy-Chemistry, Higher Mining Metallurgical Institute, Moa, Holguin 83329, Cuba
cDepartment of Chemistry, University of Reading, Whiteknights, Reading RG6 6AD, UK
dDepartment of Chemistry, University College London, 20 Gordon St., London WC1H 0AJ, UK. E-mail: d.w.lewis@ucl.ac.uk
First published on 6th June 2018
We introduce a method to allow the screening of large configurational spaces of heteroatom distributions in zeotype materials. Based on interatomic potential calculations of configurations containing up to two heteroatoms per cell, we parameterize an atomistic effective Hamiltonian to describe the energy of multiple substitutions, with consideration of both short- and long-range interactions. Then, the effective Hamiltonian is used to explore the full configurational space at other compositions, allowing the identification of the most stable structures for further analysis. We illustrate our approach with the aluminogermanate PKU-9, where we show that increasing the aluminium concentration changes the likely siting of Al, in agreement with experiment.
The incorporation of Al in germanate frameworks requires the presence of charge-compensating extra-framework cations which will impart ion-exchange and catalytic properties. Moreover, the presence of Al enhances the stability of the framework upon template removal.14 Both the amount of Al incorporated, and its location in the framework, impact the physical and chemical properties of the resulting material, similar to what happens in aluminosilicate zeolites.15–19 A large body of experimental and computational work exists aimed at identifying and attempting to explain the distribution of Al in aluminosilicate materials.20–31 In aluminogermanates the scenario is different, and computer modelling offers a valuable tool for identifying preferred siting of heteroatoms, as evidenced by previous successful applications to the investigation of Si–Al24–31 and Si–Ge distribution in zeolites.12,32–35
Heteroatoms can be distributed over the framework of zeolites with varying degrees of order, from full ordering in some cases to completely random distribution in others. Computational studies of heteroatom distributions might suffer limitations in zeolites with small concentration of heteroatoms, like those exhibiting high Si/Al ratio, or with small energy differences between the configurations of foreign atoms over distinct tetrahedral sites. In such cases, the location of heteroatoms is often random or it is directed by the synthesis conditions.30,36–40 However, the larger the energy differences between configurations, the more important thermodynamic factors will be in controlling the distribution. This explains the appearance of some – particularly naturally occurring – zeolites with ordered framework-heteroatom distributions, e.g. goosecreekite,41 edingtonite,42 thomsonite,43 or yugawaralite,44 and partial ordering in many others, such as brewsterite,45 HEU-topology clinoptilolite46 and heulandite,47 epistilbite,48 and levyne.49 Computational modelling have successfully shown the preferential ordering in goosecreekite28,29 and partial ordering in HEU-type zeolites26,49,50 in agreement with experimental results. On the other hand, the limitation of computational studies in the calculation of heteroatom distributions based on T site energetics becomes important in cases like ZSM-5 (MFI), with low aluminium content, where calculations reveal31 that there are only small energy differences between configurations. This is consistent with the experimental observation that the Al location is a function of synthetic conditions in such cases.30,36–40 Obviously, a careful scan of the configurational space of a given zeolite can provide the energetic spectrum of heteroatom configurations, and the associated occurrence probabilities on the assumption of thermodynamic equilibrium. Based on this, one can recognise those cases where the energetics of the configurations is likely to control the distribution of the heteroatoms. The role of the energetic contribution in the Si–Al distribution of a selection of zeolites was studied by Zwijnenburg and Bromley,51 and recently enlarged to 209 zeolite frameworks.52
When modelling heteroatoms distribution in zeotype frameworks, interatomic potential methods have proven particularly successful as they are low-cost and have been shown to be able to reproduce subtle structural and energetic differences: quantum-based methods remain prohibitively expensive except for considering single substitutions in a unit cell or when considering a small subset of ordered structures. But when the configurational space is as large as it is in the present study (or in similar problems in related materials) the computational expense of interatomic potential methods still remains a restriction.
The local geometry of germanate tetrahedra, besides allowing stabilization of small rings, facilitates the substitution of Al atoms in structural units rarely observed in silica and aluminosilicate zeolites, such as 3-membered rings (3MR) and spiro-5 units.53 This structural diversity might lead to a wide distribution of Al atoms in germanate frameworks, as compared to aluminosilicates. However, when the Al content is relatively high (Ge/Al ratio below ca. 4), identifying the location of Al atoms is both experimentally (for the reasons given above) and computationally challenging, the latter due to the very large size of the configurational space. As an example, we consider the aluminogermanate PKU-9 (PUN IZA topology), which exhibits a zeolite framework composed of zeolite CGS layers and spiro-5 units53 (Fig. 1). This structure has 5 distinct T sites, one located in the centre of the spiro unit (T5) and two around this centre (T1 and T2), with each of these three sites being part of a 3MR. Structure refinement of X-ray diffraction data53 suggests that the Al and Ge are randomly distributed over the five T sites. Moreover, PKU-9 has only been prepared with a Ge/Al = 3.5 and no report is found of higher Ge/Al. However, the distribution of heteroatoms in related zeolites with CGS structures varies: in phosphates ordering is found for the phosphorous and other tetrahedrally-coordinated atoms (e.g. Ga and Zn or Co) at a P/heteroatoms ratio of 1,54,55 but in CGS gallosilicates, Ga and Si atoms remain disordered even when Ga/Si is only about 2.56,57 These contrasting results suggest that we may expect some ordering in PKU-9, with significant differences in the local environment of the different T sites. Moreover, the absence of a range of PKU-9 compositions might also allude to particular topological constraints on the number and location of aluminium atoms in this particular structure. For example, if at low aluminium content the stability of the structure is strongly dependent on particular T-site occupation, hydrothermal self-assembly of such a structure is unlikely. However, in order to explore computationally such a wide range of compositions, to probe any changes in preferred occupation, we have to screen as much of the configurational space as possible, which will be limited by computational cost. This is a common issue when probing distribution configurations in materials chemistry. Therefore, we propose here an approach to model the Ge/Al distribution in complex aluminogermanates by introducing an effective Hamiltonian for a fast evaluation of configurational energies. In the 1980's, effective Hamiltonians were developed to study the siting of extra-framework cations in zeolites58,59 and Al distribution in high-symmetry zeolites were also studied shortly after.60,61 These early investigations did not consider relaxation effects and also were limited to screening relatively small configurational spaces. More recently, Monte Carlo simulations have been used to study Si–Al distribution in zeolites via sampling of non-relaxed configurations.62 However, for some zeolites relaxation effects are known to significantly affect the distribution of heteroatoms, and indeed there may be no correlation between the energy of the configurations before and after structural relaxation, which limits the applicability of such fast-sampling approaches.26 Hence, new developments are still needed to deal with large configurational spaces in complex zeolites. In the approach that we introduce here, an effective Hamiltonian suitable for fast sampling of a very large configurational space is parameterised to reproduce atomistic simulations which included full geometric relaxation. In this way, information about site-specific geometric relaxation behaviour is (indirectly) included in the model.
Common theoretical studies of heteroatom distribution in the framework of zeolites have typically considered localization of the heteroatoms as single substitutions22,24,25,30,63–66 or randomly distributed over the framework34,35,67,68 with the constraint of avoiding Al–O–Al linkages in accordance with Loewenstein's rule.69 The first scenario is only comparable to experimental situations with very low heteroatom concentrations, whereas the second scenario is only applicable to high-symmetry zeolites prepared mainly with monovalent cations. Alternatively, symmetry considerations can be used to explore heteroatom distribution in a larger configurational space.26–29,31,70 However, in zeolites with high concentrations of heteroatoms even with the use of symmetry and of some ad hoc structure-related constraints, the number of independent configurations is larger than what is tractable by fully atomistic simulations. Our selection of the aluminogermanate PKU-9 zeolite, will provide here an example on how the structural information regarding the location of the Al atoms will vary according to the approach used as the composition changes from that found experimentally to higher Ge/Al.
Our calculations give access to the energies of a very large set of configurations of heteroatom distributions, which in the first place provides a measure of the likelihood of thermodynamic control of the heteroatom distribution. As mentioned above, large energy differences suggest partial or even full ordering, whilst small energy differences may favour a random action distribution or could lead to the distribution to be strongly dependent on synthesis conditions. It is also worth recalling that during synthesis, when energetic conditions are favourable, the resulting zeolites often undergo Ostwald ripening, transforming to more stable structure.74 Indeed such transformation can also occur without changes in the zeolite topology, modifying the heteroatom distribution to gain stability.75
In order to reproduce the large polarizability of the O2− anions, the shell model of Dick and Overhauser is used,84 which is important for stabilizing low-symmetry structures.82,83,85,86 The germanate framework is modelled using the Ge–O potential of Sastre and Gale,33 while the Al–O interactions are described by the Jackson and Catlow potentials.87 In general, mixing potentials from different sources can reduce the accuracy of the calculations. However, in the present case this is not problematic as both sets of potentials were parameterized for zeolite-like materials, starting from the same O–O interaction potential as the energy reference, and therefore they are expected to be compatible. Any remaining slight inaccuracies will be masked in the context of the effective Hamiltonian approach. A potential that also considers the same energy reference and includes Ge–O and Al–O interactions has been reported by Sastre and Gale.88 To provide further confidence in the method and the potentials used we selected 60 configurations and re-minimised using the Sastre and Gale potential. These configurations were selected from the entire energy spectrum, with the same number of configurations for each decile of ordering of the lattice energies computed with the other potential. The correlation between both sets of energies is 0.991, which provide us confidence on our selected parameters.
The introduction of heteroatoms into a tetravalent framework is typically accompanied by extra-framework charge compensating cations. Under synthesis conditions these are usually the cationic template molecules and/or inorganic cations.74 However, explicit consideration of these species add even more complexity to an already large configurational space, as modelling such extra-framework species would imply an additional large set of configurations for each framework heteroatom distribution.29,50 The lower charge of Al can also be compensated by protons, forming Brønsted acid sites. However, charge compensating protons are incorporated in a second synthesis step, via calcination, after the Al atoms are already distributed, so the proton siting should not affect the Al distribution, as has been argued before by Grau-Crespo et al.27 and Sastre et al.89 In order to circumvent this problem, we have omitted the explicit consideration of the extra-framework species and the charge imbalance has been compensated by using a charge compensation background.90
(1) The energy to substitute one Ge atom with one Al in the periodic cell, ΔE1[pi], is calculated from the difference between the GULP lattice energy of the cell with one Al/Ge substitution, and E0, for each tetrahedral site pi.
(2) A second-order term, denoted as ΔE2[pi, pj] and characterizing the Al–Al interactions, is then calculated from the GULP lattice energies of the configurations with 2 periodic Al–Ge substitutions (by subtracting E0 and ΔE1[pi] and ΔE1[pj]).
(3) Then, for n Al/Ge substitutions per unit cell, the energy can be approximately calculated as:
(1) |
To validate this approach, we considered two sets of structures extracted from the full configurational space of configurations with 8 Al per unit cell. In the first (Set 1), we analyse the 10000 lowest-energy configurations given by eqn (1). Then, in order to span a wider range of energies, we also consider a second set (Set 2) of 150 000 configurations, randomly selected, but with effective-Hamiltonian energies above those of Set 1. The configurations in both sets were subject to full lattice energy minimization using GULP. The multiplicity (number of symmetry redundant configurations) of each one was determined considering the complete set of configurations (∼30 million). Then, the lattice energies in conjunction with the multiplicity of the configurations were used to determine the occurrence probability and therefore the T site occupancy. The multiplicity of sites T1 to T4 is 8, while it is 4 for T5 (Fig. 1).
Fig. 2 Fully atomistic lattice energy (from interatomic potentials using GULP) as compared to effective Hamiltonian lattice energy (eqn (1)) for the Ge–Al configurations of Set 1. The solid red line represents the best linear fit and the green dashed line represents the function y = x. |
We also find that the EH approach captures the main trends of stability across the full configurational space: Fig. 3 combines the results for all the configurations in Set 1 and Set 2. Indeed, the correlation factor is improved to 0.993, with the overall dispersion decreasing for the higher energy configurations – recall Set 2 is a random sample from the entire configurational space. This analysis would suggest that the EH successfully identifies both lower and higher energy configurations well, emphasizing its usefulness as a screening tool. The ability of the EH to determine relative energies is also accentuated by the fact that only a small number of configurations (80 of 150000) from Set 2 subsequently fall into the energy range of Set 1 when subjected to full energy minimization using interatomic potentials. Moreover, the most stable of these is still ranked only 5221 of all those structures now fully optimized.
Fig. 3 Full atomistic lattice energy (interatomic potentials using GULP) as compared to effective Hamiltonian lattice energy (eqn (1)) for the Ge–Al configurations of Sets 1 and 2. The solid red line represents the best linear fit and the yellow dashed line represents the function y = x. |
We therefore conclude that the constructed EH correctly ranks configurations and reproduce energy differences within this large configurational space, giving us confidence that it is an appropriate and effective tool in screening such materials.
Providing a quantitative description of the Al–Al interactions is complex but can be illustrated by considering first the simplest case of incorporating two Al per unit cell. In Fig. 4 are displayed the energies of the compositions with 2-Al incorporated as a function of the sum of the energies of the configurations having 1-Al substituted in each site. Hence, deviation from the line with slope 1 reflects the Al–Al interaction energy. While in ZSM-5 the values of the interaction energies are about 75 kJ mol−1, excluding non-Lowenstenian configurations,31 here we have some values up to 25 kJ mol−1 higher. This is associated with the smaller size of the unit cell, which makes the Al–Al interaction larger. Note that those configurations where the interaction energy is most significant (>125 kJ mol−1) are those that have adjacent T sites occupied by Al, suggesting (as is common) that Lowenstein's rule (formulated obviously originally for aluminosilicates) is generally obeyed for aluminium (or even other formally 3+ metals) distribution in framework materials. These results, as we show below, further suggest that even at relatively low aluminium content (Ge/Al = 17) we may expect some degree of excess Al population in T5 and conversely lower than expected content in T1 and T3.
Fig. 5 Dependence of the occurrence probabilities versus the relative energy of the configurations: (a) 1 Al atom, (b) 2 Al atoms, (c) 4 Al atoms, and (d) 8 Al atoms. |
The above results, particularly for Al atoms content up to 4 atoms per unit cell, therefore suggest that the Al average occupancy in the tetrahedral sites is not likely to be random. Site occupancy depends on the multiplicity of the sites: 8 for T1 to T4 and 4 for T5. However, we also find that the preferred location of Al varies with composition, as shown in Fig. 6. At low Al content, from 1 to 4 Al per unit cell, preference for T5 is dictated by the topology. However, when 8 Al are present in the cell, which corresponds to the experimental composition, the Al–Al repulsion (Fig. 4) becomes more relevant than the preference associated to topology. In this case, sites T1 and T2 show the higher populations, with T3 now being the least favoured. Noticeably, the differences between the various sites are now less pronounced.
In order to further quantify the heterogeneity of the site occupancy for each Al content, we use the standard deviation of the Al site occupancy probabilities as a measure, i.e. the lower this value the lower the heterogeneity, and thus the larger the likelihood of random Al siting. The calculated standard deviation values are 0.20, 0.16, 0.15, 0.24, and 0.08 for Al contents 1–4 and 8 atoms per cell, respectively. The value for 8 Al atoms per cell is at least twice lower than for the other cases, which indicates that the bias for single site occupation due to the framework topology is largely smeared by the Al–Al interactions and the local distortions accompanying the Al incorporation. However, we also note that the computed standard deviation of the average probability by Al atom for the highest Al content, 8 Al per unit cell, is about 2.7 times higher than that determined in the experimental work (0.03).53 In other words, our theoretical analysis suggests that Al siting in PKU-9, in thermodynamic equilibrium, should be less random than what was determined by XRD. This discrepancy could be due to the approximations involved in our model, or to non-thermodynamic effects (which we ignore) contributing to the actual distribution. Alternatively, it could also be a consequence of the limitations of XRD analysis in determining site occupancies for this type of system. This point therefore deserves further theoretical and experimental investigation. It may be useful to repeat a refinement starting from different occupations, such as the most stable and least stable determined here, to establish if the experimental data can indeed be used to distinguish any preferential siting. It is worth noting that the type of effective Hamiltonian introduced in this work has been recently used, along with experiments, to investigate the Si–Ge substitutional series in the chiral STW family of zeolites.98
As a case study, a detailed computational study of the Al distribution in the framework of PKU-9 aluminogermanate was presented. Al contents of 1 to 4 and 8 (experimental content) atoms per unit cell were considered. We have been able to reduce the size of the configurational space of Al distribution using the SOD code91 for Al contents up to 4 Al atoms per unit cell. For 8 Al atoms in the unit cell, a selection approach based on the stability order indicated by the effective Hamiltonian approximate energy was implemented. At Al content up to 4 atoms per cell, clear preferential Al location is predicted for T5 site. However, sites T2 and T3 are expected to show higher occupancies at the experimental composition (8 Al per cell). This reveals that the Al–Al interactions are crucial in controlling the Al distribution in aluminogermanate, and therefore information gained for lower Al incorporation cannot be simply extrapolated to higher Al content.
Conversely, the lack of significant preferential siting at the experimental compositions may provide an insight into why it is this particular composition that can be formed: for higher Ge/Al ratios, specific Ge siting may be required to form stable structures. Su et al.'s analysis of the XRD data suggests a purely random distribution,53 which our data generally supports. Nevertheless, we do predict some non-homogeneous distribution of the Al atoms over the T sites at the experimental composition. It can be argued that, without a suitable starting model, Rietveld analysis methods may not converge to anything other than a random distribution. Therefore, it would be valuable to perform further experimental studies, perhaps using Al NMR data to combine with the XRD data or using our calculated distribution as a starting point in the Rietveld analysis.
Footnote |
† Current address: Department of Physical, Chemical, and Natural Systems, University Pablo de Olavide, 41013 Seville, Spain. |
This journal is © the Owner Societies 2018 |