Time-Disordered Crystal Structure of AlPO4 ... - ACS Publications

Aug 3, 2017 - unity at 1.5 K. The elastic intensity is maximized below 100 K and drops to a ..... Government Research Training Program (RTP) Scholarsh...
0 downloads 0 Views 3MB Size
Article pubs.acs.org/JPCC

Time-Disordered Crystal Structure of AlPO4‑5 D. L Cortie,*,†,‡ B. R. McBride,† N. Narayanan,†,‡ N. R. de Souza,‡ M. Avdeev,‡ R. A. Mole,‡ G. J. McIntyre,‡ G. J. Kearley,‡ R. Withers,† D. H. Yu,*,‡ and Y. Liu*,† †

Research School of Chemistry, Australian National University, Canberra, ACT 2601, Australia Australian Nuclear Science and Technology Organisation, New Illawarra Rd., Lucas Heights NSW 2234, Australia



S Supporting Information *

ABSTRACT: Modern ab initio calculations have become increasingly accurate for predicting the symmetry of crystal structures; however, the standard methods usually only determine the 0 K static configuration and therefore misrepresent structures where dynamics play a key role. In this work, we demonstrate a clear experimental example of this phenomenon in the AlPO4-5 molecular sieve where the average high-symmetry P6cc structure emerges from local dynamics involving corner-sharing tetrahedra. Quasielastic and inelastic neutron scattering experiments were conducted to clarify the thermally activated motions between 1.5 and 300 K. Through comparison with ab initio molecular dynamics, we explain why the theoretically predicted structure is not observed in diffraction experiments. Instead, the results indicate this material is an inorganic analogue of a plastic crystal where thermal dynamics dictate the average symmetry.



INTRODUCTION An exciting scientific era has emerged in the past decade where novel materials are increasingly predicted and optimized using virtual-reality ab initio calculations. The success of this paradigm originates from the proven capability of densityfunctional theory (DFT) to produce good approximations for the underlying laws of quantum chemistry. In parallel, elaborate methods for crystal-structure prediction have evolved from data mining1 and simulated annealing2 toward genetic algorithms and particle swarm methods,3 which collectively have produced spectacular results.4−7 This has culminated in the creation of open databases, such as AFLOW,8 which currently contains ab initio predictions for over 1.6 million inorganic compounds,9 most of which have not yet been experimentally tested. Despite these advances, there are reasons to believe that current computational methods may misrepresent certain materials and, more importantly, fail to predict a wider class of potentially useful structures. It is crucial to note that thermal effects and dynamics are usually not included in standard DFT calculations since these require some form of additional approximation, and formally the method is restricted to the zero Kelvin (0 K) electronic state. For this reason, crystal structures are typically minimized to find the states of lowest internal energy using methods such as the popular conjugate gradient technique which does not necessarily reflect real dynamic processes. Although the resulting 0 K lattice approximation is usually very successful at finding a subset of stable structures, in principle, it may fail to predict the symmetry of materials where finite-temperature dynamics play © 2017 American Chemical Society

a key role in stabilizing structural motifs. In this article, we present clear experimental confirmation of this point using the example of the molecular-sieve AlPO4-5 where low-energy rigid-unit mode (RUM) dynamics are central to stabilizing the average high-symmetry P6cc structure. For decades, AlPO4-5 has been the subject of continuous research because it has an elegant weblike topology (Figure 1 a) composed of corner-sharing tetrahedra. This structure allows for entrapment of large guest molecules within extended

Figure 1. (a) A planar section of the experimentally determined average structure of AlPO4-5 viewed along the c-axis shows a large central 12-membered pore interconnected by 6- and 4-membered channels. In this schematic, the cell is cut at Z = 0.8 to visualize the connections between the downward-pointing AlO4 and upwardpointing PO4 tetrahedra. (b) Some of the bonds connecting tetrahedral planes in the experimental average structure have short lengths and angles close to 180°. Received: June 30, 2017 Revised: August 2, 2017 Published: August 3, 2017 18762

DOI: 10.1021/acs.jpcc.7b06415 J. Phys. Chem. C 2017, 121, 18762−18770

