Contact angle determination as a function of solid–liquid molecular characteristics via the Kirkwood–Buff route

Abstract

The goal of this work is to calculate contact angles as a function of the molecular characteristics of the solid and the liquid. We use Kirkwood–Buff's mechanical route and decompose the normal and transverse components of the pressure tensor into solid and liquid regions. This decomposition leads to the emergence of a new force – the repulsion induced by the liquid displaced by the solid –. We found a neutral state where liquid/wall and displaced-liquid/liquid attractions are balanced. By expanding around this state, we could calculate the density profile within the hard-core van der Waals model at low isothermal compressibility, and establish a relation between the equilibrium contact angle and the microscopic properties of the liquid and solid. This work leads to reformulate Young's equation in a novel manner and sheds new light on the physics of wetting. The density profiles and contact angles obtained theoretically are supported by molecular dynamics simulations using the Lennard-Jones potential.

---

Introduction

The first attempt to determine, from first principles, the contact angle in a solid–liquid–vapor system as a function of the molecular characteristics of the solid and the liquid was made by M. Berry.^1,2 Previous approaches^3,4 were insightful but they were not based on first principles. In Berry's work,¹ the Born–Green–Yvon (BGY) equation^5–7 was used to construct an expression for the density profile. Once this was established, a “surface tension profile” and ultimately the contact angle could be calculated. Although this work does not provide a ‘solution’ to the contact angle problem—since the density was a trial function, not an actual solution of the BGY equations—it nevertheless represents a significant attempt to relate, from first principles, the contact angle to the microscopic properties of the fluid and the solid. To the author's knowledge, this work has not been followed by subsequent developments. One inherent difficulty in these approaches is the fact that long-range interactions give rise to nonlinear integral equations, which are challenging to solve.

Later, Benner et al.⁸ addressed the problem using gradient theory.^9,10 For solid–liquid interfaces, however, a key difficulty lies in the lack of a well-defined boundary condition at the wall. Although calculations can still be carried out, the reliability of the results remains uncertain.

One often finds, in textbooks, expressions—see, for instance,¹¹—or order-of-magnitude arguments—for example,¹²—that relate interfacial tensions to the molecular characteristics of the system. The issue is that, in liquids, pressure has two components: a kinetic one, associated with the thermal motion of the molecules, and a static one, arising from intermolecular forces. These two components are of comparable magnitude. If the kinetic component is neglected—thus treating the liquid as a solid—the problem becomes greatly simplified, but a key element of the physics of the system is lost.

Young's equation¹³ offers a means to calculate the contact angle by balancing the interfacial tensions between solid–liquid, solid–vapor, and liquid–vapor interfaces. However, because no expressions derived from first principles are currently available, the equation cannot, in practice, be used to compute the contact angle directly. In the field of wetting, a theoretical expression for the contact angle would be highly valuable for structuring reasoning, interpreting experimental results, generalizing findings across different systems, and ultimately devising strategies for optimizing materials. Numerical simulations—such as molecular dynamics (MD) or density functional theory (DFT)^14–16—provide invaluable data, but they remain black boxes and offer little physical insight or guidance.

The goal of this work is to derive, analytically and from first principles, within a certain framework of approximations, the relationship between the equilibrium contact angle and the molecular characteristics of the liquid and the solid. We will make use of the Kirkwood–Buff equations,^2,6,17–19 which appear particularly well suited to this task, yet have not been thoroughly exploited for this purpose to date.

Geometry and intermolecular forces

The geometry under consideration is shown in Fig. 1A. A liquid occupies the space above a wall located at z = 0. To represent the intermolecular forces within the liquid, induced by its own environment, we consider the case of nonpolar liquids and model the interactions using a hard-sphere truncated van der Waals potential. The expressions for the potential V(r) and the corresponding force f(r), where r denotes the intermolecular separation, are given by:in which is the interaction energy, and σ the Lennard Jones (LJ) diameter. We use similar expressions for the forces exerted by the solid onto the fluid. We thus have, in the region above the wall (z > 0):with, according to Berthelot's rule , and , where and σ_S are, respectively, the solid energy and solid LJ diameter.

