Robust methods for the characterization of droplet behavior in molecular dynamics: from contact radius to contact angle

Abstract

We present a comparative analysis of methodologies for studies of wetting behavior of liquid nano-droplets. Different droplet sizes were simulated using molecular dynamics (MD) under both complete and partial wetting conditions. The results show why conventional root mean square displacement (RMSD) metrics are inadequate for capturing finite system dynamics based on internal, center of mass, and coupling contributions to its value. The z-component of the center of mass is proposed as an alternative, accurate descriptor. A new equation for the contact radius of nano-droplets is derived using the radius of gyration (R_g) of interfacial molecules. In addition, a modified sinc kernel smoother is employed to develop a new method for calculating the dynamic contact angle. The accuracy and robustness of both the new and conventional methods were evaluated through calculations of physical properties, including apparent line tension and contact angle. The apparent line tension is positive only in highly hydrophilic systems and negative otherwise. Young's contact angle shows consistent results across methods for hydrophilic cases but varies in hydrophobic systems, signifying its dependence on the droplet shape analysis technique. The proposed methods yield physically meaningful results and are more reliable than standard circle fitting techniques.

---

1 Introduction

Droplets are everywhere in our world. They form the basis for utilizing aerosols and sprays,^1,2 their interactions with surfaces govern printing and coating technologies,³ and lab-on-a-chip devices.⁴ Classical free energy models and topographical analyses remain effective approaches to the study of wetting statics.^5,6 Advancing droplet-based technologies, however, requires a clear understanding of wetting dynamics.⁷ Techniques such as environmental scanning electron microscopy are essential for probing local effects and the role of the three-phase contact line,⁸ but their limited resolution and operational complexity pose significant challenges.⁹

Molecular dynamics (MD) simulations have become a powerful tool in studies of interactions of droplets with surfaces.^10–13 In order for MD to provide reliable data on wetting dynamics, robust procedures are necessary.¹⁴ For instance, the RMSD is currently used to provide information about the droplet state,^10,15 using assumptions that may or may not be valid. Additionally, because contact angle and contact radius calculations are based on averaging, their dynamic values are inaccessible by conventional methods.^16,17 We explore two major questions; (1) What are the appropriate methods to evaluate the droplet state during spreading? (2) How can the contact radius and the contact angle be dynamically measured during an MD simulation in a physically meaningful way to obtain information about the properties of the droplet?

On an ideal atomistically smooth, homogeneous, isotropic, rigid, and non-reactive surface, the equilibrium contact angle θ_eq is defined by⁸:

where τ_SLV and r_eq are the apparent line tension and the equilibrium contact radius. The first term on the right-hand side of eqn (1) is referred to as the Young¹⁸ equation, mathematically formulated by Dupré in 1869.¹⁹ The second term is denoted the Neumann-Boruvka term²⁰ and highlights the impact of line tension at the three-phase contact line on microscopic and sub-microscopic droplets.²¹ It is common to distinguish between high-wetting (0° < θ_eq), low-wetting (90° < θ_eq < 180°), and non-wetting (θ_eq = 180°) states. Fig. 1(a) illustrates the contact angle θ and contact radius r of a droplet.

Fig. 1 (a) A liquid droplet shaped like a spherical cap on a solid substrate. Parameters θ, r, and h are contact angle, contact radius, and contact height, respectively. (b) Initial state of the 3k system (for the notation regarding systems, please see Section 2). The Lennard-Jones wall is depicted in gray. Panels (c) and (d) show the final states of the 3k system when ε_OS in eqn (2) is equal to 0.1 and 1.65 kJ mol⁻¹, respectively.

Eqn (1) shows that surface tension effects dominate the thermodynamic aspects of wetting, determining the extent of surface wettability. When investigating wetting via MD, other effects must also be considered. First, gravity can affect the shape of large droplets.²² Capillary length κ⁻¹ presents the criterion beyond which gravitational effects become important. For a liquid of density ρ, where g is the gravitational acceleration. Because κ⁻¹ is typically on the order of millimeters, gravity is not considered in nanoscale MD simulations.

Another important factor is that the interactions at the perimeter of the droplet, i.e., at the three-phase contact line, have become a major factor at nanoscale.⁸ Schimmele et al.²³ have discussed the stiffness effects due to the contact angle and its curvature related to the line tension. Using a pseudo-two-dimensional system leads to a cylindrical droplet, and eliminates curvature effects at the contact line, making apparent line tension more accessible.²⁴ On the other hand, a spherical cap is expected to form in three-dimensional systems, where, due to curvature effects, the apparent line tension is difficult to calculate.²⁵ Therefore, the value and sign of the apparent line tension have been debated in experiments and using MD.²⁶ Scocchi et al.²⁷ have also shown that the droplet size has a direct impact on the value of τ_SLV even at this scale. They observed two regimes: one for droplets smaller that ∼2000 molecules and another one for larger droplets. The value of τ_SLV is generally expected to be around 20 pN for water.²⁸ However, Zhang et al.²⁹ showed that it is scattered and spans a wide range in the literature. They also presented a modified form for the line tension based on the work of Schimmele et al.²³ Although atomistic simulations provide the highest resolution for droplet shape analysis, no single method is universally optimal at this scale. A detailed analysis of the measurement techniques, particularly for three-dimensional systems, can help resolve such controversies.

An additional effect at the nanoscale is linked to diffusion and the amount of energy accessible to the droplet. Studies have demonstrated that the self-diffusivity of a droplet exhibits a power-law increase as the size decreases.³⁰ A high self-diffusivity coefficient causes the system to attain equilibrium more rapidly. Hence, droplet mobility is critical in examining its spreading limit at small scales.

Here, we establish reliable procedures for the study of such systems¹⁴ and evaluate various conventional methods. We consider a range of wettabilities and droplet sizes in partial and complete wetting regimes using MD simulations. Typically, measurements such as total energy and a measurable parameter's root mean square deviation (RMSD) establish the droplet state. Total energy fails to offer informative data about the wetting state or the spreading limit, and RMSD can also be misleading, considering droplet mobility. We demonstrate that only one component of the RMSD is sufficient for studying this problem. A new radius of gyration analysis is introduced to monitor the evolution of the contact radius, followed by a new technique for dynamic contact angle calculation. These methods, along with established techniques, are compared with the available literature data and evaluated for their accuracy and consistency in estimating droplet physical properties.

2 Computational methods

GROMACS version 2022.5³¹ was used for the atomistic MD simulations, and molecular visualizations were done with visual molecular dynamics (VMD).³² The optimized potential for liquid simulations all-atom³³ force field was employed for the simulations. All simulations were performed in the canonical ensemble with a time step of 1 fs.

