Ab Initio Molecular Dynamics Simulation of Liquid Water with Fragment- based Quantum Mechanical Approach under Periodic Boundary Conditions†

a. Department of Basic Medicine and Clinical Pharmacy, China Pharmaceutical University, Nanjing 210009, China b. Shanghai Engineering Research Center of Molecular Therapeutics and New Drug Development, School of Chemistry and Molecular Engineering, East China Normal University, Shanghai 200062, China c. New York University-East China Normal University Center for Computational Chemistry at New York University Shanghai, Shanghai 200062, China


I. INTRODUCTION
Water is a ubiquitous and yet exceptional liquid exhibiting several anomalies. Understanding its properties is key in essentially all disciplines of the natural science and engineering [1]. Although it has been studied extensively for past years, the anomalous properties of liquid water are still not fully explored and have long been under intense scientific debate [2][3][4][5][6][7][8][9][10].
The chemical properties of liquid water are closely related to its complex hydrogen bonding network [11,12].
Ab initio molecular dynamics (AIMD) simulation allows determination of the properties of water with higher accuracy in comparison with the empirical water models [2,10]. However, theoretical predictions of the diverse properties of liquid water are still challenging. At present, most AIMD simulations of liquid water are based on the density functional theory (DFT) [28,29]. Although substantial progress has been reported in the description of the properties of liquid water by using the DFT-based AIMD simulations [29][30][31][32][33], it is still difficult to find a density functional that can provide an equally accurate description of the different types of interactions present in water, including the intramolecular covalent bonds, intermolecular hydrogen bonds, van der Waals interactions as well as the electronic and dielectric properties of the liquid water [16,34]. Consequently, the simulated properties of liquid water from DFT-based AIMD simulations were not uniformly in good agreement with the experimental observations [34]. In addition, due to the high computational cost of QM, AIMD simulations of the condensed-phase matters are usually restricted to small model systems with approximations. Even when achieving affordable computation cost, the simulations using small model systems are improper and difficult to reflect the bulk phase behavior of liquid water [35][36][37][38].
To circumvent the limitations of current AIMD simulations, various computing techniques have been developed for efficiently utilizing high-level ab initio theories on large systems [10,26,[39][40][41][42][43][44]. In our previous works [45][46][47][48][49][50][51][52] we developed the fragment-based QM method for accurate and efficient quantum calculation of large molecular systems to avoid otherwise the prohibitively large computational expenses. By using the quantum fragmentation approach, we have successfully carried out high-level wavefunction theory based AIMD simulations of ambient liquid water [10,53,54]. However, those AIMD simulations were performed in a QM/MM fashion, in which we only focused on a central water cluster (with about 140 water molecules) of the simulation box. In this study, we present AIMD simulations of ambient liquid water under the periodic boundary conditions (PBC) by using the second-order Møller-Plesset perturbation theory (MP2) with the quantum fragmentation approach. Various structural and dynamical properties of liquid water were investigated in comparison with the experimental observations.

A. The PBC-EE-GMF method
The electrostatically embedded generalized molecular fractionation (EE-GMF) method was developed in our group for accurate and efficient QM calculations of condensed-phase molecular systems [10,49,53,55,56]. In the EE-GMF scheme, each individual molecule is assigned as a fragment, and the energy of each fragment and the interaction energy between two fragments that are spatially in close contact are computed by QM, whereas the interaction energies between two distant fragments are treated via pairwise charge-charge Coulomb interactions. All QM calculations are embedded in the electrostatic field of the point charges representing the remaining system to account for the environmental effect. According to the PBC-EE-GMF scheme, the energy per unit cell of the molecular system can be expressed as the following, where n is a three-integer index of a unit cell, E i(n) is the energy of the ith molecule in the nth unit cell and E i(0)j(n) is the energy of the dimer consisting of the ith molecule in the central (0th) unit cell and the jth molecule in the nth unit cell. The first term in Eq.(1) summarizes the energy of each molecule in the central unit cell and the second term is the local twobody QM interactions if the nearest distance between any two molecules in the central unit cell is less than or equal to a predefined distance threshold λ (λ was set to 4Å in this study). For the interactions between the central unit cell and the neighboring unit cells, we also calculate the local two-body QM interactions if any water molecule in the central unit cell has close contact with other water molecules in the neighboring unit cell, as expressed in the third term of Eq.(1). In the case that the nearest distance between two molecules is larger than λ, these long-range interactions are approximately described through charge-charge Coulomb interactions. We have taken into account the background charges within the 3×3×3 supercell (i.e. of Eq.(1) takes into account the long-range interactions within the supercell of dimension as large as 11×11×11 through charge-charge Coulomb interactions as follows, where q η j(n) represents the atomic charge of ηth atom in the jth molecule of nth unit cell and R is the distance between two atomic charges.
The unit cell energy E cell can be analytically differentiated with atomic positions to obtain the atomic forces as follows, All QM calculations were carried out at the MP2/augcc-pVDZ level by using the Gaussian 09 program [57]. The atomic charges utilized for the embedding field were obtained from the flexible simple point-charge water (SPCFW) model [58].

B. AIMD simulation
A cubic unit cell containing 126 water molecules were constructed. Prior to the AIMD simulation, the system was gradually heated from 0 to 300 K in 100 ps, followed by 100 ps equilibration at 300 K with a time step of 1 fs under constant-pressure and constant-temperature (NPT) ensemble by using the classical force field of Am-ber14 [59]. The lengths of the equilibrated unit cell were 15.11, 15.94, and 15.62Å, respectively, with the density of 1.00 g/cm 3 . Subsequently, a total of 13 ps AIMD simulation was carried out under the constant-volume and constant-temperature (NVT) ensemble using a modified version of Amber14. For each step of the AIMD simulation, the energy and atomic forces of the system were calculated using the PBC-EE-GMF approach, and then passed to the MD engine (the Sander module) of Amber14. The AIMD simulation was performed at 300 K with a time step of 1 fs. The Langevin dynamics [60] was applied to regulate the temperature with a collision frequency of 5.0 ps −1 . The last 11 ps simulation trajectory was used for data analysis.

A. RDF
The structural properties of the liquid water can be inferred from the radial distribution function (RDF). FIG. 2 compares the oxygen-oxygen (g OO ), oxygenhydrogen (g OH ), and hydrogen-hydrogen (g HH ) RDFs of liquid water under ambient conditions computed through the PBC-EE-GMF-based AIMD simulation at the MP2/aug-cc-pVDZ level, as compared to the experimental results [61,62]. As can be seen from FIG. 2, this simulation yields g OO in excellent agreement with the experimental observation for the positions of the first three peaks (see FIG. 2(a)). Moreover, the simulated intensities of the first three peaks also match the experiment very well, except that the trough between the first two peaks is slightly higher than the experiment, leading to the intensity of the second peak slightly underestimated.
For the simulated g OH and g HH curves (see FIG. 2  c)), the positions of their first and third peaks are in good agreement with the corresponding experimental results, while the second peaks are slightly shifted to higher radial positions. It is noticeable that the intensities of the first peaks of g OH and g HH are significantly overestimated in comparison with the experimental results, which is mainly due to the lack of the nuclear quantum effect (NQE) in our AIMD simulation. In addition, the intensities of the second and third peaks in these two RDFs are just slightly overestimated. The inclusion of NQEs in AIMD simulation would yield g OH and g HH in closer agreement with the experiment [63], and it is not expected to significantly change g OO [64]. Overall, our AIMD simulated RDFs through the PBC-EE-GMF-based MP2/aug-cc-pVDZ method are in good accordance with the experimental observations.

B. Oxygen-oxygen-oxygen triplet angular distribution
Additional information on the local structure of water can be obtained by computing the distribution of oxygen-oxygen-oxygen angles within the first coordination shell of water molecules. Three oxygen atoms were considered as a given triplet if two of the oxygen atoms were within a prescribed cutoff distance from the third, and this cutoff (3.20Å was used in this study) was chosen to yield an average oxygen-oxygen coordination number of around 4 [29]. The calculated angular distribution, as compared to the experiment [62], is displayed in FIG. 3. As can be seen from FIG. 3, the calculated angular distribution shows a prominent peak near 105 • , which closely matches the experimental observation. In addition, both the experimental and calculated results show a small shoulder near 55 • , with the calculated probability slightly underestimated as com- FIG. 3 Oxygen-oxygen-oxygen angular distribution computed from the PBC-EE-GMF-based AIMD simulation at the MP2/aug-cc-pVDZ level, as compared to the experimental data from Ref. [59].
pared to the experimental result. The position of the prominent peak is close to the perfect tetrahedral angle of 109.5 • , which indicates a tetrahedral arrangement of the first coordination shell of water. However, the small shoulder indicates the distorted arrangement of hydrogen bonds in the first coordination shell. The tetrahedral order parameter q, which is defined as [3], where θ ij is the angle formed between the oxygen atom of a given water molecule and the oxygen atoms i and j of its four nearest neighbors. The tetrahedral order parameter q is calculated to quantify the tetrahedrality of water. A q value of 1 corresponds to a perfect tetrahedral arrangement and 0 corresponds to a completely disordered system [3]. The q value obtained in this study is 0.56, which is very close to the experimental value of 0.57. For comparison, the q value obtained from the DFT-based (using PBE0+TS-vdW(SC)) AIMD simulation was around 0.69 [29], which was significantly deviated from the experiment.

C. Coordination and hydrogen bond dynamics
The properties of liquid water are closely related to its microstructure and dynamical behavior of coordination shell and the formed hydrogen-bond networks. In this study, we quantity the coordination and hydrogen bond number in the first shell of a water molecule. The radius of the first coordination shell utilized in this study is 3.36Å, corresponding to the position of the trough between the first two peaks in g OO . From our AIMD simulation, a water molecule has an average coordination number of 4.7, which is in good agreement with the experimental estimate [61]. FIG. 4(a) plots the fluctuation of the coordination number in the first shells of three representative water molecules as a function of the simulation time. The coordination number displays a surprisingly large fluctuation, ranging from threecoordination to eight-coordination, indicating a remarkable local structure rearrangement around one water molecule. Overall, the four-and five-coordinations are dominant during the MD simulation.
Using a geometrical definition of the hydrogen bond with an O· · · O distance cutoff of 3.5Å and O· · · OH angle cutoff of 30 • , we estimated the number of hydrogen bonds per water molecule to be 3.5. This number has not been exactly determined experimentally, but estimates from previous experiments suggest values between 3 and 4 under ambient conditions [16]. However, a previous AIMD simulation at the DFT level tends to produce over-structured liquid water [29]. FIG. 4(b) shows the time-evolution of the hydrogen-bond number formed in the first shell of three representative water molecules. The hydrogen-bond number also displays a large fluctuation during the MD simulation, which is mainly resulted from the large fluctuation of the coordination number. The number of hydrogen bonds formed around one water molecule ranges from one to six, indicating the different local structures for water molecules. The dominant hydrogen bond number formed per water molecule is from three to four, indicating that the liquid water is a dynamical mixture of the tetrahedral structure and others, such as the "ring-and-chain like" structures [10,65]. where the MSD is obtained by the calculation of the squared relative displacement of the oxygen atoms as a function of time t averaged over all water molecules in the unit cell. The convergence of the calculated D value with respect to the simulation time is shown in FIG. 5. As one can see from the figure, the converged D value in this study is 0.21Å 2 /ps, which is in good agreement with the experimental value of 0.23Å 2 /ps [66]. However, the DFT-based AIMD simulations usually predict a very small D value (∼0.10Å 2 /ps) [29]. The good agreement between our theoretical prediction and experiment demonstrates the accuracy and reliability of our method using high-level wavefunction theory in describing the dynamical properties of liquid water under ambient conditions.

E. Dipole moment
Our PBC-EE-GMF-based AIMD simulation can also reveal the electronic properties of water, e.g., the distribution of the molecular dipole of water, which could not be easily measured via experiments. According to the previous reports [67][68][69][70], the estimated experimental values ranged widely from 2.6 Debye to 3.0 Debye. In this study, the dipole moment of each water molecule was calculated in the electrostatically embedded field of the point charges representing the remaining water molecules. FIG. 6 plots the distribution of the molecular dipole of water in liquid phase. As shown in FIG. 6, the molecular dipole is broadly distributed from 1.8 Debye to 3.2 Debye, reflecting the diversity of the electrostatic environment in liquid water. From this study, the prominent molecular dipole moment is around 2.6 Debye, and this value is in good agreement with the pre- vious reports [68,70,71]. In contrast, the DFT-based AIMD simulation predicted a relatively larger dipole moment (∼2.96 Debye at 300 K) due to the overestimated polarizability of DFT [29].

IV. CONCLUSION
In this work, we carried out ab initio molecular dynamics simulation of liquid water under ambient conditions at the MP2/aug-cc-pVDZ level by using the PBC-EE-GMF fragmentation approach. We studied various properties of liquid water, including the radial distribution functions, oxygen-oxygen-oxygen triplet angular distribution, coordination number and hydrogen bonds in the first shell of the water molecule, diffusion coefficient and molecular dipole moment. All simulated structural and dynamical properties, as well as the electronic properties, were uniformly in good agreement with the existing experimental observations. The results demonstrate the accuracy and reliability of the PBC-EE-GMF approach, allowing us to study the microscopic details of liquid water using high-level wavefunction theories. Therefore, the PBC-EE-GMF method can be efficiently used to assess the validity of uncertain or controversial conclusions previously made on the liquid water, and potentially provide new physical insights into other underexplored anomalous properties of liquid water.