We consider N water molecules at a constant temperature T and pressure P within a fluctuating volume V ( T , P ) . We decompose the total volume V into two components: homogeneous and heterogeneous. The homogeneous component, V iso ( T , P ) , corresponds to the isotropic contribution arising from the van der Waals interaction between the water molecules. This component is modeled using the Lennard-Jones potential in Equation 1, U r ≡ 4 ϵ r 0 r 12 − r 0 r 6 − U c , if  r < r c 0 if  r ≥ r c , (1) with a hard-core distance r 0 and a cutoff at r c ≡ 6 r 0 , where r 0 ≡ 2.9 Å is the van der Waals diameter of a single water molecule, associated to its van der Waals volume v 0 , as determined from experiments (Finney, 2001). We shift the potential by U c ≡ 4 ϵ r 0 / r c 12 − r 0 / r c 6 to avoid a discontinuity at r c . The cutoff is chosen large enough to include all significant contributions to the van der Waals interactions. We take ϵ , the characteristic energy of the Lennard-Jones interaction, as the internal unit of energy. From ab initio energy calculations (Henry, 2002), we set ϵ ≡ 5.5 kJ/mol.

The heterogeneous component of the volume reflects local fluctuations resulting from the formation of HBs under specific thermodynamic conditions. Sastry et al. demonstrated (Sastry et al., 1996) that assuming these fluctuations are proportional to the total number of HBs, N HB , is sufficient to reproduce water’s volumetric anomalies. Thus, the total volume is expressed as V ≡ V iso + v HB N HB , (2) where the proportionality factor v HB remains independent of the thermodynamic conditions ( T , P ) . This assumption is made to simplify the model and is shown to be reasonable a posteriori, at least within a limited range of T and P . We will further discuss this limitation before the conclusions.

We set v HB ≡ 0.6 v 0 . This choice stems from the volume difference per HB between high-density ice VI and VIII and low-density tetrahedral ice Ih (Bianco and Franzese, 2014; Coronas et al., 2022). It is based on the reasonable assumption that the difference between the low- and high-density ices is solely due to the open structure associated with tetrahedral HB formation and that in low-density ice, all HBs are formed, with each water molecule engaging in four HBs. In contrast, in high-density ices, all HBs are absent.

Our strategy for achieving large-scale simulation capability involves reducing the degrees of freedom of water without losing information about the HB network. To this end, we replace the atomic coordinates of the N molecules with a density field defined at the resolution of a single molecule. Therefore, we partition the homogeneous component of the total volume V iso into N cells, each sized v iso ≡ V iso / N ≥ v 0 , and ensure that N ≥ N . The case N = N , examined here, corresponds to bulk water, with each cell accommodating a single molecule and v iso representing the proper volume of a water molecule that forms no HBs. The case N > N indicates that the volume is shared by water and solutes, which will be addressed in a future publication for the 3D case (Coronas, 2023). In 2D, the case N > N has been considered already when there are vacancies in a monolayer (Franzese and de los Santos, 2009; de los Santos and Franzese, 2011; 2012) or for hydrated proteins, e.g., in (Durà-Faulí et al., 2023) and references therein.

To account for local volume fluctuations caused by the formation of HBs, we associate a heterogeneous component v HB n HB , i / 2 with each cell i ∈ [ 1 , … , N ] . Here, n HB , i represents the number of HBs formed by the molecule within cell i , with ∑ i n HB , i / 2 = N HB . The factor of 1 / 2 avoids double-counting of HBs. Consequently, each cell i has a local volume defined as v i ≡ v iso + v HB n HB , i / 2 , (3) that implies Equation 2.

To keep our coarse-graining approach (i) straightforward to implement and (ii) consistent with both low and high coordination numbers in the fluid, we partition the total volume V into cells of a cubic lattice. As a result, the relation between the van der Waals diameter r 0 and the associated volume is v 0 ≡ r 0 3 .

The partition of the volume V into a cubic grid of cells is appropriate because DFT-based Car-Parrinello molecular dynamics simulations indicate that the water coordination number does not exceed six under ambient conditions (Skarmoutsos et al., 2022). Simulations at high pressures also confirm this finding (Saitta and Datchi, 2003; Paschek et al., 2008). Therefore, each cell i has six nearest neighbors.

