Molecular dynamics (MD) simulations can provide useful insights into the atomic details of water. There is a long standing effort to simulate water using MD methods. The driving forces for MD can be obtained from either empirical force-field or ab initio calculations. Many water models have been proposed such as the TIPnP family , SPC , and MB-Pol . It is also common to use density functional theory (DFT) for ab initio molecular dynamics. However, despite the fact that the DFT formalism is exact for ground-state properties in principle, in practice approximations must be made in the exchange-correlation energy. In principle, the accuracy increases according to the "Jacob's ladder" of electronic structure theory, with the local density approximation (LDA), generalized gradient approximations (GGAs), meta-GGAs, hybrid functionals, and beyond e.g. double hybrids, random phase approximation (RPA), and correlated wave function methods, such as second-order Møller-Plesset perturbation theory (MP2). A properly chosen exchange-correlation functional has been found to predict reasonable binding energy for the water dimer, trimer, larger clusters, and ice as well as the structure of liquid water when combined with MD. LDA is unacceptable as it overbinds the water dimer by nearly a factor of 2. Unsurprisingly, LDA predicts overstructured liquid water and a very small diffusion constant. GGA functionals perform only slightly better. DFT-based MD at the GGA level, without dispersion corrections, often requires a much higher temperature than 300 K to maintain water in the liquid state. For example, one of the most popular functionals, PBE, can reproduce experimental radial distribution functions (RDFs) and density at T=440 K and P = 0.3 GPa . Even when a dispersion correction is added a higher temperature is still required to reproduce these properties. Grimme's dispersion correction to PBE does not significantly improve the structure of water. BLYP with this correction requires simulation at 360 K to reproduce the correct RDF. Among the GGA functionals, revPBE-D3 tends to perform best for RDF, density, and diffusion coefficient at ambient conditions without increasing the temperature. However, the results are very sen-sitive to the particular choice of dispersion functional. The revPBE+DRSLL functional tends to overestimate the volume of ice VIII by 20%. Recent developments in computational power and methodology have made it possible to climb higher rungs of the ladder beyond GGA, namely meta-GGA, hybrid functionals, and many-body correlated methods. Recent works show that significant improvement can be achieved by the non-empirical meta-GGA functional SCAN, as developed by Perdew and co-workers. SCAN yields the correct ordering of densities between liquid water and ice, at the same time predicting quantitative agreement with experiments at an elevated temperature of 330 K. However, the bare SCAN functional systematically overbinds water . Moreover, adding rVV10 to SCAN makes overbinding even stronger, and predicts even higher water density. Fortunately, a very recent study by Paesani and co-workers found that density-corrected SCAN improved the accuracy to the same degree as the "gold standard" coupled-cluster theory. At the same time, some empirically parametrized meta-GGA functionals such as B97M with rVV10 correction appear to perform quite well. At hybrid functional level, The revPBE0-D3 functional is able to predict the correct experimental RDF and density at room temperature. However, revPBE0-D3 still underestimates the temperature-dependent density of water by about 5% , which may be attributed to the choice of van der Waals interactions. When the hybrid functional is combined with rVV10 van der Waals interactions, the equilibrium density of water agrees better with the experiment when the parameter from the rVV10 van der Waals interaction is finetuned. Fortunately, it is possible to go beyond hybrid functionals; including virtual orbitals allows for long-range van der Waals dispersion interactions from parameter-free abinitio calculations. With tremendous computational cost, one can reach the fifth rung of the ladder with methods such as RPA and MP2. The MP2 theory gives excellent predictions of the density at room temperature, while the calculated radial distribution functions are in fair agreement with experimental data that tend to be attributed to the stronger dispersion interaction. Previous studies suggest that "MP2 water" is denser and cooler at ambient conditions compared both with experiment and with "DFT water". In addition to the underlying electronic structure theory, another issue is the accurate account of the quantum nature of the nuclei. Neglect of nuclear quantum effects (NQEs) may become one of the largest sources of error. For instance, there are subtle differences in RDFs upon the inclusion of NQEs, especially those that involve hydrogen atoms. In a classical simulation, to produce an equivalent change in the RDF, one has to increase the temperature by about 30 K. Although different quantum effects on the O-O RDF can be canceled, quantum effects on individual H-bond coordinates can be significant. In comparison with a classical simulation, proton fluctuations along a covalent bond direction increase almost 10-fold. In addition, recent work from Voth and co-workers found that higher temperatures do not accurately replicate NQEs at room temperature, which is evident in different three-body correlations as well as dynamics. Furthermore, classical molecular dynamics overestimates the vibrational frequencies of the O-H stretch region and the H-O-H bend region compared to results that take the NQEs into account. Regarding the reactivity of water, "classical" water is more basic than "quantum" water (i.e., water in nature) with a pH of 8.6, which is about a 30-fold change in the ionization constant. Similar effects have also been found at metal-water interfaces and for electron transfer in aqueous solutions . Hence, the inclusion of nuclear quantum effects is often necessary to achieve high accuracy in molecular dynamics simulations.