In this study, the standard system was a non-volatile liquid droplet on a substrate. The droplet's initial configuration was a sphere generated by Packmol³⁴ with a density of 1.0 g cm⁻³. It consisted of simple point charge (SPC)³⁵ water molecules, ranging from 3000 to 12 000, placed in a 40 × 40 × 20 nm³ box. Periodic boundary conditions were applied only along the x- and y-directions. The Nosé–Hoover thermostat^36,37 kept the temperature of water molecules at 300 K with a damping parameter of 0.1 ps. The van der Waals and the real part of the electrostatic interactions were cut off at 9.0 Å.³⁵ The long-range Coulombic interactions were calculated using the particle–particle particle–mesh (P³M) solver.³⁸ To allow for the use of P³M in a system that is periodic only along the x- and y-directions, the procedure introduced by Yeh and Berkowitz³⁹ was followed.

The ideal substrate was modelled as a structureless, flat wall that interacts via the Lennard-Jones (LJ) potential with the water. When the 12-6 LJ potential is integrated over the wall elements, the 9-3 potential for surface-water interaction is obtained. It must be noted that given the definition of the SPC water model, the water–substrate interactions only included oxygen atoms of water molecules, where^40,41

In eqn (2), ε_OS and σ_OS define the depth of the potential well and the distance of zero interaction potential for the oxygen–substrate interactions, respectively. ρ_S is the substrate number density, assumed to be unity for the simulated LJ wall. z_OS,i is the oxygen–substrate distance along the z-direction. In addition to the 9-3 LJ wall at z = 0, an additional repulsive 9-3 LJ wall at z = 20 nm acts as the mirror boundary condition for the closed system of study. The simulations included three main steps: equilibration, deposition, and spreading. During equilibration and deposition, the substrate was set to behave like a repulsive LJ-wall with a small σ_OS value in eqn (2). This ensured that no spreading occurred during the first stages of the simulations. The initial velocities were drawn from a Maxwell–Boltzmann distribution. Equilibration was performed for 20 ns for the initial droplet to assume a spherical form. The final configuration of the droplet was used as the initial state of the deposition runs. During the deposition, a constant force of 1 kJ mol⁻¹ nm⁻² was applied to the center of mass of the droplet along the z-direction toward the LJ wall. The simulations were stopped when the droplet came within 4.0 Å from the substrate. At this stage, the velocity of the droplet's center of mass was removed, the interaction parameter ε_OS in eqn (2) was adjusted, and the spreading step started. For all the spreading simulations, σ_OS in eqn (2) was set equal to 3.0 Å. Fig. 1(b) illustrates the state of a system with 3000 water molecules at the end of the deposition step. Droplet spreading was performed for 4 ns. The simulation conditions, including the thermostat, were the same during all steps.

After several preliminary simulations, eight values for the parameter ε_OS in eqn (2) were used to simulate a range of wettabilities, 0.1, 0.3, 0.5, 0.7, 1.2, 1.45, 1.65, and 2.2 kJ mol⁻¹. The value ε_OS = 0.1 is the most hydrophobic system, and ε_OS = 2.2 represents a completely wetting system. Each simulation was repeated three times from equilibration to spreading, and the averaged data was analyzed. For compact notation, the number of water molecules given in thousands is used, which is followed by the letter ‘k’ and the ε_OS value. For example, the system 3k0.1 refers to 3000 water molecules with the highest hydrophobicity, whereas 12k2.2 refers to 12 000 water molecules with the highest hydrophilicity. Fig. 1(c) and (d) depict the final states for 3k simulations with ε_OS equal to 0.1 and 1.65 kJ mol⁻¹, respectively.

In addition to the main simulations, a 5 × 5 × 5 nm³ box of 7000 water molecules was placed at the center of a 5 × 5 × 17.5 nm³ simulation box. A 10 ns equilibration and 4 ns production run were performed in the canonical ensemble to extract liquid properties. Simulations were repeated five times, and other conditions were the same as the other simulation steps for the droplet. The density, surface tension, and shear viscosity of the SPC water model were estimated to be ρ_L = 973.53 ± 3.99 kg m⁻³, γ_LV = 47.42 ± 0.16 mJ m⁻², and η = 0.41 ± 0.01 cP, respectively, showing a good agreement with the literature values.^42,43

3 Results and discussion

First, a cluster analysis was carried out on the oxygen atoms of water in order to distinguish between the molecules in the droplet and vapor. If the distance between two oxygen atoms was ≤3.4 Å, the molecules were considered to be in the same cluster. This value was chosen based on the first hydration shell of the SPC water model.³⁵

3.1 Root mean square displacement (RMSD) analysis

RMSD is calculated usingwhere N is the number of atoms in the system, and Δs_i(t) is the displacement vector of atom i. It is defined as (Δx_i, Δy_i, Δz_i), where Δx_i = x_i(t) − x_i(t₀) is the x-component of the displacement of the atom i from the reference time t₀ to time t. Changes in RMSD indicate alterations in the system's configuration and are expected to show asymptotic behavior toward equilibrium.

The RMSD for one series of simulations of both the smallest and the largest systems are shown in Fig. 2. It does not fully plateau even after the initial increase and continues to change. RMSD keeps fluctuating as the simulation continues in small systems. In larger systems, a general trend is observed where RMSD increases with ε_OS as defined in eqn (2), which provides a measure for surface wettability.

Fig. 2 Time evolution of RMSD for the systems with (a) 3000 (3k) and (b) 12 000 (12k) water molecules calculated using eqn (3). The results are shown for ε_OS range of 0.1 to 1.65 kJ mol⁻¹ as defined in eqn (2) for one series of simulations. The greater fluctuations of the RMSD value for the smaller 3k systems illustrate a different behavior to that of larger 12k systems, where a general trend is observed. It must be noted that RMSD does not fully plateau in the systems regardless of the droplet size.

3.1.1 Contributions of the different RMSD terms. The main drawback of using RMSD in finite systems, such as droplets, is that it does not accurately represent the system's dynamics, and here it fails to reveal any information about the droplet's state. Over-interpretation can lead to associating multiple conformations to the droplet, or to the conclusion that equilibrium is a stationary state without any translational motion. Different approaches have been proposed for overcoming this obstacle, espcially in the context of proteins, for example, by weighted RMSD⁴⁴ or root-mean-square fluctuations.^45,46

Here, we divide RMSD into three contributions: internal motion of the molecules inside (RMSD_INT) the droplet with respect to the center of mass (COM) of the droplet, translational motion of the COM of the droplet (RMSD_COM), and a coupling term (RMSD_CPL). The term RMSD_INT is close to the concept of deviation and self-diffusion, whereas RMSD_COM is due to the diffusion of the droplet as a whole. The derivation is provided in Appendix A. The above RMSD contributions can be written as