Article

The Journal of Physical Chemistry C micropores.10−17 So-called molecular sieve materials provide a route for targeting and filtering gas-phase pollutants and offer the potential for synthesizing novel hybrid frameworks, for example, by capturing a dye molecule.14,18,19 Despite this, many aspects of the chemical structure of AlPO4-5 still seem to be anomalous. The structure consists of alternating AlO4 and PO4 tetrahedra formed into large 12-membered rings interconnected by smaller 6- and 4-membered channels (Figure 1a). The general features have been well-established experimentally by NMR,11,20 single-crystal X-ray diffraction,13,21 neutron diffraction,12,22 and electron diffraction.23 Using diffraction, which measures the time-averaged structure, several groups have reported that AlPO4-5 crystallizes in a high-symmetry space group (P6cc or a closely related group).12,13,23 However, past theoretical calculations have indicated problems with the interplanar connectivity implicit in the aforementioned experimental structures. Specifically, the apparent bonding scheme links the ring planes with Al−O−P bonds that are far too short, with an average bond angle close to 180°23,24 (Figure 1b). Klap et al. proposed static random disorder over several partially occupied oxygen sites13 whereas Ikeda et al. modeled a similar effect using large thermal displacement parameters.12 Both approaches produce the correct average picture, but from the local perspective, if the static disorder is indeed uncorrelated between tetrahedra, this would imply an unrealistic level of flexibility in the Al−O and P−O bond lengths by necessitating distortions in the linked tetrahedra. This is at variance with the other known Al−O−P molecular sieves and inorganic crystals (e.g., berlinite). Moreover, the nearly 180° Al−O−P bond angle in the average picture contradicts calculations based on density-functional theory15 and empirical classical force field methods14,24,25 which universally predict that the crystal symmetry should relax into a lower symmetry (e.g., P21) with Al−O−P bond angles of approximately 145°. One proposal is that the order−disorder transition should occur near room temperature involving the ring planes shearing with respect to one another.14 If accurate, this would entail a large reduction in the effective pore channel of the sieve, thereby strongly influencing the ability to function as a molecular filter. However, the diffraction experiments to date do not clearly indicate any such distortion, or orthorhombic peak splitting, over a wide temperature range from 5 to 400 K. There are two viable hypotheses to explain these peculiar crystallographic features: a hidden global phase transition or a special set of local dynamics which produce a shifting time series of structural configurations. The former suggestion has been favored by every past theoretical energy minimization,14,24,25 including DFT,15 with the implicit assumption that a steady-state solution exists. However, to date, the best diffraction experiments have essentially shown no clear deviation from high symmetry.13,23 One proposal is that the transition leads to the formation of small locally disordered domains which, when averaged, obey the P6cc symmetry. However, even nanoscale probes like electron diffraction did not detect any domains.23 While local probe techniques such as NMR have found clear evidence for a more complicated site symmetry surrounding the Al nuclei, the detailed interpretation also depends on whether fast dynamics cause motional averaging leading to modified hyperfine parameters and quadrupolar splitting.11 Thus, no phase transition has ever been unambiguously detected at ambient conditions. Recent diffraction work, however, has shown high pressure causes

