# Slip and Stress From Low Strain-Rate Nonequilibrium Molecular Dynamics: The
Transient-Time Correlation Function Technique

Luca Maffioli  Department of Mathematics, School of Science, Computing and
Engineering Technologies, Swinburne University of Technology, PO Box 218,
Hawthorn, Victoria 3122, Australia  Edward R. Smith  Mechanical and Aerospace
Engineering, Brunel University London, Kingston Lane, Uxbridge UB8 3PH, United
Kingdom  James P. Ewen  Department of Mechanical Engineering, Imperial College
London, London SW7 2AZ, United Kingdom  Peter J. Daivis  School of Science,
RMIT University, GPO Bpx 2476, Victoria 3001, Australia  Daniele Dini
Department of Mechanical Engineering, Imperial College London, London SW7 2AZ,
United Kingdom  B. D. Todd<EMAIL_ADDRESS>Department of Mathematics,
School of Science, Computing and Engineering Technologies, Swinburne
University of Technology, PO Box 218, Hawthorn, Victoria 3122, Australia

###### Abstract

We derive the transient-time correlation function (TTCF) expression for the
computation of phase variables of inhomogenous confined atomistic fluids
undergoing boundary-driven planar shear (Couette) flow at constant pressure.
Using nonequilibrium molecular dynamics simulations, we then apply the TTCF
formalism to the computation of the shear stress and the slip velocity for
atomistic fluids at realistic low shear rates, in systems under constant
pressure and constant volume. We show that, compared to direct averaging of
multiple trajectories, the TTCF method dramatically improves the accuracy of
the results at low shear rates, and that it is suitable to investigate the
tribology and rheology of atomistically detailed confined fluids at realistic
flow rates.

††preprint: AIP/123-QED

##  I Introduction

Nonequilibrium molecular dynamics (NEMD) simulations have given unique
insights into nanoscale fluid behaviour in a range of applications from
tribologyEwen2018  to hydraulic fracturing and polymer processing
doi:10.1021/acs.jpclett.6b01684 ; Gartner2019 . On the other hand, NEMD
simulations are generally limited to external fields which are several orders
of magnitude larger than those encountered in experiments and applications.
For example, typical shear ratesdoi:10.1177/1350650117696181  in automotive
engines range from $10^{5}-10^{8}\;\text{s}^{-1}$, whilst those in polymer
processing and hydraulic fracturing are usually between
$10^{3}-10^{6}\;\text{s}^{-1}$ and $1-10^{3}\;\text{s}^{-1}$, respectively
Vlachopoulos2003 ; doi:10.1146/annurev-chembioeng-080615-033630 . Even using
massively-parallelized MD softwarePLIMPTON19951 ; THOMPSON2022108171  on
modern high performance computer (HPC) systems, the lowest accessible shear
rates by direct NEMDJadhao7952  are $>10^{5}\;\text{s}^{-1}$.

For atomic fluids, the viscosity remains Newtonian up to high shear
rates9943764 ; however, for molecular fluids, non-Newtonian shear thinning
behaviour can dominate the viscous response at shear rates well below those
accessible to direct NEMD simulationsJadhao7952 ; PhysRevLett.88.058302 .
Comparison of NEMD results with experimental viscosity measurements is thus
generally restricted to extrapolationJadhao7952 ; C7CP01895A  or time-
temperature superpositionJadhao7952 ; PhysRevLett.88.058302  methods.

One approach which allows the study of fluid rheology at realisitic shear
rates is the transient-time correlation function (TTCF) technique
PhysRevA.10.2461 ; Dufty1979 ; COHEN198317 ; PhysRevA.35.792 ;
PhysRevA.38.4142 , which is based on the time-correlation between the initial
rate of energy dissipation and the transient response of any arbitrary phase
variable after an external field is activated. TTCF has been applied to
investigate the rheology of a range of fluids at low shear rates, usually in
homogenous systems without confining wallsdoi:10.1080/08927020601026629 .
These studies have progressed from atomic fluidsPhysRevA.35.792 ;
PhysRevA.38.4142 ; doi:10.1080/00268970210137275  to molecular
fluidsdoi:10.1080/00268970902922625 , and even liquid metalsPhysRevB.78.184202
; doi:10.1063/1.2829872 . In addition to the viscosity, TTCF has also been
used to monitor the electrical conductivitydoi:10.1063/1.2035085 ; English2010
; doi:10.1080/00268976.2010.544263 , thermal
conductivitydoi:10.1080/08927020801930604 , colour
conductivityPhysRevE.77.027701 , and normal stress
differencesdoi:10.1080/08927020802575598 . The TTCF method has also been
extended from shear flow to elongational flowPhysRevE.56.6723 ;
PhysRevE.58.4587  as well as mixed shear and elongational
flowdoi:10.1063/1.3684753 . Relatively fewer studies have applied TTCF to
boundary-driven simulations of confined, inhomogeneous
systemsPhysRevB.72.172201 ; 22920110 ; doi:10.1080/08927022.2015.1049174 .

Boundary-driven NEMD was introduced by Ashurst and HooverPhysRevA.11.658  and
further developed by Bitsanis et al. doi:10.1063/1.453240  who included the
confining walls. In boundary-driven confined NEMD simulations, shear is
usually applied by moving the walls in opposite directions and the temperature
inside the channel is controlled by a thermostat applied either on the whole
system or only to the wall atoms, as occurs in experimentsPhysRevA.45.3706 ;
doi:10.1063/1.3450302 ; Yong2013 . This approach enables the study of new
dynamic behavior which emerges when fluids are strongly
confineddoi:10.1126/science.253.5026.1374 ; doi:10.1126/science.250.4982.792 ;
doi:10.1063/1.473692 ; Xu6560  or subjected to very high
pressuredoi:10.1063/1.3698601 ; PhysRevE.88.052406 ; 27802615 . Delhommelle
and CummingsPhysRevB.72.172201  used TTCF to study the frictional response of
a Weeks-Chandler-Anderson (WCA)doi:10.1063/1.1674820  fluid confined to a film
of about five molecular diameters over a wide range of shear rates
($10^{3}-10^{11}\;\text{s}^{-1}$). However, the wall atoms lacked thermal
motion and the fluid atoms were thermostatted, which makes an analytical
derivation of the dissipation function difficult. Bernardi et al. used TTCF to
study friction in WCA and Lennard-Jones (LJ) fluids confined between wall
atoms with thermal motion22920110 . This meant that there was no ambiguity in
the definition of the dissipation, which was analytically derived following
its mathematical definitiondoi:10.1063/1.2812241 . More recently, Bernardi et
al.doi:10.1080/08927022.2015.1049174  used TTCF to study the low shear rate
rheology of confined polymer chains of several lengths (1, 2, 4, 8 and 12
beads) represented by the finite extensible nonlinear elastic (FENE)
potentialdoi:10.1063/1.458541 .

In all previous TTCF studies of confined systemsPhysRevB.72.172201 ; 22920110
; doi:10.1080/08927022.2015.1049174  fixed channel widths have been employed
with a constant fluid volume. Experimentally, the channel width can often vary
dynamically in response to pressure, shear, and temperature. Certain important
behaviour, such as shear dilatancy26039993 , can only be captured in confined
NEMD simulations if the channel width can changeXu6560 . This can be achieved
by applying a barostat to the confining wallsdoi:10.1063/1.459524 ;
PhysRevE.90.043302 . Quantitatively different flow and friction behaviour has
been observed for systems at constant pressure rather than fixed channel
width, i.e. constant volumePhysRevE.90.043302 . In this study, we use TTCF to
investigate the friction of confined LJ and WCA fluids at low shear rates
under the conditions of both constant volume and pressure. The temperature and
pressure are controlled using a thermostat and barostat acting only on the
wall atoms to closely mimic experimental conditions. In what follows, we first
outline our methodology and establish the governing equations of motion and
their associated TTCF expressions. Next, we present results for our TTCF
computations for the shear pressure (negative of the shear stress) and slip
velocity and compare them with their standard direct-average NEMD values. We
end with some concluding remarks and present the TTCF derivation in the
Appendix.