Fig. 1 Geometry of the problem. We have M₁(x₁, y₁, z₁) and M₂(x₂, y₂, z₂). (A) The solid (grey) occupies the region z < 0, while the liquid lies above, in z > 0. (B) p_N(z) is the normal pressure acting on dS, resulting from the attractive forces exerted by the solid and the liquid located below the dashed line (i.e., for z₁ < z) on the fluid located at z₂ > z. (C) p_T(z) is the transverse pressure at dS, resulting from the attraction exerted by the solid and fluid located at y₁ < 0 (i.e. at the left of the vertical dashed line) on the solid and fluid located at y₂ > 0, i.e. at the right of the same line.

To avoid cumbersome notation, we work in a dimensionless Lennard-Jones system of units. We thus use dimensionless coordinates (and similar expressions for the other coordinates), density ρ* = ρσ³, pressure , temperature , where k_B is the Boltzmann constant, force , potential and remove the stars.

In dimensionless LJ units, for r > 1, the liquid–liquid and solid–liquid van der Waals forces are given by:

f(r) = −24r⁻⁷ and f_SL(r) = −24ζr⁻⁷ = ζf(r) (3)

in which .

Calculation of the density profiles

Owing to symmetries, the pressure tensor is diagonal. We thus have p_xx = p_yy = p_T and p_zz = p_N in which p_T and p_N are respectively, the transverse and normal components of the pressure. We use Kirkwood–Buff's equations (KB)¹⁷ which state that, in the fluid, the normal and transverse pressures are given by (in LJ units):

p_N(ρ) = ρT + p^(S)_N (4)

p_T(ρ) = ρT + p^(S)_T (5)

in which ρ is the fluid density and p^(S)_N and p^(S)_T are the static terms. The first term in the equations (the kinetic term) arises from molecular collisions, while p^(S)_N and p^(S)_T represent contributions from intermolecular forces. Importantly, these equations have been rigorously derived from statistical theory.¹⁸ Although more suited to the analysis of real gases, eqn (5), like Born–Green–Yvon (BGY) equation, are commonly used to investigate liquids. Examples include the determination of contact angles,^1,2,20 and the calculation of density profiles and surface tension at liquid–vapor interfaces.⁶

The expressions of p^(S)_N and p^(S)_T have been established by KB.¹⁸ For p_N, the principle of the calculation is sketched in Fig. 1B: it consists in summing up attractions of molecules at r₁(x₁, y₁, z₁) and r₂(x₂, y₂, z₂) below and above a plane parallel to the solid surface, for which the separation vector r = M₁M₂ = r₂ − r₁ crosses the surface element dS (see Fig. 1B).^6,18 A similar calculation hold for p_T, with a plane normal to the solid surface (see Fig. 1C). The expression of p^(S)_N and p^(S)_T, as formulated by Irving and Kirkwood,^6,19 read:

in which x₁₂, y₁₂, z₁₂ are the r coordinates and ρ⁽²⁾ is the two particle distribution functions between positions r₁ and r₂; depending on whether M₁ or M₂ are located in the liquid or in the solid, f*(r) is the liquid–liquid or solid–liquid interaction force. We will assume the following expression for ρ⁽²⁾:

ρ⁽²⁾(r, z_A, z_B) = ρ(z_A)ρ(z_B)H(r) (8)

in which ρ(z_A) and ρ(z_B) denote the mean-field densities at positions z_A and z_B, respectively, and H(r) is the Heaviside step function with a discontinuity at r = 1. This choice, which reflects the hard-sphere van der Waals potential used (see eqn (1)), inherently prevents the development of layering effects.

In eqn (6), the pressure applied on dS originates from molecules located in the solid and in the portion of fluid sandwiched between the solid and the plane at height z (see Fig. 1B). A similar remark holds for eqn (7). The idea here is to decompose p^(S)_N into two contributions: one arising from the sandwiched fluid, and the other from the solid. For that, we abandon α and return to the variable z₁, using the relation z₁ = z − αz₁₂. This leads to the following expressions of p^(S)_N:

