  The problem of mass diffusion in layered systems has relevance to
applications in different scientific disciplines, e.g., chemistry, material
science, soil science, and biomedical engineering. The mathematical challenge
in these type of model systems is to match the solutions of the time-dependent
diffusion equation in each layer, such that the boundary conditions at the
interfaces between them are satisfied. As the number of layers increases, the
solutions may become increasingly complicated. Here, we describe an alternative
computational approach to multi-layer diffusion problems, which is based on the
description of the overdamped Brownian motion of particles via the underdamped
Langevin equation. In this approach, the probability distribution function is
computed from the statistics of an ensemble of independent single particle
trajectories. To allow for simulations of Langevin dynamics in layered systems,
the numerical integrator must be supplemented with algorithms for the
transitions across the discontinuous interfaces. Algorithms for three common
types of discontinuities are presented: (i) A discontinuity in the friction
coefficient, (ii) a semi-permeable membrane, and (iii) a step-function chemical
potential. The general case of an interface where all three discontinuities are
present (Kedem-Katchalsky boundary) is also discussed. We demonstrate the
validity and accuracy of the derived algorithms by considering a simple
two-layer model system and comparing the Langevin dynamics statistics with
analytical solutions and alternative computational results.