##  II Methodology and computation details

We firstly note that all of our simulations were performed using our own in-
house code. Lupkowski and van Swoldoi:10.1063/1.459524  invented a boundary-
controlled barostat for controlling the normal pressure in molecular dynamics
simulations of liquids confined between planar fluctuating walls, where a wall
is represented by a solid block with one degree of freedom. Extensions of this
barostat were derived by Gattinoni et al.PhysRevE.90.043302  who incorporated
an atomistic description of the wall and where the atoms are tethered to an
underlying ‘virtual’ rigid lattice which plays a similar role to the Lupkowski
and van Swol barostat wall. In this work we adopt a mixed approach, where for
our 3-dimensional system the wall particles are bound to a set of virtual
lattice sites free to move along the $y$-direction (the confinement direction)
and on which an external constant force is applied. The shear is likewise
exerted on the system by applying a constant velocity to the lattice sites
along the $x$-direction. The geometry of the system is depicted in Fig. 1 and
the resulting equations of motion are:

$\begin{split}&\text{constant volume}\\\
\dot{\textbf{r}}^{f}_{i}&=\dfrac{\textbf{p}^{f}_{i}}{m_{i}}\\\
\dot{\textbf{p}}^{f}_{i}&=\textbf{F}^{2B}_{i}\\\
\dot{\textbf{r}}^{w}_{i}&=\dfrac{\textbf{p}^{w}_{i}}{m_{i}}\\\
\dot{\textbf{p}}^{w}_{i}&=\textbf{F}^{2B}_{i}-k(\textbf{r}^{w}_{i}-\textbf{r}^{l}_{i})-\alpha\textbf{p}^{w}_{i}\\\
\dot{\alpha}&=\dfrac{1}{Q}\biggl{(}\sum_{i}^{N^{w}}\textbf{p}^{w2}_{i}-3N^{w}k_{B}T\biggr{)}\\\
\dot{\textbf{r}}^{l}_{i}&=\bigl{(}\pm v\;\;,\;\;0\;\;,\;\;0\bigr{)}\\\
\dot{\textbf{p}}^{l}_{i}&=\bigl{(}0\;\;,\;\;0\;\;,\;\;0\bigr{)}.\end{split}\quad\quad\quad\begin{split}&\text{constant
pressure}\\\ \dot{\textbf{r}}^{f}_{i}&=\dfrac{\textbf{p}^{f}_{i}}{m_{i}}\\\
\dot{\textbf{p}}^{f}_{i}&=\textbf{F}^{2B}_{i}\\\
\dot{\textbf{r}}^{w}_{i}&=\dfrac{\textbf{p}^{w}_{i}}{m_{i}}\\\
\dot{\textbf{p}}^{w}_{i}&=\textbf{F}^{2B}_{i}-k(\textbf{r}^{w}_{i}-\textbf{r}^{l}_{i})-\alpha\textbf{p}^{w}_{i}\\\
\dot{\alpha}&=\dfrac{1}{Q}\biggl{(}\sum_{i}^{N^{w}}\textbf{p}^{w2}_{i}-3N^{w}k_{B}T\biggr{)}\\\
\dot{\textbf{r}}^{l}_{i}&=\bigl{(}\pm
v\;,\;\dfrac{p_{y}^{l}}{m_{i}}\;,\;0\bigr{)}\\\
\dot{\textbf{p}}^{l}_{i}&=\bigl{(}0\;,\;-F^{ext}+\dfrac{1}{N^{l}}\sum_{i}^{N^{l}}k(r_{yi}^{w}-r_{yi}^{l})\;,\;0\bigr{)}.\end{split}$
(1)

Here, the superscripts $f$, $w$ and $l$ denote, respectively, the fluid, wall,
and lattice particles. $\textbf{F}_{i}^{2B}=-\sum_{i\neq j}\nabla\phi_{ij}$ is
the interatomic two-body force, with

$\phi\left(r_{ij}\right)=\begin{cases}4\epsilon\biggl{[}\biggl{(}\dfrac{\sigma}{r_{ij}}\biggr{)}^{12}-\biggl{(}\dfrac{\sigma}{r_{ij}}\biggr{)}^{6}\biggr{]}+\phi_{c},&\text{if
}r_{ij}\leq r_{c}\\\ 0,&\text{if }r_{ij}>r_{c}\end{cases}$ (2)

and $\textbf{F}_{i}^{H}=-k(\textbf{r}_{i}^{w}-\textbf{r}_{i}^{l})$ is the
harmonic force tethering the wall particles and the lattice sites (see Fig.
2), while $\sigma$ and $\epsilon$ are the effective diameter and potential
well, respectively. The cut-off radius is $r_{c}=2^{1/6}\sigma$ for Weeks-
Chandler-Andersendoi:10.1063/1.1674820  (WCA) particles and $r_{c}=2.2\sigma$
for Lennard-Jonesdoi:10.1098/rspa.1924.0081  (LJ) ones. The term $\phi_{c}$
has the purpose to eliminate the discontinuity in the potential at $r=r_{c}$
and $\phi_{c}=\epsilon$ for the WCA potential and $\phi_{c}\simeq
0.34\epsilon$ for the LJ. As will be clear from the following results, systems
composed of simple atomistic WCA particles hardly exhibit any slip velocity at
the wall-fluid interface. In order to provide a useful example of the TTCF
formalism in the study of slip velocity, we investigated an analogous systems
with particles interacting with a Lennard-Jones potential and where the
parameters $\sigma$ and $\epsilon$ have been specifically modulated to promote
a slip at the fluid-wall interface. In what follows, all quantities are
expressed in reduced units, with $\sigma=\epsilon=1$ for the WCA systems, and
$\sigma^{f}=1$, $\sigma^{w}=0.5$, $\epsilon^{f}=\epsilon^{w}=1$ for the LJ
one. In the latter, the interaction between fluid and wall particles has been
modelled via the Lorentz-Berthelot mixing rule. For the case of constant
volume, the lattices sites are fixed, and an initial velocity $\pm v$
(positive for the upper wall, negative for the lower one) is imposed in the
$x$-direction. For the systems at constant pressure, the sites are free to
move along the $y$-direction and are subjected to the mean harmonic force and
to an additional external constant force $F^{ext}$ applied homogeneously to
each site, which must be set to match the desired pressure. Since each site is
subjected to the same force, the lattice sites move collectively as a single
rigid body.

Figure 1: Schematic of the system and location of the planes across the
channel. The $z$-direction is normal to the page. The diagram is merely
representative and the size of the particles may not be proportionate to the
actual system size.

For any generic phase variable $B\left(t\right)$, the TTCF formalism is based
on the following identitydoi:10.1063/1.2812241 ; 22920110 :

$\langle B(t)\rangle=\langle
B(0)\rangle+\int_{0}^{t}\langle\Omega(0)B(s)\rangle\text{d}s$ (3)

which relates the phase space average of $B$ at time $t$ with the time
integral of the correlation with the dissipation function
$\Omega(0)=\beta\dot{H}^{ad}$, i.e. the rate of energy dissipation without
accounting for the thermostat term at $t=0$, when the external force is
switched ondoi:10.1063/1.3684753 . For our system we find that
$\Omega=-\beta\sum_{i}k(r_{xi}^{w}-r_{xi}^{l})v$ for both the constant volume
and the constant pressure case. See the Appendix for a brief derivation of the
dissipation function in the presence of a barostat (the constant volume
derivation results in the same expression and its derivation is analogous).