The average distance between neighboring molecules is defined as r ≡ v iso 1 / 3 ∈ [ r 0 , ∞ ) . Note that r is unaffected by N HB , as the formation of HBs reduces the coordination number of water but does not alter the separation between molecules. Specifically, each water molecule minimizes the enthalpy of its local environment by forming four HBs in an almost perfect tetrahedral arrangement while excluding any ‘interstitial’ water molecule. This rearrangement leads to a decrease in local density, which corresponds with an increase in the effective volume of each molecule in the network, resulting in a variation in the total V as described by Equation 2.

The system is compressible, meaning that V fluctuates at fixed P in accordance with the equation of state. Therefore, for each i , the two components of v i ( v iso and v HB n HB , i / 2 ) vary. For the first, independent of i , it holds that v 0 / v iso = r 0 3 / r 3 ∈ [ 0,1 ] . By defining r 1 / 2 as the value where r 0 3 / r 1 / 2 3 = 0.5 , we classify as gas-like the cells with r > r 1 / 2 , i.e., those with v iso ≥ 2 v 0 , and as liquid-like the others. Since this definition is independent of i , the entire system is either gas or liquid-like, depending on the value of V iso . Nevertheless, local changes in the HB network, via n HB , i , lead to heterogeneities in local volume fluctuations.

We consider negligible the HB formation in the gas and assume that molecules within gas-like cells cannot form HBs since the average O–O distance between them exceeds the HB-breaking threshold (Luzar and Chandler, 1996). In addition, according to ab initio simulations and the Debye-Waller factor (Teixeira et al., 1990), only O O H ̂ angles (between two water molecules) within 60° result in a bonded state. Therefore, only 1 / 6 of the possible relative orientation states of two water molecules at HB distance can form a HB. To account for this, we introduce a bonding variable σ i j ∈ [ 1 , … , q ] with q = 6 . It describes the relative orientation between molecules in neighboring cells i and j . Each molecule i has six bonding variables, σ i j , one for each of the six neighboring molecules j . A HB is formed when σ i j = σ j i , in such a way that a bonded state has a probability 1 / q = 1 / 6 to occur.

As discussed in (Coronas and Franzese, 2024; Coronas et al., 2025), the model splits the HB interaction into two components: (i) covalent (pairwise directional) (Shi et al., 2018), and (ii) cooperative (many-body) (Barnes et al., 1979), with a characteristic energies J and J σ , respectively. We set J ≡ 11 kJ/mol, i.e., J / ( 4 ϵ ) ≡ 0.5 , which is consistent with the energy of a single HB and cluster analysis (Stokely et al., 2010).

The cooperative HB interactions arise from many-body effects, contributing to the tetrahedral arrangement of HBs at low temperatures (Cisneros et al., 2016). The model includes interactions up to the five-body term within the first coordination shell, as derived from polarizable models (Abella et al., 2023). Based on DFT calculations (Cobar et al., 2012), we set J σ / ( 4 ϵ ) = 0.08 to optimize the model’s accuracy in predicting the experimental equation of state and thermodynamic fluctuations (Coronas et al., 2025).

The six bonding variables σ i j of the same molecule i interact cooperatively, resulting in a reduction in energy J σ when σ i j = σ i k for j ≠ k . With a coordination number of six, this implies that the maximum cooperative energy in a cell is 15 J σ . For the selected model’s parameters (Table 1), it follows that J σ ≪ J , which aligns with the general understanding that HB cooperative reorganization occurs at a temperature significantly lower than the formation of individual HBs (Cisneros et al., 2016).

