Local temperature measurement in molecular dynamics simulations with rigid constraints

This paper presents a method for accurately calculating local temperatures in molecular dynamics simulations with rigid constraints by self-consistently evaluating the degrees of freedom, thereby correcting unphysical violations of kinetic energy equipartition and providing a sensitive indicator for numerical integration errors or insufficient equilibration.

The Gist

The Big Picture: The "Temperature" of a Dancing Crowd

Imagine you are at a crowded dance party. You want to know how hot the room is. In a normal room, you might just ask a few people, "How hot do you feel?" and take an average. In the world of computer simulations (Molecular Dynamics), scientists do the same thing: they look at how fast atoms are moving to calculate the "temperature."

But here's the catch: Atoms in these simulations are often glued together.

Think of a water molecule not as three separate dancers, but as a single rigid unit where the atoms are holding hands so tightly they can't let go. This is called a rigid constraint. It makes the computer run faster because it doesn't have to calculate every tiny wiggle of the bond.

The Problem:
When scientists try to measure the temperature of just part of this glued-together group (like just the Hydrogen atoms, or just the Carbon atoms), they run into a math problem.

If you have a rigid dumbbell (two balls connected by a stick) and you spin it, the whole thing moves as one. But if you only look at one ball, how much of that "spin energy" belongs to it?

• The Old Way: Scientists used to just split the energy evenly. "Okay, there are two atoms, so each gets half the energy."
• The Reality: This is like saying a heavy bouncer and a light dancer in a rigid group contribute equally to the group's spin. It's wrong! The heavy one contributes more to the "weight" of the spin, and the light one contributes less.

Because of this error, the computer thinks the light atoms are freezing cold and the heavy atoms are boiling hot, even though the whole system is supposed to be at a perfect, comfortable room temperature. This creates "ghost temperatures" that don't exist in real life.

The Solution: The "Inertia" Scale

The authors of this paper, Stephen Sanderson and his team, came up with a new, smarter way to slice the pie. Instead of splitting the energy evenly, they split it based on Inertia (how much an object resists moving).

The Analogy: The Merry-Go-Round
Imagine a merry-go-round (the molecule) spinning.

• The Heavy Person (Carbon): Sits near the edge. They have a lot of mass. When the ride spins, they are doing a lot of the "work" of the rotation. They "own" more of the energy.
• The Light Person (Hydrogen): Sits near the center or is just very light. They are moving, but they contribute less to the overall spin energy.

The new method says: "Don't just count heads; weigh the contribution."
They calculate exactly how much each atom contributes to the total "spin" and "slide" of the molecule. Then, they assign the temperature based on that specific contribution.

Why Does This Matter? (The "Overheating" Detector)

The paper shows that this new method isn't just about being mathematically correct; it's a canary in the coal mine for bad computer simulations.

The Metaphor: The Blurry Photo
When you take a photo with a camera, if you move the shutter too slowly, the picture gets blurry. In simulations, the "shutter speed" is the time step (how often the computer checks the position of atoms).

• If the time step is too slow (too big), the computer misses tiny details.
• Usually, the total temperature of the system looks fine.
• But, with the new method, if you look at the specific atoms (Carbon vs. Hydrogen), you see they have different temperatures.

The Discovery:
The authors found that even at a time step commonly used by scientists (2 femtoseconds, which is incredibly fast but still "slow" for this specific math), the Carbon atoms were "hotter" than the Hydrogen atoms.
This didn't mean the room was actually hot; it meant the computer simulation was "blurring" the details of how the atoms were arranged. It was a sign that the simulation was slightly inaccurate, even if the total energy looked okay.

Summary of the Takeaways

1. The Issue: When atoms are glued together in simulations, splitting the temperature evenly between them is wrong. It makes some atoms look artificially hot and others cold.
2. The Fix: The authors created a formula that divides the energy based on how much each atom actually "weighs" in the movement (inertia). This ensures every part of the molecule reports the same, correct temperature.
3. The Bonus: This new math acts like a stress test. If you use a time step that is too slow for your computer, the different atoms will start reporting different temperatures. This tells the scientist, "Hey, your simulation is getting blurry; slow down or fix your settings!"

In a nutshell: The paper teaches us how to stop a computer simulation from lying about how hot its atoms are, and gives us a new tool to catch when the simulation is running too fast to be accurate.

Technical Summary

1. Problem Statement

In Molecular Dynamics (MD) simulations, rigid constraints (e.g., fixed bond lengths and angles via algorithms like SHAKE, RATTLE, or LINCS) are commonly used to increase time steps and computational efficiency. However, these constraints reduce the system's degrees of freedom (DoF).

• The Core Issue: When calculating local temperatures (e.g., within a specific subvolume, for a specific atom type, or along a specific Cartesian direction), standard methods often fail to correctly account for the reduced DoF of constrained atoms.
• Common Pitfalls: Popular MD codes often distribute DoF evenly among atoms (e.g., subtracting 0.5 DoF per atom for a fixed bond). While this may yield correct global averages, it leads to unphysical local temperature variations in inhomogeneous systems (e.g., interfaces, flows, or aligned fields).
• Consequences:
  • Apparent violations of the equipartition theorem between neighboring subsets (e.g., Carbon vs. Hydrogen atoms in a rigid molecule).
  • Incorrect thermostatting in non-equilibrium simulations.
  • Misinterpretation of local heating or thermal conductivity.
  • Existing solutions are often case-specific and difficult to generalize to complex, multi-constraint systems.