Here, we have extended the domain of definition of ρ(z), originally limited to z > 0, i.e., from 0 to +∞, to the entire z axis (−∞, ∞), by imposing the symmetry condition ρ(−z) = ρ(z). The extended fluid density is continuous at z = 0, while its odd-order derivatives are discontinuous. The first integral thus represents the internal attractive forces inside the extended liquid, i.e. the liquid occupying the entire space.

The second integral includes two contributions: (i) the attraction exerted by the solid on the liquid located above the plane at z, and (ii) the repulsion exerted by a virtual fluid occupying the solid region z < 0 on that same liquid. We refer to this virtual fluid as the displaced liquid, intentionally using the same terminology as in Archimedes’ principle. To the best of the author's knowledge, the decomposition (9) is new.

We now take advantage of the cylindrical symmetry of the system to perform the integration over x₁₂ and y₁₂, with z₁₂ ≥ 1. This leads to the following expression for p^(S)_N:

Far from the wall, the fluid density is homogeneous, equal to ρ_L, and the pressure uniform, equal to p₀. Using eqn (6), we have:

p₀ = ρ_LT + p^(S)_N(ρ_L) = ρ_LT − aρ_L² (11)

with . The same relation holds for the transverse component. represents the contribution, at high temperature, of the attractive intermolecular interactions to the second virial coefficient, based on LJ potential.²¹

Mechanical equilibrium requires

p₀ = ρT + p^(S)_N(ρ) (12)

eqn (12) defines an integro-differential equation, which must be solved to calculate the density profile ρ(z). The key point is that when the parametervanishes, eqn (12) admits an exact solution, given by ρ = ρ_L. This can be checked by inspecting eqn (10): when Λ = 0 and ρ = ρ_L, the second integral term cancels out and the normal component of the pressure becomes p_N(ρ_L) = ρ_LT − aρ_L² = p₀, which is satisfied (see eqn (11)). Therefore, when Λ = 0, ρ(z) = ρ_L is an exact solution to eqn (12). We will call the corresponding state ‘neutral’.

The solution at small Λ can be found by expanding densities and pressures into series of increasing powers of this parameter:

At zeroth order, eqn (11) is satisfied: in this state, the solid attraction is balanced by the repulsion developed by the displaced fluid and the system behaves as an homogeneous fluid occupying the entire space: this is the neutral state. At first order, we have the following relations:

ρ_Ls(z)T + p_N⁽¹⁾(z) = 0 (15)

After calculation, we obtain:

where κ_T is the isothermal compressiblity and

We now focus on the case of weakly compressible fluids, i.e., such that κ_T ≪ 1. In practice, this condition is fulfilled by liquids well below the critical point and gases at high pressures. In the limit of small κ_T, s(z) ∼ κ_T and K(z) ∼ κ_T becomes negligible. As a consequence, the eqn (16) loses its integro-differential term and its resolution becomes trivial. This leads to:

We thus obtain, at first order in Λ and κ_T:

eqn (19) tells that, in the interfacial layer, the density excess or depletion decays as z⁻³. The density profiles (eqn (19)) are shown in Fig. 2A for ζ varying between 0 and 0.5 (with ρ_L = 0.8, ρ_S = 3.2 and κ_T = 0.7). The corresponding values of Λ are: −1, −0.5, 0, 0.5, 1.

Fig. 2 (A) Theoretical density profiles (eqn (19)) for ρ_L = 0.8, ρ_S = 3.2, κ_T = 0.7 and various ζ, varying from 0 to 0.5 (shown on the figure). The corresponding values of Λ are: −1, −0.5, 0, 0.5, 1. Insert: Density profile for Λ = −1 and κ_T = 1.5, showing drying. (B) Evolution of cos θ with ζ (eqn (26)), for ρ_S = 2.5 and three temperatures T = 1 (blue), 1.15 (red), 1.3 (black), with ρ_L = (0.7, 0.61, 0.4), and γ = (0.5, 0.15, 0.013). (C) Evolution of θ with temperature, for Λ = −1 (black line) and 1.2 (blue line), using the same ρ_L(T) and γ(T) as for (B).

