diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdFoam/Make/options b/applications/solvers/discreteMethods/molecularDynamics/mdFoam/Make/options
index 89431d0625..98f09f1b4b 100644
--- a/applications/solvers/discreteMethods/molecularDynamics/mdFoam/Make/options
+++ b/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 \
diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdFoam/createFields.H b/applications/solvers/discreteMethods/molecularDynamics/mdFoam/createFields.H
index dc88bf6e88..d5ea5ef8b6 100644
--- a/applications/solvers/discreteMethods/molecularDynamics/mdFoam/createFields.H
+++ b/applications/solvers/discreteMethods/molecularDynamics/mdFoam/createFields.H
@@ -12,6 +12,44 @@
         mesh
     );
 
-    potential pot(mesh);
+    word polyatomicCloudName("polyatomicCloud");
 
-    polyatomicCloud molecules(mesh, pot);
+    potential polyPot
+    (
+        mesh,
+        IOdictionary
+        (
+            IOobject
+            (
+                polyatomicCloudName +  "Properties",
+                mesh.time().constant(),
+                mesh,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE,
+                false
+            )
+        )
+    );
+
+    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);
diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C b/applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C
index 660c12a5b8..38510dc453 100644
--- a/applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C
+++ b/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,9 +48,13 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << endl;
 
-        molecules.evolve();
+        monoatomics.evolve();
 
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+        polyatomics.evolve();
+
+        runTime.write();
+
+        Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
     }
diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/files b/applications/solvers/discreteMethods/molecularDynamics/old/mdEquilibrationFoam/Make/files
similarity index 100%
rename from applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/files
rename to applications/solvers/discreteMethods/molecularDynamics/old/mdEquilibrationFoam/Make/files
diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/options b/applications/solvers/discreteMethods/molecularDynamics/old/mdEquilibrationFoam/Make/options
similarity index 100%
rename from applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/options
rename to applications/solvers/discreteMethods/molecularDynamics/old/mdEquilibrationFoam/Make/options
diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C b/applications/solvers/discreteMethods/molecularDynamics/old/mdEquilibrationFoam/mdEquilibrationFoam.C
similarity index 100%
rename from applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C
rename to applications/solvers/discreteMethods/molecularDynamics/old/mdEquilibrationFoam/mdEquilibrationFoam.C
diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H b/applications/solvers/discreteMethods/molecularDynamics/old/mdEquilibrationFoam/readmdEquilibrationDict.H
similarity index 100%
rename from applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H
rename to applications/solvers/discreteMethods/molecularDynamics/old/mdEquilibrationFoam/readmdEquilibrationDict.H
diff --git a/applications/utilities/mesh/manipulation/topoSet/topoSet.C b/applications/utilities/mesh/manipulation/topoSet/topoSet.C
index 2f15bdce1d..8e4e33e8af 100644
--- a/applications/utilities/mesh/manipulation/topoSet/topoSet.C
+++ b/applications/utilities/mesh/manipulation/topoSet/topoSet.C
@@ -208,7 +208,8 @@ int main(int argc, char *argv[])
                     IOobject::MUST_READ
                 );
                 Info<< "Read set " << setName << " with size "
-                    << currentSet().size() << endl;
+                    << returnReduce(currentSet().size(), sumOp<label>())
+                    << endl;
             }
 
 
@@ -292,7 +293,8 @@ int main(int argc, char *argv[])
             if (currentSet.valid())
             {
                 Info<< "    Set " << currentSet().name()
-                    << " now size " << currentSet().size()
+                    << " now size "
+                    << returnReduce(currentSet().size(), sumOp<label>())
                     << endl;
             }
         }
diff --git a/applications/utilities/preProcessing/mdInitialise/Make/options b/applications/utilities/preProcessing/mdInitialise/Make/options
index 289c0990a2..7633b07554 100644
--- a/applications/utilities/preProcessing/mdInitialise/Make/options
+++ b/applications/utilities/preProcessing/mdInitialise/Make/options
@@ -3,7 +3,6 @@ EXE_INC = \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \
     -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
-    -I$(LIB_SRC)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
@@ -13,6 +12,5 @@ EXE_LIBS = \
     -lfiniteVolume \
     -llagrangian \
     -lmolecule \
-    -lpotential \
-    -lmolecularMeasurements
+    -lpotential
 
diff --git a/applications/utilities/preProcessing/mdInitialise/mdInitialise.C b/applications/utilities/preProcessing/mdInitialise/mdInitialise.C
index 3551e46133..d29ef0a72d 100644
--- a/applications/utilities/preProcessing/mdInitialise/mdInitialise.C
+++ b/applications/utilities/preProcessing/mdInitialise/mdInitialise.C
@@ -2,7 +2,7 @@
  =========                   |
  \\      /   F ield          | OpenFOAM: The Open Source CFD Toolbox
   \\    /    O peration      |
-   \\  /     A nd            | Copyright (C) 2008-2010 OpenCFD Ltd.
+   \\  /     A nd            | Copyright (C) 2008-2011 OpenCFD Ltd.
     \\/      M anipulation   |
 -------------------------------------------------------------------------------
 License
@@ -26,8 +26,9 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#include "md.H"
 #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,26 +67,88 @@ int main(int argc, char *argv[])
         )
     );
 
-    potential pot(mesh, mdInitialiseDict, idListDict);
+    potential polyPot
+    (
+        mesh,
+        polyDict,
+        IOdictionary
+        (
+            IOobject
+            (
+                polyCloudName +  "Properties",
+                mesh.time().constant(),
+                mesh,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE,
+                false
+            )
+        ),
+        polyIdListDict
+    );
 
-    moleculeCloud molecules(mesh, pot, mdInitialiseDict);
+    polyatomicCloud poly
+    (
+        polyCloudName,
+        mesh,
+        polyPot,
+        polyDict
+    );
 
-    label totalMolecules = molecules.size();
-
-    if (Pstream::parRun())
-    {
-        reduce(totalMolecules, sumOp<label>());
-    }
-
-    Info<< nl << "Total number of molecules added: " << totalMolecules
+    Info<< nl << returnReduce(poly.size(), sumOp<label>()) << " added to "
+        << poly.name()
         << 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 moleculeCloud."
+            << "Failed writing."
             << nl << exit(FatalError);
     }
 
diff --git a/src/OSspecific/POSIX/POSIX.C b/src/OSspecific/POSIX/POSIX.C
index 48e4b3a4df..dba1e73f3a 100644
--- a/src/OSspecific/POSIX/POSIX.C
+++ b/src/OSspecific/POSIX/POSIX.C
@@ -57,6 +57,10 @@ Description
 
 #include <netinet/in.h>
 
+#ifdef USE_RANDOM
+#   include <climits>
+#endif
+
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 defineTypeNameAndDebug(Foam::POSIX, 0);
@@ -68,16 +72,19 @@ pid_t Foam::pid()
     return ::getpid();
 }
 
+
 pid_t Foam::ppid()
 {
     return ::getppid();
 }
 
+
 pid_t Foam::pgid()
 {
     return ::getpgrp();
 }
 
+
 bool Foam::env(const word& envName)
 {
     return ::getenv(envName.c_str()) != NULL;
@@ -890,7 +897,6 @@ bool Foam::mvBak(const fileName& src, const std::string& ext)
 }
 
 
-
 // Remove a file, returning true if successful otherwise false
 bool Foam::rm(const fileName& file)
 {
@@ -1221,4 +1227,34 @@ Foam::fileNameList Foam::dlLoaded()
 }
 
 
+void Foam::osRandomSeed(const label seed)
+{
+#ifdef USE_RANDOM
+    srandom((unsigned int)seed);
+#else
+    srand48(seed);
+#endif
+}
+
+
+Foam::label Foam::osRandomInteger()
+{
+#ifdef USE_RANDOM
+    return random();
+#else
+    return lrand48();
+#endif
+}
+
+
+Foam::scalar Foam::osRandomDouble()
+{
+#ifdef USE_RANDOM
+    return (scalar)random();
+#else
+    return drand48();
+#endif
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/include/OSspecific.H b/src/OpenFOAM/include/OSspecific.H
index f01ad29647..b774997339 100644
--- a/src/OpenFOAM/include/OSspecific.H
+++ b/src/OpenFOAM/include/OSspecific.H
@@ -200,6 +200,18 @@ bool dlSymFound(void* handle, const std::string& symbol);
 fileNameList dlLoaded();
 
 
+// Low level random numbers. Use Random class instead.
+
+//- Seed random number generator.
+void osRandomSeed(const label seed);
+
+//- Return random integer (uniform distribution between 0 and 2^31)
+label osRandomInteger();
+
+//- Return random double precision (uniform distribution between 0 and 1)
+scalar osRandomDouble();
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/OpenFOAM/primitives/random/Random/Random.C b/src/OpenFOAM/primitives/random/Random/Random.C
index 704e5ef118..e2496ee648 100644
--- a/src/OpenFOAM/primitives/random/Random/Random.C
+++ b/src/OpenFOAM/primitives/random/Random/Random.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "Random.H"
+#include "OSspecific.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -37,11 +38,6 @@ namespace Foam
 #    error "The random number generator may not work!"
 #endif
 
-#ifdef USE_RANDOM
-#   include <climits>
-#else
-#   include <cstdlib>
-#endif
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -57,22 +53,13 @@ Random::Random(const label seed)
         Seed = 1;
     }
 
-#   ifdef USE_RANDOM
-        srandom((unsigned int)Seed);
-#   else
-        srand48(Seed);
-#   endif
-
+    osRandomSeed(Seed);
 }
 
 
 int Random::bit()
 {
-#   ifdef USE_RANDOM
-    if (random() > INT_MAX/2)
-#   else
-    if (lrand48() > INT_MAX/2)
-#   endif
+    if (osRandomInteger() > INT_MAX/2)
     {
         return 1;
     }
@@ -85,11 +72,7 @@ int Random::bit()
 
 scalar Random::scalar01()
 {
-#   ifdef USE_RANDOM
-        return (scalar)random()/INT_MAX;
-#   else
-        return drand48();
-#   endif
+    return osRandomDouble();
 }
 
 
@@ -140,11 +123,7 @@ tensor Random::tensor01()
 
 label Random::integer(const label lower, const label upper)
 {
-#   ifdef USE_RANDOM
-        return lower + (random() % (upper+1-lower));
-#   else
-        return lower + (lrand48() % (upper+1-lower));
-#   endif
+    return lower + (osRandomInteger() % (upper+1-lower));
 }
 
 
diff --git a/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.C b/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.C
index b48c1f4e63..ac0c270503 100644
--- a/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.C
+++ b/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,7 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "cachedRandom.H"
-#include <cstdlib>
+#include "OSspecific.H"
 
 #if INT_MAX    != 2147483647
 #    error "INT_MAX    != 2147483647"
@@ -37,7 +37,7 @@ Foam::scalar Foam::cachedRandom::scalar01()
 {
     if (sampleI_ < 0)
     {
-        return drand48();
+        return osRandomDouble();
     }
 
     if (sampleI_ == samples_.size() - 1)
@@ -76,7 +76,7 @@ Foam::cachedRandom::cachedRandom(const label seed, const label count)
     }
 
     // Initialise samples
-    srand48(seed_);
+    osRandomSeed(seed_);
     forAll(samples_, i)
     {
         samples_[i] = drand48();
@@ -98,7 +98,7 @@ Foam::cachedRandom::cachedRandom(const cachedRandom& cr, const bool reset)
         )   << "Copy constructor called, but samples not being cached. "
             << "This may lead to non-repeatable behaviour" << endl;
 
-        srand48(seed_);
+        osRandomSeed(seed_);
     }
     else if (reset)
     {
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fanPressure/fanPressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/fanPressure/fanPressureFvPatchScalarField.C
index db585345f1..260fb328dc 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fanPressure/fanPressureFvPatchScalarField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fanPressure/fanPressureFvPatchScalarField.C
@@ -59,10 +59,7 @@ Foam::fanPressureFvPatchScalarField::fanPressureFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedValueFvPatchScalarField(p, iF),
-    phiName_("phi"),
-    rhoName_("rho"),
-    p0_(p.size(), 0.0),
+    totalPressureFvPatchScalarField(p, iF),
     fanCurve_(),
     direction_(ffdOut)
 {}
@@ -76,10 +73,7 @@ Foam::fanPressureFvPatchScalarField::fanPressureFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
-    phiName_(ptf.phiName_),
-    rhoName_(ptf.rhoName_),
-    p0_(ptf.p0_, mapper),
+    totalPressureFvPatchScalarField(ptf, p, iF, mapper),
     fanCurve_(ptf.fanCurve_),
     direction_(ptf.direction_)
 {}
@@ -92,10 +86,7 @@ Foam::fanPressureFvPatchScalarField::fanPressureFvPatchScalarField
     const dictionary& dict
 )
 :
-    fixedValueFvPatchScalarField(p, iF),
-    phiName_(dict.lookupOrDefault<word>("phi", "phi")),
-    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
-    p0_("p0", dict, p.size()),
+    totalPressureFvPatchScalarField(p, iF),
     fanCurve_(dict),
     direction_(fanFlowDirectionNames_.read(dict.lookup("direction")))
 {
@@ -109,10 +100,7 @@ Foam::fanPressureFvPatchScalarField::fanPressureFvPatchScalarField
     const fanPressureFvPatchScalarField& pfopsf
 )
 :
-    fixedValueFvPatchScalarField(pfopsf),
-    phiName_(pfopsf.phiName_),
-    rhoName_(pfopsf.rhoName_),
-    p0_(pfopsf.p0_),
+    totalPressureFvPatchScalarField(pfopsf),
     fanCurve_(pfopsf.fanCurve_),
     direction_(pfopsf.direction_)
 {}
@@ -124,10 +112,7 @@ Foam::fanPressureFvPatchScalarField::fanPressureFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedValueFvPatchScalarField(pfopsf, iF),
-    phiName_(pfopsf.phiName_),
-    rhoName_(pfopsf.rhoName_),
-    p0_(pfopsf.p0_),
+    totalPressureFvPatchScalarField(pfopsf, iF),
     fanCurve_(pfopsf.fanCurve_),
     direction_(pfopsf.direction_)
 {}
@@ -144,7 +129,7 @@ void Foam::fanPressureFvPatchScalarField::updateCoeffs()
 
     // Retrieve flux field
     const surfaceScalarField& phi =
-        db().lookupObject<surfaceScalarField>(phiName_);
+        db().lookupObject<surfaceScalarField>(phiName());
 
     const fvsPatchField<scalar>& phip =
         patch().patchField<surfaceScalarField, scalar>(phi);
@@ -161,7 +146,7 @@ void Foam::fanPressureFvPatchScalarField::updateCoeffs()
     else if (phi.dimensions() == dimVelocity*dimArea*dimDensity)
     {
         const scalarField& rhop =
-            patch().lookupPatchField<volScalarField, scalar>(rhoName_);
+            patch().lookupPatchField<volScalarField, scalar>(rhoName());
         aveFlowRate = dir*gSum(phip/rhop)/gSum(patch().magSf());
     }
     else
@@ -174,51 +159,23 @@ void Foam::fanPressureFvPatchScalarField::updateCoeffs()
                 << exit(FatalError);
     }
 
-    // Normal flow through fan
-    if (aveFlowRate >= 0.0)
-    {
-        // Pressure drop for this flow rate
-        const scalar pdFan = fanCurve_(aveFlowRate);
+    // Pressure drop for this flow rate
+    const scalar pdFan = fanCurve_(max(aveFlowRate, 0.0));
 
-        operator==(p0_ - dir*pdFan);
-    }
-    // Reverse flow
-    else
-    {
-        // Assume that fan has stalled if flow reversed
-        // i.e. apply dp for zero flow rate
-        const scalar pdFan = fanCurve_(0);
-
-        // Flow speed across patch
-        scalarField Up = phip/(patch().magSf());
-
-        // Pressure drop associated withback flow = dynamic pressure
-        scalarField pdBackFlow = 0.5*magSqr(Up);
-
-        if (phi.dimensions() == dimVelocity*dimArea*dimDensity)
-        {
-            const scalarField& rhop =
-                patch().lookupPatchField<volScalarField, scalar>(rhoName_);
-            pdBackFlow /= rhop;
-        }
-
-        operator==(p0_ - dir*(pdBackFlow + pdFan));
-    }
-
-    fixedValueFvPatchScalarField::updateCoeffs();
+    totalPressureFvPatchScalarField::updateCoeffs
+    (
+        p0() - dir*pdFan,
+        patch().lookupPatchField<volVectorField, vector>(UName())
+    );
 }
 
 
 void Foam::fanPressureFvPatchScalarField::write(Ostream& os) const
 {
-    fvPatchScalarField::write(os);
-    os.writeKeyword("phi") << phiName_ << token::END_STATEMENT << nl;
-    os.writeKeyword("rho") << rhoName_ << token::END_STATEMENT << nl;
+    totalPressureFvPatchScalarField::write(os);
     fanCurve_.write(os);
     os.writeKeyword("direction")
         << fanFlowDirectionNames_[direction_] << token::END_STATEMENT << nl;
-    p0_.writeEntry("p0", os);
-    writeEntry("value", os);
 }
 
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fanPressure/fanPressureFvPatchScalarField.H b/src/finiteVolume/fields/fvPatchFields/derived/fanPressure/fanPressureFvPatchScalarField.H
index 7c3d9ff411..47f5cdd49c 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fanPressure/fanPressureFvPatchScalarField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fanPressure/fanPressureFvPatchScalarField.H
@@ -25,7 +25,7 @@ Class
     Foam::fanPressureFvPatchScalarField
 
 Description
-    Assigns pressure inlet or outlet condition for a fan.
+    Assigns pressure inlet or outlet total pressure condition for a fan.
 
     User specifies:
     - pressure drop vs volumetric flow rate table (fan curve) file name;
@@ -56,8 +56,8 @@ Description
     \endverbatim
 
 See Also
-    Foam::interpolationTable and
-    Foam::timeVaryingFlowRateInletVelocityFvPatchVectorField
+    Foam::totalPressureFvPatchScalarField and
+    Foam::interpolationTable
 
 SourceFiles
     fanPressureFvPatchScalarField.C
@@ -67,8 +67,7 @@ SourceFiles
 #ifndef fanPressureFvPatchScalarField_H
 #define fanPressureFvPatchScalarField_H
 
-#include "fvPatchFields.H"
-#include "fixedValueFvPatchFields.H"
+#include "totalPressureFvPatchScalarField.H"
 #include "interpolationTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -82,19 +81,10 @@ namespace Foam
 
 class fanPressureFvPatchScalarField
 :
-    public fixedValueFvPatchScalarField
+    public totalPressureFvPatchScalarField
 {
     // Private data
 
-        //- Name of the flux transporting the field
-        word phiName_;
-
-        //- Name of the density field
-        word rhoName_;
-
-        //- Total pressure
-        scalarField p0_;
-
         //- Tabulated fan curve
         interpolationTable<scalar> fanCurve_;
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/rotatingTotalPressure/rotatingTotalPressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/rotatingTotalPressure/rotatingTotalPressureFvPatchScalarField.C
index 9f3e49576b..f4d557d694 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/rotatingTotalPressure/rotatingTotalPressureFvPatchScalarField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/rotatingTotalPressure/rotatingTotalPressureFvPatchScalarField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -112,7 +112,7 @@ void Foam::rotatingTotalPressureFvPatchScalarField::updateCoeffs()
       + rotationVelocity
     );
 
