Molecular cluster models, developed for an exfoliated kaolinite, provide a structural description comparable to that of periodic slab models for a fraction of the computational cost. These models include both the octahedral and the tetrahedral sheets of kaolinite. The first-generation model (G1) contains the inner and outer coordination sphere of the Al- and Si-honeycombs as the preferred sites for adsorption of small organic molecules. Since no experimental information is available to date at the atomic level for exfoliated kaolinite, we carried out a systematic density functional theory evaluation for establishing the most reasonable coordinates of the ions and groups. The results of molecular cluster and periodic calculations were utilized for evaluating semi-empirical Hamiltonians on larger models. Using a PM7 Hamiltonian, the structure of cluster models containing 1 + 6 (second generation) and 1 + 6 + 12 (third generation) Al- and Si-honeycombs which are out of reach for ab initio calculations, were determined. These molecular slab models offer a structural platform for adsorption, intercalation and delamination studies.

The reduction of crystalline order in kaolinite by intercalation, delamination and exfoliation increases its reactivity and thus its industrial value. The exfoliated kaolinite is a remarkable raw material that manifests rather different chemical and physical properties than the original crystalline kaolinite. The chemical composition of the exfoliated and the original crystalline kaolinite remains the same upon exfoliation; however, the morphology changes dramatically from periodic pseudo-hexagonal plates to nano-sized chips, platelets, tubes and ribbons (Singh & Gilkes, 1992; Tsunematsu & Tateyama, 1999; Gardolinski et al., 2001; Gardolinski & Lagaly, 2005; Valaskova et al., 2007; Horvath et al., 2010; Kuroda et al., 2011; Sun et al., 2011). Due to the decrease in the crystalline order, the common structural determination techniques, such as X-ray diffraction (XRD) and Raman spectroscopy, cannot be employed. Fourier-transform infrared (FTIR) spectroscopy may provide indirect information about the movements and chemical environments of groups and ions. The EXAFS analysis of X-ray absorption spectra provides some structural information in the form of radial distribution of atoms around Al and Si absorbers, or average Al–O and Si–O distances (Gualtieri et al., 2000; Shaw et al., 2009; Gates, 2013), where the fine structural details are lost due to the averaging of scattering pathways. To date, no direct experimental technique is available for obtaining atomic positions in exfoliated kaolinite. This would be of great importance, since the knowledge of structure, particularly the differences between exfoliated and platy kaolinite, hold the key for rationalized design of clay-based materials with tailored functionality.

Computational models (QM) can fill the gap in our knowledge of chemical structure by providing atomic level geometric and electronic structural information. The aim of this work is to validate QM methods and basis sets and to extend our earlier studies on molecular cluster models for non-periodic systems to larger, compositionally and chemically more realistic models. These in silico virtual chemical models have an enhanced potential to connect experiments directly with theory. Modern density functional theory (DFT)-based models (White et al., 2009; Mercier & Le Page, 2011; Dawley et al., 2012; Tunega et al., 2012) allow prediction of reactivity with regards to interaction with solvent environments and small inorganic or organic molecules without the need for an extensive empirical training set, but at a considerable computational cost. In contrast, molecular mechanics (MM) calculations (Teppen et al., 1997; Cygan et al., 2004; Heinz et al., 2013; Cygan & Tazaki, 2014) utilize classical mechanical potential energy descriptions and their parameters that are defined on the basis of molecular chemical environments or atom types without the need for any knowledge about electronic quantum states. The advanced molecular mechanical treatment provided by the INTERFACE molecular mechanical force field (Heinz & Suter, 2004; Heinz et al., 2013; Emami et al., 2014), represents considerable improvements over earlier MM force fields in which the atomic charges and other parameters reflect chemical knowledge of the distribution of the electron density, reactivity and polarizability from experiments. Regardless of how accurately a given force field parameter set is defined, the transferability of MM parameters from a crystalline, periodic system to disordered, non-periodic molecules, such as those described in this work, can be a prohibitive bottleneck.

The conceptual increase in computational accuracy for realistic size models often comes at a prohibitive computational cost. Semi-empirical Hamiltonians (Stewart, 2007; Stewart, 2008; Clark & Stewart, 2011; Stewart, 2013) offer a compromise between the accuracy of the level of theory, the computational model size and the computational cost. These methods and their parameters describe the structure of valence shells through the quantum mechanical treatment of electrons. Thus, their atomic parameters are less dependent on the molecules. This advantage of semi-empirical Hamiltonians is, however, their greatest weakness, as atoms do behave differently in different molecules. However, careful parameterization of integrals expressing electron/electron interactions can mitigate these limitations and accurate semi-empirical calculations have already been reported for minerals (Stewart, 2008).

Figure 1 (Smrcok et al., 2010) summarizes the terminology used in this work to describe the structural features of the octahedral (O) and the tetrahedral (T)-sheets or the OT-layer of kaolinite. These sheets are formed from the coordination environments of octahedral Al3+ and tetrahedral Si4+ ions, respectively. The surface hydroxide groups (s–HO) cover the O sheet. They are differentiated with respect to being above (α) the inner hydroxide group (i–HO) or adjacent (β) to it. There are six s–HO groups for a [6Al–6OH] honeycomb-like unit (Al-honeycomb, Fig. 1B), thus they need to be further divided up relative to the centroid of the six s–HO groups (XsO, upper and central pink sphere in Fig. 1A and 1B, respectively). The three s–HO groups pointing towards the XsO centroid are labelled as proximal (one αp and two βp), while the other three pointing away from the XsO are called distal groups (one αd and two βd). Analogous nomenclature applies for the bridging oxide groups (b–O2−) in the [6Si–6O] honeycomb-like unit (Si-honeycomb, Fig. 1C) with respect to their centroid (XbO) as shown for the T-sheet in the lower part of Fig. 1A and in Fig. 1C. The approximately equal distribution of the apical oxide groups (a–O2−) connecting the O- and T-sheets did not require a site-selective labelling.

