The effect of soft repulsive interactions on the diffusion of particles in quasi-one-dimensional channels: A hopping time approach

TL;DR

This study develops a transition state theory to predict hopping times for particles in quasi-one-dimensional channels, revealing how soft repulsive interactions influence diffusion crossover and how channel width affects particle dynamics.

Contribution

It introduces a quantitative transition state theory approach to predict hopping times in confined soft particle systems, accounting for particle softness and channel width effects.

Findings

Free energy barrier peaks at intermediate particle softness.

Hopping times are influenced by channel width and particle softness.

Maximum free energy barrier shifts with channel narrowing.

Abstract

Fluids confined to quasi-one-dimensional channels exhibit a dynamic crossover from single file diffusion to normal diffusion as the channel becomes wide enough for particles to hop past each other. In the crossover regime, where hopping events are rare, the diffusion coefficient in the long time limit can be related to a hopping time that measures the average time it takes a particle to escape the local cage formed by its neighbours. In this work, we show that a transition state theory that calculates the free energy barrier for two particles attempting to pass each other in the small system isobaric ensemble is able to quantitatively predict the hopping time in a system of two-dimensional soft repulsive discs [$U(r_{ij})=(\sigma/r_{ij})^{\alpha}$] confined to a hard walled channel over a range of channel radii and degrees of particle softness measured in terms of $1/\alpha$. The free energy barrier exhibits a maximum at intermediate values of $\alpha$ that moves to smaller values of $1/\alpha$ (harder particles) as the channel becomes narrower. However, the presence of the maximum is only observed in the hopping times for wide channels because the interaction potential dependence of the kinetic prefactor plays an increasingly important role for narrower channels. We also begin to explore how our transition state theory approach can be used to optimize and control dynamics in confined quasi-one-dimensional fluids.

Tables

Table 1. Table 1: Values of 1 / α 1 𝛼 1/\alpha and β ​ Δ ​ G ∗ 𝛽 Δ superscript 𝐺 \beta\Delta G^{*} at the maximum in the hopping free energy barrier.

	Global Optimization		Quadratic Fit
$R$	$1/\alpha$	$\beta \Delta G^{\ast }$	$1/\alpha$	$\beta \Delta G^{\ast }$
1.01	0.0107	10.095	0.0104	10.084
1.03	0.0123	7.978	0.0165	7.947
1.05	0.0215	6.957	0.0292	6.950
1.07	0.0429	6.294	0.0425	6.289
1.10	0.0585	5.603	0.0626	5.601

Equations

$$⟨\Delta x_t^2⟩=2F_x{t}^{1/2}\text{mbox},$$

$$D_x\propto \frac{1}{{\tau }_{hop}^{1/2}}\text{mbox},$$

$$U(r_{ij})=ϵ{\left(\frac{\sigma }{r_{ij}}\right)}^{\alpha },5\le \alpha \le 100\text{mbox},$$

