Simulating the Dehydroxylation Reaction in ... - ACS Publications

Article pubs.acs.org/JPCC

Simulating the Dehydroxylation Reaction in Smectite Models by Car−Parrinello-like−Born−Oppenheimer Molecular Dynamics and Metadynamics Daniel Muñoz-Santiburcio,*,†,‡ Alfonso Hernández-Laguna,† and Claro Ignacio Sainz-Díaz† †

Instituto Andaluz de Ciencias de la Tierra (CSIC-UGR), Avenida de las Palmeras 4, 18100 Armilla, Granada, Spain Lehrstuhl für Theoretische Chemie, Ruhr-Universität Bochum, 44780 Bochum, Germany

‡

ABSTRACT: The first step of a dehydroxylation reaction of a medium-charged dioctahedral smectite model is studied by means of ab initio molecular dynamics simulations. Combining the Car−Parrinello-like−Born−Oppenheimer molecular dynamics method with metadynamics, we find two previously proposed mechanisms (dubbed on-site and cross). We confirm the existence of an intermediate state where the tetrahedral structure is deformed, a feature that is also present in the final (semidehydroxylate) state of this reaction. This shows that the thermal transformations of phyllosilicates are more complex than often assumed, which will be relevant when elucidating the mechanisms of other transformations of these materials.

■

formed in a trans-vacant structure.4 Such an event indeed implies extreme distortions of the reactant structure; however, the proposed mechanisms for the dehydroxylation in phyllosilicates do not offer a plausible explanation for this transformation. In previous works,5,10−14 we carried out molecular dynamics (MD) simulations of the dehydroxylation reaction in different trans-vacant phyllosilicate models, describing the two alternative mechanisms (cross and on-site) proposed in the literature, as well as the corresponding back reaction (rehydroxylation) and also the subsequent water release after the formation of the product water molecule. In one of these studies,13 we described a structural distortion taking place during the dehydroxylation through the cross mechanism consisting in the breaking or weakening of a basal Si−O bond. This distortion was intimately related to the reaction itself, since it makes it possible the stabilization of the excess charge on the oxygen left behind by the migrating proton, thus making it easier for the proton to avoid being recaptured by the oxygen. Since such a structural aspect of the dehydroxylation reaction was never pointed out in other experimental or theoretical studies (even though other works point to T-sheet deformation in different processes15), several calculations were made to ensure that this distortion was indeed a realistic aspect of the dehydroxylation reaction in that model. Now, we are faced with an immediate question: was that a feature of the specif ic model employed in that study, or is it a general aspect of dehydroxylation reactions in all smectites? In that previous work, the employed model presented a high

INTRODUCTION Phyllosilicate minerals have been widely studied and now are among the best known minerals in nature. However, many of the chemical reactions in which they are involved are still poorly understood, and the precise reaction mechanisms are still a matter of debate. The low degree of crystallinity, the highly disordered cationic distribution and layer stacking sequence, and the high dispersion in particle size make the correct interpretation of the experiments difficult.1 Such is the case of the high-temperature transformations of phyllosilicates such as the dehydroxylation and rehydroxylation reactions. Smectites are widely used as part of barriers in toxic and nuclear waste deposits, and these processes may compromise their stability.2 A thermal treatment at T > 300 °C produces the reaction of the structural OH groups yielding water molecules and leaving pentacoordinated Al cations, i.e. a dehydroxylation process. The temperature range of this reaction is quite broad, from 300 to 800 °C.3 This has usually been explained in terms of the different natures of the cations linked to these OH groups, the local environment of each OH group, the cation distribution, and the polymorph type.4,5 Some authors have indicated the possibility of intermediate formation,6,7 while others have reported the importance of the particle or crystal domain size and do not accept the existence of intermediates.8,9 In these reactions, severe structural distortions must occur, but the structures of intermediates and/or transition states are not well described yet. For instance, it is known that of the two possible polymorphs of dioctahedral smectites (cis- and trans-vacant, depending on the relative positions of the OH groups with respect to the vacant sites), a cis-vacant smectite subject to successive dehydroxylation−rehydroxylation reactions is trans© XXXX American Chemical Society