The goal of this study is to provide experimentally sound and theoretically consistent estimates for the atomic positions in the exfoliated kaolinite using molecular cluster-based DFT models, where the truncation of the models is dictated by the coordination chemical environment versus the availability of computational resources. The model development was carried out in comparison with 2D and 3D periodic structure models, where only the latter model corresponds to experimental structural information. These ab initio cluster models serve as a reference for evaluating the performance of semi-empirical methods that will, in turn, allow a significant increase in the size of models by at least one order of magnitude towards realistic, nanometre-scale particles.

Computational Methods

The DFT calculations using a range of functionals were carried out using Gaussian suite programmes (Frisch, 2009). As described in a previous publication (Táborosi et al., 2014), the selection of the functionals was deliberate so as to cover the most commonly used representatives for Generalized Gradient Approximation (GGA), hybrid GGA and metaGGA functionals. While each individual functional is expected to give different numerical results, a trend was observed that all GGA functionals perform similarly and the results of the hybrid functionals depend on the amount of Hartree-Fock (HF) exchange regardless of the approximation used to construct a given functional.

The family of GGA functionals was represented by BP86 (Becke, 1988; Perdew et al., 1992) and PW91 (Perdew & Wang, 1992), while the B3LYP (Lee et al., 1988; Becke, 1993) commonly used functional was employed as a hybrid GGA functional. The TPSS was chosen as an example for metaGGA-type functional (Staroverov et al., 2003; Perdew et al., 2004; Staroverov et al., 2004). It is worth mentioning that pure GGA functionals often manifested convergence issues due to oscillations during the optimization of the electronic structure. This can be practically eliminated by embedding the computational models in a polarizable continuum model (Tomasi et al., 2002; Klamt, 2011), which is also desirable in modelling the electrostatic influence of the external chemical environment. In order to evaluate the difference between DFT- and wave functional-based methods, pure HF (Slater, 1951) calculations were also carried out. Given the strong basis set dependence of the optimized structures described earlier (Táborosi et al., 2014), only the results using the triple-ζ quality def2TZVP (Schäfer et al., 1994) basis set (abbreviated as TZVP) is reported in this study as a good approximation for the basis set saturation limit for our cluster models. The periodic 2 × 2 slab models were calculated with the SVP (Schäfer et al., 1992) basis set within the periodic boundary condition (PBC) treatment (Kudin & Scuseria, 1998; Kudin & Scuseria, 2000; Kudin et al., 2001) due to the staggering computational cost of using the def2TZVP basis set (Table 1). The PBC result at the BP86/def2TZVP level was obtained with a cumulative wall clock time of about three months for complete structural optimization on a 12-core AMD Opteron 8380 2.5GHz server (see Table 1 for more detailed timing).

All structures reported here correspond to stationary structures obtained after repeated restarts to avoid being trapped in a local minimum. Unconstrained optimizations were carried out for the periodic slab models including the unit cell parameters, while the positions of the counter ions and the atoms at the periphery or outer sphere region were kept frozen for the cluster models. The effect of the size of the periphery region was evaluated by optimizing either the entire inner sphere region or only the atoms corresponding to the Al- and Si-honeycombs. To establish a G1 reference structure, the second approach was chosen by allowing the outer sphere moieties to impose a maximal strain on the inner-sphere atoms. The set of optimized atoms formed a leaning barrel shape capped by the Al- and Si-honeycomb units. Due to the size of the model and the associated computational cost, vibrational analysis was not carried out. All calculations started from the crystalline structure of the kaolinite or its single TO-layer. In addition, any reactivity in the current study was not considered; therefore, normal modes with imaginary frequencies are not expected. The semi-empirical calculations were carried out using the MOPAC2012 package (Stewart, 2012). The results of a PM6 (Stewart, 2007) Hamiltonian with D3 dispersion (Grimme, 2006) or DH2 dispersion and hydrogen bonding (Rezac et al., 2009; Korth et al., 2010) corrections were compared to the PM7 Hamiltonian developed recently (Stewart, 2013).

The detailed steps and rationale for the development of a coordination chemistry-based cluster model of the [6Al–6OH] honeycomb were presented earlier (Táborosi et al., 2014). In the given study, the same construction strategy was employed, but simultaneously, for the Al- and Si-honeycombs, resulting in the Generation 1 (G1) cluster model of the exfoliated kaolinite (Fig. 2A). The molecular charge of the G1 model with Na+ and Mg2+ counter ions is +8. In order to neutralize G1, the peripheral region was extended with four additional Al3+-coordination environments above the two Si4+-ions that reach beyond the Al-honeycombs due to the slip of the Al- and Si-honeycombs relative to each other (Fig. 2B). Consideration of further outer-sphere coordination environments results in the Generation 2 (G2) cluster model that contains the inner Al- and Si-honeycombs, its first outer sphere environment (G1) and the next two outer spheres (Fig. 2C). Analogously, the G3 model has all the honeycombs of G2 and 12 additional from the periphery (Fig. 2D). The size of the models and number of electrons are summarized in the caption of Fig. 2. The molecular charge of the G2 and G3 models are negative, thus 2 and 12 Na+ ions were replaced symmetrically by Mg2+ ions from the outer sphere of most counter ions to get neutral computational models. Initial and optimized atomic position coordinates for all computational models are given as supporting information at http://computational.chemistry.montana.edu/clay3.

