# From Complexity to Clarity: Kolmogorov-Arnold Networks in Nuclear Binding
Energy Prediction

Hao Liu  Jin Lei<EMAIL_ADDRESS>Zhongzhou Ren  School of Physics Science
and Engineering, Tongji University, Shanghai 200092, China.

###### Abstract

This study explores the application of Kolmogorov-Arnold Networks (KANs) in
predicting nuclear binding energies, leveraging their ability to decompose
complex multi-parameter systems into simpler univariate functions. By
utilizing data from the Atomic Mass Evaluation (AME2020) and incorporating
features such as atomic number, neutron number, and shell effects, KANs
demonstrated significant predictive accuracy, surpassing traditional models.
The symbolic regression analysis yielded simplified analytical expressions for
binding energies, aligning with classical models like the liquid drop model
and the Bethe-Weizsäcker formula. These results highlight KANs’ potential in
enhancing the interpretability and understanding of nuclear phenomena, paving
the way for future applications in nuclear physics and beyond.

###### keywords:

Kolmogorov-Arnold Networks (KANs) , Nuclear Binding Energies , Symbolic
Regression

††journal: Physics Letters B

##  1 Introduction

The atomic nucleus, a quintessential quantum many-body system, exhibits
remarkable structural complexity [1]. Binding energies (BE), which are
synonymous with nuclear mass, are key characteristics of atomic nuclei and
play an indispensable role in elucidating various nuclear phenomena. These
phenomena include nuclear shapes, shell effects, pairing effects, and the
emergence and disappearance of magic numbers [2, 3, 4, 5]. Furthermore,
binding energies are crucial in the synthesis of superheavy elements and in
understanding nuclear astrophysical processes like the r-process [6] and x-ray
bursts [7]. Therefore, both theoretical predictions and experimental
measurements of nuclear binding energies are essential for advancements in
nuclear physics.

Over the years, numerous theories and methods have been developed for
predicting nuclear binding energies. These include the Bethe-Weizsäcker
formula [8, 2], the Thomas-Fermi model [9], the Hartree-Fock-Bogoliubov mean
field model [10], and ab-initio methods [11, 1, 12]. These physics-based
models have achieved commendable results within their respective applicable
ranges. Concurrently, data-driven methods such as Multi-Layer Perceptrons [13,
14, 15, 16, 17, 18], Gaussian processes [19, 20], and support vector machines
[19] have also found successful applications in the prediction of atomic
nucleus masses.

However, traditional physical models necessitate a profound understanding of
the inherent mechanisms of nuclear physical systems and face challenges when
managing complex relationships. On the other hand, while data-driven machine
learning methods can handle complex non-linear relationships, they often
require substantial data and computational resources. In the data-rich field
of nuclear physics, identifying intricate relationships among variables is a
formidable challenge. This predicament underscores the need for innovative
approaches to handle the complexity of nuclear data and to extract meaningful
insights.

Kolmogorov-Arnold Networks (KANs) offer a promising solution [21]. KANs
provide a simplified yet robust approach by decomposing complex multi-
parameter systems into manageable univariate functions. This approach is based
on the Kolmogorov-Arnold representation theorem [22, 23], which asserts that
any multivariate continuous function can be represented as a superposition of
continuous univariate functions over a compact domain. This theoretical
foundation allows KANs to handle complex, multi-parameter systems by breaking
them down into simpler, univariate components.

One of the most compelling features of KANs in this scientific context is
their capability for symbolic regression. Symbolic regression is a form of
regression analysis that seeks mathematical expressions, in symbolic form,
that best fit a given dataset. Unlike traditional numeric prediction models
that offer a fixed structure and limited insights into underlying patterns,
symbolic regression with KANs provides formulaic expressions that reveal the
inherent relationships and governing laws within the data. This significantly
enhances our understanding and offers a more profound and interpretable model.

The ability to derive symbolic equations directly from data can guide
subsequent scientific research by offering clearer insights into the
relationships and mechanics at play in nuclear phenomena. Such insights are
invaluable for advancing theories, refining experimental designs, and
developing new technologies in nuclear science. Hence, KANs do not merely
predict but also illuminate the pathway towards a deeper comprehension of the
intricate dynamics that define nuclear structures.

In the broader context of artificial intelligence (AI) applications in nuclear
physics, there is a growing interest in leveraging machine learning to study
the intricate phenomena of nuclear physics and derive formulas directly from
data. The liquid drop model, which provides a macroscopic description of the
nucleus, has been instrumental in understanding nuclear binding energies. AI
techniques, particularly those involving deep learning and symbolic
regression, have the potential to learn from vast amounts of nuclear data and
derive similar or even more refined mass formulas. These data-driven
approaches can capture subtle patterns and correlations that might be
overlooked by traditional models, thus offering a complementary perspective to
theoretical frameworks. A significant advantage of AI for physics is its
ability to uncover hidden relationships within complex datasets. In the
current study, we aim to derive mass formulas using AI, harnessing its power
to reveal insights that can enhance our understanding of nuclear binding
energies.

In this study, we aim to explore the potential of KANs in predicting nuclear
binding energies. By leveraging the capabilities of KANs, we seek to bridge
the gap between data-driven insights and theoretical understanding, ultimately
contributing to the advancement of nuclear science. The paper is organized as
follows: In Sec. 2, we provide a brief overview on the Kolmogorov Arnold
Network and data collection and feature seclection used in current study. Sec.
3 presents the application of KANs to nuclear binding enegry. Finally, we
conclude with a summary and discusses in Sec. 4.

##  2 Formalism

###  2.1 Kolmogorov Arnold Network

The Kolmogorov-Arnold theorem asserts that any continuous multivariate
function on a bounded domain can be decomposed into a finite composition of
continuous univariate functions and addition. Specifically, for a continuous
function $f:[0,1]^{n}\rightarrow\mathbb{R}$:

$f(\mathbf{x})=f(x_{1},x_{2},\ldots,x_{n})=\sum_{q=1}^{2n+1}\Phi_{q}\left(\sum_{p=1}^{n}\phi_{q,p}(x_{p})\right).$
(1)

Here, $\mathbf{x}=(x_{1},x_{2},\ldots,x_{n})$ is an $n$-dimensional vector
with each $x_{i}$ in $[0,1]$. The functions
$\phi_{q,p}:[0,1]\rightarrow\mathbb{R}$ and
$\Phi_{q}:\mathbb{R}\rightarrow\mathbb{R}$ are continuous univariate
functions. This decomposition shows that addition is the only truly
multivariate operation, as all other multivariate functions can be built using
univariate functions and sums. This is useful for identifying learnable and
explainable patterns in data. However, practical challenges in machine
learning arise due to the potential non-smooth and fractal nature of the
functions $\phi_{q,p}$ and $\Phi_{q}$, complicating decomposition and learning
processes.

Recent research by Liu et al. [21] has renewed interest in the Kolmogorov-
Arnold theorem and its use in neural networks. This has led to the creation of
Kolmogorov-Arnold Networks (KANs). Liu et al. have improved the original
framework by allowing networks to have any number of layers and widths. This
overcomes the original theorem’s limitations, which restricted it to two
layers of non-linearities and a fixed number of terms $(2n+1)$. Since most
real-world functions are usually smooth and have simple structures, these
improvements help create more practical and efficient Kolmogorov-Arnold
representations.

Figure 1:  Architecture of a shadow KAN, featuring two input features, a
hidden layer with five nodes, and a single-node output layer.

In Fig. 1, a shadow KAN is illustrated, consisting of two input variables,
$x_{1}$ and $x_{2}$, a hidden layer with five nodes, and a single output. Each
input is connected to the hidden nodes through edges characterized by
adaptable functions, replacing the traditional fixed activation functions used
in Multi-Layer Perceptrons. These adaptable functions are expanded into basis
functions, specifically B-spline functions, which dynamically adjust based on
the input data. The outputs from the hidden nodes are then processed through
additional spline-based edges, culminating in a summation at the output node.
This structure allows for increasing the network depth, enabling the analysis
of more complex physical phenomena.

###  2.2 Data Collection and Feature Selection

In this study, we used mass excess values from the Atomic Mass Evaluation
(AME2020) [24]. We focused on nuclei with both atomic number (Z) and neutron
number (N) greater than or equal to 8, covering 3456 nuclei. This selection
ensures the model’s accuracy and stability for larger, more complex nuclei.

To avoid overfitting and make the best use of our data, we applied k-fold
cross-validation [25, 26, 27], which ensures that every data point is used for
both training and validation, maximizing the use of available data. First, we
randomly shuffled the data, then split the dataset into $k$ parts (folds). In
each iteration, $k-1$ parts are used for training, and the remaining part is
used for validation. This process repeats $k$ times, with each part serving as
the validation set once. Averaging the performance across all iterations gives
a reliable measure of the model’s generalizability and helps prevent
overfitting. In our calculations, $k$ is set to 3.

We carefully selected data attributes to ensure the framework’s effectiveness
and precision. These include the atomic number (Z), neutron number (N), and
mass number (A). We also included features from previous research [19, 13,
28], such as pairing effects, shell closures, and asymmetry, which are crucial
for modeling nuclear binding energies.

Pairing effects are captured using $Z_{\text{EO}}$ and $N_{\text{EO}}$, which
can only be either 1 or 0. Here, 1 indicates odd $Z/N$ ratios, and 0 indicates
even ratios. This distinction helps the model differentiate between various
types of nuclei, each possessing unique binding energy characteristics.
Additionally, the term $A^{2/3}$ represents the surface effect, reflecting
that a nucleus’s surface area scales with the two-thirds power of its mass
number ($A$). This scaling is crucial for accurately modeling the surface
contributions to the nucleus’s overall binding energy.

The isospin asymmetry, $(N-Z)$, measures the difference between neutrons and
protons. We also consider nuclear magic numbers, with $Z_{m}$ and $N_{m}$
representing the excess numbers of protons and neutrons relative to the
nearest closed shell. Key magic numbers include
$Z_{m}(N_{m})=8,20,28,50,82,126,184$.

We also explore nuclear shells with $Z_{\text{shell}}$ and $N_{\text{shell}}$,
reflecting the orbital configurations of the last proton and neutron. These
values range from 0 to 5, depending on proton or neutron counts within
specific intervals. This helps the model capture the shell structure’s
influence on binding energy.

We tested different feature combinations to find the most informative
attributes for predicting nuclear binding energies. This approach helps
identify the best feature set for our KAN framework.

To improve model performance and interpretability while reducing overfitting,
we used L1 regularization and cross-entropy loss as described in Ref. [21].
The loss function is defined as:

$\text{Loss}=\frac{(\mathrm{BE}_{\text{pred}}-\mathrm{BE}_{\text{exp}})^{2}}{N_{p}}+\lambda\left(\mu_{1}\sum_{l=0}^{L-1}\left|\bm{\Phi}_{l}\right|_{1}+\mu_{2}\sum_{l=0}^{L-1}S\left(\bm{\Phi}_{l}\right)\right),$
(2)

