Measuring temperature by the average kinetic energy per degree of freedom to allow for a mix of linear, mono-atomic and 6dof molecules in a case.

src/lagrangian/molecularDynamics/molecule/mdTools/meanMomentumEnergyAndNMols.H

@@ -52,6 +52,8 @@ scalar singleStepTotalrDotf = 0.0;
 
 label singleStepNMols = molecules.size();
 
+label singleStepDOFs = 0;
+
 {
     IDLList<molecule>::iterator mol(molecules.begin());
 
@@ -85,9 +87,11 @@ label singleStepNMols = molecules.size();
     {
         label molId = mol().id();
 
-        scalar molMass(molecules.constProps(molId).mass());
+        const molecule::constantProperties cP(molecules.constProps(molId));
 
-        const diagTensor& molMoI(molecules.constProps(molId).momentOfInertia());
+        scalar molMass(cP.mass());
+
+        const diagTensor& molMoI(cP.momentOfInertia());
 
         const vector& molV(mol().v());
 
@@ -112,6 +116,8 @@ label singleStepNMols = molecules.size();
         singleStepTotalPE += mol().potentialEnergy();
 
         singleStepTotalrDotf += tr(mol().rf());
+
+        singleStepDOFs += cP.degreesOfFreedom();
     }
 }
 
@@ -134,46 +140,49 @@ if (Pstream::parRun())
     reduce(singleStepTotalrDotf, sumOp<scalar>());
 
     reduce(singleStepNMols, sumOp<label>());
+
+    reduce(singleStepDOFs, sumOp<label>());
 }
 
 if (singleStepNMols)
 {
-    Info<< "Number of mols in system = "
+    Info<< "Number of molecules in system = "
         << singleStepNMols << nl
         << "Overall number density = "
         << singleStepNMols/meshVolume << nl
         << "Overall mass density = "
         << singleStepTotalMass/meshVolume << nl
-        << "Average linear momentum per mol = "
+        << "Average linear momentum per molecule = "
         << singleStepTotalLinearMomentum/singleStepNMols << ' '
         << mag(singleStepTotalLinearMomentum)/singleStepNMols << nl
-        << "Average angular momentum per mol = "
+        << "Average angular momentum per molecule = "
         << singleStepTotalAngularMomentum << ' '
         << mag(singleStepTotalAngularMomentum)/singleStepNMols << nl
         << "Maximum |velocity| = "
         << singleStepMaxVelocityMag << nl
-        << "Average linear KE per mol = "
+        << "Average linear KE per molecule = "
         << singleStepTotalLinearKE/singleStepNMols << nl
-        << "Average angular KE per mol = "
+        << "Average angular KE per molecule = "
         << singleStepTotalAngularKE/singleStepNMols << nl
-        << "Average PE per mol = "
+        << "Average PE per molecule = "
         << singleStepTotalPE/singleStepNMols << nl
-        << "Average TE per mol = "
+        << "Average TE per molecule = "
         <<
         (
             singleStepTotalLinearKE
-            + singleStepTotalAngularKE
-            + singleStepTotalPE
+          + singleStepTotalAngularKE
+          + singleStepTotalPE
         )
         /singleStepNMols
         << endl;
 
-//     Info << singleStepNMols << " "
-// //     << singleStepTotalMomentum/singleStepTotalMass << " "
-//     << singleStepMaxVelocityMag << " "
-//     << singleStepTotalKE/singleStepNMols << " "
-//     << singleStepTotalPE/singleStepNMols << " "
-//     << (singleStepTotalKE + singleStepTotalPE)/singleStepNMols << endl;
+        // Info << singleStepNMols << " "
+        //     << singleStepTotalMomentum/singleStepTotalMass << " "
+        //     << singleStepMaxVelocityMag << " "
+        //     << singleStepTotalKE/singleStepNMols << " "
+        //     << singleStepTotalPE/singleStepNMols << " "
+        //     << (singleStepTotalKE + singleStepTotalPE)
+        //        /singleStepNMols << endl;
 }
 else
 {

src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H

@@ -47,17 +47,17 @@ accumulatedTotalrDotfSum += singleStepTotalrDotf;
 
 accumulatedNMols += singleStepNMols;
 
+accumulatedDOFs += singleStepDOFs;
+
 if (runTime.outputTime())
 {
-    // calculate averages
-
     if (accumulatedNMols)
     {
         Info << "calculating averages" << endl;
 
         averageTemperature =
         (
-            2.0/(6.0 * moleculeCloud::kb * accumulatedNMols)
+            2.0/(moleculeCloud::kb * accumulatedDOFs)
             *
             (
                 accumulatedTotalLinearKE + accumulatedTotalAngularKE
@@ -79,8 +79,6 @@ if (runTime.outputTime())
             meshVolume
         );
 
-        // output values
-
         Info << "----------------------------------------" << nl
             << "Averaged properties" << nl
             << "Average |velocity| = "
@@ -91,12 +89,10 @@ if (runTime.outputTime())
     }
     else
     {
-        Info << "Not averaging temperature and pressure: "
+        Info<< "Not averaging temperature and pressure: "
             << "no molecules in system" << endl;
     }
 
-    // reset counters
-
     accumulatedTotalLinearMomentum = vector::zero;
 
     accumulatedTotalMass = 0.0;
@@ -110,6 +106,8 @@ if (runTime.outputTime())
     accumulatedTotalrDotfSum = 0.0;
 
     accumulatedNMols = 0;
+
+    accumulatedDOFs = 0;
 }
 
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressureVariables.H

@@ -44,6 +44,8 @@ scalar accumulatedTotalrDotfSum = 0.0;
 
 label accumulatedNMols = 0;
 
+label accumulatedDOFs = 0;
+
 scalar averageTemperature = 0.0;
 
 scalar averagePressure = 0.0;

src/lagrangian/molecularDynamics/molecule/molecule/molecule.H

@@ -141,6 +141,8 @@ public:
 
         inline bool pointMolecule() const;
 
+        inline label degreesOfFreedom() const;
+
         inline scalar mass() const;
 
         inline label nSites() const;

src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H

@@ -423,7 +423,7 @@ Foam::molecule::constantProperties::momentOfInertia() const
 
 inline bool Foam::molecule::constantProperties::linearMolecule() const
 {
-    return (momentOfInertia_.xx() < 0);
+    return ((momentOfInertia_.xx() < 0) && (momentOfInertia_.yy() > 0));
 }
 
 
@@ -433,6 +433,23 @@ inline bool Foam::molecule::constantProperties::pointMolecule() const
 }
 
 
+inline Foam::label Foam::molecule::constantProperties::degreesOfFreedom() const
+{
+    if (linearMolecule())
+    {
+        return 5;
+    }
+    else if (pointMolecule())
+    {
+        return 3;
+    }
+    else
+    {
+        return 6;
+    }
+}
+
+
 inline Foam::scalar Foam::molecule::constantProperties::mass() const
 {
     return mass_;