Figure 2: Schematic of the external field acting on the lattice sites (black
dots). The lattice is subjected to both the tethering force (only some wall
particles are actually shown, for the sake of clarity) and the external force.
The action of the force occurs only along the $y$ (vertical) direction, and it
does not affect the dynamics in the $x$ or $z$-directions.

In summary, for the application of the TTCF method in computer simulations,
one must generate a number $N_{s}$ of uncorrelated equilibrium ($v=0$) phase
state points, which act as initial conditions of $N_{s}$ nonequilibrium
systems ($v\neq 0$) on which $B(t)$ is computed. The multiple time series are
then averaged over the $N_{s}$ measurements and integrated as per Eq. (3). To
improve the statistics, it is useful to generate several different
nonequilibrium trajectories from the same starting point, by generating
transformed, or mapped, initial conditions from each starting point. In this
work, we adopt four sets of phase-space mappings: for each particle $i$ the
mappings are generated by the following rules, taken from previous
works:PhysRevB.72.172201

$\begin{split}\bigl{(}x_{i}\;,\;y_{i}\;,\;z_{i}\;,\;p_{xi}\;,\;p_{yi}\;,\;p_{zi}\bigr{)}&\longrightarrow\bigl{(}x_{i}\;,\;y_{i}\;,\;z_{i}\;,\;-p_{xi}\;,\;-p_{yi}\;,\;-p_{zi}\bigr{)}\\\
\bigl{(}x_{i}\;,\;y_{i}\;,\;z_{i}\;,\;p_{xi}\;,\;p_{yi}\;,\;p_{zi}\bigr{)}&\longrightarrow\bigl{(}-x_{i}\;,\;y_{i}\;,\;z_{i}\;,\;-p_{xi}\;,\;p_{yi}\;,\;p_{zi}\bigr{)}\\\
\bigl{(}x_{i}\;,\;y_{i}\;,\;z_{i}\;,\;p_{xi}\;,\;p_{yi}\;,\;p_{zi}\bigr{)}&\longrightarrow\bigl{(}-x_{i}\;,\;y_{i}\;,\;z_{i}\;,\;p_{xi}\;,\;-p_{yi}\;,\;-p_{zi}\bigr{)}.\\\
\end{split}$ (4)

It is straightforward to verify that the above transformations are compatible
with the canonical ensemble, which is symmetric by any of those
transformations. This procedure ensures that $\langle\Omega(0)\rangle$ is
identically $0$, as expected theoretically, and thus the integral of the time
correlation will converge. Further details of the TTCF methodology and the use
of phase-space mappings can be found in the specialist books on nonequilibrium
statistical mechanics and nonequilibrium molecular dynamics EvansMorriss2008 ;
ToddDaivis2017 .

We computed the $P_{xy}$ component of the local pressure tensor (the negative
of the shear stress, or shear pressure) using the method of planes (MoP)
techniquePhysRevE.52.1627  at different locations across the channel, as per
Figure 1 for the WCA systems. For constant pressure simulations the planes
themselves are instantaneously adjusted since the channel width can fluctuate.
While the absolute plane positions can fluctuate in time, their relative
positions with respect to channel width remains constant. In addition, we
generated velocity and density profiles of both WCA and LJ cases. Since the
dissipation function $\Omega$ is strongly localized in the wall region, the
calculation of any of the TTCF quantities at various distances from the walls
allowed us also to verify whether the correlation
$\langle\Omega\left(0\right)B\left(t\right)\rangle$ decays when
$B\left(t\right)$ is computed far from the wall region, potentially making the
TTCF method less effective.

NEMD simulations for WCA systems composed of $3200$ fluid particles and $800$
wall particles were performed. Each wall was composed of two layers. All
systems had initial density $\rho=0.8442$, $L\simeq 16\sigma=16$ and were
thermostatted with a Nose-Hoover thermostat Nose1984 ; Hoover1985  at
temperature $T=1$ at the walls. For the case of constant pressure, we analysed
systems at relatively high pressures of $P=10$ and $P=11$. Higher pressures
than this might activate a phase transition where the fluid region becomes
solid doi:10.1063/1.3698601 ; PhysRevE.88.052406 ; 27802615 , a condition for
which the TTCF formalism fails or becomes highly inefficient. This phenomenon
might also be promoted by a small system size along the $x$-direction. For the
LJ case, we investigated systems at fluid density $\rho^{f}=0.8$ and
$\rho^{w}=8$, and a channel width of $\simeq 10\sigma$. The high density of
the wall is due to the reduced diameter of the wall particles, and corresponds
to $\rho^{w}=1$ if the quantities were normalized assuming $\sigma^{w}=1$. A
recent workPhysRevE.100.023101  pointed out that small inhomogeneous systems,
under certain conditions, tend to underestimate the solid-liquid friction.
However, in this work we have not accounted for this effect, as we focussed
solely on the derivation and application of the TTCF formalism under constant
pressure conditions for the study of highly confined fluid systems. The
reduced mass of each particle and lattice site was set to $1$, and the
stiffness of the spring was set as $k=150$ for the WCA systems and $k=75$ for
the LJ ones PhysRevE.55.4288 .

Each system underwent an equilibration $400$ time units, and subsequently
$10^{5}$ initial conditions were sampled, for a total $4\times 10^{5}$
starting points for the TTCF calculation due to the mappings, with a delay of
$2$ time units between each sample along the equilibrium simulations. Each
nonequilibrium trajectory was followed for $10$ time units. Since the TTCF
method is based on a time correlation at time $t=s$ and $t=0$, it is essential
that the equations of motion are integrated with a self-starting integrator;
this means that the popular Gear predictor-corrector would not be suitable for
the purpose. A fourth order Runge-Kutta is typically
useddoi:10.1080/08927020801930604 ; doi:10.1080/08927020802575598  which, on
the other hand, make the simulation computationally expensive. We found that
the faster velocity Verlet performs equally well. The equations of motion of
both the equilibrium and nonequilibrium systems were therefore integrated
using a velocity Verlet algorithm Verlet1967 , with time step $\delta
t=0.004$.

##  III Results

Figure 3 shows the comparison between the direct average (DAV) and the TTCF
methods in computing the shear pressure $P_{xy}\left(y\right)$ as a function
of time for WCA systems at the fluid/wall interface, for the case of a
constant volume system and reduced shear rates ranging between
$\dot{\gamma}=10^{-2}-10^{-5}$, which correspond approximately to $\simeq
5\times 10^{9}-5\times 10^{6}\;\text{s}^{-1}$ in MKS units, assuming that the
fluid is composed of monatomic argon atoms doi:10.1063/1.479848 . The DAV data
are obtained by simply averaging the time dependent relevant phase variable
(in this case, the negative of the shear stress $\sigma_{xy}$, where
$\sigma_{xy}\equiv-P_{xy}$) over multiple independent NEMD trajectories. The
number of these trajectories is simply the number of initial conditions (i.e.
$1\times 10^{5}$), whereas the number of TTCF trajectories, as described in
the previous section, is $4\times 10^{5}$ trajectories, due to the phase-space
mappings required. Even though the number of TTCF trajectories is four times
more than the number of DAV trajectories, we do not average over all four
phase-space mappings for the DAV data. The reason for this is that at short
times the DAV trajectories are highly correlated and fluctuations in
statistical error cancel due to mapping symmetry. However, Lyapunov
instability leads to a rapid de-correlation as time increases. In the long-
time limit (and certainly in the steady-state) the standard error would only
be improved by a factor of two if we were to use all $4\times 10^{5}$
trajectories for DAV. A factor of two is still insignificant compared to the
improvement in statistics we find with the TTCF results, as we will
demonstrate shortly.

