Quantum Molecular Dynamics Simulation of Condensed Hydrogens by Nuclear and Electron Wave Packet Approach

Abstract

We recently proposed a real-time simulation method of nuclear and electron wave packet molecular dynamics (NEWPMD) based on wave functions of hydrogen molecules composed of nuclear and electron WPs. Non-empirical intra- and intermolecular interactions of non-spherical hydrogen molecules were explicitly derived and a long-range dispersion force between hydrogen molecules was successfully derived. Nuclear quantum effects (NQEs) such as nuclear delocalization and zero-point energy were non-perturbatively taken into account. The developed NEWPMD method can be applied to various condensed hydrogen systems from a gas-phase to a high-pressure solid phase at feasible computational cost, providing intuitive understandings of real-time dynamics of each hydrogen molecule including H-H bond vibration, molecular orientations and librational motions. We will report the demonstrations not only in liquid hydrogens but also in solid and supercooled hydrogens which exhibited strong NQEs and thus anomalous static and dynamic properties.

1 Introduction

Although hydrogen is the simplest of all molecular species, its nuclear quantum effects (NQEs) dominate the structure and thermodynamical properties of condensed phases, making it difficult to elucidate them [1,2,3,4,5,6,7]. The NQEs of liquid hydrogen appear as broader radial distribution functions (RDFs), large self-diffusion coefficients even at low temperatures, and anomalous temperature dependence of thermal conductivity and shear viscosity [8,9,10,11,12,13,14,15]. Hydrogen is of primary importance not only for the fundamental properties but also for its potential applications as a renewable energy source without harmful by-products. It is therefore important to establish an a priori computational method to analyze and predict anomalous and unknown properties of condensed phase hydrogens.

A number of quantum mechanical and dynamical approaches have been examined to elucidate this mysterious quantum liquid. A complete treatment of the quantum dynamics for large systems using basis set methods requires storage of multidimensional wave functions and computational effort that grows exponentially with the number of degrees of freedom. Instead, semiquantum methods suitable for numerical simulations of real hydrogen ensembles have been proposed using the path integral Monte Carlo (PIMC) [8,13,16], linearized semiclassical initial value representation (LSC-IVR) [17,18], centroid molecular dynamics (CMD) [9,10,12], ring-polymer molecular dynamics [19], and the thermal Gaussian molecular dynamics [20]. These are all based on the imaginary-time PI theory, and the latter four have been implemented for calculations of time correlation functions (TCFs) accounting for particularly important NQEs such as zero-point energy effect. These previous methods unanimously adopted the Silvera-Goldman model potential [21] that assumed that the hydrogen molecules are a spherically symmetric particle without angular and intramolecular degrees of freedom.

We recently proposed a simulation method of nuclear and electron wave packet molecular dynamics (NEWPMD) based on non-empirical ab initio intra- and intermolecular interactions of non-spherical hydrogen molecules where important NQEs of a hydrogen nucleus were non-perturbatively taken into account [3,4,7,22,23]. The developed NEWPMD method can be applied to various condensed hydrogen systems from a gas-phase to a high-pressure solid phase at feasible computational cost providing intuitive understandings of real-time dynamics of each hydrogen molecule including its H-H bond vibration, molecular orientations and librational motions. In this Account, with the theoretical essence of the NEWPMD, we will report demonstrations not only in liquid hydrogens but also in solid and supercooled hydrogens which exhibit strong NQEs and thus anomalous static and dynamic properties.

2 Three-Dimensional Time-Dependent Wave Packet Approach for Many Body Electrons and Nuclei

The NEWPMD approach describes nuclei by floating and breathing Gaussian WPs [24] via the time-dependent Hartree approach, and EWPs by the perfect-pairing (PP) valence bond (VB) theory. The latter appropriately treats the Pauli exclusion energy and intermolecular long-range dispersion energy which density functional theory without Hartree-Fock-type exchange terms cannot describe [3,22].

2.1 Four-electron and four-nucleus system: two diatomic hydrogen molecules

