{'fraction_non_alphanumeric': 0.035210940043722966, 'fraction_numerical': 0.01855900274431369, 'mean_word_length': 4.572835666148263, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Molecular Mechanism of Gas Solubility in Liquid: Constant Chemical Potential Molecular Dynamics Simulations

Accurate prediction of a gas solubility in a liquid is crucial in many areas of chemistry, and a detailed understanding of the molecular mechanism of the gas solvation continues to be an active area of research. Here, we extend the idea of constant chemical potential molecular dynamics (C{\mu}MD) approach to the calculation of the gas solubility in the liquid under constant gas chemical potential conditions. As a representative example, we utilize this method to calculate the isothermal solubility of carbon dioxide in water. Additionally, we provide microscopic insight into the mechanism of solvation that preferentially occurs in areas of the surface where the hydrogen network is broken.


Introduction
Solubility measures the ability of a substance in its solid, liquid, or gaseous state to dissolve into a solvent. Accurate prediction of solubility and a detailed understanding of the molecular mechanism of solute penetration across the interfaces is crucial in many areas of chemistry, geochemistry, inorganic, physical, organic, and biochemistry [1][2][3]. For example, drug solubility is a crucial parameter in the development of active pharmaceuticals [4,5]. In the petroleum industry, design of purification process and chemical separation depends on relative solubilities [6,7]. Solubility in water is of special interest since most fundamental biological, pharmacological, industrial, and environmental functions take place in aqueous solvent. [8][9][10][11][12].
Much effort has been devoted to developing successful computational tools for solubility prediction [13][14][15][16]. Free energy methods such as free energy perturbation [17] and thermodynamic integration [18] are widely used [19][20][21]. In these methods, the solvation free energy is computed from the free energy change associated with the transfer of a solute molecule from an ideal gas to the solvent. However, the main drawback of these methods is the neglect of the solute-solute interactions both in the gas and the solvent phases.
An alternative and more realistic way of approaching this problem is to use a standard vapor-liquid equilibrium (VLE) molecular dynamics (MD) simulation. In this method, the solute molecules are allowed to diffuse into the solvent, and the solubility is obtained from the equilibrium solute density inside the liquid. However, during the solvation process, the density of the solute in the vapor phase is depleted, causing a decrease in its chemical potential. This effect can eventually lead to an underestimation of the solubility. To avoid such an artifact, one could keep the pressure constant by using an NPT simulation. Unfortunately, for the system studied here, namely the solvation of CO 2 in water, one peculiar finite-size effect takes place. That is, when the number of gas molecules is small the system makes a transition to a layered structure in which water and liquid CO 2 alternate (see Section I of the SI). When the number of gas molecules is increased this effect disappears but this forces one to use a large system. Not only that, but given the high gas compressibility, one has to deal with large volume fluctuations.
A possible way of avoiding this artifact is to use a grand-canonical (µVT) ensemble. In this approach, the simulation box is coupled to a reservoir that supplies the solute molecules to the system, and thus, one is able to control the solute density via a Monte Carlo procedure. Here we explore an alternative to this approach using CµMD [22]. This is a method recently developed in our group in order to perform constant chemical potential simulations. The method has been successfully tested in the study of crystallization from solution [23][24][25][26], and gas permeation in organic membranes [27,28]. Here, we adopt this method to simulate the solvation of the gas in the liquid under constant chemical potential conditions.
To test the effectiveness of our CµMD method, we applied it to determine the isothermal solubility of CO 2 in water in a range of pressures and temperatures. The availability of a substantial amount of experimental and theoretical result makes this system a suitable test of the method [29][30][31][32][33][34][35][36][37]. Moreover, the CO 2 -H 2 O system is one of the fundamental geochemical binary systems, and knowledge of the gas solvation process is crucial for furthering our fundamental understanding and to advance technological and environmental applications [38,39].
An additional benefit of a two-fluid simulation is that one can obtain information on how the solute molecule dissolves in the liquid. This is a subject of great relevance to which much effort has been devoted [32,40,41,[41][42][43][44][45][46][47]. In this work, by analyzing the local curvature of water interface, we find that the distribution of water molecules is not uniform and depends on surface curvature. A detailed analysis reveals that near this instantaneous interface, there are transient empty patches. These patches are mostly localized at the crest of the interface and provide an entry point for the gas molecules into the liquid.

Computational Details
In this section, we first briefly describe the CµMD method followed by the details of the thermodynamic conditions and MD simulations.

CµMD Method
In this work, we use an orthorhombic simulation box and apply periodic boundary conditions in all directions. A liquid slab sits in the middle of the simulation box ( Fig. Figure 1(a) and (b)). The two sides of the liquid slab are filled with the gas molecules. The gas chamber is divided into different regions (see Fig. Figure 1(a)). On each side of the liquid slab we define control regions (CR, shown in orange in Fig. Figure 1(a)), where the number density of the gas is controlled by an external force (see below). The region between the liquid slab and the CR is called the transition 2/12 region (TR). Beyond the force region (FR) lies a molecular reservoir (Res) that maintains the CR density constant by supplying the extra molecules. The first step in the CµMD algorithm is to calculate the solute number density c in the CR as, where V CR is the CR volume and z j the distance of a j-th particle from the box center. The f (z j ) is a continuous and differentiable switching function that counts the number of atoms in the CR. In our case, this function is defined as a product of two Fermi switching functions, as: where, z in and z out are the inner and outer CR boundaries and α is a parameter that controls the switching functions steepness ( Fig. Figure 1(a)). The function f (z j ) has a value of 1 when the solute is inside the CR and 0 when outside. From now on, we drop the atomic index for simplicity. Now, we apply a force on the instantaneous solute density (c in Eq. 1) such that it remains to the target value (c 0 ), where κ is the force constant. The G(z) in Eq. (3) is a bell-shaped function that localizes the F (z) at z F . This function (G(z)) is defined as, where σ is a broadening parameter.
The size of each region, the number density of the solute in the reservoir and the force constant are among the parameters that need to be chosen. The regions (TR, CR, FR, and the Res) should be defined in such a way that they do not overlap. Details of the CµMD parameters used in our simulations are provided in Section II of the SI.

Potential Model and Thermodynamic Conditions
Besides using an appropriate computational tool, the accuracy of the result depends on the quality of the model interaction potential [36,[48][49][50][51]. An investigation of the performance of different models has been carried out by Liu et al. using histogram-reweighting Grand Canonical Monte Carlo (GCMC) [36]. Their results show that none of the potentials tested was able to reproduce the experimental results over the range of pressure and temperature investigated. In this paper, we use for water the SPC [52] model and EPM2 [52] for CO 2 . The Lorentz-Berthelot combination rules [53] were applied to determine the CO 2 -H 2 O interaction parameters. With this choice, a good agreement is obtained for T=423 K [36]. Guided by the Liu et al. experience we chose to investigate at this temperature and the pressure P=100 K (system I), where also experimental data is available, and the system is supercritical. We also moved close to the experimental critical point (T=304.13 K, P=73.8 bar) performing a second simulation at T=323 K, P=50 bar (system II), where CO 2 is in the vapor phase.

Simulation Protocol
To prepare the initial configurations, we follow three steps: 1) we simulate the bulk liquid in a cubic box of 4055 water molecules at the chosen thermodynamic conditions, 2) if a is the equilibrium value of the initial cubic box, the simulation box is then adiabatically elongated along the z-axis until the final orthorhombic box has dimension a × a × 4a (see Table S1), 3) the space left empty by water is filled by CO 2 molecules. The number of the gas molecules is fixed to the desired value.  Figure 1(b) shows a snapshot of the simulation box. The first two steps are identical in the CµMD and NVT-MD simulations, while in the third step the number of CO 2 molecules are different. In the NVT-MD simulations, the number of CO 2 molecules in the vapor phase is chosen so as to reproduce the experimental CO 2 gas density at the desired thermodynamic conditions. In the CµMD simulations, the number of CO 2 in the CR and TR is kept close to that used in the NVT-MD simulations, while in the Res this number is ∼1.3 times higher (see Table S1). As mentioned earlier, the CµMD simulations involve several parameters such as the size of TR, CR, FR, Res, and the value of κ. The CµMD parameters are chosen after performing a small number of simulations (see Table S2). All simulations were carried out using the GROMACS-2018.3 software [54]. The CµMD method is implemented in a private version of the PLUMED2.6 code [55]. The Ewald [56] summation was used for the long-range electrostatic interactions, and a cutoff of 10Å was used for both the Coulomb and the van der Waals interactions. A time step of 2 fs was used in all simulations. A pure bulk water system was equilibrated at each pressure for 10 ns in the NPT ensemble using the Parrinello-Rahman barostat [57] and the stochastic velocity rescaling thermostat [58]. For both thermostat and barostat, the coupling constant was 1 ps. Then, the vapor-liquid systems were equilibrated for 100 ps, and statistics was accumulated for 100 ns. In order to understand the properties of the interface, a 100 ns long simulation of the vapor-liquid water system was carried out.
For the solubility calculation, both water and CO 2 equilibrated densities were measured in the region of ∼3 nm width around the center of water density profile. The analysis was performed using 10 4 configurations extracted at 10 ps intervals from the equilibrium segment of the trajectory. The Visual Molecular Dynamics (VMD) software [59] was used to visualize the trajectories and produce some of the figures.