We observe at early times the negative shear stress overshootsHeyes1980 ,
followed by damping oscillatory behaviour, similar to that observed by
Bernardi et. al 22920110 , due to the atomic vibrations in wall atoms (Fig.
3). For the highest shear rate, the statistical accuracy of the two methods is
comparable. The data show that, for the systems analysed, the standard error
of the signal $P_{xy}$ computed with the DAV method is in the order of
$10^{-2}$ or higher, and roughly constant for each level of the external
force. This means that for $\dot{\gamma}<10^{-2}$ the time evolution of
$P_{xy}$ is either inaccurate or totally indistinguishable from random noise,
as can be seen in the last two plots. On the other hand, the TTCF method
produces an extremely accurate and clean signal for arbitrarily low shear
rates, with a standard error $\epsilon_{std}$ decreasing with the shear rate
and spanning from approximately $\epsilon_{std}\approx 10^{-2}$ for
$\dot{\gamma}=10^{-2}$ to $\epsilon_{std}\approx 10^{-5}$ for
$\dot{\gamma}=10^{-5}$. The constant reduction of the uncertainty in $P_{xy}$
computed with the TTCF method indicates that an exceptionally good signal-to-
noise ratio can be achieved at shear rates lower that those adopted in this
work.

Figure 3: $P_{xy}$ at the fluid/wall interface, obtained by averaging over
planes 2 and 6 (see Fig. 1), for the WCA systems at constant volume and for
different shear rates. The size of the error bars is 4 times the standard
error. For the sake of clarity, the error bars of the DAV method for the two
lowest shear rates are not displayed and are roughly an order of magnitude
larger than the signal itself. A constant reduction of the uncertainty in the
TTCF time series indicates that this method would be highly accurate for even
lower shear rates.

In Figure 4 we see $P_{xy}$ at the various locations across the channel for
$\dot{\gamma}=10^{-5}$, computed with the TTCF method under conditions of
constant volume. The shear pressure in the fluid region (planes 3, 5 and 4)
converges to a constant value across the channel, as expected theoretically.
However, the transient region ($t<6$) shows little or no oscillatory behaviour
found close to the wall. In the wall region (planes 1, 7) the average shear
pressure is essentially null, due to the constant drag of the lattice sites;
however, the oscillatory behaviour due to the vibrating wall atoms is still
apparent. The signal is particularly stable inside the walls and at the wall-
fluid interface (planes 2, 6), due to the absence of the kinetic contribution,
which is inherently noisier. Each curve retains similar levels of uncertainty.
This indicates that a strong correlation between $\Omega\left(0\right)$ and
the fluid evolution holds arbitrarily far from the walls. We also note that
the steady-state shear pressure in the fluid region close to the wall (planes
3, 5) is reached more quickly than that near the channel centre (plane 4) due
to the time it takes for momentum transfer from the walls into the bulk of the
fluid.

Figure 4: $P_{xy}$ at various locations in the system (see Fig. 1) under
constant volume for $\dot{\gamma}=10^{-5}$, computed with the TTCF formalism.
The size of the error bars is 4 times the standard error.

Figures 5(a) and 5(b) show the comparison between the DAV and TTCF methods in
the computation of the streaming velocity profile at low shear rates
($\dot{\gamma}=10^{-5}$) for the systems at constant volume and constant
pressure, defined as

$\displaystyle\mathbf{v}\left(\mathbf{r}_{\text{bin}},t\right)=\frac{\langle\sum_{i\in{\text{bin}}}m_{i}\mathbf{v}_{i}\rangle}{\langle\sum_{i\in{\text{bin}}}m_{i}\rangle},$
(5)

where the sum ranges over all fluid atoms within a bin of finite width with
mid-point located at $\mathbf{r}_{\text{bin}}$. The standard error of the
velocity computed via DAV is in the order of $10^{-2}$ which is more than two
orders of magnitude larger than the signal itself. The velocity profile is
critical in the calculation of the slip velocity and the corresponding slip
length. In confined boundary driven systems on smooth surfaces, the velocity
profile of the fluid deviates from the linear, no-slip velocity profile
imposed by the motion of the walls. In particular, the average velocity of the
fluid particles in proximity of the wall is smaller than that of the wall and
lower than the theoretical one PhysRevA.41.6830 . Of particular significance
in these results is that TTCF has for the first time been used to determine
accurate velocity profiles for a confined fluid under very low (from an NEMD
simulation perspective) rates of strain that approach those values achievable
under laboratory conditions Neto_2005 . The direct average (DAV) of the signal
is simply too noisy and no information can be gleaned. However, the TTCF
velocity profiles, as seen from the insets in Figs. 5(a) and 5(b), demonstrate
that the signal can be clearly distinguished from the noise. This in itself is
a significant achievement and has never been seen before. It finally opens the
way for the use of NEMD to precisely compute stresses, velocities or any other
phase variable whatsoever, at physically meaningful and accessible laboratory
strain rates.

(a)

(b)

Figure 5: Velocity profile of the WCA system at constant pressure (left) at
$t=10$ and LJ system (right) at $t=7.6$, for $\dot{\gamma}=10^{-5}$. The wall
region is included, and can be recognised from the structured peaks in the
velocity at the borders of the region. The vertical black dashed lines
indicate the positions of the wall-fluid interface. The standard error of the
DAV method has been omitted and it has a magnitude in the range of $10^{-2}$.
Note the clearly defined linear velocity profile observed in the TTCF signal,
as opposed to the direct average signal, which is too noisy to allow
observation of any meaningful profile. The time $t=7.6$ for the LJ system has
been chosen to show the slip in the velocity profile (compare with Figure 7).

Figures 6 and 7 show more clearly the slip velocity, defined here as the
difference in the velocity between the inner layer of the wall particles and
the first layer of the fluid particles for WCA and LJ systems, respectively.
In our systems, the layers have a width of roughly $\sigma$. The value has
been taken as the average of the two interfaces, with the convention that the
velocities are positive, i.e. a slip velocity greater than zero indicates that
that streaming velocity of the wall is larger than that of the fluid. The data
show again the power of the TTCF method. The computation of the slip velocity
at low shear rates would be unfeasible if the direct average method was used,
and it is typically difficult to estimate directly even at high shear rates,
particularly for high-slip systems doi:10.1063/1.3675904 , due to the large
fluctuations occurring in the velocity profile. It is precisely for such
reasons that equilibrium molecular dynamics methods have been developed to
compute the friction coefficients for such systems BB1994 ; HTD2011 , and have
been shown to be more reliable than NEMD for systems of high slip, such as
water confined to graphene or carbon nanotubes doi:10.1063/1.3675904 ;
Kannam2013 . The TTCF results generated here now arm us with a powerful NEMD
technique to obtain slip velocities, slip lengths and friction coefficients
under actual nonequilibrium conditions. We further note that the oscillatory
behaviour of the slip velocity is most likely due to the vibrations of wall
atoms tethered to their lattice sites, rather than evidence of stick-slip
behaviour such as that seen previously doi:10.1126/science.250.4982.792 . A
more detailed analysis of this oscillatory behaviour will follow in a future
publication.

In Figure 8 the slip velocity, computed with the TTCF method for all the WCA
systems, is displayed. The drop in the slip velocity is approximately
proportional to the shear rate, and it is apparently not affected by the
pressure imposed on the system except for the lowest shear rate. Tables 1, 2
and 3 summarise the data. As can be seen, the error of the direct average
rapidly dominates the signal for both the quantities measured. It is useful to
compare the signal-to-noise ratio (SNR) for the DAV and the TTCF signals,
defined as the ratio between the absolute value of the average of the signal
and its standard error, at $t=10$. As can be seen in Figure 9 the SNR remains
constant for the TTCF method for each level of the external field. This is a
clear and promising evidence that accurate measurements can be achieved with
the same system size and number of samples for shear rates lower than studied
in this work, approaching values which are practical.