AlPO4-5 to transition to an incommensurate phase characterized by a modulation vector q* = (0, 0, 0.37).16 This form of modulation is virtually identical to the ambient pressure phase of the closely related SSZ24 (SiO2) compound which has the same ring topology.26 The implication is that this symmetrylowering transition is prevented at ambient pressure in AlPO4-5 by some unidentified factor. A strong candidate to explain the missing transition, as first proposed by Liu and Withers et al.,23 is that low-energy dynamics involving collective motions of the AlO and PO tetrahedra play a role in solving the local bond frustration over a wide temperature range, without necessitating a global phase transition. The first experimental indication for this form of spatial/temporal disorder was the rich diffuse patterns observed in electron diffraction.23 By analogy to the rigid-unit mode behavior known in SiO2 tridymite,27 it was proposed the underlying disorder is dynamic rather than static, as tetrahedra explore a variety of rotated configurations. A recent inelastic neutron spectroscopy (INS) study found direct evidence for corresponding dynamics over a wide temperature range.15 These initial results lend credibility to the idea that the global phase transition in AlPO4-5 is dynamically avoided; however, so far, the missing step has been to clarify the exact mechanisms responsible for the apparent average structure. The purpose of this paper is to elucidate the dynamic mechanism for the apparent crystal symmetry by identifying the source and temperature dependency of the quasielastic neutron scattering (QENS) in AlPO4-5. The QENS signal provides a fingerprint for the confined diffusion mechanics of the oxygen atoms at the corner of the rotating tetrahedra. This proves that the Al−O−P bond configurations change dynamically over a broad range of conformations on a variety of time scales. Using comparison with ab initio molecular dynamics simulations, we can show that the QENS signal is also linked to the rigid-unit modes in the network which facilitates restricted and semirandom motion of the bridging oxygen species. The relatively large oxygen jumps produce an average position of the oxygen vertices in agreement with the P6cc positions. Close analysis of the experimental QENS also indicates a transition from a purely dynamic regime above 100 K to a hidden phase of quasi-static and spatially disordered tetrahedral configurations at low temperature.



RESULTS AND DISCUSSION The neutron spectroscopy data indicate that several distinct types of low-energy dynamics are present in AlPO4-5. Figure 2 shows a constant scattering-vector (Q) slice for AlPO4-5 at 300 K using the PELICAN time-of-flight (TOF) spectrometer with an incoming wavelength of λ = 4.74 Å. The plot includes the higher energy (E) transfer regions. However, nearly all of the intensity is present below 4 meV. The low energy scale dictates that these dynamic modes will be heavily populated by thermal energy and thus can play a key role in determining the resultant dynamic structure. To separate the relative time scales of the different types of motion, the scattering intensity function (S) was fitted with two components representing the quasielastic and inelastic processes: S = SQENS + SINS

(1)

The quasielastic contribution was described using the general form 18763

DOI: 10.1021/acs.jpcc.7b06415 J. Phys. Chem. C 2017, 121, 18762−18770

Article

The Journal of Physical Chemistry C

has the intense QENS component with a width in the range |E| < 2.0 meV as illustrated in Figure 2b. Unlike the strong RUM mode intensity, which is weakly temperature dependent, the QENS intensity and width are both strongly temperature dependent in the range 100−300 K. Figures 3a−c show the data for three temperatures using λ =

Figure 2. (a) Time-of-flight neutron spectroscopy data show a mixture of inelastic and quasielastic features caused by the rigid-unit mode dynamics and confined diffusion in the system. Given the instrumental resolution and the spectral overlap of the features, it is necessary to fit both simultaneously in order to extract accurately the QENS contribution. (b) The QENS contribution consists of two distinct components corresponding to fast and slow motions with Lorentzian profiles. 2

2

SQENS = e(−Q ⟨u ⟩ /3)[A(Q )δ(ω) + (1 − A(Q ))L(Q , ω)] (2)

where A(Q) is the elastic intensity at the scattering vector Q, and L(Q,ω) is the quasielastic function described as a sum of two Lorentzians centered on zero energy, representing confined diffusion of Al, O, and P with two time scales. The meansquared vibrational displacement is represented by ⟨u2⟩. The low-energy phonon/rigid-unit mode was represented by a Lorentzian centered at ωRUM, by setting SINS = L(Q,ω,ωRUM) where L is a single Lorentzian. The latter RUM spectrum is centered at the inelastic energy of ωRUM ∼ 2.2 meV, accounting for the collective coherent motions of neighboring AlO4 and PO4 tetrahedra. The fit function was convoluted with the instrument resolution, and detailed balance was included at each temperature.28 The individual components are plotted for comparison in Figure 2a,b. Qualitatively, the clear appearance of the elastic line coexisting with the broader quasielastic intensity implies that spatially confined diffusion is present for certain chemical species around well-defined average positions.28 Meanwhile, the inelastic Lorentzian at 2.2 meV matches the feature reported in the previous study, linked to collective rigid-unit modes.15 Although the focus of the current work is the QENS, it is critical to note that a substantial “background” arises near E = 0 meV from the spectral overlap with the broad inelastic RUM signal. For this reason, all fits necessarily include the rigid-unit mode spectrum, and there is some unavoidable uncertainty in extracting the integrated intensities of the broader quasielastic component because it overlaps with the tail of the RUM signal. The two signals, however, do become clearly independent at lower temperatures. For clarity, in most subsequent plots, we only show the quasielastic window that

