Dealing with molecules that do not have all 3 diagonal components of the...

Dealing with molecules that do not have all 3 diagonal components of the principal axis inertia tensor, i.e. point masses (mono-atomics) and linear molecules (diatomics and CO2 for example).

applications/solvers/molecularDynamics/mdEquilibrationFoam/Make/options

 EXE_INC = \
     -I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \
     -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
@@ -10,4 +11,6 @@ EXE_LIBS = \
     -lfiniteVolume \
     -llagrangian \
     -lmolecule \
-    -lpotential
+    -lpotential \
+    -lmolecularMeasurements
+

applications/solvers/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C

@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    mdEquilibrationFOAM
+    mdEquilibrationFoam
 
 Description
     Equilibrates and/or preconditions MD systems
@@ -40,9 +40,9 @@ int main(int argc, char *argv[])
 #   include "createTime.H"
 #   include "createMesh.H"
 
-    moleculeCloud molecules(mesh);
+    potential pot(mesh);
 
-    molecules.removeHighEnergyOverlaps();
+    moleculeCloud molecules(mesh, pot);
 
 #   include "temperatureAndPressureVariables.H"
 
@@ -60,7 +60,7 @@ int main(int argc, char *argv[])
 
         Info << "Time = " << runTime.timeName() << endl;
 
-        molecules.integrateEquationsOfMotion();
+        molecules.evolve();
 
 #       include "meanMomentumEnergyAndNMols.H"
 

applications/solvers/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H

-Info<< "Reading MD Equilibration Dictionary" << nl << endl;
+Info<< nl << "Reading MD Equilibration Dictionary" << nl << endl;
 
 IOdictionary mdEquilibrationDict
 (

applications/solvers/molecularDynamics/mdFoam/mdFoam.C

@@ -64,15 +64,6 @@ int main(int argc, char *argv[])
 
 #       include "temperatureAndPressure.H"
 
-        if (runTime.outputTime())
-        {
-            molecules.applyConstraintsAndThermostats
-            (
-                485.4,
-                averageTemperature
-            );
-        }
-
         runTime.write();
 
         if (runTime.outputTime())

src/lagrangian/molecularDynamics/molecule/molecule/molecule.C

@@ -123,29 +123,38 @@ bool Foam::molecule::move(molecule::trackData& td)
         // but after tracking stage, i.e. rotation carried once linear motion
         // complete.
 
-        const diagTensor& momentOfInertia(constProps.momentOfInertia());
-
-        tensor R;
-
-        R = rotationTensorX(0.5*deltaT*pi_.x()/momentOfInertia.xx());
-        pi_ = pi_ & R;
-        Q_ = Q_ & R;
-
-        R = rotationTensorY(0.5*deltaT*pi_.y()/momentOfInertia.yy());
-        pi_ = pi_ & R;
-        Q_ = Q_ & R;
-
-        R = rotationTensorZ(deltaT*pi_.z()/momentOfInertia.zz());
-        pi_ = pi_ & R;
-        Q_ = Q_ & R;
-
-        R = rotationTensorY(0.5*deltaT*pi_.y()/momentOfInertia.yy());
-        pi_ = pi_ & R;
-        Q_ = Q_ & R;
-
-        R = rotationTensorX(0.5*deltaT*pi_.x()/momentOfInertia.xx());
-        pi_ = pi_ & R;
-        Q_ = Q_ & R;
+        if (!constProps.pointMolecule())
+        {
+            const diagTensor& momentOfInertia(constProps.momentOfInertia());
+
+            tensor R;
+
+            if (!constProps.linearMolecule())
+            {
+                R = rotationTensorX(0.5*deltaT*pi_.x()/momentOfInertia.xx());
+                pi_ = pi_ & R;
+                Q_ = Q_ & R;
+            }
+
+            R = rotationTensorY(0.5*deltaT*pi_.y()/momentOfInertia.yy());
+            pi_ = pi_ & R;
+            Q_ = Q_ & R;
+
+            R = rotationTensorZ(deltaT*pi_.z()/momentOfInertia.zz());
+            pi_ = pi_ & R;
+            Q_ = Q_ & R;
+
+            R = rotationTensorY(0.5*deltaT*pi_.y()/momentOfInertia.yy());
+            pi_ = pi_ & R;
+            Q_ = Q_ & R;
+
+            if (!constProps.linearMolecule())
+            {
+                R = rotationTensorX(0.5*deltaT*pi_.x()/momentOfInertia.xx());
+                pi_ = pi_ & R;
+                Q_ = Q_ & R;
+            }
+        }
 
         setSitePositions(constProps);
     }