(a)

(b)

Figure 6: Slip velocity at the interface for the WCA system at constant volume
(a) and pressure (b) for $\dot{\gamma}=10^{-5}$. The error bars for the DAV
measurements have been omitted since they are at least an order of magnitude
larger than the signal.

(a)

(b)

Figure 7: Slip velocity at the interface at constant volume for $\dot{\gamma}=10^{-2}$ (a) and $\dot{\gamma}=10^{-5}$ (b) for the LJ system. Figure 8: Slip velocity for all the systems investigated, computed with the TTCF method. $P\simeq 7.7$ represents the systems at constant volume. There is no appreciable change in the behaviour of the systems at the interface in the various systems, which means that none of them show slip velocity. Error bars are 4 times the standard error. Table 1:  Summary of the estimated slip velocity for the systems analysed. The system at $P\simeq 7.7$ is the case at constant volume. Aside from the values itself, the superiority of the TTCF method can be appreciated by the reduced uncertainty in the final estimate. The $n$ values in brackets are twice the standard error and is the error in the last $n$ digits of the slip velocities quoted. $P$ | $\simeq 7.7$ | $10$ | $11$
---|---|---|---
$\dot{\gamma}$ | DAV | TTCF | DAV | TTCF | DAV | TTCF
$0.01$ | 0.0018(44) | -0.0012(55) | 0.0005(49) | 0.0006(96) | -0.0004(47) | 0.0012(80)
$0.001$ | 0.0004(46) | 0.00008(5) | 0.0000(48) | -0.0003(4) | 0.0000(47) | -0.0003(3)
$0.0001$ | -0.0005(45) | -0.00001(6) | 0.0001(48) | 0.00001(3) | 0.0003(46) | -0.00000(2)
$0.00001$ | 0.0003(45) | 0.00000(0) | 0.0000(50) | 0.00000(0) | 0.0000(50) | 0.00000(0)
Table 2:  Summary of the estimated slip velocity for the LJ system. Aside from the values itself, the superiority of the TTCF method can be appreciated by the reduced uncertainty in the final estimate. The $n$ values in brackets are twice the standard error and is the error in the last $n$ digits of the slip velocities quoted. $\dot{\gamma}$ | DAV | TTCF
---|---|---
$0.01$ | 0.01688(3733) | 0.01584(1088)
$0.001$ | 0.00175(3817) | 0.00172(120)
$0.0001$ | 0.00015(3723) | 0.00019(12)
$0.00001$ | 0.00004(3880) | 0.00001(1)
Table 3:  Summary of the estimated shear pressure (negative shear stress) in the fluid region, as the average over the planes $2-6$. The system at $P\simeq 7.7$ is the case at constant volume. The $n$ values in brackets are twice the standard error and is the error in the last $n$ digits of the shear pressures quoted. $P$ | $\simeq 7.7$ | $10$ | $11$
---|---|---|---
$\dot{\gamma}$ | DAV | TTCF | DAV | TTCF | DAV | TTCF
$0.01$ | -0.0274(251) | -0.0205(117) | -0.0443(284) | -0.0476(134) | -0.0563(290) | -0.0475(151)
$0.001$ | -0.0027(259) | -0.0039(11) | -0.0045(283) | -0.0052(4) | -0.0057(268) | -0.0062(016)
$0.0001$ | -0.0003(247) | -0.0003(1) | -0.0004(283) | -0.0005(1) | -0.0006(286) | -0.0006(2)
$0.00001$ | 0.0000(261) | -0.00003(0) | 0.00003(274) | -0.00005(0) | 0.0000(274) | -0.00006(0)
Figure 9: Signal-to-noise ratio for the shear pressure $P_{xy}$ of the WCA
systems and for the slip velocity of the LJ system at the various shear rates.
The results of the slip velocity for the WCA systems have been omitted since
the signal itself is close to 0. The quality of the signal drops rapidly for
low shear rates for the direct average, but remains unchanged for the TTCF.

##  IV Conclusions

We have derived the analytical TTCF expression for boundary sheared systems at
constant pressure, and have shown that it is equal to that at constant volume
as long as the barostat is modelled via a conservative force. We have studied
systems at different pressure and for various driving boundary velocities, at
realistic strain rates. For all the quantities monitored, we observed an
increasing level of noise in the DAV signal, making the signal-to-noise ratio
rapidly decay to zero. On the other hand, the signal-to-noise ratio remained
constant with the TTCF method and was substantially higher, indicating that
accurate measurements can be obtained for much weaker external driving forces.
Our results show that the correlation with the dissipation function and both
the shear pressure (negative of the shear stress) and the streaming velocity
does not substantially decay when moving far from the wall. Although the
limited size of the systems along the $x$-direction could affect the final
value of the shear pressure, we have shown that the TTCF method can be applied
in the study of friction in nanochannels at realistic rates of strain. The
generality of the TTCF approach, and its simplicity, make it readily
expandable to more complex systems such as molecular fluids and mixtures and
opens up a new frontier in the use of NEMD simulation to study the underlying
physical processes that determine the tribology of highly confined systems at
realisable strain rates.

##  V Acknowledgements

The authors would like to thank Prof. D. M. Heyes (Department of Mechanical
Engineering, Imperial College London, UK) for useful discussions. We also
thank the Australian Research Council for a grant obtained through the
Discovery Projects Scheme (DP200100422) and the Royal Society for support via
International Exchanges grant $IES\ R3\ 170233$. J.P.E. and D.D. acknowledge
the financial support of the Engineering and Physical Sciences Research
Council (EPSRC) via EP/ N025954/1 and EP/P030211/1. J.P.E. was supported by
the Royal Academy of Engineering through the Research Fellowships scheme. We
acknowledge the Swinburne OzSTAR Supercomputing facility and Imperial College
London Research Computing Service for providing computational resources for
this work.

##  VI Data availability statement

The data that support the findings of this study are available from the
corresponding author or by emailing<EMAIL_ADDRESS>upon reasonable
request.

##  VII Appendix

In this Appendix, we provide a derivation of the dissipation function for
barostatted systems, and we show that it is equivalent to that of a system
with fixed walls. We recall the equations of motion

$\begin{split}\dot{\textbf{r}}^{f}_{i}&=\dfrac{\textbf{p}_{i}^{f}}{m_{i}}\\\
\dot{\textbf{p}}^{f}_{i}&=\textbf{F}_{i}^{2B}\end{split}$ (6)

$\begin{split}\dot{\textbf{r}}^{w}_{i}&=\dfrac{\textbf{p}_{i}^{w}}{m_{i}}\\\
\dot{\textbf{p}}^{w}_{i}&=\textbf{F}_{i}^{2B}+\textbf{F}_{i}^{H}-\alpha\textbf{p}^{w}_{i}\end{split}$
(7)
$\begin{split}\dot{\textbf{r}}^{l}_{i}&=\biggl{(}v\;\;,\;\;\dfrac{p_{y}^{l}}{m_{i}}\;\;,\;\;0\biggr{)}\\\
\dot{\textbf{p}}^{l}_{i}&=\biggl{(}0\;\;,\;\;-F^{ext}-\dfrac{\sum_{i}^{N^{l}}F_{yi}^{H}}{N^{l}}\;\;,\;\;0\biggr{)},\end{split}$
(8)

where, $\textbf{F}_{i}^{H}=-k(\textbf{r}_{i}^{w}-\textbf{r}_{i}^{l})$ and the
sign of the external force is negative because we are moving the upper wall
(i.e. the sign depends on the direction of the wall that is being driven). The
difference between upper and lower wall velocity is ignored, since it does not
affect the calculations.
The total energy $U=K+V$ can now be expressed as