Both experimental and computational studies indicate that bulk water molecules form four tetrahedral HBs in their lowest energy state. Excited states correspond to defects in the HB network, where molecules have either fewer or more HBs. Ab initio calculations for liquid water at ambient conditions show that under-coordinated molecules comprise approximately 45% of the HB network, while over-coordinated molecules account for less than 5% (DiStasio et al., 2014). According to neutron scattering and Raman spectroscopy studies (Giguère, 1987), bifurcated hydrogen bonds—where one acceptor interacts with two donors—could potentially enable the formation of up to five HBs. However, the energy of each bifurcated HB is approximately half that of a canonical HB. The current consensus suggests that molecules can rapidly switch their HBs between two neighboring molecules within their coordination shell, within hundreds of femtoseconds (Laage and Hynes, 2006). This dynamics, which resembles, on average, a bifurcated HB, is included in our model. Therefore, to simplify, we impose a maximum of four HBs per molecule, introducing variables η i j , which are set to 1 or 0 depending on whether the HB between molecules i and j is allowed, as described in Supplementary Material: Checkerboard partition for η variables.

Finally, the enthalpy of the system is given by: H N , P , T = U LJ − N HB J − N σ J σ + P V , (4) where U LJ is the Lennard-Jones potential in Equation 1 U LJ (r), N HB ≡ ∑ ⟨ i , j ⟩ θ ( r 1 / 2 − r i , j ) δ ( σ i , j , σ j , i ) η i , j , N σ ≡ ∑ k N ∑ ⟨ i , j ⟩ δ ( σ k , j , σ k , i ) , and V is given by Equation 2. Here, r i , j ≡ | r i ⃗ − r j ⃗ | is the distance between molecules i and j , and θ ( r ) and δ ( i , j ) are the Heaviside step and Kronecker delta functions, respectively. Thus, the formation of a macroscopic HB network leads to an increase in volume for Equation 2, a decrease in entropy due to the reduced number of accessible σ i j states, and an increase in HB enthalpy, for Equation 4.

A configuration of the CVF model is defined by the variables { r i ⃗ , σ i j , η i j } . In the present version of the model, as discussed above, we coarse-grain the molecule position r i ⃗ over the lattice cell and assign to each molecule a proper volume v i , Equation 3. Therefore, considering the definitions of V iso and η i j , the CVF configuration of the present model reduces to { V iso , σ i j , η i j } .

In each MC step, we update these variables in the following order:

1. Update V iso , keeping { σ i j , η i j } fixed.

We use the standard Metropolis algorithm to update the global variable V iso . This method involves accepting or rejecting a tentative change from V iso to V iso + Δ V iso with a probability proportional to exp [ − P Δ V iso / ( k B T ) ] , where ± Δ V iso ∝ O ( V iso / 100 ) . We update the η i j as described in Supplementary Material: Checkerboard partition for η variables. For the σ i j , we employ the parallel local Metropolis algorithm (Section 2.3) or the parallel Swendsen-Wang cluster algorithm (Section 2.4), depending on the temperature: Metropolis for T ≥ 208 K and Swendsen-Wang for lower temperatures.

The Metropolis algorithm on a regular lattice can be efficiently parallelized by dividing the space into domains for simultaneous variable updates. To maintain detailed balance, the enthalpy change from altering a variable must remain independent of other variables within the same domain. For the Ising model, common partitioning schemes include layered (Barkema and MacFarland, 1994) and checkerboard (Heermann and Burkitt, 1990) methods, with CUDA implementations available for both 2D and 3D (Hawick et al., 2011; Weigel and Yavorskii, 2011; Wojtkiewicz and Kalinowski, 2015). However, these methods are not easily applicable to the CVF model due to differing lattice topologies, so we use a layered partition that enables memory coalescing in the CVF model.

We partition the { σ i j } variables, as described in Figure 1 (Top), into six domains. Each domain contains variables interacting with six bonding indices: five on the same molecule and one on a n.n. Molecule, all from different domains. Therefore, we can update all the σ i j in the same domain simultaneously.

In CUDA applications, the main bottleneck in execution arises from data access latency (Tapia and D’Souza, 2011). Performance can be enhanced by efficiently sorting memory to exploit memory coalescing (Leist et al., 2009; Sanders and Kandrot, 2010; NVIDIA, 2022). The GPU creates, manages, schedules, and executes blocks of 32 threads simultaneously, called warps (Hawick et al., 2011). When a kernel reads (or writes) to global memory locations, it performs a single coalesced read (or write) transaction for every half-warp of 16 threads. Therefore, we are interested in sorting the vectors so that consecutive threads read (or write) consecutive memory addresses.

