diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp
index d78af46249..0fed0297b4 100644
--- a/src/RHEO/pair_rheo.cpp
+++ b/src/RHEO/pair_rheo.cpp
@@ -218,12 +218,12 @@ void PairRHEO::compute(int eflag, int vflag)
         // Thermal Evolution
         if (thermal_flag) {
           dT = dot3(dx, dWij);
-          dT *= (kappai + kappaj) * (Ti - Tj) * rinv * rinv * voli * volj;
+          dT *= (kappai + kappaj) * (Ti - Tj) * rinv * rinv * voli * volj / rho0;
           heatflow[i] += dT;
 
           if (newton_pair || j < nlocal) {
             dT = dot3(dx, dWji);
-            dT *= (kappai + kappaj) * (Tj - Ti) * rinv * rinv * voli * volj;
+            dT *= (kappai + kappaj) * (Tj - Ti) * rinv * rinv * voli * volj / rho0;
             heatflow[j] -= dT;
           }
         }