For positive Λ (ζ > 0.25), the fluid density lies above ρ_L Interestingly, at the largest Λ, densities are so large that we may expect a solid-like layer to form at the wall—a phenomenon known as van der Waals solidification.^6,11 For (Λ = 0 (ζ = 0.25)), we recover the neutral state, where the liquid density remains equal to ρ_L throughout the fluid domain. For negative Λ (ζ < 0.25), the liquid becomes depleted, and for sufficiently negative values, a vapor-like layer forms. Eventually, as shown in the inset (with κ_T = 1.5, assuming the small-κ_T approximation remains valid), drying occurs when .

Calculation of the solid–liquid interfacial tension and the contact angle

We now concentrate on eqn (7) and (10). At first order in Λ, in the small compressiblity limit, where, again, terms on the order of κ_T are neglected against terms of order unity, p^(S)_T(z), p^(S)_N(z) and the excess pressure Π(z) = p^(S)_N(z) − p^(S)_T(z) read:where

Consequently, the solid–liquid interfacial tension γ_SL reads:

where

At first order, the transverse static pressure is therefore on the order of 10% of the normal pressure. This can be understood by noting that, in the average, the distances between the fluid and solid elementary volumes involved in the calculation of p_N⁽¹⁾ are shorter than for p_T⁽¹⁾.

The interfacial tension γ_SL can be decomposed in the following manner:

γ_SL = γ_DL − γ_SL0 (24)

in which γ_DL is the displaced-liquid/liquid and γ_SL0 the ‘bare’ solid/liquid interfacial tensions.

The displaced-liquid/liquid surface tension is a novel concept. To illustrate its physical significance, we may consider a liquid initially at equilibrium, occupying all the space. Its density is uniform, equal to ρ_L. Creating a solid–liquid interface involves two steps: (i) removing part of the liquid – which depletes its density close to the newly created interface – and (ii) bringing the solid into contact with the new interface – which compresses the liquid near the wall. In this process, the displaced liquid refers to a virtual liquid occupying the volume of fluid that has been removed in order to accommodate the solid. The resulting density profile (eqn (19)) adds two contributions: one arising from the removal of the liquid (depletion), and another from its replacement by the solid (compression). They add up. In the same manner, the interfacial tension components of eqn (24), namely γ_DL and γ_SL0, add up to define the total interfacial tension γ_SL.

Eqn (22) warrants two remarks:

(i) In the case Λ > 0, γ_SL being negative, each fluid element is subjected to a compressive stress in the x, y plane. For liquid/vapor and liquid/liquid interfaces, such stresses cannot exist because the system would be unstable: a deformation of the interface generates a force which, in turn, increases the deformation. In our case, the instability is suppressed by the wall.

(ii) For γ_SL < 0, an elementary volume of fluid, placed on the surface and isolated from the surrounding fluid, will develop an internal force that causes it to spread. Λ > 0 is thus associated to wetting. In the opposite case, i.e. Λ < 0, γ_SL is positive and the internal force will cause it to dewet the surface. Λ < 0 is thus associated to non wetting.

We now calculate contact angles. For that, we first need to consider the case where the solid is in presence of the vapor. Well below the critical point, since gas densities are low, eqn (22) leads to γ_SG ≪ γ_SL, consistently with the MD simulations.^22,23 In this problem, the dry surface plays no role; it is merely a spectator in the process that determines the contact angle.

Thus, neglecting γ_SG and expressing mechanical equilibrium in the far field – as is appropriate in the derivation of Young's law²⁴ – we obtain (see Fig. 3):

γ_SL = −γ cos θ (25)

