diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
index 3c65ca9fb9..5599e861bb 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
@@ -88,8 +88,6 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-rho = thermo.rho();
-
 U -= rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
index b0a82278f0..fb89974a62 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
@@ -98,6 +98,8 @@ int main(int argc, char *argv[])
 
             turbulence->correct();
 
+            rho = thermo.rho();
+
             if (finalIter)
             {
                 mesh.data::remove("finalIteration");