Results and Discussion

Evaluation of DFT periodic boundary condition calculations for crystalline kaolinite

First, a brief section with evaluation of the performance of ab initio DFT within the framework of periodic boundary condition calculations (Kudin & Scuseria, 1998; Kudin & Scuseria, 2000; Kudin et al., 2001), as implemented in Gaussian09 (Frisch, 2009), is given. Table 2 summarizes the most important geometric parameters discussed in the study for crystalline kaolinite, both from experiments (Smrcok et al., 2010) and from GGA and metaGGA functionals. The application of hybrid GGA functionals was found to be prohibitive due to computational cost.

The agreements between any of the DFT levels and the experimental internal coordinates are within 0.03 Å and 1° for distances and angles, respectively. There is only a slight improvement when going from the SVP and TZVP basis set of 0.02 Å in root-mean-square (RMS) deviation. This means that the considerably smaller SVP basis set can be accepted as a good approximation for the basis set saturation limit for the kaolinite structure. The use of the TZVP basis set for the periodic geometry optimization would have been more expensive, computationally, without resulting in an improved molecular geometry.

The largest deviations were observed for the inner environment with regards to the orientation of the i–HO groups as this is one of the most sensitive moieties of the kaolinite structure with respect to concerted ion/dipole b–O2−…H(i–HO) interactions, H-bonding and the ionic vs. covalent nature of the Al–O(i–HO) bonds. The metaGGA level with the TPSS functional shows the best agreement, within 0.02 Å, with the experimental value of the b–O2−…H(i–HO) distance, while the same level of theory gives the worst agreement of the b–O2−…O–H(i–HO) angle with +5° deviation from the crystal structure. The trend is completely reversed for the GGA functional BP86 with deviations of 0.09 Å and −1° in distances and angles, respectively. The PW91 GGA functional offers a slightly improved agreement in the b–O2−…H(i–HO) distance with 0.06 Å, which can only be improved by 0.01 Å using the larger basis set, TZVP. To achieve an even better reproduction of the experimental structure from ab initio calculations, dispersion corrections to DFT need to be considered. It has been discussed in the literature that standard DFT functionals fail to predict structural parameters correctly when non-bonding interactions are important (Tunega et al., 2012). The effects of dispersion correction and the presence of a polarizable continuum will be discussed below for the exfoliated kaolinite models.

Establishing the reference structure for exfoliated kaolinite

The G1+ model was used to evaluate the performance of various density functionals. Figure 3 illustrates the result of the structural optimization graphically. Only those atoms that are shown in Fig. 3A were optimized. The arrows in Fig. 3 indicate the movements of the groups and atoms relative to their position in the original kaolinite (Fig. 1). Upon exfoliation, the s–HO groups cannot form H-bonds in the absence of an adjacent TO-layer, thus their positions are altered depending on the external chemical environment. In gas phase simulations, some s–HO groups remain perpendicular to the O-sheet; however, most of them fold down and became co-parallel with the O-sheet. As noted in previous work (Táborosi et al., 2014), the Al3+-ions undergo a large-scale structural rearrangement as they sink into the OT-layer, closer to its central plane. The a–O2− groups move slightly away from the b–O2− groups and they move towards the centre of the G1+ model. The Si4+-ions shift away from the centre of the TO-layer, but to a smaller extent than the Al3+-ions. The differences between the proximal and distal groups in Figures 3B and 3C became much more pronounced due to changes in the Si–O–Si angles than observed for the crystalline structure in Figures 1B and 1C.

In addition to a graphical illustration of typical structural changes in going from a crystalline (Fig. 1) to an exfoliated (Figure 3) kaolinite TO-sheet, Table 3 summarizes the numerical values for selected internal coordinates from optimized structures. Their corresponding values from the crystalline periodic structure are shown in Table 2. It has to be emphasized that the exfoliated and crystalline structures cannot be compared directly because they are different chemical species. The orientation of the s–HO groups shows a clear trend in all the molecular cluster models independent of the computational level of theory. Both βps–HO groups will be the only ones that remain pointing upward as in the crystalline structure. The presence (B3LYP/TZVP/PCM) or absence (B3LYP/TZVP) of a polarizable continuum as an electrostatic model for the external chemical environment does not affect this trend significantly in DFT calculations. The results of periodic slab calculations for a 2 × 2 cell show a slightly different arrangement, where only a few s–HO groups undergo a change in their orientations. It is important to highlight that the α–s–HO groups that are above the i–HO groups, regardless of whether they are in proximal or distal positions relative to the XsO centroid, always become parallel with the plane of the kaolinite sheet. It is also notable that the periodic calculations using the SVP basis set resulted in some disorder within the 2 × 2 slab model with respect to s–HO arrangement, which did not show up for either of the calculations employing the larger TZVP basis set. A slight expansion (0.01–0.02 Å) of the single TO-layer in the ‘a’ and ‘b’ crystallographic dimensions can be observed in the exfoliated model relative to the crystalline structure resulting in expansion of the unit cell as shown at the bottom of Table 3, while the angle of the slab unit-cell directions remain practically the same as in the original structure. As observed for the crystalline unit-cell calculations the level of computational theory does not affect the calculated structure significantly. A small difference between the TZVP and SVP basis sets here supports the basis set saturation with respect to the geometric structure already at the latter level. The deviation between the cluster and the periodic models with respect to the s–HO arrangements and its energetic importance will be discussed below.

