Samuel Alexander
Jobbins
a,
Salah Eddine
Boulfelfel
b and
Stefano
Leoni
*a
aSchool of Chemistry, Cardiff University, Cardiff, CF10 3AT, Wales, UK. E-mail: leonis@cf.ac.uk
bGeorgia Institute of Technology, School of Chemical and Biomolecular Engineering, Atlanta, GA 30332-0100, USA
First published on 4th April 2018
We introduce a novel scheme for the mechanistic investigation of solid–solid phase transitions, which we dub ‘metashooting’. Combining transition path sampling, molecular dynamics and metadynamics, this scheme allows for both a complete mechanistic analysis and detailed mapping of the free energy surface. This is illustrated by performing metashooting calculations on the pressure-induced B4/B3 → B1 phase transition in ZnO. The resulting free energy map helps to clarify the role of intermediate configurations along this activated process and the competition between different mechanistic regimes with superior accuracy. We argue that metashooting can be efficiently applied to a broader class of activated processes.
Protein folding, nucleation, crystallisation, phase transitions and chemical reactions are all important phenomena which are characterised by activated steps. Their detailed mechanistic investigation would benefit considerably from precise mapping of the underlying free energy landscape. However, despite improving numerical algorithms and increasing computational power, their systematic investigation remains extremely demanding and often completely impractical. Methods to explore multi-minima energy landscapes exist4,5 and recently there has been a pronounced effort to develop novel approaches to overcome this time-scale problem.3,6 Historically, the use of an external bias aimed at enhancing the frequency of rare events in MD simulations has been realised in methods such as umbrella sampling (US).7 This is achieved by systematically adding a bias along well-chosen coordinates, in order to sample only small but overlapping regions from which a free energy can be reconstructed. Another bias-oriented method, metadynamics (metaD),3 exploits the concept of collective variables (CVs), which span a lower dimensional coordinate space. Within the space of these CVs, a history-dependent bias is deposited, which rapidly drives the system away from its initial probable state. The definition of a good set of CVs is therefore key to the success of the method.8 Within a CV space, standard MD would linger around a minimum close to the initial state corresponding to some value of the CV. Overcoming free energy barriers by metadynamics inevitably leads to broader exploration of the CV values by crossing over to different states within a finite computational time.
Metadynamics takes place within the space of one or a set of CVs, which must be known at the beginning of the calculation. A valid set of collective variables optimally covers the entire space of the conformations of interest. This clearly requires that CVs be engineered with a global free energy map in mind. However, this is evidently problematic, as often the underlying free energy landscape is the desired final result of the simulation and it often cannot be deduced by intuition. Additionally, the requirement for such ‘global’ CVs directly clashes with the effort to keep their dimensionality as low as possible.
Recently, a novel approach has been introduced, which allows for a “bespoke bias” based on a variational approach.9,10 Within this approach, the CV problem is partially overcome by reinterpreting its meaning as a local descriptor of metastable states. Enhancing static fluctuations will amplify the probability of a crossover to another metastable state. Importantly, the bias can be deposited in such a way that partial filling can be achieved in all the basins of interest. While the success of the method still depends on an appropriate choice of CVs, this represents a novel general scheme for enhanced-sampling MD simulations.
To tackle the intrinsic disparity between the time scales of activated processes and molecular dynamics simulations, the transition path sampling method (TPS)6,11 was introduced. TPS is based upon the idea that, in a complex system, crossing a barrier between (meta-)stable states may happen in a multitude of different ways. In such a scenario, the potential energy surface of a system is unlikely to be represented by saddle points alone. Instead, numerous points on the potential energy surface become relevant, some of which may still be stationary points, but many others not. The implementation of the method consists of collecting such an ensemble of reactive trajectories, without previous knowledge of the details of the process (that is, no a priori information or guess about the reaction coordinates is required – such information should naturally result from the simulation approach). Within a transition ensemble, the relevance of a path is weighted by its probability. TPS collects true dynamical trajectories and can be plugged on top of different molecular dynamics schemes, within different simulation ensembles.2
TPS belongs to the class of pathway-orientated methods and allows for an exquisite level of detail in mechanistic investigations.12–18 Unlike metadynamics, TPS does not deposit any bias potential; similarly to metadynamics, it deploys a low-dimensional order parameter, which distinguishes between the basins of interest. TPS is designed to enhance statistics in intermediate regions between long-lived basins, however it is generally not a method of choice for free energy reconstruction as the main statistical weights are still contributed by (meta-)stable basins.
Here, we present a novel combination of TPS and metadynamics techniques, which we dub ‘metashooting’. This novel approach allows for reconstruction of the underlying free energy of a process from TPS trajectories. This combination is based on the propagation of CVs from TPS to metaD. While TPS controls the rapidity of CV variation (inherent to true dynamical trajectories), metadynamics deposits a bias around regions of the CV that vary comparatively slowly. In doing so, TPS keeps the path crossing probability ratios (A → B to B → A) unchanged, thus preventing unilateral over-biasing. By rapidly varying the controlling CVs in TPS, partial basin filling can be achieved while maintaining the intermediate regions bias-free for most of the calculations. Only in its final phases do the details of the intermediate regions appear.
In the following we illustrate the details of our approach by performing a complete metashooting analysis of the solid–solid reconstructive phase transition between the B3 zincblende phase and the B1 high-pressure rocksalt phase of zinc oxide (ZnO). ZnO has proven to be extremely challenging over the years,19–23 with several attempts to conclusively shed light on competitive mechanisms and intermediates. Saitta and Decremps first discussed the competitive nature of the ‘tetragonal’ and ‘hexagonal’ routes in the B4–B1 phase transition, and their work concluded that semiconductors involving d-electrons, such as ZnO, preferred the ‘tetragonal’ path.21 However, other work seems to contradict this – notably the experimental work of Liu et al., which utilised high resolution angular dispersive X-ray diffraction, which seemed to indicate that the ‘hexagonal’ pathway was the more favoured for ZnO at lower pressures.24 Second harmonic generation experiments pinpointed the role of the pressure medium (hydrostatic vs. non-hydrostatic) in promoting the formation of centrosymmetric iH or acentric tetragonal iT, respectively.22 In previous work,15 TPS was used to produce a detailed mechanistic pathway for the B4–B1 transition, and it was shown that the mechanism proceeds via a rich series of transition states and intermediates.12 Additionally, the two 5-coordinate intermediate motifs were visited and were shown to compete with one another in the mechanism. The analysis, based on dynamical trajectories, concluded that a hexagonal intermediate could indeed be visited, however it was not essential to the transformation. In fact, the reconstruction was found to propagate over the interfaces of local tetragonal motifs (Fig. 1).
Fig. 1 Two 5-coordinate intermediates (tetragonal iT (I4mm, left) and hexagonal iH (P63/mmc, right)) were shown to compete against one another in the B4–B1 transformation. Whilst iT was only found at the interface between B4 and B1, the iH structure sometimes dominated the entire system under study. Despite this, the iH phase was shown not to be crucial to the transformation.12 |
It is therefore apparent that the phase transitions between different polymorphs of ZnO are non-trivial. Indeed, much debate still rages over their precise nature and work is still being published to attempt to unravel the minutiae of these elusive mechanisms, with a particular emphasis on the B4–B1 transformation. In the following, after a summary of the implementation details, the results of metashooting on a challenging solid–solid phase transition in this material between the B3 and B1 phases are discussed.
Using only the 1st coordination shell means that there is no discrimination between the B3 and B4 configurations. As the system starts from a metastable B3 structure (see below), changes in the final regime can be expected. Furthermore, using solely the 1st coordination number as the order parameter exerts absolutely no bias on the calculation. The variation in the order parameter arises simply as a result of the application of the transition pressure to the system. Hence, monitoring the 1st coordination sphere does not force the system to evolve in any particular way over the course of the path sampling iterations. In TPS, both the overall mechanism and the trajectory endpoints are rectified.26 Therefore, the order parameter is chosen with minimal constraints in mind, in order to allow the process to express the most probable outcome. A higher-order CS would drive the system into a specific transformation regime within a biased TPS approach, in order to study specific, less probable transformation pathways.
The general procedure of the metashooting approach for a transition A → B is as follows:
(1) TPS step 1: Starting from a randomly chosen snapshot along the TPS trajectories, shoot off a new trajectory and time propagate.
(2) MetaD step 1: Once basin A is reached, deposit a bias using a metadynamics scheme (in this case, using the well-tempered approach), in order to locally partially fill the as-yet uncharacterised energy well. At this stage, a small number of metasteps are chosen, allowing the CV to vary about the characteristic value of the basin.
(3) TPS step 2: Generate another trajectory (as in step 1) and propagate into the now-biased basin A, to generate a true MD velocity distribution on the biased potential.
(4) TPS step 3: Similarly to TPS step 1, shoot off a new trajectory from the same snapshot as above and propagate, only this time backwards in time, to basin B.
(5) MetaD step 2: In basin B, apply the same metadynamics scheme as above for a number of steps, in order to partially fill this second basin of attraction.
(6) TPS step 4: Re-shoot a novel trajectory from a random snapshot of the previous trajectory. Propagate into the now-biased basin B.
Repeat the above iterations of metadynamics and TPS moves until the underlying free energy profile is fully converged.
The number of metadynamics steps for each metashooting move can be chosen such that only moderate variations around the current CV value are allowed, as ΔCVMax. Alternatively, as a better means to control bias deposition without imbalance, the TPS path acceptance probability in TPS can be chosen to stay the same, i.e. at no point can crossing over the activation energy barrier become easier for either direction A → B or B → A. This is realised within an adaptive scheme of rescaling momenta modifications (see below).
Fascinatingly, the final 4-coordinated product was not a pure zincblende (B3) structure, but a mixed B4/B3 system dominated by the wurtzite (B4) form. In all the successful trajectories, the final material consisted of only 10–20% B3, with the more stable hexagonal form forming the remainder of the structure. The order parameter based on the 1st coordination number permits both the wurtzite and zincblende configurations. Whilst the initial, unfavorable regime was B1 → B3, TPS rectified the mechanism into B1 → B4/B3, which contains predominantly hexagonal components.
The complexity of the system leads to several types of trajectories, with some degrees of (dis)similarity in the details of the reconstruction.
A representative B1 → B4/B3 trajectory is shown in Fig. 2. A nucleus of 4-connected ZnO (Fig. 2, 0.6 ps) forms within B1, followed by the formation of regions of 5-fold coordination at the interface with B1 (Fig. 2, 0.7 ps), which are subsequently outgrown by 4-fold coordinated regions (Fig. 2, 0.9–4.7 ps). Regions of 5-fold coordination (mostly iH) are found between the expanding 4-connected region and B1. The transformation into mixed B4/B3 entails the complete resorption of areas of 5-fold coordination (Fig. 2, 4.7–6.4 ps). The transformation of the remaining B1 motifs (4 ps) is associated with a sharp volume expansion, which causes the box volume to fluctuate, under the appearance of local 5-fold motifs (Fig. 2, 4.0–4.7 ps).
From a configuration of the type seen in Fig. 2 at 0.6 ps, a less severe compression leads to the growth of a 5-fold region at the expense of any 4-connected one. This accompanies the formation of a mostly 5-coordinated intermediate containing residual B1, denoted B1–iH. This configuration persists for as long as 500 ps in subsequent molecular dynamics simulations. The presence of an intermediate, which is not vital to the transformation and corresponds to a 5-membered intermediate along the transition pathway, was noted in previous work.12 The nature of this intermediate is the hexagonal iH form, whilst the tetragonal analogue iT only appears as a motif at interfaces during the reconstruction.
The B4/B3 → B1 mechanism slightly deviates from the reverse process (Fig. 3). In this trajectory type, the starting 4-coordinate material is first collectively distorted (Fig. 3, 0.0–3.2 ps, purple atoms). From the 5-coordinated areas of the iH motifs (Fig. 3, 3.9–5.2 ps), the B1 structure develops (Fig. 3, 5.2–6.4 ps). Thus, over half of the transformation time is spent slowly squeezing the system until it can transform into 5-coordinated iH regions, indicating the presence of a deep but fairly gently-sloping energy barrier separating the B4/B3 basin from the intermediate configurations. In all the configurations, the iH motifs are relatively long-lived. This is not the same as the reverse B1 → B4/B3 transition, in which the iH regions were less extensive and promptly eliminated at the interfaces (Fig. 2), with the majority of the simulation time being dominated by coexisting 4- and 6-coordinated motifs.
Fig. 3 Coordination-coloured snapshots of a representative B4/B3 → B1 trajectory. Mixed (B4/B3) wurtzite/zincblende (0.0 ps) is followed by a 4-coordinated distorted configuration (0.0–3.2 ps), purple atoms. Subsequent growth of 5-coordinated areas (green), mainly of type iH, is observed (3.2–5.2 ps), from which B1 (red) develops (5.2–6.0 ps) until the structure locks-in to 6-coordinate B1 (6.4 ps). Colours as in Fig. 2. |
Consistent with previous calculations,12 different intermediates are associated with different types of trajectories. The presence of B3 adds further structures to an already complex scenario. In previous work, low transition pressures (less than 10 GPa) were observed to favour the appearance of iH, which benefits from a less severe [0001]B4 compression. The results presented here appear to corroborate that view.
In metashooting, mechanistic details acquired from the transition path sampling procedure remain unspoiled. As the process is driven by TPS trajectories, metastable basins are filled in a balanced way – starting from the two limiting configurations, and finishing with the characterisation of all the relevant intermediate configurations. Repeated iterations of the metashooting algorithm (see Methods) enable the free energy landscapes corresponding only to the courses of the TPS-derived trajectories to be built up in a stepwise and equally distributed manner. As a result, regions around the transition states are considered only at the end of the analysis, avoiding any potential bias that may scramble the mechanistic details.
The algorithm should first fill up the two starting basins of attraction from the final TPS trajectories (B4/B3 and B1). Visiting these regions should become discouraged as a result of the bias potential. The progress of metashooting can be monitored by visualising the Gaussian deposition. Eventually, the metashooting algorithm should fully characterise the central (around the metastable basins) and peripheral regions of the CV space over the course of the trajectory. Convergence will be seen in a similar way to a standard metadynamics scheme – that is, the underlying free energy landscape should no longer significantly change with further iterations of the procedure. At this point, meaningful analysis of the intermediates and energy barriers associated with each process can be carried out.
Harnessing the free energy in this way should have a number of advantages. The landscape can be thought of as the surface upon which the trajectories travel, and the individual characteristics of the trajectories are determined by which regions of the surface are visited. Thus, its elucidation allows for a complete quantitative analysis of the underlying thermodynamics and kinetics of the process; measurements of the activation energy and changes in the free energy between different steps of the procedure, as well as the starting and finishing products, can trivially be taken from the plot. Moreover, the intermediates and transition states that are vitally important to the transformation will appear obvious along the reaction scheme. Such insight into trajectories using only path-based analyses of transformations is not trivial.
The free energy profile was constructed from the converged TPS trajectories. After 250 iterations of the metashooting procedure, the free energy had fully converged into a detailed landscape.
The first few iterations of the procedure (Fig. 4) accumulated Gaussian functions only within the two limiting basins (B4/B3 and B1). A small ‘shoulder-peak’ formed in the 4-coordinate basin as a result of oscillations between different forms of the low-coordination material, including B3, B4 and other forms.
After 50 iterations (Fig. 4b), a notable ‘shoulder-peak’ adjacent to the B1 basin was observed. After 150 iterations, this was resolved into a bona fide minimum, i-1. By iteration 200, a second minimum centred halfway between the two limiting basins started to appear (Fig. 4c), which was fully resolved by iteration 225 (Fig. 4d), denoted i-2. By 250 iterations (Fig. 5), the plot was fully converged – the energy barriers and positions of the minima remained unchanged with any further iterations of the metashooting procedure.
Turning point | 1st CN | 3rd CN |
---|---|---|
RS (B1) | 5.98 | 37.90 |
TS-A | 5.75 | 36.66 |
i-1 | 5.68 | 36.21 |
TS-B | 5.43 | 34.86 |
i-2 | 5.15 | 33.13 |
TS-C | 4.87 | 31.62 |
WZ-ZB (B4/B3) | 4.10 | 29.22 |
In i-1 a portion of the pristine B1 fragments contains 4- and 5-coordinated motifs. With respect to the way the square pattern is opened, the locally 5-coordinated motifs correspond to iH. This small region is extremely short-lived but is present in every successful trajectory sampled, thus explaining its place on the free energy profile. Interestingly, this mechanistic pattern cannot propagate as it is. The next free energy barrier TS-B is overcome by including an additional (111)B1 shearing mechanism, as seen in both the forward and reverse TPS trajectories. This configuration, which includes an enlarged seed region and the onset of shearing, corresponds to the free energy maximum TS-B. The TS-B configuration could easily be seen in the TPS trajectories, however the fact that this part of the transformation corresponds to two distinct events would have been completely missed without the metashooting analysis.
The second intermediate region, i-2, is a deep energy well placed halfway between the two limiting basins. Trajectories visiting configurations characterised by the 1st CN in the interval [5.0–5.3] and the 3rd CN within [32–34] can in principle contribute to this basin. Accordingly, three different types of trajectory can be distinguished:
(1) Successful B4/B3 → B1 trajectories, along which extended 5-coordinated regions are present (see Fig. 3, 3.9–5.6 ps).
(2) Successful B1 → B4/B3 trajectories, in which Zn or O is 4-, 5- and 6-coordinated, giving an overall average 1st coordination sphere of five (see Fig. 2, 0.7–0.9 ps).
(3) Trajectories where the system visits and lingers in the predominantly 5-coordinate intermediate B1–iH. The average time spent in this basin is >200 ps.
Representative trajectories of the three types above are mapped on the free energy surface in Fig. 7. Successful trajectories tend to keep a higher 3rd coordination number while they approach the large i-2 basin. On the other hand, a lower 3rd CN is associated with a trajectory entering an intermediate basin. Closer to the B1 basin, the trajectories are more similar, while the B4/B3 basin can be escaped in different ways. Around TS-B, mechanistic competition (assisted by box size fluctuations12,22) may result in very different trajectories.
It seems likely, therefore, that the depth of the intermediate basin (and hence the energy barriers associated with it) cannot be taken as representing a proper intermediate along a particular type of trajectory. Instead, numerous possible configurations with similar average coordination number values all add weight to the free energy at these values of the collective variables, even though not every configuration is realisable in every trajectory.
Unlike TS-A and TS-B, TS-C corresponds to a very ‘gentle’ maximum, with a very shallow gradient either side of it. This again corroborates what was seen in the path sampling investigations – the system spends a great deal of time in a 4-coordinate state, or a mixture of 4- and 5-coordinate states, as the result of a slow, steady process, which is mapped here onto a gentle free energy gradient. Like its associated intermediate basin i-2, the height of TS-C is influenced by a number of configurations, all of which contain some combination of 4-coordinated ZnO, iH intermediate and residual B1. Crossing this maximum implies a significant coordination number (order parameter) change, and corresponds to mixed motif configurations as at 4 ps in Fig. 2 for the B1 → B4/B3 direction, or at 4.7 ps for B4/B3 → B1 (Fig. 3).
The appearance of mixed B1–iH agrees with previous studies on ionic compounds undergoing the same B4/B3 → B1 transition,34–36 which may indicate the transferability of this mechanistic analysis.
Trajectory variety has shown that the system may proceed along different pathways, therefore crossing different features of the free energy plot. This will lead to different values for the total amount of energy overcome by the individual trajectories. The following quantitative analysis of the energy differences assumes measurement to and from the centre of each of the turning points, as described in Table 2, thereby indicating the energy associated with the minimum energy transition pathway that visits all maxima and minima. Every trajectory must surmount an energy barrier at least equivalent to the activation energy presented, as successful trajectories must overcome the highest energy barrier to cross to the other side, irrespective of the path taken.
Energy difference | kJ mol−1 | k B T/ZnO |
---|---|---|
ΔE1 | 20447.0 | 6.8 |
ΔE2 | 10483.8 | 3.5 |
ΔE3 | 13085.9 | 4.4 |
ΔE5 | 10376.6 | 3.5 |
ΔE6 | 1357.0 | 0.5 |
ΔE7 | 23726.3 | 7.9 |
ΔE4 | 35069.9 | 11.6 |
ΔE8 | 44586.7 | 14.9 |
It is also apparent from inspection that the energy barriers for each stage of the transformation are very different for the forward and reverse trajectories. This may go some way towards explaining why there are different pathways associated with the forward and reverse trajectories, and perhaps why there is such notable hysteresis present in such transformations.
The change in free energy is calculated as the difference between the energy values of the B1 and B4/B3 basin centres. The initial pressure was chosen based on the zero enthalpy difference between B3 and B1, which was dictated by the model used to commence the TPS cycles (see above).
The difference in free energy between the forward (B4/B3 → B1) and backward transformation is not zero at the chosen simulation pressure (p = 9.8 GPa). The inclusion of lattice dynamics37 in first principles calculations was shown to shift the (ΔH = 0) “equilibrium” pressure to a lower value. Consistent with those findings, Fig. 8 pinpoints the role of entropy, which is positive going from B4/B3 to B1. The disordered nature of B4/B3 reinforces this effect.
A typical B1 → B4/B3 trajectory was characterised by the formation of an initial seed, followed by the gradual transformation to a mixed 4-coordinate product via coexisting motifs. Trajectories could also end up trapped in a long-lived intermediate 5-coordinate basin, corresponding to a mixed rocksalt-5-coordinate hexagonal structure denoted B1–iH. The B4/B3 → B1 transition entailed a more stepwise process, which involved a transition to a 4-coordinate mixture via local motifs of hexagonal iH before transforming to the final B1 product.
A novel method denoted ‘metashooting’ was developed as a combination of TPS and well-tempered metadynamics, using the 1st and 3rd average coordination spheres as the collective variables. TPS ensured rapid variations of the CVs, and balanced filling of the basins of interest. The metashooting procedure is able to decipher the underlying free energy landscape of the transformation, and structural information about the nature of appropriate maxima and minima can be extracted from the resultant plot. In addition, the energy barriers and changes in free energy can be ascertained and the successful and failed trajectories can be mapped onto the final free energy plot, in order to gain a true understanding of the energetics and structural changes inherent within each individual trajectory.
Metashooting can accommodate additional collective variables, for example higher coordination sequences, in order to further distinguish between the phases. Using the 1st and 3rd coordination spheres was convenient in this case as it was an identical measurement to the order parameter used in the path sampling calculations. In general, the CV set used on the metashooting layer can be more detailed then the order parameter utilised in the TPS steps.
It is believed that the metashooting procedure, which combines transition path sampling and metadynamics, could set a novel paradigm for future investigations into condensed matter phase transitions. Such a scheme could easily be transferred to another system, under any level of theory, to ascertain the underlying thermodynamics and kinetics of a temperature- or pressure-induced phase transformation.
This journal is © The Royal Society of Chemistry 2018 |