Robust Diffusive Proton Motions in Phase IV of Solid Hydrogen - The

May 14, 2014 - The present study highlights an often overlook issue in first-principles calculations that long time MD is essential to achieve ergodic...
2 downloads 10 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCC

Robust Diffusive Proton Motions in Phase IV of Solid Hydrogen Hanyu Liu,†,‡ John Tse,†,‡ and Yanming Ma*,† †

State Key Laboratory of Superhard Materials, Jilin University, Changchun 130012, China Department of Physics and Engineering Physics, University of Saskatchewan, Saskatoon, Canada S7N 5E2

‡

S Supporting Information *

ABSTRACT: Systematic first-principles molecular dynamics (MD) simulations with long simulation times (7−13 ps) for phase IV of solid hydrogen using different supercell sizes of 96, 288, 576, and 768 atoms established that the diffusive proton motions process in the graphenelike layer is an intrinsic property and independent of the simulation cell sizes. The present study highlights an often overlook issue in first-principles calculations that long time MD is essential to achieve ergodicity, which is mandatory for a proper description of dynamics of a system. It is inappropriate to make arguments on the analysis of MD results, which are far from ergodic. Furthermore, we have simulated the vibrational density of states of phase IV based on our proton diffuse model at a pressure range of 225−300 GPa, which is qualitatively in agreement with experimental data.

■

INTRODUCTION Solid hydrogen is one of the central topics of condensed matter physics.1,2 Recent advances in experimental technique have led to the discovery of a new phase (phase IV) at room temperature.3−8 However, the structural and dynamical properties of this phase IV are under extensive debate.4,9,10 Initially, phase IV was interpreted as having an orthorhombic Pbcn symmetry (24 molecules per cell),4 which consists of alternate layers of graphene-like, three-molecule rings with elongated H2 molecules and unbound H2 molecules, which was predicted in 2007 by Pickard and Needs.11 Theoretical phonon calculations, however, show this structure is dynamically unstable. A monoclinic Pc structure that doubles that of the Pbcn structure with 96 atoms per cell and does not contain any imaginary frequency was proposed later.12 Meanwhile, our first-principles metadynamics study13 proposed an alternate Cc structure containing 96 molecules per unit cell, geometrically rather similar to Pc structure. Our calculations suggested13 for the first time that phase IV has a partially disordered structure: disorder in the weakly coupled or unbounded H2 molecules layers but ordered in the strongly coupled graphene-like layers (G-layer). More recently, a novel diffusive hydrogen process (i.e., hydrogen hoppings between two nearest molecules)7,14 with collaborative rotation of three hydrogen molecules in a ring fashion14 in the graphene-like layers at a pressure range of 250−350 GPa and temperature range 300−500 K was reported. The finding of the proton diffusion provides an explanation on the observed abrupt increase of Raman line width at the formation of phase IV and its large increase with pressure.14 Other ab initio MD simulations15,16 on phase IV found similar results.7,14 These studies show (i) the trend of the experimental Raman spectrum can be explained by the Pc or Cc structures;16 (ii) at 250 GPa and 220 K, a simulation started from phase III © XXXX American Chemical Society

