Thermal conductivity/Heat Flux from Non-Equilibrium MD simulations

GROMACS version: 2022
GROMACS modification: No
Hello Gromacs users,

I am trying to retrieve the thermal conductivity from Non-Equilibrium MD simulations using Gromacs following the paper of Römer et al. (2012) (Cookie Absent). After setting up the simulation scenario of 2 thermostatted layers and conducting the simulation in the prescribed way, one “only” needs the change in kinetic energy of the particles in the thermostatted layers (before and after the thermostatting) to get the heat flux from the continuity equation, from which one gets the thermal conductivity.

This is where my problem comes in: How does one determine the change of the kinetic energy before and after thermostatting? Is there a way for me to print that out using a non-modified version of Gromacs? And if not, what needs to be modified to achieve that goal? They also explicitly stated that they were using a modified version in their paper.

I hope I provided enough information on my problem.
I would really appreciate it if someone could help me with this issue. Thank you in advance.