Fix NPT Nose-Hoover things

[alphalm4]

Summary

• md.py : NHC chain mass (Q) initialization: tau^4 -> tau^2 (although i believe this barely effects t>0 dynamics)
• npt.py : cell force KE -> 2KE.
  The MTK cell force equation (Martyna–Tobias–Klein, 1994, eq.(2.2) + eq.(2.9)) gives (ignoring V-dependent phi and thermostat part)

$$\begin{align} \dot{p}_\epsilon &= dV(P_{\text{int}} - P_{\text{ext}}) + \frac{d}{N_f} \sum_{i=1}^{N} \frac{{\mathbf{p}_i}^2}{m_i} \\ &= dV \left( \frac{1}{dV} \left[ \sum_{i=1}^{N} \frac{{\mathbf{p}_i}^2}{m_i} + \sum_{i=1}^{N} \mathbf{r}_i \cdot \mathbf{F}_i \right] - P_{\text{ext}} \right) + \frac{d}{N_f} \sum_{i=1}^{N} \frac{{\mathbf{p}_i}^2}{m_i} \\ &= \left(1 + \frac{d}{N_f}\right) \sum_{i=1}^{N} \frac{{\mathbf{p}_i}^2}{m_i} + \sum_{i=1}^{N} \mathbf{r}_i \cdot \mathbf{F}_i - dV P_{\text{ext}} \\ &= \left(1 + \frac{d}{N_f}\right) \times (2 \text{KE}) + dV P_{\text{int}}^{\text{static}} - dV P_{\text{ext}} \\ &= \alpha \times (2 \text{KE}) - V \, {\rm tr}(\sigma) - dV P_{\text{ext}} \end{align}$$

The missing factor of 2 halved the kinetic pressure driving cell volume changes, leading to incorrect pressure coupling. Verified against jax-md box_force which uses KE2 = sum(p²/m) directly.

• npt.py : Remove redundant model eval in NHC inner step
  NHC chain half-steps only rescale momenta -> we do not need to re-eval the model(state). Just calling state.stress is ok.
• npt.py and nvt.py
  I feel the DOF should be reduced from Nd to Nd - d since the COM movement is removed, but I'm not sure (it might be reduced at some other points)

I didn't added any test scripts, but I got some internal test results that this version of NPT-NH gives nearly the same density convergence with LAMMPS NPT-NH. It'll be great to be double checked the math things by code maintainers.

Checklist

Before a pull request can be merged, the following items must be checked:

• Doc strings have been added in the Google docstring format.
• Run ruff on your code.
• Tests have been added for any new functionality or bug fixes.

[thomasloux]

Good catch, that's actually my fault arghhh

[thomasloux]

I'm almost sure this is true

[thomasloux]

I think that it is something to take into account. Technically then we should consider all dofs reductions from constraints (although constraints are not really supported for NPT). Can you implement something like NPTCRescaleState.get_number_of_degrees_of_freedom ? It considers constraints and automatically applies reduction from Fix center of mass

[alphalm4]

I found that NPTNoseHooverState already has get_number_of_degrees_of_freedom methods. is it possible to use this?

[thomasloux]

Absolutely. The idea is that you should override this method to call super().get_number_of_degrees_of_freedom - 3

[alphalm4]

Fixed it. npt_nose_hoover_init gets SimState which has DOF=3N so i subtracted 3 but npt_nose_hoover_invariant gets NPTNoseHooverState which has DOF=3N-3 so i used as-is.

I tested npt_nose_hoover_invariant and observed a relative total Hamiltonian fluctuation of ~1e-4, which seems quite good. (Though I’m not sure what the practical threshold should be.)

[thomasloux]

Same than NPT case

[thomasloux]

Good catch. I don't validate yet because we use the exact integration scheme of:
A Liouville-operator derived measure-preserving integrator for molecular dynamics simulations in the isothermal–isobaric ensemble, 2006 J. Phys. A: Math. Gen. 39 5629

But I think that there is a small error in the equation in the paper.

The equation you're correcting is equation 5.8 of the paper. Indeed we are certainly lacking a factor 2 for kinetic energy. But interestingly the equation lacks a factor 3, corresponding to your d for the term - dVP_{ext}.

I would just modify the dofs to take into account constraints as I describe in my review.
Thanks again for the contribution!

[CompRhys]

I think to merge this the values in the npt tests need to be updated. They're set to confirm the algorithm is unchanged and reproducibles but the algo changed here.

[thomasloux]

let me have a final look on the removal of a model(state) evaluation to get state stress

[thomasloux]

Good for me after testing that the extra model evaluation was not necessary

[CompRhys]

Thanks alphalm4!