Received: October 17, 2016 Revised: November 27, 2016 Published: November 29, 2016 A

DOI: 10.1021/acs.jpcc.6b10436 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 1. Smectite model. The 2 × 2 × 1 supercell is indicated by the green lines. (left) Section parallel to (010). (right) View parallel to (001). Color code: H (white), O (red), Si (yellow), Al (pink), Mg (cyan), and Na (blue).

studies involving quite different systems,20−23 the number of examples in which we can find it in combination with metadynamics is still scarce.24,25 A final technical question regards the specific collective variables (CVs) used during the metadynamics simulation. While the CVs used in ref 13 were enough to obtain the desired reactions, the resolution of the obtained FES was not perfect, making it hard to determine the position of some of the transition states. Hence, it could be useful to test new sets of CVs with the aim of obtaining a better defined FES in which the different states of the system can be easily identified. In summary, the purpose of the present study is threefold: (i) analyze the influence of the cationic substitutions on the mechanistic details such as the local deformation of the tetrahedral sheet earlier described; (ii) test new sets of CVs in order to get a better description of the underlying FES; and (iii) set up and test the CP-like−BO method in our smectite model in combination with metadynamics. The latter point is particularly important, since once we have tuned and tested this method for our particular system, it could be easily used for studying a great number of problems in similar systems, such as adsorption of organic molecules in layered silicates, phyllosilicate hydrolysis, or phase changes such as the smectite−illite transition.

number of isomorphic substitutions in the octahedral sheet, where 50% of the Al3+ ions of the corresponding nonsubstituted mineral (pyrophyllite) were replaced by Mg2+, compensating the excess negative charge with the presence of Na+ cations in the interlayer space. In the present work, we study a different model in which the number of cationic substitutions is half of that in the previous study, in order to determine whether this factor influences the details of the different mechanisms such as the distortion of the tetrahedral sheet. In this study, we also aim to improve some technical aspects of our previous works. The Born−Oppenheimer MD simulations carried out in ref 13 were computationally quite expensive. Even though the metadynamics technique16 was useful to accelerate the reactions, overcoming barriers of several tenths of kBT in a feasible amount of computing time, it is highly desirable to find more efficient methods and algorithms for this kind of study. A good choice when running the simulations in highly or massively parallel computers can be the “multiple walkers” implementation of metadynamics,17 in which several replicas of the system are used to sample the free energy surface (FES) at the same time, thus reducing the computational time for the exploration of the FES by a factor equal to the number of replicas used. On the other hand, it is also possible to achieve a high efficiency if we abandon the “standard” Born−Oppenheimer MD approach (BOMD) in favor of a novel scheme such as the “Car-Parrinello-like−BornOppenheimer” (CP-like−BO) approach proposed by Kühne et al.18,19 This method circumvents the costly full self-consistent field (SCF) minimization of the wave function at each MD step by propagating the wave function with a predictor−corrector scheme, keeping a minimal deviation from the BO surface. In complex systems, a speedup of 1−2 orders of magnitude can be achieved,18,19 although a careful fine-tuning of the several parameters involved is required. It is noticeable that, while the CP-like−BO method has now been successfully used in several

■

COMPUTATIONAL APPROACH Our Na−smectite model consists of a 2 × 2 × 1 supercell, where the unit cell has the formula Na(MgAl3)Si8O20(OH)4 (thus, it is formally a montmorillonite model). The Mg2+ cation substitutions were placed in a maximally dispersed distribution as a low energy state according to previous work26 (Figure 1). The electronic structure of the system was treated at the density functional theory (DFT) level with the GPW method27 as implemented in CP2K/Quickstep.28 The PBE exchange− correlation functional29,30 along with GTH pseudopotentials,31,32 a DZVP basis set,33 and a 700 Ry density cutoff B