The Al–O (s–HO) distances become much longer by about 0.04–0.08 Å in the exfoliated structures compared to the crystalline phase. All models reproduce this change. The periodic boundary condition models show the largest deviation of 0.08 Å. Upon changes in the Al–O distances, the corresponding Al–O–Al angles also change in the optimized structures relative to the original phase. The 2 × 2 periodic slab models show a 3–5° decrease of this angle as a result of Al3+-ions sinking toward the middle of the TO-layer, referred to as inner environment of the TO-layer. These trends were reproduced quantitatively at all levels of computational theory for the cluster model with a slight variation of 1–2°. As a further consequence of the Al3+-ion movement, the Al–O(i–HO) bond lengths pointing toward the inner environment of the TO-layer become shorter by ~0.03–0.05 Å. Although there are some numerical differences among the functionals, they all show very similar values to within 0.03 Å.

It is important to highlight that a limitation of the G1+ cluster model was found during this study; it does not describe the chemical environment of the i–HO group that points away from the Al-honeycomb fully. This can be seen from the different values obtained for the Al–O(i–HO) bond lengths for the molecular cluster model. This is not unexpected as the G1+ model is essentially the fusion of earlier models (Táborosi et al., 2014) for the [6Al–6(s–HO)] and [6Si–6(b–O2−)] honeycombs at the outer surface and does not describe the two i–HO groups adequately. A more complete description of their environment is given by the G2 model as shown below. The presence of a continuum model for the external environment only affects the Al–O(i–HO) distance (1.89 Å vs. 1.91 Å for B3LYP/TZVP) notably, where they become similar regardless of whether the chemical environment of the i–HO is complete or not.

Analogous differences can be observed for the inner chemical environment of an OT-layer amongst the crystalline, periodic slab and the molecular cluster structures. There is asymmetry in the Al–a–O2− bond lengths in the crystalline structure (Table 2). This was reproduced at every level of computational theory, but to a different extent. The slip or misfit between the Al- and Si-honeycombs induces longer Al–a–O2− bonds (2.04(1) Å vs. 1.97(1) Å, calculated values) that are closer to the i–HO group pointing inwards to the Al-honeycomb in the crystalline structure (2.02 Å vs. 1.96 Å). However, this asymmetry is essentially eliminated for the 2 × 2 periodic models (1.95(1) Å vs. 1.94(1) Å). The gas phase models, either at HF or B3LYP level, show the distortion in the Al–a–O2− bonds (2.00(2) Å versus 1.87(2) Å); however, upon using a continuum model, the differences become similar to those obtained for the slab models again without large variations amongst the functionals considered here (2.00(0) Å vs. 1.94(0) Å). These structural changes correlate with the expansion of the cluster or periodic slab in both ‘a’ and ‘b’ crystallographic directions. The diagnostic orientation of the i–HO group for the external chemical environment to the TO-layer was investigated in detail for a model that contained only the Al-honeycomb (Táborosi et al., 2014). The current results obtained for a more complete computational model in Table 3 confirm the earlier findings for internal coordinates b–O2−…H(i–HO) (δ), O(s–HO)O–H(i–HO) (α1) and b–O2−…O–H(i–HO) (α2) from a smaller model (Fig. 4 in Táborosi et al., 2014); therefore, they are not discussed in detail. In brief, the i–HO group responds to lack of an external environment at the T-sheet of the exfoliated model by changing its arrangement from parallel to approximately 45° relative to the plane of the TO-layer due to a strong H-bonding interaction with the adjacent b–O2− group. The sensitivity of this interaction was also observed at the periodic crystalline cell calculations as discussed above and is shown in Table 2.

A new structural element that was taken into account in this work is the explicit consideration of the T-sheet in the cluster model, which is often omitted in the literature. The sheet of Si4+-ions and their coordination environment are required for the cluster model according to coordination chemistry principles when considering the entire TO-sheet in adsorption intercalation, delamination and exfoliation studies. The short Si–O(a–O2−) bonds from the crystalline kaolinite becomes elongated by 0.03–0.06 Å in the periodic slab models. This distortion was well reproduced in the cluster model (0.02–0.04 Å). The use of a polarizable continuum environment reverses the elongation of the Si–O(a–O2−) bonds, as the continuum model has a similar effect on the adjacent OT-layer in the crystalline kaolinite. It is remarkable then, that despite all the changes described above, the Si–O(b–O2−) bond distances remain practically the same as in the crystalline kaolinite, to within 0.02 Å. This is only possible if there are significant changes in the Si–O–Si angles (see Table 3). In contrast to the Al–O–Al angles of the O-sheet, these angles open up by 2–11° compared to those in the original platy phase. The use of the continuum model also affects the Si–O(b–O2−) bonds in the molecular cluster model of the exfoliated kaolinite making them very similar to those of the original phase.