Figure 3. Quasielastic neutron signal is Q-dependent at (a) 300, (b) 200, and (c) 100 K and can be fitted using the two-Lorentzian model described in the text. The markers are the experimental data; the solid lines are fits to the data.

4.74 Å. Several Q values are illustrated together with the fit using the model above. The quasielastic component becomes stronger at higher Q relative to the elastic line, as expected for a confined diffusion. By comparison with the data at 300 K (Figure 3a), the relative intensity of the quasielastic contribution is lower at 200 K (Figure 3b). At 100 K (Figure 3c), the QENS contribution lies near the limits of detection of the instrument resolution (ΔR = 132 μeV) and signal-to-noise ratio (∼5 × 10−4). Indeed, the time-of-flight measurements show no detectable quasielastic signal at all below 100 K. This indicates that the relevant diffusive motions have either become slower than the instrument resolution or become confined in a real-space volume that is too tiny to measure in the accessible Q range. Normally this would be suggestive of a phase transition; however, complementary neutron powder diffraction patterns (see the Supporting Information) indicate the system remains 18764

DOI: 10.1021/acs.jpcc.7b06415 J. Phys. Chem. C 2017, 121, 18762−18770

Article

The Journal of Physical Chemistry C

Figure 4. (a) Temperature-dependent inelastic spectrum shows a decrease of inelastic intensity upon cooling. (b) The three low-energy components each show a different temperature dependency, and the QENS signal is not detectable below 100 K whereas the RUM mode persists to ∼30 K. (c) At low temperatures, the inelastic component is very small as the system enters a quasi-static state with frozen dynamics on all time scales detected within the instrument resolution (1.2 μeV).

Figure 5. (a) Collective motion of the tetrahedra produces a strong low energy inelastic feature at ∼2.2 meV which is reproduced in the S(ω,Q) calculated from the molecular dynamics. Both the calculated and experimental spectra have been integrated over Q between 0 and 2.5 Å. (b) Classical molecular dynamics calculations successfully also reproduce the main feature, and in comparison, the larger real-space simulations indicate that finite-size effects introduce a small blue-shift. (c) Confined motion of Al, P, and O species causes a weak, but clear, Q dependence in the experimentally determined EISF from the quasielastic signal at 200 and 100 K. The shaded lines indicate the range of experimental uncertainty caused by the overlapping RUM and QENS signal. The DFT dynamics reproduce the same temperature dependency, whereas the classical MD CALTROW molecular dynamics do not capture the high-temperature diffusion processes accurately. The experimental EISF was determined with an instrument resolution of ΔR = 0.132 meV.

in the average P6cc symmetry in the entire temperature range (3−300 K), in agreement with most past reports. The physical picture is therefore that the tetrahedra are dynamically disordered at high temperature and transition to a statically disordered state at low temperature. To rule out any residual diffusive dynamics below 100 K, we carried out further measurements on PELICAN with a neutron wavelength λ = 6 Å providing an energy resolution (ΔR = 65 μeV) and also on EMUthe backscattering spectrometer with an energy resolution of ΔR = 1.2 μeV which is ideal for observing ultraslow dynamics. Figure 4 summarizes the experimental results of the inelastic and quasielastic experiments. Figure 4a shows the temperature dependence of the experimental signal plotted on a linear scale. Importantly, the RUM inelastic feature systematically vanishes on the neutron energy gain side as the temperature is cooled below 50 K. This is expected based on the corresponding temperature scale with the Bose−Einstein factor set as T* = E/ kB, giving the characteristic temperature as T* = 30 K. Fits according to the model described above yielded the characteristic frequencies for the main components described in Figure 4b. Although the RUM peak loses intensity, it has a similar frequency at all temperatures, showing only a small blue-shift. This indicates that the rigid-unit modes have a strict characteristic frequency, governed by the bond strength, topology, and angular potential well for the Al−O−P bond.