Article


Figure 2. (left) Instantaneous temperature Tinst and accumulated mean temperature ⟨T⟩ during the regular MD equilibration run. (right) Kinetic energy distribution of the system during the equilibration run and fit to a Gaussian function.

indeed yielding stable propagation at the desired temperature and that the NVT ensemble was being correctly sampled (Figure 2). For accelerating the reactions and mapping the FES, we employ metadynamics16,36 (MTD) using as CVs interatomic distance differences or coordination numbers (cn), defined as

were employed. It is worth noting that the plain PBE functional is quite robust and yields reliable results for a broad range of systems, and particularly for hydrogen-bonded systems the inclusion of either D2 or D3 dispersion corrections only worsens the results.34 In consequence, given the nature of the dehydroxylation reaction, dispersion corrections were not included in our simulations. After optimizing the cell parameters of our supercell at 0 K, these were kept fixed for the rest of the calculations at the values a = 10.510 Å, b = 18.078 Å, c = 10.236 Å, α = 90.711°, β = 107.412°, and γ = 89.967°. For the ab initio molecular dynamics (AIMD) simulations35 we use the CP-like−BO scheme,18,19 taking a 0.4 fs time step. With this method, canonical sampling is established according to the Langevin-type equation

cn[A−B] =

1 NA

NA

NB

∑∑ i=1 j=1

1 − (R ij/R 0)m 1 − (R ij/R 0)n

depending on the simulation. We note that the model setup for both the “cross” and “on-site” simulations is exactly the same, and they only differ in the employed CVs. For the sake of clarity, we will employ the following notation (see also Figure 1): the specific migrating hydrogen will be always referred to as “H” (note that it is actually migrating as a proton but we omit the “+” superscript for consistency with the notation in previous works), while the rest of the hydrogens in the system are never mentioned. The oxygen to which this H is initially bonded is labeled Oi; the oxygen atom across the octahedral vacancy with respect to Oi is Ocr, and the oxygen atom linked to the same pair of octahedral cations as Oi is labeled Oos. The two apical oxygens in front of the OiH group are labeled as Oaap, and the other two apical oxygens in that same octahedral vacancy are Obap. Note that Oaap and Oi are in the same plane parallel to ab, while Obap are in the same plane as Ocr and Oos. Finally, basal O atoms in the tetrahedral sheet are labeled as Ob.

MI R̈ I = FPC − γL Ṙ I + ΞI

where ΞI is a random-noise term defined by ⟨ΞI (0) ΞI (t )⟩ = 6(γD + γL)MI kBTδ(t )

While in regular Born−Oppenheimer MD the forces acting on the nuclei are obtained after full minimization of the wave function, here the forces FPC are obtained after applying a predictor−corrector scheme to the wave functions obtained in previous steps. After several tests looking for a good balance between accuracy and efficiency, the extrapolation order was set as 7, and the number of corrector steps was adapted in each MD step by defining a threshold of 2.0 × 10−5 au (note that, in practice, the corrector steps are applied following a SCF-like procedure, and in the limit of a tight threshold a regular BOMD scheme is obtained). In our present model, this is enough to keep the number of corrector steps between 1 and 2 in all the simulations, which reduces the computing time in more than 50% with respect to a typical BOMD run. The regular Langevin parameter γL was fixed at 0.001 fs−1 and the γD parameter was adjusted on-the-fly during the preequilibration of the system until stabilizing the temperature at the desired 600 K. γD was then fixed at 5 × 10−7 fs−1 and the system was further equilibrated for 25 ps. This unusually long (compared with typical AIMD studies) equilibration time allowed us to collect enough statistics to ensure that the CP-like−BO method was

■