where $\left|\bm{\Phi}_{l}\right|_{1}$ is the L1 norm of the KAN layer,
defined as the sum of the L1 norms of all activation functions $\phi_{i,j}$,
including $n_{\text{out}}$ outputs and $n_{\text{in}}$ inputs. The L1 norm of
an activation function, averaged over $N_{p}$ inputs, is given by:

$|\phi|_{1}\equiv\frac{1}{N_{p}}\sum_{s=1}^{N_{p}}\left|\phi\left(x^{(s)}\right)\right|.$
(3)

Thus, the L1 norm of a KAN layer is:

$|\bm{\Phi}|_{1}\equiv\sum_{i=1}^{n_{\text{in}}}\sum_{j=1}^{n_{\text{out}}}\left|\phi_{i,j}\right|_{1}.$
(4)

The entropy of a KAN layer is defined as:

$S(\bm{\Phi})\equiv-\sum_{i=1}^{n_{\text{in}}}\sum_{j=1}^{n_{\text{out}}}\frac{\left|\phi_{i,j}\right|_{1}}{|\bm{\Phi}|_{1}}\log\left(\frac{\left|\phi_{i,j}\right|_{1}}{|\bm{\Phi}|_{1}}\right).$
(5)

L1 regularization helps make the model simpler by adding a penalty based on
the absolute values of the model’s parameters. This leads to models with fewer
non-zero parameters, making them easier to interpret and less likely to
overfit. Cross-entropy loss ensures the model’s predictions are well-
calibrated by measuring the difference between predicted and actual
probability distributions.

Combining these techniques, we aim to create a robust and understandable model
for predicting nuclear binding energies. We used k-fold cross-validation [25,
26, 27] to test different parts of the data, selected features carefully, and
applied advanced regularization techniques to ensure our KAN framework is
accurate and generalizable. By exploring different features and network
structures, we gain insights into what affects model performance, helping us
develop better neural network architectures for nuclear physics.

Table 1: Feature space and KANs’ structure Model | Feature
---|---
KAN-2 | N, Z
KAN-4 | N, Z, A, N-Z
KAN-9 | N, Z, A, N-Z, $A^{2/3}$, $Z_{\text{EO}}$, $N_{\text{EO}}$, $Z_{m}$, $N_{m}$
KAN-11 | N, Z, A, N-Z, $A^{2/3}$, $Z_{\text{EO}}$, $N_{\text{EO}}$, $Z_{m}$, $N_{m}$, $Z_{\text{shell}}$, $N_{\text{shell}}$

##  3 Results

Figure 2:  The absolute value of binding energy differences between KAN
predictions using different features (see Table 1) and AME2020 [24]. The RMSEs
for all the data are also provided.

Figure 3:  Experimental single neutron separation energies for $\mathrm{Ca}$,
$\mathrm{Zr}$, $\mathrm{Nd}$, and $\mathrm{Hg}$ isotopic chains in comparison
with the KAN predictions.

Figure 2 shows the absolute differences in binding energy between predictions
from various KAN models and the experimental data from AME2020 [24]. The
feature sets used in these models are listed in Table 1. Notably, the KAN-2
model, which uses only basic properties such as neutron number $N$ and proton
number $Z$, achieves a root mean square deviation error (RMSE) of just 0.92
MeV. This performance surpasses that of most microscopic models [19],
demonstrating the KAN approach’s effectiveness with minimal input features.
When additional features are incorporated, the KAN-11 model reduces the RMSE
even further to an impressive 0.26 MeV.

As more features, such as isospin asymmetry and mass number, are included, the
RMSE decreases, as illustrated in Figure 2, progressing from KAN-2 to KAN-4.
However, it is evident that in both models, the deviation error significantly
increases as nuclei approach closed shells. This error distribution
underscores the influence of shell effects on nuclear binding energy. By
incorporating features related to shell structure, models KAN-9 and KAN-11
achieve significantly improved prediction accuracy. Figure 2 demonstrates that
these models perform better for medium and heavy nuclei compared to light
nuclei. Despite efforts to enhance predictions by expanding the model space,
the challenge of accurately modeling light nuclei persists.

Figure 3 compares the separation energies for isotopic chains of four
different elements against the values in AME2020. The absence of pairing terms
in KAN-2 and KAN-4 leads to much larger prediction errors compared to KAN-9
and KAN-11, which include these terms. This difference is especially
noticeable in the color variations representing the absolute differences in
binding energy between the predictions from various KAN models and the
experimental data from AME2020, as shown in Figure 2. These variations are
present in KAN-2 and KAN-4 but not in KAN-9 and KAN-11, highlighting the
crucial role of pairing effects in detailing nuclear structure. The single-
neutron separation energy curves in Figure 3 further demonstrate the impact of
pairing. The difficulties faced by KAN-2 and KAN-4 models with abrupt changes
suggest that adding additional features like $N_{EO}$ and $Z_{EO}$ could
better model the pairing behavior in nuclear binding energies, leading to
improved performance in KAN-9 and KAN-11.

However, small differences remain in the calcium ($\mathrm{Ca}$) isotopic
chain, as shown in Figure 3. The model shows larger errors for light nuclei
compared to medium and heavy nuclei. Light nuclei are more affected by
interactions among a few particles, while heavier nuclei are influenced by
interactions among many particles. This difference makes it challenging to
create a model that works well for all types of nuclei. It suggests that we
need a more advanced KAN model to better understand the complex behavior of
nuclear systems across different mass regions.

###  3.1 Feature Analysis

Figure 4: Visualization of the trained network structure in the panel (a). The
transparency of the lines connecting nodes is proportional to $\tanh(\beta
A_{l,i,j})$, where $A_{l,i,j}$ is the mean activation, enhancing the
visibility of connections. In panel (b), the $L1$ norm for each connection is
shown on a logarithmic scale, highlighting the significance of each link.

In Fig. 4, we show the structure of the trained KAN. To clarify node
connections, we adjust the transparency of lines based on the mean value of
the activation function $\phi_{l,i,j}$ for each connection, represented by
$\tanh(\beta A_{l,i,j})$, where $A_{l,i,j}$ is the mean activation. This
technique highlights each variable’s influence on the final output, aiding in
understanding the network’s behavior. Additionally, Fig. 4(b) illustrates the
$L1$ norm for each connection, providing insights into the importance of each
link and identifying significant contributors to the model’s performance. This
dual visualization of mean activation and $L1$ norm offers a comprehensive
view of the network’s workings and the relative importance of different
features and connections.

Our analysis shows that the nucleon number $A$, representing the volume term,
has the most significant influence on the results, aligning with intuitive
expectations and the Bethe-Weizsäcker formula. The neutron number $N$ and
proton number $Z$ also play substantial roles, reflecting their fundamental
contributions to nuclear properties.

The asymmetry term $N-Z$ and the surface term $A^{2/3}$ have moderate effects,
indicating their relevance but lesser dominance compared to $A$, $N$, and $Z$.
These terms are essential for capturing the nuances of nuclear structure,
although their impact is not as pronounced as the primary variables.

Variables associated with nuclear shell structure and the pairing term show
minimal direct impact on the overall results. Despite their subtle influence,
these factors are crucial for precise predictions of nuclear binding energies,
particularly in astrophysics and nucleosynthesis, where high accuracy is
paramount. Figures 2 and 3 demonstrate that incorporating these terms, though
they contribute less directly, is essential for accurate modeling.

These findings are consistent with our calculations and corroborate results
from other studies [13, 19, 28, 15, 29], reinforcing our approach’s validity
and our model’s robustness. This analysis highlights the primary drivers of
nuclear binding energy and emphasizes considering all relevant factors, even
those with less apparent impact, to accurately capture the complexity of
nuclear interactions.

###  3.2 Symbolic regression

To create a simple model, we aim to reduce the complexity introduced by
periodic functions, such as sine functions, which can be affected by the
pairing term. To minimize this effect, we focus only on even-even nuclei.
Additionally, as discussed earlier, we limit our study to medium to heavy
nuclei, specifically those with $Z\geq 30$ and $N\geq 30$.

We initially trained a KAN network with the following architecture: 2 input
neurons, 1 hidden layer with 5 neurons, and 1 output neuron. To maintain
accuracy while increasing the regularization coefficient, we enhanced network
sparsity by pruning connections that had minimal impact. This pruning process
simplified the network to a structure with 2 input neurons, 1 hidden layer
with 2 neurons, and 1 output neuron. The final model achieved a RMSE of 0.66
MeV. Using this streamlined network, we then performed symbolic regression.

Figure 5: Panel (a) shows the original trained network structure, while Panel
(b) shows the symbolic network structure. Panels (c) and (d) illustrate the
absolute value of binding energy differences between predictions and data from
AME2020 [24].

However, it is important to note that fitting each activation function
inevitably introduces some error. These errors can accumulate and propagate
through the network, resulting in larger overall discrepancies compared to the
original errors from the activation functions. Consequently, the symbolic
regression results may exhibit greater errors relative to the initial model.
Despite this, the derived expressions can still reveal and confirm underlying
physical laws.

As shown in Fig. 5, the transition from the original network’s activation
functions (the black line in Fig. 5(a)) to the symbolic network’s activation
functions (the red line in Fig. 5(b)) results in a loss of detail in the
activation functions. These lost details propagate through the network,
leading to varying errors. A comparison of the differences between the
original and symbolic results with experimental data (Fig. 5(c) and (d))
reveals underlying patterns: near magic numbers 50, 82, and 126, the errors in
the symbolic results are notably larger than those in the original results.
This suggests that the details not fully captured by symbolic fitting in the
activation functions reflect structural information about the nuclei. Despite
not including variables related to shell structure in our input features, the
KAN network demonstrated a robust capability to identify features from the
data.

Below is the expression obtained through symbolic regression. After
simplification, it can be presented in the following form:

$\displaystyle BE(N,Z)=$ $\displaystyle A_{vol}A-A_{sym}(N-Z)-A_{coul}Z^{2}$
(6) $\displaystyle(A_{Z}(N+n)^{2/3}-A_{N}(Z+z)^{2/3}+C_{2})^{2}-C_{1}$

Table 2: Parameter and its value for the Eq. 6 Parameter | Value | Parameter | Value
---|---|---|---
$A_{vol}~{}(\mathrm{MeV})$ | 9.87917 | $A_{Z}~{}(\mathrm{MeV^{1/2}})$ | 0.53670
$A_{sym}~{}(\mathrm{MeV})$ | -1.56157 | $A_{N}~{}(\mathrm{MeV^{1/2}})$ | 2.93091
$A_{coul}~{}(\mathrm{MeV})$ | 0.060169 | $n$ | -1.06419
$C_{1}~{}(\mathrm{MeV})$ | 13.08641 | $z$ | 6.26334
$C_{2}~{}(\mathrm{MeV^{1/2}})$ | 11.76810 |  |

This expression for the binding energy achieves an RMSE of 4.2 MeV for these
even-even nuclei. Given that the binding energies for nuclei in this range lie
between 500 MeV and 2100 MeV, this suggests that the formula retains an
accuracy of over $99\%$.