@@ -172,6 +181,20 @@ bool Foam::molecule::move(molecule::trackData& td)
         v_ += 0.5*deltaT*a_;
 
         pi_ += 0.5*deltaT*tau_;
+
+        if (constProps.pointMolecule())
+        {
+            tau_ = vector::zero;
+
+            pi_ = vector::zero;
+        }
+
+        if (constProps.linearMolecule())
+        {
+            tau_.x() = 0.0;
+
+            pi_.x() = 0.0;
+        }
     }
     else
     {

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

@@ -104,6 +104,8 @@ public:
                 const List<word>& pairPotSiteIds
             );
 
+            bool linearMoleculeTest() const;
+
     public:
 
         inline constantProperties();
@@ -131,6 +133,10 @@ public:
 
         inline const diagTensor& momentOfInertia() const;
 
+        inline bool linearMolecule() const;
+
+        inline bool pointMolecule() const;
+
         inline scalar mass() const;
 
         inline label nSites() const;

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

@@ -63,11 +63,11 @@ inline Foam::molecule::constantProperties::constantProperties
 
     mass_ = sum(siteMasses);
 
+    vector centreOfMass(vector::zero);
+
     // Calculate the centre of mass of the body and subtract it from each
     // position
 
-    vector centreOfMass(vector::zero);
-
     forAll(siteReferencePositions_, i)
     {
         centreOfMass += siteReferencePositions_[i]*siteMasses[i];
@@ -77,28 +77,94 @@ inline Foam::molecule::constantProperties::constantProperties
 
     siteReferencePositions_ -= centreOfMass;
 
-    // Calculate the inertia tensor in the current orientation
+    if(siteIds_.size() == 1)
+    {
+        // Single site molecule - no rotational motion.
 
-    tensor momOfInertia(tensor::zero);
+        siteReferencePositions_[0] = vector::zero;
 
-    forAll(siteReferencePositions_, i)
+        momentOfInertia_ = diagTensor(-1, -1, -1);
+    }
+    else if(linearMoleculeTest())
     {
-        const vector& p(siteReferencePositions_[i]);
+        // Linear molecule.
 
-        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()
-        );
-    }
+        Info<< nl << "Linear molecule." << endl;
+
+        vector dir = siteReferencePositions_[1] - siteReferencePositions_[0];
+
+        dir /= mag(dir);
 
-    if(eigenValues(momOfInertia).z() > VSMALL)
+        tensor Q = rotationTensor(dir, vector(1,0,0));
+
+        siteReferencePositions_ = (Q & siteReferencePositions_);
+
+        // The rotation was around the centre of mass but remove any
+        // components that have crept in due to floating point errors
+
+        centreOfMass = vector::zero;
+
+        forAll(siteReferencePositions_, i)
+        {
+            centreOfMass += siteReferencePositions_[i]*siteMasses[i];
+        }
+
+        centreOfMass /= mass_;
+
+        siteReferencePositions_ -= centreOfMass;
+
+        diagTensor momOfInertia = diagTensor::zero;
+
+        forAll(siteReferencePositions_, i)
+        {
+            const vector& p(siteReferencePositions_[i]);
+
+            momOfInertia += siteMasses[i]*diagTensor
+            (
+                0, p.x()*p.x(), p.x()*p.x()
+            );
+        }
+
+        momentOfInertia_ = diagTensor
+         (
+             -1,
+             momOfInertia.yy(),
+             momOfInertia.zz()
+         );
+    }
+    else
     {
+        // Fully 6DOF molecule
+
+        // Calculate the inertia tensor in the current orientation
+
+        tensor momOfInertia(tensor::zero);
+
+        forAll(siteReferencePositions_, i)
+        {
+            const vector& p(siteReferencePositions_[i]);
+
+            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).x() < VSMALL)
+        {
+            FatalErrorIn("molecule::constantProperties::constantProperties")
+                << "An eigenvalue of the inertia tensor is zero, the molecule "
+                << dict.name()
+                << " is not a valid 6DOF shape."
+                << nl << abort(FatalError);
+        }
+
         // Normalise the inertia tensor magnitude to avoid SMALL numbers in the
         // components causing problems
 
-        momOfInertia /= eigenValues(momOfInertia).z();
+        momOfInertia /= eigenValues(momOfInertia).x();
 
         tensor e = eigenVectors(momOfInertia);
 
@@ -117,7 +183,7 @@ inline Foam::molecule::constantProperties::constantProperties
         // 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);