and

RMSD²_COM(t) ≡ ‖ΔS(t)‖². (6)

In eqn (4)–(6), ΔS(t) is the displacement vector of the COM at time t with respect to the reference time t₀. The tilde signifies displacement relative to the COM, and s̃_i(t) is the deviation vector of particle i from the COM at time t.

Fig. 3 shows the time evolutions of RMSD_INT and RMSD_COM for one series of 3k and 12k systems. RMSD_CPL is not shown since it has a negligible contribution to RMSD due to the impact of droplet's axisymmetric shape on term. RMSD_INT has always a slight positive slope since the molecules in the droplet never cease to move. The mean squared displacement values for the same systems are illustrated with lower opacity in Fig. 3(a) and (b). Both graphs evolve similarly, supporting the notion that RMSD_INT is close in concept to self-diffusion. Comparing the panels (a) and (b) in Fig. 3, with panels (c) and (d) clarifies that the source of RMSD_INT deviation from the mean squared displacement is due to COM contribution. Expanding eqn (4) (or A5 in Appendix A) leads to the same conclusion. Therefore, the internal term is the main reason why the RMSD never completely plateaus for similar systems.⁴⁷ However, as the droplet size becomes larger, its mobility on the surface decreases, and the RMSD_COM plateaus. Because this term captures the hindrance in the motion of the droplet as a single entity, it is expected to depend on the droplet size. Ma et al.³⁰ have shown that droplet mobility decreases with its size, following a power-law relation. This behavior is reflected in the final stage of the simulation for smaller droplets as shown in Fig. 3(c). The full decomposition of the RMSD_COM (ΔS) in x-, y-, and z-directions for the same 3000 and 12 000 molecule systems are also presented in Fig. 10 and 11 of Appendix A, respectively. This result also raises the question of whether clusters with a smaller number of molecules can even be considered droplets for the study of wettability and spreading.

Fig. 3 Time evolutions of the RMSD components for the systems with 3000 (3k), left, and 12 000 (12k), right, water molecules. The results are shown for ε_OS range of 0.1 to 1.65 kJ mol⁻¹. Panels (a) and (b): the time evolution of RMSD_INT², calculated using eqn (4) in solid colors. The mean squared displacement values for the same systems are illustrated with lower opacity. The values are shifted for clarity. RMSD_INT is close to the concept of self-diffusion and never reaches a constant value. The deviation from mean square displacement is due to the center of mass contribution, see Appendix A. Panels (c) and (d): RMSD_COM, calculated using eqn (6). RMSD_CPL (eqn (5)) is not shown as its value is close to zero. Because the droplet mobility as a whole decreases on the surface, RMSD_COM plateaus for larger systems, see Fig. 10 and 11 of Appendix A for the full decomposition of RMSD_COM along the x-, y-, and z-directions.

3.1.2 Center-of-mass analysis. Since the systems are axisymmetric, the contributions of the lateral components X_COM(t) and Y_COM(t) are negligible; the capital letters refer to the COM component. On the other hand, Z_COM(t) can provide information about the spreading limit of the droplet. Hautman and Klein⁴⁸ have shown that there is a relationship between the average Z_COM value and the contact angle. Since the contact angle is a state function, Z_COM can be defined⁴⁸ aswhere 〈Z_COM〉 is the time average of the Z-coordinate of the droplet's COM relative to the solid surface. R₀ is the radius of the spherical cap fitted to the droplet and approximated as (3N/4πρ₀)^1/3, where N is the number of molecules in the droplet, and ρ₀ is the number density of the bulk liquid. Eqn (7) has been used to calculate the contact angle in other studies,^49,50 and it has been shown that Z is a key quantity in determining the state of the droplet.

ΔZ_COM(t) provides an accessible measure of the spreading limit. Fig. 4(a) and (b) display the average behaviors of ΔZ_COM(t) for the 3k and 12k systems, respectively. Due to the greater spreading on the hydrophilic substrate, ΔZ_COM(t) increases with the droplet size. Furthermore, greater fluctuations are seen for the hydrophobic systems as compared to the hydrophilic ones. This can be explained by considering the strength of interactions. Because LJ interactions are not long-ranged, water molecules further away from the substrate are unaffected by the solid–liquid interactions. This leads to cohesive interactions such as ε_OO and the Coulombic potential playing a more dominant role in establishing the droplet's shape. Consequently, it provides a better method than the general RMSD for analyzing the state of the droplet.

Fig. 4 The evolution of ΔZ_COM(t) averaged over three simulation runs for the systems with (a) 3000 (3k) and (b) 12 000 (12k) water molecules. The results are shown for ε_OS range of 0.1 to 1.65 kJ mol⁻¹. The shaded areas represent the standard deviation. ΔZ_COM(t) is closely related to 〈Z_COM〉 in eqn (7), and provides a general measure of the equilibrium state of the droplet.

3.2 Number density-based droplet geometry

In atomic scale simulations, the inherent resolution of the system allows for identification of interfaces and their effects. When a droplet is resting on a surface, the interfacial layer can be identified using probability density profiles ρ(z) of the water molecules along the z-direction.⁵¹ To calculate ρ(z), the simulation box is divided into elements with a height of dz = 0.1 Å. The number of oxygen atoms representing water molecules is counted in each region and normalized by the region's volume.

The average values of ρ(z) for the final 1 ns of the 12k simulations are shown in Fig. 5(a). The plots do not start at zero because σ_OS > 0 as defined in eqn (2). This value approximately defines the minimum distance between the water molecules and the substrate. Almost all the droplets, with the exception of the systems with the lowest ε_OS as defined in eqn (2), show a clear first peak and a trough. A weak dependence on wettability is observed for the trough positions in the inset of Fig. 5(a). For lower ε_OS, the value is at 4.5 Å and for higher ε_OS, it shifts to 4.4 Å. Zhang et al.²⁹ showed that the value chosen for the height of such an interface can greatly affect the interpretation of the results. They chose values ranging from 1–2σ_OS. Our results confirm that the width of the interfacial layer is ∼1.5σ_OS. It is also observed that not all systems display three peaks, this is particularly the case for the more hydrophobic ones with low ε_OS. Furthermore, a depleted interfacial layer can be identified for the systems with the lowest wettability. Fig. 5(b) provides a schematic representation of the droplet.