When compared to the classical liquid drop model and the Bethe-Weizsäcker
formula [8, 2], the initial three terms in Eq. 6 appear to align with the
volume, symmetry energy, and Coulomb terms, respectively. The term with the
two-thirds power is representative of the surface energy. Any remaining terms
in the expression are likely indicative of additional factors or potential
refinements to our model.

It is also worth noting that an underlying relationship exists between $N$ and
$Z$ due to the stability and other properties of nuclei. Classical
distributions of nuclei often adhere to certain $N/Z$ ratios, as suggested by
the Bethe-Weizsäcker formula. This implies that $N$ and $Z$ are not entirely
independent variables and can, to some extent, be interchanged. Therefore, the
derived analytical expression may not be unique.

##  4  Summary and conclusions

This study has successfully demonstrated the efficacy of Kolmogorov-Arnold
Networks (KANs) in predicting nuclear binding energies, highlighting their
potential in uncovering underlying physical laws from complex data sets. KANs,
based on the Kolmogorov-Arnold representation theorem [22], offer a novel
methodological framework that decomposes multivariate functions into simpler,
univariate functions. This approach not only facilitates the handling of
complex relationships inherent in nuclear physics data but also enhances the
interpretability of the learning models.

Through the implementation of KANs, we have shown that even with minimal input
features, the networks can achieve significant predictive accuracy,
outperforming many traditional models. The addition of more sophisticated
features related to nuclear shell structure and pairing effects further
improved the accuracy, underscoring the importance of these factors in nuclear
physics.

The symbolic regression analysis conducted as part of this study yielded a
simplified, yet powerful, analytical expression for nuclear binding energies.
This expression aligns with known physical models like the classical liquid
drop model and the Bethe-Weizsäcker formula, incorporating key terms for
volume, symmetry energy, Coulomb interactions, and surface energy. The
analytical form derived not only provides a direct insight into the functional
dependencies within the nuclear data but also suggests potential areas for
further refinement of nuclear models.

Overall, the results affirm the robustness of KANs in extracting meaningful
insights from data, characterized by their ability to balance accuracy and
model complexity effectively. The study opens up new avenues for the
application of KANs in other areas of physics and complex systems where
traditional approaches might struggle to cope with the intricacy of the
relationships involved.

Future work will focus on extending the application of KANs to lighter nuclei,
where few-body dynamics become significant, and on refining the network
architecture and training processes to enhance the accuracy and applicability
of the models. Additionally, further exploration into the symbolic expressions
derived from KANs could lead to new theoretical insights and more precise
models in nuclear physics and beyond.

## Acknowledgements

This work has been supported by National Natural Science Foundation of China
(Grants No.12105204 and No.12035011), by the National Key R&D Program of China
(Contracts No. 2023YFA1606503), and by the Fundamental Research Funds for the
Central Universities.

## References

  * [1] H.-W. Hammer, S. König, U. van Kolck, Nuclear effective field theory: Status and perspectives, Rev. Mod. Phys. 92 (2020) 025004\.  doi:10.1103/RevModPhys.92.025004.
URL https://link.aps.org/doi/10.1103/RevModPhys.92.025004

  * [2] C. v. Weizsäcker, Zur theorie der kernmassen, Zeitschrift für Physik 96 (7) (1935) 431–458.
  * [3] M. Rosenbusch, P. Ascher, D. Atanasov, C. Barbieri, D. Beck, K. Blaum, C. Borgmann, M. Breitenfeldt, R. B. Cakirli, A. Cipollone, S. George, F. Herfurth, M. Kowalska, S. Kreim, D. Lunney, V. Manea, P. Navrátil, D. Neidherr, L. Schweikhard, V. Somà, J. Stanja, F. Wienholtz, R. N. Wolf, K. Zuber, Probing the $n=32$ shell closure below the magic proton number $z=20$: Mass measurements of the exotic isotopes ${}^{52,53}\mathrm{K}$, Phys. Rev. Lett. 114 (2015) 202501\.  doi:10.1103/PhysRevLett.114.202501.
URL https://link.aps.org/doi/10.1103/PhysRevLett.114.202501

  * [4] Q. Mo, M. Liu, N. Wang, Systematic study of shell gaps in nuclei, Phys. Rev. C 90 (2014) 024320.  doi:10.1103/PhysRevC.90.024320.
URL https://link.aps.org/doi/10.1103/PhysRevC.90.024320

  * [5] D. Lunney, J. M. Pearson, C. Thibault, Recent trends in the determination of nuclear masses, Rev. Mod. Phys. 75 (2003) 1021–1082.  doi:10.1103/RevModPhys.75.1021.
URL https://link.aps.org/doi/10.1103/RevModPhys.75.1021

  * [6] M. Arnould, S. Goriely, K. Takahashi, The r-process of stellar nucleosynthesis: Astrophysics and nuclear physics achievements and mysteries, Physics Reports 450 (4) (2007) 97–213.  doi:https://doi.org/10.1016/j.physrep.2007.06.002.
URL https://www.sciencedirect.com/science/article/pii/S0370157307002438

  * [7] H. Schatz, W.-J. Ong, Dependence of x-ray burst models on nuclear masses, The Astrophysical Journal 844 (2) (2017) 139.  doi:10.3847/1538-4357/aa7de9.
URL https://dx.doi.org/10.3847/1538-4357/aa7de9

  * [8] H. A. Bethe, R. F. Bacher, Nuclear physics a. stationary states of nuclei, Rev. Mod. Phys. 8 (1936) 82–229.  doi:10.1103/RevModPhys.8.82.
URL https://link.aps.org/doi/10.1103/RevModPhys.8.82

  * [9] W. Myers, W. Swiatecki, Nuclear properties according to the thomas-fermi model, Nuclear Physics A 601 (2) (1996) 141–167.  doi:https://doi.org/10.1016/0375-9474(95)00509-9.
URL https://www.sciencedirect.com/science/article/pii/0375947495005099

  * [10] S. Goriely, S. Hilaire, M. Girod, S. Péru, First gogny-hartree-fock-bogoliubov nuclear mass model, Phys. Rev. Lett. 102 (2009) 242501.  doi:10.1103/PhysRevLett.102.242501.
URL https://link.aps.org/doi/10.1103/PhysRevLett.102.242501

  * [11] M. Hadizadeh, M. Radin, K. Mohseni, A three-dimensional momentum-space calculation of three-body bound state in a relativistic faddeev scheme, Scientific Reports 10 (1) (2020) 1949.  doi:https://doi.org/10.1038/s41598-020-58577-4.
URL https://www.nature.com/articles/s41598-020-58577-4

  * [12] E. Epelbaum, H.-W. Hammer, U.-G. Meißner, Modern theory of nuclear forces, Rev. Mod. Phys. 81 (2009) 1773–1825.  doi:10.1103/RevModPhys.81.1773.
URL https://link.aps.org/doi/10.1103/RevModPhys.81.1773

  * [13] L.-X. Zeng, Y.-Y. Yin, X.-X. Dong, L.-S. Geng, Nuclear binding energies in artificial neural networks, Phys. Rev. C 109 (2024) 034318.  doi:10.1103/PhysRevC.109.034318.
URL https://link.aps.org/doi/10.1103/PhysRevC.109.034318

  * [14] C.-Q. Li, C.-N. Tong, H.-J. Du, L.-G. Pang, Deep learning approach to nuclear masses and $\alpha$-decay half-lives, Phys. Rev. C 105 (2022) 064306.  doi:10.1103/PhysRevC.105.064306.
URL https://link.aps.org/doi/10.1103/PhysRevC.105.064306

  * [15] A. E. Lovell, A. T. Mohan, T. M. Sprouse, M. R. Mumpower, Nuclear masses learned from a probabilistic neural network, Phys. Rev. C 106 (2022) 014305.  doi:10.1103/PhysRevC.106.014305.
URL https://link.aps.org/doi/10.1103/PhysRevC.106.014305

  * [16] Y. Liu, C. Su, J. Liu, P. Danielewicz, C. Xu, Z. Ren, Improved naive bayesian probability classifier in predictions of nuclear mass, Phys. Rev. C 104 (2021) 014315.  doi:10.1103/PhysRevC.104.014315.
URL https://link.aps.org/doi/10.1103/PhysRevC.104.014315

  * [17] L. Neufcourt, Y. Cao, W. Nazarewicz, F. Viens, Bayesian approach to model-based extrapolation of nuclear observables, Phys. Rev. C 98 (2018) 034318\.  doi:10.1103/PhysRevC.98.034318.
URL https://link.aps.org/doi/10.1103/PhysRevC.98.034318

  * [18] Z. M. Niu, H. Z. Liang, Nuclear mass predictions with machine learning reaching the accuracy required by $r$-process studies, Phys. Rev. C 106 (2022) L021303.  doi:10.1103/PhysRevC.106.L021303.
URL https://link.aps.org/doi/10.1103/PhysRevC.106.L021303

  * [19] E. Yüksel, D. Soydaner, H. Bahtiyar, Nuclear mass predictions using machine learning models, Phys. Rev. C 109 (2024) 064322.  doi:10.1103/PhysRevC.109.064322.
URL https://link.aps.org/doi/10.1103/PhysRevC.109.064322

  * [20] Z.-Y. Yuan, D. Bai, Z. Wang, Z.-Z. Ren, Reliable calculations of nuclear binding energies by the gaussian process of machine learning, Nuclear Science and Techniques 35 (6) (2024) 105.
  * [21] Z. Liu, Y. Wang, S. Vaidya, F. Ruehle, J. Halverson, M. Soljačić, T. Y. Hou, M. Tegmark, Kan: Kolmogorov-arnold networks, arXiv preprint arXiv:2404.19756.
  * [22] K. A. N., On the representation of continuous functions of many variables by superposition of continuous functions of one variable and addition, Dokl. Akad. Nauk. SSSR 114 (1957) 953–956.
URL https://cir.nii.ac.jp/crid/1570009749540988928

  * [23] A. B. Givental, B. A. Khesin, J. E. Marsden, A. N. Varchenko, V. A. Vassiliev, O. Y. Viro, V. M. Zakalyukin (Eds.), On the representation of functions of several variables as a superposition of functions of a smaller number of variables, Springer Berlin Heidelberg, Berlin, Heidelberg, 2009, pp. 25–46.  doi:10.1007/978-3-642-01742-1_5.
URL https://doi.org/10.1007/978-3-642-01742-1_5

  * [24] M. Wang, W. Huang, F. Kondev, G. Audi, S. Naimi, The ame 2020 atomic mass evaluation (ii). tables, graphs and references*, Chinese Physics C 45 (3) (2021) 030003.  doi:10.1088/1674-1137/abddaf.
URL https://dx.doi.org/10.1088/1674-1137/abddaf

  * [25] D. M. Allen, The relationship between variable selection and data agumentation and a method for prediction, Technometrics 16 (1) (1974) 125–127.  arXiv:https://www.tandfonline.com/doi/pdf/10.1080/00401706.1974.10489157, doi:10.1080/00401706.1974.10489157.
