Improving ab initio molecular dynamics of liquid water

Abstract

Water is arguably the most important liquid on the earth, yet the physics behind many water properties is still poorly understood. Being able to describe on the fly breaking and formation of chemical bonds, cooperative polarization effects, proton transfer, etc., ab initio molecular dynamics (AIMD) within Kohn-Sham density functional theory (DFT) is an ideal theoretical tool for investigating liquid water structure. To treat water from first-principle quantum mechanical theory, one needs both a good description of the potential energy surface for the nuclei and a quantum mechanical treatment of nuclear dynamics. The nuclear potential energy surface is derived from the quantum mechanical ground-state of the electrons following the Born-Oppenheimer adiabatic approximation of nuclei and electron dynamics.

In this thesis we focus on the improved description of the nuclear potential energy surface that can be achieved by more accurate functional approximations. So far water simulations have used the generalized gradient approximation (GGA). This approximation is affected by the so-called self-interaction error that causes an excessive delocalization of the protons in hydrogen bonds (H-bonds). Moreover, it neglects the non-local correlations responsible for van der Waals (vdW) interactions. These interactions are weaker than H-bonds, but a number of recent papers and the work discussed in this thesis show that vdW interactions play a sizeable role in the water structure.

We adopt the hybrid functional approximation which mixes a fraction of exact exchange into the GGA functional to mitigate the spurious self-interaction error, leading to a better description of H-bond and vibrational properties. In addition, we adopt a non-local functional approximation that accounts for vdW interactions, which are crucial to describe correctly the structure and the equilibrium density of water. Both hybrid functional and non-local correlation are found to be indispensable, resulting in an improved liquid water structure compared to the GGA approximation. The improvement is manifest in several properties, including the radial distribution functions, the bond angle distribution, the broken H-bond statistics, etc.

The statistical effects of quantized nuclear dynamics can be treated in an exact numerical way through path-integral (PI) AIMD simulation. But this calculation is computationally expensive. Due to the large computational cost of the improved DFT functional used in this study, we adopt a cheap approximation for quantum nuclear motion. According to this approximation, the nuclei move classically at a temperature ~ 30K higher than desired temperature. This choice was motivated by a recent work comparing PI AIMD simulation with classical AIMD simulation at ambient conditions with the same DFT functional. This work found that the oxygen-oxygen radial distribution function resulting from classical trajectories mimicked closely the corresponding quantity from the path-integral study when the temperature of the classical simulation was ~ 30K higher than that of quantum simulation. This approximate treatment of quantum nuclei restricts our investigation to the analysis of the oxygen distribution functions, as observables that depend more directly on the proton distribution are more strongly quantum mechanical. We find that the oxygen distribution resulting from our improved DFT approximations agrees well with the available experimental data.

An analysis of the O site distribution using an order parameter that differentiates between sites representative of a high density liquid (HDL) environment and sites representative of a low density liquid (LDL) environment is particularly illuminating, as it allows us to better understand the competing role of H-bonds and vdW interactions in liquid water. This analysis shows that HDL and LDL sites are simultaneously present in the liquid at ambient conditions, consistent with the idea that the inherent potential energy surface of the liquid has two dominant conformations, with low and high density respectively.