RESULTS AND DISCUSSION In this model there are two kinds of OH groups depending on the pair of octahedral cations to which they are bonded: Al− OH−Al and Al−OH−Mg. The behavior of Al−OH−Al groups has been extensively studied in previous works with a pyrophyllite model.10−12 In consequence, this work is focused on the dehydroxylation reaction involving only Al−OH−Mg groups. Cross Mechanism Simulation. The employed CVs for this simulation were the difference between the interatomic distances CV1 = d(H−Oi) − d(H−Ocr) and the coordination number of the H with the apical oxygens, CV2 = cn[H−Oap], C


Article

The Journal of Physical Chemistry C where we set m = 12, n = 24, and R0 = 1.35 Å (see Computational Approach). With the high exponents and the low cutoff radius, we define a sharp function with the aim of obtaining a clear distinction between the free energy wells of the intermediate and the reactant/product states. The reaction was successfully obtained, with a mechanism similar to the one observed in earlier studies.13 The H migrates to one of the apical oxygens Oaap, and from there it migrates to the receiving OcrH group with the subsequent formation of the water molecule. The stability of this product was checked by performing 5 ps of regular MD. As it was also observed in earlier studies,13 the just-formed water molecule, despite being rather mobile, does not leave the ditrigonal cavity but keeps interacting with the octahedral cations (via the lone electron pairs of the oxygen, which in turn prevents the octahedral cations from adopting the ideal 5-fold coordination geometry) and forming hydrogen bonds with the two Obap atoms, in both the MTD and regular MD simulations at 600 K (Figure 3, left).

Figure 5. Profiles of the minimum free energy paths for the cross and on-site mechanism simulations. Note that the free energy value of the product in the cross simulation is an upper limit for its actual value.

