diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
index 0e39c2ff33a28e47a37595aa9a13a7183178d7da..927f2199ec9ffda267a312737879e327563a883e 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
@@ -172,9 +172,11 @@ void Foam::DsmcCloud<ParcelType>::initialise
                         cP.mass()
                     );
 
-                    scalar Ei =
-                        0.5*cP.internalDegreesOfFreedom()
-                       *kb*temperature;
+                    scalar Ei = equipartitionInternalEnergy
+                    (
+                        temperature,
+                        cP.internalDegreesOfFreedom()
+                    );
 
                     U += velocity;
 
@@ -690,6 +692,47 @@ Foam::vector Foam::DsmcCloud<ParcelType>::equipartitionLinearVelocity
 }
 
 
+template<class ParcelType>
+Foam::scalar Foam::DsmcCloud<ParcelType>::equipartitionInternalEnergy
+(
+    scalar temperature,
+    scalar iDof
+)
+{
+    scalar Ei = 0.0;
+
+    if
+    (
+        iDof < 2.0 + SMALL
+     && iDof > 2.0 - SMALL
+    )
+    {
+        // Special case for iDof = 2, i.e. diatomics;
+        Ei = -log(rndGen_.scalar01())*kb*temperature;
+    }
+    else
+    {
+        scalar a = 0.5*iDof - 1;
+
+        scalar energyRatio;
+
+        scalar P;
+
+        do
+        {
+            energyRatio = 10*rndGen_.scalar01();
+
+            P = pow((energyRatio/a), a)*exp(a - energyRatio);
+
+        } while (P < rndGen_.scalar01());
+
+        Ei = energyRatio*kb*temperature;
+    }
+
+    return Ei;
+}
+
+
 template<class ParcelType>
 void Foam::DsmcCloud<ParcelType>::dumpParticlePositions() const
 {
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
index 8e8fd5cbd74a667868ac619f8cca1bf6b9242b6d..b246b0499040fa15a68ab38ceb2543b8249e43f1 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
@@ -281,8 +281,17 @@ public:
                     scalar mass
                 );
 
-                //- From the Maxwellian distribution:
+                //- Generate a random internal energy, sampled from the
+                // equilibrium distribution (Bird eqn 11.22 and 11.23 and
+                // adapting code from DSMC3.FOR)
+                scalar equipartitionInternalEnergy
+                (
+                    scalar temperature,
+                    scalar internalDegreesOfFreedom
+                );
+
 
+                // From the Maxwellian distribution:
                 //- Average particle speed
                 inline scalar maxwellianAverageSpeed
                 (
diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C
index 93058e076e411949e115cb2b088203cb145b0a4b..f9c135718947fab6e8637c6ca0b24221864f2630 100644
--- a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C
+++ b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C
@@ -109,25 +109,28 @@ void Foam::DsmcParcel<ParcelType>::hitWallPatch
 {
     const constantProperties& constProps(td.cloud().constProps(typeId_));
 
+    scalar m = constProps.mass();
+
     // pre-interaction energy
-    scalar preIE = 0.5*constProps.mass()*(U_ & U_) + Ei_;
+    scalar preIE = 0.5*m*(U_ & U_) + Ei_;
 
     // pre-interaction momentum
-    vector preIMom = constProps.mass()*U_;
+    vector preIMom = m*U_;
 
     td.cloud().wallInteraction().correct
     (
         wpp,
         this->face(),
         U_,
-        constProps.mass()
+        Ei_,
+        typeId_
     );
 
     // post-interaction energy
-    scalar postIE = 0.5*constProps.mass()*(U_ & U_) + Ei_;
+    scalar postIE = 0.5*m*(U_ & U_) + Ei_;
 
     // post-interaction momentum
-    vector postIMom = constProps.mass()*U_;
+    vector postIMom = m*U_;
 
     label wppIndex = wpp.index();
 
diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C
index 1cf6fd90dcf09d46521d1079e8f34db1de2ce8d1..ff75da32b4e3c51e4279d2dac3615ba4039f7dab 100644
--- a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C
+++ b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C
@@ -81,9 +81,7 @@ Foam::scalar Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::energyRatio
                     ChiBMinusOne
                 );
         }