We achieve this by sorting the arrays that store η i j and σ i j variables according to the index i d x = a r m ⋅ N + c e l l , where a r m ∈ { 0,1 , … , 5 } , and c e l l ∈ { 0,1 , … , N − 1 } . The index a r m represents the six possible neighbors of the cell (from 0 to 5: left, right, front, back, top, bottom), and c e l l ≡ i is the index of the cell (Figure 1 Bottom).

We implement a CUDA kernel gpu_metropolis(arm) that launches one thread per water molecule i , where a r m indicates which of the six independent domains is updated (Supplementary Algorithm 1). We define a parallel Metropolis update as six sequential calls to gpu_mertopolis(arm), where arm is chosen randomly to mimic the random selection of σ i j variables in the sequential Metropolis and to avoid the propagation of correlation waves.

We illustrate how the kernel gpu_metropolis(arm) performs coalesced memory transactions with the following example. We consider a half-warp that updates the block d i r = 0 (left domain) of the water cells c e l l ∈ { 0 , … , 15 } . Thus, idx takes values from 0 to 15. When the kernel estimates Δ N HB , it reads the right arms of the neighboring cells 1 to 16, i.e., the (consecutive) positions idx = 65 to 80. The same occurs when estimating Δ N σ , as the kernel reads memory positions in consecutive domains. We observe that an exception to this rule arises when the neighboring cell is positioned on the opposite side of the simulation box due to the periodic boundary conditions.

Local MC algorithms, such as Metropolis, experience a critical slowdown in their dynamics as the correlation length approaches the system size (as discussed in Section 3.1). In contrast, cluster MC algorithms efficiently update entire correlated regions of spins (clusters) simultaneously. Consequently, they produce statistically independent configurations at significantly lower computational costs. This efficiency is crucial in the supercooled region, for example, where the model exhibits a liquid-liquid phase transition culminating in a liquid-liquid critical point (Coronas and Franzese, 2024). In this context, we examine the Swendsen-Wang (SW) multi-cluster algorithm (Swendsen and Wang, 1987). The algorithm is defined so that, at each step, clusters of σ i j variables are formed with sizes ranging from 1 (an isolated σ i j variable) to the system’s size N . This formation follows a distribution that reproduces that of thermodynamically correlated degrees of freedom, as discussed in detail in (Bianco and Franzese, 2019) based on site-bond correlated percolation (Kasteleyn and Fortuin, 1969; Coniglio et al., 1979). The new configuration is generated by updating all the (correlated) σ i j variables within the same cluster to a new state. The sequential SW algorithm for the CVF model proceeds as follows:

1. Visit all the cells i . For each i , loop over all the pairs of variables ( σ i j , σ i k ) . If they are in the same state, place a fictitious bond between them with probability p σ = 1 − exp − J σ / k B T .

The SW algorithm performs three independent tasks. First, it places fictitious bonds between σ i j variables to generate the clusters. Second, it identifies all the clusters. Third, it updates each cluster. The first and third tasks are highly localized and can be easily parallelized. However, this is not the case for the cluster labeling operation. To tackle this challenge, we build on the work of Hawick et al., who developed various parallel labeling algorithms for arbitrary and lattice graphs using CUDA (Hawick et al., 2010). Among these, the label equivalence algorithm was refined by Kalentev et al. (Kalentev et al., 2011) and later applied by Komura and Okabe to SW simulations of the 2D Potts model (Komura and Okabe, 2012). In this context, we modify the Hawick-Kalentev label-equivalence algorithm for the CVF model.

For a given SW step, we first generate the clusters. We directly parallelize this task so that each thread works on one CVF cell. Each thread is responsible for the cooperative interactions within its cell and the covalent interactions with the left, front, and top directions. To accomplish this, we allocate the array connected of size ( 15 + 3 ) N = 18 N , which indicates whether two neighboring σ i j variables belong to the same cluster. We nest this array based on the index c o n _ i d x = l i n k N + c e l l , where l i n k ∈ { 0 , … , 17 } . l i n k = 0 , 1, and 2 represents the covalent connections between cell and its neighbors in the left, front, and top directions. The values l i n k ∈ { 3 , … , 17 } represent the 15 cooperative connections within c e l l .