Nevertheless, the associated energy cost dictates that the system becomes increasingly static at low temperatures. Meanwhile, below 100 K, the quasielastic intensities also drop toward the background, indicating a smaller real-space volume as the diffusive motion amplitude decreases. In parallel, the widths of the quasielastic components decrease at lower temperature reflecting a slowing down of the diffusional processes, approaching the TOF instrument resolution by 100 K. To clarify whether any quasielastic intensity remains below 100 K, we thus turn to the backscattering spectroscopy data which has higher energy resolution. Figure 4c shows the elastic intensity integrated over all detectors and normalized to unity at 1.5 K. The elastic intensity is maximized below 100 K and drops to a value of approximately 0.50 at high temperatures. The complementary time-of-flight and backscattering data therefore indicate that there is no quasielastic intensity below 30 K, suggesting the diffusive motions have effectively become frozen across the broad time scale spanning 300 MHz to 5 THz. To understand better how these temperature-dependent dynamics are linked to average crystal symmetry, we carried out molecular dynamics simulations using DFT and complementary classical dynamics with the well-established CALTROW force field. To extend on past work,15 the simulations were run above 100 K, and long runs were performed over 20 ps for DFT dynamics and 100 ps for the classical dynamics. Figure 5a 18765

DOI: 10.1021/acs.jpcc.7b06415 J. Phys. Chem. C 2017, 121, 18762−18770

Article

The Journal of Physical Chemistry C

Figure 6. (a) Experimentally determined average structure has Al−O−P bond angles close to 180°. (b) Stop-motion image from the DFT molecular dynamics at 300 K from 0 to 5 ps shows that the oxygen atoms explore a relatively large volume corresponding to a range of Al−O−P bond angles centered around a mean value of nearly 150°. (c) Stop-motion image is taken from the DFT molecular dynamics (0−5 ps) at 100 K and shows tightly constrained diffusion of the oxygen atoms and an Al−O−P angle with a nearly constant value. (d) Time and spatial average of the oxygen position at various temperatures visualized as an isosurface around the average bond axis. (e) Histogram showing the range of the instantaneous interplanar bond angle in the simulations at several temperatures. (f) Histograms showing the range of interplanar bond lengths at various temperatures in the simulations.

shows the S(ω,Q) calculated from the DFT molecular dynamics at 100 K, integrated over the same Q range as the experimental result. Like the experiment, a spectrum identified with the rigid-unit modes15 is observed with the strongest peak in the range 2−3 meV. The DFT molecular dynamics are clearly able to reproduce the main low-energy features. Furthermore, by decomposing the calculated S(ω,Q) into its incoherent and coherent contributions, we find that the RUM spectrum has components from both types of scattering in a ratio of approximately 1:3. The larger coherent contribution matches well with the main features of the low-energy rigid-unit mode spectrum described in previous work,15 implying collective motion of atoms in the tetrahedra. The nonzero incoherent contribution indicates that the frequencies governing collective motions are also related to self-correlations in the vertex coordinates. This simply means that, as expected, the rigid-unit mode dynamics from tilting the linked tetrahedra automatically generates same-frequency self-correlations in oxygen at the corners of the tetrahedra via librational− translational coupling. To understand the small energy shift in the calculated spectrum, we also conducted classical MD simulations which have the merit of accessing larger spatial scales (up to 512 atoms). Figure 5b shows the simulated spectrum based on the classical MD dynamics for a single unit cell and the 2 × 2 × 2 supercell. The larger unit cell gives excellent quantitative agreement with the inelastic data, whereas the single-unit-cell simulations show the presence of a finitesize effect that blue-shifts the calculated spectrum, as expected for smaller correlation lengths. Having validated the general features of the modeling, we then calculated the quasielastic intensity based on the molecular dynamics. Figure 5c shows the experimentally measured and calculated EISF for two temperatures. To make an accurate comparison between the integrated intensities in the calculation and experiment, the experimental EISF was taken as equal to the