Fig. 5 (a) Probability density of water molecules along the z-direction for systems with 12 000 (12k) water molecules. The results are shown for ε_OS range of 0.1 to 1.65 kJ mol⁻¹. The inset shows the probability density up to z = 15 Å. A slight shift in the positions of the troughs is observed, signifying an increase in the width of the interfacial layer with a decrease in ε_OS. (b) The droplet model based on the probability density profiles.

3.3 Radius of gyration (R_g²) analysis for contact radius

The spreading limit presents the point at which the droplet reaches its equilibrium. Using the obtained values for the height of the interfacial layer, a new model for contact radius is proposed. In its standard form, R_g² provides a measure of particle distribution about a reference axis. In this model, the interfacial layer is considered to be a very thin cylindrical disc with its axis of symmetry providing the reference. The radius of this disc, r, is defined as the contact radius, see Fig. 5(b). Based on the relationship between the disc radius and its radius of gyration, the contact radius can be calculated and monitored throughout the simulation. The derivation of this model is presented in Appendix B, and the final equation is written aswhere R_g,IL is defined for the interfacial layer. To assess the model's performance, the radial probability distribution of the water molecules in the interfacial layer was used to estimate r_eq. The procedure involves counting water molecules within circular shells of width 0.2 Å, centered at the interfacial layer COM. It is then normalized by the volume of the shells. A hyperbolic tangent fit to the radial probability distribution average gives the r_eq value.⁵²

Fig. 6 shows four snapshots of the interfacial water molecules, represented by the oxygen atoms, for 3k1.65 and 12k1.65 systems, respectively. The radial probability distribution estimate of r_eq for the final 1 ns of the simulations is depicted using black circles. It is evident that as the droplet reaches its spreading limit, the value (red, dashed circle) approaches r_eq (black circle). While the radial probability distribution function is reliable for calculating the equilibrium contact radius, it requires averaging over an interval. However, this method allows for direct monitoring of the contact radius during the spreading process.

Fig. 6 Top-view images of the interfacial layer for 3k (top row) and 12k (bottom row) systems with ε_OS = 1.65 kJ mol⁻¹ at different time steps from a single simulation series. Blue dots indicate oxygen atom positions, representing water molecules. The red circle denotes the contact area, with its radius computed using eqn (8). The black circle, derived from the radial probability distribution of water molecules in the interfacial layer, has a radius r_eq. As the droplet reaches its spreading limit, the two circles align with high accuracy.

3.4 Droplet profile and contact angle

The contact angle is a fundamental thermodynamic property. It sets the criterion for the equilibrium of a droplet resting on a surface and can be viewed as the wetting potential.⁸ There are many types of contact angles, but in this text, the definitions presented in ref. 8, 21 and 53 are adopted. Therefore, θ_eq is considered accessible via sampling the simulation trajectory at equilibrium, θ_Y is the Young's contact angle on an ideal surface, and the dynamic contact angle is the contact angle of a spreading droplet as it relaxes toward θ_eq.

Computing the contact angle for nano-droplets has a range of complications,²⁵ mainly due to the shape fluctuations observed at this scale. The conventional method uses the mass-density distribution of molecules^17,54–56 with modifications to this method proposed recently.⁵⁷ The mass density based methods require averaging, making the calculation of the dynamic contact angle challenging. Instead, one can perform shape analysis on the liquid–vapor interface that is made instantaneously accessible by methods such as Generalized Identification of Truly Interfacial Molecules (GITIM)^58–60 or a Euclidean 3D convex hull of the density profile.⁶¹ GITIM uses the concept of probe spheres (α-shapes^60,62,63) to identify planar and non-planar interfaces, such as external surfaces and internal voids. Once the liquid–vapor interface is extracted, circular or elliptical fitting procedures can be applied to calculate the contact angle.⁶⁴ Other methods use coarse-graining concepts, and equivalent area and volume to fit a spherical cap to the droplet molecules.⁶⁵ Such methods assume a general ellipsoidal shape for the droplet and are applicable at equilibrium when the spreading limit is reached. However, they fail during the early stages of spreading and droplet motion, especially in the presence of a precursor film. Additionally, the dimensionality of the system can impose further challenges as cylindrical droplets lead to more stable profiles making them easier to analyze than three-dimensional ones.

To overcome the above limitations, a new procedure for dynamic measurement of the contact angle is proposed. GITIM is utilized to gain direct access to the interfacial molecules. A probe radius equal to σ_OO for the SPC water model was defined to find the water molecules that reside at the liquid–vapor interface. Our main assumption is that the droplet has an axisymmetric profile, and the original interfacial points are rotated about the droplet's axis of symmetry. Fig. 7(a) illustrates the original interfacial points and the rotation procedure about the droplet's axis of symmetry that produces a 2D profile.

Fig. 7 (a) The procedure to rotate the positions of the GITIM^58,59 interfacial molecules to extract the droplet profile on the xz-plane. Panels (b) and (c) illustrate the xz-plane scatter plot of the original liquid–vapor interfacial points (gray) and rotated profile points (blue) for the system with 12 000 water molecules and ε_OS = 1.2 kJ mol⁻¹ at t = 2 ns. The black solid and dashed lines represent the surface and the solid–liquid interfacial layer, respectively. In (b), the orange and green lines illustrate the extrapolated profile to the surface and the first density peak, respectively, and the red dots show the circle fitted to the rotated profile points. In (c), the red dots show the left and right side of the rotated profile that was separately smoothed by the modified sinc smoother.⁶⁶

Fig. 7(b) shows the xz-plane profile of the interfacial molecules before (gray) and after rotation (blue), as well as the fitted circle (red) for one series of 12k1.2 simulations after droplet equilibration at t = 2.0 ns. The orange lines represent the calculated slope for the contact angle and the contact radius when the profile is extrapolated to the surface, as is common in mass density-based methods. The green line depicts the same values for the profile extrapolated to the interfacial layer's density peak, see also Fig. 5(a). The radius of gyration method estimates the value of the contact radius to be 66.4 Å for this snapshot, which is similar to the value calculated by extrapolating to the first peak of the interfacial layer.

For very low contact angles, extrapolation overestimates the contact radius, whereas underestimation of the contact radius is expected as contact angles approach 180°. Consequently, extrapolation to the first peak should be adopted, as it yields physically meaningful contact radius values at this scale. The mass density-based method estimates the contact angle to be 66.0°. As shown in Fig. 7(b), the same value by extrapolating to the surface and to the first peak of the interfacial layer is 66.4° and 67.7°, respectively. Therefore, the contact angle is not significantly affected by the extrapolation process.