in which γ is the liquid–vapor surface tension (in LJ dimensionless units). Using eqn (20), (22) and (25), we find the following expression for cos θ:

Fig. 3 According to Newton's second law, γ_SG being negligible, the forces (shown in red) exerted by the liquid surroundings on the liquid volume under consideration (shown in blue) must balance. This leads to eqn (25). (A) Λ > 0: γ_SL < 0, acute contact angle; (B) Λ < 0: γ_SL > 0, obtuse contact angle.

eqn (26) expresses, in analytical from, the contact angle in function of the molecular characteristics of the liquid and the solid. This was the goal of the paper and this is the first time that a relation of this type is derived from first principles.

To illustrate the behavior of eqn (26), Fig. 2B shows the evolution cos θ for ρ_S = 2.5 and three different temperatures. The corresponding values of ρ_L and γ are provided in the caption – in practice, to determine ρ_L, we used the 32-parameters EOS of ref. 25 and, for γ, we used ref. 26. The linearity of cos θ with ζ, in the partial wetting range (see Fig. 2B) is a consequence of eqn (26). As we approach the critical point (assuming the theory remains accurate), the curves in Fig. 2B tend toward a stepwise function, indicating that, depending on ζ (in fact on the sign of Λ), the system becomes either fully wetting or dewetting. This is critical wetting or drying.²⁷ The same phenomenon is visible in Fig. 2C, where θ is plotted as a function of temperature, for two Λ of opposite signs.

Young's equation

By using (24), Young's eqn (25) can be written in the following manner:

γ cos θ = γ_SL0 − γ_DL (27)

In this form, the equation expresses a balance between (i) the liquid/solid, (ii) the displaced-liquid/liquid interfacial tensions and (iii) the liquid/vapor surface tension. Eqn (27) represents an alternative form of Young's equation. It tells that when solid liquid interactions, represented by γ_SL0, prevails, complete wetting is favored. In the opposite case, when liquid cohesion, represented by γ_DL, dominates, dewetting is favored. The partial wetting domain is defined by the condition:

−γ < γ_SL0 − γ_DL < γ (28)

The theory thus recovers the three established wetting regimes: complete wetting, partial wetting, and dewetting. Note that eqn (27) is not the only way to express Young's equation. Alternative formulations, closer to the traditional one, can also be found. However, deriving and discussing them would go beyond the scope of this paper.

Comparison with molecular dynamics simulations performed with LJ fluids

We now compare our theoretical predictions with numerical simulations. The most relevant and well-documented studies of liquid/solid interfacial regions are based on the Lennard-Jones potential.^{22,23,27–33} Several of these studies^{22,27,28,30,31} employed the 9-3 wall potential, explored different geometries, or used Lennard-Jones truncated and shifted (LJTS) model. These are different from the assumptions of the present theory. Others^23,32,33 used the standard 12-6 Lennard-Jones potential for both the liquid and the solid and are therefore closer to the situation we consider. Focusing on a [0.7–1] range of temperature, where compressiblity is small, we build a data basis including 21 points, with contact angles varying from 30 to 120°, ζ between 0.4 and 0.9, ρ_L between 0.7 and 0.84 and ρ_S between 0.9 et 7. In such conditions, κ_T does not exceed 0.4. Note that in these studies, as the contact angles are measured on droplets larger than ten nanometers, line tension corrections³⁴ are neglected.

The numerical results are compared with the density profile and the contact angle given by formulas (19) and (26). To determine ρ_L as a function of temperature, we use the 32-parameters EOS of ref. 25. The same EOS allows to calculate κ_T. The liquid–vapor surface tension γ is determined by using ref. 26.

The theoretical and numerical density profiles are shown in Fig. 4A and B. The theory (black lines) being mean field, we averaged out the numerical profiles (dashed red lines), over one crystal period (black dots). Fig. 4A shows acceptable agreement between theory and simulation for the wetting case (ζ = 0.9).