URL https://www.tandfonline.com/doi/abs/10.1080/00401706.1974.10489157

  * [26] M. Stone, Cross-Validatory Choice and Assessment of Statistical Predictions, Journal of the Royal Statistical Society: Series B (Methodological) 36 (2) (2018) 111–133.  arXiv:https://academic.oup.com/jrsssb/article-pdf/36/2/111/49096683/jrsssb\\_36\\_2\\_111.pdf, doi:10.1111/j.2517-6161.1974.tb00994.x.
URL https://doi.org/10.1111/j.2517-6161.1974.tb00994.x

  * [27] M. Stone, An Asymptotic Equivalence of Choice of Model by Cross-Validation and Akaike’s Criterion, Journal of the Royal Statistical Society: Series B (Methodological) 39 (1) (2018) 44–47.  arXiv:https://academic.oup.com/jrsssb/article-pdf/39/1/44/49117086/jrsssb\\_39\\_1\\_44.pdf, doi:10.1111/j.2517-6161.1977.tb01603.x.
URL https://doi.org/10.1111/j.2517-6161.1977.tb01603.x

  * [28] N. Wang, M. Liu, X. Wu, J. Meng, Surface diffuseness correction in global mass formula, Physics Letters B 734 (2014) 215–219.  doi:https://doi.org/10.1016/j.physletb.2014.05.049.
URL https://www.sciencedirect.com/science/article/pii/S037026931400358X

  * [29] M. R. Mumpower, T. M. Sprouse, A. E. Lovell, A. T. Mohan, Physically interpretable machine learning for nuclear masses, Phys. Rev. C 106 (2022) L021301.  doi:10.1103/PhysRevC.106.L021301.
URL https://link.aps.org/doi/10.1103/PhysRevC.106.L021301


# Supporting Information: q-AQUA: a many-body CCSD(T) water potential,
including 4-body interactions, demonstrates the quantum nature of clusters and
the condensed phase.

Qi Yu  Department of Chemistry Yale University, New Haven, Connecticut 06520,
U.S.A<EMAIL_ADDRESS>Chen Qu  Independent researcher, Toronto, Canada
<EMAIL_ADDRESS>Paul L. Houston<EMAIL_ADDRESS>Department of Chemistry
and Chemical Biology, Cornell University, Ithaca, New York 14853, U.S.A. and
Department of Chemistry and Biochemistry, Georgia Institute of Technology,
Atlanta, Georgia 30332, U.S.A  Riccardo Conte<EMAIL_ADDRESS>Dipartimento di Chimica, Università degli Studi di Milano, via Golgi 19, 20133
Milano, Italy  Apurba Nandi  Department of Chemistry and Cherry L. Emerson
Center for Scientific Computation, Emory University, Atlanta, Georgia 30322,
U.S.A.  Joel M. Bowman<EMAIL_ADDRESS>Department of Chemistry and Cherry
L. Emerson Center for Scientific Computation, Emory University, Atlanta,
Georgia 30322, U.S.A.

## Summary

This supplemental information first contains discussion of differences between
the q-AQUA, MB-pol, and WHBB potentials. Next we provide a description of
three methods to optimize the basis set of polynomials used in calculating
n-body potential energy surfaces: Purification, Pruning, and Analytic
Gradients. Details are then provided for each of the 2-body, 3-body, and
4-body potential energy surfaces. The next section provides further details on
the hexamer results, including comparison of the predictions of the q-AQUA and
MB-pol potentials to the CCSD(T)/CBS calculations for the harmonic frequencies
for four of the water hexamers. This is followed by sections providing details
on the Diffusion Monte Carlo (DMC) calculations and on the Molecular Dynamics
(MD) simulations.

## Differences between q-AQUA, MB-pol, and WHBB potentials

It is important to distinguish the new PES based on Eq. (1) of the main text
from the MB-pol one, and the WHBB one. First, we note that MB-pol is not
strictly a many-body PES in the sense of Eq. (1) of the main text. To see
that, we give the explicit expression for the MB-pol potential for $N$
monomers. It is

$V_{\text{MB-
pol}}(1,\cdots,N)=V_{\text{TTM4-F}}(1,\cdots,N)+\sum_{i>j}^{N}\Delta{V_{2b}(i,j)}S_{2b}+\sum_{i>j>k}^{N}\Delta{V_{3b}}(i,j,k)S_{3b},$
(1)

where $V_{\text{TTM4-F}}(1,\cdots,N)$ is the TTM4-F PES1 and $\Delta
V_{2b}(i,j)$ and $\Delta V_{3b}(i,j,k)$ are correction PESs. These are PIP
fits to the difference of CCSD(T) and TTM4-F energies for 2-b and 3-b
interactions, respectively. These corrections are relatively short range and
they are switched to zero at long range using switching functions, denoted
generically by $S_{2b}$ and $S_{3b}$. (Details of this switching are given in
the MB-pol papers.2, 3) Note that $V_{\text{TTM4-F}}(1,\cdots,N)$ contains the
1-b PS monomer potential4 and also the PS dipole moment surface as part of the
electrostatic interaction. Higher-body interactions are given by polarization
and electrostatic terms within the Thole formalism (some parameters in these
used in MB-pol were optimized by Paesani and co-workers).
$V_{\text{TTM4-F}}(1,\cdots,N)$ itself cannot be written simply as a many-body
expansion. The MB-pol approach is clearly very insightful; however, it does
rely on the accuracy of the TTM4-F potential for 4-b and higher-body
interactions for which no corrections are made.

By contrast the WHBB PES, which preceded the MB-pol one, is given by two
variations. One is simply the MB expansion of Eq. (1) truncated at the 3-b
level, namely

$\begin{split}V_{\text{WHBB}}(1,\cdots,N)=&\sum_{i=1}^{N}V_{1b}(i)+\sum_{i>j}^{N}[V^{\text{CCSD(T)}}_{2b}(i,j)S_{2b}+(1-S_{2b})V^{\text{TTM3-F}}_{2b}(i,j)]+\\\
&\sum_{i>j>k}^{N}V_{3b}(i,j,k)S_{3b},\end{split}$ (2)

where the monomer potential is the PS one, the 2-b term is a PIP fit to
thousands of CCSD(T) energies, which is switched in the range 4.5–5.5 Å to the
2-b interaction from TTM3-F.5, 6 Details of this 2-b PES have been given
previously.7 Finally, a pure ab initio 3-b interaction is a PIP fit to tens of
thousands of MP2 energies. This fit is switched to zero by an exponential
damping function with a range of 8.0 bohr. Again details have been given
previously.5 An option in the WHBB PES is to add 4 and higher-body
interactions using those in TTM3-F. However, this adds considerably to the cpu
effort and makes only marginal improvement in benchmarks comparisons for the
water hexamer.

## Purification, Pruning, and Analytic Gradients

### The Purification Scheme

Consider a group of identical particles whose potential energy surface we wish
to describe using a multi-body expansion (MBE) (see Eq. (1) of the main text).
We wish to describe this potential energy surface by using a basis set of
permutationally invariant polynomials (PIPs) whose coefficients will be
determined by numerical regression so as to smoothly fit a dataset of
electronic energies and, perhaps, gradients for different geometries and whose
polynomials will be powers of, typically, Morse ($\exp(-r_{i,j}/a)$) or
reciprocal variables of the internuclear distances. Permutationally invariant
polynomials have been discussed in detail previously.8, 9, 10 It is convenient
to generate them using the monomial symmetrization approach (MSA) software,11,
12 which ensures that the polynomials have the correct symmetry with respect
to exchange of identical atoms. However, there are two potential issues for
application of MSA to the MBE problem. First, not all of these MSA PIPs will
have the correct limiting behavior for the many body expansion. Second, the
MSA software is not designed to consider the exchange of identical groups of
atoms, so we need to ensure this type of permutational symmetry by another
method.

Consider first the issue of the limiting behavior. We use “purification” to
refer to the process of identifying and setting aside PIPs that do not have
the correct limiting behavior when individual monomers or groups of monomers
are removed from the others to a great distance. Because the fits are to the
$n$-body contributions, when any of the group members are separated by large
distances, thus eliminating the n-body nature of the interaction, the n-body
contribution must go to zero. To ensure the correct limits for all
polynomials, we take an arbitrary geometry, then remove $n$-mers from one
another to large distances, and finally calculate the values of the PIPs to
see whether they go to zero. In the case of the 4-body interactions, we remove
each of the 4 water monomers from the other three and each of the six pairs of
water dimers from one another. In the case of the 3-body interactions, we
remove each water monomer from the remaining pair. In our calculations, the
distances are augmented by 100 Å, and we accept the polynomial as having the
correct behavior if its value is below $10^{-6}$.

Polynomials that do not have the correct limiting behavior cannot immediately
be removed from consideration because there may be other polynomials that, for
example, are composed of products of one with an incorrect limit and one with
a correct limit. Instead, the ones with an incorrect limit are simply given
different names than those with the correct limit. When the process of
evaluating the polynomials is finished, we then look at the definitions of all
those with the correct limits and determine which of the monomials and which
of the renamed polynomials (with incorrect limits) are required to calculate
them. We then remove those that are not required and renumber those that
remain, keeping the order of calculation so that no partial calculation that
contributes to any polynomial needs to be performed twice. The result is a set
of polynomials that all have the correct limiting behavior and that can be
calculated efficiently.

We now turn to the issue of permutational symmetry for exchange of the
identical monomers. As mentioned previously, when the monomers are groups of
atoms, this exchange is not taken into account by the MSA software. Thus, the
polynomials that we create by purification will not, in general, have
permutational symmetry with respect to exchange of identical monomers. One
method for dealing with this issue is to augment the dataset by adding all
relevant permutations of the Cartesian coordinates and assigning them the same
energy. If there are $n$ monomers, this requires a set of $n!$ geometries for
each energy. A better method is to identify those groups of polynomials that
have permutational symmetry with respect to monomer exchange and to then form
“superpolynomials” that are the sum of the polynomial members of each group.
Each group will have $n!$, not necessarily unique, members whose sum is
independent of the monomer permutations and which can then receive a single
coefficient in the energy or gradient determination.

We identify the permutationally invariant groups of polynomials by taking a
single set of $n!$ permutationally related geometries and calculating the
value of each polynomial. Taken together, the groups of polynomials for each
permutation will have the same group of values, but the values of individual
polynomials will vary from permutation to permutation. For each permutation,
one can form pairs of the polynomial identities and their values, and then
sort the pairs by their values. Looking at all pairs that have the same value
component in all permutations gives the identities of the polynomials, some of
which may be repeated, that make up a permutationally invariant group. In
general, there will be as many groups as there were original polynomials.
These groups, each with $n!$ (not necessarily unique) polynomial
contributions, are then summed to form the “superpolynomials”.

### Polynomial Pruning

