Wetting and drying of a flat, smooth surface

The wetting behaviour of a Lennard-Jones (LJ) fluid within a slit pore is simulated. The aim is to calculate the density profile \(\rho(x)\) perpendicular to the wall surfaces for a range of wetting parameters \(w\). The set-up and the choice of parameters follows approximately the work of Evans et al. [1].

Prerequisites

You need a proper installation of HAL’s MD package and, for the analysis part, Python with the h5py package. The HDF5 tools (e.g., h5ls) are not required, but of advantage.

These simulation scripts are found in the folder /share/doc/halmd/examples relative to the installation path:

• wetting/potential.lua

• wetting/wetting_equilibration.lua

• wetting/wetting_production.lua

• plotting/density_profile.py

• plotting/plot_density_profile.py

Set-up

A Lennard-Jones (LJ) fluid is confined within a slit pore delimited by two identical planar walls separated by a distance \(D=30\sigma\) along the \(x\)-direction of the simulation box, where \(\sigma\) is the interaction range of the LJ potential. Periodic boundary conditions are applied in the directions parallel to the planar walls.

The wetting behaviour of a LJ fluid is simulated using the smoothly truncated LJ potential for particle pair interactions and the planar wall potential for the particle–wall interaction.

Equilibration phase

Melting of the initial lattice configuration with an NVT simulation:

halmd wetting_equilibration.lua -v \
    --output wetting_thermalisation \
    --density 0.75 \
    --wetting 6.0 \
    --pore-width 30 \
    --temperature 0.9 \
    --time 500 \
    --sampling trajectory=1000

Production phase

The same script is then used for production, with all sampling turned on:

halmd wetting_production.lua -v \
    --output wetting_production \
    --input wetting_thermalisation.h5 \
    --wetting 6.0 \
    --pore-width 30 \
    --time 1000 \
    --sampling trajectory=1000 structure=1000

Analysis of density profiles

The density profile can be obtained from the data on the wavevectors compatible with the reciprocal space of the periodic simulation box and the particles’ positions. The Python script provided at plotting/density_profile.py extracts this data from the output file wetting_product.h5 generated during the production phase to compute the 1D density profile. Another Python script, plotting/plot_density_profile.py, is employed to visualise the computed profile. An example call of the latter script is:

python plot_density_profile.py --inputs wetting_production.h5

The calculation of the 1D density profile in plotting/density_profile.py involves the following steps:

Firstly, the complex Fourier modes are computed from the wavevectors and particle positions along the x-direction (coord_index=0).

def compute_density_modes_1D(wavevector: np.ndarray, positions: np.ndarray, coord_index: int) -> np.ndarray[complex]:
    """
    Compute time-averaged density modes in 1D.

    Args:
    - coord_index: index of coordinate of interest (coord_index ∈ {0, 1, 2})
    """

    density_modes =  np.zeros(wavevector.shape[0], dtype = complex)
    for i, k_i in enumerate(tqdm(wavevector)):
        density_modes[i] = np.sum(np.exp(1j * k_i[coord_index] * positions[:,:, coord_index]))

    return density_modes / positions.shape[0]

Next, the wavevectors are transformed to the corresponding grid indices of the simulation box in real space.

def wavevector_to_grid_indices(wavevector: np.ndarray[float], box_edges: np.ndarray[float]) -> np.ndarray[float]:
    """
    Map wave indices to the appropriate indices in a discrete Fourier grid.

    Args:
    - wavevector: the wave indices to map, given as an array of shape (n_vectors, 3)
    - box_edges: the physical dimensions of the real-space grid, given as (Lx, Ly, Lz)

    Returns:
    - grid_indices: an array containing the grid indices for all wavevectors
    """

    # Compute the grid spacing in reciprocal space
    dkx, dky, dkz =2 * np.pi / box_edges

    # Map the wavevector components to grid indices
    wx = np.array(np.round(wavevector[:,0] / dkx), dtype =int)
    wy = np.array(np.round(wavevector[:,1] / dky), dtype =int)
    wz = np.array(np.round(wavevector[:,2] / dkz), dtype =int)
    assert np.min(wx)==-np.max(wx) , "Density-modes need to be on a symmetric grid around 0, i.e. k_max = -k_min in x- direction"
    assert np.min(wy)==-np.max(wy) , "Density-modes need to be on a symmetric grid around 0, i.e. k_max = -k_min in y- direction"
    assert np.min(wz)==-np.max(wz) , "Density-modes need to be on a symmetric grid around 0, i.e. k_max = -k_min in z- direction"

    # Combine wx, wy, wz into a single array of shape (199, 3)
    grid_indices = np.column_stack((wx, wy, wz))

    return grid_indices

In the compute_1D_density function, inverse Fourier transformation is applied on the modes to obtain the density. Due to symmetry, the density at the two ends of the pore slit is set to be equal. The function returns the normalised density with respect to bulk density and the position grid.

def wavevector_to_grid_indices(wavevector: np.ndarray[float], box_edges: np.ndarray[float]) -> np.ndarray[float]:
    """
    Map wave indices to the appropriate indices in a discrete Fourier grid.

    Args:
    - wavevector: the wave indices to map, given as an array of shape (n_vectors, 3)
    - box_edges: the physical dimensions of the real-space grid, given as (Lx, Ly, Lz)

    Returns:
    - grid_indices: an array containing the grid indices for all wavevectors
    """

    # Compute the grid spacing in reciprocal space
    dkx, dky, dkz =2 * np.pi / box_edges

    # Map the wavevector components to grid indices
    wx = np.array(np.round(wavevector[:,0] / dkx), dtype =int)
    wy = np.array(np.round(wavevector[:,1] / dky), dtype =int)
    wz = np.array(np.round(wavevector[:,2] / dkz), dtype =int)
    assert np.min(wx)==-np.max(wx) , "Density-modes need to be on a symmetric grid around 0, i.e. k_max = -k_min in x- direction"
    assert np.min(wy)==-np.max(wy) , "Density-modes need to be on a symmetric grid around 0, i.e. k_max = -k_min in y- direction"
    assert np.min(wz)==-np.max(wz) , "Density-modes need to be on a symmetric grid around 0, i.e. k_max = -k_min in z- direction"

    # Combine wx, wy, wz into a single array of shape (199, 3)
    grid_indices = np.column_stack((wx, wy, wz))

    return grid_indices

To observe the normalised density profiles at vapour liquid coexistence for various wetting parameters \(w\), one should first execute the two Lua simulation scripts (wetting/wetting_equilibration.lua and wetting/wetting_production.lua) for \(w \in \{0.1, 0.5, 1.0, 1.5, 3.0, 6.0\}\). Once the output files are obtained, the profiles can be visualised in a single plot by running the following command:

python plot_density_profile.py --inputs wetting_production_*.h5

[1]

Robert Evans, Maria C. Stewart, Nigel B. Wilding, Drying and wetting transitions of a Lennard-Jones fluid: Simulations and density functional theory, J. Chem. Phys. 147, 044701 (2017) [Link]