grompp crashes with deform-init-flow = yes for energy minimization, normal-mode analysis and test-particle insertion

Summary

With deform set and deform-init-flow = yes, gmx grompp crashes with a segmentation fault for every integrator without velocities (steep, cg, l-bfgs, nm, tpi, tpic).

deformInitFlow() adds the flow profile to state->v (https://gitlab.com/gromacs/gromacs/-/blob/dd255a67520b520ed49eb4bc5779cf4b7ca28891/src/gromacs/gmxpreprocess/grompp.cpp#L2075-2098), and it is called whenever deform-init-flow = yes (https://gitlab.com/gromacs/gromacs/-/blob/dd255a67520b520ed49eb4bc5779cf4b7ca28891/src/gromacs/gmxpreprocess/grompp.cpp#L2738-2742). But the state only has velocities for dynamical integrators (https://gitlab.com/gromacs/gromacs/-/blob/dd255a67520b520ed49eb4bc5779cf4b7ca28891/src/gromacs/gmxpreprocess/grompp.cpp#L803-816). For the other integrators state.v is empty, and deformInitFlow() writes past its end.

grompp even recommends the option: with deform and deform-init-flow = no, the same steep input gives the note "Unless the velocities in the initial configuration already obey the flow profile, the flow profile should be initialized by setting deform-init-flow=yes when using the deform option" (https://gitlab.com/gromacs/gromacs/-/blob/dd255a67520b520ed49eb4bc5779cf4b7ca28891/src/gromacs/gmxpreprocess/readir.cpp#L1794-1810). With gen-vel = yes in the mdp file (which these integrators otherwise ignore) and deform-init-flow = no, grompp stops with the error "When the box is deformed and velocities are generated, the flow profile should be initialized by setting deform-init-flow=yes", so there is no setting of deform-init-flow that works. As far as I can see, box deformation is only applied by the dynamical integrators, so for these integrators the option has no effect anyway.

GROMACS version

Tested: 2026.0 (conda-forge) and main at dd255a67 (2027.0-dev-20261009-dd255a67). By code inspection, all 2024 and 2025 releases and release-2026 have the same unguarded call. It came with 28697ebf "Improve the deform option" (!3141 (merged), #4607 (closed)), first released in 2024.0.

Steps to reproduce

spc-and-methanol.gro and spc-and-methanol.top from src/testutils/simulationdatabase, with

integrator       = steep
nsteps           = 10
cutoff-scheme    = Verlet
coulombtype      = Reaction-Field
rcoulomb         = 0.9
rvdw             = 0.9
deform           = 0 0 0 0 0.01 0
deform-init-flow = yes
gmx grompp -f em.mdp -c spc-and-methanol.gro -p spc-and-methanol.top -o em.tpr

The same happens with integrator = cg, l-bfgs and nm, and with tpi and tpic (which also need tc-grps, tau-t and ref-t). With integrator = md, the same input works.

What is the current bug behavior?

grompp crashes with a segmentation fault (exit status 139) after printing the degrees of freedom. There is no error message, and no .tpr is written.

What did you expect the correct behavior to be?

grompp writes the .tpr. For integrators without velocities, deform-init-flow is ignored, with a note.

Possible fixes

  • Only call deformInitFlow() when the state has velocities (state.hasEntry(StateEntry::V)).
  • In readir, for integrators without velocities, print a note that deform-init-flow is ignored. Also skip the note (and the error for gen-vel = yes) that asks for deform-init-flow = yes.

Related ordering issues in grompp

While looking at this, I found two more ordering issues. In both, grompp uses the state from -c before cont_status() replaces it with the frame from -t (https://gitlab.com/gromacs/gromacs/-/blob/dd255a67520b520ed49eb4bc5779cf4b7ca28891/src/gromacs/gmxpreprocess/grompp.cpp#L2770-2788).

  1. -t together with deform-init-flow = yes: the flow profile is computed from the -c coordinates and box. When velocities are then read from -t, they overwrite the flow profile. When they are not read (gen-vel = yes), the flow profile no longer matches the coordinates. Example on main with deform = 0 0 0 0 0.301 0 (so v_x = 0.1 z in a 3.01 nm box), -c with zero velocities, and -t a .trr frame of the same configuration shifted by 1 nm in z with zero velocities:

    • without -t: v_x = 0.1209, 0.1082, ... for z = 1.209, 1.082, ... (correct)
    • with -t, gen-vel = no: all velocities in the .tpr are zero
    • with -t, gen-vel = yes, gen-temp = 0: z = 2.209, 2.082, ... but v_x = 0.1209, 0.1082, ... The same happens in the usual setup where -c and -t hold the same frame, for example when a deform run is started from an equilibration without deform (-c npt.gro -t npt.cpt): with gen-vel = no, the velocities in the .tpr are exactly those from -t, without the flow profile and without a note, although the documentation says that the profile is added to the initial velocities. This exists since deform-init-flow was added.

  2. NVE (integrator = md, tcoupl = no) with -t and gen-vel = no: the temperature for the Verlet buffer is computed by calc_temp() from the -c velocities (https://gitlab.com/gromacs/gromacs/-/blob/dd255a67520b520ed49eb4bc5779cf4b7ca28891/src/gromacs/gmxpreprocess/grompp.cpp#L2662), but the run starts with the -t velocities. mdrun keeps this buffer for NVE. For spc-and-methanol, with -t frames whose velocities are 2x and 0.5x those in the .gro, grompp prints "will use the initial temperature of 1046.791 K" and sets rlist = 0.935 nm in all three cases (without -t and with each of the two frames). Giving the same velocities with -c gives 4187 K / rlist 1.017 nm and 262 K / 0.904 nm. With -c and -t from the same frame the velocities are the same, so this matters when they differ, for example when the -c file has no velocities: grompp then uses the fixed 10% buffer for T = 0 (rlist 0.99 nm here), even when the -t frame is hot. The atom density for the buffer also comes from the -c coordinates and box. This order exists since the automatic NVE buffer was added in 5.0.

A fix for both: read the -t frame before the buffer is determined. Then compute the buffer temperature from the final velocities, and add the flow profile last, from the final coordinates and box. Reading -t right after new_status() would probably be simplest. Velocities from a previous deform run already contain the flow profile, so a note for -t velocities combined with deform-init-flow = yes could help.

On main, velocities of frozen atoms are also not cleared after deformInitFlow() or after reading -t. That is already reported in a comment on #5714 and addressed in !6344 (merged), so I left it out here.

I can open an MR with the fix.