ENH: Added monoatomic object and updated applications.

applications/solvers/discreteMethods/molecularDynamics/mdFoam/Make/options

@@ -1,4 +1,7 @@
+EXE_DEBUG = -DFULLDEBUG -g -O0
+
 EXE_INC = \
+    ${EXE_DEBUG} \
     -I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \
     -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \

applications/solvers/discreteMethods/molecularDynamics/mdFoam/createFields.H

@@ -12,16 +12,16 @@
         mesh
     );
 
-    word cloudName("polyatomicCloud");
+    word polyatomicCloudName("polyatomicCloud");
 
-    potential pot
+    potential polyPot
     (
         mesh,
         IOdictionary
         (
             IOobject
             (
-                cloudName +  "Properties",
+                polyatomicCloudName +  "Properties",
                 mesh.time().constant(),
                 mesh,
                 IOobject::MUST_READ_IF_MODIFIED,
@@ -31,4 +31,25 @@
         )
     );
 
-    polyatomicCloud molecules(cloudName, mesh, pot);
+    polyatomicCloud polyatomics(polyatomicCloudName, mesh, polyPot);
+
+    word monoatomicCloudName("monoatomicCloud");
+
+    potential monoPot
+    (
+        mesh,
+        IOdictionary
+        (
+            IOobject
+            (
+                monoatomicCloudName +  "Properties",
+                mesh.time().constant(),
+                mesh,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE,
+                false
+            )
+        )
+    );
+
+    monoatomicCloud monoatomics(monoatomicCloudName, mesh, monoPot);

applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C

@@ -30,6 +30,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "monoatomicCloud.H"
 #include "polyatomicCloud.H"
 
 int main(int argc, char *argv[])
@@ -47,11 +48,13 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << endl;
 
-        molecules.evolve();
+        monoatomics.evolve();
+
+        polyatomics.evolve();
 
         runTime.write();
 
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+        Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
     }

applications/utilities/preProcessing/mdInitialise/mdInitialise.C

@@ -28,6 +28,7 @@ Description
 
 #include "fvCFD.H"
 #include "polyatomicCloud.H"
+#include "monoatomicCloud.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,11 +51,15 @@ int main(int argc, char *argv[])
         )
     );
 
-    IOdictionary idListDict
+    word polyCloudName("polyatomicCloud");
+
+    const dictionary& polyDict(mdInitialiseDict.subDict(polyCloudName));
+
+    IOdictionary polyIdListDict
     (
         IOobject
         (
-            "idList",
+            polyCloudName + "_idList",
             mesh.time().constant(),
             mesh,
             IOobject::NO_READ,
@@ -62,17 +67,15 @@ int main(int argc, char *argv[])
         )
     );
 
-    word cloudName("polyatomicCloud");
+    potential polyPot
-
-    potential pot
     (
         mesh,
-        mdInitialiseDict,
+        polyDict,
         IOdictionary
         (
             IOobject
             (
-                cloudName +  "Properties",
+                polyCloudName +  "Properties",
                 mesh.time().constant(),
                 mesh,
                 IOobject::MUST_READ_IF_MODIFIED,
@@ -80,27 +83,72 @@ int main(int argc, char *argv[])
                 false
             )
         ),
-        idListDict
+        polyIdListDict
     );
 
-    polyatomicCloud molecules(cloudName, mesh, pot, mdInitialiseDict);
+    polyatomicCloud poly
+    (
+        polyCloudName,
+        mesh,
+        polyPot,
+        polyDict
+    );
 
-    label totalMolecules = molecules.size();
+    Info<< nl << returnReduce(poly.size(), sumOp<label>()) << " added to "
-
+        << poly.name()
-    if (Pstream::parRun())
-    {
-        reduce(totalMolecules, sumOp<label>());
-    }
-
-    Info<< nl << "Total number of molecules added: " << totalMolecules
         << nl << endl;
 
