Conserved quantities and ensemble measure for Martyna--Tobias--Klein barostats with restricted cell degrees of freedom

This paper derives the conserved energy-like quantity and ensemble measure for Martyna--Tobias--Klein barostats with restricted cell degrees of freedom, demonstrating that the standard formulation adapts by replacing the total degrees of freedom with the number of active axes to ensure exact conservation and correct sampling of the restricted isothermal--isobaric ensemble.

The Gist

Imagine you are a chef trying to bake the perfect loaf of bread in a magical, self-adjusting oven. In the world of computer simulations (specifically molecular dynamics), this "bread" is a collection of atoms, and the "oven" is the environment around them.

To get the bread right, you need to control two things:

1. Temperature: You need a thermostat to keep it from getting too hot or cold.
2. Pressure: You need a "barostat" (a pressure controller) to make sure the bread expands or contracts to the right size.

For a long time, scientists had a very sophisticated recipe called the Martyna–Tobias–Klein (MTK) method. This recipe allowed the oven to change shape in every possible direction at once. Imagine the oven is a stretchy, transparent box. The MTK method let that box stretch, shrink, twist, and turn in all 3 dimensions simultaneously. It was powerful, but sometimes, you don't need that much flexibility.

The Problem: The "Over-Engineered" Oven

Sometimes, you only want to control the pressure in specific directions.

• Example 1: You are simulating a thin sheet of material (like a slice of bread). You only care if it gets thicker or thinner (up and down). You don't want it to stretch sideways or twist.
• Example 2: You are simulating a crystal that can only expand along one specific line.

In the old MTK recipe, even if you only wanted to control one direction, the computer still tried to calculate how the box would twist and turn in all other directions. It was like using a giant, complex robotic arm to just turn a light switch on and off. It worked, but it was messy and hard to track the "energy" of the system.

The Solution: The "Masked" Barostat

This paper, written by Kohei Shinohara, introduces a clever trick: The Masked MTK Barostat.

Think of the "mask" as a pair of glasses with opaque lenses.

• The computer still sees the whole stretchy box.
• But the "mask" tells the computer: "Hey, only look at these specific axes (directions). Ignore everything else. Pretend the other sides of the box are made of solid steel and cannot move."

The author figured out the mathematical rules for this "masked" version. The big challenge was proving that if you use this mask, the simulation doesn't go crazy. In physics simulations, if you mess up the math, the "energy" of the system might drift, and your simulation will produce garbage results (like your bread turning into a rock).

The Three Big Discoveries

1. The "Conserved Energy" (The Unbreakable Ledger)
In physics, there's a rule that total energy must stay constant (or follow a specific pattern) for the simulation to be real.

• The Old Way: The formula for this energy was complex, counting every single way the box could wiggle (in 3D, that's 9 different ways).
• The New Way: Shinohara showed that when you put the mask on, you just need to swap the number 9 for the number of active directions (let's say 1).
• The Analogy: Imagine a bank account. The old formula said, "You have 9 different types of savings accounts." The new formula says, "You only have 1 active savings account." The math is almost identical; you just change the count. This proves the simulation is still mathematically sound.

2. The "Rule of the Room" (The Ensemble Measure)
When you simulate atoms, you aren't just watching one specific moment; you are trying to sample all possible states the system could be in (like taking millions of photos of the bread rising).

• The paper proves that even with the mask, the computer is still taking the right "photos." It correctly samples the Isothermal-Isobaric Ensemble (constant temperature and pressure).
• The Analogy: Imagine you are trying to guess the average height of people in a room. If you only look at the people standing in the center (the active axes) and ignore the people against the walls (the fixed axes), you might get the wrong average. Shinohara proved that if you follow his new math, the "average" you calculate is still perfectly correct, even though you are ignoring the walls.

3. The "Step-by-Step Dance" (The Integration Scheme)
To run the simulation, the computer has to take tiny steps forward in time.

• The author provided a new "dance routine" (an algorithm) for the computer to follow. This routine is built on the same logic as the old one but skips the steps where the "inactive" parts of the box would move.
• The Analogy: It's like a dance where the old version required you to spin, jump, and slide. The new version says, "Okay, just slide left and right. Don't spin." The author wrote down the exact steps so the computer doesn't trip over its own feet.

Why Does This Matter?

This might sound like abstract math, but it has real-world uses:

• Materials Science: Scientists can now simulate thin films (like the coating on your phone screen) much faster and more accurately because they don't have to waste computer power calculating sideways movements that don't happen.
• Efficiency: It saves computing time. By "masking" the directions you don't care about, the simulation runs faster.
• Reliability: Before this paper, people who tried to do this "masked" simulation were guessing. They didn't have a mathematical proof that their energy was conserved. Now, they have a verified recipe.

Summary

Kohei Shinohara took a complex, all-powerful tool for simulating pressure in materials and figured out how to turn off the parts you don't need without breaking the machine. He proved that if you "mask" the directions you don't want to move, the physics still works perfectly, the energy stays balanced, and the results are scientifically valid. It's like taking a Swiss Army knife and realizing you can safely lock away the saw and screwdriver to just use the blade, without the whole tool falling apart.

Technical Summary

1. Problem Statement

In molecular dynamics (MD) simulations, the Martyna–Tobias–Klein (MTK) equations of motion are the standard method for sampling the isothermal–isobaric (NPT) ensemble by coupling the simulation cell to a barostat. While the standard MTK formulation handles both isotropic and fully anisotropic cell fluctuations (where all $d^2$ components of the cell matrix fluctuate in $d$ dimensions), many practical applications require restricted cell fluctuations.

• Examples: Slab geometries often require pressure control only along the surface normal ($NP_zT$), and uniaxial stress simulations constrain lateral dimensions.
• The Gap: Existing literature provides derivations for isotropic and fully anisotropic cases but lacks a compact, rigorous derivation for the intermediate case where only a subset ($n_c$) of cell axes are active. Furthermore, there is no established derivation for the corresponding conserved energy-like quantity and ensemble measure for these "masked" barostats, making it difficult to verify energy conservation and ensemble correctness in such simulations.

2. Methodology

The author employs the generalized Liouville framework for non-Hamiltonian systems to derive the theoretical foundations for the masked MTK barostat. The methodology involves:

• Defining the System: Restricting the barostat momentum $p_g$ to a subspace spanned by $n_c$ active axes. The cell matrix $h$ is allowed to have arbitrary shape, but the active axes must be orthogonal to all other axes (a condition automatically satisfied by orthorhombic cells).
• Deriving Equations of Motion (EOM): Formulating the EOMs for particle positions, momenta, and the restricted cell degrees of freedom, coupled with Nose-Hoover chain (NHC) thermostats and barostats.
• Deriving Conserved Quantities: Constructing a conserved energy-like quantity ($H'$) and verifying its time derivative is zero ($\dot{H}' = 0$) using the derived EOMs.
• Calculating Phase-Space Compressibility: Computing the phase-space compressibility ($\kappa$) to determine the metric factor ($e^{-w}$) required for the ensemble measure.
• Ensemble Verification: Integrating the microcanonical partition function to prove that the dynamics samples the correct restricted isothermal–isobaric ensemble.
• Integration Scheme: Adapting the Liouville-operator-based symmetric Trotter splitting (used in standard MTK) to the masked variant.

3. Key Contributions

The paper provides the first rigorous derivation for MTK barostats with restricted degrees of freedom. The specific contributions are:

A. Modified Equations of Motion

The EOMs are derived for a system where only $n_c$ diagonal components of the cell matrix fluctuate.

• The barostat momentum $p_g$ is restricted to $n_c$ scalar components ($p_c$).
• The driving forces for the barostat NHC chain are modified to depend on $n_c$ rather than the full dimensionality.

B. The Conserved Energy-Like Quantity ($H'$)

A new conserved quantity is derived (Eq. 17):
$$H' = H(r, p) + \sum_{c=1}^{n_c} \frac{p_c^2}{2W_g} + P \det[h] + \sum_{j=1}^M \left( \frac{p_{\eta_j}^2}{2Q_j} + \frac{p_{\xi_j}^2}{2Q'_j} \right) + N_f kT \eta_1 + n_c kT \xi_1 + kT(\eta_c + \xi_c)$$
Key Structural Change: Compared to the fully anisotropic case, the term $d^2$ (representing total barostat DOFs) is replaced by $n_c$ in:

1. The sum of barostat kinetic energy.
2. The coupling term to the first barostat chain variable ($n_c kT \xi_1$).
3. The mass parameter $Q'_1$ for the first barostat chain ($Q'_1 = n_c kT \tau^2$).