The rotation procedure results in an accurate contact radius and contact angle for a droplet at equilibrium when it is combined with circle fitting. To extend the method to dynamic contact angle calculations, the rotated two dimensional profile is first divided into left and right segments. The profile is then expressed in polar coordinates. Radial sampling is carried out relative to the local curvature, which provides consistent sampling on both hydrophilic and hydrophobic substrates. The resulting points are linearly interpolated to generate a uniform set of profile positions. The sampled radii are then smoothed using the modified sinc filter of Schmid et al.⁶⁶ For evaluation of wetting parameters, including contact angle and contact radius, the smoothed data are mapped back to Cartesian coordinates. A local linear regression near the three phase point of each two dimensional profile yields the contact angle. The analysis requires four main parameters. The first is the number of sampling points, which is set by the largest characteristic dimension of the droplet. For hydrophilic substrates, this dimension is the equilibrium contact radius, whereas for hydrophobic substrates it is the droplet height. In practice, this is determined from the maximum extent of the atomic coordinates. The results show that one to two sampling points per unit length are sufficient to achieve uniformity in the extracted profile for shape analysis. The method remains flexible and can be adjusted to the specific characteristics of the simulated droplet. The remaining three parameters are obtained from a number density analysis, see also Fig. 5(a). They include the interfacial layer width and the positions of the first and second peaks along the z-direction. The first peak sets the lower bound for extrapolating the profile, and the second peak defines the region used for local linear regression near the three phase point. When a well defined interfacial region cannot be identified, an interfacial width of approximately ∼1.5σ_OS is recommended based on the present results and the work of Zhang et al.²⁹ In such cases, the first and the second peaks are set to one half and twice this value. Therefore, an advantage of the method is the ability to specify the peak locations to analyze local contact angles in the presence of a precursor film. Additionally, omitting the rotation step allows the method to be applied directly to pseudo-two-dimensional cylindrical droplets. In general, because the procedure is carried out independently on the left and right sides, it also enables analysis of advancing and receding fronts for contact angle hysteresis. Fig. 7(c) shows the instantaneous contact radius and contact angle calculated for the left and right profiles. The average contact radius and contact angle are 66.5 Å and 67.0, respectively. These results are in good agreement with the mass density-based and circle fitting methods.

3.5 Comparative analysis of droplet characterization methods

In this study, the radial probability distribution of the interfacial layer⁵² and the mass density distribution⁵⁴ are employed as conventional methods for determining the equilibrium contact radius and contact angle, respectively. Fig. 8(a) presents contact radius data obtained from the proposed methods for three interaction parameters in one instance of the 12 000-molecule system. Black dashed lines denote the radial probability distribution of the interfacial layer. For clarity, only one data point per five picoseconds is plotted for the circle fitting and local smoothing methods.

Fig. 8 (a) Comparison of radius of gyration, local smoothing, and circle fitting methods in estimating the contact radius of one series of the 12 000 molecule simulation. The dashed lines represent the radial probability distribution of the interfacial layer molecules used as the conventional procedure.⁵² (b) Comparison of local smoothing and circle fitting methods in calculating the contact angle from one series of the 12 000-molecule droplet. The dashed lines represent that conventional mass density⁵⁴ procedure for calculating the contact angle. The circle fitting and local smoothing data are shown in 5 ps intervals for better visualization. (c) The root mean square error (RMSE) of the extracted profiles for one series of the 12 000-molecule system.

The radius of gyration method shows excellent agreement with the reference approach⁵² at equilibrium. However, both circle fitting and local smoothing display slight overestimation and underestimation as described. The R_g,IL leads to consistent, accurate calculation of the contact radius. The local smoothing method evolves similar to R_g,IL while slightly deviating from the exact values. In contrast, the circle fitting method fails to produce contact radius data for t < 50 ps, resulting in divergent behavior during early spreading. R_g,IL is the recommended method for the study of spreading dynamics, the local smoothing method is the second best choice, and the circle fitting methods should not be employed for calculating contact radius. Contact angle measurement is inherently more complex, as it mainly relies on shape analysis. Fig. 8(b) compares results from the same three systems as in panel (a). In this case, the black dashed lines indicate contact angles computed using the conventional mass density distribution. While all methods seem to converge toward similar equilibrium values, the circle fitting method fails to detect droplet contact with the surface for t < 50 ps.

To compare the accuracy of the circle fitting and local smoothing techniques, root mean square error (RMSE) values were computed from the extracted interfacial profiles. Fig. 8(c) presents RMSE results for a 12 000-molecule system. Across the studied range of wettabilities, the local smoothing method consistently produces lower RMSE values compared to circle fitting. Although the errors are comparable in the complete wetting regime where ε_SO = 2.2 kJ mol⁻¹, local smoothing reduces the RMSE by approximately 50% on hydrophobic substrates. Similar trends were observed for other droplet sizes.

As recently confirmed by Polovinkin et al.,⁶⁹ contact angle measurements are strongly influenced by the chosen analysis method. Comparative evaluation of techniques for estimating droplet's physical properties helps to identify their respective limitations. Eqn (1) enables the calculation of apparent line tension (τ_SLV) and Young's contact angle (θ_Y) by linear regression of cos(θ_eq) versus 1/r_eq. Fig. 9(a) illustrates the computed τ_SLV values. For this fitting procedure, r_eq was calculated using the radius of gyration method. It must be emphasized that employing r_eq values obtained from local smoothing or circle fitting yields qualitatively similar results. Furthermore, the calculated values based on cylindrical slicing⁶⁷ is also displayed for comparison. Because cylindrical slicing cannot provide a complete profile of the droplet at very low contact angles, τ_SLV and θ_Y could not be calculated for the system with ε_SO = 2.2 kJ mol⁻¹ using this method.

Fig. 9 The calculated (a) apparent line tension τ_SLV and (b) Young's contact angle using mass density, local smoothing, cylindrical slicing, and circle fitting procedures. All methods except local smoothing involve fitting a circle to the droplet profile extracted from different sources. The mass-density method uses the spatial distribution of water molecules within the simulation box to find the profile.⁵⁴ The cylindrical-slicing method estimates the profile by sectioning the droplet parallel to the substrate and fitting a hyperbolic tangent to the resulting radial probability distribution.⁶⁷ The circle-fitting procedure obtains the profile using GITIM^59,60 and rotates the profile points before applying a circle to them. The local method smooths⁶⁶ the rotated profile and performs a linear regression near the three-phase contact line, extending up to the second layer of the number-density profile (see Fig. 5). The data is based on linear fits to eqn (1). The dashed line represents the calculated values for a cylindrical droplet of dispersive fluid confined between two dispersive walls.⁶⁸ The gray area shows the apparent line tension for a three dimensional droplet on aliphatic chains with hydroxyl head groups.²⁵ The cylindrical slicing method could not compute the contact angle for ε_SO = 2.2 kJ mol⁻¹.