Gas Solubility
The equilibrium solubility of CO 2 in water is obtained from the time averaged spatial density profile in the liquid slab. It is worth emphasizing that due to the penetration of gas molecules in the liquid slab, there is a decrement in the liquid density, which is prominent in the high pressure and negligible for the low pressure cases.
In Fig. Figure 2 we illustrate the density profile of CO 2 in system I using the NVT-MD approach (dotted line). As described in the introduction, due to the solvation of CO 2 into water, the number of CO 2 in the gas phase is depleted and the vapor density is reduced by as much as ∼10 % and the solubility is also reduced. Our approach on the other hand keeps the vapor density at the desired value (see solid line in Fig. Figure 2), leading to the correct estimate. In system II this depletion is even larger ∼18 % (see Fig. S2 of SI).
We obtain for the solubility of CO 2 in system I, m CO2 =0.0140, which is in excellent agreement with the experimental data [29] (0.0135). In system II, our calculated solubility m CO2 =0.008 is in a good agreement with the theoretical result (∼0.008) of Ref. 36, but in disagreement with the experiment [36].

Molecular Mechanism of Solvation
In the previous section, we discussed the solubility of CO 2 and showed how using the CµMD method one can compute solubility using a relatively small system. We devote this section to a discussion on how the local interfacial surface corrugation plays a role in controlling the gas intake into the liquid. To this end, we first briefly summarize the Willard-Chandler (WC) approach [60] that we use to define the instantaneous interface. Subsequently, we investigate the hydrogen bond (HB) network of the interfacial water molecules in the vicinity of the surface. Finally, we reveal the emergence of empty patches at the interface and discuss the gas intake into the liquid. Here, we continue studying the same thermodynamic conditions as before (namely systems I and II).