C. Ensemble Measure and Metric Factor

The phase-space compressibility is calculated as:
$$\kappa = -\frac{d}{dt}(N_f \eta_1 + n_c \xi_1 + \eta_c + \xi_c)$$
Consequently, the metric factor is $e^{-w} = \exp(N_f \eta_1 + n_c \xi_1 + \eta_c + \xi_c)$.
The paper proves that integrating out the momenta and auxiliary variables yields a partition function proportional to the standard NPT ensemble, restricted to the submanifold where the $d - n_c$ inactive cell axes are held fixed.

D. Integration Scheme

A complete, time-reversible, measure-preserving integration scheme is provided (Appendix D). It utilizes the same Liouville operator splitting as the standard MTK but modifies the barostat operator to sum only over the active axes. This ensures the numerical stability and symplectic-like properties of the simulation.

4. Results and Verification

• Exact Conservation: The author analytically verifies that $\dot{H}' = 0$, confirming that the derived energy-like quantity is exactly conserved under the new equations of motion.
• Ensemble Correctness: The derivation demonstrates that the microcanonical partition function of the extended system maps directly to the isothermal–isobaric partition function $\Delta(N, P, T)$, provided the inactive axes are fixed.
• Momentum Conservation: In the absence of external forces, a conserved momentum-like quantity $K(x)$ is identified, ensuring the dynamics remains confined to the correct submanifold even with total momentum constraints.
• Implementation: The author notes that this formulation has been implemented in the ASE (Atomic Simulation Environment) library, making it immediately available for practical use.

5. Significance and Limitations

Significance:

• Theoretical Closure: It fills a critical gap in the theoretical literature, providing the necessary tools to verify energy conservation and ensemble correctness for restricted-axis NPT simulations.
• Practical Utility: Many materials science simulations (e.g., surfaces, nanowires, 2D materials) rely on restricted barostats. This work ensures these simulations are statistically rigorous.
• Efficiency: By reducing the number of barostat degrees of freedom from $d^2$ to $n_c$, the method can be computationally more efficient and numerically stable for specific geometries.

Limitations:

• Orthogonality Constraint: The derivation strictly requires that active axes be orthogonal to all other axes.
  • Orthorhombic cells: Fully supported.
  • Hexagonal cells: Only the unique axis (c-axis) can be barostatted alone; simultaneous control of in-plane vectors requires redefining the cell as an orthohexagonal supercell.
  • Triclinic cells: Fully excluded, as no pair of axes is generally orthogonal.
• Diagonal Fluctuations Only: The current formulation addresses only diagonal cell-length fluctuations. Off-diagonal tilt factors are left for future work.

In summary, this paper provides the essential mathematical framework and integration algorithms for performing rigorous, constant-pressure molecular dynamics simulations where only specific dimensions of the simulation cell are allowed to fluctuate.