to the reactant state (i.e., before the free energy well of the product state was fully compensated by the history-dependent potential added during the MTD run), so the free energy value obtained for the product state is an upper limit for its actual value. The free energy well corresponding to the intermediate is quite shallow, and the barriers for either the recovery of the reactant or the formation of the product are 7.6 and 8.2 kcal/ mol respectively, which is on the order of only 6−7 kBT at 600 K, hence the low stability of this intermediate. This behavior can be understood by exploring the structural features of the system and comparing them to the highly charged smectite model studied in ref 13. In that system a deformation of the Tsheet was observed, consisting in the breaking of a Si−O basal bond and the formation of another bond between this Si and the Oi. The basal O atom that was undercoordinated due to the just-broken Si−O bond (labeled as Ob in Figure 1) was stabilized thanks to the coordination of the Na+ cations in the interlayer space. In that case, the Si−Oi bond just formed had a distance of 1.7 Å, on the order of all the rest of the Si−O bond distances in the system. However, in our present system a similar deformation of the T-sheet is also observed (both during the MTD run and also in the regular MD run of the

Figure 3. (left) Product for the cross mechanism simulation. (right) Intermediate for the cross mechanism simulation, optimized at 0 K.

Regarding the possible reaction intermediate in which the H is coordinated to the apical oxygen, it can be confirmed via geometry optimization at 0 K that such structure is a minimum in the PES (Figure 3, right). However, a regular MD run at 600 K (starting from a intermediate configuration extracted during the MTD simulation) reverts to the reactant structure in less than 2 ps. The FES reconstructed from the MTD simulation (Figure 4) allows the calculation of the minimum free energy path (Figure 5), observing an overall free energy barrier of 66.1 kcal/mol. However, for efficiency reasons the MTD simulation was stopped before observing the recrossing from the product

Figure 4. FES for the cross mechanism simulation. D


Article


Figure 6. (left) FES for the on-site mechanism simulation. (right) Product for the on-site mechanism.

transferred back and forth to the Obap, and eventually filling up the reactant basin in the FES and observing the H transfer from Oi to Oos. The mechanism is exactly like the one observed in our study with a highly charged smectite model:13 first, the OiH group is reoriented, crossing the plane formed by Oi, Oos, and the Al−Mg octahedral cations, so that the angle between OiH and OosH bond vectors is roughly 0. Then, the H is readily transferred to the Oos, forming the on-site product. Its stability was checked by performing a 5 ps regular MD run. Also, this time the H2O is never seen to leave the cavity in both the MTD and regular MD simulations. Regarding the deformation of the T-sheet, in this case the T-sheet is essentially unchanged. In the MD run of the product, the basal Si−O distances remain the same as in the reactant, except for a transient state in which the same deformation as the one observed in the cross mechanism simulation takes place for a brief interval of time (less than 1 ps). This is also consistent with the results observed in the highly charged smectite model.13 Regarding the reconstructed FES (Figure 6) and the minimum free energy paths (Figure 5), this simulation provides an estimate of the free energy barriers for both the cross and onsite mechanisms. In this case, the simulation was prolonged until observing the recrossing from the on-site product to the reactant state. The activation free energy for the dehydroxylation by the on-site mechanism is around 72 kcal/mol, so the on-site is less favored than the cross mechanism. Here, the employed CVs are not able to resolve the complexity of the cross mechanism as in the previous simulation, so we obtain only one TS, being the activation free energy for the cross mechanism around 61 kcal/mol. Taking into account the error derived from the neglected CVs in this last case, we can consider this value essentially in agreement with the values obtained in the cross mechanism simulation. About the free energy differences between the reactant and the products, it is clear that the cross product is preferred over the on-site product. This is also consistent with our previous study.13 It must be noted that both products are not equivalent: in the on-site product, the fully dehydroxylated mineral can be obtained by repetition of exactly the same mechanism on the rest of octahedral cation pairs. On the other hand, the cross product requires subsequent reaction steps necessarily different from

product), but it is not so pronounced as that observed in ref 13. In this case (Figure 3, right), one of the Si−Ob bonds is elongated and the Oi coordinates with the Si, but neither the just-formed Si−Oi bond stabilizes at the typical value of 1.7 Å nor the elongated Si−Ob bond is completely broken. The Si is thus pentacoordinated at the center of a trigonal bipyramid in which the Si−Oap and two Si−O bonds remain roughly in the same plane, with Si−O bond distances around 1.7 Å, while both the Si−Oi and Si−Ob bonds present a mean distance of 1.85 Å. This fact explains the low stability of the intermediate: the Si− Oi bond is not as strong as it was in the highly charged smectite model, thus the excess negative charge in the Oi is not completely neutralized and consequently it is much easier for the Oi to recapture the H while it is coordinated to the Oaap. On-Site Mechanism Simulation. For the on-site mechanism, we employed coordination numbers as CVs: CV1 = cn[Oi−H] and CV2 = cn[Oos−H]. In this case, we used R0 = 1.5 Å and the same exponents as in ref 13, m = 6 and n = 12. Interestingly, the migration of the H to the Ocr is also observed during this simulation. At the beginning of the simulation, there are several deprotonations/reprotonations of the Oi, where the H migrates back and forth from the Oi to the Oaap. However, after some of these reprotonations, the H migrates from Oi to an Obap, and from there to Ocr, then obtaining the product via the cross mechanism (the deformation of the tetrahedral sheet was also observed). The simulation was allowed to continue, and the water molecule did not leave the cavity as in the previous case. Once the free energy basin belonging to this configuration (which is exactly the same as the product obtained in the previous simulation) was filled up, the H migrated back to the Oi recovering the initial structure. It is worth noting that this last migration happened in a direct fashion from the Ocr to Oi without the H being coordinated to any of the apical oxygens. This is possible thanks to the higher mobility of the H2O compared to the structural OH groups, which makes less likely the direct proton transfer from the OiH group to Ocr, even though prior studies in pyrophyllite were able to reproduce such a mechanism with a direct proton migration.11 After the recovery of the reactant configuration, the MTD simulation continued, observing more deprotonations/reprotonations of the Oi, this time being the H E