The MSA software referred to above is generated from just two parameters, the
permutational symmetry of the atoms and the total order desired; i.e.,the
maximum sum of the polynomial powers of the transformed internuclear
distances. The permutational symmetry need not be the complete symmetry of the
system, but must include the relevant exchange possibilities. For example, in
the case of the water trimer, each water would have the symmetry designation
of 21, meaning that the two hydrogens on the water can permute with one
another and the oxygen does not permute. The symmetry of the total system
could be described as 63, but that would imply that the hydrogens or oxygen on
one water could exchange with those on another. At high enough energy, this
rearrangement can certainly happen, but for lower energies, each water stays
intact. The advantage of not using 63 symmetry is that there are more relevant
polynomials and coefficients, so the fit is generally more accurate for the
energy range in which the waters remain intact. Thus, one might consider,
instead, the symmetry 222111, for which the hydrogens on any individual water
can exchange with one another, but atoms on one water do not exchange with
those any other water. Entire water molecules can exchange, but only if this
symmetry is dealt with by one of the methods discussed in the previous
section.

For any given permutational symmetry, the maximum number of PIPs will be
determined by the maximum polynomial order desired. Of course, the processes
of purifying and grouping to achieve symmetric exchange of the waters as a
whole will reduce the numbers of polynomials from that given by the MSA
output. In general and up to a limit, the more polynomials, the more accurate
will be the fit, but the longer will be the time for calculating the energy
and/or gradients. The limit is determined by the need to have far fewer
coefficients than one has energy/gradient constraints so as to avoid
“overfitting”, a situation in which all the points are well-fit but the region
between them has strong oscillations in energy. In many cases, one would like
to have a number of polynomials/coefficients that is somewhere in between the
numbers arrived at for two different choices of order. In such cases, we
“prune” the number of polynomials provided by the larger order to the number
desired for a reasonable compromise between accuracy and speed. We use the
following method to decide which of the polynomials derived from the larger
order should be kept.

A decision on which polynomials are most important may be made on the basis of
the data set. The first step is to evaluate the maximum values of the each
transformed internuclear distance compared to the range of values among the
geometries of the data set. Recall that the internuclear distances are
transformed, usually as Morse or reciprocal values. Thus, long distances have
small transformed values (and are less important) than short distances, which
have transformed values nearer to unity (for the Morse transformation). Taking
the maximal values of the transformed internuclear distances, we then evaluate
all the polynomials available from the MSA output for the larger order. We
then eliminate those polynomials with the smaller values until we arrive at
the desired number of final polynomials. The method works both for regular
polynomials generated by the MSA software and for “superpolynomials” generated
as described in the previous section.

### Analytical Gradients

It is often very desirable to have analytical gradients of the potential so
that forces between atoms can be calculated, either for use in the fitting
process or, for example, in molecular dynamics simulations. Unfortunately,
many higher levels of electronic structure theory do not automatically provide
gradients without substantial execution cost. There are several methods for
accomplishing the goal of providing analytical gradients. One is simply to
differentiate all the terms in the energy expression with respect to each of
the Cartesian coordinates.12 However, this requires a time cost on the order
of $3Nt_{e}$ for calculation of the gradients, where $N$ is the number of
atoms and $t_{e}$ is the time for calculating the energy. A better method is
the reverse (or backward) derivative technique,13, 14 which, for molecules up
to 15 atoms, we have recently shown to have a time cost of 3–4 $t_{e}$,
independent of the number of atoms.15 We use this method for evaluating the
gradients for all the n-body potentials described below.

## Details concerning the Two-Body PES

As pointed out in the literature numerous times and using numerous metrics,
the 2-b interaction accounts for a substantial amount (around 80-90 %) of the
total interaction energy of water clusters. Thus, it is important to obtain
this interaction potential using the highest available level of electronic
structure. CCSD(T) electronic energies for the 2-b interaction have been
reported for flexible monomers previously by us,7 by Paesani and co-workers,2
and more recently by Metz and Szalewicz.16 The data sets span significantly
different energy ranges and configurations. For example our previous dataset
of 30083 energies spans an energy range predominantly from $\pm$ 7 kcal/mol
with a significant number of configurations with energies as high as 28
kcal/mol. The MB-pol data set of about 40000 configurations is predominantly
in the range $\pm$ 7 kcal/mol with a small number of configurations with
energies as high as 14 kcal/mol. Finally the one from ref. 16 consists of 4758
energies and spans the range of -5 to 10 kcal/mol.

We generated a total of 71892 CCSD(T)/CBS 2-b energies with structures
selected from our previous HBB2 data set as well as from the MB-pol data set,
using the criteria that the OO distances be in the range from 2.0 to 8.0 Å.
The 2-b energy is defined as:

$V_{2b}=V_{(\ce{H2O})_{2}}-V_{\ce{H2O},1}-V_{\ce{H2O},2}$ (3)

For each water dimer structure, we first performed CCSD(T)/aug-cc-pVTZ
calculations with an additional 3s3p2d1f basis set, following the same
procedure as MB-pol 2-b, namely, exponents equal to (0.9, 0.3, 0.1) for sp,
(0.6, 0.2) for d, and 0.3 for f. This additional basis is placed at the center
of mass (COM) of each dimer configuration. We also determined the basis set
superposition error (BSSE) correction with the counterpoise scheme. Second, we
performed CCSD(T)/aug-cc-pVQZ calculations with the same additional basis set
and with the BSSE correction. The final CCSD(T)/CBS energies were obtained by
extrapolation over the CCSD(T)/aug-cc-pVTZ and CCSD(T)/aug-cc-pVQZ 2-b
energies. All of these _ab initio_ calculations were performed using Molpro
package.17

Figure 1: *

Figure S1: Distribution of 2-b energies.

Figure S1 shows the distribution of CCSD(T)/CBS 2-b energies in the range of
-10 to 30 kcal/mol. Not shown are additional structures, numbering
approximately 0.5% of the total, that have a high energy beyond this range.
The OO distance distribution for the data set spans the range from 2 to 8 Å.
We fit a basis set of permutationally invariant polynomials to this new 2-b
data set using 7th-order, 42 symmetry. The fitting RMS error for the whole
data set is 25 cm-1. Figure S2 shows the correlation plot between the 2-b fit
and the CCSD(T)/CBS reference energies and confirms the accuracy of the final
fit. However, this 2-b data set only involves dimer structures with OO
distance between 2 and 8 Å, and the long-range 2-body interaction is also
important. To account for the long range 2-body interaction, we applied a
dipole-dipole interaction potential using a high-level dipole moment
surface.18 This long-range 2-b interaction is expressed as:

$V_{2b}^{\text{long-
range}}=\sum_{i=1}^{3}\sum_{j=1}^{3}\frac{q_{i}q_{j}}{r_{ij}}$ (4)

where $q_{i},q_{j}$ are the partial charges on $i$th atom of monomer 1 and
$j$th atom of monomer 2, and $r_{ij}$ is the distance between two atoms. The
partial charge on each atom of the water monomer is obtained from the water
monomer dipole moment surface by Lodi et al.18

Figure 2: *

Figure S2: Correlation plot between 2-b fit and CCSD(T)/CBS reference data.

Finally, a smooth switching function was employed to connect the short-range
2-b fit to the long-range dipole-dipole interaction. The expression for the
final 2-b potential energy surface is:

$V_{2b}=s(r_{\text{OO}})V_{2b}^{\text{fit}}+(1-s(r_{\text{OO}}))V_{2b}^{\text{long-
range}}$ (5)

where $s$ is calculated as

$\displaystyle s$ $\displaystyle=1.0,~{}~{}r_{\text{OO}}\leq 6.5\ \text{\AA}$
$\displaystyle=1-10\left(\frac{r_{\text{OO}}-6.5}{7.8-6.5}\right)^{3}+15\left(\frac{r_{\text{OO}}-6.5}{7.8-6.5}\right)^{4}-6\left(\frac{r_{\text{OO}}-6.5}{7.8-6.5}\right)^{5},~{}~{}6.5<r_{\text{OO}}<7.8\
\text{\AA}$ $\displaystyle=0.0,~{}~{}r_{\text{OO}}\geq 7.8\ \text{\AA}$

Figures S3 and S4 show two potential cuts for the 2-b interaction and compare
q-AQUA, MB-pol, and CCSD(T) results.

Figure 3: *

Figure S3: Comparison of the new 2-b fit and direct CCSD(T) energies for an
attractive cut.

Figure 4: *

Figure S4: Comparison of the new 2-b fit and direct CCSD(T) energies for a
repulsive cut.

## Details concerning the Three-Body PES

We selected a total of 45332 trimer structures from the previous WHBB and MB-
pol trimer configurations, using the criteria of keeping the maximum OO
distance within the 2–9.5 Å range. The 3-b energy is defined as

$V_{3b}=V_{(\ce{H2O})_{3}}-\sum_{i}^{3}V_{(\ce{H2O})_{2},i}+\sum_{i}^{3}V_{\ce{H2O},i}$
(6)

For each water trimer structure, CCSD(T)-F12a/aug-cc-pVTZ calculations were
performed and the BSSE correction was included with the counterpoise scheme.
All calculations were conducted using the Molpro package.17 The distributions
of the 3-b energies and the maximum OO distances are shown in Figure S5. For
the calculated 3-b energies, most trimer structures have 3-b interaction
within the range of -6 to 6 kcal/mol, as can be seen in the left panel of Fig.
S5. For the maximum OO distance distribution, it can be seen from the right
panel in Figure S5 that current data set spans a significantly broader region
than the MB-pol data set.

We divided the new 3-b data set into two separate sets: one with maximum OO
distance in the range from 2.0 to 7.0 Å, and another with maximum OO distance
in the range from 5.0 to 9.5 Å, resulting in two data sets with 42145 and
15282 structures, respectively. The short-range data set was fit using 4th-
order 222111-symmetry permutational invariant polynomials which are functions
of Morse variables, $\exp(-r_{ij}/a)$, where $r_{ij}$ is the internuclear
distance and $a=3$ bohr. The fitting RMS error for this data set is 9 cm-1.
The long-range data set is fit using 3rd-order 222111-symmetry permutational
invariant polynomials, which are functions of inverse of the internuclear
distance, $r_{ij}$. The fitting RMS error for this data set is 11 cm-1. The
correlation plots between two 3-b fits and BSSE corrected CCSD(T)-F12a/aug-cc-
pVTZ reference data are shown in Figs. S6 and S7, respectively.

Figure 5: *

Figure S5: Distribution of 3-b energies, left panel, and OO distances for MB-
pol and q-AQUA data sets.

Our 3-b potential energy surface is a mixed one with the two separate fits
using a switching function based on the maximum OO distance, and also smoothly
switched to zero. The final expression is :

$V_{3b}=s_{1}(r_{\text{OO}})\times\left[s_{2}(r_{\text{OO}})V_{3b}^{\text{short-
range}}+(1-s_{2}(r_{\text{OO}}))V_{3b}^{\text{long-range}}\right]$ (7)

where $s_{1}$ is calculated by

$\displaystyle s_{1}$ $\displaystyle=1.0,~{}~{}r_{\text{OO}}\leq 7.0\
\text{\AA}$
$\displaystyle=1-10\left(\frac{r_{\text{OO}}-7.0}{9.0-7.0}\right)^{3}+15\left(\frac{r_{\text{OO}}-7.0}{9.0-7.0}\right)^{4}-6\left(\frac{r_{\text{OO}}-7.0}{9.0-7.0}\right)^{5},~{}~{}7.0<r_{\text{OO}}<9.0\
\text{\AA}$ $\displaystyle=0.0,~{}~{}r_{\text{OO}}\geq 9.0\ \text{\AA}$