Once the bonds are placed, we apply the label equivalence algorithm. We allocate the l a b e l array of size 6 N , which indicates the cluster that σ i j belongs to. Thanks to Kalentev’s sophistication, this array also resolves label equivalences (Kalentev et al., 2011). The advantage is the reduction of the memory cost of the algorithm, which is significant due to the limited storage resources of the GPUs. We initialize l a b e l as l a b e l [ i d x ] = i d x , where i d x is the σ i j index defined in Metropolis. The algorithm resolves label equivalences through iterative calls to the scanning and analysis functions (Kalentev et al., 2011; Komura and Okabe, 2012). When the algorithm converges, all the σ i j variables in the same cluster will take the same l a b e l value.

The scanning function compares the label of a site i d x to the labels of all the n.n. σ i j within the cluster. For every i d x , label [ i d x ] is updated to the minimum value among all the labels of the bonded sites, including itself. In Ref. (Komura and Okabe, 2012), Komura and Okabe implemented this function using a single kernel for the 2D Potts model. However, for the CVF model, we find it more convenient to divide this function into two kernels. First, in gpu_scanning_covalent, each thread scans left, front, and top covalent interactions (Supplementary Algorithm 2). Second, gpu_scanning_cooperative scans the cooperative interactions (Supplementary Algorithm 3). An alternative implementation in a single kernel leads to race conditions when two threads attempt to update the same element of label ¹ .

Next, the analysis function updates label [ i d x ] (Supplementary Algorithm 4). This step further propagates the minimum value of label to other σ i j variables within the same cluster. Although the parallel implementation of the analysis function experiences race conditions, these collisions between threads will eventually be resolved in subsequent applications of the scanning and analysis functions (Kalentev et al., 2011). To minimize the impact of thread conflicts, we implement the gpu_analysis ( a r m ) kernel, which updates only the label of those σ i j variables in the a r m domain. We then loop through the six domains to account for all the lattice sites.

To check whether the algorithm has converged, we first store a copy of the label vector before calling the scanning and analysis functions, and then we compare it to the updated label. We parallelize this task by assigning one thread to each CVF cell. The algorithm converges when the l a b e l remains unchanged. We provide an example of label convergence after successive applications of the scanning and analysis functions in Table 2.

The CVF model predicts a liquid-liquid phase transition (LLPT) between high-density liquid (HDL) and low-density liquid (LDL) phases in the supercooled region (Coronas and Franzese, 2024). The LLPT ends in a liquid-liquid critical point (LLCP), located at P C = 174 ± 14 MPa and T C = 186 ± 4 K in the thermodynamic limit ( N → ∞ ) . This is in close agreement with finite- N estimates from iAMOEBA (Pathak et al., 2016), TIP4P/Ice (Debenedetti et al., 2020), and ML-BOP (Dhabal et al., 2024) models, as well as with a recent estimate from a collection of experimental data (Mallamace and Mallamace, 2024).

Approaching the LLCP, the correlation length ξ of the water HB network increases and ultimately diverges at the critical point. Consequently, the autocorrelation time τ of the local Metropolis MC dynamics, which is proportional to ξ , also increases approaching the LLCP. This can be demonstrated by calculating the autocorrelation function C M Δ t ≡ 〈 M i t 0 + Δ t M t 0 〉 − 〈 M 〉 2 〈 M 2 〉 − 〈 M 〉 2 , (5) where M ≡ 1 q N max { N q ′ } is an order parameter, and N q ′ is the number of σ i j variables in the state q ′ ∈ { 0 , … , q − 1 } . We define the autocorrelation time τ as the time at which C ( τ ) = 1 / e . Thus, C ( Δ t ) allows us to estimate the autocorrelation time τ of the HB network. Two CVF configurations are uncorrelated if they are sampled after a number of MC steps ≥ τ .