-    totalPressureFvPatchScalarField::updateCoeffs(Up);
+    totalPressureFvPatchScalarField::updateCoeffs(p0(), Up);
 }
 
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/rotatingTotalPressure/rotatingTotalPressureFvPatchScalarField.H b/src/finiteVolume/fields/fvPatchFields/derived/rotatingTotalPressure/rotatingTotalPressureFvPatchScalarField.H
index 087bf70b4a..01d9a73ad4 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/rotatingTotalPressure/rotatingTotalPressureFvPatchScalarField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/rotatingTotalPressure/rotatingTotalPressureFvPatchScalarField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.C
index 7defeec931..e60f802e5a 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -153,7 +153,11 @@ void Foam::totalPressureFvPatchScalarField::rmap
 }
 
 
-void Foam::totalPressureFvPatchScalarField::updateCoeffs(const vectorField& Up)
+void Foam::totalPressureFvPatchScalarField::updateCoeffs
+(
+    const scalarField& p0p,
+    const vectorField& Up
+)
 {
     if (updated())
     {
@@ -165,7 +169,7 @@ void Foam::totalPressureFvPatchScalarField::updateCoeffs(const vectorField& Up)
 
     if (psiName_ == "none" && rhoName_ == "none")
     {
-        operator==(p0_ - 0.5*(1.0 - pos(phip))*magSqr(Up));
+        operator==(p0p - 0.5*(1.0 - pos(phip))*magSqr(Up));
     }
     else if (rhoName_ == "none")
     {
@@ -178,7 +182,7 @@ void Foam::totalPressureFvPatchScalarField::updateCoeffs(const vectorField& Up)
 
             operator==
             (
-                p0_
+                p0p
                /pow
                 (
                     (1.0 + 0.5*psip*gM1ByG*(1.0 - pos(phip))*magSqr(Up)),
@@ -188,7 +192,7 @@ void Foam::totalPressureFvPatchScalarField::updateCoeffs(const vectorField& Up)
         }
         else
         {
-            operator==(p0_/(1.0 + 0.5*psip*(1.0 - pos(phip))*magSqr(Up)));
+            operator==(p0p/(1.0 + 0.5*psip*(1.0 - pos(phip))*magSqr(Up)));
         }
     }
     else if (psiName_ == "none")
@@ -196,7 +200,7 @@ void Foam::totalPressureFvPatchScalarField::updateCoeffs(const vectorField& Up)
         const fvPatchField<scalar>& rho =
             patch().lookupPatchField<volScalarField, scalar>(rhoName_);
 
-        operator==(p0_ - 0.5*rho*(1.0 - pos(phip))*magSqr(Up));
+        operator==(p0p - 0.5*rho*(1.0 - pos(phip))*magSqr(Up));
     }
     else
     {
@@ -220,7 +224,11 @@ void Foam::totalPressureFvPatchScalarField::updateCoeffs(const vectorField& Up)
 
 void Foam::totalPressureFvPatchScalarField::updateCoeffs()
 {
-    updateCoeffs(patch().lookupPatchField<volVectorField, vector>(UName_));
+    updateCoeffs
+    (
+        p0(),
+        patch().lookupPatchField<volVectorField, vector>(UName())
+    );
 }
 
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.H b/src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.H
index 9aa386fe83..1eceb87faa 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -157,6 +157,45 @@ public:
                 return UName_;
             }
 
+            //- Return the name of the flux field
+            const word& phiName() const
+            {
+                return phiName_;
+            }
+
+            //- Return reference to the name of the flux field
+            //  to allow adjustment
+            word& phiName()
+            {
+                return phiName_;
+            }
+
+            //- Return the name of the density field
+            const word& rhoName() const
+            {
+                return rhoName_;
+            }
+
+            //- Return reference to the name of the density field
+            //  to allow adjustment
+            word& rhoName()
+            {
+                return rhoName_;
+            }
+
+            //- Return the name of the compressibility field
+            const word& psiName() const
+            {
+                return psiName_;
+            }
+
+            //- Return reference to the name of the compressibility field
+            //  to allow adjustment
+            word& psiName()
+            {
+                return psiName_;
+            }
+
             //- Return the heat capacity ratio
             scalar gamma() const
             {
@@ -201,8 +240,12 @@ public:
         // Evaluation functions
 
             //- Update the coefficients associated with the patch field
-            //  using the given patch velocity field
-            virtual void updateCoeffs(const vectorField& Up);
+            //  using the given patch total pressure and velocity fields
+            virtual void updateCoeffs
+            (
+                const scalarField& p0p,
+                const vectorField& Up
+            );
 
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.C b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.C
index 812296e6f5..5580d3cb05 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.C
@@ -214,7 +214,8 @@ void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
             rHat_AB
            *(kN*pow(normalOverlapMag, b_) - etaN*(U_AB & rHat_AB));
 
-        // Cohesion force
+        // Cohesion force, energy density multiplied by the area of
+        // particle-particle overlap
         if (cohesion_)
         {
             fN_AB +=
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.C b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.C
index 9ae09dd98f..946ec3adb8 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -76,7 +76,8 @@ void Foam::WallLocalSpringSliderDashpot<CloudType>::evaluateWall
     typename CloudType::parcelType& p,
     const point& site,
     const WallSiteData<vector>& data,
-    scalar pREff
+    scalar pREff,
+    bool cohesion
 ) const
 {
     // wall patch index
@@ -88,14 +89,18 @@ void Foam::WallLocalSpringSliderDashpot<CloudType>::evaluateWall
     scalar alpha = alpha_[wPI];
     scalar b = b_[wPI];
     scalar mu = mu_[wPI];
+    scalar cohesionEnergyDensity = cohesionEnergyDensity_[wPI];
+    cohesion = cohesion && cohesion_[wPI];
 
     vector r_PW = p.position() - site;
 
     vector U_PW = p.U() - data.wallData();
 
-    scalar normalOverlapMag = max(pREff - mag(r_PW), 0.0);
+    scalar r_PW_mag = mag(r_PW);
 
-    vector rHat_PW = r_PW/(mag(r_PW) + VSMALL);
+    scalar normalOverlapMag = max(pREff - r_PW_mag, 0.0);
+
+    vector rHat_PW = r_PW/(r_PW_mag + VSMALL);
 
     scalar kN = (4.0/3.0)*sqrt(pREff)*Estar;
 
@@ -105,6 +110,16 @@ void Foam::WallLocalSpringSliderDashpot<CloudType>::evaluateWall
         rHat_PW
        *(kN*pow(normalOverlapMag, b) - etaN*(U_PW & rHat_PW));
 
+    // Cohesion force, energy density multiplied by the area of wall/particle
+    // overlap
+    if (cohesion)
+    {
+        fN_PW +=
+           -cohesionEnergyDensity
+           *mathematical::pi*(sqr(pREff) - sqr(r_PW_mag))
+           *rHat_PW;
+    }
+
     p.f() += fN_PW;
 
     vector USlip_PW =
@@ -168,6 +183,8 @@ Foam::WallLocalSpringSliderDashpot<CloudType>::WallLocalSpringSliderDashpot
     alpha_(),
     b_(),
     mu_(),
+    cohesionEnergyDensity_(),
+    cohesion_(),
     patchMap_(),
     maxEstarIndex_(-1),
     collisionResolutionSteps_
@@ -212,6 +229,8 @@ Foam::WallLocalSpringSliderDashpot<CloudType>::WallLocalSpringSliderDashpot
     alpha_.setSize(nWallPatches);
     b_.setSize(nWallPatches);
     mu_.setSize(nWallPatches);
+    cohesionEnergyDensity_.setSize(nWallPatches);
+    cohesion_.setSize(nWallPatches);
 
     scalar maxEstar = -GREAT;
 
@@ -238,6 +257,13 @@ Foam::WallLocalSpringSliderDashpot<CloudType>::WallLocalSpringSliderDashpot
 
         mu_[wPI] = readScalar(patchCoeffDict.lookup("mu"));
 
+        cohesionEnergyDensity_[wPI] = readScalar
+        (
+            patchCoeffDict.lookup("cohesionEnergyDensity")
+        );
+
+        cohesion_[wPI] = (mag(cohesionEnergyDensity_[wPI]) > VSMALL);
+
         if (Estar_[wPI] > maxEstar)
         {
             maxEstarIndex_ = wPI;
@@ -325,20 +351,22 @@ void Foam::WallLocalSpringSliderDashpot<CloudType>::evaluateWall
             p,
             flatSitePoints[siteI],
             flatSiteData[siteI],
-            pREff
+            pREff,
+            true
         );
     }
 
     forAll(sharpSitePoints, siteI)
     {
-        // Treating sharp sites like flat sites
+        // Treating sharp sites like flat sites, except suppress cohesion
 
         evaluateWall
         (
             p,
             sharpSitePoints[siteI],
             sharpSiteData[siteI],
-            pREff
+            pREff,
+            false
         );
     }
 }
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.H b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.H
index 700ce321ab..c4ad801b13 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.H
@@ -65,6 +65,12 @@ class WallLocalSpringSliderDashpot
         //- Coefficient of friction in for tangential sliding
         scalarList mu_;
 
+        //- Cohesion energy density [J/m^3]
+        scalarList cohesionEnergyDensity_;
+
+        // Switch cohesion on and off
+        boolList cohesion_;
+
         //- Mapping the patch index to the model data
         labelList patchMap_;
 
@@ -115,7 +121,8 @@ class WallLocalSpringSliderDashpot
             typename CloudType::parcelType& p,
             const point& site,
             const WallSiteData<vector>& data,
-            scalar pREff
+            scalar pREff,
+            bool cohesion
         ) const;
 
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.C b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.C
index 5f14b8e274..1f547f0d5d 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -77,16 +77,19 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
     const point& site,
     const WallSiteData<vector>& data,
     scalar pREff,
-    scalar kN
+    scalar kN,
+    bool cohesion
 ) const
 {
     vector r_PW = p.position() - site;
 
     vector U_PW = p.U() - data.wallData();
 
-    scalar normalOverlapMag = max(pREff - mag(r_PW), 0.0);
+    scalar r_PW_mag = mag(r_PW);
 
-    vector rHat_PW = r_PW/(mag(r_PW) + VSMALL);
+    scalar normalOverlapMag = max(pREff - r_PW_mag, 0.0);
+
+    vector rHat_PW = r_PW/(r_PW_mag + VSMALL);
 
     scalar etaN = alpha_*sqrt(p.mass()*kN)*pow025(normalOverlapMag);
 
@@ -94,6 +97,16 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
         rHat_PW
        *(kN*pow(normalOverlapMag, b_) - etaN*(U_PW & rHat_PW));
 
+    // Cohesion force, energy density multiplied by the area of wall/particle
+    // overlap
+    if (cohesion)
+    {
+        fN_PW +=
+           -cohesionEnergyDensity_
+           *mathematical::pi*(sqr(pREff) - sqr(r_PW_mag))
+           *rHat_PW;
+    }
+
     p.f() += fN_PW;
 
     vector USlip_PW =
@@ -157,6 +170,11 @@ Foam::WallSpringSliderDashpot<CloudType>::WallSpringSliderDashpot
     alpha_(readScalar(this->coeffDict().lookup("alpha"))),
     b_(readScalar(this->coeffDict().lookup("b"))),
     mu_(readScalar(this->coeffDict().lookup("mu"))),
+    cohesionEnergyDensity_
+    (
+        readScalar(this->coeffDict().lookup("cohesionEnergyDensity"))
+    ),
+    cohesion_(false),
     collisionResolutionSteps_
     (
         readScalar
@@ -183,6 +201,8 @@ Foam::WallSpringSliderDashpot<CloudType>::WallSpringSliderDashpot
     Estar_ = 1/((1 - sqr(pNu))/pE + (1 - sqr(nu))/E);
 
     Gstar_ = 1/(2*((2 + pNu - sqr(pNu))/pE + (2 + nu - sqr(nu))/E));
+
+    cohesion_ = (mag(cohesionEnergyDensity_) > VSMALL);
 }
 
 
@@ -266,13 +286,14 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
             flatSitePoints[siteI],
             flatSiteData[siteI],
             pREff,
-            kN
+            kN,
+            cohesion_
         );
     }
 
     forAll(sharpSitePoints, siteI)
     {
-        // Treating sharp sites like flat sites
+        // Treating sharp sites like flat sites, except suppress cohesion
 
         evaluateWall
         (
@@ -280,7 +301,8 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
             sharpSitePoints[siteI],
             sharpSiteData[siteI],
             pREff,
-            kN
+            kN,
+            false
         );
     }
 }
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.H b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.H
index 10010c60c3..e6e93ddb24 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.H
@@ -65,6 +65,12 @@ class WallSpringSliderDashpot
         //- Coefficient of friction in for tangential sliding
         scalar mu_;
 
+        //- Cohesion energy density [J/m^3]
+        scalar cohesionEnergyDensity_;
+
+        // Switch cohesion on and off
+        bool cohesion_;
+
         //- The number of steps over which to resolve the minimum
         //  harmonic approximation of the collision period
         scalar collisionResolutionSteps_;
@@ -110,7 +116,8 @@ class WallSpringSliderDashpot
             const point& site,
             const WallSiteData<vector>& data,
             scalar pREff,
-            scalar kN
+            scalar kN,
+            bool cohesion
         ) const;
 
 
diff --git a/src/lagrangian/molecularDynamics/molecule/Make/files b/src/lagrangian/molecularDynamics/molecule/Make/files
index 7f8abcb0f3..e71272b545 100644
--- a/src/lagrangian/molecularDynamics/molecule/Make/files
+++ b/src/lagrangian/molecularDynamics/molecule/Make/files
@@ -1,8 +1,14 @@
-constPropSite/constPropSite.C
+clouds/baseClasses/moleculeCloud/moleculeCloud.C
 
-polyatomic/polyatomic/polyatomic.C
-polyatomic/polyatomic/polyatomicIO.C
+molecules/constPropSite/constPropSite.C
 
-polyatomic/polyatomicCloud/polyatomicCloud.C
+molecules/monoatomic/monoatomic.C
+molecules/monoatomic/monoatomicIO.C
+molecules/polyatomic/polyatomic.C
+molecules/polyatomic/polyatomicIO.C
+
+/* controllers/basic/controllers/controllers.C
+controllers/basic/stateController/stateController.C
+controllers/basic/fluxController/fluxController.C */
 
 LIB = $(FOAM_LIBBIN)/libmolecule
diff --git a/src/lagrangian/molecularDynamics/molecule/Make/options b/src/lagrangian/molecularDynamics/molecule/Make/options
index 5698f81a65..78ecf2089f 100644
--- a/src/lagrangian/molecularDynamics/molecule/Make/options
+++ b/src/lagrangian/molecularDynamics/molecule/Make/options
@@ -2,13 +2,11 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
-    -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
-    -I$(LIB_SRC)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
     -llagrangian \
-    -lpotential \
-    -lmolecularMeasurements
+    -lpotential
 
diff --git a/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloud.C b/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.C
similarity index 80%
rename from src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloud.C
rename to src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.C
index a49d9d2319..c6d54036f7 100644
--- a/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloud.C
+++ b/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.C
@@ -23,36 +23,29 @@ License
 
 \*----------------------------------------------------------------------------*/
 