$$U_W(y_i)=\{\begin{array}{l}0\text{if}|y_i|\le R-\sigma /2,\\ \infty \text{if}|y_i|>R-\sigma /2,\end{array}$$

$$k_{TST}=\frac{{⟨|v_c^{\ast }|⟩}_{eq}}{2}{P}_0(x_c^{\ast })\text{mbox},$$

$${\tau }_{hop}(TST)=\frac{1}{k_{TST}}\text{mbox}.$$

$$\Delta ={\int }_L{Q}_{n,v,T}^{\ast }{e}^{-\beta p_{l}2RL}2RdL\text{mbox},$$

$${\Delta }_{x_c^{\prime }}d{x}_c={\int }_L{Q}_{n,v,T}^{\ast }{e}^{-\beta p_{l}2RL}\delta (x_c^{\prime }-x_c)2RdLd{x}_c\text{mbox},$$

$$\Delta ={\int }_{x_c}{\Delta }_{x_c^{\prime }}d{x}_c\text{mbox}.$$

$$P_0(x_c^{\prime })d{x}_c=\frac{{\Delta }_{x_c^{\prime }}d{x}_c}{\Delta }\text{mbox},$$

$${\int }_0^{\infty }{P}_0(x_c)d{x}_c=1\text{mbox}.$$

$$\beta \Delta G^{\ast }=-\ln P^{\ast },$$

$$P^{\ast }={\int }_0^{x_c^{\ast }}{P}_0(x_c)d{x}_c\text{mbox},$$

$$\begin{array}{rl}\text{mbox}acc& \text{mbox}(old\to \text{mbox}new)=\text{mbox}min(1,\exp \{-\beta [U^{\text{mbox}new}-U^{\text{mbox}old}]\\ & -\beta p_{l}2R[L^{\text{mbox}new}-L^{\text{mbox}old}]+(n-1)\ln [L^{\text{mbox}new}/L^{\text{mbox}old}]\})\text{mbox}.\end{array}$$

$$⟨|v_c^{\ast }|⟩=⟨\frac{|x_c(t+\Delta t)-x_c(t)|}{\Delta t}⟩\text{mbox},$$

$$Q_{1d}=\frac{2}{2!{\Lambda }^2}{\int }_0^{2R_r}d{r}_1{\int }_0^{2R_r-r_1}{e}^{-\beta U(r_{ij})}d{r}_{ij}\text{mbox},$$

$$\beta \Delta F_{1d}^{\ast }=-\ln \left(Q_{1d}/Q_{ig}\right)\text{mbox},$$

Full text

I Introduction

The movement of fluids through narrow quasi-one-dimensional channels appears in a variety of natural and engineered systems Kärger (2015); Gubbins et al. (2010) such as zeolites, Kukla et al. (1996); Kärger et al. (1992); Kumar (2014) carbon nanotubes, Das et al. (2010); Valiullin and Kärger (2010); Chen et al. (2010) metallic organic frameworks, Salles et al. (2011); Jobic (2016) confined colloids, Wei, Bechinger, and Leiderer (2000); Lutz, Kollmann, and Bechinger (2004) and in the transport of water and ions through membranes. Hodgkins and Keynes (1955); Finkelstein (1987) It also plays a fundamental role in the functioning of nano- and micro- fluidic devices Siems et al. (2012); Locatelli et al. (2016); Yang et al. (2010) and provides a basis for separating mixtures. Adhangale and Keffer (2003); Ball et al. (2009); Wanasundara, Spiteri, and Bowles (2012) When the particles are subjected to stochastic forces Levitt (1973); Mon and Percus (2003) or a Brownian background Percus (1974) and the channel is narrow enough to prevent passing, the system exhibits a form of sub-diffusion known as single file diffusion (SFD), where the mean squared displacement (MSD) increases as the square root of time in long time limit. The MSD can then be described by an Einstein-like relation,

[TABLE]

where $F_{x}$ is the mobility factor. Hahn and Kärger (1998); Lin et al. (2005) As the channel becomes wider, the particles are able to pass each other, and the system exhibits a dynamical transition or crossover from SFD to normal diffusion, Kumar (2014); Chen et al. (2010); Lucena et al. (2012); Herrera-Velarde, Pérez-Angel, and Castañeda-Priego (2016); Ooshida, Goto, and Otsuki (2018) where the MSD increases linearly in time in the long time limit. However, in the crossover regime, hopping events are rare because particles must overcome an entropic barrier caused by a restriction in configuration space as they attempt to pass. If the particles remain caged between their neighbours long enough to perform SFD before escaping, the MSD follows Eq. (1) at intermediate times before crossing over to normal diffusion at long times. Wanasundara, Spiteri, and Bowles (2014)

Percus and Mon Mon and Percus (2002) developed a simple theory that connects the intermediate time behaviour to the long time normal diffusion coefficient, $D_{x}$, as,

[TABLE]

where $\tau_{hop}$ is a phenomenological hopping time that measures the average time it takes for a particle to escape its cage. This hopping time approach to understanding diffusion in confined environments is attractive because all the factors that influence particle mobility, such as density, temperature, and particle–particle and particle–wall interactions, are captured in a single parameter that can be measured directly in simulation. A variety of theoretical and computational methods can be used to study the hopping time. For example, projection operator techniques have been used to map the diffusion of the system into a modified one-dimensional Fick–Jacobs equation Kalinay and Percus (2005); Kalinay (2007); Kalinay and Percus (2008); Mon (2008a, b, 2009) to take advantage of the quasi-one-dimensional nature of the channels. However, the activated nature of the hopping process also suggests $\tau_{hop}$ can be obtained theoretically using barrier crossing methods, and a simple transition state theory has been shown to provide the correct scaling exponent for the hopping time as a function of channel radius Bowles, Mon, and Percus (2004) for two-dimensional hard discs. It was recently shown that a transition state theory (TST) qualitatively predicts the hopping time for the same system, Ahmadi and Bowles (2017) where the height of the free energy barrier for two particles attempting to pass was obtained using the small system isobaric-isothermal ensemble. Corti and Gerardo (1998); Corti (2001)

The goals of the current paper are to show that TST provides quantitative measures of the hopping time for particles confined to quasi-one-dimensional channels when both the barrier and prefactor are calculated and to demonstrate how the hopping time approach can be used to study the role of particle–particle interaction on diffusion in confined fluids. To achieve these goals, we compare the hopping times obtained using TST with those measured directly from the simulation of large number of particles for a system of two-dimensional soft repulsive discs confined to a hard wall channel for a range of channel radii. We also examine the effect of particle softness on the hopping time, and hence diffusion, by varying the magnitude of the exponent associated with the repulsive particle-particle potential. Intuitively, we might expect the hard particle system to exhibit the longest hopping times because they restrict the configuration space of the transition state to the greatest degree. Making the potential softer should make it easier for the particles to reach the transition state. Our results show that the hopping barrier goes through a maximum as the particles become softer because of a competition between the energetic and entropic contributions to the free energy of the transition state. However, the hopping times themselves only appear to exhibit the maximum as a function of particle softness for wider channels because the potential dependence of the kinetic prefactor in the TST theory becomes increasingly important for the narrower channels. We also begin to explore how our TST method can be used to optimize the dynamic properties of single file systems.

The remainder of the paper is organized as follows. Section II outlines the model system and the methods, including outlining the general TST approach to hopping times, the hopping free energy barrier calculations, the kinetic prefactor calculations, the hopping times measurements in large systems, and the hopping time optimization. Section III describes our results and discussion, and Section IV summarizes our conclusions.

II Model and Methods

II.1 Model

We study a system of two-dimensional repulsive discs confined to a narrow two-dimensional channel that extends longitudinally along the $x$-axis and has a radius, $R$, along the $y$-axis, where the origin of the coordinate system is located at the center of the channel. The particles interact through an inverse-power-law potential,

[TABLE]

where $r_{ij}=|r_{i}-r_{j}|$ is the distance between particles $i$ and $j$, $\alpha$ is the power-law exponent describing the repulsive potential, $\sigma$ is the particle diameter, $\epsilon$ is the interaction strength, and the interaction is cut off at $r_{c}=2.5\sigma$. We use the quantity $1/\alpha$ as a measure of particle softness because it tends to zero as $\alpha\rightarrow\infty$, describing the hard sphere potential, and increases as the potential allows greater particle overlap. The particles interact with the walls of the channel as hard particles, so the particle-wall interaction is given by,

[TABLE]

where $y_{i}$ is the $y$-coordinate of particle $i$. Our simulations are carried out using reduced units. We also study the equivalent two-dimensional hard sphere system under the same conditions used to study the soft particles, so we can make a direct comparison.

II.2 Transition State Theory for Hopping times

The particles in a fluid confined to a quasi-one-dimensional channel are arranged in a single file so that each particle is caged by its (first) nearest neighbour in each direction along the channel. In order to diffuse normally in the long time limit, the particles must hop past each other. This hopping is an activated process because the excluded volume interactions of the particles and the wall restrict the configuration space available to the particles as they pass, giving rise to an entropic free energy barrier. Transition state theory, which assumes that all the trajectories that cross through the transition state are reactive, i.e., recrossing events are ignored, provides an upper estimate for the transition rate in an activated process given by, Bennet (1977); Chandler (1978); Ruiz-Montero, Frenkel, and Brey (1997); ten Wolde, Ruiz-Montero, and Frenkel (1999)

[TABLE]

where $P_{0}(x_{c}^{*})$ is the probability density of the system located at the transition state, denoted $x_{c}^{*}$, and $v_{c}^{*}$ is the velocity of the system along the reaction coordinate as it passes through the transition state. The TST hopping time is then given by,

[TABLE]

The free energy barrier associated with a particle hopping past one of its neighbours to escape its cage is calculated using the small system isobaric-isothermal ($n,p,T$) ensemble developed by Corti et al. Corti and Gerardo (1998); Corti (2001) In general, the method involves considering the properties of a small system of $n$ particles confined to a small volume, $v$, immersed in a larger system of $N-n$ particles with volume $V-v$, both at a fixed temperature, $T$. One of the $n$ small system particles, denoted the shell particle, is constrained to a region $dv$ at the boundary of the small system, defining the small system volume. This avoids an over counting of configurations associated with the fluctuation of $v$ when integrating over the degrees of freedom of the large $N-n$ system to yield a small system constant pressure ensemble that does not rely on the use of a system size dependent volume scale.

Here, we provide a brief description of the small system isobaric-isothermal ensemble in the context of its application to hopping times. Ahmadi and Bowles (2017) Figure 1(a) shows the system, consisting of $n=2$ particles, at a constant external longitudinal pressure, $p_{l}$, and fixed $T$. The small system has a length, $L$, and radius, $R$, with the center located at $r_{0}$. The shell particle (particle one) is confined to a region $2\,R\,dL$ at $x_{1}=L/2>0$ to define the volume. The partition function for the small system given by,

[TABLE]

where $Q^{*}_{n,v,T}$ is the canonical partition function of the small $n,v,T$ system with one shell particle and $\beta=1/k_{B}T$, where $k_{B}$ is the Boltzmann constant. In obtaining Eq. (7), it has been assumed that the interactions between the particles in the small and large systems are negligible.

Particle two in the system represents the caged particle, and we define a reaction coordinate for the hopping process as the axial separation between the particles, $x_{c}=x_{1}-x_{2}$, so that the transition state occurs when they are side by side in the channel, i.e., with $x_{c}=0$, as shown in Fig. 1(b). The reaction coordinate partition function is then given by,

[TABLE]

where $\delta(\cdot)$ is the Dirac delta function. Integrating Eq. (8) yields the partition function,

[TABLE]

The probability of finding the caged particle at $x_{c}^{\prime}$ is then,

[TABLE]

where the probability density, $P_{0}(x_{c}^{\prime})$, is normalized so that

[TABLE]

The Gibbs free energy barrier for hopping, which is defined as the work required to bring the caged particle from anywhere in the cage to the transition state is given by,

[TABLE]

where

[TABLE]

and the range of the reaction coordinate, $x_{c}=0$ to $x_{c}^{*}$, is the transition state region. It should also be noted that we select $x_{c}^{*}$ to be small compared to the size of the motion required to take the particle from one cage to the other so that configurations that have little or no chance of crossing the barrier in the small time limit are excluded from the transition state.

The free energy barrier for hopping is obtained from a Monte Carlo (MC) simulation as follows: Particle one, the shell particle, is located at $L/2$, and particle two is placed within the cell between $-L/2$ and $L/2$. For each MC move, a particle is randomly selected and moved by $\delta x$ and $\delta y$, up to a maximum displacement of $|\Delta x|=0.06\sigma$ and $|\Delta y|=0.12\sigma$. The move is immediately rejected if the trial displacement causes the particle to overlap with the hard wall or if particle two is moved outside the cage. If particle one is moved, the cell length is changed by $2\delta x$ and the position of particle two is re-scaled to ensure it remains within the simulation cell. Because the position of the shell particle defines the system volume, the position of the shell molecule must be positive during the simulation to ensure $V>0$. The MC probability of accepting the trial move from the old to the new configurations is given by Corti (2002); Ahmadi and Bowles (2017)

[TABLE]

For each system, the value of $\Delta G$ is obtained as an average over 20 independent runs. For each run, $2\times 10^{6}$ MC cycles are used to reach equilibrium, and the free energy barrier is calculated over the next $8\times 10^{8}$ MC cycles, where an MC cycle consists of $n=2$ MC moves. We sampled configurations every 1000 MC cycles to ensure that they are not correlated. The probability is calculated by building a normalized histogram of all configurations along the reaction coordinate, where we have used bin sizes of $\Delta x_{c}=x_{c}^{*}=0.05\sigma$, and the free energy barrier is calculated using the probabilities. The error bars are calculated as the standard deviation of $\Delta G^{*}$ over the 20 runs. Our simulations consider systems with $\alpha$ in the range 5–100 and $R/\sigma$ in the range 1.01–1.10 under conditions with $\beta p_{l}=1$ and $\beta=1$. The probability density at the transition state, $P_{0}(x_{c}^{*})$, is used in Eq. (5) to calculate the transition rate, and the hopping time is obtained by dividing the transition state probability, $P^{*}$, by $\Delta x_{c}$.

The average velocity at which the system crosses the barrier is defined by

[TABLE]

where $x_{c}(t+\Delta t)$ and $x_{c}(t)$ are the values of the reaction coordinate close to the transition state. To obtain $\left<|v_{c}^{*}|\right>$, we performed simulations in the canonical ($n,V,T$) ensemble using the standard Metropolis MC scheme Frenkel and Smit (2002) to move the particles as a simple approximation to Brownian motion. Patti and Cuetos (2012) Each MC move involves randomly selecting one of the two particles and moving it randomly by a step $\delta x$ and $\delta y$ up to a maximum displacement of $|\Delta x|=0.06\sigma$ and $|\Delta y|=0.12\sigma$. The move was rejected if the trial displacement caused any of two particles to overlap with the hard wall; otherwise the move was accepted according to the standard Metropolis MC probabilities. A Monte Carlo cycle involves $n=2$ attempted MC moves and defines the unit of time. We then used $\Delta t=1$ in Eq. (15). The starting configuration for each run, with $x_{c}(t=0)<0.05$, was taken from the free energy simulations described in the previous section, and average quantities were calculated over 1000 distinct initial configurations.

II.3 Hopping Time Measurement

In order to test the predictions of our TST, we directly measure the hopping time in a series of canonical ($N,V,T$) simulations consisting of $N=500$ particles confined to a channel of radius, $R$, and length, $L$. Periodic boundary conditions are used in the longitudinal direction, and the particle dynamics is the same as that described in the previous section. For each MC move, a particle is randomly selected and moved in a random direction by a step $\delta x$ and $\delta y$ up to a maximum displacement of $|\Delta x|=0.06\sigma$ and $|\Delta y|=0.12\sigma$. An MC cycle involves $N$ MC move attempts and defines the unit of time, $t$, in the simulation. Particles are initially placed uniformly along the channel but randomly across the width of the channel. Then $3\times 10^{7}$ MC cycles are performed to equilibrate the system before collecting data over the next $5\times 10^{7}\;\textup{to}\;3\times 10^{8}$ MC cycles in order to be long enough for the average hopping time to converge. After equilibrium is achieved, the cage for each particle is defined by their immediate right and left neighbours, and the initial hopping time is set to zero for all the particles. The number of MC cycles that it takes for each particle to escape their cage is counted as their hopping time. After each hopping event, the hopping time for the particle is reset to zero and the new cages are identified. The hopping time, $\tau_{hop}$, is calculated as an average over all hopping times recorded for all particles. We also performed simulations with $N=100-800$ to check for system size effects on the hopping time calculations, and these additional simulations confirm $N=500$ is large enough to account for system size effects.

The appropriate linear density, $\rho_{l}=N/L$ for each system ($R,\alpha$) studied is determined by performing an $N,p_{l},T$ simulation, where $N=500$, $p_{l}/k_{B}T=1$. In the limit of large system sizes, the standard constant pressure simulation method and the shell particle method yield the same results, Corti (2002) so we continued to use the shell particle method. An MC trial move follows that outlined for the free energy calculations (Eq. (14)), with the maximum step size of $|\Delta x|=0.06\sigma$ and $|\Delta y|=0.12\sigma$, and results in an MC acceptance ratio of $85-92\%$. The average $\rho_{l}$ for the system is measured over $5\times 10^{7}$ MC cycles, after the system reaches equilibrium over $10^{7}$ MC cycles, and data are sampled every $10\,000$ MC moves to ensure that they are not correlated. Figure 2 shows that the densities increase as the channels expand but remain relatively constant as a function of $\alpha$.

III Results and Discussion

III.1 Hopping times

The hopping times obtained directly from the large system simulations provide a measure by which we can qualitatively examine the influence of particle softness on diffusion. Figure 3 shows that, as soft particles become harder (decreasing $1/\alpha$), $\tau_{hop}$ initially increases, doing so more rapidly for the narrower channels. This is consistent with our expectation that particle hopping in a channel should be more difficult for harder particles. We also know that the hopping time must diverge in the limits $R\rightarrow\sigma$ and $1/\alpha\rightarrow 0$ because perfectly hard particles are unable to pass each other when the channel radius decreases below the passing threshold, causing the fluid to undergo a dynamic transition from normal diffusion to SFD. However, we also find that the hopping times exhibit an unexpected maximum, suggesting that there are cases where the softer particles diffuse more slowly than the harder particles. The maximum, which is directly observable for the wider channels, sharpens and moves to lower values of $1/\alpha$ as the channels narrow until it is no longer visible within the range of $1/\alpha$ studied. In these narrow channel cases, the hopping times grow larger than those measured for the perfect hard discs model, suggesting that the maximum may still occur at lower values of $1/\alpha$ than were studied.

We now examine the properties of the free energy barrier to hopping. Figure 4 shows the Gibbs free energy as a function of separation of the two particles along the reaction coordinate for systems with different values of $\alpha$. At large $x_{c}$, the interaction between the particles is negligible, and they are both able to move freely across the entire diameter of the channel. However, constraining the particles to a particular $x_{c}$ puts a lower limit on the fluctuations available to the cell in the $n,p_{l},T$ ensemble. The linear increase in the free energy with increasing $x_{c}$ is related to the work performed against the external pressure as the particles are separated, forcing the cage to increase in volume. The same effect is observed for two ideal gas particles. Ahmadi and Bowles (2017) The free energy goes through a minimum at $x_{c}\approx\sigma$ and begins to increase as the repulsive interaction between the particles begins to restrict their available configuration space. The minimum is sharp for the harder particles, but with decreasing $\alpha$, the softer, longer range interactions lead to a shallower and softer minimum as well as a slower increase toward the maximum near the transition state at $x_{c}=0$. When the channels are very narrow, the maximum is located very close to the geometrically defined transition state, but this is not always the case. Ahmadi and Bowles (2017) The particles studied here are circular, so the greatest degree of configurational restriction occurs when $x_{c}=0$, but as the channel becomes wider, the influence of the pressure-volume work become significant, and this can lead to the appearance of the maximum before the geometric transition state is reached. Nevertheless, it is important to recognize that it is the probability of being at the geometric transition state that contributes to the hopping time because the particles are unable to diffuse in the long time limit unless they exchange positions along the channel.

Figure 5 shows the height of the free energy barrier for hopping exhibits the same general trends as $\ln\tau_{hop}(sim)$ measured from the simulation, and moreover we are able to directly identify the location of the free energy maximum for all channel radii. Measurements of the prefactor term (see Fig. 6) also account for the disappearance of the maxima in the hopping times for narrow channels. For wider channels, the prefactor is relatively independent of the interaction potential, so the hopping times essentially follow the behaviour of the free energy. However, for narrower channels, the prefactor term becomes strongly potential dependent, and its increase outweighs the effect of the decreasing free energy. This highlights the importance of understanding the properties of the prefactor even in a regime where the free energy barriers are high and are expected to dominate the overall behaviour. Figure 6 also shows that the prefactor is relatively independent of channel width for the softer particles but is highly channel width dependent for the harder particles.

Figure 7 shows that our TST approach generally predicts the hopping times $\tau_{hop}(TST)$ to within a factor of two, underestimating the times for wide channels and overestimating them for narrow channels. Transition state theory should underestimate the hopping time because it assumes that all trajectories that cross through the transition state lead to a hopping event. Improvements that account for trajectory recrossing, where the particle re-enters its original cage, correct terms associated with the kinetic prefactor and lead to longer hopping times. However, our two particle approach makes some additional assumptions that may lead to further variations from the measured value. The construction of the two particle model, through the mapping of the large canonical ensemble onto the small system $n,p,T$ ensemble is formally exact, but when we perform our calculations, we neglect the interactions between the system and its surroundings, leading to errors in the free energy calculations. In particular, we expect the free energy calculations to be less accurate for the systems with softer, longer range potentials because we neglect the interactions outside of the cell. As the channels become narrower, the fluid structure becomes more single file, reducing the role of long range interactions and improving the accuracy of the free energy calculations. It is important to note that the prefactor calculations are also performed using just two particles; this ignores the role second neighbours may play in blocking the escape of a particle from its cage. Nevertheless, the simplicity of our two particle analysis combined with its ability to capture the general trends in the hopping times for the particles suggests that it could be useful in determining the factors that influence diffusion in highly confined fluids.

III.2 Transition State Partition Function

To gain insight into the origin of the maximum in the barrier height as a function of particle softness, we study the ensemble of states associated with the transition state, which for the current system is described by the partition function for two particles confined to a one-dimensional line (see Fig. 8). The one-dimensional transition state partition function for the system can be written as,

[TABLE]

where $\Lambda=(2mk_{B}T/h^{2})^{1/2}$ is the thermal de Broglie wavelength for a particle with mass $m$, $h$ is Planck’s constant, $R_{r}=R-\sigma/2$ is the reduced channel radius accessible to the centers of the particles due to the hard wall interaction, $r_{ij}=r_{2}-r_{1}>0$ is the particle separation, $U(r_{ij})$ is the interaction potential given by Eq. (3), and the factor of two accounts for the need to consider both possible particle orders on the line. The Helmholtz free energy for the transition state, relative to two ideal gas particles, is then given by,

[TABLE]

where $Q_{ig}=(1/2!\Lambda^{2})(2R)^{2}$ is the ideal gas partition function, noting that the ideal gas particles have no diameter and hence are able to sample the entire width of the channel. The energetic and entropic contributions to the free energy, $F=E-TS$, can then be obtained using the usual canonical ensemble expressions, $\left<E\right>=k_{B}T^{2}(\partial\ln Q/\partial T)_{N,V}$ and $S=k_{B}\ln Q+\left<E\right>/T$, respectively.

Figure 9 shows that, despite its simplicity, the one-dimensional transition state partition function captures the key features of the effect of the interaction potential on the hopping free energy barrier, including the presence of the maximum and how it evolves as a function of channel radius. However, it is important to recognize that the transition state partition function cannot be used directly in the TST expression to yield the hopping time because it does not give the probability of finding the system at the top of the barrier, but it clearly contains the relevant details on how the properties of the transition state vary as the system parameters vary.

The model also reveals that the maximum in the height of the hopping barrier as a function of particle softness results from a competition between the energy and entropy in the transition state as the interaction potential changes (see Fig. 10). The barrier for the hard disc system is entirely entropic in nature, with the excluded volume interaction between the particles and the wall restricting the accessible configurations space for the two particles. As the particles become softer, both the energy and the entropy relative to the ideal gas begin to increase because the particles can effectively overlap, but for small $1/\alpha$, the energetic cost increases more rapidly than the entropy, causing an overall increase in the height of the barrier. For larger values of $1/\alpha$ the situation is reversed. The energy essentially plateaus, and the entropy increasingly dominates, leading to an overall decrease in the hopping free energy barrier as the particles continue to become softer.

III.3 Optimization

Using the hopping time approach to study diffusion in quasi-one-dimensional systems is attractive because all the molecular details that influence mobility are contained within a single parameter, $\tau_{hop}$, suggesting it could be a useful target for the optimization of diffusion in these systems. Our TST then provides a simple and efficient way to evaluate the $\tau_{hop}$ as part of an optimization scheme.

To demonstrate the general principle, we use the pythOPT global optimization problem-solving software environment Voss (2016) to search for the value of the particle softness $(1/\alpha)$ that maximizes the hopping barrier for a given channel radius. This process corresponds to finding the maxima in Fig. 5 and approximately identifies the conditions under which diffusion is slowest. The pythOPT environment offers multiple global optimization solvers as well as a suite of benchmark problems and routines for performance analysis. The solver used in this study was the Guaranteed Convergence Particle Swarm Optimization (GCPSO). van den Bergh and Engelbrecht (2002) At each step in the search, the MC free energy calculation outlined in Section II is used to obtain the figure of merit. This obviously ignores the effects associated with the prefactor, but they can be included in the current analysis in a straightforward way. CGPSO simulations used 50 swarm particles and one million function evaluations. Each DG simulation took approximately five minutes to execute, and 40 simulations were used for each function evaluation. The optimization for each radius required approximately five CPU-days.

To compare the optimization results with estimates obtained from our systematic set of simulations, we fit a quadratic polynomial to the free energy as a function of $1/\alpha$ to the data in Fig. 5 for each channel radius for the five closest points the free energy maximum. Table 1 shows that there is excellent agreement between the maxima obtained from pythOPT and the systematic study for $\beta\Delta G^{*}$. As can be expected from the sensitivity of optimization problems Heath (2018), the agreement is not as good for the location $1/\alpha$ of the maxima, but the values are still close given the large degree of fluctuation in our free energy estimates. For example, for $R=1.01$, the fluctuation in $\beta\Delta G^{*}$ is on the order of $0.20\ kT$, and this could include values of $1/\alpha$ in the range $0-0.17$.

The example given here is relatively straightforward because we only vary one parameter and there is a well defined maximum in the free energy surface. However, it demonstrates the principle of using the TST approach, either by simply calculating the free energy or the full hopping time through the added calculation of the prefactor term, in a search for an optimal condition to control diffusion. In general, an optimization for the purposes of engineering design would involve the variation of a variety of parameters relating to channel diameter and particle–particle and particle–wall interactions.

IV Conclusions

The hopping time, $\tau_{hop}$, is a phenomenological parameter that connects a local parameter, i.e., the time it takes for a particle to escape the cage of its neighbours, to the long time diffusion coefficient of the system. It contains information about the roles system parameters, such as the density and particle-particle and particle-wall interactions, play in the dynamics of quasi-one-dimensional fluids. We have shown that a transition state theory, where the free energy barrier for two particles attempting to pass is calculated in the small system isobaric-isothermal ensemble, provides quantitative predictions for $\tau_{hop}$ and is able to reveal interesting details regarding effect of inter-particle interaction in the dynamics of the system. This suggests that the hopping time TST approach could be an effective tool for the optimization and control of diffusion in nano- and micro-fluidic devices.

Acknowledgements.

We would like to thank Natural Sciences and Engineering Research Council of Canada (NSERC) for financial support. Computing resources and support were provided by Compute Canada and WestGrid.