Release v2.13.0

Changed

• Deformation and shear stages now say what a modulus from them is worth. Every run
  reports the strain it will reach and is told that one ramp gives a modulus but not
  its uncertainty: eight independent ramps of one protocol on one network gave shear
  moduli from 1.29 to 1.74 GPa --- 10% scatter --- while each run's own fit reported an
  error an order of magnitude smaller. Take the scatter between seeds as the error bar,
  not the fit's.

  A run fast enough that little of it sits at low strain also draws a warning, since what
  gets fitted is then a secant over a large window. It does not claim a bias, because
  none is established: a single slow ramp appeared to show the default understating G
  threefold and a replicate of that same protocol refuted it. Rate dependence past yield
  is real and does replicate --- 1.22 +/- 0.06 GPa at edot = 1e-3 against
  0.96 +/- 0.02 at 1e-4. The default is unchanged.

• Deformation and shear stages now sample stress ten times more often. Both
  inherited nstenergy = 500 from the packaged npt.mdp, one frame per ps --- and
  stress is the quantity these stages exist to measure. Tested on a real shear run:
  going to a frame every 0.1 ps cut the standard error on the fitted modulus by
  sqrt(10), which is the independent-sample scaling, so the extra frames carry real
  information rather than repeating themselves. The cost is a larger edr and nothing
  else --- no change to the trajectory, the ensemble, or the physics.

Fixed

• Gromacs warnings were discarded, so nobody ever saw one. htpolynet runs every
  grompp with -maxwarn 4, so up to four warnings do not stop a run, and the command
  wrapper only echoes captured output when a command fails. Between them, a warning on
  a successful run reached neither the user nor the diagnostic log.

  That is how the Berendsen barostat deprecation --- emitted by every constant-pressure
  stage htpolynet runs, on every Gromacs since 2025 --- went unnoticed for as long as it
  has been emitted, and it is why the shear stage's off-diagonal conflict had to be found
  by running grompp by hand. Warnings are how Gromacs says a setup is questionable,
  and they are worth seeing even when the run proceeds. They are now logged at WARNING
  level, with the mdp that produced them.

• Uniaxial deformation was rejected by Gromacs 2025 before it started, so Young's
  modulus could not be measured at all. From Gromacs 2025, combining deform with
  generated velocities is a hard grompp error unless deform-init-flow = yes is set,
  because the initial velocities have to carry the flow profile the deformation implies.
  htpolynet never set it. Every deform stage therefore failed at grompp on a current
  Gromacs --- not with a bad number, but with no run. Both deform and the new shear
  now set it.

  This was invisible because nothing ran the stage: analysis/postsim.py had no tests,
  and the deformation traces on hand were produced by an older Gromacs. It was found by
  putting a real mdp through gmx grompp while trying to build a shear fixture.

Added

• The shear stage's box coupling is now a choice, because it changes the answer by
  29%. ring_cure-style shear has to be run somehow, and there are three defensible
  ways to hold the box while shearing it; measured on one network with replicates they
  give G of 1.17, 1.35 and 1.51 GPa. That spread is wider than the trajectory scatter
  of any of them, so it does not average away. coupling: clamped (the default, and
  what every measurement so far used) keeps the two unsheared directions at pressure;
  coupling: isochoric uses no barostat at all, needs no uncoupling, and is the only one
  of the three that passes grompp with no warnings whatsoever.

  The default is unchanged, and not because it is the best: the Poisson ratio read off
  the box dimensions ranks it last of the three. It is unchanged because that is one
  cross-check at one strain, and switching would silently move every G anyone has
  measured. The stage now logs which coupling it used, so a G is never quoted without
  it.

