# Tutorial: wetting and drying at 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
perpendicular to the wall surfaces for a range of wetting parameters .
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 along the -direction of the simulation box, where 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 , one should first execute the two Lua simulation scripts
(`wetting/wetting_equilibration.lua`

and `wetting/wetting_production.lua`

)
for . 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] |