Discrepancy in Heat Flux Calculation Using Dual-Thermostat Approach for Thermal Conductivity

User discussions
1

GROMACS version: 2022.3
GROMACS modification: No
We are currently working on computing thermal conductivity using the dual thermostat approach in GROMACS. In this method, the simulation box is enlarged and divided into N slabs, with two slabs equally spaced (regarding periodicity) coupled to high-temperature and low-temperature thermostats. This setup produces a linear temperature profile, allowing us to calculate the thermal conductivity coefficient by dividing the heat flux by the temperature gradient. More details on this technique can be found in the 2012 paper of Algaer and Müller-Plathe (10.1080/1539445X.2011.599699).

However, I’m encountering two issues:

  1. The heat fluxes obtained from the two thermostats are not equal, even after running long simulations (40 ns).

  2. Currently, we’re computing the heat flux using the Berendsen expression for the change in kinetic energy. We would like to switch to using the v-rescale or Nose-Hoover thermostat, but we was not able to find an explicit expression for the change in kinetic energy with these thermostats. Is there an expression available, or is there a way to retrieve it using GROMACS commands or from GROMACS code?

Any insights or suggestions would be greatly appreciated.