The local smoothing method yields negative τ_SLV values across most contact angles, except at low angles. The mass density results show similar qualitative behavior. In contrast, the circle fitting and cylindrical slicing methods give positive τ_SLV values for θ_Y > 90°. Fig. 9(b) shows that local smoothing predicts higher contact angles at lower solid–liquid interaction strengths. For θ_Y < 90°, all four methods show good agreement.

To further evaluate the accuracy of these methods, reference data from the literature are included in Fig. 9(a). The dashed line corresponds to a cylindrical droplet of dispersive fluid confined between two dispersive walls.⁶⁸ The shaded region denotes the range of τ_SLV values reported for a three-dimensional droplet on hexagonally anchored aliphatic chains terminated with hydroxyl groups.²⁵ Both studies report negative τ_SLV for hydrophobic substrates, with positive values appearing only when θ_Y → 0°.

In addition to droplet geometry and atomistic versus dispersive substrate characteristics, the choice of force field and water model affects the calculated line tension. This contributes to the variation in reported τ_SLV values. For example, Kanduč et al.²⁵ employed the united atom GROMOS force field⁷⁰ and omitted the Yeh–Berkowitz correction³⁹ for long-range electrostatics under two-dimensional periodic boundary conditions. Additionally, Kanduč⁷¹ studied the same system in an earlier work with the results suggesting that water molecules partially penetrate the aliphatic chains, leading to a slight convex-shaped contact area toward the substrate. Furthermore, droplets containing fewer than 2000 molecules have been shown to exhibit line tension behavior that deviates from the expected regime.²⁷ As discussed in Section 3.1.1, treating such small clusters as droplets in wetting studies remains a subject of debate. It is also important to note that the presented data should be interpreted as the simulation of a nano-droplet on an ideal surface as opposed to an atomistic one. These findings suggest that the wide variability in literature values for τ_SLV²⁹ arises from inconsistencies in shape analysis methods. Namely, circle fitting cannot capture the local bending near the contact line at the limit of molecular scale of MD simulations. This comparative evaluation indicates that local smoothing offers a more reliable and physically meaningful result than conventional methods.

4. Conclusion

In this study, we have presented a comprehensive investigation of conventional and new techniques for measuring dynamic properties of nano-droplets during spreading. MD simulations were conducted and repeated at least three times for a range of SPC water model droplet sizes spreading on an LJ wall with varying wettabilities.

RMSD was shown to be unsuitable for measuring system dynamics because its center of mass component plateaus only for large systems, while its internal component reflects the diffusion of molecules within the droplet. Instead, the z-component of the center of mass is recommended.

New methodologies were introduced to address limitations of conventional methods, including interval averaging and incorrect shape assumptions. A new equation based on the radius of gyration of solid–liquid interfacial molecules is proposed to compute the instantaneous contact radius. The GITIM algorithm^58,59 is employed to identify the liquid–vapor interface. Combined with a modified sinc smoother,⁶⁶ a new procedure is developed to extract a locally smooth droplet profile and calculate the dynamic contact angle.

Comparative analysis shows that both new methods provide reliable results that are in agreement with literature values. The apparent line tension is negative for hydrophobic substrates and positive only in highly hydrophilic systems. Additionally, Young's contact angle remains consistent across methods, except for strongly hydrophobic cases. The physically consistent results confirm the robustness of the methods and their suitability for future studies of dynamic wetting phenomena in MD simulations.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data is available at https://github.com/SoftSimu/droplets.

Appendices

A RMSD in terms of different contributing motions

The RMSD can be given aswith s_i(t) ≡ x(t)î + y(t)ĵ + z(t)k̂, and Δs_i(t) ≡ s_i(t) − s_i(t₀) = [x_i(t) − x_i(t₀)]î + [y_i(t) − y_i(t₀)]ĵ + [z_i(t) − z_i(t₀)]k̂. The number of molecules is denoted by N.

The RMSD above can be decomposed into two types of motion: (1) the motion of the whole droplet described by the center of mass (COM) position vector S(t), and (2) the motion of the molecules with respect to the COM position, s̃_i(t). Thus, we can use the following definition:

s̃_i(t) ≡ s_i(t) − S(t) (A2)

The displacements can then be given as as Δs̃_i(t) = s̃_i(t) − s̃_i(t₀) and ΔS(t) = S(t) − S(t₀). Thus, Δs_i(t) in eqn (A1) becomes and eqn (A1) can be re-written asand finally, separated intoIn this stage, we delineate various types of RMSD as

RMSD²_COM(t) ≡ ‖ΔS(t)‖². (A7)

With the above, eqn (4) transforms into

RMSD² = RMSD²_INT + 2 RMSD_CPL + RMSD²_COM. (A8)

eqn (A1) is the RMSD with respect to a general coordinate, whereas eqn (A4) is the RMSD of an N-particle system relative to the COM. Eqn (A5) and (A7) are in a square form, but the right-hand side of eqn (A6) can be negative. Eqn (A7) can be interpreted as the RMSD of a hypothetical single-particle system located at the COM coordinates, that is, RMSD_COM² = 6D_COMt, representing the translational diffusion of a hypothetical particle located at the COM. Similarly, RMSD_INT² = 6D_INTt and RMSD_CPL² = 6D_CPLt. D_INT describes the relative motion of the N-particle system with respect to its COM, whereas D_CPL quantifies the dynamic correlation between internal and COM motions. These equations connect the RMSD contributions to their mean squared displacement counterpart, predicting that the diffusion coefficient may also be decomposed into D = D_INT+ 2D_CPL + D_COM. It is important to notice the decomposition requires that the COM is defined in a consistent manner, the number of particles remains constant, and that there are no external constraints regarding the COM. Physically, the coupling term describes the coupling between the motion of particle i, and the overall translation of the system, i.e., how much that particle's displacement is aligned with collective shift of the system. In equilibrium systems, the coupling term vanishes.

As an example, the calculation of each component of RMSD_COM is demonstrated using the same simulation data as in Fig. 3(c) and (d). Fig. 10(a)–(g) display RMSD_COM in black for each wettability case (different ε_SO) in a 3000 molecule system. The quantities ΔX_COM, ΔY_COM, and ΔZ_COM appear as dashed, dotted, and solid colored curves. Fig. 10(h) presents the corresponding RMSD_COM data shown in Fig. 3(c) on a linear scale. Fig. 11 provides the same analysis for one series of the 12 000 molecule system, corresponding to Fig. 3(d).