(C2/c) was found to transform into phase IV by heating the system to 300 K, a direct confirmation of our previous results;14,15 and (iii) rotation of the 3-molecule ring in the Glayer of phase IV predicted in refs 13 and 14 was confirmed by another MD simulations.15 In ref 15, the authors however claimed that there are more than one G-layer in phase IV (G′ and G″) and diffusive hydrogen motions do not exist. They further argued that these discrepancies with other MD calculations7,14 were attributed to finite size effects.15 It is noteworthy that a very recent ab initio MD calculation17 also reported atomic diffusion in a molecular solid hydrogen phase (Cmca-12) when temperature goes beyond 429 K in the pressure range 260−350 GPa. Moreover, the hydrogen transfer was found to confine within one type of layer,17 similar to the diffusive hydrogen motions observed in the G-layer of phase IV.7,14 A major difference between ref 15 and the other studies7,14 is the length of the simulation time. In refs 7 and 14, MD simulation times are ranging from 7 ps14 to 90 ps.7 This is in clear contrast to the simulation time span of only 1−3 ps in refs15 and 16. It is well established that short-time MD simulation does not sample the potential landscape adequately and might not capture important events. In this study, firstprinciples MD simulations with longer simulation times (7−13 ps) at 300 GPa and 300 K for phase IV of solid hydrogen using different supercell sizes of 96, 288, 576, and 768 atoms are systematically performed to examine the finite size effects. Through carefully analysis of the atomic trajectories, it is found that the above discrepancies on the description of the hydrogen Received: April 7, 2014 Revised: May 9, 2014

A

dx.doi.org/10.1021/jp503409p | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

dynamics in the G-layers are not due to cell sizes but rather on the time span of the MD calculations. Rationally, results derived from short-time MD simulations that are not long enough to satisfy the ergodicity assumption can be unreliable. Our current results confirmed previous MD results7,14 and allow us to conclude that the diffusive proton motions in G-layers for phase IV are robust.

■

COMPUTATIONAL DETAILS Ab initio MD simulations for solid hydrogen were performed at 300 GPa and 300 K with the NPT18,19 (N = number of particles, P = pressure, T = temperature) ensemble recently implemented in Vienna Ab initio Simulation Package (VASP) code.20−22 The all-electron projector-augmented wave (PAW) method was adopted.23 Exchange and correlation effects were treated in generalized-gradient approximation (GGA).24 We have adopted simulation models containing 96, 288, 576, and 768 atoms for phase IV with simulation times ranging from 7 to 13 ps. A plane-wave cutoff of 500 eV was used. Brillouin zone sampling of 4 × 4 × 4, 2 × 2 × 4, 2 × 2 × 2, and 2 × 2 × 2 kmeshes25 for supercells consisting of 96 (Cc), 288 (Pc), 576 (Cc), and 768 (Cc) atoms, respectively, was employed. The Pc structure with 288 atoms was created by using a supercell of 2 × 2 × 1 Pc structure (96 atoms). The Cc structure with 768 atoms was created by using a supercell of 2 × 2 × 2 Cc structure (96 atoms). Note that the Cc structure with 576 atoms from the metadynamics simulation result with initial phase of phase I for solid hydrogen at 250 GPa and 300 K (disordered hexagonal structure). The chosen computational parameters were sufficient to produce converged results. The time step used in the MD simulation was 0.5 fs, and the self-consistency on the total energy was chosen to be 1 × 10−5 eV.

Figure 1. (a) Mean-square displacements as a function of time (7−13 ps) at 300 GPa and room temperature for phase IV of solid hydrogen with different simulated supercells of 96, 288, 576, and 768 atoms are shown. Layer I and II are a freely rotating hydrogen molecules layer and strongly coupled graphene-like layer (G-layer), respectively. (b) and (c) show the crystal structure of Cc phase.