Fig. 4 Theory (cos θ_theor) compared to LJ molecular dynamics simulations (cos θ_MD); (A) density profiles: red dashed line: molecular simulation³² for T = 1, ζ = 0.9; ●: same profile, averaged over a period (z error bar: ; ρ error bar: , with N the molecule number of ref. 32); full line: theoretical profile (eqn (19)); (B) same as A, with ζ = 0.5; (C) cosinus of the contact angles cos θ_theo (eqn (26)) plotted against MD simulations (cos θ_MD):^23,32,33,³³, , , ◀, ,²³ ×.³² Note that in the comparison, we increased Bo-Shi's temperatures by 10% to match the LJ critical point.

In the non-wetting case (ζ = 0.5), although the theory underestimates the depletion near the wall, it remains within the error bars. We thus observe, in the two cases, quantitative agreement.

Fig. 4C compares the contact angles predicted by the theory (see eqn (26)) with those issued from molecular dynamics simulations.^23,32,33 They are close but different. The distribution of the contact angle differences shows a mean value of −20° and a standard deviation of 10.2°. At this stage, we may conclude that theory and molecular simulations are consistent. The discrepancy may arise from several factors – for example, differences in the interaction potentials (hard-sphere van der Waals in the theory versus Lennard-Jones in the simulations), in the pair distribution functions (hard-sphere in the theory versus soft-core LJ repulsion in the simulations), or inherent limitations of the Kirkwood–Buff approach in accurately representing a liquid. More work is needed to assess the situation.

Introducing corrections for matching MD simulations and theory

From a practical perspective, and to quantify the discrepancy between numerical results and theory in a simple manner, without altering the structures of the theoretical expressions, we amend the expressions of Λ and cos θ as follows:in which q₁ and q₂ are tuning parameters. These parameters are purely empirical. With q₁ = −3T + 3.9 and q₂ = 0.55, we obtain Fig. 5.

Fig. 5 (A) Same as Fig. 4C, but plotting θ (instead of cos θ), and using formula (29), with q₁ = −3T + 3.9 and q₂ = 0.55; theoretical contact angles θ (degree) plotted against numerical data;^23,32,33 (B) and (C) density profiles: same as Fig. 4A and B, but using formula (29); (D) contact angles plotted with temperature, for three ζ: 0.817 , 0.653 , 0.490 ■;²³ full lines are the theoretical lines, obtained from eqn (29). (E) contact angles θ,²³ plotted with ζ for T = 0.7 (red squares) and T = 1 (blue squares); full lines: eqn (29) at the same temperatures.

The figure shows that, with a correction of approximately 40% on q₁ and q₂ (compared to their theoretical value q₁ = q₂ = 1), the data collapse closely onto the line of perfect agreement, with a standard deviation of 6°. Similar comments apply to the density profiles (Fig. 5B and C) and the evolution of the contact angle with respect to T and ζ (Fig. 5D and E).

Conclusion

To summarize, starting from equations derived from the statistical theory of liquids (Kirkwood–Buff equations), and using a hard-core van der Waals potential, we have solved perturbatively, in the limit of small compressibility, the mechanical equations governing the equilibrium of a liquid in contact with a solid. We obtained analytical expressions for both the density profile and the contact angle, in function of the molecular characteristics of the solid and the liquid. The expressions we found are new and are supported by MD simulations.

This study introduces a new concept: the ‘displaced liquid’. It refers to a virtual liquid that occupies the space left by the liquid pushed aside to accommodate the solid. In a manner evoking Archimedes principle, this displaced liquid develops a repulsive force onto the working fluid. The balance between this repulsive force and the attractive forces from the solid determines the density profile and the pressure field in the (working) fluid. This framework naturally leads to a reformulation of Young's equation.

This work sheds new light on the physics of wetting, offering a more intuitive description of the phenomenon than existing approaches.^35,36 From a practical standpoint, the contact angle expression we obtain could help optimize devices where wetting plays a crucial role, including applications in materials science, surface engineering, biomedical fields, and food science.