As overall benchmarks for the difference between the crystalline, periodic and cluster models for exfoliated kaolinite, root-mean-square (RMS) values were calculated for distances and angles weighted by their number of occurrences. Supporting information deposited with the Journal contains a table with numerical comparisons. The periodic slab models show 0.09–0.11 Å and 13–16° RMS deviations relative to the original structure. This can be taken as the degree of difference between the experimental structures of the exfoliated (unknown) and the original (discussed above) structures, due to the excellent agreement between the calculated and experimental structures for the original phase independent of the functional when a saturated basis set was used. The molecular cluster models display a wider range for the distances with 0.03–0.07 Å and a somewhat narrower range for the angles, 3–14°. Considering the conceptually highest level of the BP86/TZVP PBC model, at unfortunately prohibitive computational cost for routine work, the RMS deviations of the HF, the B3LYP and all the PCM models are 0.05 Å and 4°, 0.03 Å and 2°, 0.05–0.07 Å and 5–12°, respectively. These numbers reveal a trend that was emphasized previously (Táborosi et al., 2014), i.e. that the HF method can actually provide a reasonable structural description for the exfoliated kaolinite. The closest molecular representation of the molecular cluster model to the highest level of a 2 × 2 periodic structure was achieved by using the B3LYP/TZVP method without (0.03 Å and 2°) or with (0.05 Å and 5°) a PCM environment.

Given the large number of individual deviations relative to those of the original kaolinite described above and presented in Table 3, we completed an additional round of analyses for the movement of five parallel hexagons formed by six O(s–HO), Al3+, a–O2−, Si4+ and a–O2− ions by comparing the displacements of their centroids (XsO, XAl, XaO, XSi and XbO, respectively) relative to the constrained centroid positions from the periphery of the model. Table 4 compares the total and only the displacements along the ‘c’ crystallographic direction, which is approximated by the ‘z’ directional displacements. Similar analyses could not be performed for the periodic models due to the movement of all the atoms and thus the disappearance of any reference positions. However, the changes in internal coordinates in Table 3 suggest that the movement of the sheets in the 2 × 2 periodic models do follow similar trends.

With the exception of the XaO centroid for the a–O2− ions, all centroids show displacements predominantly along the ‘z’ directions perpendicular to the TO-layer. The signs in Table 4 indicate the direction of the movement. As suggested earlier (Táborosi et al., 2014), there is a considerable contraction of the layers for cluster models without a dielectric continuum model of the external environment. With respect to experimental conditions, this in vacuo calculated structure of the exfoliated kaolinite corresponds to an anaerobically heated sample below the dehydroxylation temperature at ~250–300°C. When the model was embedded into a solvent dielectric environment, the contraction of the TO-layer was reduced; however, there are still significant changes, for example, in the positions of the Al3+-ions. Figure 4 illustrates the structural contraction and/or expansion layer-by-layer for the exfoliated kaolinite graphically as a function of the external chemical environment. These individual changes may appear to be only 0.02–0.09 Å; however, collectively these will contribute to a dramatic change in the reactivity of an exfoliated surface from an original TO-layer of kaolinite. It is worth noting that the changes in the Al–O–Al and Si–O–Si angles in Table 4 are consistent with the structural modifications, due to the movements of the layers of ions and groups, shown in Fig. 4.

Semi-empirical Hamiltonians

The comparison of the molecular cluster and periodic slab calculations supports the validity of using coordination chemistry-based principles in composing computational models for exfoliated kaolinite surfaces. The G1+ model can be treated conveniently at DFT level using any modern computational facility. The molecular cluster model with the largest basis set has about a half order of magnitude advantage over a periodic model of approximately the same size at the GGA level with respect to computational cost (see Table 1) when the total number of optimization cycles needed for obtaining stationary structures is considered. In addition, the size of the G1 + model allows investigation of adsorbents up to the size of a porphin ring. It is also advantageous that practically all levels of computational theory gave similar optimized structures for the exfoliated kaolinite. Thus, the specific density functional can be selected to the specific organic/inorganic reagent without affecting the kaolinite structure considerably. The accuracy of the ab initio levels of computational theory poses a serious drawback, which is the computational cost. The employment of a more realistic cluster model, such as the G2 and above, with close to 200 atoms and 700 valence electrons may become computationally prohibitive even at the economical GGA (Table 1) level using a triple-ζ quality basis set. A compromise between model size and computational cost is given by the use of semi-empirical Hamiltonians.

Comparison of DFT and semi-empirical structures for G1+ model

Figure 5 shows the side, top and bottom views of the exfoliated kaolinite model G1+ calculated at PM7 level with the COSMO (Klamt, 2011) continuum model. It is remarkable how similar it is to the B3LYP/TZVP structure in Fig. 3 with respect to movement of the s–HO groups, displacement of the Al3+-ions, co-planar movement of the a–O2− ions, sinking of the Si4+-ions and inward movement of the b–O2− ions. Qualitatively, the PM7 Hamiltonian with a dielectric continuum model reproduces all of the structural changes observed for the exfoliated kaolinite relative to the original phase as those in the B3LYP/TZVP calculations for less than 1% of the computational cost (Table 1).

Table 5 provides a basis for the quantitative analysis of internal coordinates and selected geometric parameters. Despite the notable agreements between the DFT and semi-empirical structures in Figures 2 and 5, respectively, there are deviations that may be eliminated by improving the semi-empirical parameter set and the molecular cluster modelling approach. First, the PM6 calculations show unacceptable large structural deviations relative to the G1+ model structure calculated at the B3LYP/TZVP level. The RMS deviations in the PM6 calculations for distances and angles are 0.2 Å and 30° relative to the B3LYP/TZVP cluster model. The orientation of the s–HO groups only matches the reference structure for the a-position, while the rest is completely inverted. The Al–O(s–HO) distances are greatly elongated, by 0.04–0.08 Å, while this is considerably smaller for the Al–O(i–HO) distances. The inner environment shows a slightly better agreement with the reference B3LYP/TZVP structure due to a more complete coordination sphere; however, the position and orientation of the i–HO groups are completely reversed as they point toward the O-sheet and not. the T-sheet as in all DFT calculations. The longer Si–a–O2− distances and the shorter Si–b–O2− distances were also reversed in the PM6 structure.

