Physics Simulation – ETHermodynamics

Quickstart

- This article explains the process of simulating an ETHermodynamics experiment. You can check it out if you have questions about the underlying mathematics of ETHermodynamics.
- The processes described are ordered in chronological order, i.e. in the order they are executed in the simulation.
- The physics simulation described here is the default one used for most ETHermodynamics systems. The one-dimensional chamber uses a slightly different approach that is described in the Chambers: 1D-Chamber article.

The System

Chamber

The purpose of the ETHermodynamics software is to estimate values of thermodynamic parameters for a macroscopic system of $i = 1, \ldots, N$ point particles with masses $m_\text{i}$ that are enclosed in a rectangular container $[0, L_\text{x}] × [0, L_\text{y}] × [0, L_\text{z}] $ with volume ${V} = L_\text{x} L_\text{y} L_\text{z}$.

Particles

Particles are point particles that have two intrinsic properties, the mass $m_\text{i}$ and type (or class) $c_\text{i}$. The particle type determines their behavior in collisions with permeable walls, chemical reactions, phases, etc. To simplify the notation, corresponding type-specific parameters are treated as particle-specific, i.e., denoted as $p_\text{i}$ instead of $p(c_\text{i})$ for particle $i$.

Piston

One wall can be a movable piston with mass ${M}_\text{p}$, on which an external pressure ${P}_\text{ext}$ is exerted. The system can interact with its surroundings by heat transfer and/or mechanical work in the form of piston movement.

Multiple Chambers & Connector

The system can be extended to a system with two rectangular chambers that are separated along the vertical axis by a fixed or movable wall (connector). This can be an adiabatic wall that does not allow heat or particles to pass, a diathermal wall that allows exchange of heat but not of particles, and permeable walls that allow the transfer of particles, or heat and particles. The thickness of this wall is assumed to be negligible (zero). The sizes ${L_\text{x}}, {L_\text{y}}, {L_\text{z}}$ and particle numbers ${N}$ that characterize the geometry and content of a chamber must be specified for each chamber. Thermodynamic quantities can be obtained for each chamber.

Scaling

To make computations feasible, the time evolution of a much smaller number n of particles is simulated explicitly. Results for the macroscopic system of N particles are obtained by scaling, assuming that for each explicitly simulated particle there are $\lambda = N/n$ particles moving identically.

Typically, N is of the order of 10²³, where n is in the range of 10⁴–10⁵. The scaling is implemented by multiplying corresponding quantities by the scaling ratio $\lambda$. Note that the dimensions of the container and the time are not scaled. The same scaling ratio is used for all chambers of the system.

As a result of scaling, fluctuations of thermodynamic parameters are overestimated by approximately a factor of $\sqrt{\lambda}$. Fluctuations are therefore expected to be about 10⁹ times larger than in reality, but still small.

Setup

Particles are initially placed randomly in the container at coordinates

$$r_{\text{ij}} = L_{j}u_{\text{ij}}$$

for particles $i = 1, \ldots, N $ and dimensions $j = 1, 2, 3 (x , y , z)$, where $u_{ij}$ are random numbers distributed uniformly in the interval $[0,1]$. Velocities are initialized such that for a given initial temperature $T$ the total kinetic energy is exactly $E_{\text{kin}} = \frac{3}{2}nk_{B}T$, with $n$ the number of particles and $k_B$ the Boltzmann constant. To this end, the velocity components $v_{ij}$ are initialized to

$$ v_{ij} = \sqrt{\frac{3nk_{B}T}{Gm_{i}}}g_{ij} $$

where $g_{ij}$ are normally distributed random numbers with mean 0 and standard deviation 1, and $G = \sum_{i = 1}^{n}{\sum_{j = 1}^{3}g_{ij}^{2}}$. This distribution of velocities is not the Boltzmann distribution that one expects in thermal equilibrium, but the system will in general quickly evolve towards equilibrium.

Time evolution

The simulation proceeds in $N_{\text{steps}}$ time steps of fixed duration $\Delta t$ for a total simulation time $t_{\text{total}} = N_{\text{steps}}\Delta t$. The time step should be chosen such that, in one step, no particle $i$‘s coordinate changes by more than the corresponding container dimension, $L_j$, i.e. with $|v_{ij}|$ the particle’s speed:

$$ | \Delta x_{ij} | =|v_{ij}|\Delta t < L_j $$

for all i and j. At normal conditions, time steps can be in the range of $ \Delta t = 10^{-5}-10^{-4} \text{ s} $.