Fig. 10 (a)–(g) RMSD_COM for each value of ε_SO in black. ΔX_COM, ΔY_COM, and ΔZ_COM are shown as dashed, dotted, and solid colored lines, respectively. (h) RMSD_COM for the different ε_SO value, same as Fig. 3(c). All results are for one series of a 3000-molecule droplet. Only one data point per 20 picoseconds is shown for clarity.

Fig. 11 (a)–(g) RMSD_COM for each value of ε_SO in black. ΔX_COM, ΔY_COM, and ΔZ_COM are shown as dashed, dotted, and solid colored lines, respectively. (h) RMSD_COM for the different ε_SO value, same as Fig. 3(c). All results are for one series of a 12 000-molecule droplet. Only one data point per 20 picoseconds is shown for clarity.

The motion of a nanodroplet in the xy-plane was initially attributed to Brownian motion of a particle cluster on a smooth surface⁴⁷ and has since been observed on atomistic surfaces as well.⁷² These studies have shown that both surface wettability and droplet size affect the magnitude of lateral motion. Specifically, random lateral displacement is more pronounced on hydrophobic surfaces and for smaller droplets.⁴⁷ Furthermore, it has been shown that the lateral motion can be reinforced to suppress pinning effects. For instance, Foroutan et al.⁵¹ showed that droplet translation can offset the influence of surface heterogeneities, causing the nanodroplet to behave as if on an ideal surface. Recent work has confirmed this observation and linked the behavior to low contact angle hysteresis on self-assembled monolayers.¹³ Ma et al.³⁰ also examined how the diffusion coefficient of a water droplet depends on the amplitude of propagating ripples on a graphene sheet. A full dynamical analysis of droplet motion lies outside the scope of this work. However, for droplet diffusion on surfaces beyond the studies in ref. 13, 30, 51 and 72–74, the RMSD decomposition provides a theoretical basis for further investigation.

B Radius of gyration and contact radius

Radius of gyration analysis is generally used to measure the spatial properties of clusters. The elements of the gyration tensor (G) can be written as⁷⁵where x and y determine the directions. p_ix is the position of particle i along the x-direction, and p̄_x denotes the mean position along x-direction. The trace of G is related to the radius of gyration (R_g²) viawherein p_i and p̄ are the position and mean position vectors, respectively. If p̄ is defined as the COM vector, then G would describe the mass distribution of the particles around it. To find a definition for the radius of the disc-like interfacial layer, we invoke the equation of the moment of inertia (I) for a thin, uniform disc⁷⁶,where M is the mass, and r is the radius of the disc. In a particle system, I and M can be given asandwhere p̄ is the position vector of the axis of rotation. p_i and m_i are the position vector and mass of the particle i, respectively. Here, we consider the axis of rotation passing through the center-of-mass of the interfacial layer. With the above,for a system of identical particles that reside in the interfacial layer,Therefore, r can be defined in terms of R_g,IL as

Acknowledgements

F. E. thanks the Flight PS752 Commemorative Scholarship Program. M. K. thanks the Natural Sciences and Engineering Research Council of Canada (NSERC), the Canada Research Chairs Program, and Foundation PS. Digital Research Alliance of Canada provided the computing facilities. F. E. thanks Alireza Fazeli for his invaluable input in analyzing the data. This work was supported by the Research Council of Finland (Flagship of Advanced Mathematics for Sensing Imaging and Modelling grant 358944).

Notes and references