-#include "polyatomicCloud.H"
+#include "MoleculeCloud.H"
 #include "fvMesh.H"
 #include "mathematicalConstants.H"
 
 using namespace Foam::constant::mathematical;
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTemplateTypeNameAndDebug(Cloud<polyatomic>, 0);
-}
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::polyatomicCloud::buildConstProps()
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::buildConstProps()
 {
-    Info<< nl << "Reading polyatomicProperties dictionary." << endl;
-
     const List<word>& idList(pot_.idList());
 
     constPropList_.setSize(idList.size());
 
     const List<word>& siteIdList(pot_.siteIdList());
 
-    IOdictionary polyatomicPropertiesDict
+    IOdictionary molPropertiesDict
     (
         IOobject
         (
-            "polyatomicProperties",
+            this->name() + "Properties",
             mesh_.time().constant(),
             mesh_,
             IOobject::MUST_READ_IF_MODIFIED,
@@ -61,11 +54,13 @@ void Foam::polyatomicCloud::buildConstProps()
         )
     );
 
+    Info<< nl << "Reading " << molPropertiesDict.name() << endl;
+
     forAll(idList, i)
     {
         const word& id = idList[i];
 
-        const dictionary& molDict(polyatomicPropertiesDict.subDict(id));
+        const dictionary& molDict(molPropertiesDict.subDict(id));
 
         List<word> siteIdNames = molDict.lookup("siteIds");
 
@@ -81,17 +76,14 @@ void Foam::polyatomicCloud::buildConstProps()
             {
                 FatalErrorIn
                 (
-                    "Foam::polyatomic::constantProperties::constantProperties"
-                    "("
-                    "const dictionary& dict"
-                    ")"
+                    "void Foam::MoleculeCloud<MoleculeType>::buildConstProps()"
                 )
                     << siteId << " site not found."
                     << nl << abort(FatalError);
             }
         }
 
-        constPropList_[i] = polyatomic::constantProperties
+        constPropList_[i] = typename MoleculeType::constantProperties
         (
             molDict,
             siteIds
@@ -100,11 +92,15 @@ void Foam::polyatomicCloud::buildConstProps()
 }
 
 
-void Foam::polyatomicCloud::setSiteSizesAndPositions()
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::setSiteSizesAndPositions()
 {
-    forAllIter(polyatomicCloud, *this, mol)
+    forAllIter(typename MoleculeCloud<MoleculeType>, *this, mol)
     {
-        const polyatomic::constantProperties& cP = constProps(mol().id());
+        const typename MoleculeType::constantProperties& cP
+        (
+            constProps(mol().id())
+        );
 
         mol().setSiteSizes(cP.nSites());
 
@@ -113,14 +109,15 @@ void Foam::polyatomicCloud::setSiteSizesAndPositions()
 }
 
 
-void Foam::polyatomicCloud::buildCellOccupancy()
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::buildCellOccupancy()
 {
     forAll(cellOccupancy_, cO)
     {
         cellOccupancy_[cO].clear();
     }
 
-    forAllIter(polyatomicCloud, *this, mol)
+    forAllIter(typename MoleculeCloud<MoleculeType>, *this, mol)
     {
         cellOccupancy_[mol().cell()].append(&mol());
     }
@@ -132,15 +129,16 @@ void Foam::polyatomicCloud::buildCellOccupancy()
 }
 
 
-void Foam::polyatomicCloud::calculatePairForce()
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::calculatePairForce()
 {
     PstreamBuffers pBufs(Pstream::nonBlocking);
 
     // Start sending referred data
     il_.sendReferredData(cellOccupancy(), pBufs);
 
-    polyatomic* molI = NULL;
-    polyatomic* molJ = NULL;
+    MoleculeType* molI = NULL;
+    MoleculeType* molJ = NULL;
 
     {
         // Real-Real interactions
@@ -155,7 +153,7 @@ void Foam::polyatomicCloud::calculatePairForce()
 
                 forAll(dil[d], interactingCells)
                 {
-                    List<polyatomic*> cellJ =
+                    List<MoleculeType*> cellJ =
                         cellOccupancy_[dil[d][interactingCells]];
 
                     forAll(cellJ, cellJMols)
@@ -187,24 +185,24 @@ void Foam::polyatomicCloud::calculatePairForce()
 
         const labelListList& ril = il_.ril();
 
-        List<IDLList<polyatomic> >& referredMols = il_.referredParticles();
+        List<IDLList<MoleculeType> >& referredMols = il_.referredParticles();
 
         forAll(ril, r)
         {
             const List<label>& realCells = ril[r];
 
-            IDLList<polyatomic>& refMols = referredMols[r];
+            IDLList<MoleculeType>& refMols = referredMols[r];
 
             forAllIter
             (
-                IDLList<polyatomic>,
+                typename IDLList<MoleculeType>,
                 refMols,
                 refMol
             )
             {
                 forAll(realCells, rC)
                 {
-                    List<polyatomic*> cellI = cellOccupancy_[realCells[rC]];
+                    List<MoleculeType*> cellI = cellOccupancy_[realCells[rC]];
 
                     forAll(cellI, cellIMols)
                     {
@@ -219,11 +217,12 @@ void Foam::polyatomicCloud::calculatePairForce()
 }
 
 
-void Foam::polyatomicCloud::calculateTetherForce()
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::calculateTetherForce()
 {
     const tetherPotentialList& tetherPot(pot_.tetherPotentials());
 
-    forAllIter(polyatomicCloud, *this, mol)
+    forAllIter(typename MoleculeCloud<MoleculeType>, *this, mol)
     {
         if (mol().tethered())
         {
@@ -246,16 +245,18 @@ void Foam::polyatomicCloud::calculateTetherForce()
 }
 
 
-void Foam::polyatomicCloud::calculateExternalForce()
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::calculateExternalForce()
 {
-    forAllIter(polyatomicCloud, *this, mol)
+    forAllIter(typename MoleculeCloud<MoleculeType>, *this, mol)
     {
         mol().a() += pot_.gravity();
     }
 }
 
 
-void Foam::polyatomicCloud::removeHighEnergyOverlaps()
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::removeHighEnergyOverlaps()
 {
     Info<< nl << "Removing high energy overlaps, limit = "
         << pot_.potentialEnergyLimit()
@@ -274,23 +275,23 @@ void Foam::polyatomicCloud::removeHighEnergyOverlaps()
 
     // Real-Real interaction
 
-    polyatomic* molI = NULL;
-    polyatomic* molJ = NULL;
+    MoleculeType* molI = NULL;
+    MoleculeType* molJ = NULL;
 
     {
-        DynamicList<polyatomic*> molsToDelete;
+        DynamicList<MoleculeType*> molsToDelete;
 
         const labelListList& dil(il_.dil());
 
         forAll(dil, d)
         {
-            forAll(cellOccupancy_[d],cellIMols)
+            forAll(cellOccupancy_[d], cellIMols)
             {
                 molI = cellOccupancy_[d][cellIMols];
 
                 forAll(dil[d], interactingCells)
                 {
-                    List<polyatomic*> cellJ =
+                    List<MoleculeType*> cellJ =
                         cellOccupancy_[dil[d][interactingCells]];
 
                     forAll(cellJ, cellJMols)
@@ -322,35 +323,35 @@ void Foam::polyatomicCloud::removeHighEnergyOverlaps()
                         }
                     }
                 }
-            }
 
-            forAll(cellOccupancy_[d], cellIOtherMols)
-            {
-                molJ = cellOccupancy_[d][cellIOtherMols];
-
-                if (molJ > molI)
+                forAll(cellOccupancy_[d], cellIOtherMols)
                 {
-                    if (evaluatePotentialLimit(*molI, *molJ))
+                    molJ = cellOccupancy_[d][cellIOtherMols];
+
+                    if (molJ > molI)
                     {
-                        label idI = molI->id();
-
-                        label idJ = molJ->id();
-
-                        if
-                        (
-                            idI == idJ
-                         || findIndex(pot_.removalOrder(), idJ)
-                          < findIndex(pot_.removalOrder(), idI)
-                        )
+                        if (evaluatePotentialLimit(*molI, *molJ))
                         {
-                            if (findIndex(molsToDelete, molJ) == -1)
+                            label idI = molI->id();
+
+                            label idJ = molJ->id();
+
+                            if
+                            (
+                                idI == idJ
+                             || findIndex(pot_.removalOrder(), idJ)
+                              < findIndex(pot_.removalOrder(), idI)
+                            )
                             {
-                                molsToDelete.append(molJ);
+                                if (findIndex(molsToDelete, molJ) == -1)
+                                {
+                                    molsToDelete.append(molJ);
+                                }
+                            }
+                            else if (findIndex(molsToDelete, molI) == -1)
+                            {
+                                molsToDelete.append(molI);
                             }
-                        }
-                        else if (findIndex(molsToDelete, molI) == -1)
-                        {
-                            molsToDelete.append(molI);
                         }
                     }
                 }
@@ -376,19 +377,19 @@ void Foam::polyatomicCloud::removeHighEnergyOverlaps()
     // Real-Referred interaction
 
     {
-        DynamicList<polyatomic*> molsToDelete;
+        DynamicList<MoleculeType*> molsToDelete;
 
         const labelListList& ril(il_.ril());
 
-        List<IDLList<polyatomic> >& referredMols = il_.referredParticles();
+        List<IDLList<MoleculeType> >& referredMols = il_.referredParticles();
 
         forAll(ril, r)
         {
-            IDLList<polyatomic>& refMols = referredMols[r];
+            IDLList<MoleculeType>& refMols = referredMols[r];
 
             forAllIter
             (
-                IDLList<polyatomic>,
+                typename IDLList<MoleculeType>,
                 refMols,
                 refMol
             )
@@ -401,7 +402,7 @@ void Foam::polyatomicCloud::removeHighEnergyOverlaps()
                 {
                     label cellI = realCells[rC];
 
-                    List<polyatomic*> cellIMols = cellOccupancy_[cellI];
+                    List<MoleculeType*> cellIMols = cellOccupancy_[cellI];
 
                     forAll(cellIMols, cIM)
                     {
@@ -430,7 +431,7 @@ void Foam::polyatomicCloud::removeHighEnergyOverlaps()
                              == findIndex(pot_.removalOrder(), idJ)
                             )
                             {
-                                // Remove one of the polyatomics
+                                // Remove one of the molecules
                                 // arbitrarily, assuring that a
                                 // consistent decision is made for
                                 // both real-referred pairs.
@@ -470,17 +471,18 @@ void Foam::polyatomicCloud::removeHighEnergyOverlaps()
         reduce(molsRemoved, sumOp<label>());
     }
 
-    Info<< tab << molsRemoved << " polyatomics removed" << endl;
+    Info<< tab << molsRemoved << " molecules removed" << endl;
 }
 
 
-void Foam::polyatomicCloud::initialiseMolecules
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::initialiseMolecules
 (
-    const IOdictionary& mdInitialiseDict
+    const dictionary& mdInitialiseDict
 )
 {
     Info<< nl
-        << "Initialising polyatomics in each zone specified in "
+        << "Initialising molecules in each zone specified in "
         << mdInitialiseDict.name()
         << endl;
 
@@ -488,7 +490,10 @@ void Foam::polyatomicCloud::initialiseMolecules
 
     if (!cellZones.size())
     {
-        FatalErrorIn("void Foam::polyatomicCloud::initialiseMolecules")
+        FatalErrorIn
+        (
+            "void Foam::MoleculeCloud<MoleculeType>::initialiseMolecules"
+        )
             << "No cellZones found in the mesh."
             << abort(FatalError);
     }
@@ -528,7 +533,10 @@ void Foam::polyatomicCloud::initialiseMolecules
 
                 if (latticeIds.size() != latticePositions.size())
                 {
-                    FatalErrorIn("Foam::polyatomicCloud::initialiseMolecules")
+                    FatalErrorIn
+                    (
+                        "Foam::MoleculeCloud<MoleculeType>::initialiseMolecules"
+                    )
                         << "latticeIds and latticePositions must be the same "
                         << " size." << nl
                         << abort(FatalError);
@@ -550,7 +558,10 @@ void Foam::polyatomicCloud::initialiseMolecules
 
                     if (numberDensity < VSMALL)
                     {
-                        WarningIn("polyatomicCloud::initialiseMolecules")
+                        WarningIn
+                        (
+                            "MoleculeCloud<MoleculeType>::initialiseMolecules"
+                        )
                             << "numberDensity too small, not filling zone "
                             << zone.name() << endl;
 
@@ -571,7 +582,7 @@ void Foam::polyatomicCloud::initialiseMolecules
                     {
                         label id = findIndex(pot_.idList(), latticeIds[i]);
 
-                        const polyatomic::constantProperties& cP
+                        const typename MoleculeType::constantProperties& cP
                         (
                             constProps(id)
                         );
@@ -586,14 +597,16 @@ void Foam::polyatomicCloud::initialiseMolecules
 
                     if (massDensity < VSMALL)
                     {
-                        WarningIn("polyatomicCloud::initialiseMolecules")
+                        WarningIn
+                        (
+                            "MoleculeCloud<MoleculeType>::initialiseMolecules"
+                        )
                             << "massDensity too small, not filling zone "
                             << zone.name() << endl;
 
                         continue;
                     }
 
-
                     latticeCellScale = pow
                     (
                         unitCellMass/(det(latticeCellShape)*massDensity),
@@ -602,7 +615,10 @@ void Foam::polyatomicCloud::initialiseMolecules
                 }
                 else
                 {
-                    FatalErrorIn("Foam::polyatomicCloud::initialiseMolecules")
+                    FatalErrorIn
+                    (
+                        "Foam::MoleculeCloud<MoleculeType>::initialiseMolecules"
+                    )
                         << "massDensity or numberDensity not specified " << nl
                         << abort(FatalError);
                 }
@@ -697,8 +713,8 @@ void Foam::polyatomicCloud::initialiseMolecules
 
                 anchor += (R & (latticeCellShape & latticeAnchor));
 
-                // Continue trying to place polyatomics as long as at
-                // least one polyatomic is placed in each iteration.
+                // Continue trying to place molecule as long as at
+                // least one molecule is placed in each iteration.
                 // The "|| totalZoneMols == 0" condition means that the
                 // algorithm will continue if the origin is outside the
                 // zone.
@@ -719,10 +735,6 @@ void Foam::polyatomicCloud::initialiseMolecules
 
                     bool partOfLayerInBounds = false;
 
-                    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-                    // start of placement of polyatomics
-                    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
                     if (n == 0)
                     {
                         // Special treatment is required for the first position,
@@ -959,10 +971,6 @@ void Foam::polyatomicCloud::initialiseMolecules
                         }
                     }
 
-                    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-                    // end of placement of polyatomics
-                    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
                     if
                     (
                         totalZoneMols == 0
@@ -971,11 +979,12 @@ void Foam::polyatomicCloud::initialiseMolecules
                     {
                         WarningIn
                         (
-                            "Foam::polyatomicCloud::initialiseMolecules()"
+                            "Foam::MoleculeCloud<MoleculeType>::"
+                            "initialiseMolecules()"
                         )
                             << "A whole layer of unit cells was placed "
                             << "outside the bounds of the mesh, but no "
-                            << "polyatomics have been placed in zone '"
+                            << "molecules have been placed in zone '"
                             << zone.name()
                             << "'.  This is likely to be because the zone "
                             << "has few cells ("
@@ -1001,7 +1010,8 @@ void Foam::polyatomicCloud::initialiseMolecules
 }
 
 
-void Foam::polyatomicCloud::createMolecule
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::createMolecule
 (
     const point& position,
     label cell,
@@ -1020,7 +1030,7 @@ void Foam::polyatomicCloud::createMolecule
 
     if (cell == -1)
     {
-        FatalErrorIn("Foam::polyatomicCloud::createMolecule")
+        FatalErrorIn("Foam::MoleculeCloud<MoleculeType>::createMolecule")
             << "Position specified does not correspond to a mesh cell." << nl
             << abort(FatalError);
     }
@@ -1033,59 +1043,31 @@ void Foam::polyatomicCloud::createMolecule
     {
         specialPosition = position;
 
-        special = polyatomic::SPECIAL_TETHERED;
+        special = MoleculeType::SPECIAL_TETHERED;
     }
 
-    const polyatomic::constantProperties& cP(constProps(id));
+    const typename MoleculeType::constantProperties& cP(constProps(id));
 
-    vector v = equipartitionLinearVelocity(temperature, cP.mass());
-
-    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),
-            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)
-        );
-    }
+    typename MoleculeType::trackingData td
+    (
+        *this,
+        MoleculeType::trackingData::tpAccess
+    );
 
     addParticle
     (
-        new polyatomic
+        new MoleculeType
         (
             mesh_,
             position,
             cell,
             tetFace,
             tetPt,
-            Q,
-            v,
-            vector::zero,
-            pi,
-            vector::zero,
+            temperature,
+            bulkVelocity,
             specialPosition,
-            constProps(id),
+            cP,
+            td,
             special,
             id
         )
@@ -1093,11 +1075,12 @@ void Foam::polyatomicCloud::createMolecule
 }
 
 
-Foam::label Foam::polyatomicCloud::nSites() const
+template<class MoleculeType>
+Foam::label Foam::MoleculeCloud<MoleculeType>::nSites() const
 {
     label n = 0;
 
-    forAllConstIter(polyatomicCloud, *this, mol)
+    forAllConstIter(typename MoleculeCloud<MoleculeType>, *this, mol)
     {
         n += constProps(mol().id()).nSites();
     }
@@ -1108,24 +1091,27 @@ Foam::label Foam::polyatomicCloud::nSites() const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::polyatomicCloud::polyatomicCloud
+template<class MoleculeType>
+Foam::MoleculeCloud<MoleculeType>::MoleculeCloud
 (
+    const word& cloudName,
     const polyMesh& mesh,
     const potential& pot,
     bool readFields
 )
 :
-    Cloud<polyatomic>(mesh, "polyatomicCloud", false),
+    Cloud<MoleculeType>(mesh, cloudName, false),
+    moleculeCloud(),
     mesh_(mesh),
     pot_(pot),
     cellOccupancy_(mesh_.nCells()),
     il_(mesh_, pot_.pairPotentials().rCutMax(), false),
     constPropList_(),
-    rndGen_(clock::getTime())
+    rndGen_(label(971501) + 1526*Pstream::myProcNo())
 {
     if (readFields)
     {
-        polyatomic::readFields(*this);
+        MoleculeType::readFields(*this);
     }
 
     buildConstProps();
@@ -1138,27 +1124,30 @@ Foam::polyatomicCloud::polyatomicCloud
 }
 
 
-Foam::polyatomicCloud::polyatomicCloud
+template<class MoleculeType>
+Foam::MoleculeCloud<MoleculeType>::MoleculeCloud
 (
+    const word& cloudName,
     const polyMesh& mesh,
     const potential& pot,
-    const IOdictionary& mdInitialiseDict,
+    const dictionary& mdInitialiseDict,
     bool readFields
 )
 :
-    Cloud<polyatomic>(mesh, "polyatomicCloud", false),
+    Cloud<MoleculeType>(mesh, cloudName, false),
+    moleculeCloud(),
     mesh_(mesh),
     pot_(pot),
-    il_(mesh_, 0.0, false),
+    il_(mesh_),
     constPropList_(),
-    rndGen_(clock::getTime())
+    rndGen_(label(971501) + 1526*Pstream::myProcNo())
 {
     if (readFields)
     {
-        polyatomic::readFields(*this);
+        MoleculeType::readFields(*this);
     }
 
-    clear();
+    this->clear();
 
     buildConstProps();
 
@@ -1168,30 +1157,50 @@ Foam::polyatomicCloud::polyatomicCloud
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::polyatomicCloud::evolve()
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::evolve()
 {
-    polyatomic::trackingData td0(*this, 0);
-    Cloud<polyatomic>::move(td0, mesh_.time().deltaTValue());
+    typename MoleculeType::trackingData td0
+    (
+        *this,
+        MoleculeType::trackingData::tpFirstVelocityHalfStep
+    );
+    Cloud<MoleculeType>::move(td0, mesh_.time().deltaTValue());
 
-    polyatomic::trackingData td1(*this, 1);
-    Cloud<polyatomic>::move(td1, mesh_.time().deltaTValue());
+    typename MoleculeType::trackingData td1
+    (
+        *this,
+        MoleculeType::trackingData::tpLinearTrack
+    );
+    Cloud<MoleculeType>::move(td1, mesh_.time().deltaTValue());
 
-    polyatomic::trackingData td2(*this, 2);
-    Cloud<polyatomic>::move(td2, mesh_.time().deltaTValue());
+    typename MoleculeType::trackingData td2
+    (
+        *this,
+        MoleculeType::trackingData::tpRotationalTrack
+    );
+    Cloud<MoleculeType>::move(td2, mesh_.time().deltaTValue());
 
     calculateForce();
 
-    polyatomic::trackingData td3(*this, 3);
-    Cloud<polyatomic>::move(td3, mesh_.time().deltaTValue());
+    typename MoleculeType::trackingData td3
+    (
+        *this,
+        MoleculeType::trackingData::tpSecondVelocityHalfStep
+    );
+    Cloud<MoleculeType>::move(td3, mesh_.time().deltaTValue());
+
+    info();
 }
 
 
-void Foam::polyatomicCloud::calculateForce()
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::calculateForce()
 {
     buildCellOccupancy();
 
     // Set accumulated quantities to zero
-    forAllIter(polyatomicCloud, *this, mol)
+    forAllIter(typename MoleculeCloud<MoleculeType>, *this, mol)
     {
         mol().siteForces() = vector::zero;
 
@@ -1208,45 +1217,36 @@ void Foam::polyatomicCloud::calculateForce()
 }
 
 
-void Foam::polyatomicCloud::applyConstraintsAndThermostats
-(
-    const scalar targetTemperature,
-    const scalar measuredTemperature
-)
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::info()
 {
-    scalar temperatureCorrectionFactor =
-        sqrt(targetTemperature/max(VSMALL, measuredTemperature));
+    // Calculates and prints the mean momentum and energy in the system
+    // and the number of molecules.
 
-    Info<< "----------------------------------------" << nl
-        << "Temperature equilibration" << nl
-        << "Target temperature = "
-        << targetTemperature << nl
-        << "Measured temperature = "
-        << measuredTemperature << nl
-        << "Temperature correction factor = "
-        << temperatureCorrectionFactor << nl
-        << "----------------------------------------"
-        << endl;
+    typename MoleculeType::trackingData td
+    (
+        *this,
+        MoleculeType::trackingData::tpAccess
+    );
 
-    forAllIter(polyatomicCloud, *this, mol)
-    {
-        mol().v() *= temperatureCorrectionFactor;
-
-        mol().pi() *= temperatureCorrectionFactor;
-    }
+    MoleculeType::info(td);
 }
 
 
-void Foam::polyatomicCloud::writeXYZ(const fileName& fName) const
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::writeXYZ(const fileName& fName) const
 {
     OFstream os(fName);
 
     os  << nSites() << nl
-        << "polyatomicCloud site positions in angstroms" << nl;
+        << "MoleculeCloud<MoleculeType> site positions in angstroms" << nl;
 
-    forAllConstIter(polyatomicCloud, *this, mol)
+    forAllConstIter(typename MoleculeCloud<MoleculeType>, *this, mol)
     {
-        const polyatomic::constantProperties& cP = constProps(mol().id());
+        const typename MoleculeType::constantProperties& cP
+        (
+            constProps(mol().id())
+        );
 
         forAll(mol().sitePositions(), i)
         {
diff --git a/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloud.H b/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.H
similarity index 60%
rename from src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloud.H
rename to src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.H
index 78afbaec2f..38c912a17d 100644
--- a/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloud.H
+++ b/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.H
@@ -22,22 +22,22 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::polyatomicCloud
+    Foam::MoleculeCloud
 
 Description
 
 SourceFiles
-    polyatomicCloudI.H
-    polyatomicCloud.C
+    MoleculeCloudI.H
+    MoleculeCloud.C
 
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef polyatomicCloud_H
-#define polyatomicCloud_H
+#ifndef MoleculeCloud_H
+#define MoleculeCloud_H
 
 #include "Cloud.H"
-#include "polyatomic.H"
+#include "moleculeCloud.H"
 #include "IOdictionary.H"
 #include "potential.H"
 #include "InteractionLists.H"
@@ -51,65 +51,80 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                       Class polyatomicCloud Declaration
+                       Class MoleculeCloud Declaration
 \*---------------------------------------------------------------------------*/
 
-class polyatomicCloud
+template<class MoleculeType>
+class MoleculeCloud
 :
-    public Cloud<polyatomic>
+    public Cloud<MoleculeType>,
+    public moleculeCloud
 {
 
 private:
 
     // Private data
 
+        //-
         const polyMesh& mesh_;
 
+        //-
         const potential& pot_;
 
-        List<DynamicList<polyatomic*> > cellOccupancy_;
+        //-
+        List<DynamicList<MoleculeType*> > cellOccupancy_;
 
-        InteractionLists<polyatomic> il_;
+        //-
+        InteractionLists<MoleculeType> il_;
 
-        List<polyatomic::constantProperties> constPropList_;
+        //-
+        List<typename MoleculeType::constantProperties> constPropList_;
 
+        //-
         Random rndGen_;
 
 
     // Private Member Functions
 
+        //-
         void buildConstProps();
 
+        //-
         void setSiteSizesAndPositions();
 
-        //- Determine which polyatomics are in which cells
+        //- Determine which molecules are in which cells
         void buildCellOccupancy();
 
+        //-
         void calculatePairForce();
 
+        //-
         inline void evaluatePair
         (
-            polyatomic& molI,
-            polyatomic& molJ
+            MoleculeType& molI,
+            MoleculeType& molJ
         );
 
+        //-
         inline bool evaluatePotentialLimit
         (
-            polyatomic& molI,
-            polyatomic& molJ
+            MoleculeType& molI,
+            MoleculeType& molJ
         ) const;
 
+        //-
         void calculateTetherForce();
 
+        //-
         void calculateExternalForce();
 
+        //-
         void removeHighEnergyOverlaps();
 
-        void initialiseMolecules
-        (
-            const IOdictionary& mdInitialiseDict
-        );
+        //-
+        void initialiseMolecules(const dictionary& mdInitialiseDict);
 
+        //-
         void createMolecule
         (
             const point& position,
@@ -122,25 +137,14 @@ private:
             const vector& bulkVelocity
         );
 
+        //-
         label nSites() const;
 
-        inline vector equipartitionLinearVelocity
-        (
-            scalar temperature,
-            scalar mass
-        );
-
-        inline vector equipartitionAngularMomentum
-        (
-            scalar temperature,
-            const polyatomic::constantProperties& cP
-        );
-
         //- Disallow default bitwise copy construct
-        polyatomicCloud(const polyatomicCloud&);
+        MoleculeCloud(const MoleculeCloud&);
 
         //- Disallow default bitwise assignment
-        void operator=(const polyatomicCloud&);
+        void operator=(const MoleculeCloud&);
 
 
 public:
@@ -148,59 +152,81 @@ public:
     // Constructors
 
         //- Construct given mesh and potential references
-        polyatomicCloud
+        MoleculeCloud
         (
+            const word& cloudName,
             const polyMesh& mesh,
             const potential& pot,
             bool readFields = true
         );
 
         //- Construct given mesh, potential and mdInitialiseDict
-        polyatomicCloud
+        MoleculeCloud
         (
+            const word& cloudName,
             const polyMesh& mesh,
             const potential& pot,
-            const IOdictionary& mdInitialiseDict,
+            const dictionary& mdInitialiseDict,
             bool readFields = true
         );
 
 
     // Member Functions
 
-        //- Evolve the polyatomics (move, calculate forces, control state etc)
+        //- Evolve the molecules (move, calculate forces, control state etc)
         void evolve();
 
+        //-
         void calculateForce();
 
-        void applyConstraintsAndThermostats
-        (
-            const scalar targetTemperature,
-            const scalar measuredTemperature
-        );
+        //- Print cloud information
+        void info();
 
 
         // Access
 
+            //-
             inline const polyMesh& mesh() const;
 
+            //-
             inline const potential& pot() const;
 
-            inline const List<DynamicList<polyatomic*> >& cellOccupancy() const;
+            //-
+            inline const List<DynamicList<MoleculeType*> >&
+            cellOccupancy() const;
 
-            inline const InteractionLists<polyatomic>& il() const;
+            //-
+            inline const InteractionLists<MoleculeType>& il() const;
 
-            inline const List<polyatomic::constantProperties>
+            //-
+            inline const List<typename MoleculeType::constantProperties>
             constProps() const;
 
-            inline const polyatomic::constantProperties&
-                constProps(label id) const;
+            //-
+            inline const typename MoleculeType::constantProperties&
+            constProps(label id) const;
 
+            //-
             inline Random& rndGen();
 
+            //-
+            inline vector equipartitionLinearVelocity
+            (
+                scalar temperature,
+                scalar mass
+            );
+
+            //-
+            inline vector equipartitionAngularMomentum
+            (
+                scalar temperature,
+                const typename MoleculeType::constantProperties& cP
+            );
+
 
     // Member Operators
 
-        //- Write polyatomic sites in XYZ format
+        //- Write molecule sites in XYZ format
         void writeXYZ(const fileName& fName) const;
 };
 
@@ -211,7 +237,13 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#include "polyatomicCloudI.H"
+#include "MoleculeCloudI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "MoleculeCloud.C"
+#endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloudI.H b/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloudI.H
similarity index 78%
rename from src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloudI.H
rename to src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloudI.H
index c3207f1554..42c506024f 100644
--- a/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloudI.H
+++ b/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloudI.H
@@ -29,10 +29,11 @@ using namespace Foam::constant;
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-inline void Foam::polyatomicCloud::evaluatePair
+template<class MoleculeType>
+inline void Foam::MoleculeCloud<MoleculeType>::evaluatePair
 (
-    polyatomic& molI,
-    polyatomic& molJ
+    MoleculeType& molI,
+    MoleculeType& molJ
 )
 {
     const pairPotentialList& pairPot = pot_.pairPotentials();
@@ -43,9 +44,15 @@ inline void Foam::polyatomicCloud::evaluatePair
 
     label idJ = molJ.id();
 
-    const polyatomic::constantProperties& constPropI(constProps(idI));
+    const typename MoleculeType::constantProperties& constPropI
+    (
+        constProps(idI)
+    );
 
-    const polyatomic::constantProperties& constPropJ(constProps(idJ));
+    const typename MoleculeType::constantProperties& constPropJ
+    (
+        constProps(idJ)
+    );
 
     forAll(constPropI.pairPotSites(), pI)
     {
@@ -149,10 +156,11 @@ inline void Foam::polyatomicCloud::evaluatePair
 }
 
 
-inline bool Foam::polyatomicCloud::evaluatePotentialLimit
+template<class MoleculeType>
+inline bool Foam::MoleculeCloud<MoleculeType>::evaluatePotentialLimit
 (
-    polyatomic& molI,
-    polyatomic& molJ
+    MoleculeType& molI,
+    MoleculeType& molJ
 ) const
 {
     const pairPotentialList& pairPot = pot_.pairPotentials();
@@ -163,9 +171,15 @@ inline bool Foam::polyatomicCloud::evaluatePotentialLimit
 
     label idJ = molJ.id();
 
-    const polyatomic::constantProperties& constPropI(constProps(idI));
+    const typename MoleculeType::constantProperties& constPropI
+    (
+        constProps(idI)
+    );
 
-    const polyatomic::constantProperties& constPropJ(constProps(idJ));
+    const typename MoleculeType::constantProperties& constPropJ
+    (
+        constProps(idJ)
+    );
 
     forAll(constPropI.pairPotSites(), pI)
     {
@@ -194,14 +208,18 @@ inline bool Foam::polyatomicCloud::evaluatePotentialLimit
 
                 if (rsIsJMag < SMALL)
                 {
-                    WarningIn("polyatomicCloud::removeHighEnergyOverlaps()")
+                    WarningIn
+                    (
+                        "MoleculeCloud<MoleculeType>::"
+                        "removeHighEnergyOverlaps()"
+                    )
                         << "Molecule site pair closer than "
                         << SMALL
                         << ": mag separation = " << rsIsJMag
                         << ". These may have been placed on top of each"
                         << " other by a rounding error in mdInitialise in"
-                        << " parallel or a block filled with polyatomics"
-                        << " twice. Removing one of the polyatomics."
+                        << " parallel or a block filled with moleculess"
+                        << " twice. Removing one of the molecules."
                         << endl;
 
                     return true;
@@ -250,14 +268,18 @@ inline bool Foam::polyatomicCloud::evaluatePotentialLimit
 
                 if (rsIsJMag < SMALL)
                 {
-                    WarningIn("polyatomicCloud::removeHighEnergyOverlaps()")
+                    WarningIn
+                    (
+                        "MoleculeCloud<MoleculeType>::"
+                        "removeHighEnergyOverlaps()"
+                    )
                         << "Molecule site pair closer than "
                         << SMALL
                         << ": mag separation = " << rsIsJMag
                         << ". These may have been placed on top of each"
                         << " other by a rounding error in mdInitialise in"
-                        << " parallel or a block filled with polyatomics"
-                        << " twice. Removing one of the polyatomics."
+                        << " parallel or a block filled with molecules"
+                        << " twice. Removing one of the molecules."
                         << endl;
 
                     return true;
@@ -288,7 +310,64 @@ inline bool Foam::polyatomicCloud::evaluatePotentialLimit
 }
 
 
-inline Foam::vector Foam::polyatomicCloud::equipartitionLinearVelocity
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class MoleculeType>
+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,
     scalar mass
@@ -303,10 +382,12 @@ inline Foam::vector Foam::polyatomicCloud::equipartitionLinearVelocity
 }
 
 
-inline Foam::vector Foam::polyatomicCloud::equipartitionAngularMomentum
+template<class MoleculeType>
+inline Foam::vector
+Foam::MoleculeCloud<MoleculeType>::equipartitionAngularMomentum
 (
     scalar temperature,
-    const polyatomic::constantProperties& cP
+    const typename MoleculeType::constantProperties& cP
 )
 {
     scalar sqrtKbT = sqrt(physicoChemical::k.value()*temperature);
@@ -332,52 +413,4 @@ inline Foam::vector Foam::polyatomicCloud::equipartitionAngularMomentum
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-inline const Foam::polyMesh& Foam::polyatomicCloud::mesh() const
-{
-    return mesh_;
-}
-
-
-inline const Foam::potential& Foam::polyatomicCloud::pot() const
-{
-    return pot_;
-}
-
-
-inline const Foam::List<Foam::DynamicList<Foam::polyatomic*> >&
-    Foam::polyatomicCloud::cellOccupancy() const
-{
-    return cellOccupancy_;
-}
-
-
-inline const Foam::InteractionLists<Foam::polyatomic>&
-    Foam::polyatomicCloud::il() const
-{
-    return il_;
-}
-
-
-inline const Foam::List<Foam::polyatomic::constantProperties>
-    Foam::polyatomicCloud::constProps() const
-{
-    return constPropList_;
-}
-
-
-inline const Foam::polyatomic::constantProperties&
-    Foam::polyatomicCloud::constProps(label id) const
-{
-    return constPropList_[id];
-}
-
-
-inline Foam::Random& Foam::polyatomicCloud::rndGen()
-{
-    return rndGen_;
-}
-
-
 // ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/clouds/baseClasses/moleculeCloud/moleculeCloud.C b/src/lagrangian/molecularDynamics/molecule/clouds/baseClasses/moleculeCloud/moleculeCloud.C
new file mode 100644
index 0000000000..4fda04df7d
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/clouds/baseClasses/moleculeCloud/moleculeCloud.C
@@ -0,0 +1,48 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "moleculeCloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(moleculeCloud, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::moleculeCloud::moleculeCloud()
+{}
+
+
+// * * * * * * * * * * * * * * * *  Destructors  * * * * * * * * * * * * * * //
+
+Foam::moleculeCloud::~moleculeCloud()
+{}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/clouds/baseClasses/moleculeCloud/moleculeCloud.H b/src/lagrangian/molecularDynamics/molecule/clouds/baseClasses/moleculeCloud/moleculeCloud.H
new file mode 100644
index 0000000000..c6ba2096d3
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/clouds/baseClasses/moleculeCloud/moleculeCloud.H
@@ -0,0 +1,83 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::moleculeCloud
+
+Description
+    Virtual abstract base class for templated moleculeCloud
+
+SourceFiles
+    moleculeCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef moleculeCloud_H
+#define moleculeCloud_H
+
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class moleculeCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+class moleculeCloud
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        moleculeCloud(const moleculeCloud&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const moleculeCloud&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("moleculeCloud");
+
+    // Constructors
+
+        //- Null constructor
+        moleculeCloud();
+
+    //- Destructor
+    virtual ~moleculeCloud();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/clouds/derived/monoatomicCloud/monoatomicCloud.H b/src/lagrangian/molecularDynamics/molecule/clouds/derived/monoatomicCloud/monoatomicCloud.H
new file mode 100644
index 0000000000..727bfe226e
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/clouds/derived/polyatomicCloud/polyatomicCloud.H b/src/lagrangian/molecularDynamics/molecule/clouds/derived/polyatomicCloud/polyatomicCloud.H
new file mode 100644
index 0000000000..bb902fe3a4
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/clouds/derived/polyatomicCloud/polyatomicCloud.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::polyatomicCloud
+
+Description
+    Cloud class to simulate polyatomic molecules
+
+SourceFiles
+    polyatomicCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef polyatomicCloud_H
+#define polyatomicCloud_H
+
+#include "MoleculeCloud.H"
+#include "polyatomic.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef MoleculeCloud<polyatomic> polyatomicCloud;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.C b/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.C
new file mode 100644
index 0000000000..f0bad336ee
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.C
@@ -0,0 +1,573 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "controllers.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::controllers::controllers
+(
+    const polyMesh& mesh
+)
+:
+    time_(mesh.time()),
+    controllersDict_
+    (
+        IOobject
+        (
+            "controllersDict",
+            time_.system(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        )
+    ),
+    stateControllersList_(),
+    sCNames_(),
+    sCIds_(),
+    sCFixedPathNames_(),
+    stateControllers_(),
+    fluxControllersList_(),
+    fCNames_(),
+    fCIds_(),
+    fCFixedPathNames_(),
+    fluxControllers_()
+{}
+
+
+Foam::controllers::controllers
+(
+    const polyMesh& mesh,
+    polyatomicCloud& cloud
+)
+:
+    time_(mesh.time()),
+    controllersDict_
+    (
+        IOobject
+        (
+            "controllersDict",
+            time_.system(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    ),
+    stateControllersList_(controllersDict_.lookup("stateControllers")),
+    sCNames_(stateControllersList_.size()),
+    sCIds_(stateControllersList_.size()),
+    sCFixedPathNames_(stateControllersList_.size()),
+    stateControllers_(stateControllersList_.size()),
+    fluxControllersList_(controllersDict_.lookup("fluxControllers")),
+    fCNames_(fluxControllersList_.size()),
+    fCIds_(fluxControllersList_.size()),
+    fCFixedPathNames_(fluxControllersList_.size()),
+    fluxControllers_(fluxControllersList_.size())
+{
+
+    Info << nl << "Creating controllers" << nl << endl;
+
+    // state controllers
+
+    if (!stateControllers_.empty())
+    {
+        forAll(stateControllers_, sC)
+        {
+            const entry& controllersI = stateControllersList_[sC];
+
+            const dictionary& controllersIDict = controllersI.dict();
+
+            stateControllers_[sC] = autoPtr<stateController>
+            (
+                stateController::New(time_, cloud, controllersIDict)
+            );
+
+            sCNames_[sC] = stateControllers_[sC]->type();
+
+            sCIds_[sC] = sC;
+        }
+    }
+
+    //- flux controllers
+
+    if (!fluxControllers_.empty())
+    {
+        forAll(fluxControllers_, fC)
+        {
+            const entry& controllersI = fluxControllersList_[fC];
+
+            const dictionary& controllersIDict = controllersI.dict();
+
+            fluxControllers_[fC] = autoPtr<fluxController>
+            (
+                fluxController::New(time_, cloud, controllersIDict)
+            );
+
+            fCNames_[fC] = fluxControllers_[fC]->type();
+            fCIds_[fC] = fC;
+        }
+    }
+
+    // creating directories for state controllers
+    if (!nStateControllers_.empty())
+    {
+        // case/controllers
+        fileName controllersPath(time_.path()/"controllers");
+
+        if (!isDir(controllersPath))
+        {
+            mkDir(controllersPath);
+        }
+
+        // case/controllers/<cloudName>
+        fileName controllersPath(controllersPath/cloud.name());
+
+        if (!isDir(controllersPath))
+        {
+            mkDir(controllersPath);
+        }
+
+        // case/controllers/<cloudName>/stateControllers
+        fileName stateControllersPath(controllersPath/"stateControllers");
+
+        if (!isDir(stateControllersPath))
+        {
+            mkDir(stateControllersPath);
+        }
+
+        forAll(stateControllers_, sC)
+        {
+            if (stateControllers_[sC]->writeInCase())
+            {
+                // case/controllers/<cloudName>/
+                //     stateControllers/<stateControllerModel>
+                fileName stateControllerPath(stateControllersPath/sCNames_[sC]);
+
+                if (!isDir(stateControllerPath))
+                {
+                    mkDir(stateControllerPath);
+                }
+
+                const word& regionName = stateControllers_[sC]->regionName();
+
+                // case/controllers/<cloudName>/
+                //     stateControllers/<stateControllerModel>/<cellZoneName>
+                fileName zonePath(stateControllerPath/regionName);
+
+                if (!isDir(zonePath))
+                {
+                    mkDir(zonePath);
+                }
+
+                sCFixedPathNames_[sC] = zonePath;
+            }
+        }
+    }
+
+    // creating directories for flux controllers
+    if (nFluxControllers_ > 0)
+    {
+        // case/controllers
+        fileName controllersPath(time_.path()/"controllers");
+
+        if ( !isDir(controllersPath) )
+        {
+            mkDir(controllersPath);
+        }
+
+        // case/controllers/<cloudName>
+        fileName controllersPath(time_.path()/cloud.name());
+
+        if ( !isDir(controllersPath) )
+        {
+            mkDir(controllersPath);
+        }
+
+        // case/controllers/<cloudName>/fluxControllers
+        fileName fluxControllersPath(controllersPath/"fluxControllers");
+
+        if (!isDir(fluxControllersPath))
+        {
+            mkDir(fluxControllersPath);
+        }
+
+        forAll(fluxControllers_, fC)
+        {
+            if (fluxControllers_[fC]->writeInCase())
+            {
+                // case/controllers/<cloudName>/
+                //    fluxControllers/<fluxControllerModel>
+                fileName fluxControllerPath(fluxControllersPath/fCNames_[fC]);
+
+                if (!isDir(fluxControllerPath))
+                {
+                    mkDir(fluxControllerPath);
+                }
+
+                const word& regionName = fluxControllers_[fC]->regionName();
+
+                // case/controllers/<cloudName>/
+                //    fluxControllers/<fluxControllerModel>/<faceZoneName>
+                fileName zonePath(fluxControllerPath/regionName);
+
+                if (!isDir(zonePath))
+                {
+                    mkDir(zonePath);
+                }
+
+                fCFixedPathNames_[fC] = zonePath;
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+controllers::~controllers()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void controllers::initialConfig()
+{
+    forAll(stateControllers_, sC)
+    {
+        stateControllers_[sC]->initialConfiguration();
+    }
+
+    forAll(fluxControllers_, fC)
+    {
+        fluxControllers_[fC]->initialConfiguration();
+    }
+}
+
+
+void controllers::updateTimeInfo()
+{
+    forAll(stateControllers_, sC)
+    {
+        stateControllers_[sC]->updateTime();
+    }
+
+    forAll(fluxControllers_, fC)
+    {
+        fluxControllers_[fC]->updateTime();
+    }
+}
+
+
+void controllers::controlState()
+{
+    forAll(stateControllers_, sC)
+    {
+        stateControllers_[sC]->controlMols();
+    }
+}
+
+
+void controllers::controlVelocitiesI()
+{
+    forAll(stateControllers_, sC)
+    {
+        stateControllers_[sC]->controlMolsBeg();
+    }
+}
+
+void controllers::controlVelocitiesII()
+{
+    forAll(stateControllers_, sC)
+    {
+        stateControllers_[sC]->controlMolsEnd();
+    }
+}
+
+
+void  controllers::controlPriorToForces()
+{
+    forAll(stateControllers_, sC)
+    {
+        stateControllers_[sC]->controlBeforeForces();
+    }
+}
+
+
+void controllers::calculateStateProps()
+{
+    forAll(stateControllers_, sC)
+    {
+        stateControllers_[sC]->calculateProperties();
+    }
+
+    forAll(fluxControllers_, fC)
+    {
+        fluxControllers_[fC]->calculateProperties();
+    }
+}
+
+
+void controllers::outputStateResults()
+{
+    const Time& runTime = time_;
+
+    if (runTime.outputTime())
+    {
+        // creating a set of directories in the current time directory
+        {
+            List<fileName> timePathNames(sCFixedPathNames_.size());
+
+            if (nStateControllers_ > 0)
+            {
+                if (Pstream::master())
+                {
+                    // case/<timeDir>/uniform
+                    fileName uniformTimePath
+                    (
+                        runTime.path()/runTime.timeName()/"uniform"
+                    );
+
+                    if (!isDir(uniformTimePath))
+                    {
+                        mkDir(uniformTimePath);
+                    }
+
+                    if (!stateControllers_.empty())
+                    {
+                        // case/<timeDir>/uniform/controllers
+                        fileName controllersTimePath
+                        (
+                            uniformTimePath/"controllers"
+                        );
+
+                        if (!isDir(controllersTimePath))
+                        {
+                            mkDir(controllersTimePath);
+                        }
+
+                        // case/<timeDir>/uniform/controllers/<cloudName>
+                        fileName cloudTimePath
+                        (
+                            controllersTimePath/cloud.name()
+                        );
+
+                        if (!isDir(cloudTimePath))
+                        {
+                            mkDir(cloudTimePath);
+                        }
+
+                        // case/<timeDir>/uniform/controllers/<cloudName>/
+                        fileName stateControllersTimePath
+                        (
+                            cloudTimePath/"stateControllers"
+                        );
+
+                        if (!isDir(stateControllersTimePath))
+                        {
+                            mkDir(stateControllersTimePath);
+                        }
+
+                        forAll(stateControllers_, sC)
+                        {
+                            if (stateControllers_[sC]->writeInTimeDir())
+                            {
+                                // case/<timeDir>/uniform/controllers/
+                                //     <cloudName>/<stateControllerModel>
+                                fileName sCTimePath
+                                (
+                                    stateControllersTimePath/sCNames_[sC]
+                                );
+
+                                if (!isDir(sCTimePath))
+                                {
+                                    mkDir(sCTimePath);
+                                }
+
+                                // Creating directory for different zones but
+                                // of the same model
+                                const word& regionName =
+                                    stateControllers_[sC]->regionName();
+
+                                // case/<timeDir>/uniform/controllers/
+                                //    <cloudName>/<stateControllerModel>/
+                                //    <cellZoneName>
+                                fileName zoneTimePath(sCTimePath/regionName);
+
+                                if (!isDir(zoneTimePath))
+                                {
+                                    mkDir(zoneTimePath);
+                                }
+
+                                timePathNames[sC] = zoneTimePath;
+                            }
+                        }
+                    }
+                }
+            }
+
+            // write out data
+            forAll(stateControllers_, sC)
+            {
+                stateControllers_[sC]->output
+                (
+                    sCFixedPathNames_[sC],
+                    timePathNames[sC]
+                );
+            }
+        }
+
+        {
+            List<fileName> timePathNames(fCFixedPathNames_.size());
+
+            if (nFluxControllers_ > 0)
+            {
+                if (Pstream::master())
+                {
+                    // case/<timeDir>/uniform
+                    fileName uniformTimePath
+                    (
+                        runTime.path()/runTime.timeName()/"uniform"
+                    );
+
+                    if (!isDir(uniformTimePath))
+                    {
+                        mkDir(uniformTimePath);
+                    }
+
+                    if (!fluxControllers_.empty())
+                    {
+                        // case/<timeDir>/uniform/controllers
+                        fileName controllersTimePath
+                        (
+                            uniformTimePath/"controllers"
+                        );
+
+                        if (!isDir(controllersTimePath))
+                        {
+                            mkDir(controllersTimePath);
+                        }
+
+                        // case/<timeDir>/uniform/controllers/<cloudName>
+                        fileName cloudTimePath
+                        (
+                            controllersTimePath/cloud.name()
+                        );
+
+                        if (!isDir(cloudTimePath))
+                        {
+                            mkDir(cloudTimePath);
+                        }
+
+                        // case/<timeDir>/uniform/fluxControllers
+                        fileName controllersTimePath
+                        (
+                            cloudTimePath/"fluxControllers"
+                        );
+
+                        if (!isDir(controllersTimePath))
+                        {
+                            mkDir(controllersTimePath);
+                        }
+
+                        forAll(fluxControllers_, fC)
+                        {
+                            if (stateControllers_[fC]->writeInTimeDir())
+                            {
+                                // case/<timeDir>/uniform/controllers/
+                                //     <cloudName>/<fluxControllerModel>
+                                fileName fCTimePath
+                                (
+                                    controllersTimePath/fCNames_[fC]
+                                );
+
+                                if (!isDir(fCTimePath))
+                                {
+                                    mkDir(fCTimePath);
+                                }
+
+                                const word& regionName =
+                                    fluxControllers_[fC]->regionName();
+
+                                // case/<timeDir>/uniform/controllers/
+                                //     <cloudName>/<fluxControllerModel>/
+                                //     <faceZoneName>
+                                fileName zoneTimePath(fCTimePath/regionName);
+
+                                if (!isDir(zoneTimePath))
+                                {
+                                    mkDir(zoneTimePath);
+                                }
+
+                                timePathNames[fC] = zoneTimePath;
+                            }
+                        }
+                    }
+                }
+            }
+
+            // write out data
+            forAll(fluxControllers_, fC)
+            {
+                fluxControllers_[fC]->output
+                (
+                    fCFixedPathNames_[fC],
+                    timePathNames[fC]
+                );
+            }
+        }
+
+        // Re-read dictionaries for modified properties (run-time selection)
+        {
+            stateControllersList_.clear();
+
+            stateControllersList_ = controllersDict_.lookup("stateControllers");
+
+            forAll(stateControllers_, sC)
+            {
+                const entry& controllersI = stateControllersList_[sC];
+                const dictionary& controllersIDict = controllersI.dict();
+
+                stateControllers_[sC]->updateProperties(controllersIDict);
+            }
+        }
+
+        {
+            fluxControllersList_.clear();
+
+            fluxControllersList_ = controllersDict_.lookup("fluxControllers");
+
+            forAll(fluxControllers_, fC)
+            {
+                const entry& controllersI = fluxControllersList_[fC];
+                const dictionary& controllersIDict = controllersI.dict();
+
+                fluxControllers_[fC]->updateProperties(controllersIDict);
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.H b/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.H
new file mode 100644
index 0000000000..22e9ca7cca
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.H
@@ -0,0 +1,163 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+    controllers
+
+Description
+
+    Stores all the information for the controllers models defined within
+    the controllersDict, and selects & builds the models automatically.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef controllers_H
+#define controllers_H
+
+#include "List.H"
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "polyMesh.H"
+#include "timeData.H"
+#include "stateController.H"
+#include "fluxController.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class controllers Declaration
+\*---------------------------------------------------------------------------*/
+
+class controllers
+{
+    // Private data
+
+        Time& time_;
+
+        //- The entire dictionary (containing multiple subDictionaries)
+        IOdictionary controllersDict_;
+
+        //- state controllers
+        PtrList<entry> stateControllersList_;
+        List<word> sCNames_;
+        List<label> sCIds_;
+        List<fileName> sCFixedPathNames_;
+        List< autoPtr<stateController> > stateControllers_;
+
+        //- flux controllers
+        PtrList<entry> fluxControllersList_;
+        List<word> fCNames_;
+        List<label> fCIds_;
+        List<fileName> fCFixedPathNames_;
+        List< autoPtr<fluxController> > fluxControllers_;
+
+
+public:
+
+    // Constructors
+
+        //- Null Constructor
+        controllers
+        (
+            const polyMesh& mesh
+        );
+
+        //- Constructor for with cloud
+        controllers
+        (
+            const polyMesh& mesh,
+            polyatomicCloud& cloud
+        );
+
+    //- Destructor
+    ~controllers();
+
+    // Member Functions
+
+        //- Initial configuration call this function after the polyatomicCloud
+        //  is completely initialised
+        void initialConfig();
+
+        //- this function is to be called at the beginning of the MD time-step.
+        //  since we have placed a non-referenced time-data class in the
+        //  state-controller class.
+        void updateTimeInfo();
+
+        //- control molecular state -- call this after the intermolecular force
+        //  calulation
+        void controlState();
+
+        //-
+        void controlVelocitiesI();
+
+        //-
+        void controlVelocitiesII();
+
+        //-
+        void controlPriorToForces();
+
+        //- calculate properties -- call this at the end of the MD time-step.
+        void calculateStateProps();
+
+        //- output -- call this function at the end of the MD time-step
+        void outputStateResults();
+
+        // Access
+
+            //-
+            inline List< autoPtr<stateController> >& stateControllers();
+
+            //-
+            inline const List< autoPtr<stateController> >&
+            stateControllers() const;
+
+            //-
+            inline List< autoPtr<fluxController> >& fluxControllers();
+
+            //-
+            inline const List< autoPtr<fluxController> >&
+            fluxControllers() const;
+
+            //-
+            inline const List<word>& stateControllersNames() const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "controllersI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllersI.H b/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllersI.H
new file mode 100644
index 0000000000..66c10f0784
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllersI.H
@@ -0,0 +1,64 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::List<Foam::autoPtr<Foam::stateController> >&
+Foam::controllers::stateControllers()
+{
+    return stateControllers_;
+}
+
+
+const Foam::List<Foam::autoPtr<Foam::stateController> >&
+Foam::controllers::stateControllers() const
+{
+    return stateControllers_;
+}
+
+
+Foam::List<Foam::autoPtr<Foam::fluxController> >&
+Foam::controllers::fluxControllers()
+{
+    return fluxControllers_;
+}
+
+
+const Foam::List< autoPtr<fluxController> >&
+Foam::controllers::fluxControllers() const
+{
+    return fluxControllers_;
+}
+
+
+const Foam::List<Foam::word>& Foam::controllers::stateControllersNames() const
+{
+    return sCNames_;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.C b/src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.C
new file mode 100644
index 0000000000..711b661758
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.C
@@ -0,0 +1,524 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "waterFluxController.H"
+#include "IFstream.H"
+#include "graph.H"
+#include "polyatomicCloud.H"
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(waterFluxController, 0);
+
+defineRunTimeSelectionTable(waterFluxController, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+waterFluxController::waterFluxController
+(
+    Time& t,
+    polyatomicCloud& cloud,
+    const dictionary& dict
+)
+:
+    mesh_(refCast<const fvMesh>(cloud.mesh())),
+    cloud_(cloud),
+    rndGen_(clock::getTime()),
+    controllerDict_(dict.subDict("controllerProperties")),
+    timeDict_(controllerDict_.subDict("timeProperties")),
+    time_(t, timeDict_),
+    regionName_(controllerDict_.lookup("zoneName")),
+    regionId_(-1),
+    zoneSurfaceArea_(0.0),
+    internalFaces_(),
+    processorFaces_(),
+    control_(true),
+    readStateFromFile_(true),
+    singleValueController_(false),
+    density_(0.0),
+    velocity_(vector::zero),
+    temperature_(0.0),
+    pressure_(0.0),
+    strainRate_(tensor::zero),
+    tempGradient_(vector::zero),
+    fieldController_(false),
+    densities_(),
+    velocities_(),
+    temperatures_(),
+    pressures_(),
+    writeInTimeDir_(true),
+    writeInCase_(true)
+{
+    const faceZoneMesh& faceZones = mesh_.faceZones();
+    regionId_ = faceZones.findZoneID(regionName_);
+
+    if (regionId_ == -1)
+    {
+        FatalErrorIn("waterFluxController::waterFluxController()")
+        << "Cannot find region (faceZone): " << regionName_ << nl << "in: "
+        << t.time().system()/"controllersDict"
+        << exit(FatalError);
+    }
+
+    control_ = Switch(controllerDict_.lookup("controlSwitch"));
+    readStateFromFile_ = Switch(controllerDict_.lookup("readStateFromFile"));
+
+    setFacesInfo();
+}
+
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+autoPtr<waterFluxController> waterFluxController::New
+(
+    Time& t,
+    polyatomicCloud& cloud,
+    const dictionary& dict
+)
+{
+    word waterFluxControllerName
+    (
+        dict.lookup("fluxControllerModel")
+    );
+
+    Info<< "Selecting fluxController "
+         << waterFluxControllerName << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(waterFluxControllerName);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalError
+            << "waterFluxController::New(const dictionary&) : " << endl
+            << "    unknown waterFluxController type "
+            << waterFluxControllerName
+            << ", constructor not in hash table" << endl << endl
+            << "    Valid injector types are :" << endl;
+        Info<< dictionaryConstructorTablePtr_->toc() << abort(FatalError);
+    }
+
+    return autoPtr<waterFluxController>
+        (
+                cstrIter()(t, cloud, dict)
+        );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+waterFluxController::~waterFluxController()
+{}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+// void waterFluxController::updateTime()
+// {
+//     time_++;
+//
+//     const scalar& t = time_.time().timeOutputValue();
+//
+//     if ((t - initialTime_) < timePeriod_)
+//     {
+//         time_.controlTimeInterval().endTime() = false;
+// //         control_ = false;
+//     }
+//     else
+//     {
+// //         control_ = true;
+//     }
+// }
+
+
+void waterFluxController::setFacesInfo()
+{
+    const labelList& faces = controlZone();
+
+    if (Pstream::parRun())
+    {
+        DynamicList<label> processorFaces(0);
+
+        forAll(mesh_.boundaryMesh(), patchI)
+        {
+            const polyPatch& patch = mesh_.boundaryMesh()[patchI];
+
+            if (isA<processorPolyPatch>(patch))
+            {
+                for (label p = 0; p < patch.size(); p++)
+                {
+                    label patchFaceI = p + patch.start();
+                    label faceId = findIndex (faces, patchFaceI);
+
+                    if (faceId != -1)
+                    {
+                        processorFaces.append(patchFaceI);
+                    }
+                }
+            }
+        }
+
+        processorFaces.shrink();
+
+        processorFaces_.setSize(processorFaces.size(), -1);
+
+        forAll(processorFaces, f)
+        {
+            processorFaces_[f] = processorFaces[f];
+        }
+
+        label nInternalFaces = faces.size() - processorFaces.size();
+        internalFaces_.setSize(nInternalFaces, -1);
+
+        label counter = 0;
+
+        forAll(faces, f)
+        {
+            const label& faceI = faces[f];
+
+            if (findIndex(processorFaces, faceI) == -1)
+            {
+                internalFaces_[counter] = faceI;
+                counter++;
+            }
+        }
+
+//         Pout << "processorFaces: " << processorFaces_ << endl;
+//         Pout << "internalFaces: " << internalFaces_ << endl;
+
+        forAll(internalFaces_, f)
+        {
+            const label& faceI = internalFaces_[f];
+            zoneSurfaceArea_ += mag(mesh_.faceAreas()[faceI]);
+        }
+
+        // faces on a zone located on a processor cut belong to both processors
+        // (hence the 0.5)
+
+        forAll(processorFaces_, f)
+        {
+            const label& faceI = processorFaces_[f];
+            zoneSurfaceArea_ += 0.5*mag(mesh_.faceAreas()[faceI]);
+        }
+
+
+        if (Pstream::parRun())
+        {
+            for (int p = 0; p < Pstream::nProcs(); p++)
+            {
+                if (p != Pstream::myProcNo())
+                {
+                    const int proc = p;
+                    {
+                        OPstream toNeighbour(Pstream::blocking, proc);
+                        toNeighbour << zoneSurfaceArea_;
+                    }
+                }
+            }
+
+            //- receiving
+            for (int p = 0; p < Pstream::nProcs(); p++)
+            {
+                if (p != Pstream::myProcNo())
+                {
+                    scalar zoneSurfaceAreaProc;
+
+                    const int proc = p;
+                    {
+                        IPstream fromNeighbour(Pstream::blocking, proc);
+                        fromNeighbour >> zoneSurfaceAreaProc;
+                    }
+
+                    zoneSurfaceArea_ += zoneSurfaceAreaProc;
+                }
+            }
+        }
+    }
+    else
+    {
+        forAll(faces, f)
+        {
+            const label& faceI = faces[f];
+
+            zoneSurfaceArea_ += mag(mesh_.faceAreas()[faceI]);
+        }
+    }
+}
+
+
+
+void waterFluxController::updateTime()
+{
+    time_++;
+
+//     const scalar& t = time_.time().timeOutputValue();
+//
+//     if ((t - initialTime_) < timePeriod_)
+//     {
+//         time_.controlTimeInterval().endTime() = false;
+// //         control_ = false;
+//     }
+//     else
+//     {
+// //         control_ = true;
+//     }
+}
+
+void waterFluxController::updateFluxControllerProperties
+(
+    const dictionary& newDict
+)
+{
+    controllerDict_ = newDict.subDict("controllerProperties");
+
+    //- you can reset the controlling zone from here. This essentially
+    //  means that the coupling zone can infact move arbitrarily. To make
+    //  this happen we probably need to devise a technique for automatically
+    //  changing the cellZone else where, and then calling this function to
+    //  reset the controlling zone in which the controller operates in.
+
+    if (controllerDict_.found("controlSwitch"))
+    {
+        control_ = Switch(controllerDict_.lookup("controlSwitch"));
+    }
+
+    if (controllerDict_.found("readStateFromFile"))
+    {
+        readStateFromFile_ = Switch
+        (
+            controllerDict_.lookup("readStateFromFile")
+        );
+    }
+}
+
+const labelList& waterFluxController::controlZone() const
+{
+    return mesh_.faceZones()[regionId_];
+}
+
+label waterFluxController::isFaceOnControlZone(const label& faceI)
+{
+    const label f = findIndex(controlZone(), faceI);
+
+    return f;
+}
+
+const word& waterFluxController::regionName() const
+{
+    return regionName_;
+}
+
+const scalar& waterFluxController::density() const
+{
+    return density_;
+}
+
+scalar& waterFluxController::density()
+{
+    return density_;
+}
+
+const vector& waterFluxController::velocity() const
+{
+    return velocity_;
+}
+
+vector& waterFluxController::velocity()
+{
+    return velocity_;
+}
+
+const scalar& waterFluxController::temperature() const
+{
+    return temperature_;
+}
+
+scalar& waterFluxController::temperature()
+{
+    return temperature_;
+}
+
+const scalar& waterFluxController::pressure() const
+{
+    return pressure_;
+}
+
+scalar& waterFluxController::pressure()
+{
+    return pressure_;
+}
+
+const tensor& waterFluxController::strainRate() const
+{
+    return strainRate_;
+}
+
+tensor& waterFluxController::strainRate()
+{
+    return strainRate_;
+}
+
+const vector& waterFluxController::tempGradient() const
+{
+    return tempGradient_;
+}
+
+vector& waterFluxController::tempGradient()
+{
+    return tempGradient_;
+}
+
+
+const scalarField& waterFluxController::densityField() const
+{
+    return densities_;
+}
+
+scalarField& waterFluxController::densityField()
+{
+    return densities_;
+}
+
+const vectorField& waterFluxController::velocityField() const
+{
+    return velocities_;
+}
+vectorField& waterFluxController::velocityField()
+{
+    return velocities_;
+}
+
+const scalarField& waterFluxController::temperatureField() const
+{
+    return temperatures_;
+}
+
+scalarField& waterFluxController::temperatureField()
+{
+    return temperatures_;
+}
+
+const scalarField& waterFluxController::pressureField() const
+{
+    return pressures_;
+}
+
+scalarField& waterFluxController::pressureField()
+{
+    return pressures_;
+}
+
+
+const bool& waterFluxController::singleValueController() const
+{
+    return singleValueController_;
+}
+
+bool& waterFluxController::singleValueController()
+{
+    return singleValueController_;
+}
+
+const bool& waterFluxController::fieldController() const
+{
+    return fieldController_;
+}
+
+bool& waterFluxController::fieldController()
+{
+    return fieldController_;
+}
+
+
+const bool& waterFluxController::writeInTimeDir() const
+{
+    return writeInTimeDir_;
+}
+
+const bool& waterFluxController::writeInCase() const
+{
+    return writeInCase_;
+}
+
+
+// const scalar waterFluxController::avReqDensity() const
+// {
+//     scalar totalDensity = 0.0;
+//
+//     forAll(densities_, c)
+//     {
+//         totalDensity += densities_[c];
+//     }
+//
+//     if (cells_.size() > 0)
+//     {
+//         totalDensity /= scalar(cells_.size());
+//     }
+//
+//     return totalDensity;
+// }
+//
+// const vector waterFluxController::avReqVelocity() const
+// {
+//     vector totalVel = vector::zero;
+//
+//     forAll(velocities_, c)
+//     {
+//         totalVel += velocities_[c];
+//     }
+//
+//     if (cells_.size() > 0)
+//     {
+//         totalVel /= scalar(cells_.size());
+//     }
+//
+//     return totalVel;
+// }
+//
+// const scalar waterFluxController::avReqTemperature() const
+// {
+//     scalar totalTemp = 0.0;
+//
+//     forAll(densities_, c)
+//     {
+//         totalTemp += temperatures_[c];
+//     }
+//
+//     if (cells_.size() > 0)
+//     {
+//         totalTemp /= scalar(cells_.size());
+//     }
+//
+//     return totalTemp;
+// }
+
+
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.H b/src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.H
new file mode 100644
index 0000000000..9b8a375277
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.H
@@ -0,0 +1,272 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    waterFluxController
+
+Description
+
+SourceFiles
+    waterFluxControllerI.H
+    waterFluxController.C
+    waterFluxControllerIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef waterFluxController_H
+#define waterFluxController_H
+
+#include "IOdictionary.H"
+#include "Time.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "vector.H"
+#include "volFields.H"
+#include "Random.H"
+#include "polyatomic.H"
+#include "timeData.H"
+#include "writeTimeData.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class waterFluxController Declaration
+\*---------------------------------------------------------------------------*/
+
+class waterFluxController
+{
+
+protected:
+
+    // Protected data
+
+//         Time& time_;
+
+        const fvMesh& mesh_;
+
+        polyatomicCloud& cloud_;
+
+        Random rndGen_;
+
+        //- subDictionary containing the properties
+        dictionary controllerDict_;
+
+
+        dictionary timeDict_;
+
+        timeData time_;
+
+        //- name of face zone
+        word regionName_;
+        label regionId_;
+//         labelList faces_;
+
+        scalar zoneSurfaceArea_;
+
+        labelList internalFaces_;
+        labelList processorFaces_;
+
+        bool control_;
+        bool readStateFromFile_;
+
+        //- set all the properties below from model if required
+
+            bool singleValueController_;
+
+            // target values
+            scalar density_;
+            vector velocity_;
+            scalar temperature_;
+            scalar pressure_;
+
+            tensor strainRate_;
+            vector tempGradient_;
+
+            bool fieldController_;
+
+            //- targeted fields
+            scalarField densities_;
+            vectorField velocities_;
+            scalarField temperatures_;
+            scalarField pressures_;
+
+            bool writeInTimeDir_;
+            bool writeInCase_;
+
+
+
+
+    // Private Member Functions
+
+        void setFacesInfo();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("waterFluxController");
+
+    // Declare runtime constructor selection table
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            waterFluxController,
+            dictionary,
+            (
+                Time& t,
+                polyatomicCloud& cloud,
+                const dictionary& dict
+            ),
+            (t, cloud, dict)
+        );
+
+    // Constructors
+
+        //- Construct from components
+        waterFluxController
+        (
+                        Time& t,
+            polyatomicCloud& cloud,
+            const dictionary& dict
+        );
+
+
+    // Selectors
+
+        static autoPtr<waterFluxController> New
+        (
+                        Time& t,
+            polyatomicCloud& cloud,
+            const dictionary& dict
+        );
+
+
+    // Destructor
+
+        virtual ~waterFluxController();
+
+
+    // Member Functions
+
+        void updateTime();
+
+        //- create an initial configuration
+        virtual void initialConfiguration() = 0;
+
+
+        //- calculate any required properties
+        virtual void calculateProperties() = 0;
+
+
+        //- control the polyatomic from the tracking function
+        virtual void controlMol
+        (
+            polyatomic& mol,
+            polyatomic::trackData& td
+        ) = 0;
+
+        //- output data
+        virtual void output
+        (
+            const fileName& fixedPathName,
+            const fileName& timePath
+        ) = 0;
+
+
+        //- E. update properties from a modified dictionary
+        virtual void updateProperties(const dictionary&) = 0;
+
+
+        void updateFluxControllerProperties(const dictionary&);
+
+        // Access
+
+            //- return the control zone cells
+            const labelList& controlZone() const;
+
+
+            label isFaceOnControlZone(const label& faceI);
+
+            //- return the control zone name
+            const word& regionName() const;
+
+           //- return the targeted values
+            const scalar& density() const;
+            scalar& density();
+
+            const vector& velocity() const;
+            vector& velocity();
+
+            const scalar& temperature() const;
+            scalar& temperature();
+
+            const scalar& pressure() const;
+            scalar& pressure();
+
+            const tensor& strainRate() const;
+            tensor& strainRate();
+
+            const vector& tempGradient() const;
+            vector& tempGradient();
+
+            //- return the targeted fields
+            const scalarField& densityField() const;
+            scalarField& densityField();
+
+            const vectorField& velocityField() const;
+            vectorField& velocityField();
+
+            const scalarField& temperatureField() const;
+            scalarField& temperatureField();
+
+            const scalarField& pressureField() const;
+            scalarField& pressureField();
+
+
+            const bool& singleValueController() const;
+            bool& singleValueController();
+
+            const bool& fieldController() const;
+            bool& fieldController();
+
+            const bool& writeInTimeDir() const;
+            const bool& writeInCase() const;
+
+//             const scalar avReqDensity() const;
+//             const vector avReqVelocity() const;
+//             const scalar avReqTemperature() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.C b/src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.C
new file mode 100644
index 0000000000..b98aa126a3
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.C
@@ -0,0 +1,490 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "stateController.H"
+#include "IFstream.H"
+#include "polyatomicCloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(stateController, 0);
+
+defineRunTimeSelectionTable(stateController, dictionary);
+
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::stateController::stateController
+(
+    polyatomicCloud& cloud,
+    const dictionary& dict
+)
+:
+    mesh_(refCast<const fvMesh>(cloud.mesh())),
+    cloud_(cloud),
+    rndGen_(clock::getTime()),
+    controllerDict_(dict.subDict("controllerProperties")),
+    timeDict_(controllerDict_.subDict("timeProperties")),
+    time_(mesh_.time(), timeDict_),
+    timePeriod_(readScalar(timeDict_.lookup("initialTimePeriod"))), //temp
+    initialTime_(time_.time().startTime().value()),
+    regionName_(controllerDict_.lookup("zoneName")),
+    regionId_(-1),
+    control_(true),
+    readStateFromFile_(true),
+    singleValueController_(false),
+    density_(0.0),
+    velocity_(vector::zero),
+    temperature_(0.0),
+    pressure_(0.0),
+    strainRate_(tensor::zero),
+    tempGradient_(vector::zero),
+    fieldController_(false),
+    densities_(),
+    velocities_(),
+    temperatures_(),
+    pressures_(),
+    writeInTimeDir_(true),
+    writeInCase_(true)
+{
+    const cellZoneMesh& cellZones = mesh_.cellZones();
+
+    regionId_ = cellZones.findZoneID(regionName_);
+
+    if (regionId_ == -1)
+    {
+        FatalErrorIn("stateController::stateController()")
+            << "Cannot find region: " << regionName_ << nl << "in: "
+            << time_.time().system()/"controllersDict"
+            << exit(FatalError);
+    }
+
+    control_ = Switch(controllerDict_.lookup("controlSwitch"));
+
+    readStateFromFile_ = Switch(controllerDict_.lookup("readStateFromFile"));
+
+    const scalar& avTimeInterval = time_.averageTimeInterval().deltaT();
+
+    if ((timePeriod_ < avTimeInterval) && (timePeriod_ > 0.0))
+    {
+        timePeriod_ = avTimeInterval;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::stateController> Foam::stateController::New
+(
+    polyatomicCloud& cloud,
+    const dictionary& dict
+)
+{
+    word stateControllerName
+    (
+        dict.lookup("stateControllerModel")
+    );
+
+    Info<< "Selecting stateController "
+         << stateControllerName << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(stateControllerName);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalError
+            << "stateController::New(const dictionary&) : " << endl
+            << "    unknown stateController type "
+            << stateControllerName
+            << ", constructor not in hash table" << endl << endl
+            << "    Valid types are :" << endl;
+        Info<< dictionaryConstructorTablePtr_->toc() << abort(FatalError);
+    }
+
+    return autoPtr<stateController>
+    (
+        cstrIter()(cloud, dict)
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::stateController::~stateController()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::stateController::updateTime()
+{
+    time_++;
+
+    const scalar& t = time_.time().timeOutputValue();
+
+    if ((t - initialTime_) < timePeriod_)
+    {
+        time_.controlTimeInterval().endTime() = false;
+    }
+}
+
+
+void Foam::stateController::updateStateControllerProperties
+(
+    const dictionary& newDict
+)
+{
+    controllerDict_ = newDict.subDict("controllerProperties");
+
+    if (controllerDict_.found("controlSwitch"))
+    {
+        control_ = Switch(controllerDict_.lookup("controlSwitch"));
+    }
+
+    if (controllerDict_.found("readStateFromFile"))
+    {
+        readStateFromFile_ = Switch
+        (
+            controllerDict_.lookup("readStateFromFile")
+        );
+    }
+
+    timeDict_ = controllerDict_.subDict("timeProperties");
+
+    if (timeDict_.found("resetAtOutput"))
+    {
+        time_.resetFieldsAtOutput() = Switch(timeDict_.lookup("resetAtOutput"));
+    }
+}
+
+
+const Foam::labelList& Foam::stateController::controlZone() const
+{
+    return mesh_.cellZones()[regionId_];
+}
+
+const Foam::word& Foam::stateController::regionName() const
+{
+    return regionName_;
+}
+
+
+Foam::scalar Foam::stateController::density() const
+{
+    return density_;
+}
+
+
+Foam::scalar& Foam::stateController::density()
+{
+    return density_;
+}
+
+
+const Foam::vector& Foam::stateController::velocity() const
+{
+    return velocity_;
+}
+
+
+Foam::vector& Foam::stateController::velocity()
+{
+    return velocity_;
+}
+
+
+Foam::scalar Foam::stateController::temperature() const
+{
+    return temperature_;
+}
+
+
+Foam::scalar& Foam::stateController::temperature()
+{
+    return temperature_;
+}
+
+const Foam::scalar& Foam::stateController::pressure() const
+{
+    return pressure_;
+}
+
+
+Foam::scalar& Foam::stateController::pressure()
+{
+    return pressure_;
+}
+
+
+const Foam::tensor& Foam::stateController::strainRate() const
+{
+    return strainRate_;
+}
+
+
+Foam::tensor& Foam::stateController::strainRate()
+{
+    return strainRate_;
+}
+
+
+const Foam::vector& Foam::stateController::tempGradient() const
+{
+    return tempGradient_;
+}
+
+
+Foam::vector& Foam::stateController::tempGradient()
+{
+    return tempGradient_;
+}
+
+
+const Foam::scalarField& Foam::stateController::densityField() const
+{
+    return densities_;
+}
+
+Foam::scalarField& Foam::stateController::densityField()
+{
+    return densities_;
+}
+
+
+const Foam::vectorField& Foam::stateController::velocityField() const
+{
+    return velocities_;
+}
+
+
+Foam::vectorField& Foam::stateController::velocityField()
+{
+    return velocities_;
+}
+
+
+const Foam::scalarField& Foam::stateController::temperatureField() const
+{
+    return temperatures_;
+}
+
+
+Foam::scalarField& Foam::stateController::temperatureField()
+{
+    return temperatures_;
+}
+
+
+const Foam::scalarField& Foam::stateController::pressureField() const
+{
+    return pressures_;
+}
+
+
+Foam::scalarField& Foam::stateController::pressureField()
+{
+    return pressures_;
+}
+
+
+bool Foam::stateController::singleValueController() const
+{
+    return singleValueController_;
+}
+
+
+bool& Foam::stateController::singleValueController()
+{
+    return singleValueController_;
+}
+
+
+bool Foam::stateController::fieldController() const
+{
+    return fieldController_;
+}
+
+
+bool& Foam::stateController::fieldController()
+{
+    return fieldController_;
+}
+
+
+bool Foam::stateController::writeInTimeDir() const
+{
+    return writeInTimeDir_;
+}
+
+
+bool Foam::stateController::writeInCase() const
+{
+    return writeInCase_;
+}
+
+
+Foam::scalar Foam::stateController::avReqDensity()
+{
+    scalar totalDensity = 0.0;
+
+    if (singleValueController_)
+    {
+        totalDensity = density_;
+    }
+    else if (fieldController_)
+    {
+        label controlCells = controlZone().size();
+
+        forAll(densities_, c)
+        {
+            totalDensity += densities_[c];
+        }
+
+        if (Pstream::parRun())
+        {
+            reduce(totalDensity, sumOp<scalar>());
+
+            reduce(controlCells, sumOp<label>());
+        }
+
+        if (controlCells > 0)
+        {
+            totalDensity /= scalar(controlCells);
+        }
+    }
+
+    return totalDensity;
+}
+
+
+Foam::vector Foam::stateController::avReqVelocity()
+{
+    vector totalVel = vector::zero;
+
+    if (singleValueController_)
+    {
+        totalVel = velocity_;
+    }
+    else if (fieldController_)
+    {
+        label controlCells = controlZone().size();
+
+        forAll(velocities_, c)
+        {
+            totalVel += velocities_[c];
+        }
+
+        if (Pstream::parRun())
+        {
+            reduce(totalVel, sumOp<vector>());
+
+            reduce(controlCells, sumOp<label>());
+        }
+
+        if (controlCells > 0)
+        {
+            totalVel /= scalar(controlCells);
+        }
+    }
+
+    return totalVel;
+}
+
+
+Foam::scalar Foam::stateController::avReqTemperature()
+{
+    scalar totalTemp = 0.0;
+
+    if (singleValueController_)
+    {
+        totalTemp = temperature_;
+    }
+    else if (fieldController_)
+    {
+        label controlCells = controlZone().size();
+
+        forAll(temperatures_, c)
+        {
+            totalTemp += temperatures_[c];
+        }
+
+        if (Pstream::parRun())
+        {
+            reduce(totalTemp, sumOp<scalar>());
+
+            reduce(controlCells, sumOp<label>());
+        }
+
+        if (controlCells > 0)
+        {
+            totalTemp /= scalar(controlCells);
+        }
+    }
+
+    return totalTemp;
+}
+
+
+Foam::scalar Foam::stateController::avReqPressure()
+{
+    scalar totalPressure = 0.0;
+
+    if (singleValueController_)
+    {
+        totalPressure = pressure_;
+    }
+    else if (fieldController_)
+    {
+        label controlCells = controlZone().size();
+
+        forAll(pressures_, c)
+        {
+            totalPressure += pressures_[c];
+        }
+
+        if (Pstream::parRun())
+        {
+            reduce(totalPressure, sumOp<scalar>());
+
+            reduce(controlCells, sumOp<label>());
+        }
+
+        if (controlCells > 0)
+        {
+            totalPressure /= scalar(controlCells);
+        }
+    }
+
+    return totalPressure;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.H b/src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.H
new file mode 100644
index 0000000000..3e142c3962
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.H
@@ -0,0 +1,270 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    stateController
+
+Description
+
+    Basic/abstract class of a state controller
+
+SourceFiles
+    stateControllerI.H
+    stateController.C
+    stateControllerIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef stateController_H
+#define stateController_H
+
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "vector.H"
+#include "volFields.H"
+#include "Random.H"
+#include "polyatomic.H"
+#include "timeData.H"
+#include "writeTimeData.H"
+#include "selectIds.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class stateController Declaration
+\*---------------------------------------------------------------------------*/
+
+class stateController
+{
+
+protected:
+
+    // Protected data
+
+        //-
+        const fvMesh& mesh_;
+
+        //-
+        polyatomicCloud& cloud_;
+
+        //-
+        Random rndGen_;
+
+        //- subDictionary containing the properties
+        dictionary controllerDict_;
+
+        //-
+        dictionary timeDict_;
+
+        //-
+        timeData time_;
+
+        //-
+        scalar timePeriod_;
+
+        //-
+        scalar initialTime_;
+
+        //- name of control zone
+        word regionName_;
+
+        //-
+        label regionId_;
+
+        //-
+        bool control_;
+
+        //-
+        bool readStateFromFile_;
+
+        //- set all the properties below from model if required
+
+            //-
+            bool singleValueController_;
+
+            //- target values
+            scalar density_;
+            vector velocity_;
+            scalar temperature_;
+            scalar pressure_;
+            tensor strainRate_;
+            vector tempGradient_;
+
+            //- set this in model
+            bool fieldController_;
+
+            //- targeted fields
+            scalarField densities_;
+            vectorField velocities_;
+            scalarField temperatures_;
+            scalarField pressures_;
+
+            // set these in model
+            bool writeInTimeDir_;
+            bool writeInCase_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("stateController");
+
+    //- Declare runtime constructor selection table
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            stateController,
+            dictionary,
+            (
+                polyatomicCloud& cloud,
+                const dictionary& dict
+            ),
+            (t, cloud, dict)
+        );
+
+    // Constructors
+
+        //- Construct from components
+        stateController
+        (
+            polyatomicCloud& cloud,
+            const dictionary& dict
+        );
+
+
+    // Selectors
+
+        static autoPtr<stateController> New
+        (
+            polyatomicCloud& cloud,
+            const dictionary& dict
+        );
+
+
+    // Destructor
+
+        virtual ~stateController();
+
+
+    // Member Functions
+
+
+        void updateTime();
+
+        //- create an initial configuration
+        virtual void initialConfiguration() = 0;
+
+        //- calculate any required properties
+        virtual void calculateProperties() = 0;
+
+        //- control molecules at different stages of the integration time-step
+        virtual void controlMolsBeg() = 0;
+
+        virtual void controlBeforeForces() = 0;
+
+        virtual void controlMols() = 0;
+
+        virtual void controlMolsEnd() = 0;
+
+
+        //- output data
+        virtual void output
+        (
+            const fileName& fixedPathName,
+            const fileName& timePath
+        ) = 0;
+
+        //- update properties from a modified dictionary
+        virtual void updateProperties(const dictionary&) = 0;
+
+        void updateStateControllerProperties(const dictionary&);
+
+        // Access
+
+            //- return the control zone cells
+            const labelList& controlZone() const;
+
+            //- return the control zone name
+            const word& regionName() const;
+
+            //- return the targeted fields
+            scalar density() const;
+            scalar& density();
+
+            const vector& velocity() const;
+            vector& velocity();
+
+            scalar temperature() const;
+            scalar& temperature();
+
+            scalar pressure() const;
+            scalar& pressure();
+
+            const tensor& strainRate() const;
+            tensor& strainRate();
+
+            const vector& tempGradient() const;
+            vector& tempGradient();
+
+            //- return the targeted fields
+            const scalarField& densityField() const;
+            scalarField& densityField();
+
+            const vectorField& velocityField() const;
+            vectorField& velocityField();
+
+            const scalarField& temperatureField() const;
+            scalarField& temperatureField();
+
+            const scalarField& pressureField() const;
+            scalarField& pressureField();
+
+            bool singleValueController() const;
+            bool& singleValueController();
+
+            bool fieldController() const;
+            bool& fieldController();
+
+            bool writeInTimeDir() const;
+            bool writeInCase() const;
+
+            scalar avReqDensity();
+            vector avReqVelocity();
+            scalar avReqTemperature();
+            scalar avReqPressure();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.C b/src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.C
similarity index 66%
rename from src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.C
rename to src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.C
index 431f37aeb4..08ea263dd6 100644
--- a/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.C
+++ b/src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.C
@@ -21,12 +21,7 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Class
-    constPropSite
-
-Description
-
-\*----------------------------------------------------------------------------*/
+\*---------------------------------------------------------------------------*/
 
 #include "constPropSite.H"
 
@@ -71,6 +66,49 @@ Foam::constPropSite::~constPropSite()
 {}
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Istream& Foam::operator>>(Istream& is, constPropSite& cPS)
+{
+    is  >> cPS.siteReferencePosition_
+        >> cPS.siteMass_
+        >> cPS.siteCharge_
+        >> cPS.siteId_
+        >> cPS.name_
+        >> cPS.pairPotentialSite_
+        >> cPS.electrostaticSite_;
+
+    // Check state of Istream
+    is.check
+    (
+        "Foam::Istream& Foam::operator>>"
+        "(Foam::Istream&, Foam::constPropSite&)"
+    );
+
+    return is;
+}
+
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const constPropSite& cPS)
+{
+
+    os  << token::SPACE << cPS.siteReferencePosition()
+        << token::SPACE << cPS.siteMass()
+        << token::SPACE << cPS.siteCharge()
+        << token::SPACE << cPS.siteId()
+        << token::SPACE << cPS.name()
+        << token::SPACE << cPS.pairPotentialSite()
+        << token::SPACE << cPS.electrostaticSite();
+
+    // Check state of Ostream
+    os.check
+    (
+        "Foam::Ostream& Foam::operator<<(Foam::Ostream&, "
+        "const Foam::constPropSite&)"
+    );
+
+    return os;
+}
+
 
 // ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.H b/src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.H
similarity index 86%
rename from src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.H
rename to src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.H
index 0ffb63b0a8..910c0ccaef 100644
--- a/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.H
@@ -20,11 +20,13 @@ License
 
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
 Class
 
 
 Description
 
+
 SourceFiles
     constPropSiteI.H
     constPropSite.C
@@ -35,14 +37,23 @@ SourceFiles
 #define constPropSite_H
 
 #include "vector.H"
-#include "IFstream.H"
-#include "OFstream.H"
+#include "IOstreams.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+// Forward declaration of classes
+class Istream;
+class Ostream;
+
+// Forward declaration of friend functions and operators
+class constPropSite;
+Istream& operator>>(Istream&, constPropSite&);
+Ostream& operator<<(Ostream&, const constPropSite&);
+
+
 /*---------------------------------------------------------------------------*\
                          Class constPropSite Declaration
 \*---------------------------------------------------------------------------*/
@@ -144,6 +155,18 @@ public:
 
             //-
             inline bool& electrostaticSite();
+
+
+    // Member Operators
+
+        inline bool operator==(const constPropSite&) const;
+        inline bool operator!=(const constPropSite&) const;
+
+
+    // IOstream Operators
+
+        friend Istream& operator>>(Istream&, constPropSite&);
+        friend Ostream& operator<<(Ostream&, const constPropSite&);
 };
 
 
diff --git a/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSiteI.H b/src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSiteI.H
similarity index 81%
rename from src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSiteI.H
rename to src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSiteI.H
index dbeb395920..9c648c94c8 100644
--- a/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSiteI.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSiteI.H
@@ -21,11 +21,6 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Class
-
-
-Description
-
 \*---------------------------------------------------------------------------*/
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -114,4 +109,31 @@ inline bool& Foam::constPropSite::electrostaticSite()
 }
 
 
+// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
+
+bool Foam::constPropSite::operator==
+(
+    const constPropSite& rhs
+) const
+{
+    return
+       siteReferencePosition_ == rhs.siteReferencePosition_
+    && siteMass_ == rhs.siteMass_
+    && siteCharge_ == rhs.siteCharge_
+    && siteId_ == rhs.siteId_
+    && name_ == rhs.name_
+    && pairPotentialSite_ == rhs.pairPotentialSite_
+    && electrostaticSite_ == rhs.electrostaticSite_;
+}
+
+
+bool Foam::constPropSite::operator!=
+(
+    const constPropSite& rhs
+) const
+{
+    return !(*this == rhs);
+}
+
+
 // ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomic.C b/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomic.C
new file mode 100644
index 0000000000..f5b0316487
--- /dev/null
+++ b/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;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomic.H b/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomic.H
new file mode 100644
index 0000000000..d1e999435b
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicI.H b/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicI.H
new file mode 100644
index 0000000000..4485d53b8d
--- /dev/null
+++ b/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_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicIO.C b/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicIO.C
new file mode 100644
index 0000000000..a9a27bb992
--- /dev/null
+++ b/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;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomic.C b/src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.C
similarity index 93%
rename from src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomic.C
rename to src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.C
index dc22e42cf1..18fbfbb5c9 100644
--- a/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomic.C
+++ b/src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.C
@@ -23,11 +23,19 @@ License
 
 \*----------------------------------------------------------------------------*/
 
-#include "polyatomicCloud.H"
 #include "polyatomic.H"
 #include "Random.H"
 #include "Time.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(polyatomic, 0);
+    defineTemplateTypeNameAndDebug(Cloud<polyatomic>, 0);
+}
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 Foam::tensor Foam::polyatomic::rotationTensorX(scalar phi) const
@@ -81,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
@@ -90,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
 
@@ -108,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
@@ -149,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
@@ -285,7 +293,7 @@ void Foam::polyatomic::hitWallPatch
 )
 {
     // Use of the normal from tetIs is not required as
-    // hasWallImpactDistance for a polyatomicCloud is false.
+    // hasWallImpactDistance is false.
     vector nw = normal();
     nw /= mag(nw);
 
diff --git a/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomic.H b/src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.H
similarity index 82%
rename from src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomic.H
rename to src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.H
index 3b314ff2ac..0a7797516b 100644
--- a/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomic.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.H
@@ -42,15 +42,13 @@ SourceFiles
 #include "autoPtr.H"
 #include "diagTensor.H"
 #include "constPropSite.H"
+#include "MoleculeCloud.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-// Class forward declarations
-class polyatomicCloud;
-
 /*---------------------------------------------------------------------------*\
                            Class polyatomic Declaration
 \*---------------------------------------------------------------------------*/
@@ -91,22 +89,19 @@ public:
             //- Moment of intertia (in principal axis configiration)
             diagTensor momentOfInertia_;
 
+            //-
             scalar mass_;
 
 
         // Private Member Functions
 
-            void setInteracionSiteBools
-            (
-                const List<word>& siteIds,
-                const List<word>& pairPotSiteIds
-            );
-
+            //-
             bool linearMoleculeTest() const;
 
 
     public:
 
+        //-
         inline constantProperties();
 
         //- Construct from dictionary
@@ -118,22 +113,31 @@ public:
 
         // Member functions
 
+        //-
         inline const List<constPropSite>& sites() const;
 
+        //-
         inline const List<label>& pairPotSites() const;
 
+        //-
         inline const List<label>& electrostaticSites() const;
 
+        //-
         inline const diagTensor& momentOfInertia() const;
 
+        //-
         inline bool linearMolecule() const;
 
+        //-
         inline bool pointMolecule() const;
 
+        //-
         inline label degreesOfFreedom() const;
 
+        //-
         inline scalar mass() const;
 
+        //-
         inline label nSites() const;
     };
 
@@ -141,8 +145,22 @@ public:
     //- Class used to pass tracking data to the trackToFace function
     class trackingData
     :
-        public particle::TrackingData<polyatomicCloud>
+        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_;
 
@@ -151,9 +169,9 @@ public:
 
         // Constructors
 
-            trackingData(polyatomicCloud& cloud, label part)
+            trackingData(MoleculeCloud<polyatomic>& cloud, trackPart part)
             :
-                particle::TrackingData<polyatomicCloud>(cloud),
+                particle::TrackingData<MoleculeCloud<polyatomic> >(cloud),
                 part_(part)
             {}
 
@@ -176,10 +194,10 @@ private:
         //  bodyLocalVector = Q_.T() & globalVector
         tensor Q_;
 
-        // Linear velocity of polyatomic
+        //- Linear velocity of polyatomic
         vector v_;
 
-        // Total linear acceleration of polyatomic
+        //- Total linear acceleration of polyatomic
         vector a_;
 
         //- Angular momentum of polyatomic, in body local reference frame
@@ -188,8 +206,10 @@ private:
         //- Total torque on polyatomic, in body local reference frame
         vector tau_;
 
+        //-
         vector specialPosition_;
 
+        //-
         scalar potentialEnergy_;
 
         // - r_ij f_ij, stress dyad
@@ -198,31 +218,59 @@ private:
 //         // - r_ij outer product f_ij: virial contribution
 //         tensor rDotf_;
 
+        //-
         label special_;
 
+        //-
         label id_;
 
+        //-
         List<vector> siteForces_;
 
+        //-
         List<vector> sitePositions_;
 
 
     // Private Member Functions
 
+
+        //-
         tensor rotationTensorX(scalar deltaT) const;
 
+        //-
         tensor rotationTensorY(scalar deltaT) const;
 
+        //-
         tensor rotationTensorZ(scalar deltaT) const;
 
 
 public:
 
+    //- Runtime type information
+    TypeName("polyatomic");
+
     friend class Cloud<polyatomic>;
 
     // Constructors
 
-        //- Construct from components
+        //- 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
         (
             const polyMesh& mesh,
@@ -279,53 +327,91 @@ public:
 
         // 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 tensor& Q() const;
+
+            //-
             inline tensor& Q();
 
+            //-
             inline const vector& v() const;
+
+            //-
             inline vector& v();
 
+            //-
             inline const vector& a() const;
+
+            //-
             inline vector& a();
 
+            //-
             inline const vector& pi() const;
+
+            //-
             inline vector& pi();
 
+            //-
             inline const vector& tau() const;
+
+            //-
             inline vector& tau();
 
+            //-
             inline const List<vector>& siteForces() const;
+
+            //-
             inline List<vector>& siteForces();
 
+            //-
             inline const List<vector>& sitePositions() const;
+
+            //-
             inline List<vector>& 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;
 
 
@@ -367,10 +453,15 @@ public:
 
     // I-O
 
+        //- Read
         static void readFields(Cloud<polyatomic>& mC);
 
+        //- Write
         static void writeFields(const Cloud<polyatomic>& mC);
 
+        //- Show info
+        static void info(trackingData& td);
+
 
     // IOstream Operators
 
diff --git a/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomicI.H b/src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomicI.H
similarity index 85%
rename from src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomicI.H
rename to src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomicI.H
index faf24d1f1c..57149e5fa9 100644
--- a/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomicI.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomicI.H
@@ -57,11 +57,12 @@ inline Foam::polyatomic::constantProperties::constantProperties
 
     List<word> pairPotentialIds(dict.lookup("pairPotentialSiteIds"));
 
+    sites_.setSize(siteReferencePositions.size());
+
     if
     (
-        ( siteIdNames.size() != sites_.size() )
-     || ( siteReferencePositions.size() != sites_.size() )
-     || ( siteCharges.size() != sites_.size() )
+        (siteIdNames.size() != sites_.size())
+     || (siteCharges.size() != sites_.size())
     )
     {
             FatalErrorIn
@@ -72,12 +73,13 @@ inline Foam::polyatomic::constantProperties::constantProperties
                 ")"
             )
             << "Sizes of site id, charge and "
-            << "referencePositions are not the same: " << sites_.size()
-            << nl << abort(FatalError);
+            << "referencePositions are not the same: " << nl
+            << siteIdNames  << nl
+            << siteReferencePositions << nl
+            << siteCharges << nl
+            << abort(FatalError);
     }
 
-    sites_.setSize(siteIdNames.size());
-
     electrostaticSites_.setSize(sites_.size(), -1);
     pairPotSites_.setSize(sites_.size(), -1);
 
@@ -107,7 +109,7 @@ inline Foam::polyatomic::constantProperties::constantProperties
             electrostaticSites_[electrostaticI++] = sI;
         }
 
-        if (sites_[sI].pairPotentialSite() && !sites_[sI].electrostaticSite())
+        if (!sites_[sI].pairPotentialSite() && !sites_[sI].electrostaticSite())
         {
             WarningIn
             (
@@ -303,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,
@@ -315,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),
@@ -332,10 +393,10 @@ inline Foam::polyatomic::polyatomic
     rf_(tensor::zero),
     special_(special),
     id_(id),
-    siteForces_(constProps.nSites(), vector::zero),
-    sitePositions_(constProps.nSites())
+    siteForces_(cP.nSites(), vector::zero),
+    sitePositions_(cP.nSites())
 {
-    setSitePositions(constProps);
+    setSitePositions(cP);
 }
 
 
@@ -413,7 +474,7 @@ inline bool Foam::polyatomic::constantProperties::linearMolecule() const
 }
 
 
-inline bool Foam::polyatomic::constantProperties::pointMolecule() const
+inline  bool Foam::polyatomic::constantProperties::pointMolecule() const
 {
     return (momentOfInertia_.zz() < 0);
 }
diff --git a/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomicIO.C b/src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomicIO.C
similarity index 67%
rename from src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomicIO.C
rename to src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomicIO.C
index 0d93a071c0..86dcbfe49f 100644
--- a/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomicIO.C
+++ b/src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomicIO.C
@@ -25,7 +25,7 @@ License
 
 #include "polyatomic.H"
 #include "IOstreams.H"
-#include "polyatomicCloud.H"
+#include "Cloud.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -92,7 +92,11 @@ Foam::polyatomic::polyatomic
     is.check
     (
         "Foam::polyatomic::polyatomic"
-        "(const Cloud<polyatomic>& cloud, Foam::Istream&), bool"
+        "("
+            "const polyMesh& mesh,"
+            "Istream& is,"
+            "bool readFields"
+        ")"
     );
 }
 
@@ -134,7 +138,8 @@ void Foam::polyatomic::readFields(Cloud<polyatomic>& mC)
     mC.checkFieldIOobject(mC, id);
 
     label i = 0;
-    forAllIter(polyatomicCloud, mC, iter)
+
+    forAllIter(typename Cloud<polyatomic>, mC, iter)
     {
         polyatomic& mol = iter();
 
@@ -203,7 +208,7 @@ void Foam::polyatomic::writeFields(const Cloud<polyatomic>& mC)
     );
 
     label i = 0;
-    forAllConstIter(polyatomicCloud, mC, iter)
+    forAllConstIter(typename Cloud<polyatomic>, mC, iter)
     {
         const polyatomic& mol = iter();
 
@@ -241,17 +246,114 @@ void Foam::polyatomic::writeFields(const Cloud<polyatomic>& mC)
     orientation1.write();
     orientation2.write();
     orientation3.write();
+}
 
-    Info<< "writeFields " << mC.name() << endl;
 
-    if (isA<polyatomicCloud>(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 polyatomicCloud& m = dynamic_cast<const polyatomicCloud&>(mC);
+        const label molId = mol().id();
+        scalar molMass(td.cloud().constProps(molId).mass());
+        totalMass += molMass;
+        //CentreOfMass += mol().position()*molMass;
+    }
 
-        m.writeXYZ
+    // if (nMols)
+    // {
+    //     CentreOfMass /= totalMass;
+    // }
+
+    forAllConstIter(typename Cloud<polyatomic>, td.cloud(), mol)
+    {
+        const label molId = mol().id();
+        const polyatomic::constantProperties cP
         (
-            m.mesh().time().timePath()/cloud::prefix/"polyatomicCloud.xmol"
+            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;
     }
 }
 
diff --git a/src/lagrangian/molecularDynamics/potential/potential/potential.C b/src/lagrangian/molecularDynamics/potential/potential/potential.C
index 4d1f0df703..76a624a652 100644
--- a/src/lagrangian/molecularDynamics/potential/potential/potential.C
+++ b/src/lagrangian/molecularDynamics/potential/potential/potential.C
@@ -93,7 +93,10 @@ void Foam::potential::setSiteIdList(const dictionary& moleculePropertiesDict)
 }
 
 
-void Foam::potential::potential::readPotentialDict()
+void Foam::potential::potential::readPotentialDict
+(
+    const dictionary& moleculePropertiesDict
+)
 {
     Info<< nl <<  "Reading potential dictionary:" << endl;
 
@@ -111,31 +114,28 @@ void Foam::potential::potential::readPotentialDict()
 
     idList_ = List<word>(idListDict.lookup("idList"));
 
-    setSiteIdList
-    (
-        IOdictionary
-        (
-            IOobject
-            (
-                "moleculeProperties",
-                mesh_.time().constant(),
-                mesh_,
-                IOobject::MUST_READ_IF_MODIFIED,
-                IOobject::NO_WRITE,
-                false
-            )
-        )
-    );
+    setSiteIdList(moleculePropertiesDict);
 
     List<word> pairPotentialSiteIdList
     (
         SubList<word>(siteIdList_, nPairPotIds_)
     );
 
-    Info<< nl << "Unique site ids found: " << siteIdList_
-        << nl << "Site Ids requiring a pair potential: "
-        << pairPotentialSiteIdList
-        << endl;
+    Info<< nl << "Unique site ids found:";
+
+    forAll(siteIdList_, i)
+    {
+        Info<< " " << siteIdList_[i];
+    }
+
+    Info << nl << "Site Ids requiring a pair potential:";
+
+    forAll(pairPotentialSiteIdList, i)
+    {
+        Info<< " " << pairPotentialSiteIdList[i];
+    }
+
+    Info<< nl;
 
     List<word> tetherSiteIdList(0);
 
@@ -229,37 +229,25 @@ void Foam::potential::potential::readPotentialDict()
 
     if (potentialDict.found("external"))
     {
-        Info<< nl << "Reading external forces:" << endl;
+        Info<< nl << "Reading external forces: ";
 
         const dictionary& externalDict = potentialDict.subDict("external");
 
         // gravity
         externalDict.readIfPresent("gravity", gravity_);
-    }
 
-    Info<< nl << tab << "gravity = " << gravity_ << endl;
+        Info<< "gravity = " << gravity_ << nl << endl;
+    }
 }
 
 
 void Foam::potential::potential::readMdInitialiseDict
 (
-    const IOdictionary& mdInitialiseDict,
-    IOdictionary& idListDict
+    const dictionary& mdInitialiseDict,
+    const dictionary& moleculePropertiesDict,
+    dictionary& idListDict
 )
 {
-    IOdictionary moleculePropertiesDict
-    (
-        IOobject
-        (
-            "moleculeProperties",
-            mesh_.time().constant(),
-            mesh_,
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE,
-            false
-        )
-    );
-
     DynamicList<word> idList;
 
     DynamicList<word> tetherSiteIdList;
@@ -286,11 +274,12 @@ void Foam::potential::potential::readMdInitialiseDict
                 (
                     "potential::readMdInitialiseDict"
                     "("
-                        "const IOdictionary&, "
-                        "IOdictionary&"
+                        "const dictionary&, "
+                        "dictionary&"
                     ")"
                 )   << "Molecule type " << id
-                    << " not found in moleculeProperties dictionary." << nl
+                    << " not found in " << moleculePropertiesDict.name()
+                    << " dictionary." << nl
                     << abort(FatalError);
             }
 
@@ -341,8 +330,8 @@ void Foam::potential::potential::readMdInitialiseDict
                 (
                     "potential::readMdInitialiseDict"
                     "("
-                        "const IOdictionary&, "
-                        "IOdictionary&"
+                        "const dictionary&, "
+                        "dictionary&"
                     ")"
                 )   << "Tether id  " << tetherSiteId
                     << " not found as a site of any molecule in zone." << nl
@@ -364,24 +353,34 @@ void Foam::potential::potential::readMdInitialiseDict
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::potential::potential(const polyMesh& mesh)
+Foam::potential::potential
+(
+    const polyMesh& mesh,
+    const dictionary& moleculePropertiesDict
+)
 :
     mesh_(mesh)
 {
-    readPotentialDict();
+    readPotentialDict(moleculePropertiesDict);
 }
 
 
 Foam::potential::potential
 (
     const polyMesh& mesh,
-    const IOdictionary& mdInitialiseDict,
-    IOdictionary& idListDict
+    const dictionary& mdInitialiseDict,
+    const dictionary& moleculePropertiesDict,
+    dictionary& idListDict
 )
 :
     mesh_(mesh)
 {
-    readMdInitialiseDict(mdInitialiseDict, idListDict);
+    readMdInitialiseDict
+    (
+        mdInitialiseDict,
+        moleculePropertiesDict,
+        idListDict
+    );
 }
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
diff --git a/src/lagrangian/molecularDynamics/potential/potential/potential.H b/src/lagrangian/molecularDynamics/potential/potential/potential.H
index 8f24458558..02e08c65d0 100644
--- a/src/lagrangian/molecularDynamics/potential/potential/potential.H
+++ b/src/lagrangian/molecularDynamics/potential/potential/potential.H
@@ -55,35 +55,49 @@ class potential
 {
     // Private data
 
+        //-
         const polyMesh& mesh_;
 
+        //-
         List<word> idList_;
 
+        //-
         List<word> siteIdList_;
 
+        //-
         label nPairPotIds_;
 
+        //-
         scalar potentialEnergyLimit_;
 
+        //-
         labelList removalOrder_;
 
+        //-
         pairPotentialList pairPotentials_;
 
+        //-
         tetherPotentialList tetherPotentials_;
 
+        //-
         vector gravity_;
 
 
     // Private Member Functions
 
+
+        //-
         void setSiteIdList(const dictionary& moleculePropertiesDict);
 
-        void readPotentialDict();
+        //-
+        void readPotentialDict(const dictionary& moleculePropertiesDict);
 
+        //-
         void readMdInitialiseDict
         (
-            const IOdictionary& mdInitialiseDict,
-            IOdictionary& idListDict
+            const dictionary& mdInitialiseDict,
+            const dictionary& moleculePropertiesDict,
+            dictionary& idListDict
         );
 
         //- Disallow default bitwise copy construct
@@ -98,14 +112,19 @@ public:
     // Constructors
 
         //- Construct from mesh reference
-        potential(const polyMesh& mesh);
+        potential
+        (
+            const polyMesh& mesh,
+            const dictionary& moleculePropertiesDict
+        );
 
         //- Construct from mdInitialiseDict
         potential
         (
             const polyMesh& mesh,
-            const IOdictionary& mdInitialiseDict,
-            IOdictionary& idListDict
+            const dictionary& mdInitialiseDict,
+            const dictionary& moleculePropertiesDict,
+            dictionary& idListDict
         );
 
 
@@ -117,22 +136,31 @@ public:
 
         // Access
 
+            //-
             inline label nIds() const;
 
+            //-
             inline const List<word>& idList() const;
 
+            //-
             inline const List<word>& siteIdList() const;
 
+            //-
             inline scalar potentialEnergyLimit() const;
 
+            //-
             inline label nPairPotentials() const;
 
+            //-
             inline const labelList& removalOrder() const;
 
+            //-
             inline const pairPotentialList& pairPotentials() const;
 
+            //-
             inline const tetherPotentialList& tetherPotentials() const;
 
+            //-
             inline const vector& gravity() const;
 };
 
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index fa208c73e4..73e0ed3d43 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -160,9 +160,19 @@ Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment
 
     if (closestSpokeSurface == -1)
     {
-        FatalErrorIn("conformalVoronoiMesh::requiredAlignment")
-            << "No secondary surface hit found in spoke search."
-            << nl << exit(FatalError);
+        WarningIn
+        (
+            "conformalVoronoiMesh::requiredAlignment"
+            "("
+                "const Foam::point& pt"
+            ") const"
+        )
+            << "No secondary surface hit found in spoke search "
+            << "using " << s
+            << " spokes, try increasing alignmentSearchSpokes."
+            << endl;
+
+        return I;
     }
 
     // Auxiliary alignment generated by spoke intersection normal.
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
index b7fc0b3a1c..8d491429e4 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
@@ -1110,7 +1110,7 @@ Foam::conformalVoronoiMesh::collapseFace
 
         // Output face and collapse axis for visualisation
 
-        Info<< "# Aspect ratio = " << aspectRatio << nl
+        Pout<< "# Aspect ratio = " << aspectRatio << nl
             << "# inertia = " << J << nl
             << "# determinant = " << detJ << nl
             << "# eigenvalues = " << eigenValues(J) << nl
@@ -1120,17 +1120,17 @@ Foam::conformalVoronoiMesh::collapseFace
 
         forAll(f, fPtI)
         {
-            meshTools::writeOBJ(Info, pts[f[fPtI]]);
+            meshTools::writeOBJ(Pout, pts[f[fPtI]]);
         }
 
-        Info<< "f";
+        Pout<< "f";
 
         forAll(f, fPtI)
         {
-            Info << " " << fPtI + 1;
+            Pout << " " << fPtI + 1;
         }
 
-        Info<< endl;
+        Pout<< endl;
 
         return fcmNone;
     }
@@ -1402,7 +1402,7 @@ Foam::conformalVoronoiMesh::collapseFace
     // {
     //     // Output face and collapse axis for visualisation
 
-    //     Info<< "# Aspect ratio = " << aspectRatio << nl
+    //     Pout<< "# Aspect ratio = " << aspectRatio << nl
     //         << "# determinant = " << detJ << nl
     //         << "# collapseAxis = " << collapseAxis << nl
     //         << "# mode = " << mode << nl
@@ -1412,33 +1412,33 @@ Foam::conformalVoronoiMesh::collapseFace
 
     //     scalar scale = 2.0*mag(fC - pts[f[0]]);
 
-    //     meshTools::writeOBJ(Info, fC);
-    //     meshTools::writeOBJ(Info, fC + scale*collapseAxis);
+    //     meshTools::writeOBJ(Pout, fC);
+    //     meshTools::writeOBJ(Pout, fC + scale*collapseAxis);
 
-    //     Info<< "f 1 2" << endl;
+    //     Pout<< "f 1 2" << endl;
 
     //     forAll(f, fPtI)
     //     {
-    //         meshTools::writeOBJ(Info, pts[f[fPtI]]);
+    //         meshTools::writeOBJ(Pout, pts[f[fPtI]]);
     //     }
 
-    //     Info<< "f";
+    //     Pout<< "f";
 
     //     forAll(f, fPtI)
     //     {
-    //         Info << " " << fPtI + 3;
+    //         Pout << " " << fPtI + 3;
     //     }
 
-    //     Info<< nl << "# " << d << endl;
+    //     Pout<< nl << "# " << d << endl;
 
-    //     Info<< "# " << d.first() << " " << d.last() << endl;
+    //     Pout<< "# " << d.first() << " " << d.last() << endl;
 
     //     forAll(d, dI)
     //     {
-    //         meshTools::writeOBJ(Info, fC + (d[dI] - dShift)*collapseAxis);
+    //         meshTools::writeOBJ(Pout, fC + (d[dI] - dShift)*collapseAxis);
     //     }
 
-    //     Info<< endl;
+    //     Pout<< endl;
     // }
 
     return mode;
@@ -1490,7 +1490,7 @@ void Foam::conformalVoronoiMesh::deferredCollapseFaceSet
         }
     }
 
-    Info<< "facesToCollapse" << nl << faceLabels << endl;
+    Pout<< "facesToCollapse" << nl << faceLabels << endl;
 }
 
 
@@ -1632,7 +1632,8 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
 
     if (checkFaces.size() < fAreas.size())
     {
-        Info<< "Excluding " << fAreas.size() - checkFaces.size()
+        Info<< "Excluding "
+            << returnReduce(fAreas.size() - checkFaces.size(), sumOp<label>())
             << " faces from check, < " << faceAreaLimit << " area" << endl;
     }
 
@@ -1655,7 +1656,7 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
         {
             if (cells[cI].size() < 4 && cells[cI].size() > 0)
             {
-                // Info<< "cell " << cI << " " << cells[cI]
+                // Pout<< "cell " << cI << " " << cells[cI]
                 //     << " has " << cells[cI].size() << " faces."
                 //     << endl;
 
@@ -1669,7 +1670,8 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
         }
 
         Info<< "    cells with more than 1 but fewer than 4 faces          : "
-            << nInvalidPolyhedra << endl;
+            << returnReduce(nInvalidPolyhedra, sumOp<label>())
+            << endl;
 
         // Check for cells with one internal face only
 
@@ -1712,7 +1714,8 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
         }
 
         Info<< "    cells with with zero or one non-boundary face          : "
-            << oneInternalFaceCells << endl;
+            << returnReduce(oneInternalFaceCells, sumOp<label>())
+            << endl;
     }
 
 
@@ -1795,7 +1798,7 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
     }
 
     Info<< nl << "Maximum number of filter limits applied: "
-        << maxFilterCount << endl;
+        << returnReduce(maxFilterCount, maxOp<label>()) << endl;
 
     return wrongFaces;
 }
@@ -2489,7 +2492,9 @@ void Foam::conformalVoronoiMesh::removeUnusedPoints
 
     inplaceReorder(oldToNew, pts);
 
-    Info<< "    Removing " << pts.size() - pointI << " unused points"
+    Info<< "    Removing "
+        << returnReduce(pts.size() - pointI, sumOp<label>())
+        << " unused points"
         << endl;
 
     pts.setSize(pointI);
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
index 9ab849e7d1..de3a0ee99f 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
@@ -696,9 +696,10 @@ void Foam::conformalVoronoiMesh::findRemainingProtrusionSet
         }
     }
 
-    if (!protrudingCells.empty())
+    if (returnReduce(protrudingCells.size(), sumOp<label>()) > 0)
     {
-        Pout<< nl << "Found " << protrudingCells.size()
+        Info<< nl << "Found "
+            << returnReduce(protrudingCells.size(), sumOp<label>())
             << " cells protruding from the surface, writing cellSet "
             << protrudingCells.name()
             << endl;
diff --git a/src/postProcessing/functionObjects/IO/partialWrite/partialWrite.C b/src/postProcessing/functionObjects/IO/partialWrite/partialWrite.C
index 327c7b804a..1ea7ae253d 100644
--- a/src/postProcessing/functionObjects/IO/partialWrite/partialWrite.C
+++ b/src/postProcessing/functionObjects/IO/partialWrite/partialWrite.C
@@ -27,6 +27,7 @@ License
 #include "dictionary.H"
 #include "Time.H"
 #include "IOobjectList.H"
+#include "polyMesh.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -123,14 +124,26 @@ void Foam::partialWrite::write()
         else
         {
             // Delete all but marked objects
+            fileName dbDir;
+            if (isA<polyMesh>(obr_))
+            {
+                dbDir = dynamic_cast<const polyMesh&>(obr_).dbDir();
+            }
+
             IOobjectList objects(obr_, obr_.time().timeName());
 
             forAllConstIter(HashPtrTable<IOobject>, objects, iter)
             {
                 if (!objectNames_.found(iter()->name()))
                 {
-                    const fileName f = obr_.time().timePath()/iter()->name();
-                    //Pout<< "   rm " << f << endl;
+                    const fileName f =
+                        obr_.time().timePath()
+                       /dbDir
+                       /iter()->name();
+                    if (debug)
+                    {
+                        Pout<< "   rm " << f << endl;
+                    }
                     rm(f);
                 }
             }
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/Allrun b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/Allrun
index ae53653faa..026f614b53 100755
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/Allrun
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/Allrun
@@ -32,4 +32,8 @@ decomposePar -region panelRegion > log.decomposParPanelRegion.log 2>&1
 
 runParallel `getApplication` 6
 
+paraFoam -touch
+paraFoam -touch -region panelRegion
+
+
 # ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties
index 7422e40c44..ccd5d89ef8 100644
--- a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties
+++ b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties
@@ -112,6 +112,7 @@ subModels
                 alpha           0.12;
                 b               1.5;
                 mu              0.43;
+                cohesionEnergyDensity 0;
             }
             frontAndBack
             {
@@ -120,6 +121,7 @@ subModels
                 alpha           0.12;
                 b               1.5;
                 mu              0.1;
+                cohesionEnergyDensity 0;
             }
         };
     }
diff --git a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties
index 36448f7b30..3adb6ab771 100644
--- a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties
+++ b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties
@@ -121,6 +121,7 @@ subModels
                 alpha           0.12;
                 b               1.5;
                 mu              0.43;
+                cohesionEnergyDensity 0;
             }
             frontAndBack
             {
@@ -129,6 +130,7 @@ subModels
                 alpha           0.12;
                 b               1.5;
                 mu              0.1;
+                cohesionEnergyDensity 0;
             }
         };
     }
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/Allrun.pre b/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/Allrun.pre
index 747433b184..bd6fc4f90f 100755
--- a/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/Allrun.pre
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/Allrun.pre
@@ -8,3 +8,6 @@ runApplication setSet -batch wallFilmRegion.setSet
 mv log.setSet log.wallFilmRegion.setSet
 
 runApplication extrudeToRegionMesh -overwrite
+
+paraFoam -touch
+paraFoam -touch -region wallFilmRegion
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/Allrun.pre b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/Allrun.pre
index 5583402201..b5423180b6 100755
--- a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/Allrun.pre
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/Allrun.pre
@@ -32,3 +32,6 @@ cp -r system/wallFilmRegion.org system/wallFilmRegion
 
 find ./0 -maxdepth 1 -type f -exec \
     sed -i "s/wallFilm/\"(region0_to.*)\"/g" {} \;
+
+paraFoam -touch
+paraFoam -touch -region wallFilmRegion
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/Allrun.pre b/tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/Allrun.pre
index 0e3c1209c5..d9b17ad097 100755
--- a/tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/Allrun.pre
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/Allrun.pre
@@ -10,3 +10,5 @@ mv log.setSet log.wallFilmRegion.setSet
 
 runApplication extrudeToRegionMesh -overwrite
 
+paraFoam -touch
+paraFoam -touch -region wallFilmRegion
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/Allrun.pre b/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/Allrun.pre
index 14da52453e..bbae3db32c 100755
--- a/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/Allrun.pre
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/Allrun.pre
@@ -14,3 +14,6 @@ runApplication setSet -region wallFilmRegion -batch createWallFilmRegionPatches.
 mv log.setSet log.createWallFilmRegionPatches.setSet
 
 runApplication createPatch -region wallFilmRegion -overwrite
+
+paraFoam -touch
+paraFoam -touch -region wallFilmRegion