This code computes the density profile, surface excess, and wall tension, for a fluid of pure dissipative particle dynamics (DPD) beads against a wall, using a mean-field density functional theory (DFT).
The codes are:
wall_dft_single.py: solve a single state point;wall_dft_awall_scan.py: solve for an array of wall repulsion amplitudes;wall_dft_rhob_scan.py: solve for an array of bulk densities;wall_dft_zero.py: solve for zero surface excess or wall tension;wall_dft_minim.py: solve for minimum perturbation to bulk;wallDFT.py: underlying python module implementing functionality.
All the python scripts can be run without command line parameters, using the built-in defaults. They were developed to support a forthcoming publication in the area.
Any other scripts in this repository are experimental; you are welcome to investigate them and use them, but they are not particularly well documented ! See:
lle_interface.py: solve liquid-liquid phase equilibria.
Consider a fluid of particles at a non-uniform density ρ(r) interacting with a pair potential U(r), in the presence of an external potential Uext(r). A mean-field density functional free energy for this system is
- ∫ ρ(r) [ln ρ(r) − 1] + ∫ ρ(r) Uext(r) + ½ ∫ ρ(r) ρ(r') U(|r−r'|) .
If we consider the case where Uext(r) = Uwall(z) represents a wall in the plane normal to the z-direction, then this reduces to
- ∫ ρ(z) [ln ρ(z) − 1] + ∫ ρ(z) Uwall(z) + ½ ∫ ρ(z) ρ(z') U(z−z') ,
where U(z) is a partial integral of U(r) corresponding to the interaction between two parallel sheets of particles at unit density, separated by a distance z.
The corresponding grand potential (per unit area) is given by subtracting μ ∫ ρ(z) from the above expression, where μ is the chemical potential. By eliminating μ in favour of the bulk density ρb, and minimising the grand potential, one arrives at the following condition for the density profile in the wall potential,
- ρ(z) = ρb exp[ − Uwall(z) − ∫ dz' Δρ(z') U(z−z') ] ,
where Δρ(z) = ρ(z) − ρb. This has to be solved self-consistently, and this is what the code in this repository does. As input, one needs to specify the bulk density and the two potential functions U(z) and Uwall(z). The method is to start with the initial guess ρ(z) = ρb exp[ − Uwall(z) ], and iterate the above, mixing in a fraction of the new solution at each iteration step to assure convergence (Picard method). The integral in the exponential here can be evaluated as a convolution, using a standard numerical library routine.
Given a solution ρ(z), one can compute the wall tension γ and the surface excess Γ. The former is just the excess grand potential per unit area, and I define the latter as ∫ dz Δρ(z) outside the hard wall (see below) (z ≥ 0). The grand potential per unit area Ω / A = ∫ dz ρ(z) ωN(z) where ωN(z) = − 1 − ½ ∫ dz' ρ(z') U(z−z'). The bulk grand potential per unit volume needed to calculate γ is ωb = − p where p = ρb + ½ ρb2 ∫ dz U(z) is the pressure.
It follows from classical thermodynamics that dγ = − Γ dμ (Gibbs isotherm). It is also true that dp = ρ dμ (Gibbs-Duhem relation). Hence the dγ/dp = − Γ / ρ , and this can be verified numerically. Note that Γ / ρ, which can be positive or negative, can be interpreted as an adsorption length.
To quantify the perturbation in the bulk caused by the wall, the codes also report the integral of |Δρ(z)| outside of the wall potential (z > 1 in the wall models defined below).
In DPD the pair potential is
- U(r) = A (1−r)2 / 2
for r < 1, and U(r) = 0 for r > 1. The length scale is set by the range rc = 1, and the energy scale is set by the choice kBT = 1. The DPD model is then characterised by the bulk density ρb and repulsion amplitude A. For example, the standard DPD water model has ρb = 3 and A = 25.
This choice of pair potential gives the following partial integral to be used in the wall DFT calculations,
- U(z) = π A (1 − z)3 (1 + 3z) / 12
for 0 ≤ z ≤ 1, together with U(z) = 0 for z > 1, and U(z) = U(|z|) for z < 0.
I consider two different wall models, a 'standard' one where
- Uwall(z) = Awall (1 − z)2 / 2
for 0 ≤ z ≤ 1, together with Uwall(z) = 0 for z > 1, and Uwall(z) = ∞ for z < 0 (ie a hard repulsive barrier at z = 0); and a minimally perturbative 'continuum' wall model with
- Uwall(z) = π Awall ρwall (1 − z)4 (2 + 3z) / 60
for 0 ≤ z ≤ 1, together with similarly Uwall(z) = 0 for z > 1, and Uwall(z) = ∞ for z < 0. Note that for z > 0 the corresponding wall force − ∂Uwall / ∂z = ρb U(z) for U(z) given above.
One can show that in the latter case, for Awall ρwall = A ρb, the mean-field DFT density profile is flat up to the wall, with ρ(z) = ρb for z > 0. Under these conditions, the wall is minimally perturbative in the sense that Δρ(z) = 0 outside the wall. Trivially, the surface excess Γ = 0. One can show that the wall tension calculated from the excess mean-field grand potential per unit area in this case is γ = π A ρb2 / 240.
The continuum wall model derives its name because it also corresponds to a DPD particle interacting with a continuum of wall particles at a density ρwall in the z < 0 half-space. This is a model first introduced by Goicochea and Alarcón in J. Chem. Phys. 134, 014703 (2011).
For early work on minimally-pertubative wall models in a more generic setting see Henderson in Mol. Phys. 74, 1125 (1991). I am grateful to Bob Evans for drawing my attention to this.
The code lle_interface.py calculates the density profiles and
interfacial properties (surface excesses, surface tension) for
coexisting phases in a binary DPD mixture. The methods are similar to
those above, except there is no wall potential and there are two
density fields. The trickiest part is to solve for the phase
equilibrium. This is done assuming NPT, where the pressure is set by
specifying a baseline (default ρ = 3, A = 25). For fixed pressure
there are two coexisting phase compositions x1 and
x1. These are obtained by solving numerically for
equality of chemical potentials in the coexisting phases. The initial
part of the code computes these compositions, and the associated
chemical potentials of the two components (the same in both phases).
These are then used in the density profile calculation.
This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.
This program is copyright © 2024 Patrick B Warren (STFC).
Send email to patrick.warren{at}stfc.ac.uk.