fraction of the elastic component (A(Q)) to the overall QENS intensity from the fits using eq 2. In contrast to the initial study,15 the RUM spectrum in this work was carefully excluded from the EISF since it is not strictly part of the quasielastic scattering (it has nonzero energy transfer). The range of experimental values is a shaded plot reflecting the choice of constraining the RUM amplitude in the fitting model to minimum values or maximum values which result in statistically good fits (figure-of-merit χ2 between 2.0 and 2.5). The theoretical EISF calculated with the DFT-MD falls well within the experimental range for 200 and 100 K. In common with experiment, the ab initio dynamics predict a weak but observable quasielastic intensity (corresponding to EISF < 1). By decomposing the theoretical DFT EISF into elemental contributions, we found that oxygen makes the dominant contribution to the quasielastic intensities. Both the theoretical and experimental EISF show a weak Q dependence over the measured scattering vector range (0.4 < Q < 2.5 Å−1), indicating that the confined diffusion for the oxygen occurs within small spatial volumes around the average lattice points with an effective radius of ∼0.8 Å. As in the experiment, the quasielastic signal effectively vanishes for low temperatures (below 100 K) in the DFT simulations (corresponding to EISF = 1) over the accessible Q range and time range. This shows that the simulations correctly indicate that some species have insufficient energy to cross an energy barrier in the system. While the temperature dependence of the signal is reproduced well in the ab initio microcanonical simulations, surprisingly the classical force fields fail to reproduce the high-temperature dynamics showing an EISF that remains close to 1 at all temperatures. This indicates that the activation barriers for the diffusive motions are overestimated in the classical force field in contrast to those calculated with the ab initio method. The solid agreement with experiment validates the DFT-MD results and motivates a closer analysis of the individual atomic 18766

DOI: 10.1021/acs.jpcc.7b06415 J. Phys. Chem. C 2017, 121, 18762−18770

Article

The Journal of Physical Chemistry C

To understand better the time scales associated with the sporadic jump motion, and to move beyond the average picture, we performed further analysis of certain individual trajectories focusing on “problematic” O2 atoms at the bridge between AlO and PO planes of tetrahedra. Figure 7 shows a representative trajectory for one such oxygen at 300 K decomposed into Cartesian coordinates where Z is taken to point along the c-direction. The oxygen mostly undergoes displacements primarily in the XY plane whereas the Z

trajectories to identify the specific motions and barriers to explain the emergence of the average structure. Figure 6a shows the average structure of AlPO4-5 reported by Klap13 which contains the interplanar bond angle close to 180°. For comparison, Figures 6b,c show the dynamic structures at 300 and 100 K visualized as stop-motion images of the chemical species for 5 ps from the DFT calculations, in 10 fs frame steps. In contrast to the average structure determined in diffraction, the dynamic picture shows the instantaneous Al−O−P angle is almost always much smaller than 180°; however, the precise orientation varies as a function of time. In particular, there is considerable freedom to rotate the interplanar Al−O−P bond around the c-axis, so that in practice the time-averaged central position does indeed lie along the 180° position above 100 K, as the oxygen move within considerable volumes. By taking the spatial and time average of all interplanar Al−O−P bonds in the system, Figure 6d shows the probability isosurfaces of finding the interplanar oxygen within a spatial region occupied at least 10% of the time at 300, 100, and 1.5 K. At higher temperatures, this corresponds to a pancake volume with a radius in the ab plane of about 0.8 Å at 300 K and c-axis height 0.3 Å centered on the 180° bond axis. However, at lower temperatures in the regime of 100 K, the motion of the bridging oxygen occupies a doughnut-shaped isosurface with a slightly reduced radius (0.6 Å). The interplanar oxygen undergoes dramatic precessions around the energetically unfavorable bond angle, thereby allowing for a longer instantaneous bond length, and this is mediated by some smaller motions for the in-plane oxygen positions. Yet the average position remains close to the 180° bond at all temperatures. The simulations at very low temperature (1.5 K) indicate that eventually the motion of the oxygen becomes suppressed, and instead, these atoms occupy sites with increasing levels of discreteness near six energy minima positioned on a ring around the 180° bond position. To quantify these traits, Figure 6e shows a histogram of all the bond angles for all bridging bonds in the 20 ps trajectories. It is worth noting that most bond angles lie close to a mean of 150° at all temperatures in agreement with the general expectation for Al−O−P materials. However, the width of the bond angles distribution depends on temperature. Figure 6f shows the matching histogram for the distribution of bond lengths in the molecular dynamics simulation. The simulations indicate the average instantaneous bond length of the P−O interplanar bond (1.55 Å) is considerably larger than the bond length from the time-averaged structure (1.45 Å), whereas the difference is far less noticeable for the Al−O bond. This result essentially confirms that rigid-unit modes are indeed the dominant dynamic mechanism in AlPO4-5 because the Al−O and P−O lengths are found to vary less than 0.2 and 0.1 Å, respectively, in the molecular dynamics. The latter bond flexing is far too small to account for the larger range of spatial positions explored by the dynamic oxygen. However, the combination of bond flexing and small translations allows a certain degree of randomness in the tetrahedral network, so that not all rotations necessarily involve a perfectly matching counter rotation in the neighboring tetrahedra. Consequently, these dynamics should be considered “quasi” rigid-unit modes rather than perfect rigid modes. One interesting corollary is that some of the Al−O−P bond angles at 300 K sporadically have values close to 180°, indicating a rare type of jumping motion of the oxygen which is a key factor that enlarges the regions of confined diffusion.

