Lees-Edwards Boundary Conditions

In the realm of computational physics and chemistry, molecular dynamics simulations are a powerful tool for understanding how matter behaves at the atomic scale. By simulating a small number of particles in a box, we can predict macroscopic properties. For static materials, standard periodic boundary conditions allow a small, finite system to effectively represent an infinite one. However, this simple approach breaks down when we want to model systems in motion, such as a fluid under shear. How can we simulate the continuous, uniform flow of a river or a polymer melt within the confines of a tiny, fixed box? This is the fundamental challenge addressed by the Lees-Edwards boundary conditions.

This article delves into this ingenious computational method. The first chapter, "Principles and Mechanisms," will unpack the core idea behind the Lees-Edwards technique, explaining how it creates a seamless shearing world. We will explore the necessary adjustments to particle dynamics, such as the SLLOD equations of motion, and the critical importance of separating ordered flow from random thermal motion. Subsequently, the "Applications and Interdisciplinary Connections" chapter will demonstrate the method's power, from its primary role in calculating the viscosity of simple and complex fluids to its deeper connections with the fundamental principles of non-equilibrium thermodynamics.

To understand the world, we often build miniature versions of it. In computational physics, our "worlds" are simulation boxes, tiny cubes filled with particles that obey the laws of physics. But how can a tiny, finite box tell us anything about the vast, continuous flow of a river or the slow ooze of honey? The challenge is to make our box behave as if it were a small patch of an infinitely large system. For a static material, the solution is simple: periodic boundary conditions. When a particle exits one side of the box, it magically reappears on the opposite side. It’s as if our box is tiled endlessly in all directions.

But what if the world we want to model is not static? What if it's shearing?

Imagine a simple shear flow, like a river flowing faster at the surface than at the bottom. Or, for a more hands-on picture, think of a deck of cards. If you push the top card sideways, the whole deck deforms. Each card slides a little relative to the one below it. This is the essence of shear. We can describe this motion with a simple, elegant mathematical rule called the ​​streaming velocity​​, where the velocity of the fluid in one direction (say, $x$) increases linearly with the position in a perpendicular direction (say, $y$). Mathematically, we write this as $\mathbf{u}(\mathbf{r}) = \dot{\gamma} y \hat{\mathbf{x}}$, where $\dot{\gamma}$ is a constant called the ​​shear rate​​—it tells us how rapidly the velocity changes with position.

Now, try to put this in our periodic box. Standard periodic boundaries assume the boxes in our infinite tiling are all sitting still relative to each other. But in a shear flow, the layer of fluid just above our box is sliding past it! The layer below is sliding the other way. The entire infinite lattice of periodic images is deforming. Our simple rule of "reappearing on the opposite side" is no longer enough. We need a boundary condition that understands this continuous, sliding motion. This is the brilliant insight of Arthur Lees and S.F. Edwards.

The ​​Lees-Edwards boundary conditions​​ (LEBC) are a marvel of physical intuition and geometric cleverness. The idea is to make the periodic boundaries themselves move in a way that perfectly mimics the shear flow. Imagine again our infinite stack of simulation boxes. At time $t=0$, they are all perfectly aligned. As the shear flow starts, the entire plane of boxes at a height $L_y$ above our central box starts sliding in the $x$-direction. The plane at $2L_y$ slides twice as fast, and so on.

The rule for a particle is simple: when it moves out of the top of the box (crossing the $y=L_y$ boundary), it re-enters at the bottom ($y=0$), but its $x$-position is shifted by exactly the amount the top layer of boxes has slid past the bottom layer in the elapsed time $t$. This displacement is $\Delta x(t) = \dot{\gamma} t L_y$. The same logic applies in reverse for a particle crossing the bottom boundary.

