# Langevin thermostat with fixcm=True: Issue #1285 fix introduced new bug in ASE 3.26.0 ## Summary Issue #1285 reported that the Langevin thermostat with `fixcm=True` produced incorrect temperature distributions in ASE 3.24. ASE 3.26.0 attempted to fix this, but the fix introduced a new bug: while the mean temperature is now correct, **the distribution is too wide with extreme unphysical outliers**. ## Background **Previous bug (ASE 3.24):** https://gitlab.com/ase/ase/-/issues/1285 - Langevin with `fixcm=True` **reported** incorrect temperatures (mean \~150K for 300K target) - **However, the actual distribution of states was correct** - bond distance distributions matched the target temperature - This was a **temperature calculation/reporting bug**, not a dynamics bug **Current bug (ASE 3.26.0):** - Temperature reporting is now correct ✓ - **But the actual distribution of states is now wrong** ✗ - The dynamics now sample an ensemble that is too hot - Bond distance distribution is 47% broader than expected ## Expected Behavior ### Temperature Distribution For a canonical ensemble (NVT), the instantaneous temperature should follow a **Gamma distribution**: ``` P(T) ∝ T^(f/2-1) * exp(-f*T/(2*T_target)) ``` Where f = degrees of freedom: - **fixcm=False**: f=6 (3N for N=2 atoms) → σ\_T = 173 K - **fixcm=True**: f=3 (1 vib + 2 rot, COM removed) → σ\_T = 245 K ### Bond Distance Distribution For a diatomic molecule in thermal equilibrium, the O-O distance should follow a **Boltzmann distribution**: ``` P(r) ∝ exp(-V(r)/(k_B*T)) ``` Where V(r) is the Morse potential. At T=300K with the Morse parameters used (ε=5.2 eV, r₀=1.21 Å, ρ₀=2.7): - **Expected mean:** \~1.212 Å - **Expected std dev:** \~0.022 Å The bond distance distribution should be **independent of the fixcm setting** - removing center-of-mass motion doesn't change the internal vibrational dynamics. Therefore: - **fixcm=True and fixcm=False should give identical O-O distance distributions** ## Actual Behavior ### ASE 3.24 with fixcm=True (original bug from #1285) - **Reported** mean temperature: **149.7 ± 136.2 K** (target: 300 K) - **BUT: O-O distance distribution matches Boltzmann theory at 300K!** - Measured: Mean = 1.2124 Å, σ = 0.0213 Å - Theory (Boltzmann at 300K): Mean ≈ 1.212 Å, σ ≈ 0.022 Å ✓ - Matches fixcm=OFF: 1.2124 Å, σ = 0.0221 Å ✓ - This proves the **actual dynamics were at \~300K** - **Conclusion: Temperature calculation bug**, not a thermostat bug ### ASE 3.26.0 with fixcm=True (new bug after fix) - **Reported** mean temperature: **309.3 ± 268.1 K** (target: 300 K) ✓ Reporting fixed! - **BUT: O-O distance distribution NO LONGER matches Boltzmann theory** - Measured: Mean = 1.2150 Å, **σ = 0.0313 Å** - Theory (Boltzmann at 300K): Mean ≈ 1.212 Å, σ ≈ 0.022 Å - fixcm=OFF: Mean = 1.2123 Å, σ = 0.0218 Å ✓ - **σ is 47% broader than theory** - system is sampling a hotter ensemble - **σ is 44% broader than fixcm=OFF** - the two should be identical! - **Conclusion: The fix corrected temperature reporting but broke the actual dynamics** ### ASE 3.24/3.26.0 with fixcm=False (both work correctly) - **ASE 3.24**: 278.6 ± 164.7 K (σ matches theory: 173 K ✓, mean \~7% low) - **ASE 3.26.0**: 306.4 ± 169.6 K (σ matches theory: 173 K ✓, mean perfect ✓) ## Impact **ASE 3.24:** Temperature reporting bug - simulations were actually correct but reported wrong temperatures. **ASE 3.26.0:** The fix for #1285 introduced a **worse bug** - the reported temperatures are now correct, but the **actual dynamics are wrong**. Simulations now suffer from: 1. Incorrect ensemble sampling - system is too hot 2. Bond distance distributions 47% broader than expected **This is more serious than the original bug** because users will see "correct" reported temperatures (\~300K mean) and not realize their simulations are sampling the wrong ensemble. At least in 3.24, the wrong reported temperature (\~150K) was an obvious red flag. ## Minimal Reproducible Example ```python import numpy as np from ase import Atoms, units from ase.md.langevin import Langevin from ase.calculators.morse import MorsePotential from ase.md.velocitydistribution import MaxwellBoltzmannDistribution # Create O2 molecule atoms = Atoms('O2', positions=[[0, 0, 0], [0, 0, 1.21]]) atoms.center(vacuum=10.0) atoms.set_pbc([True, True, True]) # Morse potential calc = MorsePotential(epsilon=5.2, r0=1.21, rho0=2.7) atoms.calc = calc # Initialize velocities MaxwellBoltzmannDistribution(atoms, temperature_K=300) # Langevin dynamics with DEFAULT fixcm=True dyn = Langevin( atoms, timestep=0.5 * units.fs, temperature_K=300, friction=0.01, fixcm=True # DEFAULT ) # Collect temperature data temperatures = [] def collect(): temperatures.append(atoms.get_temperature()) dyn.attach(collect, interval=1) dyn.run(200000) # 100 ps print(f"Mean: {np.mean(temperatures):.1f} K (target: 300 K)") # Check O-O distance distribution distances = [np.linalg.norm(atoms.get_positions()[1] - atoms.get_positions()[0]) for _ in range(len(temperatures))] print(f"O-O distance σ: {np.std(distances):.4f} Å (expected: ~0.022 Å)") ``` **Expected output (ASE 3.26.0):** ``` Mean: ~300 K ✓ O-O distance σ: ~0.022 Å (same as fixcm=OFF) ``` **Actual output (ASE 3.26.0):** ``` Mean: 309.3 K ✓ O-O distance σ: 0.0313 Å ✗ (47% broader than fixcm=OFF) ``` ## Comparison Plot {width=900 height=334} **Left panel:** Temperature distributions compared to theoretical Gamma distributions - Solid black: Theory for f=6 (fixcm=OFF), σ=173 K - Dashed black: Theory for f=3 (fixcm=ON), σ=245 K - Cyan/Orange (fixcm=OFF): Match solid theory line perfectly ✓ - Blue (ASE 3.24 fixcm=ON): Wrong mean, wrong shape (original bug #1285) - Red (ASE 3.26.0 fixcm=ON): Roughly follows dashed theory line **Right panel:** O-O distance distributions compared to Boltzmann theory at 300K - Black curve: Theoretical Boltzmann distribution P(r) ∝ exp(-V(r)/(k_B×300K)) - Blue/Cyan/Orange: Match theory perfectly ✓ - Red (ASE 3.26.0 fixcm=ON): **Significantly broader than theory** ✗ - smoking gun evidence that the dynamics are wrong ## System Information - **ASE versions tested:** 3.24.0, 3.26.0 - **Python version:** 3.13 - **Platform:** macOS Darwin 24.6.0 - **Test system:** O2 diatomic molecule with Morse potential - **Simulation parameters:** - Duration: 100 ps (200,000 steps @ 0.5 fs timestep) - Target temperature: 300 K - Friction: 0.01 fs⁻¹ ## Statistical Summary ### Temperature Statistics (Reported) | Condition | Mean (K) | Std Dev (K) | Expected σ | Status | |-----------|----------|-------------|------------|--------| | ASE 3.24 fixcm=ON | 149.7 | 136.2 | 245 | :x: Reporting wrong (but dynamics OK) | | ASE 3.26.0 fixcm=ON | 309.3 | 268.1 | 245 | ✓ Reporting fixed | | ASE 3.24 fixcm=OFF | 278.6 | 164.7 | 173 | ✓ σ matches, mean \~7% low | | ASE 3.26.0 fixcm=OFF | 306.4 | 169.6 | 173 | :white_check_mark: Perfect! | ### O-O Distance Statistics (Actual Dynamics) | Condition | Mean (Å) | Std Dev (Å) | vs Theory | vs fixcm=OFF | Status | |-----------|----------|-------------|-----------|--------------|--------| | **Theory (Boltzmann, 300K)** | \~1.212 | \~0.022 | \- | \- | \- | | ASE 3.24 fixcm=ON | 1.2124 | 0.0213 | ✓ Match | Same | :white_check_mark: Dynamics correct! | | ASE 3.26.0 fixcm=ON | 1.2150 | **0.0313** | ✗ **+42%** | ✗ **+44%** | :x: Too broad - system too hot | | ASE 3.24 fixcm=OFF | 1.2124 | 0.0221 | ✓ Match | \- | :white_check_mark: | | ASE 3.26.0 fixcm=OFF | 1.2123 | 0.0218 | ✓ Match | \- | :white_check_mark: | ## Workaround / Suggested Immediate Fix Setting `fixcm=False` works correctly in both ASE 3.24 and 3.26.0. **Suggested immediate fix:** Change the default value of `fixcm` from `True` to `False` in the Langevin thermostat until the underlying bug is resolved. While this contradicts the original intent for isolated molecules/clusters, it prevents users from unknowingly sampling incorrect ensembles. ## Suggested Investigation The fix for issue #1285 may have overcorrected the energy balance. The interaction between `fixcm` and the Langevin stochastic forces appears to inject too much energy into the system, causing the bond distance distribution to broaden beyond the theoretical Boltzmann prediction. ## Full Reproduction Code Complete test scripts, simulation data (100 ps trajectories), and analysis code available if needed. --- **Related issue:** #1285 (original fixcm=True temperature bug in ASE 3.24)