We started by considering an H₂ dimer and denote the four electrons by a, b, c, and d, and the four nuclei by A, B, C, and D. The nuclei (A, B) and electrons (a, b) form one hydrogen molecule, while the remaining (C, D) and (c, d) compose another hydrogen molecule. The time-dependent total wave function ψ(t), where the EWP pairs (a, b) and (c, d) are coupled in the singlet configuration, was introduced as   

$\begin{array}{c}\psi \left(t\right)=A\left[{\varphi }_{\text{a}}\right(q_1\left){\varphi }_{\text{b}}\right(q_2\left){\varphi }_{\text{c}}\right(q_3\left){\varphi }_{\text{d}}\right(q_4\left)\Theta \right(1,2,3,4\left)\right]\\ \times {\Phi }_{\text{A}}\left(Q_1\right){\Phi }_{\text{B}}\left(Q_2\right){\Phi }_{\text{C}}\left(Q_3\right){\Phi }_{\text{D}}\left(Q_4\right),\end{array}$(1)

where A is an antisymmetrizer for the four electrons [25,26]. For simplicity, $\hslash =1$, an electron charge is unity and all the coordinates are mass scaled. The spin function Θ(1,2,3,4) expresses the PP VB, $\Theta (1,2,3,4)=\left(\alpha \right(1\left)\beta \right(2)-\beta (1\left)\alpha \right(2\left)\right)/\sqrt{2}\times \left(\alpha \right(3\left)\beta \right(4)-\beta (3\left)\alpha \right(4\left)\right)/\sqrt{2}$Each term in ψ(t) is composed of the three-dimensional symmetric Gaussian WP basis functions. The nuclear wave function part was introduced by the Hartree product of the four NWPs each of which has the form   

$\begin{array}{c}{\Phi }_K\left(Q_i\right)\equiv {\left(\frac{1}{2\pi {\Omega }_K^2\left(t\right)}\right)}^{\frac{3}{4}}\times \\ exp\left[X_K{(Q_i-R_K(t\left)\right)}^2+i{P}_K\left(t\right)\cdot (Q_i-R_K(t\left)\right)\right]\end{array}$(2)

with $X_K=(-1+2i{\Pi }_K(t\left){\Omega }_K\right(t\left)\right)/4{\Omega }_K^2\left(t\right)$. The Gaussian NWP is specified by the WP center R_K(t), its width Ω_K(t), and their conjugate momenta P_K(t) and Π_K(t), respectively. Each EWP was defined by   

${\varphi }_k\left(q_i\right)\equiv {\left(\frac{1}{2\pi {\rho }_k^2\left(t\right)}\right)}^{\frac{3}{4}}exp\left[-\frac{{(q_i-r_k(t\left)\right)}^2}{4{\rho }_k^2\left(t\right)}\right],$(3)

whose WP width ρ_k(t) can change depending on the H-H bond length r_HH(t); 0.398+0.360r_HH(t) Å for the large EWP width, and 0.176+0.244r_HH(t) Å for the small EWP width. These linear dependences of ρ_k(t) on r_HH(t) were deduced from the NEWP calculation on a H₂ molecule [22]. We adopted these thawed EWPs to calculate all of the intramolecular potentials which depend only on intramolecular EWPs: ${\overline{J}}_2$, which will be introduced in eq.(10), is composed of such intramolecular potentials. The widths of EWPs were fixed for the intermolecular potentials as ρ_k = 0.654 Å for a and c and 0.370 Å for b and d, respectively, which were determined from the optimized EWPs in a H₂ molecule [22]. Further, from the same optimization of the EWPs in a H₂ molecule, we let the EWP centers r_k(t) depend on the hydrogen nuclear coordinates at each moment as r_a(t)=r_b(t)=(R_A(t)+R_B(t))/2 and r_c(t)=r_d(t)=(R_C(t)+R_D(t))/2. If a hydrogen molecule was highly asymmetrically deformed by a strong external field as seen in a case of polarization, extended formulation would be needed. The EWPs are not evolved variationally, but instead can instantly adjust their widths and positions to NWP dynamics at each moment because timescales of the EWP dynamics are much shorter than timescales of the NWP dynamics. All these approximations greatly simplified the derivations of the equations of motion (EOM) for hydrogen molecules written below.