-    }
-
-    while (P < rndGen.scalar01());
+    } while (P < rndGen.scalar01());
 
     return energyRatio;
 }
diff --git a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
index a69f74aa994d88459e395570adecbe65be66700e..1367198efac558a052321a53321b1d01b022f00c 100644
--- a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
+++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
@@ -237,9 +237,11 @@ void Foam::FreeStream<CloudType>::inflow()
 
                 U += velocity_;
 
-                scalar Ei =
-                    0.5*cloud.constProps(typeId).internalDegreesOfFreedom()
-                   *CloudType::kb*temperature_;
+                scalar Ei = cloud.equipartitionInternalEnergy
+                (
+                    temperature_,
+                    cloud.constProps(typeId).internalDegreesOfFreedom()
+                );
 
                 cloud.addNewParcel
                 (
diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
index a6ed8e5718381c05491a1cf7c1aa859fa1472753..d8973fc5001601154b781516e05bb4f7a16a0f9d 100644
--- a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
+++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
@@ -54,7 +54,8 @@ void Foam::MaxwellianThermal<CloudType>::correct
     const wallPolyPatch& wpp,
     const label faceId,
     vector& U,
-    scalar mass
+    scalar& Ei,
+    label typeId
 )
 {
     label wppIndex = wpp.index();
@@ -102,6 +103,10 @@ void Foam::MaxwellianThermal<CloudType>::correct
 
     scalar T = cloud.T().boundaryField()[wppIndex][wppLocalFace];
 
+    scalar mass = cloud.constProps(typeId).mass();
+
+    scalar iDof = cloud.constProps(typeId).internalDegreesOfFreedom();
+
     U =
         sqrt(CloudType::kb*T/mass)
        *(
@@ -111,6 +116,8 @@ void Foam::MaxwellianThermal<CloudType>::correct
         );
 
     U += cloud.U().boundaryField()[wppIndex][wppLocalFace];
+
+    Ei = cloud.equipartitionInternalEnergy(T, iDof);
 }
 
 
diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.H b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.H
index 1542d20ff051b9e9f16afecf36b2e406a127c23b..a52a066aec21685e52947116f10f069b69b89b64 100644
--- a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.H
+++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.H
@@ -78,7 +78,8 @@ public:
             const wallPolyPatch& wpp,
             const label faceId,
             vector& U,
-            scalar mass
+            scalar& Ei,
+            label typeId
         );
 };
 
diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C b/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C
index 38e6621d30bcdfb697f461645cd5a686712ccce0..4de5507bf00f07b9df5d646cefb4f889e029d265 100644
--- a/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C
+++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C
@@ -56,7 +56,8 @@ void Foam::SpecularReflection<CloudType>::correct
     const wallPolyPatch& wpp,
     const label faceId,
     vector& U,
-    scalar mass
+    scalar& Ei,
+    label typeId
 )
 {
     vector nw = wpp.faceAreas()[wpp.whichFace(faceId)];
diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.H b/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.H
index 67d0adf8b074f6f737878f136130f5cb152797d5..d2208bd3503d29b6ca40a4805718a9f655882134 100644
--- a/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.H
+++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.H
@@ -76,7 +76,8 @@ public:
             const wallPolyPatch& wpp,
             const label faceId,
             vector& U,
-            scalar mass
+            scalar& Ei,
+            label typeId
         );
 };
 
diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/WallInteractionModel/WallInteractionModel.H b/src/lagrangian/dsmc/submodels/WallInteractionModel/WallInteractionModel/WallInteractionModel.H
index d8667201a7ab52a3a214ffceebf8daaf767c4d43..582bf2b1065ce26ec0113398357fed1a75d92a69 100644
--- a/src/lagrangian/dsmc/submodels/WallInteractionModel/WallInteractionModel/WallInteractionModel.H
+++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/WallInteractionModel/WallInteractionModel.H
@@ -130,7 +130,8 @@ public:
             const wallPolyPatch& wpp,
             const label faceId,
             vector& U,
-            scalar mass
+            scalar& Ei,
+            label typeId
         ) = 0;
 };