Article

The Journal of Physical Chemistry C this first step in order to obtain the fully dehydroxylated mineral. Further Discussion. An immediate conclusion of this study is that the cross mechanism for proton migration is preferred over the on-site mechanism. This is in agreement with our previous study.13 However, the free energy barriers reported in the present work are slightly higher than in that one. This may be a consequence of the lower charge in the interlayer space in this model in comparison with the previous model. As we discussed previously, once the proton leaves the Oi, the excess negative charge on this atom can be stabilized by coordinating with a basal Si atom. For this, a Si−Ob bond must be broken or weakened, which is facilitated by stabilization of the just-formed excess charge on the Ob by the interlayer Na+ cations. Since in the present model there are fewer interlayer cations, the sequence of interactions that allows the stabilization of the Oi is diminished, thus making it more difficult for the leaving proton to escape from the attraction of the Oi. This work is focused mainly on the first step of the dehydroxylation process, i.e. the formation of the hydrated semidehydroxylate intermediate. This is the limiting step of the whole process according to previous works.14 Previous AIMD simulations on the dehydroxylation of pyrophyllite showed12,14 that the following steps of the process, namely the release of the water from the tetrahedral cavity to the interlayer space (19−24 kcal/mol), the dehydroxylation of the semidehydroxylate state (21−31 kcal/mol), and migration of the water outside of the interlayer space (15 kcal/mol) require much less activation free energy than the first step (48−50 kcal/mol). Although in this smectite model there are interlayer cations to which the water can coordinate during the water release process, this effect will most likely have a small impact on the free energy barrier for the release of the water molecule considering the high temperatures at which this reaction takes place. The dehydroxylation in phyllosilicates takes place over a wide range of temperatures. This is usually explained in terms of the different strength of the bonds between the structural OH groups and the octahedral cations. However, in the present study the migrating proton is released from a Al−OH−Mg group and received by a Al−OH−Mg group, just as in the previous study, and the difference between the systems resides not in the octahedral cations to which the reacting OH groups are bonded but in the environment of these reacting groups. Hence, it is possible that the wide range of temperatures experimentally observed for the dehydroxylation is due not only to the different strength of the bonds between the OH groups and the octahedral cations, but also to the different ability of the environment to stabilize the residual charge left behind by the migrating proton.

observed between the activation free energies in two different models where the reacting OH groups reside in equivalent pairs of octahedral cations, thus giving another plausible explanation for the wide range of temperatures for this reaction observed in experiments. This peculiar mechanistic detail of the dehydroxylation reaction in phyllosilicates had never been suggested until our previous investigation of a highly charged smectite model.13 Its confirmation in the present study demonstrates that the mechanisms of thermal transformations of phyllosilicates are more complex than what has hitherto been assumed, which may be very relevant in future investigations of other elusive mechanisms such as that of the transformation from cis- to trans-vacant smectite upon de- and rehydroxylation.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] or [email protected]. de. ORCID

Daniel Muñoz-Santiburcio: 0000-0001-9490-5975 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors are thankful to the Centro Técnico de Informática (CSIC), and the Centro de Servicios de Informática y Redes de Comunicaciones (CSIRC), Universidad de Granada, for providing the computing time. This work was supported by the Spanish MCYT and European FEDER grants CGL200502681, CGL2008-02850, CGL2008-06245-CO2-01, and FIS2013-48444-C2-2P, and Andalusian grants RNM363, RNM264, and RNM1897.

■

REFERENCES