Figure 7. (a) Relative positions of a single interplanar oxygen from its time average, extracted from the DFT molecular dynamic trajectory at 300 K. The (XY) plane components (green and red) show much larger deviations than the Z component (blue), indicating the motion is approximately confined within a circle. (b) Radial distance in the XY plane of the same oxygen from its time average position. Sporadic jumps across the circle correspond to low values of the instantaneous radius as indicated by the ∗ symbols. (c) Top-down view of the entire trajectory at 100 K with the line shading reflecting the chronological sequence. At low temperature, the oxygen is confined within one region of a circle. (d) At 300 K, the trajectory indicates multiple crossing events as the oxygen jumps across the circle. (e) Schematic representation of the circular motion with conserved radius. (f) Schematic illustration of the jump motion where the radius varies. (g) Fourier transform of the correlation functions indicating different time scales for the two types of motions. 18767

DOI: 10.1021/acs.jpcc.7b06415 J. Phys. Chem. C 2017, 121, 18762−18770

Article

The Journal of Physical Chemistry C

analogous to a glass transition. It is somewhat surprising that the average structure appears to be very similar for the timedisordered and statically disordered states. On the basis of the lack of a long-range structural transition, we propose that the material can be cooled into a variety of metastable configurations that are explored by its high-temperature dynamics. The resulting configurational randomness provides a natural explanation for the residual entropy observed in heat capacity measurements of AlPO4-5.15 The flexibility of this compound provides a compelling example of how thermal excitations can be key to stabilizing an otherwise unlikely crystal structure. More broadly, it is reasonable to expect similar effects in many other materials, e.g., metal−organic framework (MOF) and other members of the zeolite family. Unveiling the dynamic mechanisms for bond-frustrated structures is a key first step in the grand challenge of improving computer-based prediction of functional materials which defy the 0 K approximation.

coordinate varies very little (Figure 7a). Furthermore, by taking the average across the whole trajectory, it is clear the individual O have average positions at the center of a circle, corresponding to the 180° Al−O−P bond. However, in general, the O atom spends almost no time at the central position. Instead, the majority of the time is spent at the circumference except for a few sporadic jumps when the O atom approaches the central position for a brief instant as it attempts to jump across the ring. To quantify the distance of the oxygen from the average position, we define a radius r(t) as the distance between the instantaneous position and the average position. Figure 7b shows that the radius has a well-defined nonzero average value except at a few times denotated by the symbols (∗). The circular motions which preserve the radius have a higher frequency (>1 THz) compared to the sporadic jump motions which modify the radius (