diff --git a/applications/solvers/molecularDynamics/mdFoam/mdFoam.C b/applications/solvers/molecularDynamics/mdFoam/mdFoam.C
index 4357c5c0b1..f65424b437 100644
--- a/applications/solvers/molecularDynamics/mdFoam/mdFoam.C
+++ b/applications/solvers/molecularDynamics/mdFoam/mdFoam.C
@@ -46,20 +46,40 @@ int main(int argc, char *argv[])
 
 #   include "temperatureAndPressureVariables.H"
 
+    label nAveragingSteps = 0;
+
     Info << "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
         runTime++;
 
+        nAveragingSteps++;
+
         Info << "Time = " << runTime.timeName() << endl;
 
         molecules.evolve();
 
 #       include "meanMomentumEnergyAndNMols.H"
 
+#       include "temperatureAndPressure.H"
+
+        if (runTime.outputTime())
+        {
+            molecules.applyConstraintsAndThermostats
+            (
+                485.4,
+                averageTemperature
+            );
+        }
+
         runTime.write();
 
+        if (runTime.outputTime())
+        {
+            nAveragingSteps = 0;
+        }
+
         Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H
index 24d6fb88e5..ddf99c30d5 100644
--- a/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H
@@ -33,11 +33,13 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-accumulatedTotalMomentum += singleStepTotalMomentum;
+accumulatedTotalLinearMomentum += singleStepTotalLinearMomentum;
 
 accumulatedTotalMass += singleStepTotalMass;
 
-accumulatedTotalKE += singleStepTotalKE;
+accumulatedTotalLinearKE += singleStepTotalLinearKE;
+
+accumulatedTotalAngularKE += singleStepTotalAngularKE;
 
 accumulatedTotalPE += singleStepTotalPE;
 
@@ -55,12 +57,12 @@ if (runTime.outputTime())
 
         averageTemperature =
         (
-            2.0/(3.0 * moleculeCloud::kb * accumulatedNMols)
+            2.0/(6.0 * moleculeCloud::kb * accumulatedNMols)
             *
             (
-                accumulatedTotalKE
+                accumulatedTotalLinearKE + accumulatedTotalAngularKE
                 -
-                0.5*magSqr(accumulatedTotalMomentum)/accumulatedTotalMass
+                0.5*magSqr(accumulatedTotalLinearMomentum)/accumulatedTotalMass
             )
         );
 
@@ -82,7 +84,7 @@ if (runTime.outputTime())
         Info << "----------------------------------------" << nl
             << "Averaged properties" << nl
             << "Average |velocity| = "
-                << mag(accumulatedTotalMomentum)/accumulatedTotalMass
+                << mag(accumulatedTotalLinearMomentum)/accumulatedTotalMass
                 << " m/s" << nl
             << "Average temperature = "
                 << averageTemperature << " K" << nl
@@ -98,11 +100,13 @@ if (runTime.outputTime())
 
     // reset counters
 
-    accumulatedTotalMomentum = vector::zero;
+    accumulatedTotalLinearMomentum = vector::zero;
 
     accumulatedTotalMass = 0.0;
 
-    accumulatedTotalKE = 0.0;
+    accumulatedTotalLinearKE = 0.0;
+
+    accumulatedTotalAngularKE = 0.0;
 
     accumulatedTotalPE = 0.0;
 
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressureVariables.H b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressureVariables.H
index 283535c7a4..45b5c0d221 100644
--- a/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressureVariables.H
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressureVariables.H
@@ -30,11 +30,13 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-vector accumulatedTotalMomentum(vector::zero);
+vector accumulatedTotalLinearMomentum(vector::zero);
 
 scalar accumulatedTotalMass = 0.0;
 
-scalar accumulatedTotalKE = 0.0;
+scalar accumulatedTotalAngularKE = 0.0;
+
+scalar accumulatedTotalLinearKE = 0.0;
 
 scalar accumulatedTotalPE = 0.0;
 
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
index 2d8cd8f00b..5f90b7856e 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
@@ -246,38 +246,12 @@ void Foam::molecule::hitWallPatch
     nw /= mag(nw);
 
     scalar vn = v_ & nw;