Vapor-Liquid Equilibrium Interface
In the WC approach, the atomic density is broadened by placing a Gaussian at each atomic position. This broadened density is then collocated on a three-dimensional grid of points. The set of points from which the density is half of its bulk value density defines the interface. Here we constructed the interface by an array of 21 × 21 × 101 grid points which amounts to choosing a coarse-graining length to 2.4Å. For this analysis, we used the code implemented in Ref. 61.
In Fig. Figure 3(a), we show one such surface. The ondulation of this surface originates from capillary wave fluctuations. We shall show below that these ondulations have a significant effect on the local water structure and the gas adsorption process. In order to carry out this analysis, we need to identify the crests and troughs of the interface. These are defined by their position relative to the instantaneous average z-coordinate (Fig. Figure 3(a)).
In addition, we measure the thickness of this discretized surface as the difference in heights between the largest and smallest z-values. The interface area is obtained by first performing a Delaunay triangulation between the interface points and summing the triangles facet areas. In pure water the surface thickness increases by ∼50 % in going from system II where it is ∼8.5Å to system I that has a surface thickness of ∼12.5Å. In the case of the CO 2 -H 2 O system, the CO 2 molecules crowd at the interface (see Fig. Figure 1) and increase the surface thickness by ∼3Å in both systems.

HB Network in Vicinity of Surface Ondulation
Now we investigate the influence of water surface ondulation on the interfacial water network. The instantaneous density profile of water along the z-axis, depicted in Fig. S3, shows clearly the layering of water molecules, a well-known interfacial property amply discussed in the literature [46,60]. Furthermore, it has been shown that the density and orientational parameters of water molecules at the first layer deviate significantly from that of the bulk [62]. In order to analyze this effect, we identify the surface water molecules (W int ) as those whose distance from the interfacial grid points is less than ∼3.5Å. Subsequently, the surface water molecules are classified into the crest and trough molecules using the criteria described in the SI (Section III of the SI). Here the densities of vapor CO 2 and liquid water are at T=423 K and P=100 bar.

