Adding pressure measurement via r dot f and temperature measurement and control including rotational energy. Adding random orientation on initialisation. Modifying constantProperties construction to detect point masses. Tidy up hitWallPatch function to remove commented out stochastic wall code.

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;

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;
 

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;
 

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;
+    }
 }
 
 

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
+        centreOfMass += siteReferencePositions_[i]*siteMasses[i];
-        // position
+    }
 
-        vector centreOfMass(vector::zero);
+    centreOfMass /= mass_;
 
-        forAll(siteReferencePositions_, i)
+    siteReferencePositions_ -= centreOfMass;
-        {
-            centreOfMass += siteReferencePositions_[i]*siteMasses[i];
-        }
 
-        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);
+        momOfInertia += siteMasses[i]*tensor
-
+        (
-        forAll(siteReferencePositions_, i)
+            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(),
-            const vector& p(siteReferencePositions_[i]);
+            -p.z()*p.x(), -p.z()*p.y(), p.x()*p.x() + p.y()*p.y()
-
+        );
-            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()
-            );
-        }
 
+    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
+        // Calculate the transformation between the principle axes to the
-        // axes
+        // 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
+        // The rotation was around the centre of mass but remove any
-        // that have crept in due to floating point errors
+        // 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
+        // Calculate the moment of inertia in the principle component
-        // frame
+        // 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);
+    }
 }
 
 

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;
     }
 }
 

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() += rsIsJ * fsIsJ;
-
-                    // molJ->rf() += rIJ * fIJ;
                 }
             }
 
@@ -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() += rsIsJ * fsIsJ;
-
-                    // molJ->rf() += rIJ * fIJ;
                 }
             }
         }
@@ -200,9 +196,7 @@ inline void Foam::moleculeCloud::evaluatePair
 
                     molReal->potentialEnergy() += 0.5*potentialEnergy;
 
-                    // What to do here?
+                    molReal->rf() += rsRealsRef * fsRealsRef;
-
-                    // molReal->rf() += rRealRef * fRealRef;
 
                 }
             }
@@ -233,9 +227,7 @@ inline void Foam::moleculeCloud::evaluatePair
 
                     molReal->potentialEnergy() += 0.5*potentialEnergy;
 
-                    // What to do here?
+                    molReal->rf() += rsRealsRef * fsRealsRef;
-
-                    // molReal->rf() += rRealRef * fRealRef;
                 }
             }
         }