(1) Geiger, C. A. Silicate garnet: A micro to macroscopic (re) view. Am. Mineral. 2008, 93, 360−372. (2) Pusch, R.; Kasbohm, J.; Knutsson, S.; Yang, T.; Nguyen-Thanh, L. The role of smectite clay barriers for isolating high-level radioactive waste (HLW) in shallow and deep repositories. Procedia Earth Planet. Sci. 2015, 15, 680−687. (3) Derkowski, A.; Drits, V. A.; McCarty, D. K. Nature of rehydroxylation in dioctahedral 2:1 layer clay minerals. Am. Mineral. 2012, 97, 610−629. (4) Drits, V.; Besson, G.; Muller, F. An improved model for structural transformation of heat-treated aluminous dioctahedral 2:1 layer silicates. Clays Clay Miner. 1995, 43, 718−731. (5) Molina-Montes, E.; Timón, V.; Hernández-Laguna, A.; SainzDíaz, C. I. Dehydroxylation mechanisms in Al3+/Fe3+ dioctahedral phyllosilicates by quantum mechanical methods with cluster models. Geochim. Cosmochim. Acta 2008, 72, 3929−3938. (6) Kloprogge, J. T.; Komarneni, S.; Yanagisawa, K.; Fry, R.; Frost, R. L. Infrared emission spectroscopic study of the dehydroxylation via surface silanol groups of synthetic and natural beidellite. J. Colloid Interface Sci. 1999, 212, 562−569. (7) Zhang, M.; Redfern, S. A.; Salje, E. K.; Carpenter, M. A.; Hayward, C. L. Thermal behavior of vibrational phonons and hydroxyls of muscovite in dehydroxylation: In situ high-temperature infrared spectroscopic investigations. Am. Mineral. 2010, 95, 1444− 1457. (8) Drits, V. A.; Derkowski, A.; McCarty, D. K. New insight into the structural transformation of partially dehydroxylated pyrophyllite. Am. Mineral. 2011, 96, 153−171.

■

CONCLUSIONS The dehydroxylation of a medium-charge smectite model has been studied by means of ab initio simulations. The cross mechanism was found to be the most probable one; however, the activation free energy for the on-site mechanism is close enough that both can happen at high temperature. We confirmed the involvement of the tetrahedral layer, which plays an important role in the reaction by establishing a weak bond between a basal Si and the O that releases the proton, thus stabilizing the residual negative charge left by the migrating proton. This deformation, which is in turn made possible by the interlayer cations that stabilize the basal O atom that is left undercoordinated, is responsible for the differences F


Article