-//     vector vt = v_ - vn*nw;
-
-//     Random rand(clock::getTime());
-
-//     scalar tmac = 0.8;
-
-//     scalar wallTemp = 2.5;
-
-//     if (rand.scalar01() < tmac)
-//     {
-//         // Diffuse reflection
-//
-//         vector tw1 = vt/mag(vt);
-//
-//         vector tw2 = nw ^ tw1;
-//
-//         V_ = sqrt(wallTemp/mass_)*rand.GaussNormal()*tw1
-//                 + sqrt(wallTemp/mass_)*rand.GaussNormal()*tw2
-//                 - mag(sqrt(wallTemp/mass_)*rand.GaussNormal())*nw;
-//     }
-
-//     else
-//     {
-        // Specular reflection
-
-        if (vn > 0)
-        {
-            v_ -= 2*vn*nw;
-        }
-
-//     }
 
+    // Specular reflection
+    if (vn > 0)
+    {
+        v_ -= 2*vn*nw;
+    }
 }
 
 
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
index 7d3f5c3aa7..c1d69a87c4 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
@@ -63,46 +63,47 @@ inline Foam::molecule::constantProperties::constantProperties
 
     mass_ = sum(siteMasses);
 
+    // Calculate the centre of mass of the body and subtract it from each
+    // position
+
+    vector centreOfMass(vector::zero);
+
+    forAll(siteReferencePositions_, i)
     {
-        // Calculate the centre of mass of the body and subtract it from each
-        // position
+        centreOfMass += siteReferencePositions_[i]*siteMasses[i];
+    }
 
-        vector centreOfMass(vector::zero);
+    centreOfMass /= mass_;
 
-        forAll(siteReferencePositions_, i)
-        {
-            centreOfMass += siteReferencePositions_[i]*siteMasses[i];
-        }
+    siteReferencePositions_ -= centreOfMass;
 
-        centreOfMass /= mass_;
+    // Calculate the inertia tensor in the current orientation
 
-        siteReferencePositions_ -= centreOfMass;
+    tensor momOfInertia(tensor::zero);
 
-        // Calculate the inertia tensor in the current orientation
+    forAll(siteReferencePositions_, i)
+    {
+        const vector& p(siteReferencePositions_[i]);
 
-        tensor I(tensor::zero);
-
-        forAll(siteReferencePositions_, i)
-        {
-            const vector& p(siteReferencePositions_[i]);
-
-            I += siteMasses[i]*tensor
-            (
-                p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(),
-                -p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(),
-                -p.z()*p.x(), -p.z()*p.y(), p.x()*p.x() + p.y()*p.y()
-            );
-        }
+        momOfInertia += siteMasses[i]*tensor
+        (
+            p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(),
+            -p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(),
+            -p.z()*p.x(), -p.z()*p.y(), p.x()*p.x() + p.y()*p.y()
+        );
+    }
 
+    if(eigenValues(momOfInertia).z() > VSMALL)
+    {
         // Normalise the inertia tensor magnitude to avoid SMALL numbers in the
         // components causing problems
 
-        I /= eigenValues(I).z();
+        momOfInertia /= eigenValues(momOfInertia).z();
 
-        tensor e = eigenVectors(I);
+        tensor e = eigenVectors(momOfInertia);
 
-        // Calculate the transformation between the principle axes to the global
-        // axes
+        // Calculate the transformation between the principle axes to the
+        // global axes
 
         tensor Q =
         vector(1,0,0)*e.x() + vector(0,1,0)*e.y() + vector(0,0,1)*e.z();
@@ -110,13 +111,11 @@ inline Foam::molecule::constantProperties::constantProperties
         // Transform the site positions
 
         siteReferencePositions_ = (Q & siteReferencePositions_);
-    }
 