$\begin{split}K&=\sum_{i}^{N^{f}}\dfrac{\textbf{p}_{i}^{2f}}{2m_{i}}+\sum_{i}^{N^{w}}\dfrac{\textbf{p}_{i}^{2w}}{2m_{i}}+N^{l}\dfrac{p_{y}^{2l}}{2m_{i}}\\\
V&=\sum_{i}^{N^{t}}\sum^{N^{t}}_{j>i}\phi_{ij}^{2B}+\sum_{i}^{N^{w}}\phi_{i}^{H}+F^{ext}\sum_{i}^{N^{l}}(r_{yi}^{l}-r_{yi0}^{l})\end{split}$
(9)

with $N^{t}=N^{f}+N^{w}$,
$\phi_{i}^{H}=\frac{1}{2}k(\textbf{r}_{i}^{w}-\textbf{r}_{i}^{l})^{2}$ and
$r_{yi0}^{l}$ the arbitrary reference position of the lattice site $i$. The
total time derivative of the energy $U$, which is proportional to the
dissipation function we require ($\Omega$), is

$\begin{split}\dfrac{\text{d}U}{\text{d}t}=\dfrac{\partial U}{\partial
t}+&\sum_{i}^{Nf}\biggl{(}\dfrac{\partial}{\partial\textbf{r}^{f}_{i}}\cdot\dot{\textbf{r}}^{f}_{i}+\dfrac{\partial}{\partial\textbf{p}^{f}_{i}}\cdot\dot{\textbf{p}}^{f}_{i}\biggr{)}U+\\\
&\sum_{i}^{Nw}\biggl{(}\dfrac{\partial}{\partial\textbf{r}^{w}_{i}}\cdot\dot{\textbf{r}}^{w}_{i}+\dfrac{\partial}{\partial\textbf{p}^{w}_{i}}\cdot\dot{\textbf{p}}^{w}_{i}\biggr{)}U+\sum_{i}^{Nl}\biggl{(}\dfrac{\partial}{\partial\textbf{r}^{l}_{i}}\cdot\dot{\textbf{r}}^{l}_{i}+\dfrac{\partial}{\partial\textbf{p}^{l}_{i}}\cdot\dot{\textbf{p}}^{l}_{i}\biggr{)}U.\end{split}$
(10)

Recalling that for any generic particle $k$ we have

$\begin{split}\dfrac{\partial}{\partial\textbf{r}_{k}}\biggl{(}\sum_{i}^{Nt}\sum^{Nt}_{j>i}\phi_{ij}^{2B}\biggr{)}&=-\textbf{F}_{k}^{2B}\\\
\dfrac{\partial}{\partial\textbf{r}^{w}_{k}}\biggl{(}\sum_{i}^{Nw}\phi_{i}^{H}\biggr{)}&=-\textbf{F}_{k}^{H}\\\
\dfrac{\partial}{\partial\textbf{r}^{l}_{l}}\biggl{(}\sum_{i}^{Nw}\phi_{i}^{H}\biggr{)}&=\textbf{F}_{k}^{H},\end{split}$
(11)

then the total time derivatives of the kinetic and potential energy are

$\begin{split}\dfrac{\text{d}K}{\text{d}t}&=\sum_{i}\biggl{(}\dfrac{\partial}{\partial\textbf{p}_{i}}\cdot\dot{\textbf{p}}_{i}\biggr{)}K=\sum_{i}^{Nf}\textbf{p}^{f}_{i}\cdot(\textbf{F}_{i}^{2B})+\sum_{i}^{Nw}\textbf{p}^{w}_{i}\cdot(\textbf{F}_{i}^{2B}+\textbf{F}_{i}^{H})+\sum_{i}^{N^{l}}\dfrac{p^{l}_{y}}{m_{i}}(-F^{ext}-\dfrac{\sum_{i}F_{yi}^{H}}{N^{l}})\end{split}$
(12)

and

$\begin{split}\dfrac{\text{d}V}{\text{d}t}&=\sum_{i}\biggl{(}\dfrac{\partial}{\partial\textbf{r}_{i}}\cdot\dot{\textbf{r}}_{i}\biggr{)}V=\sum_{i}^{Nf}\dfrac{\textbf{p}^{f}_{i}}{m_{i}}\cdot(-\textbf{F}_{i}^{2B})+\sum_{i}^{Nw}\dfrac{\textbf{p}^{w}_{i}}{m_{i}}\cdot(-\textbf{F}_{i}^{2B}-\textbf{F}_{i}^{H})+\sum_{i}^{Nl}\dot{\textbf{r}}_{i}^{l}\cdot(\textbf{F}_{i}^{H})+F^{ext}\sum_{i}^{N^{l}}\dfrac{p^{l}_{y}}{m_{i}}.\end{split}$
(13)

With the assumption that all the masses $m_{i}$ are identical and equal to
$m$, Eq. (13) can be simplified and combined with the kinetic contribution
into

$\dfrac{\text{d}U}{\text{d}t}=\dfrac{N^{l}p_{y}^{l}}{m}\biggl{(}-\sum_{i}^{N^{l}}\dfrac{F_{yi}^{H}}{N^{l}}\biggr{)}+\sum_{i}^{N^{l}}\dot{\textbf{r}}_{i}^{l}\cdot\textbf{F}_{i}^{H}.$
(14)

Recalling that for a particle $k$ of the lattice

$\dot{\textbf{r}}_{k}^{l}\cdot\textbf{F}_{k}^{H}=\biggl{(}v\;\;,\;\;\dfrac{p_{y}^{l}}{M}\;\;,\;\;0\biggr{)}\cdot\biggl{(}F^{H}_{xi}\;\;,\;\;F^{H}_{yi}\;\;,\;\;F^{H}_{zi}\biggr{)}^{T}=F^{H}_{xi}v+F^{H}_{yi}\dfrac{p_{y}^{l}}{m}$
(15)

we finally have

$\dfrac{\text{d}U}{\text{d}t}=-\dfrac{p_{y}^{l}}{m}\sum_{i}^{N^{l}}F^{H}_{yi}+\dfrac{p_{y}^{l}}{m}\sum_{i}^{N^{l}}F^{H}_{yi}+v\sum_{i}^{N^{l}}F^{H}_{xi}=-\sum_{i}^{Nl}k(r_{xi}^{w}-r_{xi}^{l})v,$
(16)

hence

$\Omega=\beta\frac{dU}{dt}=\beta\sum_{i}^{Nl}k(r_{xi}^{w}-r_{xi}^{l})v.$ (17)

We note that the dissipation function derived is identical to that of a system
with fixed walls. This property follows immediately from the conservative
nature of the forces involved in barostatting the system, and hence its
generality is quite wide and immediately applicable to several other methods
of pressure control for inhomogeneous systems.