These time steps are much longer than in conventional molecular dynamics simulation, where strong interactions between atoms dictate time steps of about $2\text{ fs}=2\cdot10^{-15}\text{ s}$. These much longer time steps allow the simulation of a thermodynamic system over several seconds on a laptop computer.

Each time step $ t \rightarrow t + \Delta t$ comprises several actions, which are executed in the following order:

1. Stage transition
2. Free time evolution of particles
3. Collisions of particles with container walls
4. Heat exchange with surrounding reservoir
5. Energy exchange between particles
6. Heat exchange through the connector
7. Particle exchange through the connector
8. Chemical reactions
9. Storage of quantities for analysis

A complete simulation can be composed of multiple stages. In each stage, the external system parameters have prescribed fixed values (or, in the case of the external pressure, follow a prescribed schedule; see below).

Stage transition

Changes of the temperature of the external thermal bath are applied immediately since the system in general adapts quickly to the new condition.

On the other hand, abrupt changes of the external pressure on the piston can lead to undesired oscillations of the piston position. To avoid this, the external pressure can be changed smoothly from a given initial value $P_\text{ext}(0)$ to a new value ${P}_\text{ext}$ during a simulation stage.

In this case, the external pressure is adapted quadratically during the first three quarters of the simulation stage, and stays constant at ${P}_\text{ext}$ during the last quarter:

$$ P_{\text{ext}}\left( t \right) = \left\{ \begin{matrix}
P_{\text{ext}} – (P_{\text{ext}} – P_{\text{ext}}\left( 0 \right))\left( 1 – \frac{t}{\frac{3}{4}t_{\text{tot}}} \right)^{2} && \ \ \ \ \ \ \ \text{if}\ 0 \leq t < \frac{3}{4}t_{\text{tot}} \\
P_{\text{ext}} && \text{if}\ t \geq \frac{3}{4}t_{\text{tot}} \\
\end{matrix} \right. $$

Free time evolution of particles

In the absence of forces, particles move linearly:

$${\widetilde{r}}_{\text{ij}}\left( t + \mathrm{\Delta}t \right) = r_{\text{ij}}\left( t \right) + v_{\text{ij}}\left( t \right)\mathrm{\Delta}t$$

Velocities remain unchanged, and particles may temporarily move outside the container before they are returned into the container in the next action.

Collisions of particles with container walls

This action handles all collisions of particles with the fixed container walls and with the moving piston within a time step. The collisions are elastic. For efficiency, the collisions are not handled in the time order in which they occur during the time step but in a random order of the particles.

This treatment is exact for collisions of particles with fixed walls but incurs an approximation for collisions of particles with the piston that changes its velocity during the time step due to the impact of particles.

For each particle ${i} = \pi_\text{1}, \ldots, \pi_\text{n}$, where ${\pi}$ is a random permutation of $1,\ldots, n$, it is checked for each coordinate dimension ${j} = 1, 2, 3 $, whether the particle has left the container during the time step.

Particle ${i}$ has collided with a wall in dimension ${j}$ if ${\widetilde{r}}_{\text{ij}}\left( t + \mathrm{\Delta}t \right) < 0$ or ${\widetilde{r}}_{\text{ij}}\left( t + \mathrm{\Delta}t \right) > L_{j}$.

Its position then changes as

$$ r_{\text{ij}}\left( t + \mathrm{\Delta}t \right) = \left\{ \begin{matrix}
– {\widetilde{r}}_{\text{ij}}\left( t + \mathrm{\Delta}t \right) & \text{if}{\ \widetilde{r}}_{\text{ij}}\left( t + \mathrm{\Delta}t \right) < 0\ \\
2L_{j} – {\widetilde{r}}_{\text{ij}}\left( t + \mathrm{\Delta}t \right) & \text{if}{\ \widetilde{r}}_{\text{ij}}\left( t + \mathrm{\Delta}t \right) > L_{j} \\
\end{matrix} \right.\ $$

If the wall is fixed, the velocity of the particle is reversed:

$$ v_{\text{ij}}\left( t + \mathrm{\Delta}t \right) = – v_{\text{ij}}\left( t \right) $$

Note that within a time step a particle may undergo up to 3 collisions with walls in different dimensions, but not multiple collisions with two walls of the same dimension because we assume that the time step $\Delta t$ has been chosen short enough to exclude this possibility, i.e. $ \Delta t < L_\text{j}/|v_\text{ij}|$ for all ${i}$ and ${j}$.

If particle ${i}$ hits the movable piston at position ${z} = {L}_\text{z}({t})$, the ${z}$-components of the velocities of the particle ${v_\text{iz}}$ and the piston ${v}_\text{p}$ change according to the laws of momentum and energy conservation:

$$ v_{\text{iz}}\left( t + \mathrm{\Delta}t \right) = \frac{M_{i}v_{\text{iz}}\left( t \right) + M_{p}(2{\widetilde{v}}_{p} – v_{\text{iz}}\left( t \right))}{M_{i} + M_{p}} $$

$$ v_{p} = \frac{M_{p}{\widetilde{v}}_{p} + M_{i}(2{v_{\text{iz}}\left( t \right) – \widetilde{v}}_{p})}{M_{i} + M_{p}} $$

where $ {\widetilde{v}}_{p}$ and $v_{p} $ are the velocities of the piston before and after its collision with particle ${i}$, ${M_\text{i}} = {\lambda m_\text{i}}$ is the scaled-up mass of particle ${i}$, and ${M}_\text{p}$ the mass of the piston.

For the first collision in the time step, ${\widetilde{v}}_{p}$ is initialized to ${\widetilde{v}}_{p} = v_{p}(t)$. Then all collisions are treated in random order, iteratively updating $v_{p}$.

The pressure exerted on the container walls by the particles is obtained from the momentum change in the collisions:

$$ P = \frac{1}{A\mathrm{\Delta}t}\sum_{i = 1}^{n}{M_{i}\left| v_{\text{ij}}\left( t + \mathrm{\Delta}t \right) – v_{\text{ij}}\left( t \right) \right|}$$

where ${A} = 2({L_\text{x}L_\text{y}} + {L_\text{x}L_\text{z}} + {L_\text{y}L_\text{z}})$ is the inner surface of the container.

The external pressure ${P}_\text{ext}$ leads to a force $ {F}_\text{ext} = {P}_\text{ext} L_\text{x} L_\text{y} $ on the piston, and hence to an acceleration ${a} = -{F}_\text{ext~}/M_{p}$ of the piston. The negative sign results from the fact that the external force acts against an increase of the piston position ${L_\text{z}}$.

The position of the piston is changed only at the end of the time step as a result of the collisions and the external pressure:

$$ L_{z}\left( t + \mathrm{\Delta}t \right) = L_{z}\left( t \right) + v_{p}\mathrm{\Delta}t + \frac{1}{2}a{\mathrm{\Delta}t}^{2} $$

where $v_{p}$ is the velocity obtained by the last collision. The final velocity of the piston at the end of the time step is

$$ v_{p}\left( t + \mathrm{\Delta}t \right) = v_{p} + a\mathrm{\Delta}t $$

The work done by the system in the time step results from the motion of the piston against the external pressure and the change of the kinetic energy of the piston:

$$ \mathrm{\Delta}W = {- F}_{\text{ext}}{(L}_{z}\left( t + \mathrm{\Delta}t \right) – L_{z}\left( t \right)) – \frac{1}{2}\text{\ M}_{p}\left( {v_{p}\left( t + \mathrm{\Delta}t \right)}^{2} – {v_{p}\left( t \right)}^{2} \right) $$

The sign follows the convention that energy entering the system is counted positive.

Heat exchange with surrounding reservoir

Coupling to a thermal bath is implemented by adapting the kinetic energy of particles in the vicinity of the container walls (including the piston) to a Boltzmann distribution corresponding to the temperature ${T}_\text{ext}$ of the thermal bath.

It is important to apply the energy exchange to all particles in a small layer next to the container walls rather than only to the particles that collide with the walls because the former acts uniformly on particles of any velocity whereas the latter would affect fast moving particles more often than slow ones (since fast particles collide more often with the walls).

The strength of coupling with the thermal bath is controlled by the parameter ${r}_\text{heat}$ (with unit $s^{-1}$) that specifies how often, on average per unit time, a particle undergoes an energy exchange with the thermal bath.

In a time step, all particles ${i}$ that are closer than $\Delta {L} = {r}_\text{heat} \Delta {t V}/{A} $ to a wall (a good approximation as long as $\Delta {L} \ll {L_\text{x}}, {L_\text{y}}, {L_\text{z}}$), have their velocities scaled

$$ v_{\text{ij}} \longmapsto \sqrt{\frac{E_{kin,i}^{(new)}}{E_{kin,i}}}v_{\text{ij}} $$

where $E_{kin,i} = \frac{1}{2}m_{i}\sum_{j = 1}^{3}v_{\text{ij}}^{2} $ is the kinetic energy of particle ${i}$ before scaling, and $E_{kin,i}^{(new)} = \frac{3}{2}k_{B}T_{\text{ext}}\left( – \log u_{i} \right) $, with $u_\text{i}$ a random number distributed uniformly in the interval $(0,1]$, is the new kinetic energy.

Since the random numbers $-\log{u_\text{i}}$ are exponentially distributed in $ [0,\infty)$, $ E_{kin,i}^{(new)}$ follows a Boltzmann distribution corresponding to the temperature ${T}_\text{ext}$ of the thermal bath.

The heat transferred into the system in a time step is

$$ \mathrm{\Delta}Q = \frac{N}{n}\sum_{i}^{}{E_{kin,i}^{(new)} -}E_{kin,i} $$

where the sum runs over all particles ${i}$ that have exchanged energy with the thermal bath.

Energy exchange between particles

Even in an ideal gas, a (weak) mechanism that allows energy to be exchanged between particles is necessary to establish thermodynamic equilibrium. In the ETHermodynamics program, energy exchange between particles is implemented by selecting random pairs of particles whose velocities are reset randomly but such that the total kinetic energy of each particle pair is preserved. This approach is unphysical but requires far less computation time than, for instance, a pair potential.

The rate of energy exchange between particles is controlled by the parameter ${r}_\text{mix}$ (with unit $\text{s}^\text{-1})$ that specifies how often, on average per unit time, a particle undergoes an energy exchange with another particle. In one time step, ${r}_\text{mix} \Delta t/2$ random particle pairs are selected.

The velocities of the particles ${i}$ and ${k}$ in such a pair are then reset to random values such that

their kinetic energies are consistent with a Boltzmann distribution and
the total kinetic energy of the pair of particles is not changed.

To fulfil the first condition, new velocities for the two particles ${l} = {i}, {k}$ are initialized (on an arbitrary scale) as

$$ {\widetilde{v}}_{\text{lj}} = \sqrt{\frac{- \log u_{l}}{m_{l}g_{l}^{2}}}g_{\text{lj}} $$

where ${u_\text{l}}$ are random numbers distributed uniformly in the interval $(0,1]$, ${g_\text{lj}}$ are normally distributed random numbers with mean 0 and standard deviation 1, and $g_{l}^{2} = \sum_{j = 1}^{3}g_{\text{lj}}^{2}$. To keep the total energy for the particle pair constant, the velocities of the two particles ${l} = {i}, {k}$ are then scaled by

$$ v_{\text{lj}} \longmapsto \sqrt{\frac{m_{i}v_{i}^{2} + m_{k}v_{k}^{2}}{m_{i}{\widetilde{v}}_{i}^{2} + m_{k}{\widetilde{v}}_{k}^{2}}}{\ \widetilde{v}}_{\text{lj}} $$

where $ v_{i}^{2} = \sum_{j = 1}^{3}v_{\text{ij}}^{2} $ etc.

Heat Exchange through the connector

Heat exchange through the connector is implemented in the same way as the energy exchange between particles treated above, except that it is applied only to particle pairs close to the connector, such that the two particles reside on different sides of the wall.

In analogy to the heat exchange with a surrounding reservoir, the particles that exchange energy are chosen such that they are closer than $\Delta_{L} = r_\text{heat} \Delta{t V}/A$ to the connector, where $V$ is the volume and ${A}$ the inner surface of the chamber that holds the particle ($\Delta {L}$ varies for chambers with different ${V}/{A}$ ratio).

The parameter ${r}_\text{heat}$ (with unit $\text{s}^{-1}$), introduced above for the heat exchange with a surrounding reservoir, specifies how often, on average per unit time, a particle undergoes an energy exchange with a thermal bath or a diathermal wall if all walls were heat-conducting.

If ${n}_\text{1} ≤ n_\text{2}$ are the numbers of particles in the two chambers that are close enough to the diathermal wall to exchange energy, then ${n}_\text{1}$ pairs ($i, \pi _\text{i}$) are formed, where ${i} = 1,\ldots, {n}_\text{1}$ denotes the index of a particle near the diathermal wall of the one chamber, and $\pi_\text{i}$, which is the ${i}$-th element of a random permutation of $1,\ldots, {n}_{2}$, denotes the index of a particle near the diathermal wall of the other chamber.

The velocities of the two particles in each of these pairs are then changed assuming a hard spheres collision as described for the energy exchange between particles.

The heat transferred into a given chamber in a time step is

$$ \mathrm{\Delta}Q = \frac{N}{n}\sum_{i}^{}{E_{kin,i}^{(new)} -}E_{kin,i} $$

where the sum runs over all particles $i$ in the chamber that have exchanged energy through the diathermal wall.

Particle exchange through the connector

A permeable connector is characterized by a particle-specific reflection probability ${p_\text{i}} = p \cdot p_{\text{modifier}, i}$, where $p$ is the connector-permeability (identical for all particles, but variable over stages) and $p_{\text{modifier}, i}$ is a particle-type specific parameter.

A particle $i$ that collides with a permeable wall is randomly either reflected with probability ${p_\text{i}}$ as described above for fixed or movable walls, or continues with probability $1-{p_\text{i}}$ its linear trajectory as if no wall would have been present.

Chemical reactions

Introduction

Chemical reactions change the type of particles. A general chemical reaction

$$ R: A_1 + A_2 + \dots \ \rightarrow B_1 + B_2 + \dots $$

converts reactants $A_1, A_2, \dots$ into products $B_1, B_2, \dots$ Without loss of generality, stochiometric factors can be handled by multiple identical reactants or products.

Implementation

The implementation of reactions consists of two parts: First, to determine which particles take part in a reaction in a given time step. Second, to carry out the reaction by replacing the reactants with the products and computing their positions and velocities.

Step 1

rate law

The first part can be implemented phenomenologically by assuming that the concentrations of reactants $[A_i]$ and products $[B_j]$ change according to a rate law of the form

$$ – \frac{d\left\lbrack A_{i} \right\rbrack}{dt} = \frac{d\left\lbrack B_{j} \right\rbrack}{dt} = k_{R}\prod_{k}^{}\left\lbrack A_{k} \right\rbrack^{\alpha_{k}} $$

The concentrations are defined as $[A_i] = A_i/(N_AV)$, where N_i is the (scaled) number of particles of type $[A_i]$ (or $[B_i]$) in the chamber of volume $V$ in which the reaction takes place, and $N_A$ the Avogadro constant.

$k_R$ is the reaction rate, and the exponents $\alpha_k$ specify the dependence of the reaction rate on the corresponding reactant concentration.

Typically, the exponents $a_k$ are small integer or fractional numbers (e.g., 0, ½, 1, 2). The index $k$ of the product runs over all reactants.

Number of particles

The (unscaled) number of simulated particles $n_i$ of type $[A_i]$ is related to the molar concentration by $[A_i] = n_i/(n_AV)$, where$n_A = N_A\lambda$ is Avogadro’s constant divided by the overall scaling ratio $\lambda = N/n$.

From the rate law, the number of reactions $\Delta n_i$. that take place in one time step of duration $\Delta t$ can be obtained:

In one time step, $\Delta n_i$ particle tuples $A_1, A_2, \dots$ are selected at random. For simplicity, it is only required that the reactants are in the same chamber but not that they are close together in space.

Microscopic rate law

Alternatively, in a more microscopic view, the macroscopic rate law can be replaced by the assumption that the reaction takes place with an intrinsic reaction rate (probability per unit time) $r_\text{R}$ whenever at least one molecule of each reactant is within a small “reaction region of volume $V_\text{R}$.

To this end, the chamber of volume $V$ is assumed to be composed of non-overlapping small reaction volumes $V_\text{R}$ in the form of a grid. In each time step, each grid cell is checked for the simultaneous presence of all reactants.

This approach corresponds to a special case of the above phenomenological rate law, which can be seen as follows. If the particles are distributed uniformly within the chamber, the probability $p_\text{i}$ to find at least one reactant particle of type $A_\text{i}$ in a given region of volume $V_\text{R}$ is

$$ p_{i} = 1 – \left( \frac{V – V_{R}}{V} \right)^{n_{i}} \approx 1 – \left( 1 – n_{i}\frac{V_{R}}{V} \right) = n_{i}\frac{V_{R}}{V} = \frac{V_{R}}{n_{A}}\left\lbrack A_{i} \right\rbrack $$

The probability that a reaction takes place within the given region of volume $V_\text{R}$ during a time $\Delta t$ then becomes

$$ p = r_{R}\Delta t\prod_{k}^{}p_{k} = r_{R}\Delta t\left( \frac{V_{R}}{n_{A}} \right)^{r}\prod_{k}^{}\left\lbrack A_{k} \right\rbrack $$

where $r$ is the number of reactants and the product runs over all reactants. Since there are $V/V_\text{R}$ cells of volume $V_\text{R}$, the total number of reactions that take place in one time step in the chamber of volume $V$ is

$$ \Delta n_{i} = {p\frac{V}{V_{R}} = r_{R}\Delta t\frac{V}{V_{R}}\left( \frac{V_{R}}{n_{A}} \right)}^{r}\prod_{k}^{}\left\lbrack A_{k} \right\rbrack $$

Macroscopic Rate Law

Comparing this with the phenomenological rate law relates the macroscopic reaction rate constant $k_\text{R}$ to the microscopic intrinsic reaction rate $r_\text{R}$ and the reaction volume $V_\text{R}$,

$$ k_{R} = r_{R}\frac{V_{R}^{r – 1}}{n_{A}^{r + 1}} $$

and yields $\alpha_\text{k} = 1 $ for all exponents in the phenomenological rate law.

Step 2

When a reaction takes place, the reactants are replaced by products of different type.

Each product particle is placed at the position of a former reactant particle (a reactant position may be used multiple times if there are more product particles than reactant particles), and its velocity is reset using the same method as for the energy exchange between particles (see above).

The total kinetic energy of the product particles becomes $E_\text{kin,B} = E_\text{kin,A} + E_\text{R}$, the sum of the total kinetic energies of the reactants and the energy released by the reaction, $E_\text{R}$, which is positive for an exothermal reaction and negative for an endothermal reaction.

The required generalization of the velocity reset method to more than two particles and for $E_\text{R} ≠ 0$ is straightforward.

Storage of quantities for analysis

Position-dependent entropy

An entropy ${S_\text{q}}$ of the spatial distribution of particles can be defined using a spatial grid by

$$ S_{q} = k_{B}\log{\Omega_{q}} \ \ \ \ \ \ \ \text{with} \ \ \ \ \ \ \ \Omega_{q} = \frac{N!}{N_{1}! N_{2}! \ldots} $$

$\Omega _\text{q}$ is the multiplicity of the macrostate characterized by the grid cell occupancies ${N}_\text{1}, {N}_\text{2},\ldots $, i.e. the number of ways in which ${N}$ particles can be distributed in the system such that grid cell ${i}$ is occupied exactly ${N_\text{i}}$ times.

${S}_\text{q}$ is maximal if the particles are distributed evenly over all accessible grid cells. The entropy ${S}_\text{q}$ is a state function that can be calculated instantaneously for any given microscopic state of the system by counting the occupancies of the grid cells.

With the approximation $ log {n}! ≈ {n} log {n} – {n} $ that is good for large ${n}$, the entropy can be calculated as

$$ S_{q} = – Nk_{B}\sum_{i}^{}p_{i}\log p_{i} $$

where the sum runs over all occupied grid cells, ${p_\text{i}} = {N_\text{i}}/{N}$ is the fractional occupancy of grid cell ${i}$ with ${N_\text{i}}$ the number of particles in cell ${i}$ and ${N}$ the total number of particles.

This formula may be generalized to cells of different size by redefining ${p_\text{i}}$ as ${p_\text{i}} = {\rho _\text{i}}/{\rho} / {M}$, where ${\rho_\text{i}} = {N_\text{i}}/{V_\text{i}}$ and ${\rho} = {N}/{V}$ are the particle densities for grid cell ${i}$ with volume ${V_\text{i}}$ and the entire system with volume ${V}$, and ${M}$ denotes the number of (accessible) cells.

The two definitions of ${p_\text{i}}$ are equivalent if all cells are of the same size (${V_\text{i}} = {V}/{M}$). Cells with different size can occur, for instance, at the boundaries of chambers.

Since the entire grid is imaginary, the cells and their sizes are in principle arbitrary, except that we assume them to be macroscopic in the sense that the number of particles in a given cell can, at least in principle, be measured. Scaling the volumes of all cells by a constant factor ${\gamma}$ results in an entropy

$$ S_{q}^{(\gamma)} = – Nk_{B}\sum_{j = 1}^{M/\gamma}{{\gamma p}_{j}\log{\gamma p}_{j}} \approx {- Nk_{B}\sum_{i = 1}^{M}{p_{i}\log{\gamma p}_{i}} = S}_{q} – Nk_{B}\log\gamma $$

that differs from the entropy ${S}_\text{q}$ based on the original grid only by an additive constant. Therefore, it is also permitted to calculate the entropy with ${p_\text{i}} = {\rho_\text{i}}/{\rho} $ instead of ${p_\text{i}} = {\rho_\text{i}}/{\rho} /{M}$. For the approximations to hold well, the cells should be large enough to contain $N_{i} \gg 1$ particles.

As a simple sum over the grid cells, the entropy ${S}_\text{q}$ is additive in the sense that the contributions from different spatial regions (for instance, different chambers, or parts of chambers), simply sum up. Likewise, the entropy of a system with particles of different types is the sum of the entropies calculated for each type of particles.

The entropy ${S}_\text{q}$ can change in two ways.

In the first case, a closed system with fixed thermodynamic properties, such as chamber size (${V}$), particle number (${N}$), total energy (${E}$) etc., that is in an arbitrary microscopic state will evolve towards thermodynamic equilibrium, a macroscopic state in which all thermodynamic quantities (but not the microscopic degrees of freedom of the particles) remain stationary.

Assuming that the system is with equal (or similar) probability in any of its accessible microscopic states, it will eventually stay almost always in or close to the macroscopic state with the largest multiplicity and thus largest entropy. The maximum of the multiplicity and entropy is exceedingly sharp for a system with a macroscopic number of particles.

In this first case, the given thermodynamic quantities ${N}, {V}, {E},\ldots{}$ have the same values for all accessible microscopic states, and changes of the spatial entropy ${S}_\text{q}$ do not reflect (non-existent) changes of the fixed thermodynamic quantities ${N}, {V}, {E},\ldots{}$ but evaluate a different macroscopic quantity, namely the spatial distribution of the particles, which is expressed by the set of occupation numbers of the imaginary cells.

In the second case, we consider changes of the entropy between equilibrium states with different thermodynamic properties.

For example, in an adiabatic expansion starting from thermodynamic equilibrium, the volume ${V}$ will increase and the internal energy ${E}$ will decrease, leading to a new equilibrium state with different values of the thermodynamic parameters. Entropy changes for such processes are described, for instance, by the Sackur-Tetrode formula for the entropy of an ideal gas

$$ S\left( E,V,N \right) = Nk_{B}\left( \log\frac{V}{N} + {\frac{3}{2}\log}\frac{E}{N} + X \right) $$

where ${X}$ is an arbitrary constant. The first term, $ {Nk}_\text{B}\log{{V}/{N}}$, corresponds to the spatial entropy ${S}_\text{q}$.

For processes connecting two equilibrium states, $A \rightarrow B$, the change of the microscopically defined entropy ${S}_\text{q}$ (plus the analogous momentum-dependent entropy ${S}_\text{m}$ that will be introduced below) should be equivalent to the change of an entropy defined in macroscopic thermodynamics on the basis of the heat transfer in a reversible process, $\Delta {Q}_\text{rev}$,

$$ \mathrm{\Delta}S = S_{B} – S_{A} = \int_{A}^{B}\frac{{\mathrm{\Delta}Q}_{\text{rev}}}{T} $$

The equivalence of the microscopic and thermodynamic definitions of entropy can be confirmed in a simulation experiment that records $ \Delta {Q}/{T}$ for a process connecting two equilibrium states ${A}$ and ${B}$, in which the external thermodynamic parameters ${N}, {V}, {E},\ldots{}$ change sufficiently slowly, such that the system proceeds from ${A}$ to ${B}$ through a series of equilibrium states for the current values of the external thermodynamic parameters ${N}, {V}, {E},\ldots{}$

Energy-dependent entropy

In analogy to the entropy of the spatial distribution of the particles, an entropy of the energy distribution of the individual particles can be defined using a one-dimensional energy grid by

$$ S_{e} = k_{B}\log{\Omega_{e}} \ \ \ \ \ \ \ \text{with} \ \ \ \ \ \ \ \Omega_{e} = \frac{N!}{N_{1}!\ N_{2}!\ldots} $$

$ \Omega_\text{e} $ is the multiplicity of the macrostate characterized by the energy grid cell occupancies ${N}_{1}, {N}_{2},\ldots$, i.e. the number of ways in which ${N}$ particles can be distributed in the system such that energy grid cell ${i}$ is occupied exactly ${N_\text{i}}$ times.

The entropy ${S}_\text{e}$ is a state function that can be calculated instantaneously for any given microscopic state of the system by counting the occupancies of the energy grid cells.

With the approximation $ \log{{n}!} ≈ {n} \log{n} – {n}$ that is good for large ${n}$, the entropy can be calculated as

$$ S_{e} = – Nk_{B}\sum_{i}^{}p_{i}\log p_{i} $$

where the sum runs over all occupied energy grid cells, ${p_\text{i}} = {N_\text{i}}/{N}$ is the fractional occupancy of energy grid cell ${i $ with ${N_\text{i}}$ the number of particles in cell ${i}$ and ${N}$ the total number of particles.

The energy distribution is subject to the constraint that the energies of the particles must sum up to the total energy ${E}$.

Under this constraint, ${S}_\text{e}$ will be maximal if the particles are distributed over the accessible grid cells approximately according to a Boltzmann distribution (approximate because the occupancy of grid cells with energy larger than the total energy ${E}$ is strictly zero).

The energy grid may be equidistant with spacing $\Delta {E}$, such that grid cell ${i} = 1,2,\ldots{}$ covers the energy interval $[({i} — 1)\Delta {E}, {i}\Delta {E})$.

The energy spacing $ \Delta {E}$ must be held constant in order to obtain comparable entropy values. As in the case of the positional entropy, a change in $\Delta {E}$ results in an additive constant added to ${S_\text{e}}$.

In principle, at least ${E}/\Delta{E}$ grid cells are necessary to accommodate all possible energy distributions, including the most extreme distribution with all but one particle at rest. However, the higher energy cells are very unlikely to be occupied by “normal” distributions.

In good approximation, the energy grid cells can therefore be limited to cover a fixed multiple (for instance, 10) of the average energy per particle, $3/2{k}_\text{B} {T}_\text{max}$, where ${T}_\text{max}$ is the highest temperature expected in the experiment. The few particles with higher energy will be counted in the highest energy cell.

Momentum-dependent entropy

The purely energy-based entropy ${S_\text{e}}$ does not take into account the directions in which particles fly.

For instance, the very unlikely situation that all particles fly parallel to the ${x}$-axis would have the same energy-based entropy ${S}_\text{e}$ as the natural case of random flight directions.

This aspect can be captured by defining an entropy ${S_\text{m}}$ that is based on the distribution of momenta (proportional to velocities) instead of particle energies. This entropy is implemented in analogy to the spatial entropy using a three-dimensional momentum grid.

$$ S_{m} = – Nk_{B}\sum_{i}^{}p_{i}\log p_{i} $$

where the sum runs over all occupied momentum grid cells, ${p_\text{i}} = {N_\text{i}}/{N}$ is the fractional occupancy of grid cell ${i}$ with ${N_\text{i}}$ the number of particles in cell ${i}$ and ${N}$ the total number of particles.

To keep particles of different mass on the same scale, the grid is not applied directly to the momenta ${p} = {mv}$ but to $\widetilde{p} = \sqrt{m/2}{v}$, which is on the scale of $\sqrt{E}$ and therefore similar for particles of any mass if one expects equipartition of energy.

The grid size of the this scaled momenta grid is set to a fixed value that corresponds to a temperature difference of 50 K.

Quantity storage

At the end of each time step, quantities of interest are recorded for analysis. These include, but are not limited to:

Volume ${V} = {L_\text{x} L_\text{y} L_\text{z}}$
Temperature $T = \frac{\sum_{i = 1}^{n}{\frac{1}{2}m_{i}v_{i}^{2}}}{\left( \frac{3}{2}nk_{B} \right)}$
Pressure ${P}$
External pressure ${P}_\text{ext}$
Piston velocity ${v}_\text{p}$
Internal energy $U = \frac{N}{n}\sum_{i = 1}^{n}{\frac{1}{2}m_{i}v_{i}^{2}} = \frac{3}{2}Nk_{B}T$
Heat transfer into the system $\Delta {Q}$
Total heat transferred into the system $Q \longmapsto Q + \Delta Q$
Work done by the system ${W}$
Total work done by the system $W \longmapsto W + \Delta W$
Entropy, composed of spatial and momentum-related terms $ {S} = {S}_\text{q} + {S}_\text{m} $
Ideality: agreement with the ideal gas law ${PV}/({Nk}_\text{B}{T})$ (=1 for an ideal gas)

Cumulative quantities are initialized at the beginning of the simulation: ${Q} = 0, {W} = 0$. Recording of heat and work can be switched on and off for individual stages.

For reporting, most quantities ${X}({t})$ are averaged over a time period that is equal to the interval between experiment snapshots (this setting can be changed in the Scene tabs).

$$ \overline{X}\left( t \right) = \frac{1}{n_{\text{ave}}}\sum_{i = 0}^{n_{\text{ave}} – 1}{X(t – i\mathrm{\Delta}t)} $$

Not averaged are external settings that can change instantaneously, e.g. external pressure