and $s_{2}$ is calculated by

$\displaystyle s_{2}$ $\displaystyle=1.0,~{}~{}r_{\text{OO}}\leq 5.0\
\text{\AA}$
$\displaystyle=1-10\left(\frac{r_{\text{OO}}-5.0}{6.0-5.0}\right)^{3}+15\left(\frac{r_{\text{OO}}-5.0}{6.0-5.0}\right)^{4}-6\left(\frac{r_{\text{OO}}-5.0}{6.0-5.0}\right)^{5},~{}~{}5.0<r_{\text{OO}}<6.0\
\text{\AA}$ $\displaystyle=0.0,~{}~{}r_{\text{OO}}\geq 6.0\ \text{\AA}$

Figure 6: *

Figure S6: Correlation plot between short-range 3-b fit and BSSE-corrected
CCSD(T)-F12a/aVTZ reference data.

Figure 7: *

Figure S7: Correlation plot between long-range 3-b fit and BSSE-corrected
CCSD(T)/CBS reference data.

Figure 8: *

Figure S8: Comparison of the new 3-b fit and direct CCSD(T) energies for an
attractive cut.

Figure S8 shows an attractive cut through the potential as one monomer is
moved relative to the remaining dimer. The results of the q-AQUA and MB-pol
potentials are shown along with the calculated CCSD(T)-F12 results. Both
potentials do well in the range down to about 2.5 Å, while the MB-pol
potential slightly underestimates and the q-AQUA sightly overestimates the
potential at shorter distances. Results for a repulsive cut have been shown in
the main text.

## Details concerning the Four-Body PES

The four-body potential energy surface was a new potential which improved on
the one previously reported.19 It was based on 4th-order 222221111 symmetry.
The basis set was purified to ensure the correct limiting behavior at long
distances between each water monomer and the remaining trimer and between each
pair of water dimers. Compaction reduced the number of monomials, polynomials,
and “renamed” polynomials to 930, 1648, and 279, respectively, from the
original basis set of 2910 monomials and 10,736 polynomials. Determining the
groups of polynomials that were permutationally symmetric with respect to
water exchange further reduced the number of polynomials to 872 while
increasing the number of renamed polynomials to 1927. The advantages of the
grouping are that an unreplicated data set can be used in the fit and that the
fitting time is reduced. In order to achieve a good balance between speed and
accuracy, the number of polynomials was further reduced to 200 by the pruning
method described above. The Fortran program produced for evaluating energies
also used the reverse derivative method for calculating gradients.

The coefficients were determined by fitting the basis set to a set of 3692
geometric points calculated at the CCSD(T)-F12 level. Figure S9 shows the
distributions of the 4-b energy and the OO distances from the data base.
Figure S10 provides a correlation plot between the CCSD(T) energies and the
fit energies. The RMS fitting error is 7.2 cm-1.

Figure 9: *

Figure S9: Distribution of 4-b energies and OO distances.

Figure 10: *

Figure S10: Correlation plot between 4-b fit and CCSD(T)-F12/haTZ references.

One might expect that MB-pol 4-body, which uses the TTM4-F, should be precise
in the long range, but this is not the case, as is shown in Fig. S11. In this
figure, the difference between CCSD(T)-F12 4-body energy and the MB-pol/TTM4-F
4-body energy is plotted against the maximum OO distance in the tetramer, for
all the 3692 configurations in our 4-body data set. We can see that TTM4-F has
large errors even when the OO distance is around 7 Å. The RMSE of the TTM4-F
4-body is 21.2 cm-1; for comparison, our 4-body fit has an RMSE of 7.2 cm-1.

Figure 11: *

Figure S11: Difference between CCSD(T) 4-body energy and TTM4-F 4-body energy,
as a function of the maximum OO distance in the tetramer.

##  1 Details concerning water hexamers

###  1.1 Table of 2-b, 3-b, 4-b, and total dissociation energies for water
hexamer isomers

Table S1 shows results for the water hexamers (a plot of these data is given
in the main text).

Table 1: *

Table S1: 2-b, 3-b, 4-b and total dissociation energies (kcal/mol) for various
water hexamer isomers.  De 2b energy Isomer CCSD(T)/CBSa q-AQUA MB-pol
CCSD(T)/CBSb q-AQUA MB-pol Prism 45.92 46.00 45.73 -38.94 -39.07 -38.93 Cage
45.67 45.67 45.46 -38.47 -38.49 -38.48 Book 1 45.20 44.99 44.59 -36.02 -35.97
-35.81 Book 2 44.90 44.71 44.36 -36.13 -36.04 -35.92 Bag 44.30 44.31 43.71
-35.28 -35.38 -35.24 Ring 44.12 43.46 43.30 -32.71 -32.59 -32.48 Boat 1 43.13
42.71 42.45 -32.30 -32.27 -32.02 Boat 2 43.07 42.66 42.48 -32.24 -32.21 -32.02
MAE / 0.25 0.53 / 0.07 0.15 3-b energy 4-b energy Isomer CCSD(T)/CBSb q-AQUA
MB-pol CCSD(T)/CBSb q-AQUA MB-pol Prism -8.70 -8.71 -8.77 -0.66 -0.65 -0.52
Cage -8.97 -9.08 -8.93 -0.53 -0.49 -0.47 Book 1 -10.38 -10.36 -10.26 -1.08
-1.07 -0.92 Book 2 -10.11 -10.06 -10.03 -1.00 -1.06 -0.85 Bag -10.35 -10.34
-10.15 -1.16 -1.18 -0.90 Ring -11.78 -11.61 -11.60 -1.78 -1.59 -1.44 Boat 1
-11.34 -11.23 -11.30 -1.63 -1.52 -1.35 Boat 2 -11.34 -11.25 -11.29 -1.61 -1.49
-1.35 MAE / 0.07 0.10 / 0.07 0.21

  * a

From Ref. 20

  * b

Ref. 21

###  1.2 Harmonic Frequencies of Water Hexamers

Tables S2 and S3 show a comparison of the harmonic frequencies for four water
hexamers (Prism, Cage, Book1, and Ring) for the q-AQUA and MP-pol potentials
as compared with the benchmark CCSD(T)/CBS calculations of Ref. 22. As can be
seen, the results for q-AQUA and MB-pol are both fairly accurate, with a
slight improvement in the mean absolute error (MAE) from the q-AQUA potential.
For both q-AQUA and MB-pol, the discrepancies between the potential energy
fits and the CCSD(T)/CBS calculations are dominated by the high frequency
modes.

Table 2: *

Table S2: Comparison of harmonic frequencies (in cm-1) of Prism and Cage
hexamers from q-AQUA and MB-pol with CCSD(T)/CBS frequencies.a

Isomer | Prism | Cage
---|---|---
Method | CCSD(T) | q-AQUA | MB-pol | CCSD(T) | q-AQUA | MB-pol
mode 1 | 61 | 61 | 62 | 42 | 42 | 43
mode 2 | 70 | 71 | 73 | 56 | 58 | 57
mode 3 | 74 | 78 | 75 | 74 | 75 | 75
mode 4 | 98 | 91 | 110 | 100 | 99 | 98
mode 5 | 112 | 103 | 120 | 127 | 120 | 129
mode 6 | 149 | 145 | 150 | 153 | 149 | 153
mode 7 | 173 | 168 | 173 | 185 | 173 | 178
mode 8 | 178 | 175 | 179 | 194 | 186 | 191
mode 9 | 212 | 211 | 204 | 210 | 209 | 207
mode 10 | 217 | 225 | 218 | 223 | 218 | 214
mode 11 | 238 | 234 | 230 | 231 | 227 | 223
mode 12 | 246 | 241 | 241 | 234 | 234 | 230
mode 13 | 275 | 279 | 272 | 242 | 244 | 238
mode 14 | 284 | 282 | 277 | 253 | 249 | 247
mode 15 | 287 | 287 | 281 | 281 | 277 | 267
mode 16 | 357 | 346 | 347 | 293 | 295 | 290
mode 17 | 367 | 362 | 366 | 383 | 386 | 378
mode 18 | 420 | 417 | 407 | 395 | 399 | 391
mode 19 | 427 | 422 | 424 | 437 | 437 | 426
mode 20 | 462 | 462 | 447 | 453 | 456 | 441
mode 21 | 491 | 488 | 490 | 467 | 463 | 461
mode 22 | 530 | 530 | 527 | 534 | 528 | 523
mode 23 | 547 | 541 | 536 | 553 | 539 | 542
mode 24 | 612 | 604 | 590 | 620 | 618 | 603
mode 25 | 638 | 636 | 625 | 682 | 683 | 671
mode 26 | 675 | 669 | 661 | 717 | 720 | 700
mode 27 | 711 | 709 | 696 | 774 | 761 | 752
mode 28 | 823 | 821 | 808 | 790 | 785 | 772
mode 29 | 868 | 859 | 851 | 852 | 842 | 833
mode 30 | 1001 | 998 | 977 | 975 | 973 | 952
mode 31 | 1663 | 1658 | 1652 | 1666 | 1660 | 1663
mode 32 | 1674 | 1667 | 1671 | 1673 | 1667 | 1672
mode 33 | 1683 | 1678 | 1680 | 1684 | 1681 | 1681
mode 34 | 1699 | 1696 | 1689 | 1698 | 1694 | 1691
mode 35 | 1716 | 1709 | 1711 | 1707 | 1704 | 1704
mode 36 | 1733 | 1730 | 1729 | 1722 | 1716 | 1712
mode 37 | 3301 | 3317 | 3306 | 3324 | 3339 | 3310
mode 38 | 3509 | 3519 | 3512 | 3517 | 3530 | 3530
mode 39 | 3601 | 3620 | 3603 | 3556 | 3556 | 3563
mode 40 | 3620 | 3651 | 3627 | 3604 | 3632 | 3616
mode 41 | 3717 | 3753 | 3730 | 3650 | 3675 | 3671
mode 42 | 3735 | 3762 | 3748 | 3718 | 3759 | 3736
mode 43 | 3784 | 3812 | 3792 | 3757 | 3784 | 3768
mode 44 | 3799 | 3824 | 3807 | 3792 | 3819 | 3806
mode 45 | 3821 | 3835 | 3819 | 3895 | 3894 | 3901
mode 46 | 3898 | 3898 | 3909 | 3896 | 3905 | 3905
mode 47 | 3899 | 3908 | 3911 | 3899 | 3909 | 3912
mode 48 | 3901 | 3913 | 3912 | 3908 | 3928 | 3916
MAE | / | 7.9 | 7.8 | / | 7.7 | 8.9

  * a

From Ref. 22

Table 3: *

Table S3: Comparison of harmonic frequencies (in cm-1) of book-1 and ring
hexamers from q-AQUA and MB-pol with CCSD(T)/CBS frequencies.a

Isomer | Book-1 | Ring
---|---|---
Method | CCSD(T) | q-AQUA | MB-pol | CCSD(T) | q-AQUA | MB-pol
mode 1 | 27 | 24 | 25 | 28 | 24 | 27
mode 2 | 37 | 37 | 36 | 28 | 24 | 27
mode 3 | 53 | 52 | 51 | 45 | 41 | 43
mode 4 | 67 | 67 | 69 | 45 | 41 | 43
mode 5 | 85 | 86 | 86 | 50 | 45 | 48
mode 6 | 156 | 155 | 155 | 82 | 77 | 77
mode 7 | 179 | 176 | 172 | 156 | 155 | 152
mode 8 | 189 | 183 | 181 | 172 | 161 | 167
mode 9 | 195 | 194 | 189 | 195 | 184 | 182
mode 10 | 225 | 219 | 214 | 195 | 184 | 182
mode 11 | 233 | 229 | 223 | 211 | 205 | 200
mode 12 | 245 | 238 | 232 | 211 | 205 | 200
mode 13 | 250 | 247 | 236 | 254 | 241 | 236
mode 14 | 271 | 269 | 257 | 254 | 241 | 236
mode 15 | 282 | 275 | 266 | 282 | 268 | 263
mode 16 | 291 | 287 | 277 | 292 | 286 | 276
mode 17 | 302 | 301 | 289 | 292 | 286 | 276
mode 18 | 377 | 369 | 367 | 323 | 315 | 305
mode 19 | 393 | 383 | 381 | 407 | 387 | 386
mode 20 | 432 | 422 | 426 | 426 | 404 | 409
mode 21 | 443 | 437 | 437 | 426 | 404 | 409
mode 22 | 467 | 458 | 458 | 441 | 432 | 431
mode 23 | 533 | 525 | 522 | 450 | 434 | 438
mode 24 | 601 | 592 | 586 | 450 | 434 | 438
mode 25 | 708 | 698 | 678 | 757 | 737 | 728
mode 26 | 735 | 731 | 711 | 776 | 753 | 738
mode 27 | 811 | 798 | 777 | 776 | 753 | 738
mode 28 | 829 | 820 | 800 | 867 | 842 | 818
mode 29 | 874 | 866 | 835 | 867 | 842 | 818
mode 30 | 989 | 980 | 953 | 941 | 916 | 892
mode 31 | 1661 | 1654 | 1655 | 1665 | 1665 | 1659
mode 32 | 1673 | 1669 | 1667 | 1676 | 1670 | 1668
mode 33 | 1675 | 1672 | 1670 | 1676 | 1670 | 1668
mode 34 | 1691 | 1682 | 1683 | 1701 | 1691 | 1691
mode 35 | 1702 | 1696 | 1695 | 1701 | 1691 | 1691
mode 36 | 1730 | 1721 | 1721 | 1716 | 1706 | 1705
mode 37 | 3386 | 3397 | 3413 | 3440 | 3469 | 3488
mode 38 | 3455 | 3474 | 3463 | 3505 | 3534 | 3531
mode 39 | 3503 | 3540 | 3521 | 3505 | 3534 | 3531
mode 40 | 3587 | 3598 | 3612 | 3554 | 3595 | 3577
mode 41 | 3637 | 3661 | 3653 | 3554 | 3595 | 3577
mode 42 | 3651 | 3680 | 3667 | 3570 | 3622 | 3596
mode 43 | 3768 | 3791 | 3781 | 3901 | 3901 | 3908
mode 44 | 3893 | 3896 | 3899 | 3901 | 3901 | 3908
mode 45 | 3898 | 3905 | 3908 | 3901 | 3901 | 3908
mode 46 | 3900 | 3907 | 3909 | 3901 | 3902 | 3911
mode 47 | 3900 | 3909 | 3912 | 3901 | 3902 | 3911
mode 48 | 3903 | 3916 | 3913 | 3901 | 3907 | 3912
MAE | / | 8.2 | 12.6 | / | 13.5 | 16.5

  * a

From Ref. 22

##  2 Details of diffusion Monte Carlo calculations

The use of the diffusion Monte Carlo (DMC) method to estimate the vibrational
ground-state energy (i.e., zero-point energy) is based on the similarity
between the diffusion equation and the imaginary-time Schrödinger equation
with an energy shift $E_{\text{ref}}$

$\frac{\partial\psi(\bm{x},\tau)}{\partial\tau}=\sum_{i=1}^{N}\frac{\hbar^{2}}{2m_{i}}\nabla_{i}^{2}\psi(\bm{x},\tau)-\left[V(\bm{x})-E_{\text{ref}}\right]\psi(\bm{x},\tau)$
(8)

The reference energy $E_{\text{ref}}$ in the above equation is used to
stabilize the diffusion system in its ground state and thus is the estimator
of the zero-point energy.23 We employed the unbiased, unconstrained
implementation of DMC, 24 in which the DMC calculation starts from an initial
guess of the ground-state wave function, represented by a population of $N(0)$
equally weighted Gaussian random walkers. These walkers then diffuse randomly
in imaginary time according to a Gaussian distribution and the population is
controlled by birth-death processes, given as

$\displaystyle P_{\text{birth}}=\exp[-(E_{i}-E_{\text{ref}})\Delta\tau]-1{\ \
\ }(\text{if}{\ }E_{i}<E_{\text{ref}}),$ (9) $\displaystyle
P_{\text{birth}}=1-\exp[-(E_{i}-E_{\text{ref}})\Delta\tau]{\ \ \ }(\text{if}{\
}E_{i}>E_{\text{ref}}),$ (10)

where $E_{i}$ is the energy of the $i$th walker. To maintain the number of
random walkers at about the initial value $N(0)$, $E_{\text{ref}}$ is adjusted
at the end of each time step according to

$E_{\text{ref}}(\tau)=\langle V(\tau)\rangle-\alpha\frac{N(\tau)-N(0)}{N(0)}$
(11)

where $N(\tau)$ is the number of walkers at the time step $\tau$ and $\langle
V(\tau)\rangle$ represents the average potential energy of all of the walkers
at that step.

In this study, the imaginary time step $\Delta\tau$ = 10 a.u. and $\alpha$ =
10. The random walkers are initiated at the Prism, Cage, and Book isomers
respectively, and for each isomer, three DMC calculations are performed. In
each DMC calculation, the number of walkers is 50,000, and these walkers are
equilibrated for 10,000 time steps followed by 50,000 propagation steps.
Previous studies show that this number of walkers may produce systematic
uncertainties in terms of the absolute values of the ZPEs due to finite number
of walkers; however, the relative ZPEs between different isomers are reliable,
based on previous studies.25, 26 In addition to the systematic uncertainties,
there are statistical uncertainties because the reference energy oscillates
around the ZPE. The statistical uncertainty is estimated as the standard
deviation of the 3 DMC runs for the same isomer.

## Details of the molecular dynamics simulations

We interfaced the q-AQUA water potential with the i-PI software27 to enable
both classical and path integral molecular dynamics simulations to be
performed for bulk water. Both classical and path integral molecular
simulations were conducted with 256 water molecules in a periodically
replicated simulation box with the experimental density set to be that at 300
K. For the classical MD simulation, at each temperature, we first propagated
the trajectory for 100 ps to reach thermal equilibrium and then ran three
independent trajectories for 100 ps each. The static and dynamical properties
were calculated as an average over the three trajectories. The stochastic
velocity rescaling thermostat is employed for the NVT classical MD
simulations. For the path integral and also the ring polymer MD simulations
(PIMD and RPMD), the trajectories were propagated in the NVT ensemble with a
path-integral white noise Langevin thermostat (PILE-g), and 8 beads were
included. For each temperature, we ran a single trajectory for 50 ps to
achieve thermal equilibrium and then ran three independent trajectories for 50
ps with time step of 0.25 fs.

### Mean square displacement and self-diffusion constants

The self diffusion coefficient, $D$, of liquid water can be calculated from:

$D=\frac{1}{3}\int_{0}^{\infty}\langle\textbf{v}(0)\cdot\textbf{v}(t)\rangle
dt=\frac{1}{6}\lim_{t\rightarrow\infty}\frac{d\langle\parallel\textbf{r}(t)-\textbf{r}(0)\parallel^{2}\rangle}{dt}$
(12)

where $\langle\parallel\textbf{r}(t)-\textbf{r}(0)\parallel^{2}\rangle$ is the
mean square displacement (MSD). For each trajectory, we used the center of
mass of each water molecule to calculate the MSDs and conducted linear fits to
obtain the slope of the MSD curve. The self-diffusion constant D is simply
$1/6$ of the MSD slope and the final values reported in the main text are from
the averaged values over different trajectories.

Figure 12: *

Figure S12: Mean square displacements for classical MD (blue) and PIMD (red)
at the temperature of 298 K as a function of time.

Figure S12 shows two typical MSD curves from classical and ring polymer MD
simulation at 298 K. As seen, both of them show linear behavior with the
propagation of time and it is reliable to perform a linear fit for obtaining
the diffusion coefficient. However, it should be noted that to obtain more
accurate and converged diffusion constants of liquid water, longer time
trajectories are required. We will conduct more extensive MD simulations in
the future.

### Second-order rotational time correlation functions and orientational
relaxation time

The second-order rotational time correlation function for OH bond is
calculated as:

$C_{2}(t)=\big{\langle}P_{2}(\textbf{u}_{\text{OH}}(t)\cdot\textbf{u}_{\text{OH}}(0))\big{\rangle}$
(13)

where $\textbf{u}_{\text{OH}}(t)$ is the unit vector along each OH bond, and
$P_{2}$ is the second-order Legendre polynomial.

Figure 13: *

Figure S13: Typical second-order rotational time correlation functions,
$C_{2}(t)$, from classical and ring polymer MD simulations at 298 K

Figure S13 shows two typical second-order rotational time correlation
functions at 298 K from classical MD and RPMD simulations. To predict the
orientational relaxation time, $\tau_{2}$, for both classical MD and RPMD
simulations, we fit the correlation function $C_{2}(t)$ to a single
exponential function first. The fitted exponential function was used to
describe the long-time decay of $C_{2}(t)$ for times longer than 10 ps. We
then numerically integrated the $C_{2}(t)$ function covering both short-time
and long-time regions to obtain $\tau_{2}$.

### Many body effects in liquid water

As mentioned above, the q-AQUA water potential is constructed from 1-b,
intrinsic 2-b, 3-b and 4-b interactions. It becomes natural to investigate how
different components of potential influence the properties of liquid water. In
Figure S14, we show both the OO radial distribution function and the MSD
curves for liquid water at 298 K from classical MD simulation with different
levels of potentials. As seen, using only 1-b + 2-b, both the static
structural properties and dynamic properties of liquid water deviate
considerably from the experiment. In particular, the peak positions of the OO
radial distribution function do not agree with the experimental results. The
slope of MSD curve is around 6 times of that when using full q-AQUA potential,
where the latter is proven to have a diffusion coefficient close to the
experiment. Thus, the use of only the 1-b + 2-b interactions cannot describe
the important properties of liquid water – higher body interactions are
needed. This is verified by the fact that when the 3-b interaction is added
(see blue lines in Figure S14), both the radial distribution function and the
MSD curve reach reasonable agreement with experiment. The self diffusion
coefficients calculated from the 1-b + 2-b + 3-b interactions with classical
MD simulations at different temperatures are listed in Table S4. It can be
seen that when the PES is truncated at three-body level, it can already
provide qualitatively reasonable diffusion coefficients comparing with the
experiment.