6/12
We are particularly interested in knowing the behavior of these surface water molecules and their HB network patterns across the crest and trough regions. Recently, Pezzotti et al. [46] identified a two dimensional (2D) extended HB network in this first layer of water near the interface with an average of 2.9 HB per water at 315 K and P=1 bar. We obtain an average of 2.4 and 2.7 HB per molecule for systems I and II, respectively. Here, we consider as hydrogen-bonded two water molecules that satisfy the geometrical criteria of O-O distance < 3.5Å and H-O...O angle <30 • . Fig. Figure 4(a) illustrates that among all water molecules of the 2D HB network, the number of water molecules that participate in this network is on average higher on the trough than the crest. At T=423 K, this difference is around 30 % (see Fig. Figure 5a), while at T=323 K this is around 10 % (See Fig. S5 in the SI). In other words, hydrogen bonding is more favorable on trough than in crest. This is due to the fact that the trough having a relatively higher surface area than crest can accommodate an extended HB network (See Fig S6 in the SI). This result affirms a non-uniform distribution of water molecules at the interface, and as we will discuss later, this kind of heterogeneity leads to the introduction of the surface empty patches concept.

Surface Empty Patches and Gas Solvation Mechanism
The concept of empty patch or cavity is well-defined in the literature and extensively used for examining the density fluctuations in bulk liquids [63][64][65][66][67][68][69][70][71][72]. In this work, for the first time, we characterize transient empty patches on the instantaneous surface.

7/12
In order to identify these patches, we draw on the interface circles of radius R P centered at the interface grid points. The union of all such circles that do not contain any water oxygen defines the surface empty patches (Fig. Figure 3(b)).
To do this, we use the graph algorithm implemented in Matlab [73]. A similar approach was used previously in the context of void analysis in the supercooled liquid water [72]. To define the surface empty patches area, we perform the Delaunay triangulation on a finer grid. For this analysis, we use R P =2.5Å. Fig. Figure 5(a) shows some representative small and large transient empty patches on the pure vapor-liquid water interface. As illustrated in Fig. Figure 5(a), a major portion of the empty patch area lies on the crest. There are only a few small empty patches on the trough. We now express this observation in a quantitative way. To this purpose, we calculate the fraction of free area (sum of all empty patches area of the interface divided by the total area of the instantaneous interface) that belongs to the crest and trough. Figs. Figure 4(b) and S5(b) reveal that in the vapor-liquid water, around 70 % of the total free area is localized on the crest, and only 30 % sites on the trough. This result is consistent with the presence of a larger number of HBs on the trough (Figs. Figure 5 and S5) than on the crest. Moreover, a similar percentage of the free area in CO 2 -H 2 O system (see dotted lines in Fig. Figure 4 panel (b)), reveals that the presence of gas molecules near the interface has no significant effect on the organization of water molecules in the vicinity of the instantaneous interface, which is consistent with the result of Ref. 44. Figs. Figure 4(c) and S5(b) show that also the HB percentage on the crest and trough remains almost unaltered, especially for the lower temperature T=323 K.
The distribution of CO 2 molecules within 2Å inside and outside of the interface (see panels (c) and (d) of Figs. Figure 4 and S5), shows that in both cases, the percentage of CO 2 molecules in the vicinity of the crest is higher than that of the trough. This suggests that the penetration of the gas molecules is more likely to take place from the crest via the interface. Panel (b) of Fig. Figure 5 shows two representative CO 2 molecules that penetrate the interface through the empty patches of the crest.

Conclusion
Scheme 1. Illustration of the diffusion mechanism of gas molecule through the interface.
In this study, we have applied the recently developed CµMD method to calculate the isothermal solubility of CO 2 in water at two thermodynamic conditions around the critical point of CO 2 . Proper control of the gas pressure near the liquid results in an accurate calculation of this important physical property. An excellent agreement of the calculated solubility with the experimental value is obtained at temperature 423 K and pressure 100 bar. The underestimation of the solubility data at lower temperature and pressure (T=323 K and P=50 bar) is due to the force-field's inability in reproducing experimental thermodynamic conditions. More accurate potentials, such as those derived from the first principle many-body potentials or accurate machine learning-based force fields, could provide a better estimate of the calculated solubility.
Additionally, we propose the mechanism of gas molecules solvation in the aqueous medium. The schematic in Scheme I illustrates a possible diffusion process of the gas molecules through the transient empty patches of the interface. An extended 2D-HB network on the trough region of the interface makes gas molecules difficult to penetrate through them to the interface. On the contrary, the large available free area on the crest region of the interface allows a facile diffusion of the gas molecules to the bulk liquid. As far as we know, this is the first time that we quantify the diffusion mechanism of gas molecules based on the empty patches on the local curvature of the instantaneous interface. Such microscopic details are important to understand the solvation mechanism of gases, in general solutes, in different solvents.