Superprotonic phase transition of $\mathrm{Cs}\mathrm{H}\mathrm{S}{\mathrm{O}}_4$: A molecular dynamics simulation study

The superprotonic phase transition (phase II → phase I; 414 K) of cesium hydrogen sulfate, $\mathrm{Cs}\mathrm{H}\mathrm{S}{\mathrm{O}}_4$, was simulated using molecular dynamics with the “first principles” MSXX force field (FF). The structure, binding energy, and vibrational frequencies of the $\mathrm{Cs}\mathrm{H}\mathrm{S}{\mathrm{O}}_4$ monomer, the binding energy of the ${\left({\mathrm{H}}_2\mathrm{S}{\mathrm{O}}_4\right)}_2$ dimer, and the torsion barrier of the $\mathrm{H}\mathrm{S}{\mathrm{O}}_4^{-}$ ion were determined from quantum mechanical calculations, and the parameters of the Dreiding FF for Cs, S, O, and H adjusted to reproduce these quantities. Each hydrogen atom was treated as bonded exclusively to a single oxygen atom (proton donor), but allowed to form hydrogen bonds to various second nearest oxygen atoms (proton acceptors). Fixed temperature-pressure (NPT) dynamics were employed to study the structure as a function of temperature from 298 to 723 K. In addition, the influence of several force field parameters, including the hydrogen torsional barrier height, hydrogen bond strength, and oxygen charge distribution, on the structural behavior of $\mathrm{Cs}\mathrm{H}\mathrm{S}{\mathrm{O}}_4$ was probed. Although the FF does not allow proton migration (i.e., proton jumps) between oxygen atoms, a clear phase transition occurred as demonstrated by a discrete change of unit cell symmetry (monoclinic to tetragonal), cell volume, and molar enthalpy. The dynamics of the $\mathrm{H}\mathrm{S}{\mathrm{O}}_4^{-}$ group reorientational motion also changed dramatically at the transition. The observation of a transition to the expected tetragonal phase using a FF in which protons cannot migrate indicates that proton diffusion does not drive the transition to the superprotonic phase. Rather, high conductivity is a consequence of the rapid reorientations that occur in the high temperature phase. Furthermore, because no input from the superprotonic phase was employed in these simulations, it may be possible to employ MD to hypothesize superprotonic materials.

Figure 1

(Color online) Crystal structure of $\mathrm{Cs}\mathrm{H}\mathrm{S}{\mathrm{O}}_4$: (a) monoclinic phase II (Ref. 10) and (b) tetragonal phase I as proposed by Jirak et al. (Ref. 11). In (b), each oxygen position has half occupancy and the hydrogen atoms are placed in the middle of the disordered hydrogen bonds (dashed lines).

Figure 2

(Color online) Possible configurations of the sulfate tetrahedra in the superprotonic phase: the structures by Jirak et al. (a), (Ref. 11), Merinov et al. (b), (Ref. 13), and Belushkin et al. (c), (Ref. 14) have two, two, and four orientations, respectively, which transform into each other by rotations of 32°, 30°, and 30°, respectively. The * designates one possible arrangement for the oxygen atoms of a tetrahedra.

Figure 3

(Color online) Structures used to adjust the FF parameters: (a) $\mathrm{Cs}\mathrm{H}\mathrm{S}{\mathrm{O}}_4$ monomer, (b) ${\left({\mathrm{H}}_2\mathrm{S}{\mathrm{O}}_4\right)}_2$ dimer, and (c) $\mathrm{H}\mathrm{S}{\mathrm{O}}_4^{-}$ ion [projected down the S-O(H) bond]. The $\mathrm{Cs}\mathrm{H}\mathrm{S}{\mathrm{O}}_4$ monomer (a) was used to determine the vdW parameters for Cs and all $\mathrm{H}\mathrm{S}{\mathrm{O}}_4^{-}$ FF parameters except for the hydrogen bond and the O-S-O-H torsional terms, which were adjusted with the results for (b) and (c), respectively.