2.2 Time-dependent variational principle on a pair of hydrogen molecules

Once we introduced the time-dependent total wave function ψ(t), the action integral for the four-electron and four-nucleus system can be defined as   

$\Gamma \equiv \int Ldt=\int dt⟨\psi \left(t\right)|i\frac{\partial }{\partial t}-\hat{H}|\psi \left(t\right)⟩,$(4)

with the Lagrangian L [27,28]. The Hamiltonian operator Ĥ is given by   

$\begin{array}{c}\hat{H}=\sum _{i=1}^4-\frac{1}{2M_{\text{nuc}}}\frac{{\partial }^2}{\partial Q_i^2}\\ +V(q_1,q_2,q_3,q_4;Q_1,Q_2,Q_3,Q_4),\end{array}$(5)

where M_nuc is a relative mass of a proton to an electron. V(q₁, q₂, q₃, q₄; Q₁, Q₂, Q₃, Q₄) is a sum of kinetic energy of electrons V_ke,elec, and three electrostatic potentials of electron-electron, nucleus-nucleus, and nucleus-electron, V_ee, V_nn and V_ne, respectively. The explicit Lagrangian form can be obtained by performing analytical integrations with ψ(t) and Ĥ as   

$\begin{array}{c}L=\sum _K^{\text{A},\text{B},\text{C},\text{D}}[P_K\cdot {\dot{R}}_K+\frac{3}{2}{\Pi }_K{\dot{\Omega }}_K-\frac{3}{2}{\Omega }_K{\dot{\Pi }}_K\\ -\frac{3}{8M_{\text{nuc}}{\Omega }_K^2}-\frac{3{\Pi }_K^2}{2M_{\text{nuc}}}-\frac{P_K^2}{2M_{\text{nuc}}}\\ -⟨V(q_1,q_2,q_3,q_4;Q_1,Q_2,Q_3,Q_4)⟩.\end{array}$(6)

Here, the expectation is denoted by $<\cdots >\equiv <\psi \left(t\right)|\cdots |\psi \left(t\right)>$ The resulting expectation   

$\begin{array}{c}⟨V(q_1,q_2,q_3,q_4;Q_1,Q_2,Q_3,Q_4)⟩=\\ ⟨V_{\text{ke},\text{elec}}⟩+⟨V_{\text{ee}}⟩+⟨V_{\text{nn}}⟩+⟨V_{\text{ne}}⟩,\end{array}$(7)

which is a function of R_K(t) and Ω_K(t), was explicitly and non-perturbatively derived as follows.

The electrostatic energy for the nucleus-nucleus interaction was obtained as   

$⟨V_{\text{nn}}⟩=\sum _{I>J}^{\text{A},\text{B},\text{C},\text{D}}\frac{1}{|R_I-R_J|}\text{erf}\left(\frac{|R_I-R_J|}{2^{\frac{1}{2}}{\left({\Omega }_I^2+{\Omega }_J^2\right)}^{\frac{1}{2}}}\right).$(8)

The other kinetic and potential functions related to the EWPs were given by   

$\begin{array}{c}{V}^{ab,cd}\equiv ⟨V_{\text{ke},\text{elec}}⟩+⟨V_{\text{ee}}⟩+⟨V_{\text{ne}}⟩\\ =\frac{1}{\Delta }(J_0+J_2+J_3+J_4),\end{array}$(9)

where J_n represents an n-electron exchange integral and Δ is a renormalization factor both of which were given in Ref [3]. The electron potential function in eq.(9) includes too many electron exchange integrals to solve the final EOM with reasonable computational cost. We thus removed unimportant terms from the electron exchange integrals as   

$V^{ab,cd}=\frac{1}{\Delta }(J_0+{\overline{J}}_2+{\overline{J}}_3).$(10)

We neglected all the terms in J₄. ${\overline{J}}_2$ includes exchanges only between the intramolecular electron pairs (a,b) and (c,d). ${\overline{J}}_3$ accounts for only the exchange terms related to the electrons a and c of larger Gaussian EWP width in different molecules. The current approximation led to the comparable result to the full result which included all the exchange integrals up to J₄ [3]. Further graphical comparisons of the important intra- and intermolecular potential forms can be seen in Ref [22].

The time-dependent variational principle on the action integral, δΓ/δR_K = 0, etc., resulted in the EOM   

$\begin{array}{c}{\dot{R}}_K=\frac{\partial H_{\text{ext}}}{\partial P_K},{\dot{P}}_K=-\frac{\partial H_{\text{ext}}}{\partial R_K},\\ {\dot{\Omega }}_K=\frac{1}{3}\frac{\partial H_{\text{ext}}}{\partial {\Pi }_K},{\dot{\Pi }}_K=-\frac{1}{3}\frac{\partial H_{\text{ext}}}{\partial {\Omega }_K},\end{array}$(11)

with the extended Hamiltonian,   

$\begin{array}{c}{H}_{\text{ext}}\equiv \sum _K^{\text{A},\text{B},\text{C},\text{D}}\left[\frac{P_K^2}{2M_{\text{nuc}}}+\frac{3{\Pi }_K^2}{2M_{\text{nuc}}}+\frac{3{ħ}^2}{8M_{\text{nuc}}{\Omega }_K^2}\right]\\ +<V_{\text{nn}}>+V^{ab,cd}.\end{array}$(12)

Note that ħ has been retrieved here and it will be hereafter explicitly shown. {R_K, P_K} and {Ω_K, Π_K} can be regarded as conjugate coordinate and momentum pairs. Now dynamics can be described with the potential concept in this extended Hamiltonian. In order to solve the EOM (11), we explicitly and non-perturbatively derived the potential expectation appearing in eq.(12) with the total wave function ψ(t). It is thus distinguished from most of the previous NWP approaches in which a potential surface was given in advance by a separate modeling and, in many cases, expanded quadratically around the moving NWP centers to perturbatively take into account the NQEs [20,27,28,29,30,31,32,33,34,35].

2.3 Many-body hydrogen system

Finally, we extended the above formulation for the four-electron and four-nucleus system to a system composed of many-body hydrogen molecules. The extended Hamiltonian appearing in the EOM (11) for the $N_{\text{mol}}$-body system was derived as   

$\begin{array}{c}{H}_{\text{ext}}\left(N_{\text{mol}}\right)\equiv \sum _K^{N_{\text{nuc}}}\left[\frac{P_K^2}{2M_{\text{nuc}}}+\frac{3{\Pi }_K^2}{2M_{\text{nuc}}}+\frac{3{ħ}^2}{8M_{\text{nuc}}{\Omega }_K^2}\right]\\ +\sum _{I>J}^{N_{\text{nuc}}}\frac{1}{|R_I-R_J|}\text{erf}\left(\frac{|R_I-R_J|}{2^{\frac{1}{2}}{\left({\Omega }_I^2+{\Omega }_J^2\right)}^{\frac{1}{2}}}\right)\\ +\sum _{ab>cd}^{N_{\text{mol}}}{V}^{ab,cd}-\sum _{ab}^{N_{\text{mol}}}(N_{\text{mol}}-2)v^{ab},\end{array}$(13)

where N_nuc and N_mol( = N_nuc/2) are total numbers of nuclei and molecules, respectively. Two electrons a and b(or c and d) constituting one hydrogen molecule should be calculated as a pair in V^ab,cd, leading to the maximum number of summation N_mol. Since the summation of V^ab,cd by all molecular pairs caused multiple counting of the intramolecular electron energies, we needed to subtract an energy of a single hydrogen molecule v^ab whose explicit form was written in Ref [3]. Real-time microscopic trajectories of condensed-phase hydrogen molecules can be simulated using the EOM (11) with H_ext(N_mol) where no empirical parameter for the intra- and intermolecular interactions appears. Since the current EOM representing quantum hydrogens include only the auxiliary coordinates and momenta of the NWPs, many of standard numerical tools for classical MD simulations like thermostats, boundary conditions, parallel computations and sampling methods can be applied. The computational cost for the current NEWPMD is quite reasonable; propagating 576 molecules for 1 ps requires approximately 18 minutes on a single core of a Xeon CPU with 2.67 GHz. The truncations of the electron exchange integrals explained in eq.(10) contributed to reduce the computational cost.

3 Simulation Details

In order to check whether the NEWPMD method describes condensed hydrogen systems correctly, we performed real-time simulations of 1200, 640, 576, and 252 hydrogen molecules in a simulation box with a periodic boundary condition. Figure 1 shows an example of liquid hydrogens taken from the real-time dynamics simulation. We started cooling and equilibration runs from a hexagonal close-packed (h.c.p.) crystal structure where each hydrogen molecule had random orientation. The density was set at the saturated vapor pressure; e.g., the molar volume is 23.06 × 10⁻⁶m³/mol and 32.00 × 10⁻⁶m³/mol at 2.5 K solid and 25 K liquid, respectively [9,36]. The density of the supercooled state was set by extrapolating the saturated vapor-pressure line toward the lower temperature region [9,36]; e.g., the molar volume is 24.20 × 10⁻⁶m³/mol at 2.5 K [23].

Figure 1.

 Snap shot from real-time liquid dynamics of 576 molecules at 25 K.

All integrations of the EOM were performed by the velocity-verlet method with the time step 0.1 fs for the cooling and equilibration runs of hundreds of picoseconds. In the cooling and equilibration runs, we made only the atomic center momentum degrees of freedom, P_K(t), influenced by a heat bath set by a velocity scaling thermostat and Berendsen method with desired temperature. Other degrees of freedom were freely time-evolved by the EOM (11). After the careful cooling and equilibration runs, the whole phase space reached the thermal equilibrium owing to heat conduction between the degrees of freedom controlled by the heat bath and the other degrees of freedom free from the heat bath. After the equilibration runs, we carried out the NVE (microcanonical) simulations for hundres of picoseconds. The increased time step, 0.5 fs, for the NVE simulations gave good total energy conservation; e.g., the ratios of the root-mean-squared (rms) fluctuation of the total energy over the rms fluctuation of the total kinetic energy were 0.04 and 0.02 for 576 molecules at 25 K and 14 K, respectively.

4 Digest of Results

The non-empirical NEWPMD method for the first time enabled systematic investigations and quantitative comparisons of various hydrogen phases. We calculated structural and dynamical properties of condensed-phase hydrogens including their real-time microscopic dynamics of H-H bond vibrations, molecular orientations and librations, and demonstrated the ability of the NEWPMD.

4.1 Liquid hydrogens

The NEWPMD gave the correct structures and transport properties of liquid hydrogen under vapor pressure without any empirical parameters [3,4]. Figure 2 shows RDFs, g(r), of liquid hydrogen at 14 K and 25 K defined by,   

$g\left(r\right)=\frac{⟨n\left(r\right)⟩}{4\pi r^{2}dr{n}_0},$(14)

where n(r) represents the number of atoms between shells of radii r and r+dr, and n₀ is the atomic number density of the whole system [27,28]. The whole shape of the RDF including all the peak positions, height and width almost agreed well with the previous PIMC [8,13,16], LSC-IVR [17] and CMD [9,10,12] calculations and the neutron scattering experiments [11,12].

Figure 2.

 Basic properties of liquid hydrogens at 14 K and 25 K. (a) RDFs obtained by the NEWPMD in comparison with the PIMC results. (b) Time-dependent diffusion properties. Self-diffusion coefficients were obtained by lim_t→∞D(t) according to Einstein's relation.

The shoulder appearing in the first RDF peak at 14 K reflected the diatomic structure of a hydrogen molecule in the first solvation shell; the most stable solvation structure was the T-shape configuration, and the nearest hydrogen atoms around the center hydrogen contributed to the shoulder around 3.25Å. This diatomic shoulder was further supported by the two-dimensional RDF of radial distance and relative molecular orientation θ defined in Figure 2 [4]. The current hydrogen molecules closely interacting each other were not a free spherical molecule due to the asymmetric external interaction potential.

Comparison of the RDFs at 14 K and 25 K indicated that the peaks and valleys of the RDF became broader and deviated toward the larger distance at the higher temperature, while the position of the first peak did not change, in agreement with the previous PIMC and CMD results [9,10,12]. In addition, the shoulder of the first RDF peak was weakened at 25 K, reflecting the less structured liquid.

Figure 2 (b) shows the time-dependent diffusion properties calculated from the slope of the mean-square displacement (MSD) as a function of time,   

$D\left(t\right)=\frac{⟨\left|R\right(t)-R(0){|}^2⟩}{6t},$(15)

where R(t) denotes the position vector of a hydrogen NWP center at time t [27,28]. We calculated the average MSD by averaging over different initial times and hydrogen WP centers with systems of different size. The self-diffusion coefficients of the infinite system obtained by extrapolating the NEWPMD data to N_mol→∞ were D(∞) = 0.345Ų/ps at 14 K and 1.65 Ų/ps at 25 K, in agreement with the corresponding experimental values of 0.35~0.37 Ų/ps at 14 K and 1.4~1.6 Ų/ps at 25 K [3,19,20]. The self-diffusion coefficients became about five times larger by increasing the temperature only by less than twice. This can be explained by combined and correlated effects of the less structured liquid and the shorter H-H bond length at the higher temperature as reported in Ref [4]; since the most stable solvation structure was the T-shape configuration, a H-H bond was more elongated by attraction from electron distributions in bonding region of a paired hydrogen molecule as the liquid structure became more condensed. The average H-H bond lengths, 0.7659 Å at 14 K and 0.7655 Å at 25 K, were reasonably compared to the experimental result, 0.755 Å [22,37]. Reflecting the bond elongation, the H-H bond frequency at 14 K, 4604.7 cm⁻¹, was red-shifted from the frequency at 25 K, 4611.9 cm⁻¹. The current H-H vibrational frequency overestimated the experimental data at 14.2 K, 4151.5 cm⁻¹ [36] by 10.9% which stemmed from the limitation of the Gaussian-type EWPs to describe electrons in each hydrogen molecule. This discrepancy however would be almost irrelevant to the discussions in this Account.

The NEWPMD included the intramolecular degrees of freedom of a non-spherical hydrogen molecule, and can describe the detailed microscopic structures of each molecule even in the condensed phases. In fact, we for the first time discussed distributions, TCFs and frequencies of the H-H bond vibration, intra- and intermolecular orientations, and NWP width in the liquid phase [3,4].

4.2 Solid hydrogens

We successfully reproduced the stable vapor-pressure solid of the h.c.p. lattice structure with the reasonable freezing temperature and lattice phonon frequencies [7]. Figure 3 shows typical h.c.p. structures taken from the vapor-pressure solid dynamics. In Ref [7], we characterized dynamical and structural characters of solid hydrogens under vapor pressure, demonstrating the difference from liquid and high-pressure solid hydrogens [38]. The most stable adjacent configuration exhibited a zigzag structure, in contrast with the T-shape liquid configuration. The periodic angular distributions of hydrogen molecules indicated that molecules were not a completely free rotor in the vapor-pressure solid reflecting the asymmetric potentials from surrounding molecules on adjacent lattice sites. The discrete jumps of librational and H-H vibrational frequencies as well as H-H bond length caused by the structural rearrangements effectively discriminated the liquid and solid phases.

Figure 3.

 Snap shots from real-time solid dynamics of 211 molecules at 2.5 K.

The current thermodynamic points, vapor pressure and temperatures below 25 K, are quite different from thermodynamics states where high-pressure phase I solids have been discussed.(See Figure 8 in Ref [39].) The vapor-pressure hydrogen solid simply connected to the vapor-pressure hydrogen liquid: The blueshift of the H-H vibrational frequency with increasing the temperature simply linked to the density decrease along the saturated vapor pressure line as shownin Figure 5 of Ref [36]. The similar blueshift of the H-H vibrational frequency was observed in the experiments on the vapor-pressure solids and liquids [36]. The high-pressure phase I solid exhibited the completely opposite tendency: the H-H bond length decreased with increasing the pressure up to 30 GPa, indicating simple compression of each hydrogen molecule [37].

4.3 Supercooled hydrogens

Supercooled hydrogen liquid as well as superfluid have continued to elude experimental observation due to rapid crystallization. In Ref [23], we computationally realized and investigated the supercooled hydrogens by the NEWPMD and concluded that the hydrogen supercooled liquid was not a simply cooled liquid but rather exhibited intrinsic structural and dynamical characters including a precursor of tunneling and superfluidity which neither normal hydrogen liquid nor solid possessed. In fact, we demonstrated that the supercooled hydrogen liquid had the distinct structural and dynamical properties such as closer first peak of the RDF than the solid RDF, on-line ar MSD, red-shifting boson peaks with increasing temperature, diffusion caused by mesoscale cluster deformation, and enhanced localization and delocalization of hydrogen nuclei. All of these insights and information we obtained will provide a milestone for planning experiments focusing on hydrogen metastable states and for identifying and characterizing various unknown hydrogen phases.

5 Summary

We reported the first computational study on real-time dynamics of vapor-pressure hydrogen liquids, solids, and supercooled liquids exhibiting strong NQEs which have been hardly accessible by use of previous computational and theoretical methods like density functional theory and the semiquantum molecular dynamics simulations with path integrals. The real-time NEWPMD method is based on the non-empirical intra- and intermolecular interactions of non-spherical hydrogen molecules and thus can be called the quantum molecular dynamics simulation method where any electronic model potentials for the interactions do not need to be given in advance. We confirmed that the NEWPMD with the reasonable computational cost successfully reproduced the experimental basic properties of hydrogen liquid and solid such as structures, transport coefficients, and intra- and intermolecular vibrations, demonstrating that the NEWPMD can be a powerful and efficient simulation method to study nonequilibrium and dynamical processes of various condensed-phase hydrogen systems. These studies will open a new avenue of hydrogen material research and become an essential step toward elucidating mysteries of hydrogen systems.

References

• [1]   H. Saito, H. Nagao, K. Nishikawa, K. Kinugawa, J. Chem. Phys., 119, 953 (2003).
• [2]   M. Zoppi, J. Phys. Cond. Mett., 15, R1047 (2003).
• [3]   K. Hyeon-Deuk, K. Ando, J. Chem. Phys.Commun., 140, 171101 (2014a).
• [4]   K. Hyeon-Deuk, K. Ando, Phys. Rev. B, 90, 165132 (2014b).
• [5]   M. Kühnel, J. M. Fernández, G. Tejeda, A. Kalinin, S. Montero, R. E. Grisenti, Phys. Rev. Lett., 106, 245301 (2011).
• [6]   M. Kuhnel, J. M. Fernandez, F. Tramonto, G. Tejeda, E. Moreno, A. Kalinin, M. Nava, D. E. Galli, S. Montero, R. E. Grisenti, Phys. Rev. B, 89, 180201 (2014).
• [7]   K. Hyeon-Deuk, K. Ando, J. Chem. Phys.Commun., 140, 171102 (2015).
• [8]   D. Scharf, G. J. Martyna, M. L. Klein, Low Temp. Phys., 19, 364 (1993).
• [9]   Y. Yonetani, K. Kinugawa, J. Chem. Phys., 120, 10624 (2004).
• [10]   Y. Yonetani, K. Kinugawa, J. Chem. Phys., 119, 9651 (2003).
• [11]   F. J. Bermejo, K. Kinugawa, C. Cabrillo, S. M. Bennington, B. Fak, M. T. Fernandez-Diaz, P. Verkerk, J. Dawidowski, R. Fernandez-Perea, Phys. Rev. Lett., 84, 5359 (2000).
• [12]   F. J. Bermejo, B. Fak, S. M. Bennington, K. Kinugawa, J. Dawidowski, M. T. Fernandez-Diaz, C. Cabrillo, R. Fernandez-Perea, Phys. Rev. B, 66, 212202 (2002).
• [13]   M. Zoppi, M. Neumann, M. Celli, Phys. Rev. B, 65, 092204 (2002).
• [14]   F. J. Bermejo, B. Fak, S. M. Bennington, R. Fernandez-Perea, C. Cabrillo, J. Dawidowski, M. T. Fernandez-Diaz, P. Verkerk, Phys. Rev. B, 60, 15154 (1999).
• [15]   D. Colognesi, M. Celli, M. Neumann, M. Zoppi, Phys. Rev. E Stat. Nonlin. Soft Matter Phys., 70, 061202 (2004).
• [16]   A. Nakayama, N. Makri, J. Chem. Phys., 119, 8592 (2003).
• [17]   J. A. Poulsen, G. Nyman, P. J. Rossky, J. Phys. Chem. B, 108, 19799 (2004).
• [18]   J. Liu, W. H. Miller, J. Chem. Phys., 128, 144511 (2008).
• [19]   T. F. Miller, D. E. Manolopoulos, J. Chem. Phys., 122, 184503 (2005).
• [20]   I. Georgescu, J. Deckman, L. J. Fredrickson, V. A. Mandelshtam, J. Chem. Phys., 134, 174109 (2011).
• [21]   I. F. Silvera, V. V. Goldman, J. Chem. Phys., 69, 4209 (1978).
• [22]   K. Hyeon-Deuk, K. Ando, Chem. Phys. Lett., 532, 124 (2012).
• [23]   K. Hyeon-Deuk, K. Ando, Phys. Chem. Chem. Phys. Commun., 18, 2314 (2016).
• [24]   A. K. Pattanayak, W. C. Schieve, Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics, 50, 3601 (1994).
• [25]   K. Ando, Bull. Chem. Soc. Jpn., 82, 975 (2009).
• [26]   K. Ando, Chem. Phys. Lett., 523, 134 (2012).
• [27]   K. Hyeon-Deuk, K. Ando, J. Chem. Phys., 132, 164507 (2010).
• [28]   K. Hyeon-Deuk, K. Ando, J. Chem. Phys., 131, 064501 (2009).
• [29]   F. Grossmann, Chem. Phys. Lett., 262, 470 (1996).
• [30]   F. Grossmann, T. Kramer, J. Phys. A, 44, 445309 (2011).
• [31]   E. Heatwole, O. V. Prezhdo, J. Chem. Phys., 126, 204108 (2007).
• [32]   S.-Y. Lee, E. J. Heller, J. Chem. Phys., 76, 3035 (1982).
• [33]   M. H. Beck, A. Jackle, G. A. Worth, H.-D. Meyer, Phys. Rep., 324, 1 (2000).
• [34]   I. Burghardt, K. Giri, G. A. Worth, J. Chem. Phys., 129, 174104 (2008).
• [35]   D. V. Shalashilin, M. S. Child, J. Chem. Phys., 128, 054102 (2008).
• [36]   R. Sliter, A. F. Vilesov, J. Chem. Phys., 131, 074502 (2009).
• [37]   P. Loubeyre, M. Jean-Louis, I. F. Silvera, Phys. Rev. B, 43, 10191 (1991).
• [38]   H. Kitamura, S. Tsuneyuki, T. Ogitsu, T. Miyake, Nature, 404, 259 (2000).
• [39]   J. M. McMahon, M. A. Morales, C. Pierleoni, D. M. Ceperley, Rev. Mod. Phys., 84, 1607 (2012).