## References

  * (1) J. P. Ewen, D. M. Heyes, and D. Dini, “Advances in nonequilibrium molecular dynamics simulations of lubricants and additives,” Friction, vol. 6, no. 4, pp. 349–386.
  * (2) A. Obliger, R. Pellenq, F.-J. Ulm, and B. Coasne, “Free volume theory of hydrocarbon mixture transport in nanoporous materials,” The Journal of Physical Chemistry Letters, vol. 7, no. 19, pp. 3712–3717, 2016.  PMID: 27570884.
  * (3) T. E. Gartner III and A. Jayaraman, “Modeling and simulations of polymers: A roadmap,” Macromolecules, vol. 52, p. 755, 2019.
  * (4) R. Taylor and B. de Kraker, “Shear rates in engines and implications for lubricant design,” Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology, vol. 231, no. 9, pp. 1106–1116, 2017.
  * (5) J. Vlachopoulos and D. Strutt, “Polymer processing,” Materials Science and Technology, vol. 19, pp. 1161–1169, 2019.
  * (6) A. C. Barbati, J. Desroches, A. Robisson, and G. H. McKinley, “Complex fluids and hydraulic fracturing,” Annual Review of Chemical and Biomolecular Engineering, vol. 7, no. 1, pp. 415–453, 2016.  PMID: 27070765.
  * (7) S. Plimpton, “Fast parallel algorithms for short-range molecular dynamics,” Journal of Computational Physics, vol. 117, no. 1, pp. 1–19, 1995.
  * (8) A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, “Lammps \- a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales,” Computer Physics Communications, vol. 271, p. 108171, 2022.
  * (9) V. Jadhao and M. O. Robbins, “Probing large viscosities in glass-formers with nonequilibrium simulations,” Proceedings of the National Academy of Sciences, vol. 114, no. 30, pp. 7952–7957, 2017.
  * (10) D. M. Heyes, “Transport coefficients of Lennard-Jones fluids: A molecular-dynamics and effective-hard-sphere treatment.,” Phys. Rev. B Condens Matter, vol. 37, no. 10, pp. 5677–5696, 1988.
  * (11) S. Bair, C. McCabe, and P. T. Cummings, “Comparison of nonequilibrium molecular dynamics with experimental measurements in the nonlinear shear-thinning regime,” Phys. Rev. Lett., vol. 88, p. 058302, Jan 2002\.
  * (12) J. P. Ewen, C. Gattinoni, J. Zhang, D. M. Heyes, H. A. Spikes, and D. Dini, “On the effect of confined fluid molecular structure on nonequilibrium phase behaviour and friction,” Phys. Chem. Chem. Phys., vol. 19, pp. 17883–17894, 2017.
  * (13) W. M. Visscher, “Transport processes in solids and linear-response theory,” Phys. Rev. A, vol. 10, pp. 2461–2472, Dec 1974.
  * (14) J. W. Dufty and M. J. Lindenfeld, “Nonlinear transport in the Boltzmann limit,” Journal of Statistical Physics, vol. 20, pp. 259–301, 1979.
  * (15) E. Cohen, “Kinetic theory of non-equilibrium fluids,” Physica A: Statistical Mechanics and its Applications, vol. 118, no. 1, pp. 17–42, 1983\.
  * (16) G. P. Morriss and D. J. Evans, “Application of transient correlation functions to shear flow far from equilibrium,” Phys. Rev. A, vol. 35, pp. 792–797, Jan 1987.
  * (17) D. J. Evans and G. P. Morriss, “Transient-time-correlation functions and the rheology of fluids,” Phys. Rev. A, vol. 38, pp. 4142–4148, Oct 1988.
  * (18) B. D. Todd and P. J. Daivis, “Homogeneous non-equilibrium molecular dynamics simulations of viscous flow: techniques and applications,” Molecular Simulation, vol. 33, no. 3, pp. 189–229, 2007.
  * (19) I. Borsak, P. T. Cummings, and D. J. Evans, “Shear viscosity of a simple fluid over a wide range of strain rates,” Molecular Physics, vol. 100, no. 16, pp. 2735–2738, 2002.
  * (20) O. A. Mazyar, G. Pan, and C. McCabe, “Transient time correlation function calculation of the viscosity of a molecular fluid at low shear rates: a comparison of stress tensors,” Molecular Physics, vol. 107, no. 14, pp. 1423–1429, 2009.
  * (21) C. Desgranges and J. Delhommelle, “Rheology of liquid fcc metals: Equilibrium and transient-time correlation-function nonequilibrium molecular dynamics simulations,” Phys. Rev. B, vol. 78, p. 184202, Nov 2008.
  * (22) C. Desgranges and J. Delhommelle, “Shear viscosity of liquid copper at experimentally accessible shear rates: Application of the transient-time correlation function formalism,” The Journal of Chemical Physics, vol. 128, no. 8, p. 084506, 2008.
  * (23) J. Delhommelle, P. T. Cummings, and J. Petravic, “Conductivity of molten sodium chloride in an arbitrarily weak dc electric field,” The Journal of Chemical Physics, vol. 123, no. 11, p. 114505, 2005.
  * (24) N. J. English, D. A. Mooney, and S. W. O’Brien, “Electrical conductivity and dipolar relaxation of binary dimethylimidazolium chloride–water solutions: A molecular dynamics study,” Journal of Molecular Liquids, vol. 157, no. 2-3, pp. 163–167, 2010.
  * (25) N. J. English, D. A. Mooney, and S. O’Brien, “Ionic liquids in external electric and electromagnetic fields: a molecular dynamics study,” Molecular Physics, vol. 109, no. 4, pp. 625–638, 2011.
  * (26) C. Desgranges and J. Delhommelle, “Estimating the conductivity of a nanoconfined liquid subjected to an experimentally accessible external field,” Molecular Simulation, vol. 34, no. 2, pp. 177–181, 2008.
  * (27) C. Desgranges and J. Delhommelle, “Molecular simulation of transport in nanopores: Application of the transient-time correlation function formalism,” Phys. Rev. E, vol. 77, p. 027701, Feb 2008.
  * (28) C. Desgranges and J. Delhommelle, “Accurate determination of normal stress differences via transient-time correlation function – non-equilibrium molecular dynamics (ttcf–nemd) simulations,” Molecular Simulation, vol. 35, no. 5, pp. 405–408, 2009.
  * (29) B. D. Todd, “Application of transient-time correlation functions to nonequilibrium molecular-dynamics simulations of elongational flow,” Phys. Rev. E, vol. 56, pp. 6723–6728, Dec 1997.
  * (30) B. D. Todd, “Nonlinear response theory for time-periodic elongational flows,” Phys. Rev. E, vol. 58, pp. 4587–4593, Oct 1998.
  * (31) R. Hartkamp, S. Bernardi, and B. D. Todd, “Transient-time correlation function applied to mixed shear and elongational flows,” The Journal of Chemical Physics, vol. 136, no. 6, p. 064105, 2012.
  * (32) J. Delhommelle and P. T. Cummings, “Simulation of friction in nanoconfined fluids for an arbitrarily low shear rate,” Phys. Rev. B, vol. 72, p. 172201, Nov 2005.
  * (33) S. Bernardi, S. J. Brookes, and D. J. Searles, “Response theory for confined systems,” J. Chem. Phys., vol. 137, p. 074114, Aug 2012.
  * (34) S. Bernardi and D. J. Searles, “Local response in nanopores,” Molecular Simulation, vol. 42, no. 6-7, pp. 463–473, 2016.
  * (35) W. T. Ashurst and W. G. Hoover, “Dense-fluid shear viscosity via nonequilibrium molecular dynamics,” Phys. Rev. A, vol. 11, pp. 658–678, Feb 1975.
  * (36) I. Bitsanis, J. J. Magda, M. Tirrell, and H. T. Davis, “Molecular dynamics of flow in micropores,” The Journal of Chemical Physics, vol. 87, no. 3, pp. 1733–1750, 1987.
  * (37) S. Y. Liem, D. Brown, and J. H. R. Clarke, “Investigation of the homogeneous-shear nonequilibrium-molecular-dynamics method,” Phys. Rev. A, vol. 45, pp. 3706–3713, Mar 1992.
  * (38) S. Bernardi, B. D. Todd, and D. J. Searles, “Thermostating highly confined fluids,” The Journal of Chemical Physics, vol. 132, no. 24, p. 244706, 2010\.
  * (39) X. Yong and L. T. Zhang, “Slip in nanoscale shear flow: mechanisms of interfacial friction,” Microfluidics and Nanofluidics, vol. 14, pp. 299–308, Jan 2013.
  * (40) S. Granick, “Motions and relaxations of confined liquids,” Science, vol. 253, no. 5026, pp. 1374–1379, 1991.
  * (41) P. A. Thompson and M. O. Robbins, “Origin of stick-slip motion in boundary lubrication,” Science, vol. 250, no. 4982, pp. 792–794, 1990.
  * (42) M. J. Stevens, M. Mondello, G. S. Grest, S. T. Cui, H. D. Cochran, and P. T. Cummings, “Comparison of shear flow of hexadecane in a confined geometry and in bulk,” The Journal of Chemical Physics, vol. 106, no. 17, pp. 7303–7314, 1997.
  * (43) R.-G. Xu and Y. Leng, “Squeezing and stick–slip friction behaviors of lubricants in boundary lubrication,” Proceedings of the National Academy of Sciences, vol. 115, no. 26, pp. 6560–6565, 2018.
  * (44) D. M. Heyes, E. R. Smith, D. Dini, H. A. Spikes, and T. A. Zaki, “Pressure dependence of confined liquid behavior subjected to boundary-driven shear,” The Journal of Chemical Physics, vol. 136, no. 13, p. 134705, 2012.
  * (45) C. Gattinoni, D. M. Heyes, C. D. Lorenz, and D. Dini, “Traction and nonequilibrium phase behavior of confined sheared liquids at high pressure,” Phys. Rev. E, vol. 88, p. 052406, Nov 2013.
  * (46) S. Maćkowiak, D. M. Heyes, D. Dini, and A. C. Brańka, “Non-equilibrium phase behavior and friction of confined molecular films under shear: A non-equilibrium molecular dynamics study.,” J. Chem. Phys., vol. 145, p. 164704, Oct 2016.
  * (47) J. D. Weeks, D. Chandler, and H. C. Andersen, “Role of repulsive forces in determining the equilibrium structure of simple liquids,” The Journal of Chemical Physics, vol. 54, no. 12, pp. 5237–5247, 1971.
  * (48) D. J. Evans, D. J. Searles, and S. R. Williams, “On the fluctuation theorem for the dissipation function and its connection with response theory,” The Journal of Chemical Physics, vol. 128, no. 1, p. 014504, 2008.
  * (49) K. Kremer and G. S. Grest, “Dynamics of entangled linear polymer melts: A molecular dynamics simulation,” The Journal of Chemical Physics, vol. 92, no. 8, pp. 5057–5086, 1990.
  * (50) I. Rosenhek-Goldian, N. Kampf, A. Yeredor, and J. Klein, “On the question of whether lubricants fluidize in stick-slip friction,” Proceedings of the National Academy of Sciences of the United States of America, vol. 112, no. 23, pp. 7117–7122, 2015.
  * (51) M. Lupkowski and F. van Swol, “Computer simulation of fluids interacting with fluctuating walls,” The Journal of Chemical Physics, vol. 93, no. 1, pp. 737–745, 1990.
  * (52) C. Gattinoni, S. Maćkowiak, D. M. Heyes, A. C. Brańka, and D. Dini, “Boundary-controlled barostats for slab geometries in molecular dynamics simulations,” Phys. Rev. E, vol. 90, p. 043302, Oct 2014.
  * (53) J. E. Jones, “On the determination of molecular fields. - II. From the variation of the viscosity of a gas with temperature,” Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, vol. 106, no. 738, pp. 441–462, 1924.
  * (54) D. J. Evans and G. Morriss, Statistical Mechanics of Nonequilibrium Liquids.  Cambridge University Press, 2008.
  * (55) B. D. Todd and P. J. Daivis, Nonequilibrium Molecular Dynamics: Theory, Algorithms and Applications.  Cambridge University Press, 2017.
  * (56) B. D. Todd, D. J. Evans, and P. J. Daivis, “Pressure tensor for inhomogeneous fluids,” Phys. Rev. E, vol. 52, pp. 1627–1638, Aug 1995.
  * (57) S. Nose, “A unified formulation of the constant temperature molecular dynamics methods,” J. Chem. Phys., vol. 81, p. 511, 1984.
  * (58) W. G. Hoover, “Canonical dynamics: Equilibrium phase-space distributions,” Phys. Rev. A, vol. 31, pp. 1695–1697, 1985.
  * (59) K. Ogawa, H. Oga, H. Kusudo, Y. Yamaguchi, T. Omori, S. Merabia, and L. Joly, “Large effect of lateral box size in molecular dynamics simulations of liquid-solid friction,” Phys. Rev. E, vol. 100, p. 023101, Aug 2019.
  * (60) K. P. Travis, B. D. Todd, and D. J. Evans, “Departure from navier-stokes hydrodynamics in confined liquids,” Phys. Rev. E, vol. 55, pp. 4288–4295, Apr 1997.
  * (61) L. Verlet, “Computer experiments on classical fluids. i. thermodynamical properties of lennard-jones molecules,” Phys. Rev., vol. 159, p. 98, 1967\.
  * (62) J. A. White, “Lennard-jones as a model for argon and test of extended renormalization group calculations,” The Journal of Chemical Physics, vol. 111, no. 20, pp. 9352–9356, 1999.
  * (63) D. M. Heyes, J. J. Kim, C. J. Montrose, and T. A. Litovitz, “Time dependent nonlinear shear stress effects in simple liquids: A molecular dynamics study,” J. Chem. Phys., vol. 73, p. 3987, 1980.
  * (64) P. A. Thompson and M. O. Robbins, “Shear flow near solids: Epitaxial order and flow boundary conditions,” Phys. Rev. A, vol. 41, pp. 6830–6837, Jun 1990\.
  * (65) C. Neto, D. R. Evans, E. Bonaccurso, H.-J. Butt, and V. S. J. Craig, “Boundary slip in newtonian liquids: a review of experimental studies,” Reports on Progress in Physics, vol. 68, pp. 2859–2897, oct 2005.
  * (66) S. K. Kannam, B. D. Todd, J. S. Hansen, and P. J. Daivis, “Slip length of water on graphene: Limitations of non-equilibrium molecular dynamics simulations,” The Journal of Chemical Physics, vol. 136, no. 2, p. 024705, 2012.
  * (67) L. Bocquet and J.-L. Barrat, “Hydrodynamic boundary-conditions, correlation functions, and kubo relations for confined fluids,” Phys. Rev. E, vol. 49, p. 3079, 1994.
  * (68) J. S. Hansen, B. D. Todd, and P. J. Daivis, “Prediction of fluid velocity slip and solid surfaces,” Phys. Rev. E, vol. 84, p. 016313, 2011.
  * (69) S. K. Kannam, B. D. Todd, J. S. Hansen, and P. J. Daivis, “How fast does water flow in carbon nanotubes?,” J. Chem. Phys., vol. 138, p. 094701, 2013.


* [37] O. Sigmund and J. S. Jensen. Systematic design of phononic band-gap materials and structures by topology optimization, Phil. Trans. R. Soc. London, A  361:1001–1019, 2003.
  * [38] V. P. Smyshlyaev. Propagation and localization of elastic waves in highly anisotropic periodic composites via two-scale homogenization, Mechanics of Materials, 41:434–447, 2009
  * [39] J. Vasseur, P. Deymier, G. Frantziskonis, G. Hong, B. Djafari-Rouhani, and L. Dobrzynski. Experimental evidence for the existence of absolute acoustic band gaps in two-dimensional periodic composite media, . Phys. Condens. Matter, 10:6051–6064, 1998.
  * [40] J. Vondrejc, E. Rohan, J. Heczko. Shape optimization of phononic band gap structures using the homogenization approach. International Journal of Solids and Structures. Vol.113-114:147-168, 2017.
  * [41] V.V. Zhikov. On Spectrum Gaps of Some Divergent Elliptic Operators with Periodic Coefficients. St. Petersburg Math. J., 16:773-790, 2005.