  In the context of forthcoming galaxy surveys, to ensure unbiased constraints
on cosmology and gravity when using non-linear structure information,
percent-level accuracy is required when modelling the power spectrum. This
calls for frameworks that can accurately capture the relevant physical effects,
while allowing for deviations from $\Lambda$CDM. Massive neutrino and baryonic
physics are two of the most relevant such effects. We present an integration of
the halo model reaction frameworks for massive neutrinos and
beyond-$\Lambda$CDM cosmologies. The integrated halo model reaction, combined
with a pseudo power spectrum modelled by HMCode2020 is then compared against
$N$-body simulations that include both massive neutrinos and an $f(R)$
modification to gravity. We find that the framework is 4% accurate down to at
least $k\approx 3 \, h/{\rm Mpc}$ for a modification to gravity of $|f_{\rm
R0}|\leq 10^{-5}$ and for the total neutrino mass $M_\nu \equiv \sum m_\nu \leq
0.15$ eV. We also find that the framework is 4% consistent with EuclidEmulator2
as well as the Bacco emulator for \B{most of the considered} $\nu w$CDM
cosmologies down to at least $k \approx 3 \, h$/Mpc. Finally, we compare
against hydrodynamical simulations employing HMCode2020's baryonic feedback
modelling on top of the halo model reaction. For $\nu \Lambda$CDM cosmologies
we find 2% accuracy for $M_\nu \leq 0.48$eV down to at least $k\approx 5h$/Mpc.
Similar accuracy is found when comparing to $\nu w$CDM hydrodynamical
simulations with $M_\nu = 0.06$eV. This offers the first non-linear,
theoretically general means of accurately including massive neutrinos for
beyond-$\Lambda$CDM cosmologies, and further suggests that baryonic, massive
neutrino and dark energy physics can be reliably modelled independently.