So, the full remapping rule for a particle at an arbitrary position $(x, y, z)$ back into the primary simulation box (let's say, from $0$ to $L$ in each direction) is a two-step process. First, we find which periodic image cell the particle is in. The integer $n_y = \lfloor y/L_y \rfloor$ tells us how many boxes "up" the particle is. We then apply the corresponding shifts to bring it back to the primary cell's frame of reference. The velocity must also be corrected to account for jumping between layers moving at different speeds. The complete transformation looks like this:

• The new position is $\mathbf{r}' = \begin{pmatrix} \left( x - \dot{\gamma} t L_y n_y \right) \pmod{L_x} y \pmod{L_y} z \pmod{L_z} \end{pmatrix}$.
• The new velocity is $\mathbf{v}' = \begin{pmatrix} v_x - \dot{\gamma} L_y n_y v_y v_z \end{pmatrix}$.

This elegant trick creates a seamless, infinitely shearing medium from a single, finite box. There are no walls, no artificial surfaces—just the pure, bulk behavior of a fluid in motion.

In a flowing river, is the water "hot" because the river is moving at several meters per second? Of course not. The temperature of the water relates to the random, chaotic jiggling of individual water molecules, not their collective, ordered motion downstream. This distinction is absolutely critical in a sheared system.

We must decompose a particle's total velocity, $\mathbf{v}_i$, into two parts: the ordered ​​streaming velocity​​, $\mathbf{u}(\mathbf{r}_i)$, and the random, thermal part, which we call the ​​peculiar velocity​​, $\mathbf{c}_i$.

$\mathbf{c}_i = \mathbf{v}_i - \mathbf{u}(\mathbf{r}_i)$

This is not just a mathematical convenience; it is a profound physical statement. The internal energy, the temperature, and the pressure (or more generally, the ​​stress tensor​​) of the fluid are all defined by the statistics of these peculiar velocities. They represent the chaos. The streaming velocity represents the order. The kinetic energy associated with the streaming velocity, $\frac{1}{2}m\mathbf{u}^2$, is macroscopic, ordered energy, not thermal energy.

What would happen if we got this wrong? Suppose we used a thermostat—a device for adding or removing heat to keep the temperature constant—but we mistakenly programmed it to control the temperature based on the total velocity $\mathbf{v}_i$. The streaming velocity $\mathbf{u}(\mathbf{r}_i) = \dot{\gamma} y \hat{\mathbf{x}}$ is largest near the top and bottom of the box ($y=\pm L_y/2$) and zero at the center ($y=0$). Our faulty thermostat would see the high-velocity particles near the boundaries as "too hot" and would aggressively cool them. It would see the particles at the center as "correct." The result? A completely unphysical temperature profile, with the fluid being artificially cold at the boundaries and hot in the middle. The thermostat would be fighting the very flow we are trying to simulate, introducing spurious forces that would even distort the measured stress in the fluid. This "what-if" scenario beautifully illustrates why separating ordered from peculiar motion is not just an option, but a necessity.

To calculate the force on a particle, we need to find all its neighbors within a certain cutoff distance, $r_c$. In a standard periodic box, this is easy. But in our sliding-brick world, a neighboring particle might be in a periodic image that has been shifted significantly sideways. The shortest distance between two particles is no longer found by just looking in the adjacent cells of a fixed grid. How do we solve this geometric puzzle efficiently?

There are two equally beautiful solutions:

1. ​​The Deshearing Trick:​​ Imagine you have a skewed picture. You can apply a transformation in your favorite image editor to "un-skew" it and make it rectangular. We can do the same for our simulation box! At any given moment $t$, our sheared box is a parallelogram. We can apply a simple mathematical transformation (an affine map) that squashes it back into a rectangle. In this "desheared" coordinate system, everything is simple again, and we can find the shortest distance using the standard minimum image convention. Once we have this vector, we apply the inverse transformation to get the true separation vector back in the real, sheared world. It's a clever change of perspective that turns a hard problem into an easy one.

2. ​​The Abstract Matrix View:​​ A more powerful way to think about the problem is to embrace the sheared geometry. The simulation box at any time $t$ is a deforming parallelogram (or, in 3D, a triclinic cell). Any such cell can be described by a matrix, $\mathbf{H}(t)$, whose columns are the vectors that define the cell's edges. For our simple shear, this matrix is: $\mathbf{H}(t) = \begin{pmatrix} L_x \dot{\gamma} t L_y 0 \\ 0 L_y 0 \\ 0 0 L_z \end{pmatrix}$ With this matrix, all geometric operations become elegant matrix-vector algebra. To find the shortest distance between two particles, we transform their separation vector into "fractional coordinates" (coordinates relative to the box vectors) using the inverse matrix, $\mathbf{H}^{-1}(t)$. In this fractional space, we apply the simple nearest-image rule, and then transform back to the real world using $\mathbf{H}(t)$. This approach is not only clean but also general—it can handle any kind of box deformation, not just simple shear.

We have the boundary conditions and the geometry. But what are the actual equations of motion we need to solve for each particle? We can't just use Newton's $m\ddot{\mathbf{r}}_i = \mathbf{F}_i$, because that applies to an inertial (non-accelerating, non-deforming) frame. Our Lees-Edwards world is a non-inertial, deforming frame.

The correct equations of motion are known as the ​​SLLOD​​ equations. They look like Newton's laws but with an extra term that accounts for the deforming reference frame. When written for the peculiar momentum, $\mathbf{p}_i = m\mathbf{c}_i$, the equation for particle $i$ takes the form:

$\dot{\mathbf{p}}_i = \mathbf{F}_i - \dot{\gamma} p_{iy} \hat{\mathbf{x}} - \alpha\mathbf{p}_i$

Here, $\mathbf{F}_i$ is the true inter-particle force, and $\alpha\mathbf{p}_i$ is a thermostatting force to remove viscous heat. The fascinating new term is $-\dot{\gamma} p_{iy} \hat{\mathbf{x}}$. This is a "fictitious force," much like the Coriolis force you feel on a merry-go-round. It arises purely because we are describing the physics relative to a deforming (shearing) background. It is the very term that drives the flow and sustains the shear.

The precise form of this term is deeply significant. One could imagine an alternative, seemingly plausible form, like $-\dot{\gamma} p_{ix} \hat{\mathbf{y}}$ (known as the Doll's tensor formulation). However, this seemingly minor change has a disastrous physical consequence: it violates the conservation of angular momentum. A simulation using this incorrect equation would produce a fluid that has an unphysical internal rotation, resulting in a stress tensor that is not symmetric ($\sigma_{xy} \neq \sigma_{yx}$), which is impossible for a simple fluid. The SLLOD formulation gets the physics of rotation correct, preserving the fundamental symmetries of the underlying mechanics. It is a beautiful example of how the mathematical structure of an algorithm must respect the deep principles of physics.

Stirring a cup of thick honey requires effort. The work you do with the spoon is converted into heat, warming the honey slightly. This phenomenon, called ​​viscous heating​​, is also present in our simulation. The shear flow continuously does work on the fluid, and this work increases the fluid's internal energy—its temperature.

The SLLOD equations and Lees-Edwards boundaries allow us to calculate this process precisely. The rate at which the shear does work on the system, $\dot{W}$, and thus increases its internal energy, is given by a beautifully simple macroscopic-looking formula built from microscopic quantities:

$\dot{W} = -V \sigma_{xy} \dot{\gamma}$

Here, $V$ is the volume of the box, $\dot{\gamma}$ is the shear rate, and $\sigma_{xy}$ is the shear stress—the measure of the fluid's internal friction. Microscopically, this stress is composed of a kinetic part from the motion of particles between layers and a potential part from the forces acting between particles in different layers. This equation is the heart of non-equilibrium molecular dynamics. It connects the microscopic details—the positions, forces, and peculiar momenta of individual particles—to a macroscopic, measurable material property: ​​viscosity​​. By running the simulation, we can measure the stress $\sigma_{xy}$ that results from an imposed shear rate $\dot{\gamma}$, and from their ratio, we can compute the viscosity of our simulated fluid.

From a simple geometric trick of sliding periodic images, we have built a complete, self-consistent world that allows us to probe the fundamental properties of matter far from equilibrium. This journey, from a deck of cards to the calculation of viscosity, showcases the power and beauty of theoretical physics in connecting the microscopic to the macroscopic.

Having grasped the elegant machinery of Lees-Edwards boundary conditions, we can now ask the most important question for any scientific tool: What is it good for? The answer, it turns out, is wonderfully broad. This clever trick of wrapping space in a shearing motion is not merely a computational convenience; it is a veritable portal into the non-equilibrium world, allowing us to probe the behavior of matter under flow with unprecedented clarity. From calculating the simple "stickiness" of water to unraveling the complex thermodynamics of systems far from equilibrium, the applications are a journey through physics, chemistry, and engineering.

At its heart, the Lees-Edwards method provides a perfect, idealized setup to measure a fluid's most familiar transport property: its shear viscosity, the quantity we intuitively understand as thickness or resistance to flow. Imagine trying to measure the viscosity of honey. You might slide a plate over a layer of it and measure the required force. But what about the edges of the plate? What about the interface with the air? These messy details complicate the measurement.

The Lees-Edwards condition is the physicist's dream of this experiment. It creates a state of pure, uniform shear throughout the entire bulk of the simulated fluid, with no walls, no surfaces, and no edges to worry about. By imposing a known, constant shear rate, $\dot{\gamma}$, we can simply "listen" to the fluid's response. This response is the shear stress, $\sigma_{xy}$, which is the internal friction the fluid generates against the imposed flow. In a computer simulation, we can precisely calculate this stress by tracking the motions and interactions of every single particle. The stress arises from two microscopic sources: the physical transport of momentum by particles moving between imaginary layers of the fluid (the kinetic contribution), and the sum of all the intermolecular forces stretching across those layers (the configurational or virial contribution).

The ratio of the measured stress to the imposed shear rate gives us the viscosity, $\eta = -\langle \sigma_{xy} \rangle / \dot{\gamma}$. This method is so powerful and direct that it has become a gold standard for computing the viscosity of everything from simple Lennard-Jones fluids—the physicist's model "atom"—to complex and vital substances like water, whose flow properties are critical to countless biological and chemical processes.

Of course, the world is more interesting than simple honey. Many fluids have a viscosity that changes depending on how fast you try to make them flow. Think of ketchup: it's thick in the bottle, but flows easily once you get it moving. This behavior is called "shear thinning," and it is a hallmark of what we call complex or non-Newtonian fluids.

Here, the Lees-Edwards method truly comes into its own. By systematically varying the shear rate $\dot{\gamma}$ over many orders of magnitude, we can use a series of simulations to map out the entire viscosity function, $\eta(\dot{\gamma})$. This allows us to precisely identify the onset of shear thinning and characterize the full rheological, or flow, behavior of a material. Such studies are not just academic exercises; they are essential in materials science and chemical engineering for designing paints, coatings, foods, and polymer melts. To do this accurately requires great care, controlling for simulation artifacts related to system size and the thermostat used to siphon off heat, but the reward is a complete picture of a material's complex response to stress.

Furthermore, the method is not confined to simple shear. The underlying mathematical framework, known as the SLLOD algorithm, can be generalized to impose any kind of uniform velocity gradient, described by a tensor $\boldsymbol{\kappa}$. This allows us to simulate not just shearing flows, but also elongational flows—the kind of stretching that occurs when you pull on a strand of mozzarella cheese—or any combination thereof. This versatility is crucial for understanding processes like polymer extrusion or the behavior of biological cells in complex blood flow environments.

The Lees-Edwards method does more than just help us engineer materials; it provides a window into the fundamental laws of thermodynamics and statistical mechanics.

First, let's consider the energy. Maintaining a fluid in a state of shear requires a continuous input of work. This work is dissipated by the fluid's internal friction and converted into heat. In a simulation, this "viscous heating" would quickly raise the temperature, so a thermostat is needed to constantly remove the heat and maintain a steady state. This process has a profound thermodynamic consequence: it continuously generates entropy. The Lees-Edwards setup allows us to see this principle in action with beautiful clarity. The rate of mechanical work done on the system, given by $-V \sigma_{xy} \dot{\gamma}$, must exactly equal the heat flow $\dot{Q}$ removed by the thermostat in the steady state. According to the laws of stochastic thermodynamics, this heat flow results in an entropy production rate in the surrounding environment of $\dot{S} = \dot{Q}/T$. Thus, we find a direct, elegant link between the mechanical properties of the flow and a cornerstone of thermodynamics: $\dot{S} = (-V \sigma_{xy} \dot{\gamma}) / T$. The simple act of shearing a fluid is an entropy-generating machine.

Second, the very nature of the SLLOD equations used with Lees-Edwards boundaries is revealing. One might draw an analogy to the motion of a charged particle in a magnetic field, where the equations of motion also contain terms that mix position and velocity components. However, this analogy, while tempting, breaks down in a crucial way. The dynamics of a particle in a magnetic field are Hamiltonian; they conserve energy and phase-space volume. The SLLOD equations are fundamentally non-Hamiltonian and dissipative. Even without interparticle forces, the shear term $-\boldsymbol{\kappa} \cdot \mathbf{p}$ in the momentum equation causes the peculiar kinetic energy to change, a direct manifestation of viscous heating. The phase space of the system continuously contracts, a signature of dissipation that is balanced by the action of the thermostat. This teaches us that simulating non-equilibrium flow requires us to work with a different class of dynamics, one where the fundamental conservation laws of the equilibrium world no longer hold in the same way.

Finally, this non-equilibrium method connects back to the equilibrium world through one of the deepest principles in physics: the Fluctuation-Dissipation Theorem. This theorem states that the way a system responds to an external kick (dissipation) is intimately related to the natural, spontaneous fluctuations it exhibits at rest. For viscosity, this means that the value we measure by forcing the fluid to flow with Lees-Edwards boundaries (a non-equilibrium measurement) can also, in principle, be calculated from the correlations of random, fluctuating stresses in a fluid at complete rest (an equilibrium measurement using the Green-Kubo formula). Comparing these two methods provides a powerful consistency check and a profound demonstration of the unity of statistical physics.

Lastly, the Lees-Edwards method serves as a bridge between the microscopic world of atoms and the macroscopic world of continuum hydrodynamics. A simulation box is finite, and this has consequences. The allowed patterns of fluid motion—the hydrodynamic modes—are quantized by the box's finite dimensions, much like the sound waves on a guitar string are limited to specific harmonics. This discretization of modes introduces a subtle but important finite-size effect on the calculated viscosity. By carefully analyzing how the measured viscosity changes with the size and shape of the simulation box, we can test and validate theories of fluctuating hydrodynamics, learning how the continuum description of a fluid emerges from its discrete, atomic constituents.

In a sense, the Lees-Edwards boundary condition is the perfect microscope for the physics of flow. It strips away all the non-essential, complicating factors of a real-world experiment and allows us to focus on the intrinsic, bulk behavior of matter. By doing so, it not only gives us the power to compute practical numbers for engineering but also to see the beautiful and deep connections that link mechanics, thermodynamics, and statistical physics in the dynamic world far from equilibrium.