We calculate, with the Metropolis algorithm, C M ( Δ t ) for a system with N = 32   768 water molecules (Figure 2a). At a pressure of P = 160 MPa, which is close to the critical pressure in the thermodynamic limit ( P = 174 ± 14 MPa), τ displays non-monotonic behavior with fast dynamics at high T = 208 K and low T = 188 K, alongside an apparent divergence at T = 193 K. This behavior is linked to a structural change between HDL-like and LDL-like forms of water, characterized by the Widom line (the locus of maxima of ξ ) emerging from the LLCP at higher pressure (Coronas and Franzese, 2024). As an approximate estimate of the Widom line, we present the locus of extrema of the specific heat C P , namely, the maxima of the enthalpy fluctuations (Figure 2b). At extremely low pressure, P = − 300 MPa, far from the critical region, τ exhibits a similar non-monotonic behavior, but without any apparent divergence upon crossing the Widom line. This is consistent with a decrease in the maxima of the correlation length ξ as the distance between the LLCP and the point along the Widom line increases.

Cluster MC algorithms are suitable for efficiently sampling the critical region, as they bypass the critical slowdown of the dynamics by updating regions of correlated HBs simultaneously. We compare the autocorrelation function (Equation 5) computed with local Metropolis (Figure 3a) and cluster SW (Figure 3b) algorithms. As discussed in Section 3.1, we observe slow dynamics of the system at low P = − 300 MPa and T = 205 K, near the Widom line. With increasing pressure ( P = 0.1 MPa) at constant temperature, the system remains in a metastable supercooled liquid state, exhibiting rapid decorrelation. Finally, at T = 195 K and P = 160 MPa, close to the LLCP (174 ± 14 MPa, 186 ± 4 K) (Coronas and Franzese, 2024), the autocorrelation time appears to diverge. The comparison with SW illustrates that cluster MC circumvents the critical slowdown of the dynamics in all cases, even near the LLCP.

Caution should be taken for the interpretation of time autocorrelation functions computed through MC simulations. Molecular Dynamics (MD) and MC fundamentally differ in how they explore the configurational phase space. While MD solves the time evolution of the system, MC proposes random updates of the system that are accepted or rejected with a probability given by the Boltzmann factor. Consequently, the MC timestep does not have a direct physical meaning. It is a computational step used to explore the configuration space, but it does not correspond to any real elapsed time. If MD is used, time correlation functions are a physical property of the system and describe how quickly a property, such as the HB lifetime, decorrelates. Conversely, when employing MC, time correlation functions reflect the algorithm’s ability to propose (and accept) configurations that vary with respect to a specific property. In the MC context, the time correlation function is not a property of the system but a property of the algorithm. Indeed, Metropolis and SW algorithms explore the same free energy landscape with identical equilibrium configurations; however, they differ in the number of steps needed to reach other equilibrium microstates. Hence, differences in MC time autocorrelation functions between Metropolis and SW do not reflect different physical properties of the simulated system but rather different algorithm efficiencies.

As discussed above, the autocorrelation time τ for the SW algorithm is significantly shorter than that for the Metropolis MC. However, SW cluster MC is notably more computationally expensive than Metropolis. Therefore, to determine which MC dynamics is more efficient in generating uncorrelated configurations, one must compare the time each algorithm takes to produce τ MC steps.

First, we analyze the computational cost of the parallel Metropolis algorithm for different system sizes, N ≤ 17   576   000 . The hardware and software specifications of the workstation are detailed in Supplementary Material: Workstation. Depending on the system size, we perform between 2 and 10 independent simulations of 5   000 MC steps. We find that the results are robust against changes in thermodynamic conditions; that is, changes in T and P do not affect the computational cost of the algorithm (Supplementary Algorithm 1).

Our results show that the time necessary (cost) for a parallel Metropolis update scales linearly for 32   768 ≤ N ≤ 2   097   152 (Figure 4). For these systems, the GPU resources are neither saturated (large N ) nor under-exploited (small N ); thus, the time spent on data accessing scales linearly with N . For small N ≤ 8   000 , the computational resources of the GPU are not optimized. We find that in this range, the time cost of a Metropolis step remains approximately constant ( ∼ 0.1 ms, Figure 4: inset). For large N ≥ 2   097   152 , the size of the arrays of random numbers must be reduced to fit within the GPU global memory (Supplementary Material: Generation and usage of random numbers). The additional time cost arises from both the increasing number of executions of the kernels for generating random numbers and the time involved in memory transactions. In particular, we benchmark accessible size-systems up to N = 17   576   000 water molecules with a time cost of 280 m per Metropolis update (Figure 4), which corresponds to a cubic simulation box of 75 × 75 × 75 nm 3 .

