^{1,a)}and Koji Ando

^{1}

### Abstract

A semiquantal (SQ) molecular dynamics (MD) simulation method based on an extended Hamiltonian formulation has been developed using multi-dimensional thawed Gaussian wave packets (WPs), and applied to an analysis of hydrogen-bond (H-bond)dynamics in liquid water. A set of Hamilton's equations of motion in an extended phase space, which includes variance-covariance matrix elements as auxiliary coordinates representing anisotropic delocalization of the WPs, is derived from the time-dependent variational principle. The present theory allows us to perform real-time and real-space SQMD simulations and analyze nuclear quantum effects on dynamics in large molecular systems in terms of anisotropicfluctuations of the WPs. Introducing the Liouville operator formalism in the extended phase space, we have also developed an explicit symplectic algorithm for the numerical integration, which can provide greater stability in the long-time SQMD simulations. The application of the present theory to H-bonddynamics in liquid water is carried out under a single-particle approximation in which the variance-covariance matrix and the corresponding canonically conjugate matrix are reduced to block-diagonal structures by neglecting the interparticle correlations. As a result, it is found that the anisotropy of the WPs is indispensable for reproducing the disordered H-bond network compared to the classical counterpart with the use of the potential model providing competing quantum effects between intra- and intermolecular zero-point fluctuations. In addition, the significant WP delocalization along the out-of-plane direction of the jumping hydrogen atom associated with the concerted breaking and forming of H-bonds has been detected in the H-bond exchange mechanism. The relevance of the dynamical WP broadening to the relaxation of H-bond number fluctuations has also been discussed. The present SQ method provides the novel framework for investigating nuclear quantum dynamics in the many-body molecular systems in which the local anisotropicfluctuations of nuclear WPs play an essential role.

The authors would be grateful to Dr. Kjartan T. Wikfeldt for supplying the RDFs obtained from the RMC method. J.O. would like to thank Professor Shinji Saito for providing helpful comments. K.A. acknowledges supports from KAKENHI Nos. 20108017 (“π-space”) and 22550012.

I. INTRODUCTION

II. GENERALIZED SEMIQUANTAL TIME-DEPENDENT THEORY

A. Correlated multi-dimensional Gaussian wave packet dynamics: Hamiltonian formulation

B. Expectation values for potential functions

III. SYMPLECTIC AND TIME-REVERSIBLE INTEGRATOR

IV. COMPUTATIONAL DETAILS

A. Single-particle approximation

B. Flexible SPC water model

C. Semiquantal molecular dynamics simulations

V. RESULTS AND DISCUSSION

A. Static equilibrium properties

B. H-bond exchange dynamics

C. H-bond number fluctuation

VI. CONCLUSIONS

### Key Topics

- Hydrogen bonding
- 57.0
- Quantum effects
- 26.0
- Molecular dynamics
- 24.0
- Anisotropy
- 23.0
- Intermolecular forces
- 17.0

## Figures

Radial distribution functions (RDFs) for (a) O–O, (b) O–H, and (c) H–H atom pairs of liquid water calculated from the expectation values using the SQ models. Also shown are the RDFs obtained from the classical (CL) model and those fitted to experimental x-ray (Ref. 108) and neutron (Ref. 109) diffraction data in reciprocal space using the reverse Monte Carlo (RMC) method without H-bond constraints (Ref. 110). The insets show the magnifications of the corresponding first peaks.

Radial distribution functions (RDFs) for (a) O–O, (b) O–H, and (c) H–H atom pairs of liquid water calculated from the expectation values using the SQ models. Also shown are the RDFs obtained from the classical (CL) model and those fitted to experimental x-ray (Ref. 108) and neutron (Ref. 109) diffraction data in reciprocal space using the reverse Monte Carlo (RMC) method without H-bond constraints (Ref. 110). The insets show the magnifications of the corresponding first peaks.

Same as Fig. 1 except that the RDFs for the SQ models are calculated with respect to the WP centers.

Same as Fig. 1 except that the RDFs for the SQ models are calculated with respect to the WP centers.

(a) Optimized energies of the water dimer as a function of the O–O distance obtained from the SQ and classical (CL) simulations without intramolecular constraints and (b) those with the fixed intramolecular geometry (see text), where the energy is defined as the deviation from the energy of two isolated water molecules and the O–O distance is calculated with respect to the WP centers. The insets show the magnifications of the corresponding minimum regions.