1. M. F. Ruiz-Lopez, J. S. Francisco, M. T. Martins-Costa and J. M. Anglada, Nat. Rev. Chem., 2020, 4, 459–475.
2. N. S. Ha, M. de Raad, L. Z. Han, A. Golini, C. J. Petzold and T. R. Northen, RSC Chem. Biol., 2021, 2, 1331–1351.
3. Y. Jiang and C.-H. Choi, Adv. Mater. Interfaces, 2021, 8, 2001205.
4. Y. Wang, Z. Chen, F. Bian, L. Shang, K. Zhu and Y. Zhao, Expert Opin. Drug Discovery, 2020, 15, 969–979.
5. Y. Howell, R. Blaikie and S. Lowrey, J. Colloid Interface Sci., 2025, 677, 1014–1021.
6. K. M. Al Balushi, K. Sefiane and D. Orejon, J. Colloid Interface Sci., 2022, 612, 792–805.
7. P. G. de Gennes, Rev. Mod. Phys., 1985, 57, 827–863.
8. E. Bormashenko, Adv. Colloid Interface Sci., 2015, 222, 92–103.
9. D. B. Peckys, E. Macías-Sánchez and N. de Jonge, MRS Bull., 2020, 45, 754–760.
10. M. Foroutan, F. Esmaeilian and M. T. Rad, Phys. Fluids, 2021, 33, 032017.
11. C. Xie, J. Shi, Y. Luo, G. Chu and H. Li, Langmuir, 2022, 38, 9822–9832.
12. D. Niu, H. Gao, G. Tang and Y. Yan, Langmuir, 2021, 37, 9009–9016.
13. S. Lepikko, Y. M. Jaques, M. Junaid, M. Backholm, J. Lahtinen, J. Julin, V. Jokinen, T. Sajavaara, M. Sammalkorpi, A. S. Foster and R. H. A. Ras, Nat. Chem., 2024, 16, 506–513.
14. M. Lbadaoui-Darvas, G. Garberoglio, K. S. Karadima, M. N. D. S. Cordeiro, A. Nenes and S. Takahama, Mol. Simul., 2023, 49, 1229–1266.
15. J. Cordeiro and S. Desai, J. Micro Nano-Manuf., 2017, 5, 031008.
16. D. Sergi, G. Scocchi and A. Ortona, Fluid Phase Equilib., 2012, 332, 173–177.
17. L. Muñoz-Rugeles, B. A. Arenas-Blanco, J. M. del Campo and E. Mejía-Ospino, Phys. Chem. Chem. Phys., 2022, 24, 11412–11419.
18. T. Young, Philos. Trans. R. Soc. London, 1805, 95, 65–87.
19. A. Dupré, Théorie mécanique de la chaleur, Gauthier-Villars, Paris, 1869.
20. L. Boruvka and A. W. Neumann, J. Chem. Phys., 1977, 66, 5464–5476.
21. J. W. Drelich, Adv. Colloid Interface Sci., 2019, 267, 1–14.
22. E. Bormashenko, Colloids Surf., A, 2009, 345, 163–165.
23. L. Schimmele, M. Napiórkowski and S. Dietrich, J. Chem. Phys., 2007, 127, 164715.
24. J. Włoch, A. P. Terzyk and P. Kowalczyk, Chem. Phys. Lett., 2017, 674, 98–102.
25. M. Kanduč, L. Eixeres, S. Liese and R. R. Netz, Phys. Rev. E, 2018, 98, 032804.
26. B. Zhao, S. Luo, E. Bonaccurso, G. K. Auernhammer, X. Deng, Z. Li and L. Chen, Phys. Rev. Lett., 2019, 123, 094501.
27. G. Scocchi, D. Sergi, C. D'Angelo and A. Ortona, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 061602.
28. L. Guillemot, T. Biben, A. Galarneau, G. Vigier and E. Charlaix, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 19557–19562.
29. J. Zhang, P. Wang, M. K. Borg, J. M. Reese and D. Wen, Phys. Fluids, 2018, 30, 082003.
30. M. Ma, G. Tocci, A. Michaelides and G. Aeppli, Nat. Mater., 2016, 15, 66–71.
31. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1–2, 19–25.
32. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
33. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
34. L. Martínez, R. Andrade, E. G. Birgin and J. M. Martínez, J. Comput. Chem., 2009, 30, 2157–2164.
35. H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, J. Phys. Chem., 1987, 91, 6269–6271.
36. S. Nosé, Mol. Phys., 1984, 52, 255–268.
37. W. G. Hoover, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31, 1695–1697.
38. R. Hockney and J. Eastwood, Computer Simulation Using Particles, CRC Press, New York, 1st edn, 1988, pp. 120–165.
39. I. C. Yeh and M. L. Berkowitz, J. Chem. Phys., 1999, 111, 3155–3162.
40. J. A. Nieminen, D. B. Abraham, M. Karttunen and K. Kaski, Phys. Rev. Lett., 1992, 69, 124–127.
41. F. F. Abraham and Y. Singh, J. Chem. Phys., 1977, 67, 2384–2385.
42. S. P. Kadaoluwa Pathirannahalage, N. Meftahi, A. Elbourne, A. C. G. Weiss, C. F. McConville, A. Padua, D. A. Winkler, M. Costa Gomes, T. L. Greaves, T. C. Le, Q. A. Besford and A. J. Christofferson, J. Chem. Inf. Model., 2021, 61, 4521–4536.
43. Y. Wu, H. L. Tepper and G. A. Voth, J. Chem. Phys., 2006, 124, 024503.
44. A. J. W. Orry and R. Abagyan, Methods of Protein Structure Comparison, Humana Press, Totowa, NJ, 2012, pp. 231–257.
45. A. Kuzmanic and B. Zagrovic, Biophys. J., 2010, 98, 861–871.
46. J. W. Pitera, J. Phys. Chem. B, 2014, 118, 6526–6530.
47. Z.-X. Niu, T. Huang and Y. Chen, Front. Phys., 2018, 13, 137804.
48. J. Hautman and M. L. Klein, Phys. Rev. Lett., 1991, 67, 1763–1766.
49. M. Khalkhali, N. Kazemi, H. Zhang and Q. Liu, J. Chem. Phys., 2017, 146, 114704.
50. C. F. Fan and T. Caǧin, J. Chem. Phys., 1995, 103, 9053–9061.
51. M. Foroutan, S. M. Fatemi, F. Esmaeilian and V. F. Naeini, Langmuir, 2018, 34, 14085–14095.
52. E. Bertrand, T. D. Blake and J. D. Coninck, J. Condens. Matter Phys., 2009, 21, 464124.
53. T. Blake and J. Haynes, J. Colloid Interface Sci., 1969, 30, 421–423.
54. T. Werder, J. H. Walther, R. L. Jaffe, T. Halicioglu and P. Koumoutsakos, J. Phys. Chem. B, 2003, 107, 1345–1352.
55. C. Li, J. Huang and Z. Li, Sci. Rep., 2016, 6, 26488.
56. M. Foroutan, S. M. Fatemi, F. Esmaeilian, V. Fadaei Naeini and M. Baniassadi, Phys. Fluids, 2018, 30, 052101.
57. M. Foroutan, M. Torabi Rad, A. Boudaghi and H. Ataeizadeh, J. Iran. Chem. Soc., 2022, 19, 423–433.
58. M. Sega, G. Hantal, B. Fábián and P. Jedlovszky, J. Comput. Chem., 2018, 39, 2118–2125.
59. M. Sega, S. S. Kantorovich, P. Jedlovszky and M. Jorge, J. Chem. Phys., 2013, 138, 044110.
60. L. B. Pártay, G. Hantal, P. Jedlovszky, A. Vincze and G. Horvai, J. Comput. Chem., 2008, 29, 945–956.
61. A. Maslov, M. K. Tiwari and D. Seveno, Adv. Eng. Mater., 2024, 2401463.
62. H. Edelsbrunner, D. Kirkpatrick and R. Seidel, IEEE Trans. Inf. Theory, 1983, 29, 551–559.
63. H. Edelsbrunner and E. P. Mücke, ACM Trans. Graphics, 1994, 13, 43–72.
64. A. K. Giri and M. Sega, J. Mol. Liq., 2024, 407, 125155.
65. Y. Wang, A. Kiziltas, P. Blanchard and T. R. Walsh, J. Chem. Inf. Model., 2022, 62, 6302–6308.
66. M. Schmid, D. Rath and U. Diebold, ACS Meas. Sci. Au, 2022, 2, 185–196.
67. M. J. de Ruijter, T. D. Blake and J. De Coninck, Langmuir, 1999, 15, 7836–7847.
68. R. Bey, B. Coasne and C. Picard, J. Chem. Phys., 2020, 152, 094707.
69. M. S. Polovinkin, N. A. Volkov, D. V. Tatyanenko and A. K. Shchekin, Colloids Surf., A, 2024, 702, 134932.
70. C. Oostenbrink, A. Villa, A. E. Mark and W. F. Van Gunsteren, J. Comput. Chem., 2004, 25, 1656–1676.
71. M. Kanduč, J. Chem. Phys., 2017, 147, 174701.
72. A. Boudaghi and M. Foroutan, J. Mol. Liq., 2022, 348, 118017.
73. T. A. Ho, D. V. Papavassiliou, L. L. Lee and A. Striolo, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 16170–16175.
74. M. Darvishi and M. Foroutan, J. Mol. Liq., 2017, 225, 1–10.
75. D. C. Rapaport, The Art of Molecular Dynamics Simulation, Cambridge University Press, Cambridge, UK, 2nd edn, 2004, pp. 245–266.
76. F. P. Beer, E. R. Johnston, D. F. Mazurek, P. J. Cornwell and B. P. Self, Vector mechanics for engineers, McGraw-Hill, New York, US, 12th edn, 2019, pp. 536–537.

---