and (ii) reconstruction of all molecules with time in layer II. This leads to MSD of hydrogen atoms in layer II looks like step behavior. It is noteworthy that the reconstruction process of all molecules in layer II occurs at each 2−3 ps. From the analysis of the averaged atomic positions of the Pc structure,15 Magdau and Ackland proposed that phase IV of solid hydrogen adopts a dynamic, simple hexagonal structure with three layers: a freely rotating hydrogen molecules layer, a static hexagonal trimers layer, and a rotating hexagonal trimers layer. It is correct to interpret the crystal structure with the averaged atomic positions if the atoms are not diffusive. However, it is not applicable in phase IV. As already shown above, the hydrogen atoms in G-layer diffuse substantially. Therefore, characterization of this “liquid-like” structure by averaging the atomic positions may not be appropriate. To demonstrate such a failure, we have computed the averaged atomic positions of phase IV for different running time durations at 300 GPa and 300 K (Figure 2). The averaged atomic positions derived from a short-MD run (1.5 ps) generated indeed structures with two rotating 3-molecule Glayers (Figures 2a,c) and two static hexagonal 3-molecule Glayers (Figures 2b,d). On the contrary, the averaged atomic positions from a longer MD simulation (5 ps) show the presence of rotating 3-molecule rings in all G-layers (Figure 2e−h). Eventually, the averaged atomic positions from the trajectory of longest MD run (13 ps) gave featureless structure due to the extensive diffusive hydrogen motions in G-layers (Figure 2i−l) as expected. Figure 3a shows the trajectory of the G-layer derived from the last 7.5 ps run of a MD simulation using 768 atoms at 300 K and 300 GPa. It failed to show the 3-molecule ring character. In contrast, the short time trajectory did reveal the molecular nature, as shown in Figure 3b−d, where the G-layers reconstructed with time. The analysis presented here clearly shows the occurrence of diffusive hydrogen motions. It is

■

RESULTS AND DISCUSSION The atomic mean-square displacement (MSD) is an important method to investigate the diffusion of a system. This can be done by calculating MSD as a function of time from MSD = (1/N)∑Ni=0⟨[ri(t) − ri(0)]2⟩, where N is the total number of total atoms in the system, ri(t) is the atomic position at a time of t, ri(0) is the initial position of atom i, and ⟨...⟩ represents the ensemble average. We computed the MSD of hydrogen for phase IV as derived from the present MD simulations at 300 GPa and 300 K for different supercell sizes. The results summarized in Figure 1 show unambiguously that the hydrogen MSD in the G-layer (or layer II) grows with time, even in the largest supercell (768 atoms). This observation demonstrates unequivocally the proton diffusion at 300 GPa and 300 K. There is a weak dependence on the MSD values with the size of the simulation cells, showing that the MSD decreases with increasing cell size at smaller size regime. However, this value is well converged when the cell size reaches 288 atoms or beyond (Figure 1). From the slopes of the MSD curves, the averaged diffusion coefficients are estimated to be ∼0.7 × 10−5 cm2/s for supercells of 288, 576, and 768 atoms. This diffusion magnitude is comparable to that for a typical liquid, e.g., ∼2.3 × 10−5 cm2/ s in liquid water at 298.16 K.26 In this manner, hydrogen in the G-layers may be considered as a sublattice melting because the value of MSD is similar to that of typical liquid water. However, unlike water, hydrogen diffusion in phase IV of the solid is not a random process but follows a well-defined pattern through a series of hopping motions. These hopping motions consist of two important processes: (i) three-molecule cluster rotations B

dx.doi.org/10.1021/jp503409p | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

Figure 2. Averaged atomic positions of different G-layers for phase IV using supercells of 576 atoms at 300 GPa and 300 K from 1.5 ps (top panel), 5 ps (middle panel), and 13 ps (bottom panel) NPT-MD simulation.

inappropriate to characterize the structure of phase IV from the averaged atomic positions via short-time MD running slices.15 To further investigate the proton motion in G-layers, a randomly chosen hydrogen atom, marked by red color, was traced for 10 ps. As shown in Figure 3a, it is evident that this marked hydrogen atom has migrated from one site to another site in the G-layer (or layer II). A supporting movie of the atomic trajectories of the MD simulation has been prepared and deposited as Supporting Information.27 The present study shows unambiguously that diffusive hydrogen motion is ubiquitous and an intrinsic property of phase IV in solid hydrogen. The remarkable feature of phase IV is the emergence of two different high vibrational frequencies above ∼220 GPa.3,4,28 Moreover, one vibrational frequency (ν1 mode) above 4000 cm−1 remains invariant with increasing pressure, while the other one (ν2 mode) is drastically softening. To investigate Raman spectroscopic features of phase IV, we have extracted the vibrational density of states from MD simulations at a pressure