Figure 4

Potential energy curves for a $\mathrm{H}\mathrm{S}{\mathrm{O}}_4^{-}$ ion with fixed $\mathrm{O}\left(1\right)\text{−}\mathrm{S}\text{−}{\mathrm{O}}_{\mathrm{D}}\text{−}\mathrm{H}$ torsional angles: (a) optimized by QM and (b) minimized with the adjusted FF. The calculations were carried out from 0° to 60° and extended over a full 360° using symmetry.

Figure 5

(Color online) Calculated IR spectra for MD simulations: For phases II and I of $\mathrm{Cs}\mathrm{H}\mathrm{S}{\mathrm{O}}_4$, (a) and (b), respectively. Peaks with strongest absorption are marked with an asterisk and the frequencies at the peak maximums are compared to those of the strongest peaks found by two separate infrared spectroscopy experiments on polycrystalline $\mathrm{Cs}\mathrm{H}\mathrm{S}{\mathrm{O}}_4$ (Refs. 20, 36). In both phase II (c) and phase I (d) the simulation and experimental peak positions are very similar (for perfect match ⇒ slope$=1$). A plot of the experimental values below versus those above the superprotonic transition (not shown) results in a slope of 0.98(1), $R^2=0.999$, as the measured IR spectra changes little over the phase II-I transition.

Figure 6

(Color online) Cell parameters as a function of temperature (MD simulations): average values calculated from the final 150 ps of each 300 ps simulation. Lines indicate the two transitions at 598 and 698 K.

Figure 7

Calculated x-ray powder diffraction patterns (time averaged from 150 ps of MD, 1500 frames). (a) Phase II $\mathrm{Cs}\mathrm{H}\mathrm{S}{\mathrm{O}}_4$ calculated from the experimental structure (Ref. 10). (b) and (c) Phase II calculated from MD simulations at 298 and 573 K, respectively. (d) Phase I $\mathrm{Cs}\mathrm{H}\mathrm{S}{\mathrm{O}}_4$ calculated from the Jirak et al. experimental structure (Ref. 11). (e) Phase I calculated from MD simulations at 623 K. (d) X-ray diffraction pattern calculated from MD simulations at 723 K, above the second (melting) transition.

Figure 8

(Color online) (a) Potential energy and (b) volume as a function of temperature from MD simulations: 298 k to 723 k.

Figure 9

Laboratory (fixed) frame and equations used to calculate the orientation of a tetrahedron’s $\mathrm{S}\text{−}{\mathrm{O}}_{\mathrm{D}}$ vector.

Figure 10

Probability distribution functions for the S-O vectors: (a), (b), (c), (d), and (f) were created using only the orientations of the S85-O86 (donor oxygen) at 298, 498, 573, 598, and 723 K, respectively, whereas (e) is a combination of the orientations for all four S-O vectors of the S85 tetrahedron in the temperature range 598–673 K (phase I). The directions of the S-O vectors were converted into the polar coordinates $(\theta ,\varphi )$ and then mapped onto the 2 D $(\theta ,\varphi )$ space. A probability distribution function was created by partitioning the full ranges of each variable $(0⩽\theta <\pi ,0⩽\varphi <2\pi )$ into a $20\times 20$ matrix and assigning each $(\theta ,\varphi )$ pair to a particular cell.

Figure 11

(Color online) Probability distribution functions for all S-O vectors in phase I (from the 598–673 K simulations). The four distinct directions correspond fairly well to S-O directions derived from the published structures (Refs. 11, 13, 14).

Figure 12

Orientation/reorientation of $\mathrm{H}\mathrm{S}{\mathrm{O}}_4^{-}$ ion defined by its $\mathrm{S}\text{−}{\mathrm{O}}_{\mathrm{D}}$ vector: (a) an $\mathrm{H}\mathrm{S}{\mathrm{O}}_4^{-}$ ion has the vector pointing from its S to ${\mathrm{O}}_{\mathrm{D}}$ atom initially aligned vertically; (b) after a time $t^{\prime }$, the $\mathrm{H}\mathrm{S}{\mathrm{O}}_4^{-}$ ion has changed its orientation and the S to ${\mathrm{O}}_{\mathrm{D}}$ vector now forms an angle $\theta \left(t^{\prime }\right)$ with the vertical.