-    IOstream::defaultPrecision(15);
+    word monoCloudName("monoatomicCloud");
+
+    const dictionary& monoDict(mdInitialiseDict.subDict(monoCloudName));
+
+    IOdictionary monoIdListDict
+    (
+        IOobject
+        (
+            monoCloudName + "_idList",
+            mesh.time().constant(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        )
+    );
+
+    potential monoPot
+    (
+        mesh,
+        monoDict,
+        IOdictionary
+        (
+            IOobject
+            (
+                monoCloudName +  "Properties",
+                mesh.time().constant(),
+                mesh,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE,
+                false
+            )
+        ),
+        monoIdListDict
+    );
+
+    monoatomicCloud mono
+    (
+        monoCloudName,
+        mesh,
+        monoPot,
+        monoDict
+    );
+
+    Info<< nl << returnReduce(mono.size(), sumOp<label>()) << " added to "
+        << mono.name()
+        << nl << endl;
 
     if (!mesh.write())
     {
         FatalErrorIn(args.executable())
-            << "Failed writing polyatomicCloud."
+            << "Failed writing."
             << nl << exit(FatalError);
     }
 

src/lagrangian/molecularDynamics/molecule/Make/files

@@ -2,6 +2,8 @@ clouds/baseClasses/moleculeCloud/moleculeCloud.C
 
 molecules/constPropSite/constPropSite.C
 
+molecules/monoatomic/monoatomic.C
+molecules/monoatomic/monoatomicIO.C
 molecules/polyatomic/polyatomic.C
 molecules/polyatomic/polyatomicIO.C
 

src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.C

@@ -478,7 +478,7 @@ void Foam::MoleculeCloud<MoleculeType>::removeHighEnergyOverlaps()
 template<class MoleculeType>
 void Foam::MoleculeCloud<MoleculeType>::initialiseMolecules
 (
-    const IOdictionary& mdInitialiseDict
+    const dictionary& mdInitialiseDict
 )
 {
     Info<< nl
@@ -607,7 +607,6 @@ void Foam::MoleculeCloud<MoleculeType>::initialiseMolecules
                         continue;
                     }
 
-
                     latticeCellScale = pow
                     (
                         unitCellMass/(det(latticeCellShape)*massDensity),
@@ -1049,37 +1048,11 @@ void Foam::MoleculeCloud<MoleculeType>::createMolecule
 
     const typename MoleculeType::constantProperties& cP(constProps(id));
 
-    vector v = equipartitionLinearVelocity(temperature, cP.mass());
+    typename MoleculeType::trackingData td
-
-    v += bulkVelocity;
-
-    vector pi = vector::zero;
-
-    tensor Q = I;
-
-    if (!cP.pointMolecule())
-    {
-        pi = equipartitionAngularMomentum(temperature, cP);
-
-        scalar phi(rndGen_.scalar01()*twoPi);
-
-        scalar theta(rndGen_.scalar01()*twoPi);
-
-        scalar psi(rndGen_.scalar01()*twoPi);
-
-        Q = tensor
     (
-            cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi),
+        *this,
-            cos(psi)*sin(phi) + cos(theta)*cos(phi)*sin(psi),
+        MoleculeType::trackingData::tpAccess
-            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
     (
@@ -1090,13 +1063,11 @@ void Foam::MoleculeCloud<MoleculeType>::createMolecule
             cell,
             tetFace,
             tetPt,
-            Q,
+            temperature,
-            v,
+            bulkVelocity,
-            vector::zero,
-            pi,
-            vector::zero,
             specialPosition,
-            constProps(id),
+            cP,
+            td,
             special,
             id
         )
@@ -1136,7 +1107,7 @@ Foam::MoleculeCloud<MoleculeType>::MoleculeCloud
     cellOccupancy_(mesh_.nCells()),
     il_(mesh_, pot_.pairPotentials().rCutMax(), false),
     constPropList_(),
-    rndGen_(clock::getTime())
+    rndGen_(label(971501) + 1526*Pstream::myProcNo())
 {
     if (readFields)
     {
@@ -1159,7 +1130,7 @@ Foam::MoleculeCloud<MoleculeType>::MoleculeCloud
     const word& cloudName,
     const polyMesh& mesh,
     const potential& pot,
-    const IOdictionary& mdInitialiseDict,
+    const dictionary& mdInitialiseDict,
     bool readFields
 )
 :
@@ -1167,9 +1138,9 @@ Foam::MoleculeCloud<MoleculeType>::MoleculeCloud
     moleculeCloud(),
     mesh_(mesh),
     pot_(pot),
-    il_(mesh_, 0.0, false),
+    il_(mesh_),
     constPropList_(),
-    rndGen_(clock::getTime())
+    rndGen_(label(971501) + 1526*Pstream::myProcNo())
 {
     if (readFields)
     {
@@ -1189,18 +1160,34 @@ Foam::MoleculeCloud<MoleculeType>::MoleculeCloud
 template<class MoleculeType>
 void Foam::MoleculeCloud<MoleculeType>::evolve()
 {
-    typename MoleculeType::trackingData td0(*this, 0);
+    typename MoleculeType::trackingData td0
+    (
+        *this,
+        MoleculeType::trackingData::tpFirstVelocityHalfStep
+    );
     Cloud<MoleculeType>::move(td0, mesh_.time().deltaTValue());
 
-    typename MoleculeType::trackingData td1(*this, 1);
+    typename MoleculeType::trackingData td1
+    (
+        *this,
+        MoleculeType::trackingData::tpLinearTrack
+    );
     Cloud<MoleculeType>::move(td1, mesh_.time().deltaTValue());
 
-    typename MoleculeType::trackingData td2(*this, 2);
+    typename MoleculeType::trackingData td2
+    (
+        *this,
+        MoleculeType::trackingData::tpRotationalTrack
+    );
     Cloud<MoleculeType>::move(td2, mesh_.time().deltaTValue());
 
     calculateForce();
 
-    typename MoleculeType::trackingData td3(*this, 3);
+    typename MoleculeType::trackingData td3
+    (
+        *this,
+        MoleculeType::trackingData::tpSecondVelocityHalfStep
+    );
     Cloud<MoleculeType>::move(td3, mesh_.time().deltaTValue());
 
     info();
@@ -1231,144 +1218,18 @@ void Foam::MoleculeCloud<MoleculeType>::calculateForce()
 
 
 template<class MoleculeType>
-void Foam::MoleculeCloud<MoleculeType>::info() const
+void Foam::MoleculeCloud<MoleculeType>::info()
 {
     // Calculates and prints the mean momentum and energy in the system
     // and the number of molecules.
 
-    vector totalLinearMomentum(vector::zero);
+    typename MoleculeType::trackingData td
-
-    vector totalAngularMomentum(vector::zero);
-
-    scalar maxVelocityMag = 0.0;
-
-    scalar totalMass = 0.0;
-
-    scalar totalLinearKE = 0.0;
-
-    scalar totalAngularKE = 0.0;
-
-    scalar totalPE = 0.0;
-
-    scalar totalrDotf = 0.0;
-
-    //vector CentreOfMass(vector::zero);
-
-    label nMols = this->size();
-
-    label dofs = 0;
-
-    {
-        forAllConstIter(typename MoleculeCloud<MoleculeType>, *this, mol)
-        {
-            const label molId = mol().id();
-
-            scalar molMass(this->constProps(molId).mass());
-
-            totalMass += molMass;
-
-            //CentreOfMass += mol().position()*molMass;
-        }
-
-        // if (nMols)
-        // {
-        //     CentreOfMass /= totalMass;
-        // }
-
-        forAllConstIter(typename MoleculeCloud<MoleculeType>, *this, mol)
-        {
-            const label molId = mol().id();
-
-            const typename MoleculeType::constantProperties cP
     (
-                this->constProps(molId)
+        *this,
+        MoleculeType::trackingData::tpAccess
     );
 
-            scalar molMass(cP.mass());
+    MoleculeType::info(td);
-
-            const diagTensor& molMoI(cP.momentOfInertia());
-
-            const vector& molV(mol().v());
-
-            const vector& molOmega(inv(molMoI) & mol().pi());
-
-            vector molPiGlobal = mol().Q() & mol().pi();
-
-            totalLinearMomentum += molV * molMass;
-
-            totalAngularMomentum += molPiGlobal;
-            //+((mol().position() - CentreOfMass) ^ (molV * molMass));
-
-            if (mag(molV) > maxVelocityMag)
-            {
-                maxVelocityMag = mag(molV);
-            }
-
-            totalLinearKE += 0.5*molMass*magSqr(molV);
-
-            totalAngularKE += 0.5*(molOmega & molMoI & molOmega);
-
-            totalPE += mol().potentialEnergy();
-
-            totalrDotf += tr(mol().rf());
-
-            dofs += cP.degreesOfFreedom();
-        }
-    }
-
-    scalar meshVolume = sum(mesh_.cellVolumes());
-
-    if (Pstream::parRun())
-    {
-        reduce(totalLinearMomentum, sumOp<vector>());
-        reduce(totalAngularMomentum, sumOp<vector>());
-        reduce(maxVelocityMag, maxOp<scalar>());
-        reduce(totalMass, sumOp<scalar>());
-        reduce(totalLinearKE, sumOp<scalar>());
-        reduce(totalAngularKE, sumOp<scalar>());
-        reduce(totalPE, sumOp<scalar>());
-        reduce(totalrDotf, sumOp<scalar>());
-        reduce(nMols, sumOp<label>());
-        reduce(dofs, sumOp<label>());
-        reduce(meshVolume, sumOp<scalar>());
-    }
-
-    if (nMols)
-    {
-        Info<< nl << "Number of molecules in " << this->name() << " = "
-            << nMols << nl
-            << "    Overall number density = "
-            << nMols/meshVolume << nl
-            << "    Overall mass density = "
-            << totalMass/meshVolume << nl
-            << "    Average linear momentum per molecule = "
-            << totalLinearMomentum/nMols << ' '
-            << mag(totalLinearMomentum)/nMols << nl
-            << "    Average angular momentum per molecule = "
-            << totalAngularMomentum << ' '
-            << mag(totalAngularMomentum)/nMols << nl
-            << "    maximum |velocity| = "
-            << maxVelocityMag << nl
-            << "    Average linear KE per molecule = "
-            << totalLinearKE/nMols << nl
-            << "    Average angular KE per molecule = "
-            << totalAngularKE/nMols << nl
-            << "    Average PE per molecule = "
-            << totalPE/nMols << nl
-            << "    Average TE per molecule = "
-            <<
-            (
-                  totalLinearKE
-                + totalAngularKE
-                + totalPE
-            )
-            /nMols
-            << nl << endl;
-    }
-    else
-    {
-        Info<< nl << "No molecules in " << this->name() << nl << endl;
-    }
 }
 
 

src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.H

@@ -122,10 +122,7 @@ private:
         void removeHighEnergyOverlaps();
 
         //-
-        void initialiseMolecules
+        void initialiseMolecules(const dictionary& mdInitialiseDict);
-        (
-            const IOdictionary& mdInitialiseDict
-        );
 
         //-
         void createMolecule
@@ -143,20 +140,6 @@ private:
         //-
         label nSites() const;
 
-        //-
-        inline vector equipartitionLinearVelocity
-        (
-            scalar temperature,
-            scalar mass
-        );
-
-        //-
-        inline vector equipartitionAngularMomentum
-        (
-            scalar temperature,
-            const typename MoleculeType::constantProperties& cP
-        );
-
         //- Disallow default bitwise copy construct
         MoleculeCloud(const MoleculeCloud&);
 
@@ -183,7 +166,7 @@ public:
             const word& cloudName,
             const polyMesh& mesh,
             const potential& pot,
-            const IOdictionary& mdInitialiseDict,
+            const dictionary& mdInitialiseDict,
             bool readFields = true
         );
 
@@ -197,7 +180,7 @@ public:
         void calculateForce();
 
         //- Print cloud information
-        void info() const;
+        void info();
 
 
         // Access
@@ -226,6 +209,20 @@ public:
             //-
             inline Random& rndGen();
 
+            //-
+            inline vector equipartitionLinearVelocity
+            (
+                scalar temperature,
+                scalar mass
+            );
+
+            //-
+            inline vector equipartitionAngularMomentum
+            (
+                scalar temperature,
+                const typename MoleculeType::constantProperties& cP
+            );
+
 
     // Member Operators
 

src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloudI.H

@@ -30,7 +30,7 @@ using namespace Foam::constant;
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class MoleculeType>
-void Foam::MoleculeCloud<MoleculeType>::evaluatePair
+inline void Foam::MoleculeCloud<MoleculeType>::evaluatePair
 (
     MoleculeType& molI,
     MoleculeType& molJ
@@ -157,7 +157,7 @@ void Foam::MoleculeCloud<MoleculeType>::evaluatePair
 
 
 template<class MoleculeType>
-bool Foam::MoleculeCloud<MoleculeType>::evaluatePotentialLimit
+inline bool Foam::MoleculeCloud<MoleculeType>::evaluatePotentialLimit
 (
     MoleculeType& molI,
     MoleculeType& molJ
@@ -310,8 +310,63 @@ bool Foam::MoleculeCloud<MoleculeType>::evaluatePotentialLimit
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
 template<class MoleculeType>
-Foam::vector
+inline const Foam::polyMesh& Foam::MoleculeCloud<MoleculeType>::mesh() const
+{
+    return mesh_;
+}
+
+
+template<class MoleculeType>
+inline const Foam::potential& Foam::MoleculeCloud<MoleculeType>::pot() const
+{
+    return pot_;
+}
+
+
+template<class MoleculeType>
+inline const Foam::List<Foam::DynamicList<MoleculeType*> >&
+Foam::MoleculeCloud<MoleculeType>::cellOccupancy() const
+{
+    return cellOccupancy_;
+}
+
+
+template<class MoleculeType>
+inline const Foam::InteractionLists<MoleculeType>&
+Foam::MoleculeCloud<MoleculeType>::il() const
+{
+    return il_;
+}
+
+
+template<class MoleculeType>
+inline const Foam::List<typename MoleculeType::constantProperties>
+Foam::MoleculeCloud<MoleculeType>::constProps() const
+{
+    return constPropList_;
+}
+
+
+template<class MoleculeType>
+inline const typename MoleculeType::constantProperties&
+Foam::MoleculeCloud<MoleculeType>::constProps(label id) const
+{
+    return constPropList_[id];
+}
+
+
+template<class MoleculeType>
+inline Foam::Random& Foam::MoleculeCloud<MoleculeType>::rndGen()
+{
+    return rndGen_;
+}
+
+
+template<class MoleculeType>
+inline Foam::vector
 Foam::MoleculeCloud<MoleculeType>::equipartitionLinearVelocity
 (
     scalar temperature,
@@ -328,7 +383,7 @@ Foam::MoleculeCloud<MoleculeType>::equipartitionLinearVelocity
 
 
 template<class MoleculeType>
-Foam::vector
+inline Foam::vector
 Foam::MoleculeCloud<MoleculeType>::equipartitionAngularMomentum
 (
     scalar temperature,
@@ -358,59 +413,4 @@ Foam::MoleculeCloud<MoleculeType>::equipartitionAngularMomentum
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class MoleculeType>
-const Foam::polyMesh& Foam::MoleculeCloud<MoleculeType>::mesh() const
-{
-    return mesh_;
-}
-
-
-template<class MoleculeType>
-const Foam::potential& Foam::MoleculeCloud<MoleculeType>::pot() const
-{
-    return pot_;
-}
-
-
-template<class MoleculeType>
-const Foam::List<Foam::DynamicList<MoleculeType*> >&
-Foam::MoleculeCloud<MoleculeType>::cellOccupancy() const
-{
-    return cellOccupancy_;
-}
-
-
-template<class MoleculeType>
-const Foam::InteractionLists<MoleculeType>&
-Foam::MoleculeCloud<MoleculeType>::il() const
-{
-    return il_;
-}
-
-
-template<class MoleculeType>
-const Foam::List<typename MoleculeType::constantProperties>
-Foam::MoleculeCloud<MoleculeType>::constProps() const
-{
-    return constPropList_;
-}
-
-
-template<class MoleculeType>
-const typename MoleculeType::constantProperties&
-Foam::MoleculeCloud<MoleculeType>::constProps(label id) const
-{
-    return constPropList_[id];
-}
-
-
-template<class MoleculeType>
-Foam::Random& Foam::MoleculeCloud<MoleculeType>::rndGen()
-{
-    return rndGen_;
-}
-
-
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/clouds/derived/monoatomicCloud/monoatomicCloud.H

@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::monoatomicCloud
+
+Description
+    Cloud class to simulate monoatomic molecules
+
+SourceFiles
+    monoatomicCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef monoatomicCloud_H
+#define monoatomicCloud_H
+
+#include "MoleculeCloud.H"
+#include "monoatomic.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef MoleculeCloud<monoatomic> monoatomicCloud;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomic.C

@@ -0,0 +1,199 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*----------------------------------------------------------------------------*/
+
+#include "monoatomic.H"
+#include "Random.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(monoatomic, 0);
+    defineTemplateTypeNameAndDebug(Cloud<monoatomic>, 0);
+}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::monoatomic::move
+(
+    monoatomic::trackingData& td,
+    const scalar trackTime
+)
+{
+    td.switchProcessor = false;
+    td.keepParticle = true;
+
+    if (special_ != SPECIAL_FROZEN)
+    {
+        return td.keepParticle;
+    }
+
+    const constantProperties& constProps(td.cloud().constProps(id_));
+
+    if (td.part() == trackingData::tpFirstVelocityHalfStep)
+    {
+        // First leapfrog velocity adjust part, required before tracking+force
+        // part
+
+        v_ += 0.5*trackTime*a_;
+    }
+    else if (td.part() == trackingData::tpLinearTrack)
+    {
+        // Leapfrog tracking part
+
+        scalar tEnd = (1.0 - stepFraction())*trackTime;
+        scalar dtMax = tEnd;
+
+        while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
+        {
+            // set the lagrangian time-step
+            scalar dt = min(dtMax, tEnd);
+
+            dt *= trackToFace(position() + dt*v_, td);
+
+            tEnd -= dt;
+            stepFraction() = 1.0 - tEnd/trackTime;
+        }
+
+        setSitePositions(constProps);
+    }
+    else if (td.part() == trackingData::tpSecondVelocityHalfStep)
+    {
+        // Second leapfrog velocity adjust part, required after tracking+force
+        // part
+
+        a_ = siteForces_[0]/constProps.mass();
+
+        v_ += 0.5*trackTime*a_;
+    }
+    else if (td.part() != trackingData::tpRotationalTrack)
+    {
+        FatalErrorIn("monoatomic::move(trackingData&, const scalar)") << nl
+            << td.part() << " is an invalid part of the integration method."
+            << abort(FatalError);
+    }
+
+    return td.keepParticle;
+}
+
+
+void Foam::monoatomic::transformProperties(const tensor& T)
+{
+    particle::transformProperties(T);
+
+    v_ = transform(T, v_);
+
+    a_ = transform(T, a_);
+
+    rf_ = transform(T, rf_);
+
+    sitePositions_[0] = position_ + (T & (sitePositions_[0] - position_));
+
+    siteForces_[0] = T & siteForces_[0];
+}
+
+
+void Foam::monoatomic::transformProperties(const vector& separation)
+{
+    particle::transformProperties(separation);
+
+    if (special_ == SPECIAL_TETHERED)
+    {
+        specialPosition_ += separation;
+    }
+
+    sitePositions_[0] += separation;
+}
+
+
+void Foam::monoatomic::setSitePositions(const constantProperties& constProps)
+{
+    sitePositions_[0] = position_;
+}
+
+
+void Foam::monoatomic::setSiteSizes(label size)
+{
+    // Nothing required, size controlled internally
+}
+
+
+bool Foam::monoatomic::hitPatch
+(
+    const polyPatch&,
+    trackingData&,
+    const label,
+    const scalar,
+    const tetIndices&
+)
+{
+    return false;
+}
+
+
+void Foam::monoatomic::hitProcessorPatch
+(
+    const processorPolyPatch&,
+    trackingData& td
+)
+{
+    td.switchProcessor = true;
+}
+
+
+void Foam::monoatomic::hitWallPatch
+(
+    const wallPolyPatch& wpp,
+    trackingData& td,
+    const tetIndices& tetIs
+)
+{
+    // Use of the normal from tetIs is not required as
+    // hasWallImpactDistance is false.
+    vector nw = normal();
+    nw /= mag(nw);
+
+    scalar vn = v_ & nw;
+
+    // Specular reflection
+    if (vn > 0)
+    {
+        v_ -= 2*vn*nw;
+    }
+}
+
+
+void Foam::monoatomic::hitPatch
+(
+    const polyPatch&,
+    trackingData& td
+)
+{
+    td.keepParticle = false;
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomic.H

@@ -0,0 +1,418 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::monoatomic
+
+Description
+    Foam::monoatomic
+
+SourceFiles
+    monoatomicI.H
+    monoatomic.C
+    monoatomicIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef monoatomic_H
+#define monoatomic_H
+
+#include "particle.H"
+#include "IOstream.H"
+#include "autoPtr.H"
+#include "diagTensor.H"
+#include "constPropSite.H"
+#include "MoleculeCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class monoatomic Declaration
+\*---------------------------------------------------------------------------*/
+
+class monoatomic
+:
+    public particle
+{
+
+public:
+
+    // Values of special that are less than zero are for built-in functionality.
+    // Values greater than zero are user specifiable/expandable (i.e. test
+    // special_ >= SPECIAL_USER)
+
+    enum specialTypes
+    {
+        SPECIAL_TETHERED = -1,
+        SPECIAL_FROZEN   = -2,
+        NOT_SPECIAL      = 0,
+        SPECIAL_USER     = 1
+    };
+
+    //- Class to hold monoatomic constant properties
+    class constantProperties
+    {
+        // Private data
+
+            //- Sites of mass, charge or interaction
+            FixedList<constPropSite, 1> sites_;
+
+            //- Which sites require electrostatic interactions
+            FixedList<label, 1> electrostaticSites_;
+
+            //- Which sites require pair interactions
+            FixedList<label, 1> pairPotSites_;
+
+            //-
+            scalar mass_;
+
+
+    public:
+
+        //-
+        inline constantProperties();
+
+        //- Construct from dictionary
+        inline constantProperties
+        (
+            const dictionary& dict,
+            const List<label>& siteIds
+        );
+
+        // Member functions
+
+        //-
+        inline const FixedList<constPropSite, 1>& sites() const;
+
+        //-
+        inline const FixedList<label, 1>& pairPotSites() const;
+
+        //-
+        inline const FixedList<label, 1>& electrostaticSites() const;
+
+        //-
+        inline label degreesOfFreedom() const;
+
+        //-
+        inline scalar mass() const;
+
+        //-
+        inline label nSites() const;
+    };
+
+
+    //- Class used to pass tracking data to the trackToFace function
+    class trackingData
+    :
+        public particle::TrackingData<MoleculeCloud<monoatomic> >
+    {
+    public:
+
+        enum trackPart
+        {
+            tpFirstVelocityHalfStep,
+            tpLinearTrack,
+            tpRotationalTrack,
+            tpSecondVelocityHalfStep,
+            tpAccess
+        };
+
+    private:
+
+        // label specifying which part of the integration algorithm is taking
+        label part_;
+
+
+    public:
+
+        // Constructors
+
+            trackingData(MoleculeCloud<monoatomic>& cloud, trackPart part)
+            :
+                particle::TrackingData<MoleculeCloud<monoatomic> >(cloud),
+                part_(part)
+            {}
+
+        // Member functions
+
+            inline label part() const
+            {
+                return part_;
+            }
+    };
+
+
+private:
+
+    // Private data
+
+          //- Linear velocity of monoatomic
+        vector v_;
+
+        //- Total linear acceleration of monoatomic
+        vector a_;
+
+            //-
+        vector specialPosition_;
+
+        //-
+        scalar potentialEnergy_;
+
+        // - r_ij f_ij, stress dyad
+        tensor rf_;
+
+//         // - r_ij outer product f_ij: virial contribution
+//         tensor rDotf_;
+
+        //-
+        label special_;
+
+        //-
+        label id_;
+
+        //-
+        FixedList<vector, 1> siteForces_;
+
+        //-
+        FixedList<vector, 1> sitePositions_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("monoatomic");
+
+    friend class Cloud<monoatomic>;
+
+    // Constructors
+
+        //- Construct with macroscopic description
+        inline monoatomic
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label cellI,
+            const label tetFaceI,
+            const label tetPtI,
+            const scalar temperature,
+            const vector& bulkVelocity,
+            const vector& specialPosition,
+            const constantProperties& constProps,
+            trackingData& td,
+            const label special,
+            const label id
+        );
+
+        //- Construct from all components
+        inline monoatomic
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label cellI,
+            const label tetFaceI,
+            const label tetPtI,
+            const vector& v,
+            const vector& a,
+            const vector& specialPosition,
+            const constantProperties& constProps,
+            const label special,
+            const label id
+        );
+
+        //- Construct from Istream
+        monoatomic
+        (
+            const polyMesh& mesh,
+            Istream& is,
+            bool readFields = true
+        );
+
+        //- Construct and return a clone
+        autoPtr<particle> clone() const
+        {
+            return autoPtr<particle>(new monoatomic(*this));
+        }
+
+        //- Factory class to read-construct particles used for
+        //  parallel transfer
+        class iNew
+        {
+            const polyMesh& mesh_;
+
+        public:
+
+            iNew(const polyMesh& mesh)
+            :
+                mesh_(mesh)
+            {}
+
+            autoPtr<monoatomic> operator()(Istream& is) const
+            {
+                return autoPtr<monoatomic>(new monoatomic(mesh_, is, true));
+            }
+        };
+
+
+    // Member Functions
+
+        // Tracking
+
+            //-
+            bool move(trackingData&, const scalar trackTime);
+
+            //-
+            virtual void transformProperties(const tensor& T);
+
+            //-
+            virtual void transformProperties(const vector& separation);
+
+            //-
+            void setSitePositions(const constantProperties& constProps);
+
+            //-
+            void setSiteSizes(label size);
+
+
+        // Access
+
+            //-
+            inline const vector& v() const;
+
+            //-
+            inline vector& v();
+
+            //-
+            inline const vector& a() const;
+
+            //-
+            inline vector& a();
+
+            //-
+            inline const FixedList<vector, 1>& siteForces() const;
+
+            //-
+            inline FixedList<vector, 1>& siteForces();
+
+            //-
+            inline const FixedList<vector, 1>& sitePositions() const;
+
+            //-
+            inline FixedList<vector, 1>& sitePositions();
+
+            //-
+            inline const vector& specialPosition() const;
+
+            //-
+            inline vector& specialPosition();
+
+            //-
+            inline scalar potentialEnergy() const;
+
+            //-
+            inline scalar& potentialEnergy();
+
+            //-
+            inline const tensor& rf() const;
+
+            //-
+            inline tensor& rf();
+
+            //-
+            inline label special() const;
+
+            //-
+            inline bool tethered() const;
+
+            //-
+            inline label id() const;
+
+
+    // Member Operators
+
+        //- Overridable function to handle the particle hitting a patch
+        //  Executed before other patch-hitting functions
+        bool hitPatch
+        (
+            const polyPatch&,
+            trackingData& td,
+            const label patchI,
+            const scalar trackFraction,
+            const tetIndices& tetIs
+        );
+
+        //- Overridable function to handle the particle hitting a processorPatch
+        void hitProcessorPatch
+        (
+            const processorPolyPatch&,
+            trackingData& td
+        );
+
+        //- Overridable function to handle the particle hitting a wallPatch
+        void hitWallPatch
+        (
+            const wallPolyPatch&,
+            trackingData& td,
+            const tetIndices&
+        );
+
+        //- Overridable function to handle the particle hitting a polyPatch
+        void hitPatch
+        (
+            const polyPatch&,
+            trackingData& td
+        );
+
+
+    // I-O
+
+        //- Read
+        static void readFields(Cloud<monoatomic>& mC);
+
+        //- Write
+        static void writeFields(const Cloud<monoatomic>& mC);
+
+        //- Show info
+        static void info(trackingData& td);
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<<(Ostream&, const monoatomic&);
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "monoatomicI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicI.H

@@ -0,0 +1,328 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+inline Foam::monoatomic::constantProperties::constantProperties()
+:
+    sites_(),
+    electrostaticSites_(-1),
+    pairPotSites_(-1),
+    mass_(0)
+{}
+
+
+inline Foam::monoatomic::constantProperties::constantProperties
+(
+    const dictionary& dict,
+    const List<label>& siteIds
+)
+:
+    sites_(),
+    electrostaticSites_(-1),
+    pairPotSites_(-1),
+    mass_(0)
+{
+    if (siteIds.size() != 1)
+    {
+        FatalErrorIn
+        (
+            "Foam::monoatomic::constantProperties::constantProperties"
+            "("
+                "const dictionary& dict, "
+                "const List<label>& siteIds"
+            ")"
+        )
+            << "monoatomic, single site only, given: " << dict
+            << nl << abort(FatalError);
+    }
+
+    FixedList<word, 1> siteIdNames = dict.lookup("siteIds");
+
+    FixedList<vector, 1> siteReferencePositions
+    (
+        dict.lookup("siteReferencePositions")
+    );
+
+    FixedList<scalar, 1> siteMasses(dict.lookup("siteMasses"));
+
+    FixedList<scalar, 1> siteCharges(dict.lookup("siteCharges"));
+
+    FixedList<word, 1> pairPotentialIds(dict.lookup("pairPotentialSiteIds"));
+
+    constPropSite site = sites_[0];
+
+    site = constPropSite
+    (
+        siteReferencePositions[0],
+        siteMasses[0],
+        siteCharges[0],
+        siteIds[0],
+        siteIdNames[0],
+        (findIndex(pairPotentialIds, siteIdNames[0]) != -1),    // pair
+        (mag(siteCharges[0]) > VSMALL)                          // charge
+    );
+
+    mass_ = site.siteMass();
+
+    if (site.pairPotentialSite())
+    {
+        pairPotSites_[0] = 0;
+    }
+
+    if (site.electrostaticSite())
+    {
+        electrostaticSites_[0] = 0;
+    }
+
+    if (!site.pairPotentialSite() && !site.electrostaticSite())
+    {
+        WarningIn
+        (
+            "Foam::monoatomic::constantProperties::constantProperties"
+            "("
+            "const dictionary& dict"
+            ")"
+        )
+        << siteIdNames[0] << " is a non-interacting site." << endl;
+    }
+
+    // Single site monoatomic - no rotational motion.
+    site.siteReferencePosition() = vector::zero;
+}
+
+
+inline Foam::monoatomic::monoatomic
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label cellI,
+    const label tetFaceI,
+    const label tetPtI,
+    const scalar temperature,
+    const vector& bulkVelocity,
+    const vector& specialPosition,
+    const constantProperties& cP,
+    monoatomic::trackingData& td,
+    const label special,
+    const label id
+)
+:
+    particle(mesh, position, cellI, tetFaceI, tetPtI),
+    v_(bulkVelocity),
+    a_(vector::zero),
+    specialPosition_(specialPosition),
+    potentialEnergy_(0.0),
+    rf_(tensor::zero),
+    special_(special),
+    id_(id),
+    siteForces_(vector::zero),
+    sitePositions_()
+{
+    setSitePositions(cP);
+
+    v_ += td.cloud().equipartitionLinearVelocity(temperature, cP.mass());
+}
+
+
+inline Foam::monoatomic::monoatomic
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label cellI,
+    const label tetFaceI,
+    const label tetPtI,
+    const vector& v,
+    const vector& a,
+    const vector& specialPosition,
+    const constantProperties& cP,
+    const label special,
+    const label id
+
+)
+:
+    particle(mesh, position, cellI, tetFaceI, tetPtI),
+    v_(v),
+    a_(a),
+    specialPosition_(specialPosition),
+    potentialEnergy_(0.0),
+    rf_(tensor::zero),
+    special_(special),
+    id_(id),
+    siteForces_(vector::zero),
+    sitePositions_()
+{
+    setSitePositions(cP);
+}
+
+// * * * * * * * constantProperties Member Functions * * * * * * * * * * * * //
+
+inline const Foam::FixedList<Foam::constPropSite, 1>&
+Foam::monoatomic::constantProperties::sites() const
+{
+    return sites_;
+}
+
+
+inline const Foam::FixedList<Foam::label, 1>&
+Foam::monoatomic::constantProperties::pairPotSites() const
+{
+    return pairPotSites_;
+}
+
+
+inline const Foam::FixedList<Foam::label, 1>&
+Foam::monoatomic::constantProperties::electrostaticSites() const
+{
+    return electrostaticSites_;
+}
+
+
+inline Foam::label
+Foam::monoatomic::constantProperties::degreesOfFreedom() const
+{
+    return 3;
+}
+
+
+inline Foam::scalar Foam::monoatomic::constantProperties::mass() const
+{
+    return mass_;
+}
+
+
+inline Foam::label Foam::monoatomic::constantProperties::nSites() const
+{
+    return 1;
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+inline const Foam::vector& Foam::monoatomic::v() const
+{
+    return v_;
+}
+
+
+inline Foam::vector& Foam::monoatomic::v()
+{
+    return v_;
+}
+
+
+inline const Foam::vector& Foam::monoatomic::a() const
+{
+    return a_;
+}
+
+
+inline Foam::vector& Foam::monoatomic::a()
+{
+    return a_;
+}
+
+
+inline const Foam::FixedList<Foam::vector, 1>&
+Foam::monoatomic::siteForces() const
+{
+    return siteForces_;
+}
+
+
+inline Foam::FixedList<Foam::vector, 1>& Foam::monoatomic::siteForces()
+{
+    return siteForces_;
+}
+
+
+inline const Foam::FixedList<Foam::vector, 1>&
+Foam::monoatomic::sitePositions() const
+{
+    return sitePositions_;
+}
+
+
+inline Foam::FixedList<Foam::vector, 1>& Foam::monoatomic::sitePositions()
+{
+    return sitePositions_;
+}
+
+
+inline const Foam::vector& Foam::monoatomic::specialPosition() const
+{
+    return specialPosition_;
+}
+
+
+inline Foam::vector& Foam::monoatomic::specialPosition()
+{
+    return specialPosition_;
+}
+
+
+inline Foam::scalar Foam::monoatomic::potentialEnergy() const
+{
+    return potentialEnergy_;
+}
+
+
+inline Foam::scalar& Foam::monoatomic::potentialEnergy()
+{
+    return potentialEnergy_;
+}
+
+
+inline const Foam::tensor& Foam::monoatomic::rf() const
+{
+    return rf_;
+}
+
+
+inline Foam::tensor& Foam::monoatomic::rf()
+{
+    return rf_;
+}
+
+
+inline Foam::label Foam::monoatomic::special() const
+{
+    return special_;
+}
+
+
+inline bool Foam::monoatomic::tethered() const
+{
+    return special_ == SPECIAL_TETHERED;
+}
+
+
+inline Foam::label Foam::monoatomic::id() const
+{
+    return id_;
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicIO.C

@@ -0,0 +1,305 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "monoatomic.H"
+#include "IOstreams.H"
+#include "Cloud.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::monoatomic::monoatomic
+(
+    const polyMesh& mesh,
+    Istream& is,
+    bool readFields
+)
+:
+    particle(mesh, is, readFields),
+    v_(vector::zero),
+    a_(vector::zero),
+    specialPosition_(vector::zero),
+    potentialEnergy_(0.0),
+    rf_(tensor::zero),
+    special_(0),
+    id_(0),
+    siteForces_(),
+    sitePositions_()
+{
+    if (readFields)
+    {
+        if (is.format() == IOstream::ASCII)
+        {
+            is  >> v_;
+            is  >> a_;
+            is  >> siteForces_;
+            potentialEnergy_ = readScalar(is);
+            is  >> rf_;
+            special_ = readLabel(is);
+            id_ = readLabel(is);
+            is  >> sitePositions_;
+            is  >> specialPosition_;
+        }
+        else
+        {
+            is.read
+            (
+                reinterpret_cast<char*>(&v_),
+                sizeof(v_)
+              + sizeof(a_)
+              + sizeof(specialPosition_)
+              + sizeof(potentialEnergy_)
+              + sizeof(rf_)
+              + sizeof(special_)
+              + sizeof(id_)
+            );
+
+            is  >> siteForces_ >> sitePositions_;
+        }
+    }
+
+    // Check state of Istream
+    is.check
+    (
+        "Foam::monoatomic::monoatomic"
+        "("
+            "const polyMesh& mesh,"
+            "Istream& is,"
+            "bool readFields"
+        ")"
+    );
+}
+
+
+void Foam::monoatomic::readFields(Cloud<monoatomic>& mC)
+{
+    if (!mC.size())
+    {
+        return;
+    }
+
+    particle::readFields(mC);
+
+    IOField<vector> v(mC.fieldIOobject("v", IOobject::MUST_READ));
+    mC.checkFieldIOobject(mC, v);
+
+    IOField<vector> a(mC.fieldIOobject("a", IOobject::MUST_READ));
+    mC.checkFieldIOobject(mC, a);
+
+    IOField<vector> specialPosition
+    (
+        mC.fieldIOobject("specialPosition", IOobject::MUST_READ)
+    );
+    mC.checkFieldIOobject(mC, specialPosition);
+
+    IOField<label> special(mC.fieldIOobject("special", IOobject::MUST_READ));
+    mC.checkFieldIOobject(mC, special);
+
+    IOField<label> id(mC.fieldIOobject("id", IOobject::MUST_READ));
+    mC.checkFieldIOobject(mC, id);
+
+    label i = 0;
+
+    forAllIter(typename Cloud<monoatomic>, mC, iter)
+    {
+        monoatomic& mol = iter();
+
+        mol.v_ = v[i];
+        mol.a_ = a[i];
+        mol.specialPosition_ = specialPosition[i];
+        mol.special_ = special[i];
+        mol.id_ = id[i];
+        i++;
+    }
+}
+
+
+void Foam::monoatomic::writeFields(const Cloud<monoatomic>& mC)
+{
+    particle::writeFields(mC);
+
+    label np = mC.size();
+
+    IOField<vector> v(mC.fieldIOobject("v", IOobject::NO_READ), np);
+    IOField<vector> a(mC.fieldIOobject("a", IOobject::NO_READ), np);
+    IOField<vector> specialPosition
+    (
+        mC.fieldIOobject("specialPosition", IOobject::NO_READ),
+        np
+    );
+    IOField<label> special(mC.fieldIOobject("special", IOobject::NO_READ), np);
+    IOField<label> id(mC.fieldIOobject("id", IOobject::NO_READ), np);
+
+    label i = 0;
+    forAllConstIter(typename Cloud<monoatomic>, mC, iter)
+    {
+        const monoatomic& mol = iter();
+
+        v[i] = mol.v_;
+        a[i] = mol.a_;
+        specialPosition[i] = mol.specialPosition_;
+        special[i] = mol.special_;
+        id[i] = mol.id_;
+        i++;
+    }
+
+    v.write();
+    a.write();
+    specialPosition.write();
+    special.write();
+    id.write();
+}
+
+
+void Foam::monoatomic::info(monoatomic::trackingData& td)
+{
+    vector totalLinearMomentum(vector::zero);
+    scalar maxVelocityMag = 0.0;
+    scalar totalMass = 0.0;
+    scalar totalLinearKE = 0.0;
+    scalar totalPE = 0.0;
+    scalar totalrDotf = 0.0;
+    label nMols = td.cloud().size();
+
+    forAllConstIter(typename Cloud<monoatomic>, td.cloud(), mol)
+    {
+        const label molId = mol().id();
+        scalar molMass(td.cloud().constProps(molId).mass());
+        totalMass += molMass;
+    }
+
+    forAllConstIter(typename Cloud<monoatomic>, td.cloud(), mol)
+    {
+        const label molId = mol().id();
+        const monoatomic::constantProperties cP
+        (
+            td.cloud().constProps(molId)
+        );
+        scalar molMass(cP.mass());
+        const vector& molV(mol().v());
+        totalLinearMomentum += molV * molMass;
+        if (mag(molV) > maxVelocityMag)
+        {
+            maxVelocityMag = mag(molV);
+        }
+        totalLinearKE += 0.5*molMass*magSqr(molV);
+        totalPE += mol().potentialEnergy();
+        totalrDotf += tr(mol().rf());
+    }
+
+    scalar meshVolume = sum(td.cloud().mesh().cellVolumes());
+
+    if (Pstream::parRun())
+    {
+        reduce(totalLinearMomentum, sumOp<vector>());
+        reduce(maxVelocityMag, maxOp<scalar>());
+        reduce(totalMass, sumOp<scalar>());
+        reduce(totalLinearKE, sumOp<scalar>());
+        reduce(totalPE, sumOp<scalar>());
+        reduce(totalrDotf, sumOp<scalar>());
+        reduce(nMols, sumOp<label>());
+        reduce(meshVolume, sumOp<scalar>());
+    }
+
+    if (nMols)
+    {
+        Info<< nl << "Number of molecules in " << td.cloud().name() << " = "
+            << nMols << nl
+            << "    Overall number density = "
+            << nMols/meshVolume << nl
+            << "    Overall mass density = "
+            << totalMass/meshVolume << nl
+            << "    Average linear momentum per molecule = "
+            << totalLinearMomentum/nMols << ' '
+            << mag(totalLinearMomentum)/nMols << nl
+            << "    maximum |velocity| = "
+            << maxVelocityMag << nl
+            << "    Average linear KE per molecule = "
+            << totalLinearKE/nMols << nl
+            << "    Average angular KE per molecule = "
+            << totalPE/nMols << nl
+            << "    Average TE per molecule = "
+            <<
+            (
+                  totalLinearKE
+                + totalPE
+            )
+            /nMols
+            << nl << endl;
+    }
+    else
+    {
+        Info<< nl << "No molecules in " << td.cloud().name() << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const monoatomic& mol)
+{
+    if (os.format() == IOstream::ASCII)
+    {
+        os  << token::SPACE << static_cast<const particle&>(mol)
+            << token::SPACE << mol.face()
+            << token::SPACE << mol.stepFraction()
+            << token::SPACE << mol.v_
+            << token::SPACE << mol.a_
+            << token::SPACE << mol.specialPosition_
+            << token::SPACE << mol.potentialEnergy_
+            << token::SPACE << mol.rf_
+            << token::SPACE << mol.special_
+            << token::SPACE << mol.id_
+            << token::SPACE << mol.siteForces_
+            << token::SPACE << mol.sitePositions_;
+    }
+    else
+    {
+        os  << static_cast<const particle&>(mol);
+        os.write
+        (
+            reinterpret_cast<const char*>(&mol.v_),
+            sizeof(mol.v_)
+          + sizeof(mol.a_)
+          + sizeof(mol.specialPosition_)
+          + sizeof(mol.potentialEnergy_)
+          + sizeof(mol.rf_)
+          + sizeof(mol.special_)
+          + sizeof(mol.id_)
+        );
+        os  << mol.siteForces_ << mol.sitePositions_;
+    }
+
+    // Check state of Ostream
+    os.check
+    (
+        "Foam::Ostream& Foam::operator<<"
+        "(Foam::Ostream&, const Foam::monoatomic&)"
+    );
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.C

@@ -89,7 +89,7 @@ bool Foam::polyatomic::move
 
     const constantProperties& constProps(td.cloud().constProps(id_));
 
-    if (td.part() == 0)
+    if (td.part() == trackingData::tpFirstVelocityHalfStep)
     {
         // First leapfrog velocity adjust part, required before tracking+force
         // part
@@ -98,7 +98,7 @@ bool Foam::polyatomic::move
 
         pi_ += 0.5*trackTime*tau_;
     }
-    else if (td.part() == 1)
+    else if (td.part() == trackingData::tpLinearTrack)
     {
         // Leapfrog tracking part
 
@@ -116,7 +116,7 @@ bool Foam::polyatomic::move
             stepFraction() = 1.0 - tEnd/trackTime;
         }
     }
-    else if (td.part() == 2)
+    else if (td.part() == trackingData::tpRotationalTrack)
     {
         // Leapfrog orientation adjustment, carried out before force calculation
         // but after tracking stage, i.e. rotation carried once linear motion
@@ -157,7 +157,7 @@ bool Foam::polyatomic::move
 
         setSitePositions(constProps);
     }
-    else if (td.part() == 3)
+    else if (td.part() == trackingData::tpSecondVelocityHalfStep)
     {
         // Second leapfrog velocity adjust part, required after tracking+force
         // part

src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.H

@@ -49,9 +49,6 @@ SourceFiles
 namespace Foam
 {
 
-template<class ParcelType>
-class DsmcParcel;
-
 /*---------------------------------------------------------------------------*\
                            Class polyatomic Declaration
 \*---------------------------------------------------------------------------*/
@@ -150,6 +147,20 @@ public:
     :
         public particle::TrackingData<MoleculeCloud<polyatomic> >
     {
+    public:
+
+        enum trackPart
+        {
+            tpFirstVelocityHalfStep,
+            tpLinearTrack,
+            tpRotationalTrack,
+            tpSecondVelocityHalfStep,
+            tpAccess
+        };
+
+
+    private:
+
         // label specifying which part of the integration algorithm is taking
         label part_;
 
@@ -158,7 +169,7 @@ public:
 
         // Constructors
 
-            trackingData(MoleculeCloud<polyatomic>& cloud, label part)
+            trackingData(MoleculeCloud<polyatomic>& cloud, trackPart part)
             :
                 particle::TrackingData<MoleculeCloud<polyatomic> >(cloud),
                 part_(part)
@@ -242,6 +253,23 @@ public:
 
     // Constructors
 
+        //- Construct with macroscopic description
+        inline polyatomic
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label cellI,
+            const label tetFaceI,
+            const label tetPtI,
+            const scalar temperature,
+            const vector& bulkVelocity,
+            const vector& specialPosition,
+            const constantProperties& constProps,
+            polyatomic::trackingData& td,
+            const label special,
+            const label id
+        );
+
         //- Construct from all components
         inline polyatomic
         (
@@ -431,6 +459,9 @@ public:
         //- Write
         static void writeFields(const Cloud<polyatomic>& mC);
 
+        //- Show info
+        static void info(trackingData& td);
+
 
     // IOstream Operators
 

src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomicI.H

@@ -305,6 +305,66 @@ inline Foam::polyatomic::constantProperties::constantProperties
 
 
 inline Foam::polyatomic::polyatomic
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label cellI,
+    const label tetFaceI,
+    const label tetPtI,
+    const scalar temperature,
+    const vector& bulkVelocity,
+    const vector& specialPosition,
+    const constantProperties& cP,
+    polyatomic::trackingData& td,
+    const label special,
+    const label id
+)
+:
+    particle(mesh, position, cellI, tetFaceI, tetPtI),
+    Q_(I),
+    v_(bulkVelocity),
+    a_(vector::zero),
+    pi_(vector::zero),
+    tau_(vector::zero),
+    specialPosition_(specialPosition),
+    potentialEnergy_(0.0),
+    rf_(tensor::zero),
+    special_(special),
+    id_(id),
+    siteForces_(cP.nSites(), vector::zero),
+    sitePositions_(cP.nSites())
+{
+    setSitePositions(cP);
+
+    v_ += td.cloud().equipartitionLinearVelocity(temperature, cP.mass());
+
+    if (!cP.pointMolecule())
+    {
+        pi_ = td.cloud().equipartitionAngularMomentum(temperature, cP);
+
+        scalar phi(td.cloud().rndGen().scalar01()*twoPi);
+
+        scalar theta(td.cloud().rndGen().scalar01()*twoPi);
+
+        scalar psi(td.cloud().rndGen().scalar01()*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)
+        );
+    }
+}
+
+
+Foam::polyatomic::polyatomic
 (
     const polyMesh& mesh,
     const vector& position,
@@ -317,10 +377,9 @@ inline Foam::polyatomic::polyatomic
     const vector& pi,
     const vector& tau,
     const vector& specialPosition,
-    const constantProperties& constProps,
+    const constantProperties& cP,
     const label special,
     const label id
-
 )
 :
     particle(mesh, position, cellI, tetFaceI, tetPtI),
@@ -334,10 +393,10 @@ inline Foam::polyatomic::polyatomic
     rf_(tensor::zero),
     special_(special),
     id_(id),
-    siteForces_(constProps.nSites(), vector::zero),
+    siteForces_(cP.nSites(), vector::zero),
-    sitePositions_(constProps.nSites())
+    sitePositions_(cP.nSites())
 {
-    setSitePositions(constProps);
+    setSitePositions(cP);
 }
 
 

src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomicIO.C

@@ -249,6 +249,115 @@ void Foam::polyatomic::writeFields(const Cloud<polyatomic>& mC)
 }
 
 
+void Foam::polyatomic::info(polyatomic::trackingData& td)
+{
+    vector totalLinearMomentum(vector::zero);
+    vector totalAngularMomentum(vector::zero);
+    scalar maxVelocityMag = 0.0;
+    scalar totalMass = 0.0;
+    scalar totalLinearKE = 0.0;
+    scalar totalAngularKE = 0.0;
+    scalar totalPE = 0.0;
+    scalar totalrDotf = 0.0;
+    //vector CentreOfMass(vector::zero);
+    label nMols = td.cloud().size();
+    label dofs = 0;
+
+    forAllConstIter(typename Cloud<polyatomic>, td.cloud(), mol)
+    {
+        const label molId = mol().id();
+        scalar molMass(td.cloud().constProps(molId).mass());
+        totalMass += molMass;
+        //CentreOfMass += mol().position()*molMass;
+    }
+
+    // if (nMols)
+    // {
+    //     CentreOfMass /= totalMass;
+    // }
+
+    forAllConstIter(typename Cloud<polyatomic>, td.cloud(), mol)
+    {
+        const label molId = mol().id();
+        const polyatomic::constantProperties cP
+        (
+            td.cloud().constProps(molId)
+        );
+        scalar molMass(cP.mass());
+        const diagTensor& molMoI(cP.momentOfInertia());
+        const vector& molV(mol().v());
+        const vector& molOmega(inv(molMoI) & mol().pi());
+        vector molPiGlobal = mol().Q() & mol().pi();
+        totalLinearMomentum += molV * molMass;
+        totalAngularMomentum += molPiGlobal;
+            //+((mol().position() - CentreOfMass) ^ (molV * molMass));
+        if (mag(molV) > maxVelocityMag)
+        {
+            maxVelocityMag = mag(molV);
+        }
+        totalLinearKE += 0.5*molMass*magSqr(molV);
+        totalAngularKE += 0.5*(molOmega & molMoI & molOmega);
+        totalPE += mol().potentialEnergy();
+        totalrDotf += tr(mol().rf());
+        dofs += cP.degreesOfFreedom();
+    }
+
+    scalar meshVolume = sum(td.cloud().mesh().cellVolumes());
+
+    if (Pstream::parRun())
+    {
+        reduce(totalLinearMomentum, sumOp<vector>());
+        reduce(totalAngularMomentum, sumOp<vector>());
+        reduce(maxVelocityMag, maxOp<scalar>());
+        reduce(totalMass, sumOp<scalar>());
+        reduce(totalLinearKE, sumOp<scalar>());
+        reduce(totalAngularKE, sumOp<scalar>());
+        reduce(totalPE, sumOp<scalar>());
+        reduce(totalrDotf, sumOp<scalar>());
+        reduce(nMols, sumOp<label>());
+        reduce(dofs, sumOp<label>());
+        reduce(meshVolume, sumOp<scalar>());
+    }
+
+    if (nMols)
+    {
+        Info<< nl << "Number of molecules in " << td.cloud().name() << " = "
+            << nMols << nl
+            << "    Overall number density = "
+            << nMols/meshVolume << nl
+            << "    Overall mass density = "
+            << totalMass/meshVolume << nl
+            << "    Average linear momentum per molecule = "
+            << totalLinearMomentum/nMols << ' '
+            << mag(totalLinearMomentum)/nMols << nl
+            << "    Average angular momentum per molecule = "
+            << totalAngularMomentum << ' '
+            << mag(totalAngularMomentum)/nMols << nl
+            << "    maximum |velocity| = "
+            << maxVelocityMag << nl
+            << "    Average linear KE per molecule = "
+            << totalLinearKE/nMols << nl
+            << "    Average angular KE per molecule = "
+            << totalAngularKE/nMols << nl
+            << "    Average PE per molecule = "
+            << totalPE/nMols << nl
+            << "    Average TE per molecule = "
+            <<
+            (
+                  totalLinearKE
+                + totalAngularKE
+                + totalPE
+            )
+            /nMols
+            << nl << endl;
+    }
+    else
+    {
+        Info<< nl << "No molecules in " << td.cloud().name() << endl;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 Foam::Ostream& Foam::operator<<(Ostream& os, const polyatomic& mol)

src/lagrangian/molecularDynamics/potential/potential/potential.C

@@ -243,9 +243,9 @@ void Foam::potential::potential::readPotentialDict
 
 void Foam::potential::potential::readMdInitialiseDict
 (
-    const IOdictionary& mdInitialiseDict,
+    const dictionary& mdInitialiseDict,
     const dictionary& moleculePropertiesDict,
-    IOdictionary& idListDict
+    dictionary& idListDict
 )
 {
     DynamicList<word> idList;
@@ -274,11 +274,12 @@ void Foam::potential::potential::readMdInitialiseDict
                 (
                     "potential::readMdInitialiseDict"
                     "("
-                        "const IOdictionary&, "
+                        "const dictionary&, "
-                        "IOdictionary&"
+                        "dictionary&"
                     ")"
                 )   << "Molecule type " << id
-                    << " not found in mdProperties dictionary." << nl
+                    << " not found in " << moleculePropertiesDict.name()
+                    << " dictionary." << nl
                     << abort(FatalError);
             }
 
@@ -329,8 +330,8 @@ void Foam::potential::potential::readMdInitialiseDict
                 (
                     "potential::readMdInitialiseDict"
                     "("
-                        "const IOdictionary&, "
+                        "const dictionary&, "
-                        "IOdictionary&"
+                        "dictionary&"
                     ")"
                 )   << "Tether id  " << tetherSiteId
                     << " not found as a site of any molecule in zone." << nl
@@ -367,9 +368,9 @@ Foam::potential::potential
 Foam::potential::potential
 (
     const polyMesh& mesh,
-    const IOdictionary& mdInitialiseDict,
+    const dictionary& mdInitialiseDict,
     const dictionary& moleculePropertiesDict,
-    IOdictionary& idListDict
+    dictionary& idListDict
 )
 :
     mesh_(mesh)

src/lagrangian/molecularDynamics/potential/potential/potential.H

@@ -95,9 +95,9 @@ class potential
         //-
         void readMdInitialiseDict
         (
-            const IOdictionary& mdInitialiseDict,
+            const dictionary& mdInitialiseDict,
             const dictionary& moleculePropertiesDict,
-            IOdictionary& idListDict
+            dictionary& idListDict
         );
 
         //- Disallow default bitwise copy construct
@@ -122,9 +122,9 @@ public:
         potential
         (
             const polyMesh& mesh,
-            const IOdictionary& mdInitialiseDict,
+            const dictionary& mdInitialiseDict,
             const dictionary& moleculePropertiesDict,
-            IOdictionary& idListDict
+            dictionary& idListDict
         );