`NPT` bug

[jla-gardner]

I think there is a bug in the npy_langevin integrator.

To show this, I use a graph-pes model with no interatomic interactions, i.e. a PerfectGasModel. For a given pressure, number of particles, and temperature, we know what the volume of the cell should converge to: $V = NkT / p$.

I set up the experiment as follows:

import torch
import numpy as np
import ase
from graph_pes.models import FixedOffset
import torch_sim

# starting conditions
p = 0.001  # in eV/A^3
T = 10  # in K
N = 50

k = 8.313e-5  # in eV/K
expected_V = N * k * T / p  # in A^3

# start with an expanded cell
starting_V = expected_V * 1.5
l = starting_V ** (1 / 3)
structure = ase.Atoms(
    symbols=["Pb"] * N,
    positions=np.random.rand(N, 3) * l,
    cell=np.eye(3) * l,
)

perfect_gas = torch_sim.models.GraphPESWrapper(FixedOffset())
# prove that this model generates 0 stress contributions
stress_pred = perfect_gas(
    torch_sim.state.initialize_state(
        [structure], torch.device("cpu"), torch.float64
    )
)["stress"]
assert (stress_pred == 0).all()

torch_sim.integrate(
    system=structure,
    model=perfect_gas,
    n_steps=20_000,
    timestep=0.001,
    temperature=T,
    integrator=torch_sim.npt_langevin,
    external_pressure=-1.0 * p, # apply compression
    trajectory_reporter=torch_sim.TrajectoryReporter(
        "perfect-gas.h5md", state_frequency=50
    ),
)

I expect this to converge to the expected volume given by a perfect gas (~40Å$^3$) in this case.

However, this is what I observe instead:

import matplotlib.pyplot as plt
%config InlineBackend.figure_format = 'retina'

trajectory = torch_sim.TorchSimTrajectory("perfect-gas.h5md")
cells = trajectory.get_array("cell")[:150]
volumes = np.linalg.det(cells)


plt.plot(volumes)
plt.axhline(expected_V, color="k", linestyle="--", lw=1)
plt.text(len(cells), expected_V, f"$NkT/p$ = {expected_V:.2f}", ha="right", va="bottom")
plt.ylim(0, starting_V)
plt.xlabel("Time")
plt.ylabel("Volume (Å$^3$)")
plt.savefig("perfect-gas-bug.png", dpi=300, bbox_inches="tight")

[janosh]

fixed in #153