(a) Optimized energies of the water dimer as a function of the O–O distance obtained from the SQ and classical (CL) simulations without intramolecular constraints and (b) those with the fixed intramolecular geometry (see text), where the energy is defined as the deviation from the energy of two isolated water molecules and the O–O distance is calculated with respect to the WP centers. The insets show the magnifications of the corresponding minimum regions.

Averaged H-bond exchange trajectories of (a) the jump angle θ and (b) the O^{a}O*O^{b} angle ψ. The standard deviations are also plotted as error bars. The insets show the corresponding magnifications.

Averaged H-bond exchange trajectories of (a) the jump angle θ and (b) the O^{a}O*O^{b} angle ψ. The standard deviations are also plotted as error bars. The insets show the corresponding magnifications.

Averaged H-bond exchange trajectories of (a) the O*O^{a} and O*O^{b}, and (b) the O^{a}O^{b} distances. The standard deviations are also plotted as error bars.

Averaged H-bond exchange trajectories of (a) the O*O^{a} and O*O^{b}, and (b) the O^{a}O^{b} distances. The standard deviations are also plotted as error bars.

Averaged H-bond exchange trajectories of the WP width fluctuations of the jumping hydrogen atom H* in the body-fixed axes obtained from the SQ simulation with the EGWP model. Also shown is the graphical definition of the body-fixed axes for the H* atom.

Averaged H-bond exchange trajectories of the WP width fluctuations of the jumping hydrogen atom H* in the body-fixed axes obtained from the SQ simulation with the EGWP model. Also shown is the graphical definition of the body-fixed axes for the H* atom.

(a) Normalized time correlation functions of the H-bond number fluctuation obtained from the SQ simulations with the thawed GWPs (TGWPs) and (b) those with the frozen GWPs (FGWPs) in which only the rotational and translational degrees of freedom of the WPs are allowed to fluctuate. Also shown is the classical (CL) result. The insets show the corresponding power spectra in arbitrary units. Notice the logarithmic scales of the axes.

(a) Normalized time correlation functions of the H-bond number fluctuation obtained from the SQ simulations with the thawed GWPs (TGWPs) and (b) those with the frozen GWPs (FGWPs) in which only the rotational and translational degrees of freedom of the WPs are allowed to fluctuate. Also shown is the classical (CL) result. The insets show the corresponding power spectra in arbitrary units. Notice the logarithmic scales of the axes.

Time series of the deviation of the total Hamiltonian from its initial value in liquid water calculated from the EGWP model using the time-reversible second-order explicit symplectic integrator (SI2) and the fourth-order Runge-Kutta method (RK4) with a time step of Δ*t* = 0.02 fs. The inset shows the time step versus the standard deviation (SD) of the total Hamiltonian in kcal/mol using SI2; here, notice the logarithmic scales of the axes. Also shown in the inset is the linear fitting with a slope of 2.00 (dotted line).

Time series of the deviation of the total Hamiltonian from its initial value in liquid water calculated from the EGWP model using the time-reversible second-order explicit symplectic integrator (SI2) and the fourth-order Runge-Kutta method (RK4) with a time step of Δ*t* = 0.02 fs. The inset shows the time step versus the standard deviation (SD) of the total Hamiltonian in kcal/mol using SI2; here, notice the logarithmic scales of the axes. Also shown in the inset is the linear fitting with a slope of 2.00 (dotted line).

Time development of the intra- and intermolecular potential energies obtained from the EGWP model. The average values (the SDs) of the intra- and intermolecular potential energies are 15.0 (0.23) and −10.5 (0.17) kcal/mol, respectively.

Time development of the intra- and intermolecular potential energies obtained from the EGWP model. The average values (the SDs) of the intra- and intermolecular potential energies are 15.0 (0.23) and −10.5 (0.17) kcal/mol, respectively.

## Tables

Ensemble-averaged static monomer properties in liquid water obtained from the SQ and classical simulations with the f-SPC model. Optimized monomer properties of the isolated water molecule obtained from the SQ simulations are also listed alongside the classical geometric parameters for the f-SPC model. The number in parentheses represents the standard errors in the final digits.

Ensemble-averaged static monomer properties in liquid water obtained from the SQ and classical simulations with the f-SPC model. Optimized monomer properties of the isolated water molecule obtained from the SQ simulations are also listed alongside the classical geometric parameters for the f-SPC model. The number in parentheses represents the standard errors in the final digits.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content