• A shear postsim stage, for the shear modulus. The counterpart of deform:
  Gromacs' deform drives one off-diagonal element of the box matrix at a constant
  rate, and G is the slope of shear stress against engineering shear strain, exactly
  as E is the slope of tensile stress against tensile strain. direction names the
  plane --- xy, xz or yz --- and edot is the shear rate in ps^-1.

  The normal directions stay pressure-coupled at P while the off-diagonal is driven,
  because the anisotropic barostat is given zero compressibility off the diagonal and so
  does not oppose the deformation. The stage writes the driven box element and the
  matching off-diagonal of the pressure tensor, with derived -strain and -stress
  columns; analysis.utils.compute_G fits those, as compute_E does for tension.

  One thing that differs from the tensile case and is easy to get wrong: the engineering
  shear strain is the driven box element over the length of the sheared face and starts
  at zero, so it carries no -1.

  A known disagreement, stated because it is not resolved in this release. Shearing
  by deform requires the barostat to be switched off for the Cartesian component the
  shear moves, since Gromacs refuses the run otherwise. G from that configuration
  comes out 1.51 +/- 0.05 GPa where a less constrained variant gives 1.17 +/- 0.02 on the
  same network --- and only the latter is consistent with that network's own E and
  Poisson ratio. So G here disagrees with E by about a quarter.

  The tempting explanation, that the constraint inflates G, does not survive its own
  test: a constant-volume variant fixes all three box dimensions rather than one, so it
  is more constrained still, and its mean sits below the clamped value. The cause is
  open. Treat G from this stage as disagreeing with E by roughly a quarter, and do
  not assume you know why.

  The mdp is accepted by Gromacs; no G has been measured. Putting it through
  gmx grompp found two defects the unit tests could not: the barostat must be
  uncoupled for the Cartesian component the shear moves, not merely off the diagonal ---
  grompp warns of "spurious periodicity effects" otherwise --- and deform-init-flow
  must be set, which is a hard error without. Both are fixed and grompp now accepts
  every plane. The energy terms it traces (Box-YX, Box-ZX, Box-ZY, Pres-XY,
  Pres-XZ, Pres-ZY) were confirmed present in Gromacs 2025.4.

  First real run, on a threading-filtered BADCy network of 12,600 atoms: both stages
  completed, the stress-strain curve has the right shape for a glassy thermoset --- an
  elastic rise to about gamma 0.08, a roll-over, then a yield plateau near 160-175 MPa
  --- and G comes out 1.3 +/- 0.2 GPa over gamma 0.08-0.10. deform with
  deform-init-flow was also verified on Gromacs 2026.3, not only the 2025.4 used here.

  That G is one window's slope, not a zero-strain tangent, and should not be quoted as
  a converged modulus. Below gamma 0.03 the fit has no power at all --- the standard
  error exceeds the estimate and the sign flips --- because per-frame virial noise is
  456 bar and a 300 ps ramp at edot = 1e-3 puts only 20 frames below that strain.
  Above it the fitted G falls from 2.08 to 1.31 GPa as the window grows, the same way
  the uniaxial fixture's E falls from 4.74 to 2.72. A slower ramp, or more frequent
  energy output at the same rate, is what a modulus fixture needs. ROADMAP.md carries
  this, and the related finding that the uniaxial trace has no resolved elastic region
  either.

• Real simulation traces as test fixtures, for Tg and Young's modulus. Until now
  every test of those fits used synthetic data, where a clean straight line recovers its
  slope to eight digits. tests/unit/fixtures/ now carries a real 61-point temperature
  ladder and a real 1001-frame uniaxial deformation, the latter written by
  PostSimDeform itself.

  They show what synthetic data cannot. On the real ladder, Tg spans 377 to 406 K
  across fit windows a reasonable person would choose; on the real deformation, E spans
  2420 to 4398 MPa, nearly a factor of two. The R-squared of a perfectly sensible
  modulus fit is under 0.5, because single-frame pressure in MD scatters hugely --- so
  anyone tempted to add an R-squared quality gate should look at that first. All of
  this is now pinned by tests, and ROADMAP.md carries what to do about it.

  There is deliberately no shear fixture, and a test asserts its absence so that
  adding one forces the accompanying tests to be written. No shear simulation has been
  run.

• Tests for the two numbers postsim exists to produce. analysis/postsim.py was a
  zero-coverage module, and compute_E and compute_tg --- the fits behind Young's
  modulus and the glass transition --- had no tests at all, while the shear stage added
  hours earlier did. Now covered: both fits against known slopes and a known
  intersection, the mdp construction for every deform and shear direction, that the
  driven axis is not pressure-coupled while the others are, and that an unknown direction
  is refused rather than guessed.

  Two of those tests document behavior rather than assert it is good. compute_E fitted
  across a rolled-over stress curve understates the modulus by more than a factor of two,
  which is what the fit window exists to prevent. And compute_tg trusts its windows: a
  hot window reaching about 20 K below the transition biases Tg roughly 7 K low, with no
  check and no warning. ROADMAP.md carries what a fix would look like.