2. Methodology

The authors propose a general, self-consistent framework to calculate the correct DoF for any subset of atoms ($S$) within a constrained system, applicable to both rigid bodies and semi-rigid fragments.

A. Theoretical Framework

The method relies on the principle that at equilibrium, the kinetic temperature calculated for any non-empty subset of a rigid body must be identical.

1. Rigid Bodies:

   • The kinetic energy of a rigid body is decomposed into translational and rotational modes.
   • Translational DoF: Partitioned based on the atom's mass contribution to the total mass of the rigid body ($d_{tr} \propto m_j / M_{total}$).
   • Rotational DoF: Partitioned based on the atom's contribution to the modal inertia. By diagonalizing the inertia tensor, the DoF for a specific rotational mode $\beta$ is assigned to atom $j$ proportional to its contribution to the moment of inertia for that mode ($d_{rot} \propto I_{j,\beta} / I_{total,\beta}$).
   • Directional Decomposition: The framework extends to calculate DoF for specific Cartesian directions (e.g., $x, y, z$), allowing for the measurement of "directional temperatures" (crucial for interfaces or shear flows).

2. Semi-Rigid Fragments:

   • For systems with some flexible degrees of freedom (e.g., fixed bonds but flexible angles), the authors utilize an internal coordinate system.
   • They construct a Jacobian matrix ($J$) that maps internal coordinate velocities to laboratory frame velocities.
   • The DoF is calculated by projecting the generalized inertia onto the modal frame, ensuring the equipartition theorem holds for the internal modes.

3. Volumetric Bodies:

   • The theory is generalized to continuous mass distributions (volumetric rigid bodies) using mass density functions, allowing for temperature gradients across a single rigid object.

B. Simulation Validation

The theory was validated using LAMMPS on several test systems:

• Ethane/Graphene Interface: A system of rigid ethane molecules in contact with a flexible graphene membrane.
• Confined Water: Rigid SPC/E water molecules confined between copper walls.
• Rigid Dumbbells: A system with a temperature gradient across a rigid dumbbell to test heat transfer.
• Bulk Ethane: Semi-rigid ethane (fixed C-H bonds, flexible angles) to test time-step dependence.

3. Key Contributions

1. Generalized DoF Partitioning: A mathematically rigorous method to assign fractional DoF to individual atoms based on their mass and geometric contribution to the system's inertia, rather than arbitrary averaging.
2. Directional Temperature Calculation: A novel approach to calculate kinetic temperatures along specific Cartesian axes, essential for non-equilibrium simulations where streaming velocities or anisotropic ordering exist.
3. Semi-Rigid Extension: An extension of the rigid body formalism to semi-rigid fragments using Jacobian matrices, making the method applicable to standard organic force fields.
4. Discretization Error Detection: The identification of local temperature divergence as a sensitive metric for detecting numerical integration errors (discretization artifacts) that do not affect global energy conservation.

4. Results

• Interface Validation (Ethane/Graphene):
  • Naive/Uniform Partitioning: Produced artificial temperature "ringing" at the interface. Carbon atoms appeared too hot, and Hydrogen atoms too cold, due to incorrect DoF assignment.
  • Inertia-Based Partitioning: Yielded a uniform temperature profile of 300 K across the entire system, correctly satisfying equipartition.
• Directional Temperature (Confined Water):
  • Near the wall, water molecules are ordered, causing kinetic energy to be unequally distributed among Cartesian directions.
  • Using the new directional DoF calculation, the local temperature in every direction ($x, y, z$) converged to the set-point (300 K), whereas naive methods showed significant deviations.
• Temperature Gradients (Rigid Dumbbells):
  • Demonstrated that a rigid body can sustain a temperature gradient across its length when one end is in a hot bath and the other in a cold bath. The method correctly captured the translational and rotational mode correlations required to maintain this gradient.
• Time-Step Sensitivity (Bulk Ethane):
  • At small time steps (0.5 fs), C and H atoms in semi-rigid ethane converged to the correct 300 K.
  • At larger, commonly used time steps (2.0 fs), a steady-state violation of equipartition emerged: C atoms became significantly hotter (~4.5 K above set-point) while H atoms were cooler.
  • This divergence was linked to configurational overheating (broadening of bond angle/dihedral distributions) caused by discretization errors in the integrator, which are masked when only global kinetic energy is monitored.

5. Significance

• Accuracy in Non-Equilibrium MD: The method provides the necessary tools to correctly thermostat and analyze systems with interfaces, flows, and external fields, preventing artifacts in thermal conductivity and local heating studies.
• Force Field Development: The discovery that rigid constraints can lead to configurational overheating at standard time steps (2 fs) suggests that many existing force fields may be implicitly parameterized to compensate for these numerical errors.
• Computational Efficiency: The proposed method allows researchers to detect integration errors using only structure and momentum data (local temperatures) without needing expensive potential energy analysis or comparative simulations with smaller time steps.
• Universality: Unlike previous approaches, this framework is applicable to arbitrary geometric constraints, volumetric bodies, and semi-rigid fragments, making it a robust standard for future MD software implementations.