The Journal of Physical Chemistry C (9) Drits, V. A.; Derkowski, A.; McCarty, D. K. Kinetics of partial dehydroxylation in dioctahedral 2:1 layer clay minerals. Am. Mineral. 2012, 97, 930−950. (10) Molina-Montes, E.; Donadio, D.; Hernández-Laguna, A.; SainzDíaz, C. I. DFT research on the dehydroxylation reaction of pyrophyllite 2. Characterization of reactants, intermediates, and transition states along the reaction path. J. Phys. Chem. A 2008, 112, 6373−6383. (11) Molina-Montes, E.; Donadio, D.; Hernández-Laguna, A.; SainzDíaz, C. I.; Parrinello, M. DFT research on the dehydroxylation reaction of pyrophyllite 1. First-principle molecular dynamics simulations. J. Phys. Chem. B 2008, 112, 7051−7060. (12) Molina-Montes, E.; Donadio, D.; Hernández-Laguna, A.; SainzDíaz, C. I. Exploring the rehydroxylation reaction of pyrophyllite by ab initio molecular dynamics. J. Phys. Chem. B 2010, 114, 7593−7601. (13) Muñoz Santiburcio, D.; Kosa, M.; Hernández-Laguna, A.; SainzDíaz, C. I.; Parrinello, M. Ab initio molecular dynamics study of the dehydroxylation reaction in a smectite model. J. Phys. Chem. C 2012, 116, 12203−12211. (14) Molina-Montes, E.; Donadio, D.; Hernández-Laguna, A.; Parrinello, M.; Sainz-Díaz, C. I. Water release from pyrophyllite during the dehydroxylation process explored by quantum mechanical simulations. J. Phys. Chem. C 2013, 117, 7526−7532. ́ (15) Szczerba, M.; Derkowski, A.; Kalinichev, A. G.; Srodoń , J. Molecular modeling of the effects of 40Ar recoil in Illite particles on their K-Ar isotope dating. Geochim. Cosmochim. Acta 2015, 159, 162− 176. (16) Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562−12566. (17) Raiteri, P.; Laio, A.; Gervasio, F. L.; Micheletti, C.; Parrinello, M. Efficient reconstruction of complex free energy landscapes by multiple walkers metadynamics. J. Phys. Chem. B 2006, 110, 3533−3539. (18) Kühne, T. D.; Krack, M.; Mohamed, F. R.; Parrinello, M. Efficient and Accurate Car-Parrinello-like Approach to BornOppenheimer Molecular Dynamics. Phys. Rev. Lett. 2007, 98, 066401. (19) Kühne, T. D. Second generation Car-Parrinello molecular dynamics. WIREs Comput. Mol. Sci. 2014, 4, 391−406. (20) Caravati, S.; Bernasconi, M.; Kühne, T.; Krack, M.; Parrinello, M. Unravelling the mechanism of pressure induced amorphization of phase change materials. Phys. Rev. Lett. 2009, 102, 205502. (21) Cucinotta, C. S.; Miceli, G.; Raiteri, P.; Krack, M.; Kühne, T. D.; Bernasconi, M.; Parrinello, M. Superionic conduction in substoichiometric LiAl alloy: an ab initio study. Phys. Rev. Lett. 2009, 103, 125901. (22) Kühne, T. D.; Khaliullin, R. Z. Electronic signature of the instantaneous asymmetry in the first coordination shell of liquid water. Nat. Commun. 2013, 4, 1450. (23) Raty, J. Y.; Zhang, W.; Luckas, J.; Chen, C.; Mazzarello, R.; Bichara, C.; Wuttig, M. Aging mechanisms in amorphous phasechange materials. Nat. Commun. 2015, 6, 7467. (24) Caravati, S.; Sosso, G. C.; Bernasconi, M.; Parrinello, M. Density functional simulations of hexagonal Ge2Sb2Te5 at high pressure. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 094117. (25) Ronneberger, I.; Zhang, W.; Eshet, H.; Mazzarello, R. Crystallization properties of the Ge2Sb2Te5 phase-change compound from advanced simulations. Adv. Funct. Mater. 2015, 25, 6407−6413. (26) Sainz-Díaz, C. I.; Palin, E. J.; Dove, M. T.; Hernández-Laguna, A. Monte Carlo simulations of ordering of Al, Fe, and Mg cations in the octahedral sheet of smectites and illites. Am. Mineral. 2003, 88, 1033−1045. (27) Lippert, G.; Hutter, J.; Parrinello, M. A hybrid Gaussian and plane wave density functional scheme. Mol. Phys. 1997, 92, 477−488. (28) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. Comput. Phys. Commun. 2005, 167, 103−128. (29) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation made simple. Phys. Rev. Lett. 1996, 77, 3865−3868.

(30) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation made simple [Phys. Rev. Lett. 77, 3865 (1996)]. Phys. Rev. Lett. 1997, 78, 1396−1396. (31) Goedecker, S.; Teter, M.; Hutter, J. Separable dual-space Gaussian pseudopotentials. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 1703−1710. (32) Krack, M. Pseudopotentials for H to Kr optimized for gradientcorrected exchange-correlation functionals. Theor. Chem. Acc. 2005, 114, 145−152. (33) VandeVondele, J.; Hutter, J. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. J. Chem. Phys. 2007, 127, 114105. (34) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (35) Marx, D.; Hutter, J. Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods; Cambridge University Press: 2009. (36) Laio, A.; Gervasio, F. L. Metadynamics: a method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science. Rep. Prog. Phys. 2008, 71, 126601.

G


Simulating the Dehydroxylation Reaction in ... - ACS Publications

Recommend Documents