Coupled Cluster Molecular Dynamics of Condensed Phase Systems Enabled by Machine Learning Potentials: Liquid Water Benchmark

Coupled cluster theory is a general and systematic electronic structure method, but in particular the highly accurate “gold standard” coupled cluster singles, doubles and perturbative triples, CCSD(T), can only be applied to small systems. To overcome this limitation, we introduce a framework to transfer CCSD(T) accuracy of finite molecular clusters to extended condensed phase systems using a high-dimensional neural network potential. This approach, which is automated, allows one to perform high-quality coupled cluster molecular dynamics, CCMD, as we demonstrate for liquid water including nuclear quantum effects. The machine learning strategy is very efficient, generic, can be systematically improved, and is applicable to a variety of complex systems.

Figure 1

Radial distribution functions of water including nuclear quantum effects from path integral CCMD simulations at 298 K (red dashed lines) compared to x-ray diffraction (a) and neutron scattering (b), (c) data (black solid lines) with their error bars (gray-shaded areas) [94, 95]; see SM, Sec. S-V.A [40] for details on the presented OO RDF including information on the error bars. Only the intramolecular OH peak in panel (b) at $\approx 1\AA$ (that is also presented in its inset) is given by the $\mathrm{\Delta }{G}_{\mathrm{OX}}(r)$ data in Fig. 2 obtained from O isotope substitution neutron diffraction [97] experiments at $300.6\pm 0.5\mathrm{K}$ (black solid line without error bars, digitally extracted) which approximates $g_{\mathrm{OH}}(r)$ as explained therein, whereas the intermolecular OH and HH RDFs are from [95]. The RDFs have been computed using a radial bin size of 0.01 Å as depicted and are not smoothened.

Figure 2

Mean-square displacement (a) and $C_2^{\mathrm{OH}}$ orientational correlation function (b) of water including nuclear quantum effects from path integral CCMD simulations at 298 K. The insets magnify the long (a) and short (b) time behavior from zero up to 50 and 0.5 ps, respectively. The functions have been computed using a time interval of 1 fs as depicted and are not smoothened.