The PM7 calculations produced inward-folded s–HO groups, which can be explained by the influence of intramolecular chemical interactions as an explicit external environment is absent in gas phase calculations. However, the orientations follow the same qualitative trend as in the B3LYP reference structure. The Al–O distances of the O-sheet are closer to those of the reference structure, but they are still elongated. The inner environment now looks more reasonable, including the position of the i–HO group, than in the PM6 model. Similarly, the differences in the Si–O distances of the T-sheet follow the same trend as in the reference structure. These improvements can be gauged by the reduced RMS deviations of 0.04 Å and 5° for distances and angles, respectively, compared to the order of magnitude larger values of the PM6 results. This suggests that the PM7 Hamiltonian and the corresponding parameter set are much improved relative to those of the PM6 for the kaolinite structures.

In order to tackle the inherent problem of the semi-empirical calculations with respect to the s–HO positions in the absence of an external environment, the PM6 and PM7 models were embedded into a polarizable continuum described by the COSMO model. Furthermore, an empirical dispersion and hydrogen bond correction (DH2) (Rezac et al., 2009; Korth et al., 2010) was also utilized as implemented in MOPAC2012. As expected, these affect the orientation of the s–HO groups, since they can be involved in ion/dipole interaction with the surface point charges of the COSMO cavity. Yet, the PM6 Hamiltonians with or without dispersion corrections show structures that correspond to too many folded s–HO groups relative to the reference DFT structure, while the PM7 model has predicted only one type of s–HO group that moved correctly. The RMS deviations for PM6 Hamiltonians improve to 0.10–0.11 Å and 20–22° when using the COSMO continuum model, while it gets slightly worse for the PM7, 0.09 Å and 14°, relative to the in vacuo PM7 calculation.

In summary, comparisons of the performance of semi-empirical Hamiltonians do not provide a ubiquitous favourite method with respect to the reproduction of the B3LYP/TZVP reference structure; however, the PM7 method without a polarizable continuum appears to best describe the inner chemical environment; while application of a polarizable continuum model provides an improved description of surface sites. Since the main interest is in the surface reactivity of kaolinite, all further calculations for the G1+, G2 and G3 models were carried out using the PM7/COSMO model. As described later, there will be unexpected advantages for larger models at the PM7/COSMO level, which can be related to the result of employing a more complete computational model.

Mapping of s–HO orientation for the G1+ model

A critical aspect of the models for the exfoliated kaolinite is the orientation of the s–HO groups that will determine whether the surface of the O-sheet will behave as dominantly electrophilic (as in the original phase) or nucleophilic by folding all s–HO groups downward, or amphoteric. The importance of the position of s–HO groups was realized after obtaining different optimized structures as a function of the initial structure when using the PM7 Hamiltonians (Table 6). Optimizations with a s–HO arrangement of ‘UUUUUU’, as in the original structure, resulted in different structures and energies than from the B3LYP reference structure with a ‘PPUPUP’ arrangement, where the ‘U’ and ‘P’ labels indicate upward and planar (or folded in) arrangement of the H–O bond relative to the O-sheet surface, respectively. The energetic difference of 33 and 482 kJ mol−1 are considerable for the PM7 and PM7/COSMO models, respectively. Gradual change of the arrangements of the s–HO groups from the crystalline ‘UUUUUU’ to the DFT ‘PPUPUP’ configuration allowed localization to lower energy structures than any of the initial attempts for the PM7/COSMO level. The variability of the optimized structures draws attention to a considerable conformational problem that had to be addressed to evaluate whether the presence of many local minima correspond to significant energetic preferences for a particular arrangement of s–HO groups.

The low computation cost of the PM7/COSMO calculations, allowed construction of all 64 possible combinations for the upward and parallel arrangements of the 6 s–HO groups of the Al-honeycomb in the G1+ model. Figure 6 summarizes the results of the initial structural optimizations, where the arrangement of the s–HO groups was locked by fixing the ‘z’ coordinates of the central s–HO groups within the Al-honeycomb. Upon relaxation of the ‘z’ coordinates, the s–HO groups start to fold in and result in seven stationary structures as summarized on the right hand side of Fig. 6. In all cases, the s–HO groups folded slightly below the plane of the TO-sheet; however, this essentially corresponds to the parallel orientation in the B3LYP reference structure to within a few degrees. In a similar manner to the internal coordinates in Table 5, a good quantitative agreement with the DFT calculations was not obtained; however, five of the final G1+ model structures at the PM7/COSMO level can be considered as indistinguishable energetically. This suggests that without an external environment the s–HO groups can adopt random orientations. Furthermore, the energy range of 43 kJ mol−1 for the conformation analysis in Fig. 6 is considerably lower than the energy deviation of more than 100 kJ mol−1 (Table 6), suggesting the presence of numerous local minima that the MOPAC optimizer may locate. This brings up the question about the possibility that the reference DFT structure may not correspond to a global minimum either. The detailed analysis of the potential energy surface at DFT level is the subject of on-going investigations.

