Create fluid in a slit pore

A fluid confined to a slit pore is set up by placing fluid particles in between two planar surfaces. The latter are implemented as flat walls by an external wall potential. The fluid may consist of particles of different species, and the walls may have different interaction parameters for each species. For simplicity, only a single species is used here.

Note

The simulation box is still defined with periodic boundary conditions at all its faces. Therefore, it is crucial that the box size in the direction perpendicular to the walls is larger than the pore width to avoid spurious pair interactions across the walls.

The setup procedure has the following steps:

• define geometry: pore width and box extents, fraction of the pore

• place particles on an fcc lattice in the pore part of the box

• assign random velocities according to a Maxwell–Boltzmann distribution

• define interactions: pair potential, external wall potential

local length = {pore_width + 10, box_length, box_length}   -- add extra space "behind" the walls
local slab = {pore_width / length[1], 1, 1}
-- compute particle number from pore volume and mean density
local nparticle = math.floor(args.density * numeric.prod(slab) * numeric.prod(length))

-- create system state
local box = mdsim.box({length = length})
local particle = mdsim.particle({dimension = #length, particles = nparticle})

-- set initial particle positions
mdsim.positions.lattice({box = box, particle = particle, slab = slab}):set()

-- set initial particle velocities
mdsim.velocities.boltzmann({particle = particle, temperature = args.temperature}):set()

-- define Lennard-Jones pair potential (with parameters ε=1 and σ=1 for a single species)
-- and register computation of pair forces
definitions.lennard_jones.create_pair_force({
    box = box, particle = particle, cutoff = args.pair_cutoff, smoothing = args.smoothing
})

-- define Lennard-Jones wall potential to form a slit pore
-- and register computation of this external force
definitions.slit_pore.create_wall_force({
    box = box, particle = particle
  , pore_width = pore_width
  , epsilon = args.wall_epsilon, wetting = args.wetting
  , cutoff = args.wall_cutoff, smoothing = args.smoothing
})
-- end of usage in doc/recipes/create_slit_pore.rst.in

We have made use of the following helper functions, which may be found in separate definitions file (see path examples/definitions):

function M.create_wall_force(args)
    local box = utility.assert_kwarg(args, "box")
    local particle = utility.assert_kwarg(args, "particle")
    local pore_width = utility.assert_kwarg(args, "pore_width")
    local sigma = args.sigma or 1
    local epsilon = args.epsilon or 1
    local wetting = args.wetting or 1
    local cutoff = args.cutoff or 2.5
    local smoothing = args.smoothing or 0.005

    local potential = mdsim.potentials.external.planar_wall({
        -- place two walls perpendicular to the x-axis
        surface_normal = {{-1, 0, 0}, {1, 0, 0}}
        -- add a space of .5σ between the nominal pore space and each of the walls
      , offset = {(pore_width + 1) / 2, (pore_width + 1) / 2}
        -- wall interaction parameters
      , sigma = sigma, epsilon = epsilon, wetting = wetting
      , cutoff = cutoff,  smoothing = smoothing
      , species = (type(sigma) == 'number') and 1 or nil
    })

    -- register computation of wall forces
    local force = mdsim.forces.external({
        box = box, particle = particle, potential = potential
    })

    return force, potential
end

function M.create_pair_force(args)
    local box = utility.assert_kwarg(args, "box")
    local particle = utility.assert_kwarg(args, "particle")
    local cutoff = args.cutoff
    local smoothing = args.smoothing or 0.005

    -- define Lennard-Jones pair potential,
    local potential = mdsim.potentials.pair.lennard_jones()
    -- optionally, apply interaction cutoff
    if cutoff then
        -- use smooth truncation
        if smoothing > 0 then
            potential = potential:truncate({"smooth_r4", cutoff = cutoff, h = smoothing})
        else
            potential = potential:truncate({cutoff = cutoff})
        end
    end

    -- register computation of pair forces
    local force = mdsim.forces.pair({
        box = box, particle = particle, potential = potential
    })

    return force, potential
end