The present work concerns the mechanics of the interfacial layer and it would be of interest to see how it relates to the energetic approach extensively developed in the literature since Gibbs' seminal paper.³⁷ In the future, it would be interesting to extend the theory to water, where polar interactions dominate. The next step would be to investigate slippage, which strongly depends on wetting.³⁸ The present approach could also be extended more complex situations, for example involving two immiscible liquids or the presence of films.

Conflicts of interest

There are no conflicts to declare.

Data availability

This study was carried out using publicly available data from ref. 22, 23, 27, 28 and 30–33.

Acknowledgements

The author is much indebted to J. F. Joanny and L. Bocquet for their comments. He thanks H. Stone, F. Brochard, D. Bonn, P. Nghe, Y. Pomeau, B. Andreotti, J. O. Fossum, A. Colin and J. Treiner for valuable inputs. Part of this work has been performed during visit in NTNU (Trondheim) (Prof. Fossum) and in University of Napoli, ITT (Prof. Torino and Prof. Netti). Their hospitality is warmly acknowledged. The author is grateful to IPGG, ESPCI PSL and CNRS for their support.

References

1. M. V. Berry, Simple fluids near rigid solids: statistical mechanics of density and contact angle, J. Phys. A: Math., Nucl. Gen., 1974, 7, 231.
2. G. G. Navascues, Liquid surfaces: theory of surface tension, Rep. Prog. Phys., 1979, 42, 1131.
3. L. A. Girifalco and J. G. Robert, A theory for the estimation of surface and interfacial energies. I. Derivation and application to interfacial tension, J. Phys. Chem., 1957, 61, 904.
4. F. M. Fowkes, Attractive Forces at Interfaces, Ind. Eng. Chem., 1964, 56, 40.
5. J. Yvon, La theorie statistique des fluides et l’equation d’etat, Hermann, Paris, 1937.
6. J. S. Rowlinson and B. Widom, Molecular Theory of Capillarity, Oxford University Press, 1989.
7. M. Green and H. S. Green, A General Kinetic Theory of Liquids I: The Molecular Distribution Functions, Proc. R. Soc. London, Ser. A, 1956, 188(1012), 10.
8. L. Benner, L. Scriven and H. Davis, Structure and stress in the gas–liquid–solid contact region, Faraday Symposia of the Chemical Society, Royal Society of Chemistry, 1981, vol. 16.
9. J. D. van der Waals, Over de Continuiteit van Den Gas En Vloeistoftoestand. (On the Continuity of the Gas and Liquid State), PhD thesis, Leiden University, 1873.
10. J. W. Cahn and J. E. Hilliard, Free Energy of a Nonuniform System. I. Interfacial Free Energy, J. Chem. Phys., 1958, 28, 258.
11. J. N. Israelachvili, Intermolecular and Surface Forces, Academic Press, London, 1985.
12. P. G. De Gennes, F. Brochard-Wyart and D. Quere, Capillarity and Wetting Phenomena: Drops, Bubbles, Pearls, Waves, Springer, 2004.
13. T. Young, An essay on the cohesion of fluids, Philos. Trans. R. Soc. London, 1805, 95, 65.
14. M. E. Tuckerman, Recent Advances in First-Principles Based Molecular Dynamics, Acc. Chem. Res., 2021, 54(1), 2.
15. D. Sibanda, S. Temitope Oyinbo and T.-C. Jen, A review of atomic layer deposition modelling and simulation methodologies: Density functional theory and molecular dynamics, Nanotechnol. Rev., 2022, 11(1), 1332.
16. Y. Yant, A. K. Narayanon Nair, S. Sun and D. Laul, Estimating fluid-solid interfacial free energies for wettabilities: A review of molecular simulation methods, Adv. Colloid Interface Sci., 2025, 103482.
17. J. G. Kirkwood, The statistical mechanical theory of transport processes I. General theory, J. Chem. Phys., 1946, 14, 180.
18. J. G. Kirkwood and F. P. Buff, The statistical mechanical theory of surface tension, J. Chem. Phys., 1949, 17, 338.
19. J. H. Irving and J. G. Kirkwood, The Statistical Mechanical Theory of Transport Processes. The Equations of Hydrodynamics, J. Chem. Phys., 1950, 18, 817.
20. G. Navascues and M. V. Berry, The statistical mechanics of wetting, Mol. Phys., 1977, 34(3), 649.
21. J. O. Hirschfelder, C. F. Curtiss and R. B. Bird, Molecular Theory of Gases and Liquids, Wiley, New York, 1954.
22. M. Nijmeijer, C. Bruin, A. Bakker and J. van Leeuwen, Molecular dynamics of the wetting and drying of a wall with a long-ranged wall-fluid interaction, J. Phys.: Condens. Matter, 1992, 4, 15.
23. S. Bo and V. K. Dhir, Molecular dynamics simulation of the contact angle of liquids on solid surfaces, J. Chem. Phys., 2009, 130, 034705.
24. P. G. de Gennes, Wetting: statics and dynamics, Rev. Mod. Phys., 1985, 57, 827.
25. S. Pieprzyk, A. C. Branka, S. Mackowiak and D. M. Heyes, Comprehensive representation of the Lennard-Jones equation of state based on molecular dynamics simulation, J. Chem. Phys., 2018, 148, 114505.
26. S. Stephan, M. Thol, J. Vrabec and H. Hasse, Thermophysical Properties of the Lennard-Jones Fluid: Database and Data Assessment, J. Chem. Inf. Model., 2019, 59(10), 4248.
27. S. Becker, H. M. Urbassek, M. Horsch and H. Hasse, Contact angle of sessile drops in Lennard-Jones systems, Langmuir, 2014, 30(45), 13606.
28. J. H. Sikkenk, O. Indekeu, M. J. van Leeuwen, E. O. Vossnaek and A. F. Bakker, Simulation of Wetting and Drying at Solid-Fluid Interfaces on the Delft Molecular Dynamics Processor, J. Stat. Phys., 1988, 52, 23.
29. S. Ravipati, B. Aymard, S. Kalliadasis and A. Galindo, On the equilibrium contact angle of sessile liquid drops from molecular dynamics simulations, J. Chem. Phys., 2018, 148, 164704.
30. K. Bucior, L. Yelash and K. Binder, Molecular-dynamics simulation of evaporation processes of fluid bridges confined in slitlike pores, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2009, 79, 031604.
31. R. Evans, M. C. Stewarta and N. B. Wilding, A unified description of hydrophilic and superhydrophobic surfaces in terms of the wetting and drying transitions of liquids, Proc. Natl. Acad. Sci. U. S. A., 2019, 116(48), 23901.
32. L. Bocquet and J. L. Barrat, Large Slip Effect at a Nonwetting Fluid-Solid Interface, Soft Matter, 2007, 3, 685.
33. Q. Li, B. Wang, Y. Chen and Z. Zhao, Wetting and evaporation of argon nanodroplets on smooth and rough substrates: Molecular dynamics simulations, Chem. Phys. Lett., 2016, 662, 73.
34. T. Werder, J. H. Walther, R. L. Jaffe, T. Halicioglu and P. Koumoutsakos, On the water-carbon interaction for use in molecular dynamics simulations of graphite and carbon nanotubes, J. Phys. Chem. B, 2003, 107(6), 1345.
35. M. Berry, The molecular mechanism of surface tension, Phys. Educ., 1971, 6, 79.
36. A. Marchand, J. H. Weijs, J. H. Snoeijer and B. Andreotti, Why is surface tension a force parallel to the interface?, Am. J. Phys., 2011, 79(10), 999.
37. J. W. Gibbs, in The Collected Works of J. Willard Gibbs: Thermodynamics, ed. W. R. Longley and R. G. Van Name, Longmans, Green and Co., New York, 1928, vol. I (reprinted by Yale University Press, 1948).
38. P. Tabeling, Introduction to Microfluidics, Oxford Univ. Press, 2nd edn, 2023.

---