+        centreOfMass = vector::zero;
 
         forAll(siteReferencePositions_, i)
         {
@@ -131,7 +197,7 @@ inline Foam::molecule::constantProperties::constantProperties
         // Calculate the moment of inertia in the principle component
         // reference frame
 
-        tensor momOfInertia(tensor::zero);
+        momOfInertia = tensor::zero;
 
         forAll(siteReferencePositions_, i)
         {
@@ -152,12 +218,6 @@ inline Foam::molecule::constantProperties::constantProperties
             momOfInertia.zz()
         );
     }
-    else
-    {
-        Info<< nl << "Molecule " << dict.name()  << " is a point mass." << endl;
-
-        momentOfInertia_ = diagTensor(0, 0, 0);
-    }
 }
 
 
@@ -236,6 +296,38 @@ inline void Foam::molecule::constantProperties::setInteracionSiteBools
 }
 
 
+inline bool Foam::molecule::constantProperties::linearMoleculeTest() const
+{
+    if (siteIds_.size() == 2)
+    {
+        return true;
+    }
+
+    vector refDir = siteReferencePositions_[1] - siteReferencePositions_[0];
+
+    refDir /= mag(refDir);
+
+    for
+    (
+        label i = 2;
+        i < siteReferencePositions_.size();
+        i++
+    )
+    {
+        vector dir = siteReferencePositions_[i] - siteReferencePositions_[i-1];
+
+        dir /= mag(dir);
+
+        if (mag(refDir & dir) < 1 - SMALL)
+        {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+
 // * * * * * * * constantProperties Member Functions * * * * * * * * * * * * //
 
 inline const Foam::Field<Foam::vector>&
@@ -322,6 +414,18 @@ Foam::molecule::constantProperties::momentOfInertia() const
 }
 
 
+inline bool Foam::molecule::constantProperties::linearMolecule() const
+{
+    return (momentOfInertia_.xx() < 0);
+}
+
+
+inline bool Foam::molecule::constantProperties::pointMolecule() const
+{
+    return (momentOfInertia_.zz() < 0);
+}
+
+
 inline Foam::scalar Foam::molecule::constantProperties::mass() const
 {
     return mass_;

src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C

@@ -976,30 +976,33 @@ void Foam::moleculeCloud::createMolecule
 
     v += bulkVelocity;
 
-    vector pi = equipartitionAngularMomentum
-    (
-        temperature,
-        cP.momentOfInertia()
-    );
+    vector pi = vector::zero;
 
-    scalar phi(rndGen_.scalar01()*mathematicalConstant::twoPi);
+    tensor Q = I;
 
-    scalar theta(rndGen_.scalar01()*mathematicalConstant::twoPi);
+    if (!cP.pointMolecule())
+    {
+        pi = equipartitionAngularMomentum(temperature, cP);
 
-    scalar psi(rndGen_.scalar01()*mathematicalConstant::twoPi);
+        scalar phi(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)
-    );
+        scalar theta(rndGen_.scalar01()*mathematicalConstant::twoPi);
+
+        scalar psi(rndGen_.scalar01()*mathematicalConstant::twoPi);
+
+        Q = tensor
+        (
+            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
     (

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

@@ -140,7 +140,7 @@ private:
         inline vector equipartitionAngularMomentum
         (
             scalar temperature,
-            const diagTensor& momentOfInertia
+            const molecule::constantProperties& cP
         );
 
         //- Disallow default bitwise copy construct

src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H

@@ -534,15 +534,29 @@ inline Foam::vector Foam::moleculeCloud::equipartitionLinearVelocity
 inline Foam::vector Foam::moleculeCloud::equipartitionAngularMomentum
 (
     scalar temperature,
-    const diagTensor& momentOfInertia
+    const molecule::constantProperties& cP
 )
 {
-    return vector
-    (
-        sqrt(kb*temperature*momentOfInertia.xx())*rndGen_.GaussNormal(),
-        sqrt(kb*temperature*momentOfInertia.yy())*rndGen_.GaussNormal(),
-        sqrt(kb*temperature*momentOfInertia.zz())*rndGen_.GaussNormal()
-    );
+    scalar sqrtKbT = sqrt(kb*temperature);
+
+    if (cP.linearMolecule())
+    {
+        return sqrtKbT*vector
+        (
+            0.0,
+            sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal(),
+            sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal()
+        );
+    }
+    else
+    {
+        return sqrtKbT*vector
+        (
+            sqrt(cP.momentOfInertia().xx())*rndGen_.GaussNormal(),
+            sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal(),
+            sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal()
+        );
+    }
 }
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //