# Topology and complexity of the hydrogen bond network in classical models of
water

Fausto Martelli  IBM Research Europe, Hartree Centre, Daresbury, WA4 4AD,
United Kingdom<EMAIL_ADDRESS>Department of Physics and CNR
Institute of Complex Systems, Sapienza University of Rome, P.le Aldo Moro 5,
00185 Roma, Italy<EMAIL_ADDRESS>
###### Abstract

Over the years, plenty of classical interaction potentials for water have been
developed and tested against structural, dynamical and thermodynamic
properties. On the other hands, it has been recently observed (F. Martelli et.
al, ACS Nano, 14, 8616–8623, 2020) that the topology of the hydrogen bond
network (HBN) is a very sensitive measure that should be considered when
developing new interaction potentials. Here we report a thorough comparison of
11 popular non polarizable classical water models against their HBN, which is
at the root of water properties. We probe the topology of the HBN using the
ring statistics and we evaluate the quality of the network inspecting the
percentage of broken and intact HBs. For each water model, we assess the
tendency to develop hexagonal rings (that promote crystallization at low
temperatures) and pentagonal rings (known to frustrate against crystallization
at low temperatures). We then introduce the _network complexity index_ , a
general descriptor to quantify how much the topology of a given network
deviates from that of the ground state, namely of hexagonal or cubic ice.
Remarkably, we find that the network complexity index allows us to relate, for
the first time, the dynamical properties of different water models with their
underlying topology of the HBN. Our study provides a benchmark against which
the performances of new models should be tested against, and introduces a
general way to quantify the complexity of a network which can be transferred
to other materials and that links the topology of the HBN with dynamical
properties. Finally, our study introduces a new perspective that can help in
rationalizing the transformations among the different phases of water and of
other materials.

Water models, Classical potentials, Network topology, Network complexity,
Hydrogen bonds

††preprint: AIP/123-QED

##  I Introduction

On our planet, water is the only substance that can be found co-existing in
the solid, liquid and vapour phases outside research laboratories. Its
molecular simplicity hides a remarkably wide list of anomalous behaviors that
stretch over the the most complex phase diagram of any pure substance Salzmann
(2019), and whose origin lies in a critical point located at low temperatures
and low pressures Palmer  _et al._ (2014); Sellberg  _et al._ (2014);
Debenedetti, Sciortino, and Zerze (2020); Kringle  _et al._ (2020); Kim  _et
al._ (2020). Nonetheless, water plays a primary role in many industrial,
biological and geological processes. Therefore, there is a great interest in
developing classical interaction potentials able to embrace water’s complex
nature. This intent is aggravated by the wide span of thermodynamic conditions
at which water exists (from low temperatures and low pressures of the
interstellar medium to high pressures and high temperatures in the core of
planets), and by the broad range of timescales required for several processes
to occur (from heterogeneous nucleation to protein folding to geological
processes).

After a long-lasting debate Poole  _et al._ (1992); Liu  _et al._ (2010);
Limmer and Chandler (2011); Wikfeldt, Nilsson, and Pettersson (2011); Palmer,
Car, and Debenedetti (2013); Limmer and Chandler (2013); Palmer  _et al._
(2014); Limmer and Chandler (2014); Chandler (2016); Palmer  _et al._ (2016a,
2018, b), it is now accepted that liquid water is a mixture Palmer  _et al._
(2014); Sellberg  _et al._ (2014); Palmer  _et al._ (2018); Debenedetti,
Sciortino, and Zerze (2020); Kim  _et al._ (2020); Kringle  _et al._ (2020) of
molecules whose local neighborhood constantly change between an ordered
tetrahedral state with local lower density, and a more distorted tetrahedral
state with local higher density Martelli  _et al._ (2020); Shi and Tanaka
(2020a, b); Shi, Russo, and Tanaka (2018a); Russo and Tanaka (2014a); Akahane
and Tanaka (2018); Santra  _et al._ (2015); Huang  _et al._ (2009); Nilsson
and Pettersson (2015); Wikfeldt, Nilsson, and Pettersson (2011); De Marzio
_et al._ (2017); Martelli (2019). These two local structures continuously
interchange with each other, giving rise to a complex network of bonds that
actively determines the properties of water at the macro scale Martelli
(2019). The percentage of such local environments depends on the thermodynamic
conditions Martelli (2019); Wikfeldt, Nilsson, and Pettersson (2011), and
eventually liquid water becomes a 1:1 mixture at low temperatures and
pressures Palmer  _et al._ (2014); Debenedetti, Sciortino, and Zerze (2020).

Classical interaction potentials have been developed over the years to
reproduce experimental structural, dynamical and/or thermodynamic properties
of water (over a relatively reduced set of thermodynamic points), and several
studies have compared them against a plethora of observables Palmer  _et al._
(2016b); Pekka and Lennart (2001); Mao and Zhang (2012); Lee and Kim (2019);
Jorgensen  _et al._ (1983); Harrach and Drossel (2014); González  _et al._
(2010); Zielkiewicz (2005); Steinczinger, Jóvári, and Pusztai (2017); Dix,
Lue, and Carbone (2018). On the other hand, even though Bernal and Fowler
first recognized almost 100 years ago that water molecules build a complex
network of bonds that is at the heart of water’s anomalous behavior Bernal and
Fowler (1933), the pivotal role of the hydrogen bond network (HBN) has so far
been showcased only in a handful of important cases involving transformations
between complex phases of water Tse  _et al._ (1999); Marton̆ák, Donadio, and
Parrinello (2004, 2005); Palmer  _et al._ (2014); Shephard  _et al._ (2017);
Martelli  _et al._ (2018); Martelli (2019), the mutual interactions between
water and biological membranes Martelli, Crain, and Franzese (2020), and
between water and graphene sheets for technological purposes Chiricotto  _et
al._ . In particular, the inspection of the topology of the HBN in biological
environments is opening new avenues in enhancing the efficacy of new
drugs/vaccines Martelli, Calero, and Franzese (2021). The HBN has never been
used as a metric to test water models. Two reasons for this lack of comparison
are: (i) it is very hard to have a direct experimental description of the HBN
in liquid phases, and (ii) it is hard to probe the HBN from a computational
perspective.

Here we fill the gap of the point (ii) above by comparing the HBN of 11
popular and widely adopted classical water models at ambient conditions. We
study the TIP3P Jorgensen  _et al._ (1983) the SPC Berendsen  _et al._ (1981),
the SPC/E Berendsen, Grigera, and Straatsma (1987) and the flexible SPC Toukan
and Rahman (1985); Amira, Spángberg, and Hermansson (2004) as 3-points models;
the TIP4P Jorgensen  _et al._ (1983), the TIP4P-Ice Abascal  _et al._ (2005),
the TIP4P/2005 Abascal and Vega (2005), the flexible TIP4P/2005 González and
Abascala (2011) and the TIP4P-Ew Horn  _et al._ (2004) as 4-points models; the
TIP5P Mahoney and Jorgensen (2000) and the TIP5P-Ew Rick (2004) as 5-points
models. We probe the topology of the HBN using the ring statistics, a
theoretical tool that has seen increasingly high relevance in determining the
properties of bulk water Martelli  _et al._ (2016, 2018); Martelli (2019);
Formanek and Martelli (2020); Palmer  _et al._ (2014); Santra  _et al._
(2015); Leoni  _et al._ (2019); Camisasca  _et al._ (2019); Marton̆ák,
Donadio, and Parrinello (2004, 2005); Russo and Tanaka (2014b); Palmer  _et
al._ (2014); Fitzner  _et al._ (2019); Shi and Tanaka (2018), of aqueous
solutions Bakó  _et al._ (2017); Pothoczki, Pusztai, and Bakó (2018, 2019); Li
_et al._ (2020a) and of water under confinement Martelli, Crain, and Franzese
(2020); Chiricotto  _et al._ . We then measure the quality of the HBN for each
water model in terms of broken and intact HBs, a measure intimately linked to
the fluidity and tetrahedrality of water DiStasio Jr.  _et al._ (2014);
Martelli, Crain, and Franzese (2020).
We also introduce the concept of _network complexity index_ $\xi$ that allows
to quantify how much a given HBN deviates from the HBN of water in the
crystalline phase at low temperatures and ambient pressure, i.e., hexagonal or
cubic ice. As a showcase, we apply this index $\xi$ to one of the rings
definition and counting scheme here adopted. Remarkably, the index $\xi$
allows us to link dynamical properties of water with the topology of the
underlying HBN and the corresponding structural properties. Finally, we show
that the HBN topology is more sensitive to the size of the simulation box with
respect to, e.g., structural properties measured _via_ the two-bodies pair
correlation function. Therefore, the inspections of the network topology
should always be considered when simulating network-forming materials. This
issue becomes particularly relevant when dealing with _ab initio_ molecular
dynamics simulations which are restricted to small simulation cells.

The article is organized as follows. In Section II we report the details of
the numerical simulations, the ring definitions and counting schemes, and we
introduce the network complexity index. In Section III we report our main
findings for all 11 water models here inspected. Conclusions and final remarks
are reported in Section IV.

##  II Computational details

In this section we describe the numerical setup, the protocols implemented to
count rings and we introduce the definition of the network complexity index.

###  II.1 Numerical simulations

Our study is based on classical molecular dynamics simulations of systems
composed of $N=1100$ water molecules described by different 11 interaction
potentials in the isobaric ($N$$p$$T$) ensemble. We have employed Nosé-Hoover
thermostat Nosé (1984); Hoover (1985) with 0.2 ps relaxation time to maintain
constant temperature at $T=300$ K, and Parrinello-Rahman barostat Parrinello
and Rahman (1981) with 2 ps relaxation time to maintain constant pressure at 1
bar. We have performed simulations with the GROMACS 18.0.1 package Abraham
_et al._ (2015). All simulations have been equilibrated for 1 ns. The
production runs achieved 3 ns. For each water model we have averaged over 10
independent trajectories. Our analysis investigates the properties of liquid
water in the presence of thermal noise.

###  II.2 Ring statistics and network complexity index

In order to compute the ring statistics it is necessary to follow two steps.
First of all, it is necessary to define the link between atoms/molecules.
Possible definitions can be based on the formation of bonds, interaction
energies, geometric distances, etc.. The second step is the definition of ring
and the corresponding counting scheme. This task is of particular relevance in
directional networks, like water or silica, were the donor/acceptor nature of
the bonds breaks the symmetry in the linker search path. Several definitions
of rings and counting schemes have been reported in the literature King
(1967); Rahman and Stillinger (1973); Guttman (1990); Franzblau (1991); Wooten
(2002); Yuan and Cormack (2002); Roux and Jund (2010). Such different
definitions have yielded to different interpretations of numerical results
even for the simplest crystalline structures, not to mention more complicated
networks such as amorphous silicate structures Jin  _et al._ (1994); Hobbs
_et al._ (1998); Guttman (1990); Marians and Hobbs (1988, 1990); Yuan and
Cormack (2002) and water Fitzner  _et al._ (2019); Camisasca  _et al._ (2019);
Martelli (2019); Santra  _et al._ (2015). In the case of water, these
inconsistencies have been recently reconciled by Formanek and Martelli showing
that they were caused by different counting schemes Formanek and Martelli
(2020).
Here, we report a thorough analysis of the HBN for the 11 classical non-
polarizable models of water. This benchmark study adopts the three ring
definitions and counting schemes reported in Ref. Formanek and Martelli
(2020). The definition of HB follows the geometric construction described in
Ref. Luzar and Chandler (1996). In this regard, any quantitative measure of
HBs in liquid water is somewhat ambiguous, since the notion of an HB itself is
not uniquely defined. However, qualitative agreement between the definition
here adopted and several other proposed definitions have been deemed
satisfactory over a wide range of thermodynamic conditions, including the one
here investigated Prada-Gracia, Shevchuk, and Rao (2013); Shi, Russo, and
Tanaka (2018b). We construct rings by starting from a tagged water molecule
(marked as 1 in fig. 1) and recursively traversing the HBN until the starting
point is reached again or the path exceeds the maximal ring size considered
(12 water molecules in our case). In the first scheme sketched in fig. 1 a),
molecule 1 _donates_ a HB emphasized by the red arrow, and we restrict our
counting only to the shortest rings King (1967), i.e., to rings that can not
be further decomposed into smaller ones (the 6-folded ring in this case). We
will refer to this scheme to as d1. This scheme emphasizes the directional
nature of the HBs resulting in an enhanced hexagonal character of the HBN
Formanek and Martelli (2020). It is therefore well suited for investigating
phenomena such as nucleation, where hexagonal rings are the elemental building
blocks of the final, crystal network. In the second scheme sketched in fig. 1
b), we consider only rings formed when molecule 1 _accepts_ a HB and we
restrict our counting only to the shortest ring King (1967) (the 5-folded ring
in this case). We will refer to this scheme to as d2. This scheme emphasizes
the formation of pentagonal rings in the network Formanek and Martelli (2020),
which are known to play an important role at supercooled conditions
frustrating against crystallization Russo and Tanaka (2014a); Shi, Russo, and
Tanaka (2018a); Martelli (2019), as well as in promoting the crystallization
of clathrate structures Li  _et al._ (2020b). These two counting schemes d1
and d2 provide drastically different distributions Formanek and Martelli
(2020) because of the different energies involved in accepting and donating a
HB. As we will show in the forthcoming discussion, such difference translates
into distinct percentages of coordination defects of the kind
$\textit{A}_{2}\textit{D}_{1}$ and $\textit{A}_{1}\textit{D}_{2}$, where
$\textit{A}_{x}\textit{D}_{y}$ indicates that a water molecules accepts $x$
and donates $y$ HBs. In the third scheme sketched in fig. 1 c) we loosen the
previous restrictions. Molecule 1 can now either accept or donate a HB, hence
counting both the 6-folded and the 5-folded ring. In the following, we will
refer to this counting scheme to as d3.
As shown in Ref. Formanek and Martelli (2020), these three different counting
schemes carry different, but complementary physical information and,
therefore, allow us to make a proper comparison on the topology of the HBN
generated by different classical models of water.

Figure 1:  Schematic representation of the three counting schemes. Red filled
circles represent oxygen atoms, while red empty circles represent hydrogen
atoms. The network of HBs is represented by the arrows Water molecule labeled
as 1 is the starting molecule from which rings are counted. a): Molecule
donates a HB (red arrow) and the shortest ring is the hexagonal one. b):
Molecule 1 can accept a HB (one of the two red arrows) and the shortest ring
is the pentagonal one c): Molecule 1 can either donate one HB or accept a HB,
generating both the short hexagonal and pentagonal ring.

When studying the topology of a network via the ring statistics, the
probability distribution P(n) of having a n-folded ring is a normalized
quantity that does not reflect the overall number of each n-folded ring. On
the other hand, the actual number of rings is important to understand the
degree of complexity in a network (for a given number of molecules). The
stable crystalline phase of water at ambient pressure is the hexagonal(cubic)
ice, Ih(c), characterized by an hexagonal network of bonds. In a sample of
liquid water at higher temperatures, the thermal noise allows water molecules
to explore a larger configurational space and, hence, the underlying HBN hosts
also shorter and longer rings. Intuitively, such network is more ”complex”
with respect to the HBN of the ground state, and its fluctuations are related
to the dynamical properties (translational and rotational diffusion) of water
molecules. Knowing the symmetry of the HBN at the ground state, we can measure
the deviation from it. We introduce the _network complexity index_ $\xi$
defined as the ratio between the number of 6-folded rings and the total number
of rings:

$\xi=\frac{n_{6}}{\sum_{i=3}^{12}n_{i}}$ (1)

where $n_{i}$ is the number of the $i-$th ring. The sum on the denominator of
eq. 1 runs from n=3, the shortest ring length possible, to n=12, the longest
ring length here considered. The network complexity index $\xi$ encodes the
number of rings in a given network. For the ground state Ih(c), the network
complexity index is trivially $\xi=1$. We can define the case of maximal
disorder in the network when each ring length occur with the same frequency in
the network. This case corresponds to a flat distribution in P(n) and, since
in our case we consider 10 possible ring lengths (from n=3 to n=12),
$\xi=0.1$. Here, we compute $\xi$ for the counting scheme d3, scheme c) in
fig. 1 which, compared to the other schemes, accounts for the presence in the
network of longer rings.

The definition of network complexity index can be extended to other networks
with different ground states (i.e., non-hexagonal networks), in which case the
nominator in eq. 1 will be $n_{j}$, $j$ being the length of the characteristic
ring length at the ground state. However, when comparing the index $\xi$ in
different simulations, care must be taken in order to ensure that (i) the same
definition of ring and counting schemes are used, (ii) the same maximum search
path (or maximal ring size) is implemented.

It is worthy to remark, at this point, that our analysis occurs in the
presence of thermal noise. As recently shown by Montes de Oca et al., the
percentage of broken and intact bonds drastically changes upon removal of the
thermal noise de Oca  _et al._ (2020), i.e., in correspondence with the
inherent potential energy surfaces (IPES). While we do expect quantitative
differences in correspondence with the IPES with respect to the results that
we will show in Section III, we are confident that the overall trend should be
preserved.

##  III Results

In this section, we report and discuss our results in terms of topology and
complexity of the HBNs, as well as their quality in terms of broken and intact
HBs for the water 11 models. We then compare the network complexity indices
and we relate them to the dynamical properties of each water model. Finally,
we show and discuss the effects of the simulation box size on the topology of
the HBN.

###  III.1 3 points models

we start our analysis by comparing, in fig. 2, the oxygen-oxygen two-bodies
pair correlation functions g2(r) for 3-points models with the g2(r) obtained
from various scattering experiments Skinner  _et al._ (2013); Soper and
Benmore (2008) (open circles and squares) and _ab initio_ molecular dynamics
(AIMD) simulations DiStasio Jr.  _et al._ (2014) (open triangles) at the PBE0
level of theory accounting for non-local van der Waals/dispersion interactions
DiStasio Jr.  _et al._ (2014). It is worth to mention, at this point, that the
PBE0+vdw level of theory gives an accurate g2(r) compared to the experimental
one, but does not capture the correct density difference between liquid water
at ambient conditions and Ih.
We can observe that the TIP3P model (black line) has the first peak at $\sim
2.8$ nm with intensity comparable with that of scattering experiments and _ab
initio_ simulations. On the other hand, after the first minimum the
distribution is almost flat with no peaks at larger distances. This is
indicative of a lack of structurization at larger distances. The SPC model
(red line) is able to perform better than the TIP3 model, with a g2(r) showing
hints of a second and a third peak, though shifted compared to the
experimental and _ab initio_ g2(r). On the other hand, the intensity of the
first peak overcomes the experimental and AIMD first peak. The SPC/E model
(green line) results in a better g2(r) in terms of intensity and position of
the second and third peak with respect to the experimental and the _ab initio_
g2(r), reflecting the better performances in density and diffusion constant
than the SPC model Berendsen, Grigera, and Straatsma (1987). On the other
hand, the intensity of the first peak is further enhanced. The flexible SPC
water model (blue line) is a re-parametrization of the SPC water model in
which the O–H stretching is made anharmonic, and thus the dynamical behavior
is well described and bulk density and permettivity are correct Praprotnik,
Janežič, and Mavri (2004). The g2(r) of the flexible SPC model is
characterized by a high intensity first peak and a deeper first minimum with
respect to the rigid SPC model, and correctly captures the profile of the
g2(r) at larger distances.

Overall, the sequence TIP3, SPC, SPC/E and flexible SPC is characterized by an
increment in the height of the first peak that also shifts slightly towards
lower distances, a corresponding deepening of the first minimum and a gradual
appearance of a second and third peak that tend to overlap with the
experimental and with AIMD ones for the SPC/E and the flexible SPC models.

Figure 2:  The oxygen-oxygen two-bodies pair correlation, g2(r), of liquid
water for the TIP3P (black), SPC (red), SPC/E (green) and the flexible SPC
model (blue). The g2(r) obtained from various scattering experiments Skinner
_et al._ (2013); Soper and Benmore (2008) and _ab initio_ molecular dynamics
simulations DiStasio Jr.  _et al._ (2014) are reported for comparison with
open symbols.

In fig. 3 we report the probability distribution P(n) of having a n-folded
ring, with n$\in[3,12]$ for the classical 3 points models. The distribution
for the TIP3P model is reported as open circles, SPC as open squares, SPC/E as
open diamonds and flexible SPC as open triangles. It is worth to remark, at
this point, that the P(n) is a normalized distribution and, therefore, it does
not reflect the actual number of each ring.
Panel a) reports the P(n) according to the definition d1 sketched in fig. 1
a). According to this counting scheme that emphasizes the directionality of
the HBs, all networks have a dominating hexagonal character. Such character is
milder in the TIP3P model (black open circles), which shows a very broad P(n)
and whose network accommodates also 10- and 11-folded rings. The hexagonal
character of the network grows moving to the SPC model (red open squares),
with a corresponding reduction of longer (n$>$8) rings. The reduction of
longer rings develops along with an enhancement of 5- and 7-folded rings,
which are comparable to 6-folded rings in terms of energy Camisasca  _et al._
(2019). In particular, the HBN of the SPC model described with this counting
scheme is almost completely deprived of rings with n$>$10\. The hexagonal
character of the network further grows in the SPC/E (green open diamonds) and
in the flexible SPC (open blue triangles). These models show almost identical
P(n)s. In particular, besides the enhanced hexagonal character, the HBNs of
these models have an enhanced pentagonal and heptagonal character.
Consequently, the contribution from longer rings, namely rings with n=8 and
n=9, is less marked. Therefore, loosening the holonomic constraints and
allowing the (re-parametrized) SPC model to vibrate has a major effect not
just on the structural properties (as shown from the g2(r), fig. 2) but also
on the topology of the HBN.
Panel b) reports the P(n) according to the definition d2 sketched in fig. 1
b). As shown in Ref. Formanek and Martelli (2020), this counting scheme
emphasizes the pentagonal character of the HBN, as the starting water molecule
must accept one HB instead of donating. Interestingly, the TIP3P model shows,
according to this counting scheme, an almost equal pentagonal and hexagonal
character with a tail populating configurations up to n=10. The HBN of the SPC
model, on the other hand, shows an improved pentagonal character with a slight
increase in the hexagonal character and a reduction of longer rings which
mostly disappear at n$>$9\. The pentagonal character of the HBN over the
hexagonal one become particularly dominant in the SPC/E and in the flexible
SPC models, for which rings characterized by n$>$8 are mostly absent.
Panel c) reports the P(n) according to the definition d3 sketched in fig. 1
c). Since this counting scheme does not discern among the donor/acceptor
character of the HBs and does not implement the shortest path criterion King
(1967), the resulting topology is more complex with respect to the previous
counting schemes. Therefore, for this scheme we also compute the network
topology index $\xi$ (eq. 1) and report them in table 1. According to this
counting scheme, the HBN for all the inspected models show a similar character
in terms of hexagonal and heptagonal rings. The TIP3P model shows a quite
broad distribution with a considerable contribution of longer (n$>$7) rings.
The presence of rings with n=12 reflects the more complex topology of the
network with respect to the previous counting schemes. Interestingly, such
distribution is comparable with that of the TIP3P model optimized MacKerell
Jr.  _et al._ (1998) for biological simulations Martelli, Crain, and Franzese
(2020). Because of the very broad character of P(n), the TIP3P model has a low
value of network complexity index $\xi=0.1496$, not too far from the value
$\xi=0.1$ that characterizes (as described in Section II.2) a flat
distribution P(n) with equal populations. Moving to the SPC model, we observe
a slight enhancement of n=5, and a more pronounced enhancement of n=6 and n=7
which, as for the TIP3P model, equally characterize the HBN. Contrarily, the
network is deprived of longer rings, namely of n=10, n=11 and n=12 rings. The
slight increase in n=6 causes a small increment in the network complexity
index to $\xi=0.1665$. Moving to the SPC/E model, we observe a mild increment
in n=5, with a relevant increment in n=6 and n=7, as well as a marked
depletion of longer rings. The complexity index for the SPC/E model further
increases to $\xi=0.1863$. A similar distribution characterizes the HBN of the
flexible SPC model, which almost overlaps with the P(n) of the SPC/E model.
The distribution of the flexible SPC model suggests that, as mentioned above,
the introduction of flexibility plays a very important role in shaping the
network. The complexity index for the flexible SPC model rises to
$\xi=0.1922$.

Figure 3:  Probability distributions of the hydrogen-bonded n-folded rings, P(n), for liquid water at ambient conditions described by the TIP3P (black open circles), the SPC (red open squares), the SPC/E (green open diamonds) and the flexible SPC (blue open triangles) models. Panel a) reports the P(n) according to the definition sketched in fig. 1 a); panel b) reports the p(n) according to the definition sketched in fig. 1 b); panel c) reports the P(n) according to the definition sketched in fig. 1 c). | TIP3P | SPC | SPC/E | SPC-Flex
---|---|---|---|---
$\xi$ | 0.1496 | 0.1665 | 0.1863 | 0.1922
Table 1: Values of the network complexity index $\xi$ computed for the rings
counting scheme d3 for the 3 points models of water.

In fig. 4, we report the percentage of broken and intact HBs for the 3 point
models of water and, for comparison, the _ab initio_ liquid water at the PBE0
level of theory with vdW long range interactions at a temperature of 330K to
account for nuclear quantum effects DiStasio Jr.  _et al._ (2014). We adopt
the following syntax: $\textit{A}_{x}\textit{D}_{y}$ indicates the number of
acceptors ($\textit{A}_{x}$) and donors ($\textit{D}_{y}$) HBs. The network in
_ab initio_ liquid water (black open circles) is dominated by a intact HBs
($\textit{A}_{2}\textit{D}_{2}$), which account for $\sim 48\%$. The second
highest configuration is the $\textit{A}_{1}\textit{D}_{2}$, with $\sim 20\%$,
followed by $\textit{A}_{2}\textit{D}_{1}$ ($12\%$),
$\textit{A}_{2}\textit{D}_{2}$ ($\sim 10\%$) and
$\textit{A}_{3}\textit{D}_{2}$ ($\sim 5\%$). The percentage of this last
configuration, i.e., the configuration $\textit{A}_{3}\textit{D}_{2}$, is
shared with all 3 points models. The percentage of broken and intact HBs in
the TIP3P model (red open squares) quantitatively differs from the
distribution of _ab initio_ liquid water, but qualitatively shows a similar
behaviour. In particular, the distribution for the TIP3P model is dominated by
a markedly reduced amount of intact HBs with a percentage
$\textit{A}_{2}\textit{D}_{2}\sim 33\%$, followed by
$\textit{A}_{1}\textit{D}_{2}$ with $\sim 23\%$,
$\textit{A}_{2}\textit{D}_{1}$ with $\sim 12\%$ and an almost equal percentage
of $\textit{A}_{1}\textit{D}_{1}$. The low percentage of intact HBs explains
the very broad distribution of rings (fig. 3) as well as the absence of
structurization beyond the first hydration shell in the g2(r) (fig. 2). Moving
from the TIP3 model to the SPC model (open green diamonds) we observe a slight
decrease in the configuration $\textit{A}_{1}\textit{D}_{1}$ to $\sim 10\%$
and an increment in the percentage of intact HBs to $\sim 37\%$. Therefore,
the SPC model is characterized by a slightly enhanced percentage of fully
coordinated configurations which cause an enhancement in 5-, 6-, and 7-folded
rings, as reported in fig. 3, and are responsible for the appearance of the
second peak in the g2(r) (fig. 2). The modified SPC models, namely the SPC/E
(open blue triangles) and the flexible SPC model (open orange left triangles),
are both characterized by a distribution that mostly overlap with the _ab
initio_ liquid water potential. The high percentage of intact HBs reaches
values in the order of $\sim 45-48\%$ for both SPC/E and the flexible SPC
models, and this enhanced tetra-coordinated character explains the enhanced
number of 5-, 6-, and 7-folded rings and the corresponding decrease in longer
rings, as reported in fig. 3, as well as an enhancement in the intensity of
the first peak in the g2(r) and a corresponding deepening of the following
minimum.

Figure 4:  Percentage-wise decomposition of the intact HBs per water molecule
into acceptor-(A) and donor-(D) for _ab initio_ liquid water at T=330 K as
black open circles, and for the 3 points models. The TIP3P model is reported
as red open squares, the SPC as green open diamonds, the SPC/E as blue open
triangles and the flexible SPC as orange open left triangles. The $x$-axis
labels $\textit{A}_{x}\textit{D}_{y}$ indicate the number of acceptor
($\textit{A}_{x}$) and donor ($\textit{D}_{y}$) HBs. For clarity we omit
combinations with minor contributions, e.g., $\textit{A}_{3}\textit{D}_{1}$,
$\textit{A}_{0}\textit{D}_{y}$, $\textit{A}_{x}\textit{D}_{0}$, etc.

###  III.2 4 points models

In fig. 5 we report the g2(r) for the 4 points models, namely the TIP4 (black
line), the TIP4P-Ew (red line), the TIP4P/2005 (green line), the flexible
TIP4P/2005 (blue line) and the TIP4P-ice (orange line) models. As for the 3
points models, the open symbols represent the g2(r) from experimental data
Skinner  _et al._ (2013); Soper and Benmore (2008) and from _ab initio_
molecular dynamics simulations DiStasio Jr.  _et al._ (2014), and serve as a
comparison. The TIP4P model is characterized by the intensity of the first
peak as high as $\sim 3$, the closest to the experimental/_ab initio_ ones.
The depth of the following minimum is slightly more pronounced than the
experimental/AIMD. At larger distances, on the other hand, the g2(r) almost
overlap with the experimental/AIMD ones. The TIP4P-Ew (red line) has been
optimized for simulating water in biological environments Horn  _et al._
(2004). With respect to the TIP4P model, the TIP4P-Ew shows a higher intensity
in the first peak of the g2(r), reaching $\sim 3.2$. Other features of the
g2(r) are qualitatively overlapping with the experimental/AIMD distribution
functions. A similar trend characterizes the g2(r) for the TIP4P/2005 (green
line) and for the flexible TIP4P/2005 (blue line) which are mostly
indistinguishable from the distribution function of the TIP4P-Ew function. The
TIP4P-Ice, on the other hand, shows a g2(r) more structured with respect to
the previous ones. The intensity of the first peak reaches values as close as
$\sim 3.5$, and the following minimum shows a deeper depth. Finally, the
intensity of the second peak is slightly more pronounced with respect to the
experimental/AIMD ones. Such ”over”-structurization is not surprising
considering that this water model has been developed to study crystalline and
non crystalline solid forms of ice Abascal  _et al._ (2005).

Overall, moving from the TIP4P model to the TIP4P-Ew, TIP4P/2005, flexible
TIP4P/2005 and the TIP4P-Ice we observe a systematic enhancement of the first
peak which also slightly shifts towards larger values. All models well capture
the second and third peaks.

Figure 5:  The oxygen-oxygen two-bodies pair correlation, g2(r), of liquid
water for the TIP4P (black), TIP4P-Ew (red), TIP4P/2005 (green), the flexible
TIP4P/2005 (blue), and the TIP4P-Ice (orange) models. The g2(r) obtained from
various scattering experiments Skinner  _et al._ (2013); Soper and Benmore
(2008) and _ab initio_ molecular dynamics simulations DiStasio Jr.  _et al._
(2014) are reported for comparison with open symbols.

In fig. 6 we report the probability distribution P(n) of having a n-folded
ring, with n$\in[3,12]$ for the TIP4P (black open circles), the TIP4P-Ew (red
open squares), the TIP4P/2005 (green open diamonds), the flexible TIP4P/2005
(blue open triangles) and the TIP4P-Ice (orange open left triangles) models.
In panel a) we report the P(n) according to the definition d1 sketched in fig.
1 a). All classical models of water show a P(n) maximized at n=6. In
particular, the TIP4P model shows, among all the 4 points models here
inspected, the broader distribution, with an almost equal contribution of
5-folded and 7-folded rings. Interestingly, the distribution for the TIP4P
model is less broad compared with the distribution for the TIP3P model
reported in fig. 3 a), with almost no rings longer than n=9. The HBN of the
TIP4P-Ew model shows a slight increase in the pentagonal and in the hexagonal
character and a corresponding lower character in rings longer than n=7. Very
similar distributions occur for the TIP4P/2005 and for the flexible TIP4P/2005
models. The hexagonal character of the HBN is further emphasized in the
TIP4P-Ice model for which we can observe also a slight increment in the
pentagonal character and a depletion of rings longer than n=7. This enhanced
hexagonal character of the HBN is a consequence of the parametrization of this
model which has been optimized to reproduce crystalline and non-crystalline
solid forms of ice.
In panel b) of fig. 6 we report the P(n) computed using the rings and counting
scheme definition d2 sketched in fig. 1 b). According to such counting scheme,
all distributions are maximized over n=5. In particular, the P(n) for the
TIP4P model shows the broader distribution, with some contribution from n=8,
while the contributions to the HBN from longer rings are mostly negligible.
The TIP4P-Ew model shows an enhanced pentagonal character and a slightly
increment of the hexagonal character as well, with a reduction of the
contribution coming from rings longer than n=6. Such tendency becomes more
pronounced moving to the TIP4P/2005 and to the flexible TIP4P/2005 models,
while differences among these three models are minimal. On the other hand, the
pentagonal character of the HBN is particularly enhanced in the TIP4P-Ice
model whose HBN, according to this counting scheme, dose not account for rings
longer than n=7.
In panel c) we report the P(n) computed using the counting scheme d3 sketched
in fig. 1 c). For this scheme, we also compute the complexity indices $\xi$,
reported in table 2. We can observe that the distribution P(n) for the TIP4P
model is fairly broad and mostly dominated by an almost equal amount of 6- and
7-folded rings followed by 5-folded rings, while the amount of longer rings
decreases with increasing the rings lengths. The complexity index for the
TIP4P model is $\xi=0.1721$, which is larger compared to the TIP3 and the SPC
models, but lower compared to the other 3 points models here inspected. The
network of TIP4P-Ew model is characterized by a consistent increment in the
hexagonal character, followed by the heptagonal character and a smaller
increment in the pentagonal character. Longer rings weight less with respect
to the TIP4P model. The complexity index for the TIP4P-Ew is $\xi=0.1887$,
larger compared to the TIP4P model. The distribution P(n) for the TIP4P/2005
and the flexible TIP4P/2005 models mostly overlap with the P(n) for the
TIP4P-Ew, and have a similar complexity index, namely $\xi=0.1915$ for the
TIP4P/2005 model and $\xi=0.1984$ for the flexible TIP4P/2005 model.
Interestingly, the P(n) for the TIP4P-Ew, for the TIP4P/2005 and for the
flexible TIP4P/2005 models show a breaking in the equal hexagonal and
heptagonal character of the network, with a slight predominance of the
hexagonal character which justifies the higher values of the index $\xi$. The
network of the TIP4P-Ice model is characterized by a further enhancement of
the hexagonal character, followed by the heptagonal and hexagonal characters,
with a depletion of longer rings. Such enhancement in the hexagonal character
is reflected by the complexity index $\xi=0.2155$, the closest to the value of
ice among all models inspected in this work.

The similar distributions of rings for the TIP4P-Ew, TIP4P/2005 and flexible
TIP4P/2005 in all counting schemes suggest that such models perform equally in
terms of HBN, and this reflects the very similar g2(r) (fig. 5). On the other
hand, the TIP4P-Ice, which shows a more structured g2(r) (fig. 5), tends to
favour a network with a shorter connectivity.

Figure 6:  Probability distributions of the hydrogen-bonded n-folded rings, P(n), for liquid water at ambient conditions described by the TIP4P (black open circles), the TIP4P-EW (red open squares), the TIP4P/2005 (green open diamonds), the flexible TIP4P/2005 (blue open triangles), and the TIP4P-Ice (open orange triangles) models. | TIP4P | TIP4P-Ew | TIP4P/2005 | TIP4P/2005-Flex | TIP4P-Ice
---|---|---|---|---|---
$\xi$ | 0.1721 | 0.1887 | 0.1915 | 0.1984 | 0.2155
Table 2: Values of the network complexity index $\xi$ computed for the rings
counting scheme d3 for the 4 points models of water.

In fig. 7 we report the percentage of broken and intact HBs for the 4 points
classical models and for _ab initio_ liquid water. We can observe that the
percentage of broken and intact HBs for all 4 points models qualitatively
resembles that of _ab initio_ liquid water (open black circles). With respect
to AIMD water, the TIP4P model (open red squares) underestimates the
percentage of $\textit{A}_{2}\textit{D}_{2}$ configurations to $\sim 43\%$ and
a slightly higher percentage of $\textit{A}_{1}\textit{D}_{2}$ configurations
($\sim 22\%$). All other configurations mostly overlap with the configurations
of AIMD water. The percentage of broken and intact HBs for the TIP4P-Ew (open
green diamonds), for the TIP4P/2005 (open blue triangles) and the flexible
TIP4P/2005 (open orange left triangles) mostly overlap with each other,
reflecting the almost overlapping distributions of rings and the similar
values of the complexity index $\xi$ (fig. 6). In particular, we observe
almost no differences between the distribution for the TIP4P/2005 and the
flexible TIP4P/2005 models. With respect to the TIP4P model, they recover the
percentage of $\textit{A}_{2}\textit{D}_{2}$ for the AIMD liquid water, while
showing a small reduction of $\textit{A}_{1}\textit{D}_{1}$ defects ($\sim
7\%$). Such increment in the percentage of $\textit{A}_{2}\textit{D}_{2}$ to
$\sim 48\%$ reflects the stronger 5-, 6- and 7-folded character of the HBN
described in fig. 6. The TIP4P-Ice model, on the other hand, overestimate the
percentage of intact HBs with respect to AIMD water, reaching $\sim 55\%$ of
the total configurations. This enhanced 4-folded coordination explains the
distribution of rings (fig.6) describing an HBN particularly enriched in 6-
and 7-folded rings, as well as the higher intensity in the first peak of the
g2(r) (fig 5).

Overall, we can state that, overall, the 4 points models are characterized by
similar HBNs, with the exception of the TIP4P model which shows a broad
distribution of rings caused by a lower percentage of intact HBs, and the
TIP4P-Ice model, whose high percentage of intact HBs causes generates a
network characterized by a strong hexagonal and heptagonal character with a
marked reduction of longer rings.

Figure 7:  Percentage-wise decomposition of the intact HBs per water molecule
into acceptor-(A) and donor-(D) for _ab initio_ liquid water at T=330 K as
black open circles, and for the 4 points models. The TIP4P model is reported
as red open squares, the TIP4P-Ew model as green open diamonds, the TIP4P/2005
model as blue open triangles, the flexible TIP4P/2005 as orange open left
triangles, and the TIP4P-Ice model as brown open lower triangles. The
distribution for the TIP4P/2005 model almost perfectly overlaps with the
distribution for the flexible TIP4P/2005 model.

###  III.3 5 points models

In fig. 8 we compare the g2(r) of two 5 points models, namely the TIP5P (black
line) and the TIP5P-E (red line) models with the g2(r) obtained from various
scattering experiments Skinner  _et al._ (2013); Soper and Benmore (2008) and
_ab initio_ molecular dynamics simulations DiStasio Jr.  _et al._ (2014)
reported as open symbols. We can observe that the g2(r) for both 5 points
models mostly overlap. With respect to the 3 points and to the 4 points
models, the g2(r) of both 5 points models are more closer to the experimental
and to the AIMD g2(r). The intensity of the first peak for the TIP5P and the
TIP5P-E models is below 3, namely $\sim$2.9. The depth of the first minimum is
slightly more pronounced with respect to the experimental and the AIMD g2(r),
while experimental and AIMD peaks at longer distances are well captured.

Figure 8:  The oxygen-oxygen two-bodies pair correlation, g2(r), of liquid
water for the TIP5P (black) and the TIP5P-EW (red) models. The g2(r) obtained
from various scattering experiments Skinner  _et al._ (2013); Soper and
Benmore (2008) and _ab initio_ molecular dynamics simulations DiStasio Jr.
_et al._ (2014) are reported for comparison with open symbols.

In fig. 9 we report the probability distribution P(n) of having a n-folded
ring, with n$\in[3,12]$ for the TIP5P model (open black circles) and for the
TIP5P-E model (open red squares). In panel a) we report the P(n) computed
using the counting scheme d1 sketched in fig. 1 a) and which emphasizes the
directionality of the HBN. Both the TIP5P and the TIP5P-E models provide
similar distributions which are slightly maximized at n=6, and contributions
of rings up to n=11.
In panel b) we report the P(n) computed using, as a counting scheme, the
definition d2 sketched in fig. 1 b). As for the previous case, also in this
case the two distributions are mostly indistinguishable and with an almost
equal pentagonal and hexagonal character. Overall, the pentagonal character of
both HBNs is less pronounced with respect to the four points models, and
qualitatively resemble the P(n) for the three points model TIP3P (fig. 3 b)).
In panel c) we report the P(n) computed using the definition d3 and sketched
in fig. 1 c). For this counting scheme, we also compute the complexity indices
reported in table 3. As for the previous counting schemes, the differences
between the two P(n)’s are minimal, indicating that the two model provide
similar HBNs. Such similarity can be quantified observing that both models
have almost the same complexity index $\xi$, i.e., $\xi=0.1774$ for the TIP5P
model and $\xi=0.1754$ for the TIP4P-Ew model. The values of $\xi$ for both
networks are comparable with that of the TIP4P model but, with respect to the
TIP4P model, the hexagonal character of the HBN is more pronounced for both
the TIP5P and the TIP5P-E models. Besides the TIP4P-Ice model, among all other
interaction potentials here inspected, the TIP5P and the TIP5P-E models are
the only one for which the counting scheme d3 is (slightly) maximized towards
n=6.

Figure 9:  Probability distributions of the hydrogen-bonded n-folded rings, P(n), for liquid water at ambient conditions described by the TIP5P (black open circles) and the TIP5P-EW (red open squares). | TIP5P | TIP5P-E
---|---|---
$\xi$ | 0.1774 | 0.1754
Table 3: Values of the network complexity index $\xi$ computed for the rings
counting scheme d3 for the 5 points models of water.

In fig. 10 we report the percentage of broken and intact HBs for the 5 points
classical models and for _ab initio_ liquid water. With respect to the AIMD
liquid water (open black circles), both the TIP5P (open red squares) and the
TIP5P-E (open green diamonds) models are characterized by a markedly lower
percentage of intact HBs ($\sim 36\%$), while low-coordinated defects occur in
higher percentages. The $\textit{A}_{1}\textit{D}_{2}$ configuration account
for the $\sim 22\%$ of the total configurations, followed by
$\textit{A}_{2}\textit{D}_{1}$ configurations with a percentage of $\sim 16\%$
and $\textit{A}_{1}\textit{D}_{1}$ configurations with $\sim 12\%$ of the
total configurations. The low amount of intact HBs (compared with the AIMD
water), reflects the broad distribution of rings and the corresponding low
value of the network complexity indices $\xi$ (fig. 9).

Figure 10:  Percentage-wise decomposition of the intact HBs per water molecule
into acceptor-(A) and donor-(D) for _ab initio_ liquid water at T=330 K as
black open circles, and for the 5 points models. The TIP5P model is reported
as red open squares and the TIP5P-EW model as green open diamonds.

###  III.4 Relation between network complexity and dynamical properties

We here show how the network complexity is linked to dynamical properties,
namely the translational and rotational diffusion. This link comes from the
observation that bonding can be viewed as a competition between the energy
gained from the formation of a bond, and the entropy loss due to the reduction
in configurational volume that occurs when two particles are constrained to
stay close relative to each other. The establishment of an extended network of
bonds occurs when the energy gain (that controls the lifetime of bonds)
balances the entropy loss.
For each water model we have computed the diffusion coefficient and the
rotational relaxation time ($\tau_{rot}$), and we have reported them against
the network complexity index $\xi$. We have computed the diffusion coefficient
from the mean squared displacement, and the rotational relaxation time
$\tau_{rot}$ from the integral of the rotational autocorrelation function as
reported in Refs Calero, Stanley, and Franzese (2016); Martelli, Calero, and
Franzese (2021) $C_{rot}(t)=\left<\textit{OH}(t)\cdot\textit{OH}(0)\right>$,
i.e., $\tau_{rot}=\int_{0}^{+\infty}C_{rot}(t)dt$. Fig. 11 show the values of
these three observables in a three dimensional plot, with projections on the
corresponding two dimensional spaces. The values for 3-points models are
reported in red, the values for 4-points models are reported in green, and the
values of 5-points models are reported in blue. We can observe a clear
correlation between the complexity of the HBN and the dynamical properties of
water molecules. The models with the highest diffusion coefficients and the
fastest rotational relaxation times are the TIP3P and the SPC models, which
are also characterized by the lowest values of the index $\xi$. Contrarily,
the TIP4P-Ice model is the model with the lowest diffusion coefficient and the
slowest rotational relaxation time, and the highest index $\xi$. Overall, we
can observe that for all models of water the higher the value of $\xi$
(reported in ascending order in the tables 1, 2, 3), the slower the diffusion
coefficient and the rotational relaxation time. Therefore, we can assert that
there is a clear correlation between the complexity of the HBN and dynamical
properties, a relation never observed before. Such correlation suggests that
faster diffusion and rotations allow water molecules to increase the possible
connections between each other, hence increasing the configurational space
that the network can explore resulting in a more complex topology able to host
a larger amount of longer rings.

Figure 11:  Three dimensional plot reporting the projection on the
corresponding two dimensions of the values acquired by the 11 classical models
of water. We report the values for the network complexity index $\xi$, the
diffusion coefficient and the rotational relaxation time $\tau_{rot}$. Data
for 3-points models are reported in red, while data for the 4-points models
are reported in green and data for the 5-points models are reported in blue.

###  III.5 System size dependence

We now turn our attention to the study of finite size effects, commonly
inspected when computing physical quantities to check whether a system suffers
from periodicity artifacts. In fig. 12 a) we report the oxygen-oxygen g2(r)
for a simulation box containing 500 water molecules (black continuous line),
1000 water molecules (red dashed line) and 1500 water molecules (green dotted-
dashed line) interacting with the TIP4P classical interaction potential. We
can observe that the three g2(r) perfectly overlap. In fig. 12 b) we report
the probability distribution P(n) of having a n-folded ring for the three
cases inspected above and computed according to the ring definition and
counting scheme d3. We can observe that the distribution computed for the
smaller simulation box with N=500 molecules is remarkably different from the
distributions with N=1000 and N=1500 molecules. In particular, we observe a
strong enhancement of n=12 rings which causes a reduced contribution of
shorter rings to the P(n). This result indicates that a search path of n=12
water molecules is too long for a small simulation box with only N=500 water
molecules, and the increment in n=12 is caused by periodicity in the
simulation box. It is worth to mention that the definitions d1 and d2 do not
show such behavior (data not reported), as the network investigated with these
definitions does not host rings as long as n=12 (see fig. 6 upper and middle
panels). Therefore, although the g2(r) for N=500 is the same as the g2(r) for
larger simulation boxes, care must be taken when inspecting the network
topology and in how such inspection is performed.

Figure 12:  Panel a): system size dependence on the oxygen-oxygen two bodies
pair correlation function for the TIP4P model using a box containing 500 water
molecules (black continuous line), 1000 water molecules (red dashed line) and
1500 water molecules (green dotted-dashed line). Panel b): ring distribution
for a system described by the TIP4P model in a box containing 500 water
molecules (black pluses), 1000 water molecules (red crosses) and 1500 water
molecules (green stars).

##  IV Conclusions

In this article we have tested 11 popular non polarizable classical
interaction potentials for water against their hydrogen bond networks (HBNs).
We have probed the topology of the HBN using three schemes that emphasize
different physical features. We have evaluated the quality of the HBNs in
terms of broken and intact HBs, and we have linked our results to structural
properties measured _via_ the two bodies pair correlation function g2(r). We
have then introduced the network complexity index $\xi$ that measures how much
the topology of a HBN deviates from that of the ground state, and we have
tested it to one of the three rings counting schemes. We have shown that the
index $\xi$ is directly related to dynamical properties, hence establishing a
clear cause-effect relationship between molecular motions and network
connectivity. Finally, we have inspected how periodicity artifacts can
influence the topology of the HBN. We have performed all studies at ambient
conditions, i.e., T=300 K and p=1 bar. Although different water models have
(very) different melting points, their network topology –and hence their
network complexity– remain roughly unchanged away from the limit of
supercooling Formanek and Martelli (2020). On the other hand, when water is
under confinement, water molecules in the proximity of the surfaces undergo a
drastic change in the dynamics Samatas  _et al._ (2018); Martelli, Crain, and
Franzese (2020); Chiricotto  _et al._ ; Gallo, Rovere, and Chen (2010);
Camisasca, Marzio, and Gallo (2020); Iorio, Camisasca, and Gallo (2019); Iorio
_et al._ (2020); Tenuzzo, Camisasca, and Gallo (2020); Calero and Franzese
(2020) and in network topology Martelli, Crain, and Franzese (2020);
Chiricotto  _et al._ . Therefore, the choice of a given interaction potential
becomes of particular relevance.
In the class of 3 points models, we have tested the TIP3P, the SPC, the SPC/E
and the flexible SPC models. We have found that the TIP3P model is
characterized by the less structured network, with broad distributions or
rings and a network rich in coordination defects which allows the network to
arrange in long rings. In particular, the counting scheme d3 gives a low value
of complexity index $\xi$ which reflects the large distance from the HBN of
crystalline ice for which $\xi=1.0$. The broad rings distribution and the low
percentage of intact HBs explain the absence of a second hydration peak in the
g2(r). The percentage of intact HBs increases in the SPC model which is,
therefore, characterized by an HBN with fewer longer rings and by a g2(r) with
signatures of a second hydration peak. The SPC/E and the flexible SPC models
are characterized by a further increase in the percentage of intact HBs,
comparable with that of _ab initio_ liquid water. The resulting HBNs
accommodate an even lower percentage of longer rings and, therefore, the
network complexity index is higher compared to the previous models, indicating
a closer (but still very far) HBN to the HBN of Ih(c), in agreement with the
more structured g2(r). It is of particular interest to observe how the
introduction of flexibility in the SPC model drastically affects the topology
of the HBN.
In the class of 4 points models, we have tested the TIP4P, the TIP4P-Ew, the
TIP4P/2005, the flexible TIP4P/2005 and the TIP4P-Ice models. The TIP4P model
is the only model whose network accommodates a lower percentage of intact HBs
with respect to _ab initio_ liquid water. The topology of the corresponding
HBN is hence the most complex, i.e., with a low complexity index. The
TIP4P-Ew, the TIP4P/2005 and the flexible TIP4P/2005 models are characterized
by a similar percentage of intact HBs, comparable with that of _ab initio_
liquid water. The corresponding HBNs have comparable topologies and values of
complexity indices. The TIP4P-Ice model, finally, shows a higher percentage of
intact HBs with respect to _ab initio_ liquid water. The topology of the
corresponding HBN is the less complex among all models here studied, with
small contributions of longer rings and the highest value of $\xi$ index. Such
results explain the over structured g2(r) with respect to both _ab initio_
water and experimental results.
In the class of 5 points models, we have tested the TIP5P and the TIP5P-E
models. Both models have similar HBNs, characterized by a lower percentage of
broken HBs with respect to _ab initio_ liquid water. Both networks are fairly
complex, accommodating longer rings causing low values of the network
complexity index. Overall, the balance between intact HBs and network topology
allows the 5 points models to be the better models in reproducing the _ab
initio_ and experimental water g2(r).

Overall, we have shown that water models endowed with the fastest dynamics are
able to establish more complex networks, while models with the slowest
dynamics establishes networks more closely related to that of cubic or
hexagonal ice. In particular, among the 11 models here inspected the TIP3P and
the SPC models are the ones with the fastest dynamics and with the network
deviating the most from that of ice. On the other hand, the TIP4P-Ice model is
the one with the slowest dynamics and, hence, with a network of HB the closer
to that of ice.

Finally, we have shown that the topology of the HBN might be affected by
finite size effects when other observables such as, e.g., the two body pair
correlation function, do not show such sensitivity.

In conclusion, the topology of the HBN and its quality in terms of broken and
intact HBs are more sensitive quantities than other physical observables
Martelli, Crain, and Franzese (2020); Chiricotto  _et al._ . Therefore, the
properties of the HBN should be inspected along with all other properties such
as, e.g., structural, dynamical and thermodynamic properties when developing
new interaction potentials. Our study provides a benchmark evaluated following
three different ring definitions and counting schemes. New interaction
potentials should be tested against such results. Nonetheless, the network
complexity index provides a direct quantitative measure of how much a HBN is
complex and far from the HBN at the ground state, and a direct link to the
dynamical properties. Such quantity can be transferred to other materials. The
effect of polarization on the topology of the HBN and on its quality should be
investigated.

###### Acknowledgements.

We acknowledge support from the STFC Hartree Centre’s Innovation Return on
Research programme, funded by the Department for Business, Energy and
Industrial Strategy.

## References

  * Salzmann (2019) C. G. Salzmann, J. Chem. Phys. 150, 060901 (2019).
  * Palmer  _et al._ (2014) J. C. Palmer, F. Martelli, Y. Liu, R. Car, A. Z. Panagiotopoulos, and P. G. Debenedetti, Nature 510, 385 (2014).
  * Sellberg  _et al._ (2014) J. A. Sellberg, C. Huang, T. A. McQueen, N. D. Loh, H. Laksmono, D. Sclesinger, R. G. Sierra, D. Nordlund, C. Y. Hampton, D. Starodub, D. P. DePonte, M. Beye, C. Chen, A. V. Martin, A. Barty, K. T. Wikfeldt, T. M. Weiss, C. Caronna, J. Feldkamp, L. B. Skinner, M. M. Seibert, M. Messerschmidt, G. J. Williams, S. Boutet, L. G. M. Pettersson, M. J. Bogan, and A. Nilsson, Nature 510, 381 (2014).
  * Debenedetti, Sciortino, and Zerze (2020) P. G. Debenedetti, F. Sciortino, and G. H. Zerze, Science 369, 289 (2020).
  * Kringle  _et al._ (2020) L. Kringle, W. A. Thornley, B. D. Kay, and G. A. Kimmel, Science 369, 1490 (2020).
  * Kim  _et al._ (2020) K. H. Kim, K. Amann-Winkel, N. Giovambattista, A. Spah, F. Perakis, H. Pathak, M. L. Parada, C. Yang, D. Mariedahl, T. Eklund, T. J. Lane, S. You, S. Jeong, M. Weston, J. H. Lee, I. Eom, M. Kim, J. Park, P. P. S.H. Chun, and A. Nilsson, Science 370, 978 (2020).
  * Poole  _et al._ (1992) P. H. Poole, F. Sciortino, U. Essmann, and H. E. Stanley, Nature 360, 324 (1992).
  * Liu  _et al._ (2010) Y. Liu, J. C. Palmer, A. Z. Panagiotopoulos, and P. G. Debenedetti, J. Chem. Phys. 137, 214505 (2010).
  * Limmer and Chandler (2011) D. T. Limmer and D. Chandler, J. Chem. Phys. 135, 134503 (2011).
  * Wikfeldt, Nilsson, and Pettersson (2011) K. T. Wikfeldt, A. Nilsson, and L. G. M. Pettersson, Phys. Chem. Chem. Phys. 13, 19918 (2011).
  * Palmer, Car, and Debenedetti (2013) J. C. Palmer, R. Car, and P. G. Debenedetti, Faraday Discuss. 167, 77 (2013).
  * Limmer and Chandler (2013) D. T. Limmer and D. Chandler, J. Chem. Phys. 138, 214504 (2013).
  * Limmer and Chandler (2014) D. T. Limmer and D. Chandler, Proc. Natl. Acad. Sci. USA 111, 9413 (2014).
  * Chandler (2016) D. Chandler, Nature 531, E1 (2016).
  * Palmer  _et al._ (2016a) J. C. Palmer, F. Martelli, Y. Liu, R. Car, A. Z. Panagiotopoulos, and G. D. P, Nature 531, E2 (2016a).
  * Palmer  _et al._ (2018) J. C. Palmer, A. Haji-Akbari, R. S. Singh, F. Martelli, R. Car, A. Z. Panagiotopoulos, and P. G. Debenedetti, J. Chem. Phys. 148, 137101 (2018).
  * Palmer  _et al._ (2016b) J. C. Palmer, R. S. Singh, R. Chen, F. Martelli, and P. G. Debenedetti, Mol. Phys. 114, 2580 (2016b).
  * Martelli  _et al._ (2020) F. Martelli, F. Leoni, F. Sciortino, and J. Russo, J. Chem. Phys. 153, 104503 (2020).
  * Shi and Tanaka (2020a) R. Shi and H. Tanaka, Proc. Natl. Ac. Sci. 117, 26591 (2020a).
  * Shi and Tanaka (2020b) R. Shi and H. Tanaka, J. Am. Chem. Soc. 142, 2868 (2020b).
  * Shi, Russo, and Tanaka (2018a) R. Shi, J. Russo, and H. Tanaka, J. Chem. Phys. 149, 224502 (2018a).
  * Russo and Tanaka (2014a) J. Russo and H. Tanaka, Nat. Commun. 5, 3556 (2014a).
  * Akahane and Tanaka (2018) J. R. K. Akahane and H. Tanaka, Proc. Natl. Ac. Sci. 115, E3333 (2018).
  * Santra  _et al._ (2015) B. Santra, R. A. DiStasio Jr., F. Martelli, and R. Car, Mol. Phys. 113, 2829 (2015).
  * Huang  _et al._ (2009) C. Huang, K. T. Wikfeldt, T. Tokushima, D. Nordlund, Y. Harada, U. Bergmann, M. Niebuhr, T. M. Weiss, Y. Horikawa, M. Leetmaa, M. P. Ljungberg, O. Takahashi, A. Lenz, L. Ojamäe, A. P. Lyubartsev, S. Shin, L. G. M. Pettersson, and A. Nilsson, Proc. Natl. Acad. Sci. USA 106, 15214 (2009).
  * Nilsson and Pettersson (2015) A. Nilsson and L. G. M. Pettersson, Nat. Commun. 6, 8998 (2015).
  * De Marzio  _et al._ (2017) M. De Marzio, G. Camicasca, M. Rovere, and P. Gallo, J. Chem. Phys. 146, 084502 (2017).
  * Martelli (2019) F. Martelli, J. Chem. Phys. 150, 094506 (2019).
  * Pekka and Lennart (2001) M. Pekka and N. Lennart, J. Phys. Chem. A 105, 9954 (2001).
  * Mao and Zhang (2012) Y. Mao and Y. Zhang, Chem. Phys. Lett. 542, 37 (2012).
  * Lee and Kim (2019) S. H. Lee and J. Kim, Mol. Phys. 117, 1926 (2019).
  * Jorgensen  _et al._ (1983) W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 (1983).
  * Harrach and Drossel (2014) M. F. Harrach and B. Drossel, J. Chem. Phys. 140, 174501 (2014).
  * González  _et al._ (2010) B. S. González, E. G. Noya, C. Vega, and L. M. Sesé, J. Phys. Chem. B 114, 2484 (2010).
  * Zielkiewicz (2005) J. Zielkiewicz, J. Chem. Phys. 123, 104501 (2005).
  * Steinczinger, Jóvári, and Pusztai (2017) Z. Steinczinger, P. Jóvári, and L. Pusztai, J. Mol. Phys. 228, 19 (2017).
  * Dix, Lue, and Carbone (2018) J. Dix, L. Lue, and P. Carbone, J. Comp. Chem. 39, 2051 (2018).
  * Bernal and Fowler (1933) J. D. Bernal and R. H. Fowler, J. Chem. Phys. 1, 515 (1933).
  * Tse  _et al._ (1999) J. S. Tse, D. D. Klug, C. A. Tulk, I. Swainson, E. C. Svensson, C.-K. Loong, V. Shpakov, V. R. Belosludov, R. V. Belosludov, and Y. Kawazoe, Nature 400, 647 (1999).
  * Marton̆ák, Donadio, and Parrinello (2004) R. Marton̆ák, D. Donadio, and M. Parrinello, Phys. Rev. Lett. 92, 225702 (2004).
  * Marton̆ák, Donadio, and Parrinello (2005) R. Marton̆ák, D. Donadio, and M. Parrinello, J. Chem. Phys. 122, 134501 (2005).
  * Shephard  _et al._ (2017) J. J. Shephard, S. Ling, G. Sosso, A. Michaelides, B. Slater, and C. G. Salzmann, J. Phys. Chem. Lett. 8, 1645 (2017).
  * Martelli  _et al._ (2018) F. Martelli, N. Giovambattista, S. Torquato, and R. Car, Phys. Rev. Materials 2, 075601 (2018).
  * Martelli, Crain, and Franzese (2020) F. Martelli, J. Crain, and G. Franzese, ACS Nano 14, 8616 (2020).
  * (45) M. Chiricotto, F. Martelli, G. Giunta, and P. Carbone, J. Phys. Chem. C Under review.
  * Martelli, Calero, and Franzese (2021) F. Martelli, C. Calero, and G. Franzese, “Re-defining the concept of hydration water near soft interfaces,” (2021), arXiv:2101.06136 [cond-mat.soft] .
  * Berendsen  _et al._ (1981) H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, and J. Hermans, in  _Intermolecular Forces_ , edited by B. Pullman (Springer, Dordrecht, 1981) pp. 331–342.
  * Berendsen, Grigera, and Straatsma (1987) H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 (1987).
  * Toukan and Rahman (1985) K. Toukan and A. Rahman, Phys. Rev. B 31, 2643 (1985).
  * Amira, Spángberg, and Hermansson (2004) S. Amira, D. Spángberg, and K. Hermansson, Chem. Phys. 303, 327 (2004).
  * Abascal  _et al._ (2005) J. L. F. Abascal, E. Sanz, R. G. Fernández, and C. Vega, J. Chem. Phys. 122, 234511 (2005).
  * Abascal and Vega (2005) J. L. F. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005).
  * González and Abascala (2011) M. A. González and J. L. F. Abascala, J. Chem. Phys. 135, 224516 (2011).
  * Horn  _et al._ (2004) H. W. Horn, W. C. Swope, J. W. Pitera, J. D. Madura, T. J. Dick, G. L. Hura, and T. Head-Gordon, J. Chem. Phys. 120, 9665 (2004).
  * Mahoney and Jorgensen (2000) M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 122, 8910 (2000).
  * Rick (2004) S. W. Rick, J. Chem. Phys. 120, 6085 (2004).
  * Martelli  _et al._ (2016) F. Martelli, H.-Y. Ko, E. C. Oğuz, and R. Car, Phys. Rev. B 97, 064105 (2016).
  * Formanek and Martelli (2020) M. Formanek and F. Martelli, AIP Adv. 10, 055205 (2020).
  * Leoni  _et al._ (2019) F. Leoni, R. Shi, H. Tanaka, and J. Russo, J. Chem. Phys. 151, 044505 (2019).
  * Camisasca  _et al._ (2019) G. Camisasca, D. Schlesinger, I. Zhovtobriukh, G. Pitsevich, and L. G. M. Pettersson, J. Chem. Phys. 151, 034508 (2019).
  * Russo and Tanaka (2014b) J. Russo and H. Tanaka, Nat. Commun. 5, 3556 (2014b).
  * Fitzner  _et al._ (2019) M. Fitzner, G. C. Sosso, S. J. Cox, and A. Michaelides, Proc. Natl. Acad. Sci. USA 116, 2009 (2019).
  * Shi and Tanaka (2018) R. Shi and H. Tanaka, Proc. Natl. Acad. Sci. USA 115, 1980 (2018).
  * Bakó  _et al._ (2017) I. Bakó, J. Oláh, A. Lábas, S. Bálint, L. Pusztai, and M. C. B. Funel, J. Mol. Liq. 228, 25 (2017).
  * Pothoczki, Pusztai, and Bakó (2018) S. Pothoczki, L. Pusztai, and I. Bakó, J. Phys. Chem. B 122, 6790 (2018).
  * Pothoczki, Pusztai, and Bakó (2019) S. Pothoczki, L. Pusztai, and I. Bakó, J. Phys. Chem. B 123, 7599 (2019).
  * Li  _et al._ (2020a) L. Li, J. Zhong, Y. Yan, J. Zhang, J. Xu, J. S. Francisco, and X. C. Zeng, Proc. Natl. Acad. Sci. USA 117, 24701 (2020a).
  * DiStasio Jr.  _et al._ (2014) R. A. DiStasio Jr., B. Santra, Z. Li, X. Wu, and R. Car, J. Chem. Phys. 141, 084502 (2014).
  * Nosé (1984) S. Nosé, Mol. Phys. 52, 255 (1984).
  * Hoover (1985) W. G. Hoover, Phys. Rev. A 31, 1695 (1985).
  * Parrinello and Rahman (1981) M. Parrinello and A. Rahman, J. Appl. Phys. 52, 7182 (1981).
  * Abraham  _et al._ (2015) M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess, and E. Lindahl, SoftwareX 1, 19 (2015).
  * King (1967) S. V. King, Nature 213, 1112 (1967).
  * Rahman and Stillinger (1973) A. Rahman and F. H. Stillinger, J. Am. Chem. Soc. 95, 7943 (1973).
  * Guttman (1990) L. Guttman, J. Non-Cryst. Solids 116, 145 (1990).
  * Franzblau (1991) D. S. Franzblau, Phys. Rev. B 44, 4925 (1991).
  * Wooten (2002) F. Wooten, Acta Cryst. A 58, 346 (2002).
  * Yuan and Cormack (2002) X. Yuan and A. N. Cormack, Comp. Mater. Sci. 24, 343 (2002).
  * Roux and Jund (2010) S. L. Roux and P. Jund, Comp. Mater. Sci. 49, 70 (2010).
  * Jin  _et al._ (1994) W. Jin, R. K. Kalia, P. Vashishta, and J. P. Rino, Phys. Rev. B 118, 50 (1994).
  * Hobbs  _et al._ (1998) L. W. Hobbs, C. E. Jesurum, V. Pulim, and B. Berger, Philos. Mag. A 68, 679 (1998).
  * Marians and Hobbs (1988) C. S. Marians and L. W. Hobbs, J. Non-Cryst. Solids 160, 317 (1988).
  * Marians and Hobbs (1990) C. S. Marians and L. W. Hobbs, J. Non-Cryst. Solids 119, 269 (1990).
  * Luzar and Chandler (1996) A. Luzar and D. Chandler, Nature 379, 55 (1996).
  * Prada-Gracia, Shevchuk, and Rao (2013) D. Prada-Gracia, R. Shevchuk, and F. Rao, J. Chem. Phys. 139, 084501 (2013).
  * Shi, Russo, and Tanaka (2018b) R. Shi, J. Russo, and H. Tanaka, J. Chem. Phys. 149, 224502 (2018b).
  * Li  _et al._ (2020b) L. Li, J. Zhong, Y. Yan, J. Zhang, J. Xu, J. S. Francisco, and X. C. Zeng, Proc. Natl. Acad. Sci. USA 117, 24701 (2020b).
  * de Oca  _et al._ (2020) J. M. M. de Oca, F. Sciortino, , and G. A. Appignanesi, J. Chem. Phys. 152, 244503 (2020).
  * Skinner  _et al._ (2013) L. B. Skinner, C. Huang, D. Schlesinger, L. G. M. Pettersson, A. Nilsson, and C. J. Benmore, J. Chem. Phys. 138, 074506 (2013).
  * Soper and Benmore (2008) A. K. Soper and C. J. Benmore, Phys. Rev. Lett. 101, 065502 (2008).
  * Praprotnik, Janežič, and Mavri (2004) M. Praprotnik, D. Janežič, and J. Mavri, J. Phys. Chem. A 108, 11056 (2004).
  * MacKerell Jr.  _et al._ (1998) A. D. MacKerell Jr., D. Bashford, M. Bellott, R. L. D. Jr, J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiórkiewicz-Kuczera, D. Yin, and M. Karplus, J. Phys. Chem. B 102, 3586 (1998).
  * Calero, Stanley, and Franzese (2016) C. Calero, E. H. Stanley, and G. Franzese, Materials 9, 319 (2016).
  * Samatas  _et al._ (2018) S. Samatas, C. Calero, F. Martelli, and G. Franzese, arXiv:1811.01911 [cond-mat.soft] (2018).
  * Gallo, Rovere, and Chen (2010) P. Gallo, M. Rovere, and S.-H. Chen, J. Phys. Chem. Lett. 1, 729 (2010).
  * Camisasca, Marzio, and Gallo (2020) G. Camisasca, M. D. Marzio, and P. Gallo, J. Chem. Phys. 153, 224503 (2020).
  * Iorio, Camisasca, and Gallo (2019) A. Iorio, G. Camisasca, and P. Gallo, Sci. China Phys. Mech. Astron. 62, 107011 (2019).
  * Iorio  _et al._ (2020) A. Iorio, M. Minozzi, G. Camisasca, M. Rovere, and P. Gallo, Philosophical Magazine 100, 2582 (2020).
  * Tenuzzo, Camisasca, and Gallo (2020) L. Tenuzzo, G. Camisasca, and P. Gallo, Molecules 25, 4570 (2020).
  * Calero and Franzese (2020) C. Calero and G. Franzese, J. Mol. Liq. 317, 114027 (2020).