Finally, when the full q-AQUA potential is employed with the inclusion of 4-b
interactions, Figure S14 shows that classical MD simulations give a more
localized OO radial distribution function with higher distribution around 4-5
Å. The slope of the MSD is lowered compared with that calculated with the 1-b
+ 2-b + 3-b interactions, and the diffusion constant at 298 K is also further
decreased. These observations call into question of the role of the 4-b
interaction and the adequacy of a classical MD simulation for accurately
capturing all the static and dynamic properties.

Figure 14: *

Figure S14: Effects of 2-b, 3-b, and 4-b interactions on the static structural
and dynamical properties of liquid water at 298 K from classical MD
simulations. Left: the OO radial distribution function (Experimental data are
from Ref. 28), Right: the mean square displacement (MSD) as a function of
time. MD (1-2 body) is the MD simulation using only 1-b + 2-b interactions. MD
(1-3 body) is the MD simulation using 1-b + 2-b + 3-b interactions. MD (1-4
body) is the MD simulation using full q-AQUA potential with all 1-b, 2-b, 3-b
and 4-b interactions included

### Nuclear quantum effects in liquid water

Figure S15 shows the effect of using 1- through 3-b interactions (left panel)
or 1- through 4-b interactions (right panel) and using either classical
molecular dynamics (MD) or quantum path-integral molecular dynamics (PIMD).
For both the 1- through 3-body and 1- through 4-body calculations, the quantum
PIMD calculation is in substantially better agreement with experiment than
that with classical MD.

Figure 15: *

Figure S15: Nuclear quantum effects on the OO radial distribution function for
liquid water at 298 K from classical MD and PIMD simulations. Left: using
potential with 1-b + 2-b + 3-b interactions. Right: using full q-AQUA
potential with all 1-b, 2-b, 3-b and 4-b interactions included. Experimental
data are from Ref. 28

Table 4: *

Table S4: Self-diffusion constant (Å2/ps) of liquid water at different
temperatures using different components of the potential and different methods

Use full q-AQUA potential (1-b + 2-b + 3-b + 4-b)
---
Temperature (K) | Classical | RPMD | Expt.a
278 | 0.080 $\pm$ 0.016 | 0.130 $\pm$ 0.015 | 0.131
288 | 0.102 $\pm$ 0.017 | 0.177 $\pm$ 0.015 | 0.177
298 | 0.145 $\pm$ 0.012 | 0.226 $\pm$ 0.020 | 0.230
320 | 0.248 $\pm$ 0.011 | 0.331 $\pm$ 0.016 | 0.360
Use partial q-AQUA potential (1-b + 2-b + 3-b)
Temperature (K) | Classical | RPMD | Expt.a
278 | 0.088 $\pm$ 0.006 | 0.170 $\pm$ 0.014 | 0.131
288 | 0.131 $\pm$ 0.011 | 0.216 $\pm$ 0.010 | 0.177
298 | 0.188 $\pm$ 0.012 | 0.242 $\pm$ 0.012 | 0.230
320 | 0.332 $\pm$ 0.031 | 0.344 $\pm$ 0.020 | 0.360

  * a

from Ref. 29 and Ref. 30

### Computational cost of q-AQUA water potential

Table S5 shows the computation cost of the q-AQUA model for both energy and
gradient calculations for a 256-water system, using single core and multicores
of 2.4 GHz Intel Xeon processor. Two calculations are shown, one with periodic
boundary conditions (PBC) and one without. The number of n-body configurations
considered in the calculation is shown in the second column and is limited by
a cut-off for the maximum OO distance. For the 2-b interaction, in gas-phase
cluster calculations there is no cut-off because the potential is switched in
the long range to the dipole-dipole interaction. In the MD simulation, a 14 Å
cut-off was used. Figs. S3 and S4 show that the potential is essentially zero
at this distance. In the 3-b calculation, the cut-off limit was 9 Å, at which
distance the potential is also exceedingly small, as seen from Fig S8. For the
4-b interaction, the switching range to zero for the MD simulations was set s
5.0–6.0 Å in order to limit the number of calculations needed. A longer range
of 5.8–7.5 Å was used in for the gas-phase cluster work. In the MD
simulations, the same cut-off maximum distances were used in calculations both
with and without PBCs. Reverse derivative methods were used in all
calculations with the exception of long-range 2-b ones that are based on the
dipole-dipole interactions.

Table 5: *

Table S5: The computation cost of the q-AQUA method for energy and gradient
calculations of a 256 water system  No PBC Time for energy (s) Time for
energy+gradient (s) Component Number 1 core 8 core 1 core 8 core 1-b 256 0.002
0.002a 0.003 0.003a 2-b 32640 0.23 0.02 0.72 0.08 3-b 84051 0.42 0.05 1.94
0.26 4-b 115922 1.26 0.17 4.34 0.52 Total 2.00 0.35 7.12 0.97 With PBC Time
for energy (s) Time for energy+gradient (s) Component Number 1 core 8 core 1
core 8 core 1-b 256 0.002 0.002a 0.003 0.003a 2-b 48919 0.38 0.04 1.18 0.18
3-b 195058 0.84 0.12 5.30 0.70 4-b 268304 2.91 0.37 11.55 1.31 Total 4.52 0.83
18.05 2.49

  * a

1-b terms are not parallelized

## References

  * Burnham et al. 2008 Burnham, C. J.; Anick, D. J.; Mankoo, P. K.; Reiter, G. F. _J. Chem. Phys._ 2008, _128_ , 154519
  * Babin et al. 2013 Babin, V.; Leforestier, C.; Paesani, F. _J. Chem. Theory Comput._ 2013, _9_ , 5395
  * Babin et al. 2014 Babin, V.; Medders, G. R.; Paesani, F. _J. Chem. Theory Comput._ 2014, _10_ , 1599
  * Partridge and Schwenke 1997 Partridge, H.; Schwenke, D. W. _J. Chem. Phys._ 1997, _106_ , 4618
  * Wang et al. 2009 Wang, Y. M.; Shepler, B. C.; Braams, B. J.; Bowman, J. M. _J. Chem. Phys._ 2009, _131_ , 054511
  * Fanourgakis and Xantheas 2008 Fanourgakis, G. S.; Xantheas, S. S. _J. Chem. Phys._ 2008, _128_ , 074506
  * Shank et al. 2009 Shank, A.; Wang, Y.; Kaledin, A.; Braams, B. J.; Bowman, J. M. _J. Chem. Phys._ 2009, _130_ , 144314
  * Braams and Bowman 2009 Braams, B. J.; Bowman, J. M. _Int. Rev. Phys. Chem._ 2009, _28_ , 577–606
  * Bowman et al. 2010 Bowman, J. M.; Braams, B. J.; Carter, S.; Chen, C.; Czaḱo, G.; Fu, B.; Huang, X.; Kamarchik, E.; Sharma, A. R.; Shepler, B. C.; Wang, Y.; Xie, Z. _J. Phys. Chem. Lett._ 2010, _1_ , 1866–1874
  * Qu et al. 2018 Qu, C.; Yu, Q.; Bowman, J. M. _Annu. Rev. Phys. Chem._ 2018, _69_ , 6.1–6.25
  * Xie 2019 Original MSA Software. 2019; https://www.mcs.anl.gov/research/projects/msa/, Accessed: 2019-12-20
  * msa 2019 MSA Software with Gradients. https://github.com/szquchen/MSA-2.0, 2019; Accessed: 2019-01-20
  * Baydin and Pearlmutter 2014 Baydin, A. G.; Pearlmutter, B. A. _JMLR: Workshop and Conference Proceedings, ICML 2014 AutoML Workshop_ 2014, 1–7
  * Baydin et al. 2017 Baydin, A. G.; Pearlmutter, B. A.; Radul, A. A.; Siskind, J. M. _The Journal of Machine Learning Research_ 2017, _18_ , 5595–5637
  * Houston et al. 2022 Houston, P. L.; Qu, C.; Nandi, A.; Conte, R.; Yu, Q.; Bowman, J. M. _J. Chem. Phys._ 2022, _156_ , 044120
  * Metz and Szalewicz 2020 Metz, M. P.; Szalewicz, K. _J. Chem. Theory Comput._ 2020, _16_ , 2317–2339
  * Werner et al. 2015 Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M. MOLPRO, version 2015.1, a package of ab initio programs. 2015; see http://www.molpro.net
  * Lodi et al. 2011 Lodi, L.; Tennyson, J.; Polyansky, O. L. _J. Chem. Phys._ 2011, _135_ , 034113
  * Nandi et al. 2021 Nandi, A.; Qu, C.; Houston, P. L.; Conte, R.; Yu, Q.; Bowman, J. M. _J. Phys. Chem. Lett._ 2021, _12_ , 10318–10324
  * Bates and Tschumper 2009 Bates, D. M.; Tschumper, G. S. _J. Phys. Chem. A_ 2009, _113_ , 3555–3559
  * Reddy et al. 2016 Reddy, S. K.; Straight, S. C.; Bajaj, P.; Huy Pham, C.; Riera, M.; Moberg, D. R.; Morales, M. A.; Knight, C.; Götz, A. W.; Paesani, F. _J. Chem. Phys._ 2016, _145_ , 194504
  * Howard and Tschumper 2015 Howard, J. C.; Tschumper, G. S. _J. Chem. Theory Comput._ 2015, _11_ , 2126
  * Anderson 1975 Anderson, J. B. _J. Chem. Phys._ 1975, _63_ , 1499–1503
  * Kosztin et al. 1996 Kosztin, I.; Faber, B.; Schulten, K. _Am. J. Phys._ 1996, _64_ , 633–644
  * Mallory et al. 2015 Mallory, J. D.; Brown, S. E.; Mandelshtam, V. A. _J. Phys. Chem. A_ 2015, _119_ , 6504–6515
  * Wang et al. 2012 Wang, Y. M.; Babin, V.; Bowman, J. M.; Paesani, F. _J. Am. Chem. Soc._ 2012, _134_ , 11116
  * Kapil et al. 2019 Kapil, V. et al. _Comput. Phys. Commun._ 2019, _236_ , 214–223
  * Skinner et al. 2013 Skinner, L. B.; Huang, C.; Schlesinger, D.; Pettersson, L. G. M.; Nilsson, A.; Benmore, C. J. _J. Chem. Phys._ 2013, _138_ , 074506
  * Mills 1973 Mills, R. _J. Phys. Chem._ 1973, _77_ , 685–688
  * Holz et al. 2000 Holz, M.; Heil, S. R.; Sacco, A. _Phys. Chem. Chem. Phys._ 2000, _2_ , 4740–4742