-    {
         // Recalculating the moment of inertia with the new site positions
 
-        // The rotation was around the centre of mass but remove any components
-        // that have crept in due to floating point errors
+        // The rotation was around the centre of mass but remove any
+        // components that have crept in due to floating point errors
 
         vector centreOfMass(vector::zero);
 
@@ -129,16 +128,16 @@ inline Foam::molecule::constantProperties::constantProperties
 
         siteReferencePositions_ -= centreOfMass;
 
-        // Calculate the moment of inertia in the principle component reference
-        // frame
+        // Calculate the moment of inertia in the principle component
+        // reference frame
 
-        tensor I(tensor::zero);
+        tensor momOfInertia(tensor::zero);
 
         forAll(siteReferencePositions_, i)
         {
             const vector& p(siteReferencePositions_[i]);
 
-            I += siteMasses[i]*tensor
+            momOfInertia += siteMasses[i]*tensor
             (
                 p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(),
                 -p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(),
@@ -146,9 +145,19 @@ inline Foam::molecule::constantProperties::constantProperties
             );
         }
 
-        momentOfInertia_ = diagTensor(I.xx(), I.yy(), I.zz());
+        momentOfInertia_ = diagTensor
+        (
+            momOfInertia.xx(),
+            momOfInertia.yy(),
+            momOfInertia.zz()
+        );
     }
+    else
+    {
+        Info<< nl << "Molecule " << dict.name()  << " is a point mass." << endl;
 
+        momentOfInertia_ = diagTensor(0, 0, 0);
+    }
 }
 
 
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
index c15a629417..1d0e2da8d1 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
@@ -982,6 +982,25 @@ void Foam::moleculeCloud::createMolecule
         cP.momentOfInertia()
     );
 
+    scalar phi(rndGen_.scalar01()*mathematicalConstant::twoPi);
+
+    scalar theta(rndGen_.scalar01()*mathematicalConstant::twoPi);
+
+    scalar psi(rndGen_.scalar01()*mathematicalConstant::twoPi);
+
+    const tensor Q
+    (
+        cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi),
+        cos(psi)*sin(phi) + cos(theta)*cos(phi)*sin(psi),
+        sin(psi)*sin(theta),
+        - sin(psi)*cos(phi) - cos(theta)*sin(phi)*cos(psi),
+        - sin(psi)*sin(phi) + cos(theta)*cos(phi)*cos(psi),
+        cos(psi)*sin(theta),
+        sin(theta)*sin(phi),
+        - sin(theta)*cos(phi),
+        cos(theta)
+    );
+
     addParticle
     (
         new molecule
@@ -989,7 +1008,7 @@ void Foam::moleculeCloud::createMolecule
             cloud,
             position,
             cell,
-            I,
+            Q,
             v,
             vector::zero,
             pi,
@@ -1124,6 +1143,8 @@ void Foam::moleculeCloud::applyConstraintsAndThermostats
     for (mol = this->begin(); mol != this->end(); ++mol)
     {
         mol().v() *= temperatureCorrectionFactor;
+
+        mol().pi() *= temperatureCorrectionFactor;
     }
 }
 
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
index ca095562cd..72e5617a8f 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
@@ -91,11 +91,9 @@ inline void Foam::moleculeCloud::evaluatePair
 
                     molJ->potentialEnergy() += 0.5*potentialEnergy;
 
-                    // What to do here?
+                    molI->rf() += rsIsJ * fsIsJ;
 
-                    // molI->rf() += rIJ * fIJ;
-
-                    // molJ->rf() += rIJ * fIJ;
+                    molJ->rf() += rsIsJ * fsIsJ;
                 }
             }
 
@@ -128,11 +126,9 @@ inline void Foam::moleculeCloud::evaluatePair
 
                     molJ->potentialEnergy() += 0.5*potentialEnergy;
 
-                    // What to do here?
+                    molI->rf() += rsIsJ * fsIsJ;
 
-                    // molI->rf() += rIJ * fIJ;
-
-                    // molJ->rf() += rIJ * fIJ;
+                    molJ->rf() += rsIsJ * fsIsJ;
                 }
             }
         }
@@ -200,9 +196,7 @@ inline void Foam::moleculeCloud::evaluatePair
 
                     molReal->potentialEnergy() += 0.5*potentialEnergy;
 
-                    // What to do here?
-
-                    // molReal->rf() += rRealRef * fRealRef;
+                    molReal->rf() += rsRealsRef * fsRealsRef;
 
                 }
             }
@@ -233,9 +227,7 @@ inline void Foam::moleculeCloud::evaluatePair
 
                     molReal->potentialEnergy() += 0.5*potentialEnergy;
 
-                    // What to do here?
-
-                    // molReal->rf() += rRealRef * fRealRef;
+                    molReal->rf() += rsRealsRef * fsRealsRef;
                 }
             }
         }