Figure 13

Autocorrelation functions for all 32 tetrahedra evaluated using their average cosine values between the reference orientations (at 150 ps) and subsequent orientations: (a) average cosine value versus temperature. Note the decrease in the average cosine value across the phases II to I transition and that the average cosine values continue to decrease above the transition; (b) time dependence of the average cosine value at 573 K which equilibrates around 0.73, (c) same for 598, (d) same for 623 K, which drops to zero in 150 ps.

Figure 14

Average angular velocity of all 32 $\mathrm{S}-{\mathrm{O}}_D$ vectors versus temperature. Note that the jump in this value at the phase II-I transition is only $\sim 25\%$. Therefore the dramatic increase in protonic conductivity across the transition (100 000–1 000 000 %) is due to a change in the nature of the reorientations (e.g., oscillation to libration/rotation) rather than the degree of the reorientations.

Figure 15

Pictorial representation of the difference between a three-fold, axial rotation around the $\mathrm{S}\text{−}{\mathrm{O}}_D$ vector and an extensive reorientation involving the $\mathrm{S}\text{−}{\mathrm{O}}_D$ vector.

Figure 16

(Color online) Rearrangements of the sulfate tetrahedra below the superprotonic phase transition at 598 K. Note the significant number of axial rotations (even near ambient temperatures), while the number of extensive reorientations remains at or very near zero.

Figure 17

(Color online) Sulfate group reorientations involving the donor oxygen for simulations exhibiting the superprotonic phase of $\mathrm{Cs}\mathrm{H}\mathrm{S}{\mathrm{O}}_4$ (598–673 K). Numbers were calculated using two sets of FWHM values for the regions of high orientational density in Fig. 11: 0.47 rad in $\theta$ and 0.68 rad in $\varphi$ (the average values) and 0.52 rad in $\theta$ and 0.80 rad in $\varphi$ (the averages increased by their respective standard deviations).

Figure 18

(Color online) Unit cell parameters as a function of temperature, as obtained from MD simulations using equal charges on all four oxygens. The superprotonic transition initiates at 423 K and is complete (i.e., $\beta =120°$ and $\mid \mathrm{a}\mid =\mid \mathrm{b}\mid$) at 523 K.

Figure 19

(Color online) Potential energy as a function of temperature, as obtained from MD simulations using equal charges on all four oxygens. Discontinuities in the potential energy function are apparent at 523 K, at the completion of the superprotonic transition, and at 673 K, at the melting transition.

Figure 20

(Color online) Calculated x-ray diffraction patterns: (b), (c), and (e) from MD simulations using equal oxygen charges at the temperatures indicated, (a) from the structure determined by single crystal x-ray diffraction at 298 K (Ref. 11) and (d) from structure determined by powder x-ray diffraction at 415 K (Ref. 10). Comparison of the patterns in (a), (b), and (c) show that the simulations at 398 and 423 K are in phase II $\mathrm{Cs}\mathrm{H}\mathrm{S}{\mathrm{O}}_4$, whereas the diffraction pattern at 523 K, (e) matches the phase I pattern calculated from experiment (d).

Figure 21

(Color online) Rearrangements of the sulfate tetrahedra across the low temperature transitions of the simulations with increased hydrogen bond strength: (a) 1, (b) 1.5, and (c) 2 times the original hydrogen bond energies. Note the jump in axial rotations at precisely the same temperatures the lattice constants were observed to change, (Ref. 51), while the number of (extensive) reorientations remains at or very near zero.

Figure 22

Autocorrelation functions for simulations with (a) decreased and (b) increased torsional barrier heights. The values are lower in (a) and higher in (b) when compared to the original simulations Fig. 13a.