We estimate the size-dependent performance gain, or speedup factor, SF ( N ) ≡ t CPU / t GPU , defined as the ratio between the time required for a parallel and a sequential update of an entire system of N molecules, as shown in Table 3. The results indicate that, for the smallest system considered ( N = 64 molecules), the parallel algorithm is less efficient than the sequential one. This is not surprising, as a sufficiently large number of threads must be executed to fully utilize the GPU resources (Hall et al., 2014). Indeed, Wojtkiewicz and Kalinowski also find S F < 1 for small systems (Wojtkiewicz and Kalinowski, 2015). The large SF = 136.72 measured for N = 2   097   152 is attributed to the significant increase in the time cost of the sequential implementation compared to the parallel approach. More specifically, we find that the large σ and η arrays exceed the RAM storage capacity, necessitating that they be loaded in portions, which delays the sequential computation. We could not measure the SF for 2   097   152 < N ≤ 17   576   000 due to the excessive cost of computing t CPU . Instead, we extrapolated SF = 208 for N = 10   077   696 and SF = 245.8 for N = 17   576   000 from a power-law fit in the range 32   768 ≤ N ≤ 2   097   152 (see Supplementary Figure S3). Further details on the computation of the SF are provided in Supplementary Material: Speedup factors.

Next, we estimate the performance of the parallel SW algorithm. Unlike the Metropolis case, the time cost of a SW update depends on the cluster size distribution, which in turn is influenced by the thermodynamic conditions (Bianco and Franzese, 2019). Close to the Widom line, the system undergoes a transition from non-percolation to percolation upon isobaric cooling. Thus, we consider two temperatures on either side of the transition: T = 195 K (percolation) and 210 K (not percolation), with P = 0.1  MPa. For every system size N , we perform between 5 and 10 independent simulations of 1 0 3 MC steps.

We find that the time cost of the parallel SW algorithm increases linearly for N ≤ 262,144 , although data at small values of N are noisy (Figure 5). As with the parallel Metropolis algorithm, we attribute this to the suboptimal usage of GPU resources.

At larger values of N , we observe additional time costs compared to linearity. As in the Metropolis case, we attribute this to limited resources for storing large arrays, such as those used for random numbers. However, the cost value reached at N = 2   097   152 for SW is approximately 1.6 times larger than in the Metropolis case for the same size, which limits our ability to explore sizes with tens of millions of water molecules.

We note that the parallel SW update is faster under percolating conditions than in the absence of percolation. Although a better performance for a cluster algorithm is expected when the correlation length is large, because larger clusters lead to fewer in number, this result is not obvious. One might expect that the total time cost of the update is governed by the time cost of labeling the largest cluster, as seen in the sequential implementation.

A possible explanation is that the analysis function converges rapidly irrespective of the cluster size, making the size of the largest cluster less relevant. Therefore, the difference in time cost between percolation and non-percolation likely stems from less efficient memory readings of the label array in smaller clusters by the scanning and labeling functions.

A further consequence of this feature of the parallel implementation on GPUs is that the speedup factor relative to the sequential implementation on CPUs is greater under percolation conditions (Table 4; Supplementary Figure S4). In particular, we find that for N ≥ 32   768 , the SF under percolating conditions is nearly twice that under non-percolating conditions.

Based on our results, we observe that the time cost of a single MC update for a CVF water sample of N = 32   768 molecules (as discussed in Sections 3.1, 3.2) is 0.15 m using Metropolis. In the case of SW, the time cost is 1.9 m for percolating clusters, i.e., approaching the LLCP, and 2.0 m for non-percolating clusters, away from the critical region.

Therefore, the SW algorithm is approximately ten times more costly than Metropolis for N = 32   768 . This result suggests that the SW algorithm should be employed whenever the autocorrelation time τ obtained with SW is at least ten times shorter than the τ obtained with Metropolis. As we have discussed above, this occurs when the system approaches the Widom line (the maximum of the correlation length ξ ) or the region of maxima of specific heat (the maximum of enthalpy fluctuations).