Figure 3. Collected trajectories of long-run MD simulations from the last 7.5 ps (a), 2.5−4.5 ps (b), 6−7 ps (c), and 7.5−8.5 ps (d). Parts b−d clearly show the G-layer reconstructed with time.

Figure 4. (a) Simulated vibrational density of states extracted from MD simulations at a pressure range of 225−300 GPa. (b) It shows the change of two major high vibrational frequencies with creasing pressure, and experimental data are from refs 3 and 4. C

dx.doi.org/10.1021/jp503409p | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

range of 225−300 GPa. It is clearly seen that the ν1 mode frequency softens very rapidly and the ν2 mode frequency remains invariant with pressure (Figure 4). Meanwhile, the full width at half-maximum (fwhm) of ν1 mode was observed by experiment3,4,27 to increase drastically, while the fwhm of ν2 mode remains unchanged. Our current simulations can reproduce this experimental observation and qualitatively agree with previous theoretical results.7,14 The change of the fwhm of ν1 mode strongly correlates the fact that molecules in G-layers are short-lived due to the diffusive proton motions in G-layers of phase IV at high pressures. In contrast, molecules in layer I are ultrastable, which can explain the little-affected Raman fwhm of ν2 mode under high pressure.

(7) Goncharov, A. F.; Tse, J. S.; Wang, H.; Yang, J.; Struzhkin, V. V.; Howie, R. T.; Gregoryanz, E. Bonding, Structures, and Band Gap Closure of Hydrogen at High Pressures. Phys. Rev. B 2013, 87, 024101. (8) Loubeyre, P.; Occelli, F.; Dumas, P. Hydrogen Phase IV Revisited via Synchrotron Infrared Measurements in H2 and D2 up to 290 GPa at 296 K. Phys. Rev. B 2013, 87, 134101. (9) Nellis, W. J.; Ruoff, A. L.; Silvera, I. F. Has Metallic Hydrogen Been Made in a Diamond Anvil Cell? arXiv 2012, 1201, 0407. (10) Eremets, M.; Troyan, I.; Lerch, P.; Drozdov, A. Infrared Study of Hydrogen up to 310 GPa at Room Temperature. High Pressure Res. 2013, 33, 377−380. (11) Pickard, C. J.; Needs, R. J. Structure of Phase III of Solid Hydrogen. Nat. Phys. 2007, 3, 473−476. (12) Pickard, C. J.; Martinez-Canales, M.; Needs, R. J. Density Functional Theory Study of Phase IV of Solid Hydrogen. Phys. Rev. B 2012, 85, 214114. (13) Liu, H.; Zhu, L.; Cui, W.; Ma, Y. Room-Temperature Structures of Solid Hydrogen at High Pressures. J. Chem. Phys. 2012, 137, 074501. (14) Liu, H.; Ma, Y. Proton or Deuteron Transfer in Phase IV of Solid Hydrogen and Deuterium. Phys. Rev. Lett. 2013, 110, 025903. (15) Magdau, I. B.; Ackland, G. J. Identification of High-Pressure Phases III and IV in Hydrogen: Simulating Raman Spectra Using Molecular Dynamics. Phys. Rev. B 2013, 87, 174110. (16) Magdau, I. B.; Ackland, G. J. High Temperature Raman Analysis of Hydrogen Phase IV from Molecular Dynamics. J. Phys.: Conf. Ser. 2014, 500, 032012. (17) Belonoshko, A. B.; Ramzan, M.; Mao, H.-k.; Ahuja, R. Atomic Diffusion in Solid Molecular Hydrogen. Sci. Rep. 2013, 3, 2340. (18) Liu, H.; Hernandez, E. R.; Yan, J.; Ma, Y. Anomalous Melting Behavior of Solid Hydrogen at High Pressures. J. Phys. Chem. C 2013, 117, 11873−11877. (19) Hernandez, E. R. Metric-Tensor Flexile-Cell Algorithm for Isothermal-Isobaric Molecular Dynamics Simulations. J. Chem. Phys. 2001, 115, 10282. (20) Hernandez, E. R.; Rodriguez-Prieto, A.; Bergara, A.; Alfe, D. First-Principles Simulations of Lithium Melting: Stability of the bcc Phase Close to Melting. Phys. Rev. Lett. 2010, 104, 185701. (21) Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for Liquid Metals. Phys. Rev. B 1993, 47, 558. (22) Kresse, G.; Furthmuller, J. Efficient Iterative Schemes for Ab Initio Total-energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B 1996, 54, 11169. (23) Blochl, P. E. Projector Augmented-Wave Method. Phys. Rev. B 1994, 50, 17953. (24) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (25) Monkhorst, H. J.; Pack, J. D. Special Points for Brillouin-Zone Integrations. Phys. Rev. B 1976, 13, 5188. (26) Krynicki, K.; Green, C. D.; Sawyer, D. W. Pressure and Temperature Dependence of Self-diffusion in Water. Faraday Discuss. Chem. Soc. 1978, 66, 199−208. (27) See Supporting Information for atomic diffusion movie. (28) Zha, C.-s.; Cohen, R.; Mao, H.-k.; Hemley, R. J. Raman Measurements of Phase Transitions in Dense Solid Hydrogen and Deuterium to 325 GPa. Proc. Natl. Acad. Sci. U. S. A. 2014, 111, 4792.