Several of the relaxed structures in Fig. 5 show a s–HO group arrangement with downward instead of a parallel H–O bond orientation. Unexpectedly, the ‘DDDDDD’ arrangement was one of the lowest energy conformations, where the label ‘D’ indicates a downward orientation of the surface O–H bond. This suggests that the H-bonding or ion/dipole interactions in PM7 appear to be parameterized to be overly dominant, while ligand/ligand repulsion interactions are underestimated. This can be improved by increasing the covalent bond character of the Al–O bonds, which can result in parallel arrangement for the s–HO with the O-sheet as seen for DFT calculations. Furthermore, a big deviation between DFT and semi-empirical results is the orientation of the i–HO groups. In both models, they move from a parallel position with the O- or T-sheet; however, they move upwards in the semi-empirical calculations and downward in the DFT models. This can also be corrected potentially by further parameterization of the PM7 Coulomb and resonance integrals for the Al-centres. Regardless of the numerical differences between the DFT and semi-empirical exfoliated kaolinite models, the arrangements of the s–HO groups are determined predominantly by the interactions between the solvent, adsorbent/solvent and pure adsorbent with the surface and not just the surface itself. The nature of these interactions with explicit molecular solvent environments is also the subject of on-going investigations.

Structural optimizations of the G2 model

Another possible source of deviation between the DFT and semi-empirical models in the surface reactivity of exfoliated kaolinite is the size of the computational model itself. The terminating peripheral atoms, the counter ions, are within 5–7 Å of the optimized region of the G1+ model. In order to evaluate the effect of outer-sphere groups and ions, the PM7/COSMO method was utilized and structural optimizations for the G2 model were carried out. Figure 7 shows the optimized structures for two different regions of the G2 model, where the s–HO groups, unexpectedly show the ‘DUUDUU’ arrangements, which is now close to the arrangement obtained for periodic slab models, but not for the G1 model in either the DFT cluster or semi-empirical calculations.

The bottom right hand corner in Fig. 7A illustrates the tilted, but still parallel arrangements of the α–s–HO and i–HO groups. This is a shortcoming of the PM7 Hamiltonian in comparison to the DFT results, where the α–s–HO should be parallel to the surface and the i–HO should have a H-bonding interaction with the adjacent b–O2− ion located diagonally in the Si-honeycomb. The possibility that with parameterization of the covalent Al–O and Si–O interactions, the structural deviations between the DFT and PM7 optimized structures can be minimized is under investigation. Figure 7B hints at the appearance of a structural ordering at the kaolinite surface, where the four β–s–HO groups remain in a H-donor/nucleophilic position, while the two α–s–HO groups become H-acceptor/electrophilic sites for an Al-honeycomb as suggested by all the DFT models.

The orientation of the s–HO groups also shows some variability as a function of the relaxed and constrained regions (Table 6). When the G1 and G2 regions are relaxed, there are considerable energetic differences (~300 kJ mol−1 per honeycomb unit) due to the different number of optimized atoms and also one of the s–HO groups turns downwards for the G1 structure while the same configuration is obtained for the G2 structure as in the B3LYP/TZVP model. In comparison to the above energetic differences, we observe little difference (to within 35 kJ mol−1) whether the G2 structure was obtained from G1-only or G2-only optimized initial structures.

Structural differences within the G1 region between the G1 only and the entire G2 (G1+G2 regions) relaxed optimizations are summarized in Table 7. With the exception of the orientation of the s–HO groups there are no significant differences between the internal coordinates for these optimizations. This suggests that orientation of the s–HO groups has a negligible effect on the internal structure of the TO-layer of exfoliated kaolinite.

Structural optimization of the G3 model

The G3 model of ~4 nm diameter with more than 700 atoms and about 3,000 valence electrons was optimized for three different regions of G1, G2 including G1, G3 including G1 and G2 regions and without any constraints for the periphery and the outer sphere of neutralizing counter-ions. The G1 optimization was completed rapidly, while the remainder showed some minor convergence issues as the number of optimized atomic positions increased. Further convergence issues emerged only if the rings of the G2-only and G3-only regions were optimized initially and then the additional regions were allowed to relax. The computational cost for structural optimizations of the entire G3 model without any constraints became comparable to those costs for a G1 model with DFT methods (Table 1), but this was mainly due to the low parallelization efficiency of the current MOPAC2012 implementation and the large number of optimization cycles needed to achieve stable stationary structures. The optimized structures are shown in Fig. 8. Figure 8A compares the side views of various optimized regions for the G3 model showing how well the structural integrity of the model remained. The i–HO groups point upward as seen for the G2 and G1 models; however, the s–HO group arrangement remains similar to those calculated for periodic models. An important observation is the pronounced curvature of the unconstrained optimized G3 model, which is consistent with the trends observed experimentally, as exfoliated kaolinite tends to adopt tubular or ribbon-like morphology due to its crystallinity. The formation of nanotubes and ribbons can be promoted by both the methoxylation of the O-sheet and the swelling of the layers in methanol suspension (Kuroda et al., 2011). The change in shape is the result of the considerable size differences between the [AlO6] and the [SiO4] coordination environments, with the latter having considerably shorter Si–O distances than the Al–O and thus a more contracted coordination environment. In the crystalline kaolinite, the strong interlayer H-bonding relaxes the structural differences between the [AlO6] and the [SiO4] moieties by rendering the Si–O (bridging) bonds longer by ~0.02–0.03 Å and the Al–O (surface) shorter by 0.05 Å, compared to the exfoliated counterparts. This morphology difference between the exfoliated and original platy phases should already be detectable experimentally for the ~5 nm particles by high-resolution transmission electron microscopy or atomic force microscopy.

It was unexpected that as the model size was increased to G3 the variability of the s–HO group arrangements virtually disappeared (Table 6 bottom section) and all calculations provide results that are comparable to the reference B3LYP/TZVP level. The similarity of the G1 only regions of the G3 models in Table 7 emphasizes that the molecular cluster model for the G1 region converges at the G2 level and consideration of G3 will not provide superior geometric structural results.

The electrostatic contour surfaces (Fig. 8B) calculated from the PM7/COSMO atomic point charges for the G2 region optimized structure of the G3 model, demonstrate how the relaxation of the s–HO groups can influence the surface reactivity. The isotropic electrophilic O-sheet of a kaolinite (as seen at the periphery of the top surface in Fig. 8B, marked blue) will disappear for the exfoliated kaolinite particles. Instead, stripes of nucleophilic regions appear as the α–s–HO groups above the i–HO groups fold in to become parallel with the surface (white and light pink coloured regions in Fig. 8B). This is a remarkable large-scale change that is expected to contribute to an amphoteric nucleophilic/electrophilic property of the surface with structural order as has also been shown for the G2 model in Fig. 7B. The T-sheet of the exfoliated kaolinite layer shows a homogeneous, nucleophilic environment (Fig. 7B bottom and right hand side insets) with large enough cavities within the Si-honeycombs that can accommodate cations of inorganic or CH3 groups of electrophilic organic intercalating, delaminating, or exfoliating agents. According to interpretations of previous DFT calculations (Táborosi et al., 2014), signs of anisotropy also appear for the T-sheet due to the strong and directed H-bonding interactions between the straightened i–HO groups and the adjacent b–O2− ions.

Summary and Conclusions

In this work, we presented evidence for the viability of molecular cluster models as accurate and economic alternatives to two-dimensional periodic boundary condition models for exfoliated kaolinite. In a systematic model-building approach, pitfalls in modelling were identified and work-arounds were provided. Using coordination chemistry-based principles for model construction and density functional theory for calculating chemical interactions with sufficiently large basis sets, chemically reasonable atomic positional coordinates can be obtained for up to 100-atom models (neutralized First Generation model, G1+) for exfoliated kaolinite. This can be achieved for about an order of magnitude less computational cost for the cluster models than for the corresponding periodic models at the same level. It is also remarkable that the HF method provides comparable structural results to the DFT calculations. Earlier it had been emphasized that there is a need for a saturated basis set of at least triple-ζ quality together with a polarization function (def2TZVP). In this work results were discussed that suggest the reliability of a smaller, and thus computationally more efficient, SVP basis set. This opens up the possibility for the ab initio DFT treatment of G2 size computational models.

Periodic and cluster model calculations at the DFT level were used to evaluate the accuracy of PM6 and PM7 Hamiltonians. PM7 was superior to PM6 even when using dispersion and H-bonding corrections for the latter. Polarizable continuum modelling according to COSMO formalism improves the description of the surface groups, which is expected to be critical for understanding electrophilic and nucleophilic reactivity of the surface towards adsorbents and intercalating reagents. The considerably lower computational cost in using PM7 in the GGA DFT calculations allowed the evaluation of Generation 2 and 3 models that correspond to a radial expansion of the G1 model with Al- and Si-honeycombs. Despite some quantitative differences between the semi-empirical G1 model and the reference DFT structure, the G2 and G3 models provide a chemically reasonable description for the surfaces of the exfoliated kaolinite. From comparison of the optimized G1 region from G2 and G3 models it was found that the G2 model is already converged and using a G3 model does not result in any different geometric structures for the G1 region. The G3 model, with ~4 nm diameter, shows a dome-shaped curvature, which parallels the morphologies observed experimentally.

The G1 model was reasonable for DFT modelling where the covalent and ionic interactions are described adequately; however, for semi-empirical models at least a G2 size model is required for modelling central Al- and/or Si-honeycombs, as adsorption sites. A notable advantage of using a G2 model is the possibility of employing proton termination for the peripheral dangling oxide ions attached to Al3+- or Si4+-ions, or bridging between them. The completion of a series of modelling studies would be necessary for understanding the behaviour of surface hydroxides at DFT level and their response to the presence of an explicit aqueous environment in both DFT and PM7 calculations. It would also be necessary to understand the importance and the magnitude of the curvature of exfoliated kaolinite, which may exclude the possibility of using slab-like periodic boundary conditions and the use of peripheral constraints in the molecular cluster models. Upon completion of these studies adsorption and intercalation processes might be simulated at quantum chemical level with explicit solvent environments for the exfoliated kaolinite sheet for the structure of which little information is available.

This research was supported by the European Union and the State of Hungary, co-financed by the European Social Fund in the framework of TÁMOP-4.2.2.A-11/1/KONV-2012-0071 (AT) and TÁMOP 4.2.4.A/2-11-1-2012-0001 ‘National Excellence Program’ (RKS). We acknowledge the use of the computational facilities at Montana State University, Bozeman, MT, USA, and the National Information Infrastructure Development (NIIF) Program of Hungary.

Electronic supporting information containing detailed model creation steps, molecular structures, formatted checkpoint files, electronic spreadsheet with comparison of experimental and various calculated data are available at http://computational.chemistry.montana.edu/clay3.