■

CONCLUSION The influence of system size effects on diffusive hydrogen motions in phase IV of solid hydrogen at 300 GPa and 300 K has been extensively examined by first-principles MD simulations with various simulation cells of 96−768 atoms. It is demonstrated that long simulation time is essential to reveal the diffusive events. Our results thus confirmed previous MD results.7,14 It is concluded that long-range diffusive proton motion in phase IV of solid hydrogen is robust. The present study highlights the importance of time span in MD simulation.

■

ASSOCIATED CONTENT

S Supporting Information *

Atomic diffusion movie. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] (Y.M.). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors acknowledge the funding support from China 973 Program under Grant 2011CB808200, National Natural Science Foundation of China under Grants 11274136, 11025418, and 91022029, 2012 Changjiang Scholar of Ministry of Education, and Changjiang Scholar and Innovative Research Team in Jilin University (No. IRT1132). Part of the calculations was performed in high performance computing center of Jilin University.

■

REFERENCES

(1) Mao, H. K.; Hemley, R. J. Ultrahigh-Pressure Transitions in Solid Hydrogen. Rev. Mod. Phys. 1994, 66, 671−692. (2) McMahon, J. M.; Morales, M. A.; Pierleoni, C.; Ceperley, D. M. The Properties of Hydrogen and Helium under Extreme Conditions. Rev. Mod. Phys. 2012, 84, 1607−1653. (3) Eremets, M. I.; Troyan, I. A. Conductive Dense Hydrogen. Nat. Mater. 2011, 10, 927−931. (4) Howie, R. T.; Guillaume, C. L.; Scheler, T.; Goncharov, A. F.; Gregoryanz, E. Mixed Molecular and Atomic Phase of Dense Hydrogen. Phys. Rev. Lett. 2012, 108, 125501. (5) Howie, R. T.; Scheler, T.; Guillaume, C. L.; Gregoryanz, E. Proton Tunneling in Phase IV of Hydrogen and Deuterium. Phys. Rev. B 2012, 86, 214104. (6) Zha, C.; Liu, Z.; Ahart, M.; Boehler, R.; Hemley, R. J. HighPressure Measurements of Hydrogen Phase IV Using Synchrotron Infrared Spectroscopy. Phys. Rev. Lett. 2013, 110, 217402. D

dx.doi.org/10.1021/jp503409p | J. Phys. Chem. C XXXX, XXX, XXX−XXX