diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/files b/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/files
new file mode 100644
index 0000000000..c3a050bb68
--- /dev/null
+++ b/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/files
@@ -0,0 +1,3 @@
+mdEquilibrationFoam.C
+
+EXE = $(FOAM_APPBIN)/mdEquilibrationFoam
diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/options b/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/options
new file mode 100644
index 0000000000..91cab4609b
--- /dev/null
+++ b/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/options
@@ -0,0 +1,15 @@
+EXE_INC = \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools \
+    -lfiniteVolume \
+    -llagrangian \
+    -lmolecule \
+    -lpotential \
+    -lmolecularMeasurements
diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C b/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C
new file mode 100644
index 0000000000..96c07a1756
--- /dev/null
+++ b/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C
@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+ =========                   |
+ \\      /   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/>.
+
+Application
+    mdEquilibrationFoam
+
+Description
+    Equilibrates and/or preconditions molecular dynamics systems
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "md.H"
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nReading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    potential pot(mesh);
+
+    moleculeCloud molecules(mesh, pot);
+
+    #include "temperatureAndPressureVariables.H"
+
+    #include "readmdEquilibrationDict.H"
+
+    label nAveragingSteps = 0;
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        nAveragingSteps++;
+
+        Info<< "Time = " << runTime.timeName() << endl;
+
+        molecules.evolve();
+
+        #include "meanMomentumEnergyAndNMols.H"
+
+        #include "temperatureAndPressure.H"
+
+        #include "temperatureEquilibration.H"
+
+        runTime.write();
+
+        if (runTime.outputTime())
+        {
+            nAveragingSteps = 0;
+        }
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H b/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H
new file mode 100644
index 0000000000..3f43981907
--- /dev/null
+++ b/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H
@@ -0,0 +1,18 @@
+Info<< nl << "Reading MD Equilibration Dictionary" << nl << endl;
+
+IOdictionary mdEquilibrationDict
+(
+    IOobject
+    (
+        "mdEquilibrationDict",
+        runTime.system(),
+        mesh,
+        IOobject::MUST_READ_IF_MODIFIED,
+        IOobject::NO_WRITE
+    )
+);
+
+scalar targetTemperature = readScalar
+(
+    mdEquilibrationDict.lookup("targetTemperature")
+);
diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdFoam/Make/options b/applications/solvers/discreteMethods/molecularDynamics/mdFoam/Make/options
index 98f09f1b4b..91cab4609b 100644
--- a/applications/solvers/discreteMethods/molecularDynamics/mdFoam/Make/options
+++ b/applications/solvers/discreteMethods/molecularDynamics/mdFoam/Make/options
@@ -1,9 +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)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
@@ -13,4 +11,5 @@ EXE_LIBS = \
     -lfiniteVolume \
     -llagrangian \
     -lmolecule \
-    -lpotential
+    -lpotential \
+    -lmolecularMeasurements
diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C b/applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C
index 38510dc453..779172803d 100644
--- a/applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C
+++ b/applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C
@@ -30,31 +30,60 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "monoatomicCloud.H"
-#include "polyatomicCloud.H"
+#include "md.H"
 
 int main(int argc, char *argv[])
 {
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "createFields.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+    Info<< "\nReading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    potential pot(mesh);
+
+    moleculeCloud molecules(mesh, pot);
+
+    #include "temperatureAndPressureVariables.H"
+
+    label nAveragingSteps = 0;
+
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.loop())
     {
+        nAveragingSteps++;
+
         Info<< "Time = " << runTime.timeName() << endl;
 
-        monoatomics.evolve();
+        molecules.evolve();
 
-        polyatomics.evolve();
+        #include "meanMomentumEnergyAndNMols.H"
+
+        #include "temperatureAndPressure.H"
 
         runTime.write();
 
-        Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+        if (runTime.outputTime())
+        {
+            nAveragingSteps = 0;
+        }
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
     }
diff --git a/applications/utilities/preProcessing/mdInitialise/Make/options b/applications/utilities/preProcessing/mdInitialise/Make/options
index 7633b07554..289c0990a2 100644
--- a/applications/utilities/preProcessing/mdInitialise/Make/options
+++ b/applications/utilities/preProcessing/mdInitialise/Make/options
@@ -3,6 +3,7 @@ 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
 
@@ -12,5 +13,6 @@ EXE_LIBS = \
     -lfiniteVolume \
     -llagrangian \
     -lmolecule \
-    -lpotential
+    -lpotential \
+    -lmolecularMeasurements
 
diff --git a/applications/utilities/preProcessing/mdInitialise/mdInitialise.C b/applications/utilities/preProcessing/mdInitialise/mdInitialise.C
index d29ef0a72d..3551e46133 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-2011 OpenCFD Ltd.
+   \\  /     A nd            | Copyright (C) 2008-2010 OpenCFD Ltd.
     \\/      M anipulation   |
 -------------------------------------------------------------------------------
 License
@@ -26,9 +26,8 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
+#include "md.H"
 #include "fvCFD.H"
-#include "polyatomicCloud.H"
-#include "monoatomicCloud.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,15 +50,11 @@ int main(int argc, char *argv[])
         )
     );
 
-    word polyCloudName("polyatomicCloud");
-
-    const dictionary& polyDict(mdInitialiseDict.subDict(polyCloudName));
-
-    IOdictionary polyIdListDict
+    IOdictionary idListDict
     (
         IOobject
         (
-            polyCloudName + "_idList",
+            "idList",
             mesh.time().constant(),
             mesh,
             IOobject::NO_READ,
@@ -67,88 +62,26 @@ int main(int argc, char *argv[])
         )
     );
 
-    potential polyPot
-    (
-        mesh,
-        polyDict,
-        IOdictionary
-        (
-            IOobject
-            (
-                polyCloudName +  "Properties",
-                mesh.time().constant(),
-                mesh,
-                IOobject::MUST_READ_IF_MODIFIED,
-                IOobject::NO_WRITE,
-                false
-            )
-        ),
-        polyIdListDict
-    );
+    potential pot(mesh, mdInitialiseDict, idListDict);
 
-    polyatomicCloud poly
-    (
-        polyCloudName,
-        mesh,
-        polyPot,
-        polyDict
-    );
+    moleculeCloud molecules(mesh, pot, mdInitialiseDict);
 
-    Info<< nl << returnReduce(poly.size(), sumOp<label>()) << " added to "
-        << poly.name()
+    label totalMolecules = molecules.size();
+
+    if (Pstream::parRun())
+    {
+        reduce(totalMolecules, sumOp<label>());
+    }
+
+    Info<< nl << "Total number of molecules added: " << totalMolecules
         << nl << endl;
 
-    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;
+    IOstream::defaultPrecision(15);
 
     if (!mesh.write())
     {
         FatalErrorIn(args.executable())
-            << "Failed writing."
+            << "Failed writing moleculeCloud."
             << nl << exit(FatalError);
     }
 
diff --git a/src/lagrangian/molecularDynamics/Allwmake b/src/lagrangian/molecularDynamics/Allwmake
index 118ba05f1a..7c17819ad5 100755
--- a/src/lagrangian/molecularDynamics/Allwmake
+++ b/src/lagrangian/molecularDynamics/Allwmake
@@ -4,6 +4,7 @@ makeType=${1:-libso}
 set -x
 
 wmake $makeType potential
+wmake $makeType molecularMeasurements
 wmake $makeType molecule
 
 # ----------------------------------------------------------------- end-of-file
diff --git a/src/lagrangian/molecularDynamics/molecularMeasurements/Make/files b/src/lagrangian/molecularDynamics/molecularMeasurements/Make/files
new file mode 100644
index 0000000000..8450dc8360
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecularMeasurements/Make/files
@@ -0,0 +1,3 @@
+distribution/distribution.C
+
+LIB = $(FOAM_LIBBIN)/libmolecularMeasurements
diff --git a/src/lagrangian/molecularDynamics/molecularMeasurements/Make/options b/src/lagrangian/molecularDynamics/molecularMeasurements/Make/options
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/src/lagrangian/molecularDynamics/molecularMeasurements/bufferedAccumulator/bufferedAccumulator.C b/src/lagrangian/molecularDynamics/molecularMeasurements/bufferedAccumulator/bufferedAccumulator.C
new file mode 100644
index 0000000000..1d9c964d4f
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecularMeasurements/bufferedAccumulator/bufferedAccumulator.C
@@ -0,0 +1,236 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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 "bufferedAccumulator.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<class Type>
+const char* const
+    Foam::bufferedAccumulator<Type>::typeName("bufferedAccumulator");
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+void Foam::bufferedAccumulator<Type>::accumulateAndResetBuffer(const label b)
+{
+    accumulationBuffer() += (*this)[b];
+
+    averagesTaken_++;
+
+    (*this)[b] = Field<Type>(bufferLength(), pTraits<Type>::zero);
+
+    bufferOffsets_[b] = 0;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::bufferedAccumulator<Type>::bufferedAccumulator()
+:
+    List< Field<Type> >(),
+    averagesTaken_(),
+    bufferOffsets_()
+{}
+
+
+template<class Type>
+Foam::bufferedAccumulator<Type>::bufferedAccumulator
+(
+    const label nBuffers,
+    const label bufferLength,
+    const label bufferingInterval
+)
+:
+    List< Field<Type> >(),
+    averagesTaken_(),
+    bufferOffsets_()
+{
+    setSizes
+    (
+        nBuffers,
+        bufferLength,
+        bufferingInterval
+    );
+}
+
+
+template<class Type>
+Foam::bufferedAccumulator<Type>::bufferedAccumulator
+(
+    const bufferedAccumulator<Type>& bA
+)
+:
+    List< Field<Type> >(static_cast< List< Field<Type> > >(bA)),
+    averagesTaken_(bA.averagesTaken()),
+    bufferOffsets_(bA.bufferOffsets())
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::bufferedAccumulator<Type>::~bufferedAccumulator()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::bufferedAccumulator<Type>::setSizes
+(
+    const label nBuffers,
+    const label bufferLength,
+    const label bufferingInterval
+)
+{
+    (*this).setSize(nBuffers + 1);
+
+    forAll((*this), b)
+    {
+        (*this)[b] = Field<Type>(bufferLength, pTraits<Type>::zero);
+    }
+
+    averagesTaken_ = 0;
+
+    bufferOffsets_.setSize(nBuffers);
+
+    forAll(bufferOffsets_, bO)
+    {
+        bufferOffsets_[bO] = -bufferingInterval * bO - 1;
+    }
+}
+
+
+template<class Type>
+Foam::label Foam::bufferedAccumulator<Type>::addToBuffers
+(
+    const List<Type>& valuesToAdd
+)
+{
+    label bufferToRefill = -1;
+
+    for (label b = 0; b < nBuffers(); b++)
+    {
+        Field<Type>& buf((*this)[b]);
+
+        label& bO = bufferOffsets_[b];
+
+        if (bO >= 0)
+        {
+            buf[bO] = valuesToAdd[b];
+        }
+
+        bO++;
+
+        if (bO == bufferLength())
+        {
+            accumulateAndResetBuffer(b);
+        }
+
+        if (bO == 0)
+        {
+            if (bufferToRefill != -1)
+            {
+                FatalErrorIn("bufferedAccumulator<Type>::addToBuffers ")
+                    << "More than one bufferedAccumulator accumulation "
+                    << "buffer filled at once, this is considered an error."
+                    << abort(FatalError);
+            }
+
+            bufferToRefill = b;
+        }
+    }
+
+    return bufferToRefill;
+}
+
+
+template<class Type>
+Foam::Field<Type> Foam::bufferedAccumulator<Type>::averaged() const
+{
+    if (averagesTaken_)
+    {
+        Field<Type> bA = accumulationBuffer()/averagesTaken_;
+
+        return bA;
+    }
+    else
+    {
+        WarningIn
+        (
+            "bufferedAccumulator<Type>::averagedbufferedAccumulator() const"
+        )   << "Averaged correlation function requested but averagesTaken = "
+            << averagesTaken_
+            << ". Returning empty field."
+            << endl;
+
+        return Field<Type>(bufferLength(), pTraits<Type>::zero);
+    }
+}
+
+
+template<class Type>
+void Foam::bufferedAccumulator<Type>::resetAveraging()
+{
+    accumulationBuffer() = Field<Type>(bufferLength(), pTraits<Type>::zero);
+
+    averagesTaken_ = 0;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::bufferedAccumulator<Type>::operator=
+(
+    const bufferedAccumulator<Type>& rhs
+)
+{
+    // Check for assignment to self
+    if (this == &rhs)
+    {
+        FatalErrorIn
+        (
+            "bufferedAccumulator<Type>::operator=(const bufferedAccumulator&)"
+        )   << "Attempted assignment to self"
+            << abort(FatalError);
+    }
+
+    List< Field<Type> >::operator=(rhs);
+
+    averagesTaken_ = rhs.averagesTaken();
+
+    bufferOffsets_ = rhs.bufferOffsets();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#   include "bufferedAccumulatorIO.C"
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecularMeasurements/bufferedAccumulator/bufferedAccumulator.H b/src/lagrangian/molecularDynamics/molecularMeasurements/bufferedAccumulator/bufferedAccumulator.H
new file mode 100644
index 0000000000..0311ad3096
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecularMeasurements/bufferedAccumulator/bufferedAccumulator.H
@@ -0,0 +1,176 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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::bufferedAccumulator
+
+Description
+
+SourceFiles
+    bufferedAccumulatorI.H
+    bufferedAccumulator.C
+    bufferedAccumulatorIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef bufferedAccumulator_H
+#define bufferedAccumulator_H
+
+#include "Field.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class Type>
+class bufferedAccumulator;
+
+template<class Type>
+Ostream& operator<<
+(
+    Ostream&,
+    const bufferedAccumulator<Type>&
+);
+
+/*---------------------------------------------------------------------------*\
+                      Class bufferedAccumulator Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class bufferedAccumulator
+:
+    public List< Field<Type> >
+{
+    // Private data
+
+        label averagesTaken_;
+
+        List<label> bufferOffsets_;
+
+
+    // Private Member Functions
+
+        inline Field<Type>& accumulationBuffer();
+
+        inline const Field<Type>& accumulationBuffer() const;
+
+        void accumulateAndResetBuffer(const label b);
+
+
+public:
+
+    //- Component type
+    typedef typename pTraits<Type>::cmptType cmptType;
+
+
+    // Static data members
+
+        static const char* const typeName;
+
+
+    // Constructors
+
+        //- Construct null
+        bufferedAccumulator();
+
+        //- Construct from components
+        bufferedAccumulator
+        (
+            const label nBuffers,
+            const label bufferLength,
+            const label bufferingInterval
+        );
+
+        //- Construct as copy
+        bufferedAccumulator(const bufferedAccumulator<Type>&);
+
+
+    //- Destructor
+    ~bufferedAccumulator();
+
+
+    // Member Functions
+
+        label addToBuffers(const List<Type>& valuesToAdd);
+
+        Field<Type> averaged() const;
+
+        void resetAveraging();
+
+
+        // Access
+
+            inline label averagesTaken() const;
+
+            inline label nBuffers() const;
+
+            inline label bufferLength() const;
+
+            inline const List<label>& bufferOffsets() const;
+
+
+        // Edit
+
+            void setSizes
+            (
+                const label nBuffers,
+                const label bufferLength,
+                const label bufferingInterval
+            );
+
+
+    // Member Operators
+
+        void operator=(const bufferedAccumulator<Type>&);
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<< <Type>
+        (
+            Ostream&,
+            const bufferedAccumulator<Type>&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "bufferedAccumulatorI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "bufferedAccumulator.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecularMeasurements/bufferedAccumulator/bufferedAccumulatorI.H b/src/lagrangian/molecularDynamics/molecularMeasurements/bufferedAccumulator/bufferedAccumulatorI.H
new file mode 100644
index 0000000000..42d0cec63d
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecularMeasurements/bufferedAccumulator/bufferedAccumulatorI.H
@@ -0,0 +1,79 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+inline Field<Type>& bufferedAccumulator<Type>::accumulationBuffer()
+{
+    return (*this)[nBuffers()];
+}
+
+
+template<class Type>
+inline const Field<Type>& bufferedAccumulator<Type>::accumulationBuffer() const
+{
+    return (*this)[nBuffers()];
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+inline label bufferedAccumulator<Type>::averagesTaken() const
+{
+    return averagesTaken_;
+}
+
+
+template<class Type>
+inline label bufferedAccumulator<Type>::nBuffers() const
+{
+    return bufferOffsets_.size();
+}
+
+
+template<class Type>
+inline label bufferedAccumulator<Type>::bufferLength() const
+{
+    return (*this)[0].size();
+}
+
+
+template<class Type>
+inline const List<label>& bufferedAccumulator<Type>::bufferOffsets() const
+{
+    return bufferOffsets_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecularMeasurements/bufferedAccumulator/bufferedAccumulatorIO.C b/src/lagrangian/molecularDynamics/molecularMeasurements/bufferedAccumulator/bufferedAccumulatorIO.C
new file mode 100644
index 0000000000..1e891fcb36
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecularMeasurements/bufferedAccumulator/bufferedAccumulatorIO.C
@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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 "bufferedAccumulator.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+template<class Type>
+Foam::Ostream&
+Foam::operator<<(Ostream& os, const bufferedAccumulator<Type>& bA)
+{
+
+    os  << bA.averagesTaken_
+        << static_cast<const List< Field<Type> >&>(bA)
+        << bA.bufferOffsets();
+
+    // Check state of Ostream
+    os.check
+    (
+        "Foam::Ostream& Foam::operator<<(Foam::Ostream&, "
+        "const Foam::bufferedAccumulator&)"
+    );
+
+    return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunction.C b/src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunction.C
new file mode 100644
index 0000000000..1eedef0e4c
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunction.C
@@ -0,0 +1,223 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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 "correlationFunction.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<class Type>
+const char* const
+    Foam::correlationFunction<Type>::typeName("correlationFunction");
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+void Foam::correlationFunction<Type>::setTimesAndSizes
+(
+    const label tZeroBufferSize
+)
+{
+    sampleSteps_  = ceil(sampleInterval_/mesh_.time().deltaTValue());
+
+    sampleInterval_ = sampleSteps_*mesh_.time().deltaTValue();
+
+    label bufferLength(ceil(duration_/sampleInterval_));
+
+    duration_ = bufferLength*sampleInterval_;
+
+    label bufferingInterval(ceil(averagingInterval_/sampleInterval_));
+
+    averagingInterval_ = bufferingInterval*sampleInterval_;
+
+    label nBuffers(ceil(duration_/averagingInterval_));
+
+    this->setSizes
+    (
+        nBuffers,
+        bufferLength,
+        bufferingInterval
+    );
+
+    tZeroBuffers_ =
+        Field< Field<Type> >
+        (
+            nBuffers,
+            Field<Type>
+            (
+                tZeroBufferSize,
+                pTraits<Type>::zero
+            )
+        );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::correlationFunction<Type>::correlationFunction
+(
+    const polyMesh& mesh,
+    const dictionary& cfDict,
+    const label tZeroBufferSize
+)
+:
+    bufferedAccumulator<scalar>(),
+    mesh_(mesh)
+{
+    duration_ = readScalar(cfDict.lookup("duration"));
+
+    sampleInterval_ = readScalar(cfDict.lookup("sampleInterval"));
+
+    averagingInterval_ = readScalar(cfDict.lookup("averagingInterval"));
+
+    setTimesAndSizes(tZeroBufferSize);
+}
+
+
+template<class Type>
+Foam::correlationFunction<Type>::correlationFunction
+(
+    const polyMesh& mesh,
+    const label tZeroBufferSize,
+    const scalar duration,
+    const scalar sampleInterval,
+    const scalar averagingInterval
+)
+:
+    bufferedAccumulator<scalar>(),
+    mesh_(mesh),
+    duration_(duration),
+    sampleInterval_(sampleInterval),
+    averagingInterval_(averagingInterval)
+{
+    setTimesAndSizes(tZeroBufferSize);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::correlationFunction<Type>::~correlationFunction()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::correlationFunction<Type>::calculateCorrelationFunction
+(
+    const Field<Type>& currentValues
+)
+{
+    if (measurandFieldSize() != currentValues.size())
+    {
+        FatalErrorIn("correlationFunction<Type>::calculateCorrelationFunction")
+            << "Trying to supply a Field of length"
+            << currentValues.size()
+            << " to calculate the correlation function. "
+            << "Expecting a Field of length "
+            << measurandFieldSize() << nl
+            << abort(FatalError);
+    }
+
+    List<scalar> cFSums(nBuffers(),0.0);
+
+    forAll(tZeroBuffers_, tZB)
+    {
+        scalar& cFSum = cFSums[tZB];
+
+        const Field<Type>& tZeroBuffer = tZeroBuffers_[tZB];
+
+        forAll(currentValues, cV)
+        {
+            const Type& tZeroBufferValue = tZeroBuffer[cV];
+
+            const Type& currentValue = currentValues[cV];
+
+            forAll(currentValue, component)
+            {
+                cFSum +=
+                (
+                    tZeroBufferValue[component]*currentValue[component]
+                );
+            }
+        }
+
+        cFSum /= (measurandFieldSize()*currentValues[0].size());
+    }
+
+    label bufferToRefill = addToBuffers(cFSums);
+
+    if (bufferToRefill != -1)
+    {
+        tZeroBuffers_[bufferToRefill] = currentValues;
+    }
+}
+
+
+template<class Type>
+void Foam::correlationFunction<Type>::calculateCorrelationFunction
+(
+    const Type& currentValue
+)
+{
+    if (measurandFieldSize() != 1)
+    {
+        FatalErrorIn("correlationFunction<Type>::calculateCorrelationFunction")
+            << "Trying to supply a single value to calculate the correlation "
+            << "function.  Expecting a Field of length "
+            << measurandFieldSize()
+            << abort(FatalError);
+    }
+
+    calculateCorrelationFunction(Field<Type>(1, currentValue));
+}
+
+
+template<class Type>
+Foam::scalar Foam::correlationFunction<Type>::integral() const
+{
+    Field<scalar> averageCF(averaged());
+
+    scalar cFIntegral = 0.0;
+
+    for (label v = 0; v < averageCF.size() - 1; v++)
+    {
+        cFIntegral +=
+            0.5
+           *sampleInterval_
+           *(averageCF[v+1] + averageCF[v]);
+    }
+
+    return cFIntegral;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#   include "correlationFunctionIO.C"
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunction.H b/src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunction.H
new file mode 100644
index 0000000000..f0c26f9b1b
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunction.H
@@ -0,0 +1,181 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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::correlationFunction
+
+Description
+
+SourceFiles
+    correlationFunctionI.H
+    correlationFunction.C
+    correlationFunctionIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef correlationFunction_H
+#define correlationFunction_H
+
+#include "bufferedAccumulator.H"
+#include "dictionary.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class Type>
+class correlationFunction;
+
+template<class Type>
+Ostream& operator<<
+(
+    Ostream&,
+    const correlationFunction<Type>&
+);
+
+/*---------------------------------------------------------------------------*\
+                    Class correlationFunction Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class correlationFunction
+:
+    public bufferedAccumulator<scalar>
+{
+    // Private data
+
+        const polyMesh& mesh_;
+
+        Field< Field<Type> > tZeroBuffers_;
+
+        scalar duration_;
+        scalar sampleInterval_;
+        scalar averagingInterval_;
+
+        label sampleSteps_;
+
+
+    // Private Member Functions
+
+        void setTimesAndSizes(const label);
+
+        //- Disallow default bitwise copy construct
+        correlationFunction(const correlationFunction<Type>&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const correlationFunction<Type>&);
+
+
+public:
+
+    //- Component type
+    typedef typename pTraits<Type>::cmptType cmptType;
+
+
+    // Static data members
+
+        static const char* const typeName;
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        correlationFunction
+        (
+            const polyMesh&,
+            const dictionary&,
+            const label tZeroBufferSize
+        );
+
+        //- Construct from components
+        correlationFunction
+        (
+            const polyMesh&,
+            const label tZeroBufferSize,
+            const scalar duration,
+            const scalar sampleInterval,
+            const scalar averagingInterval
+
+        );
+
+
+    //- Destructor
+    ~correlationFunction();
+
+
+    // Member Functions
+
+        void calculateCorrelationFunction(const Field<Type>&);
+
+        void calculateCorrelationFunction(const Type&);
+
+        scalar integral() const;
+
+        bool writeAveraged(Ostream&) const;
+
+
+    // Access
+
+        inline const Field< Field<Type> >& tZeroBuffers() const;
+
+        inline scalar duration() const;
+
+        inline scalar sampleInterval() const;
+
+        inline scalar averagingInterval() const;
+
+        inline label sampleSteps() const;
+
+        inline label measurandFieldSize() const;
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<< <Type>
+        (
+            Ostream&,
+            const correlationFunction<Type>&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "correlationFunctionI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "correlationFunction.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunctionI.H b/src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunctionI.H
new file mode 100644
index 0000000000..752facb6f9
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunctionI.H
@@ -0,0 +1,69 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+inline const Foam::Field< Foam::Field<Type> >&
+Foam::correlationFunction<Type>::tZeroBuffers() const
+{
+    return tZeroBuffers_;
+}
+
+
+template<class Type>
+inline Foam::scalar Foam::correlationFunction<Type>::duration() const
+{
+    return duration_;
+}
+
+
+template<class Type>
+inline Foam::scalar Foam::correlationFunction<Type>::sampleInterval() const
+{
+    return sampleInterval_;
+}
+
+
+template<class Type>
+inline Foam::scalar Foam::correlationFunction<Type>::averagingInterval() const
+{
+    return averagingInterval_;
+}
+
+
+template<class Type>
+inline Foam::label Foam::correlationFunction<Type>::sampleSteps() const
+{
+    return sampleSteps_;
+}
+
+
+template<class Type>
+inline Foam::label Foam::correlationFunction<Type>::measurandFieldSize() const
+{
+    return tZeroBuffers_[0].size();
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunctionIO.C b/src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunctionIO.C
new file mode 100644
index 0000000000..48b2d74228
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunctionIO.C
@@ -0,0 +1,71 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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 "correlationFunction.H"
+#include "IOstreams.H"
+
+template<class Type>
+bool Foam::correlationFunction<Type>::writeAveraged(Ostream& os) const
+{
+    Field<scalar> averageCF(averaged());
+
+    forAll(averageCF, v)
+    {
+        os  << v*sampleInterval()
+            << token::SPACE
+            << averageCF[v]
+            << nl;
+    }
+
+    return os.good();
+}
+
+
+template<class Type>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const correlationFunction<Type>& cF
+)
+{
+    os  << cF.duration()
+        << nl << cF.sampleInterval()
+        << nl << cF.averagingInterval()
+        << nl << cF.sampleSteps()
+        << nl << cF.tZeroBuffers()
+        << nl << static_cast<const bufferedAccumulator<scalar>&>(cF);
+
+    // Check state of Ostream
+    os.check
+    (
+        "Foam::Ostream& Foam::operator<<"
+        "(Ostream&, const correlationFunction<Type>&)"
+    );
+
+    return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.C b/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.C
new file mode 100644
index 0000000000..1110415d10
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.C
@@ -0,0 +1,478 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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 "distribution.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(distribution, 0);
+}
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+void Foam::distribution::write
+(
+    const fileName& file,
+    const List<Pair<scalar> >& pairs
+)
+{
+    OFstream os(file);
+
+    forAll(pairs, i)
+    {
+        os  << pairs[i].first() << ' ' << pairs[i].second() << nl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::distribution::distribution()
+:
+    Map<label>(),
+    binWidth_(1)
+{}
+
+
+Foam::distribution::distribution(const scalar binWidth)
+:
+    Map<label>(),
+    binWidth_(binWidth)
+{}
+
+
+Foam::distribution::distribution(const distribution& d)
+:
+    Map<label>(static_cast< Map<label> >(d)),
+    binWidth_(d.binWidth())
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::distribution::~distribution()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::label Foam::distribution::totalEntries() const
+{
+    label sumOfEntries = 0;
+
+    forAllConstIter(Map<label>, *this, iter)
+    {
+        sumOfEntries += iter();
+
+        if (sumOfEntries < 0)
+        {
+            WarningIn("label distribution::totalEntries()")
+                << "Accumulated distribution values total has become negative: "
+                << "sumOfEntries = " << sumOfEntries
+                << ". This is most likely to be because too many samples "
+                << "have been added to the bins and the label has 'rolled "
+                << "round'. Try distribution::approxTotalEntries which "
+                << "returns a scalar." << endl;
+
+            sumOfEntries = -1;
+
+            break;
+        }
+    }
+
+    return sumOfEntries;
+}
+
+
+Foam::scalar Foam::distribution::approxTotalEntries() const
+{
+    scalar sumOfEntries = 0;
+
+    forAllConstIter(Map<label>, *this, iter)
+    {
+        sumOfEntries += scalar(iter());
+    }
+
+    return sumOfEntries;
+}
+
+
+Foam::scalar Foam::distribution::mean() const
+{
+    scalar runningSum = 0;
+
+    scalar totEnt = approxTotalEntries();
+
+    List<label> keys = toc();
+
+    forAll(keys,k)
+    {
+        label key = keys[k];
+
+        runningSum +=
+            (0.5 + scalar(key))
+           *binWidth_
+           *scalar((*this)[key])
+           /totEnt;
+    }
+
+    return runningSum;
+}
+
+
+Foam::scalar Foam::distribution::median()
+{
+    // From:
+    // http://mathworld.wolfram.com/StatisticalMedian.html
+    // The statistical median is the value of the distribution variable
+    // where the cumulative distribution = 0.5.
+
+    scalar median = 0.0;
+
+    scalar runningSum = 0.0;
+
+    List<Pair<scalar> > normDist(normalised());
+
+    if (normDist.size())
+    {
+        if (normDist.size() == 1)
+        {
+            median = normDist[0].first();
+        }
+        else if
+        (
+            normDist.size() > 1
+         && normDist[0].second()*binWidth_ > 0.5
+        )
+        {
+            scalar xk = normDist[1].first();
+            scalar xkm1 = normDist[0].first();
+            scalar Sk =
+                (normDist[0].second() + normDist[1].second())*binWidth_;
+            scalar Skm1 = normDist[0].second()*binWidth_;
+
+            median = (0.5 - Skm1)*(xk - xkm1)/(Sk - Skm1) + xkm1;
+        }
+        else
+        {
+            label lastNonZeroIndex = 0;
+
+            forAll(normDist,nD)
+            {
+                if (runningSum + (normDist[nD].second()*binWidth_) > 0.5)
+                {
+                    scalar xk = normDist[nD].first();
+                    scalar xkm1 = normDist[lastNonZeroIndex].first();
+                    scalar Sk = runningSum + (normDist[nD].second()*binWidth_);
+                    scalar Skm1 = runningSum;
+
+                    median = (0.5 - Skm1)*(xk - xkm1)/(Sk - Skm1) + xkm1;
+
+                    break;
+                }
+                else if (normDist[nD].second() > 0.0)
+                {
+                    runningSum += normDist[nD].second()*binWidth_;
+
+                    lastNonZeroIndex = nD;
+                }
+            }
+        }
+    }
+
+    return median;
+}
+
+
+void Foam::distribution::add(const scalar valueToAdd)
+{
+    iterator iter(this->begin());
+
+    label n = label(valueToAdd/binWidth_) - label(neg(valueToAdd/binWidth_));
+
+    iter = find(n);
+
+    if (iter == this->end())
+    {
+        this->insert(n,1);
+    }
+    else
+    {
+        (*this)[n]++;
+    }
+
+    if ((*this)[n] < 0)
+    {
+        FatalErrorIn("distribution::add(const scalar valueToAdd)")
+            << "Accumulated distribution value has become negative: "
+            << "bin = " << (0.5 + scalar(n)) * binWidth_
+            << ", value = " << (*this)[n]
+            << ". This is most likely to be because too many samples "
+            << "have been added to a bin and the label has 'rolled round'"
+            << abort(FatalError);
+    }
+}
+
+
+void Foam::distribution::add(const label valueToAdd)
+{
+    add(scalar(valueToAdd));
+}
+
+
+void Foam::distribution::insertMissingKeys()
+{
+    iterator iter(this->begin());
+
+    List<label> keys = toc();
+
+    sort(keys);
+
+    if (keys.size())
+    {
+        for (label k = keys[0]; k < keys.last(); k++)
+        {
+            iter = find(k);
+
+            if (iter == this->end())
+            {
+                this->insert(k,0);
+            }
+        }
+    }
+}
+
+
+Foam::List<Foam::Pair<Foam::scalar> > Foam::distribution::normalised()
+{
+    scalar totEnt = approxTotalEntries();
+
+    insertMissingKeys();
+
+    List<label> keys = toc();
+
+    sort(keys);
+
+    List<Pair<scalar> > normDist(size());
+
+    forAll(keys,k)
+    {
+        label key = keys[k];
+
+        normDist[k].first() = (0.5 + scalar(key))*binWidth_;
+
+        normDist[k].second() = scalar((*this)[key])/totEnt/binWidth_;
+    }
+
+    if (debug)
+    {
+        Info<< "totEnt: " << totEnt << endl;
+    }
+
+    return normDist;
+}
+
+
+Foam::List<Foam::Pair<Foam::scalar> > Foam::distribution::normalisedMinusMean()
+{
+    return normalisedShifted(mean());
+}
+
+
+Foam::List<Foam::Pair<Foam::scalar> > Foam::distribution::normalisedShifted
+(
+    scalar shiftValue
+)
+{
+    List<Pair<scalar> > oldDist(normalised());
+
+    List<Pair<scalar> > newDist(oldDist.size());
+
+    forAll(oldDist,u)
+    {
+        oldDist[u].first() -= shiftValue;
+    }
+
+    scalar lowestOldBin = oldDist[0].first()/binWidth_ - 0.5;
+
+    label lowestNewKey = label
+    (
+        lowestOldBin + 0.5*sign(lowestOldBin)
+    );
+
+    scalar interpolationStartDirection =
+        sign(scalar(lowestNewKey) - lowestOldBin);
+
+    label newKey = lowestNewKey;
+
+    if (debug)
+    {
+        Info<< shiftValue
+            << nl << lowestOldBin
+            << nl << lowestNewKey
+            << nl << interpolationStartDirection
+            << endl;
+
+        scalar checkNormalisation = 0;
+
+        forAll(oldDist, oD)
+        {
+            checkNormalisation += oldDist[oD].second()*binWidth_;
+        }
+
+        Info<< "Initial normalisation = " << checkNormalisation << endl;
+    }
+
+    forAll(oldDist,u)
+    {
+        newDist[u].first() = (0.5 + scalar(newKey)) * binWidth_;
+
+        if (interpolationStartDirection < 0)
+        {
+            if (u == 0)
+            {
+                newDist[u].second() =
+                    (0.5 + scalar(newKey))*oldDist[u].second()
+                  - oldDist[u].second()
+                        *(oldDist[u].first() - binWidth_)/ binWidth_;
+            }
+            else
+            {
+                newDist[u].second() =
+                    (0.5 + scalar(newKey))
+                   *(oldDist[u].second() - oldDist[u-1].second())
+                  +
+                    (
+                        oldDist[u-1].second()*oldDist[u].first()
+                      - oldDist[u].second()*oldDist[u-1].first()
+                    )
+                    /binWidth_;
+            }
+        }
+        else
+        {
+            if (u == oldDist.size() - 1)
+            {
+                newDist[u].second() =
+                    (0.5 + scalar(newKey))*-oldDist[u].second()
+                  + oldDist[u].second()*(oldDist[u].first() + binWidth_)
+                   /binWidth_;
+            }
+            else
+            {
+                newDist[u].second() =
+                    (0.5 + scalar(newKey))
+                   *(oldDist[u+1].second() - oldDist[u].second())
+                  +
+                    (
+                        oldDist[u].second()*oldDist[u+1].first()
+                      - oldDist[u+1].second()*oldDist[u].first()
+                    )
+                   /binWidth_;
+            }
+        }
+
+        newKey++;
+    }
+
+    if (debug)
+    {
+        scalar checkNormalisation = 0;
+
+        forAll(newDist, nD)
+        {
+            checkNormalisation += newDist[nD].second()*binWidth_;
+        }
+
+        Info<< "Shifted normalisation = " << checkNormalisation << endl;
+    }
+
+    return newDist;
+}
+
+
+Foam::List<Foam::Pair<Foam::scalar> > Foam::distribution::raw()
+{
+    insertMissingKeys();
+
+    List<label> keys = toc();
+
+    sort(keys);
+
+    List<Pair<scalar> > rawDist(size());
+
+    forAll(keys,k)
+    {
+        label key = keys[k];
+
+        rawDist[k].first() = (0.5 + scalar(key))*binWidth_;
+
+        rawDist[k].second() = scalar((*this)[key]);
+    }
+
+    return rawDist;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+void Foam::distribution::operator=(const distribution& rhs)
+{
+    // Check for assignment to self
+    if (this == &rhs)
+    {
+        FatalErrorIn("distribution::operator=(const distribution&)")
+            << "Attempted assignment to self"
+            << abort(FatalError);
+    }
+
+    Map<label>::operator=(rhs);
+
+    binWidth_ = rhs.binWidth();
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const distribution& d)
+{
+    os  << d.binWidth_
+        << static_cast<const Map<label>&>(d);
+
+    // Check state of Ostream
+    os.check
+    (
+        "Ostream& operator<<(Ostream&, "
+        "const distribution&)"
+    );
+
+    return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.H b/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.H
new file mode 100644
index 0000000000..02f347d170
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.H
@@ -0,0 +1,148 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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::distribution
+
+Description
+    Accumulating histogram of values.  Specified bin resolution
+    automatic generation of bins.
+
+SourceFiles
+    distributionI.H
+    distribution.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef distribution_H
+#define distribution_H
+
+#include "Map.H"
+#include "Pair.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class distribution Declaration
+\*---------------------------------------------------------------------------*/
+
+class distribution
+:
+    public Map<label>
+{
+    // Private data
+
+        scalar binWidth_;
+
+
+public:
+
+    //- Runtime type information
+
+        TypeName("distribution");
+
+    // Static functions
+
+        //- write to file
+
+            static void write
+            (
+                const fileName& file,
+                const List<Pair<scalar> >& pairs
+            );
+
+
+    // Constructors
+
+        //- Construct null
+        distribution();
+
+        //- Construct from binWidth
+        distribution(const scalar binWidth);
+
+        //- Construct as copy
+        distribution(const distribution&);
+
+
+    //- Destructor
+    virtual ~distribution();
+
+
+    // Member Functions
+
+        label totalEntries() const;
+
+        scalar approxTotalEntries() const;
+
+        scalar mean() const;
+
+        scalar median();
+
+        //- Add a value to the appropriate bin of the distribution.
+        void add(const scalar valueToAdd);
+
+        void add(const label valueToAdd);
+
+        void insertMissingKeys();
+
+        List<Pair<scalar> > normalised();
+
+        List<Pair<scalar> > normalisedMinusMean();
+
+        List<Pair<scalar> > normalisedShifted(scalar shiftValue);
+
+        List<Pair<scalar> > raw();
+
+
+        // Access
+
+            inline scalar binWidth() const;
+
+
+    // Member Operators
+
+        void operator=(const distribution&);
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<<(Ostream&, const distribution&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "distributionI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distributionI.H b/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distributionI.H
new file mode 100644
index 0000000000..344567114b
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distributionI.H
@@ -0,0 +1,34 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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  * * * * * * * * * * * * * //
+
+inline Foam::scalar Foam::distribution::binWidth() const
+{
+    return binWidth_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/Make/files b/src/lagrangian/molecularDynamics/molecule/Make/files
index e71272b545..9086ba7ffa 100644
--- a/src/lagrangian/molecularDynamics/molecule/Make/files
+++ b/src/lagrangian/molecularDynamics/molecule/Make/files
@@ -1,14 +1,9 @@
-clouds/baseClasses/moleculeCloud/moleculeCloud.C
+reducedUnits/reducedUnits.C
+reducedUnits/reducedUnitsIO.C
 
-molecules/constPropSite/constPropSite.C
+molecule/molecule.C
+molecule/moleculeIO.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 */
+moleculeCloud/moleculeCloud.C
 
 LIB = $(FOAM_LIBBIN)/libmolecule
diff --git a/src/lagrangian/molecularDynamics/molecule/Make/options b/src/lagrangian/molecularDynamics/molecule/Make/options
index 78ecf2089f..5698f81a65 100644
--- a/src/lagrangian/molecularDynamics/molecule/Make/options
+++ b/src/lagrangian/molecularDynamics/molecule/Make/options
@@ -2,11 +2,13 @@ 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/potential/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
     -llagrangian \
-    -lpotential
+    -lpotential \
+    -lmolecularMeasurements
 
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/averageMDFields.H b/src/lagrangian/molecularDynamics/molecule/mdTools/averageMDFields.H
new file mode 100644
index 0000000000..215d6bbef8
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/averageMDFields.H
@@ -0,0 +1,222 @@
+if (runTime.outputTime())
+{
+    /*-----------------------------------------------------------------------*\
+        Number density
+    \*-----------------------------------------------------------------------*/
+
+    scalarField totalRhoN_sum(mesh.nCells(), 0.0);
+
+    forAll(allSpeciesRhoN, rN)
+    {
+        allSpeciesRhoN[rN].internalField() =
+            allSpeciesN_RU[rN]
+            /mesh.cellVolumes()
+            /nAveragingSteps;
+
+        totalRhoN_sum += allSpeciesRhoN[rN].internalField();
+    }
+
+    totalRhoN.internalField() = totalRhoN_sum;
+
+
+    /*-----------------------------------------------------------------------*\
+        Mass density
+    \*-----------------------------------------------------------------------*/
+
+    scalarField totalRhoM_sum(mesh.nCells(), 0.0);
+
+    forAll(allSpeciesRhoM, rM)
+    {
+        allSpeciesRhoM[rM].internalField() =
+            allSpeciesM_RU[rM]
+            /mesh.cellVolumes()
+            /nAveragingSteps;
+
+        totalRhoM_sum += allSpeciesRhoM[rM].internalField();
+    }
+
+    totalRhoM.internalField() = totalRhoM_sum;
+
+    /*-----------------------------------------------------------------------*\
+        Bulk velocity
+    \*-----------------------------------------------------------------------*/
+
+    vectorField totalMomentum_sum(mesh.nCells(), vector::zero);
+
+    scalarField totalMass_sum(mesh.nCells(), 0.0);
+
+    forAll(allSpeciesVelocity, v)
+    {
+        // A check for 1/0 molecules is required.
+
+        vectorField& singleSpeciesVelocity
+        (
+            allSpeciesVelocity[v].internalField()
+        );
+
+        forAll(singleSpeciesVelocity, sSV)
+        {
+            if (allSpeciesN_RU[v][sSV])
+            {
+                singleSpeciesVelocity[sSV] =
+                    allSpeciesVelocitySum_RU[v][sSV]
+                    /allSpeciesN_RU[v][sSV];
+
+                totalMomentum_sum[sSV] +=
+                    allSpeciesM_RU[v][sSV]
+                    /allSpeciesN_RU[v][sSV]
+                    *allSpeciesVelocitySum_RU[v][sSV];
+
+                totalMass_sum[sSV] += allSpeciesM_RU[v][sSV];
+            }
+            else
+            {
+                singleSpeciesVelocity[sSV] = vector::zero;
+            }
+        }
+    }
+
+    forAll(totalVelocity.internalField(), tV)
+    {
+        if (totalMass_sum[tV] > VSMALL)
+        {
+            totalVelocity.internalField()[tV] =
+                totalMomentum_sum[tV]
+                /totalMass_sum[tV];
+        }
+        else
+        {
+            totalVelocity.internalField()[tV] =
+            vector::zero;
+        }
+    }
+
+    /*-----------------------------------------------------------------------*\
+        Kinetic temperature
+    \*-----------------------------------------------------------------------*/
+
+    scalarField totalTemperatureVTerms_sum(mesh.nCells(), 0.0);
+
+    scalarField totalN_sum(mesh.nCells(), 0.0);
+
+    forAll(allSpeciesTemperature, t)
+    {
+        // A check for 1/0 molecules is required.
+
+        scalarField& singleSpeciesTemp
+        (
+            allSpeciesTemperature[t].internalField()
+        );
+
+        forAll(singleSpeciesTemp, sST)
+        {
+            if (allSpeciesN_RU[t][sST])
+            {
+                singleSpeciesTemp[sST] =
+                    allSpeciesM_RU[t][sST]
+                    /allSpeciesN_RU[t][sST]
+                    /(3.0 * moleculeCloud::kb * allSpeciesN_RU[t][sST])
+                   *(
+                        allSpeciesVelocityMagSquaredSum_RU[t][sST]
+                        -
+                        (
+                            allSpeciesVelocitySum_RU[t][sST]
+                            &
+                            allSpeciesVelocitySum_RU[t][sST]
+                        )
+                        /allSpeciesN_RU[t][sST]
+                    );
+
+                totalTemperatureVTerms_sum[sST] +=
+                    allSpeciesM_RU[t][sST]
+                   /allSpeciesN_RU[t][sST]
+                   *(
+                        allSpeciesVelocityMagSquaredSum_RU[t][sST]
+                      -
+                        (
+                            allSpeciesVelocitySum_RU[t][sST]
+                            &
+                            allSpeciesVelocitySum_RU[t][sST]
+                        )
+                        /allSpeciesN_RU[t][sST]
+                    );
+
+                totalN_sum[sST] += allSpeciesN_RU[t][sST];
+            }
+            else
+            {
+                singleSpeciesTemp[sST] = 0.0;
+            }
+        }
+    }
+
+    forAll(totalTemperature.internalField(), tT)
+    {
+        if (totalN_sum[tT] > 0)
+        {
+            totalTemperature.internalField()[tT] =
+                totalTemperatureVTerms_sum[tT]
+                /(3.0 * moleculeCloud::kb * totalN_sum[tT]);
+        }
+        else
+        {
+            totalTemperature.internalField()[tT] = 0.0;
+        }
+    }
+
+    /*-----------------------------------------------------------------------*\
+        Mean kinetic energy
+    \*-----------------------------------------------------------------------*/
+
+    scalarField totalKE_sum(mesh.nCells(), 0.0);
+
+    forAll(allSpeciesMeanKE, mKE)
+    {
+        // A check for 1/0 molecules is required.
+
+        scalarField& singleSpeciesMeanKE
+        (
+            allSpeciesMeanKE[mKE].internalField()
+        );
+
+        forAll(singleSpeciesMeanKE, sSMKE)
+        {
+            if (allSpeciesN_RU[mKE][sSMKE])
+            {
+                singleSpeciesMeanKE[sSMKE] =
+                    allSpeciesM_RU[mKE][sSMKE]
+                   /allSpeciesN_RU[mKE][sSMKE]
+                   /(2.0*allSpeciesN_RU[mKE][sSMKE])
+                   *(
+                        allSpeciesVelocityMagSquaredSum_RU[mKE][sSMKE]
+                    );
+
+                totalKE_sum[sSMKE] +=
+                    allSpeciesM_RU[mKE][sSMKE]
+                    /allSpeciesN_RU[mKE][sSMKE]
+                    /2.0
+                   *(
+                        allSpeciesVelocityMagSquaredSum_RU[mKE][sSMKE]
+                    );
+            }
+            else
+            {
+                singleSpeciesMeanKE[sSMKE] = 0.0;
+            }
+        }
+    }
+
+    forAll(totalMeanKE.internalField(), tMKE)
+    {
+        if (totalN_sum[tMKE] > 0)
+        {
+            totalMeanKE.internalField()[tMKE] =
+                totalKE_sum[tMKE]
+                /totalN_sum[tMKE];
+        }
+        else
+        {
+            totalMeanKE.internalField()[tMKE] = 0.0;
+        }
+    }
+}
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/calculateAutoCorrelationFunctions.H b/src/lagrangian/molecularDynamics/molecule/mdTools/calculateAutoCorrelationFunctions.H
new file mode 100644
index 0000000000..eda25d00f7
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/calculateAutoCorrelationFunctions.H
@@ -0,0 +1,77 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+if (mesh.time().timeIndex() % vacf.sampleSteps() == 0)
+{
+    Field<vector> uVals(molecules.size());
+
+    label uV = 0;
+
+    forAllConstIter(IDLList<molecule>, molecules, mol)
+    {
+        uVals[uV++] = mol().U();
+    }
+
+    vacf.calculateCorrelationFunction(uVals);
+}
+
+if (mesh.time().timeIndex() % pacf.sampleSteps() == 0)
+{
+    vector p = vector::zero;
+
+    forAllConstIter(IDLList<molecule>, molecules, mol)
+    {
+        p.x() +=
+            mol().mass() * mol().U().y() * mol().U().z()
+          + 0.5*mol().rf().yz();
+
+        p.y() +=
+            mol().mass() * mol().U().z() * mol().U().x()
+          + 0.5*mol().rf().zx();
+
+        p.z() +=
+            mol().mass() * mol().U().x() * mol().U().y()
+          + 0.5*mol().rf().xy();
+    }
+
+    pacf.calculateCorrelationFunction(p);
+}
+
+if (mesh.time().timeIndex() % hfacf.sampleSteps() == 0)
+{
+    vector s = vector::zero;
+
+    forAllConstIter(IDLList<molecule>, molecules, mol)
+    {
+        s +=
+        (
+            0.5*mol().mass()*magSqr(mol().U())
+          + mol().potentialEnergy()
+        )*mol().U()
+      + 0.5*(mol().rf() & mol().U());
+    }
+
+    hfacf.calculateCorrelationFunction(s);
+}
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/calculateMDFields.H b/src/lagrangian/molecularDynamics/molecule/mdTools/calculateMDFields.H
new file mode 100644
index 0000000000..af5c938da4
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/calculateMDFields.H
@@ -0,0 +1,23 @@
+const List<DynamicList<molecule*> >& cellOccupancy = molecules.cellOccupancy();
+
+forAll(cellOccupancy, cell)
+{
+    const List<molecule*>& molsInCell = cellOccupancy[cell];
+
+    forAll(molsInCell, mIC)
+    {
+        molecule* mol = molsInCell[mIC];
+
+        const label molId = mol->id();
+
+        const vector& molU = mol->U();
+
+        allSpeciesN_RU[molId][cell]++;
+
+        allSpeciesM_RU[molId][cell] += mol->mass();
+
+        allSpeciesVelocitySum_RU[molId][cell] += molU;
+
+        allSpeciesVelocityMagSquaredSum_RU[molId][cell] += molU & molU;
+    }
+}
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/calculateTransportProperties.H b/src/lagrangian/molecularDynamics/molecule/mdTools/calculateTransportProperties.H
new file mode 100644
index 0000000000..0311e077dd
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/calculateTransportProperties.H
@@ -0,0 +1,81 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+if (writeVacf)
+{
+    OFstream vacfFile(runTime.path()/"vacf");
+
+    if (!vacf.writeAveraged(vacfFile))
+    {
+        FatalErrorIn(args.executable())
+            << "Failed writing to "
+            << vacfFile.name()
+            << abort(FatalError);
+    }
+}
+
+Info<< "Diffusion coefficient = "
+    << vacf.integral() << endl;
+
+if (writePacf)
+{
+    OFstream pacfFile(runTime.path()/"pacf");
+
+    if (!pacf.writeAveraged(pacfFile))
+    {
+        FatalErrorIn(args.executable())
+            << "Failed writing to "
+            << pacfFile.name()
+            << abort(FatalError);
+    }
+}
+
+Info<< "Viscosity = "
+    << pacf.integral()/averageTemperature/moleculeCloud::kb/meshVolume
+    << endl;
+
+if (writeHFacf)
+{
+    OFstream hfacfFile
+    (
+        runTime.path()/ + "hfacf"
+    );
+
+    if (!hfacf.writeAveraged(hfacfFile))
+    {
+        FatalErrorIn(args.executable())
+            << "Failed writing to "
+            << hfacfFile.name()
+            << abort(FatalError);
+    }
+}
+
+Info<< "Thermal conductivity = "
+    << hfacf.integral()
+        /averageTemperature
+        /averageTemperature
+        /moleculeCloud::kb
+        / meshVolume
+    << endl;
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/createAutoCorrelationFunctions.H b/src/lagrangian/molecularDynamics/molecule/mdTools/createAutoCorrelationFunctions.H
new file mode 100644
index 0000000000..b257b7eece
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/createAutoCorrelationFunctions.H
@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+Info << nl << "Creating autocorrelation functions." << endl;
+
+IOdictionary mdTransportProperitesDict
+(
+    IOobject
+    (
+        "mdTransportProperitesDict",
+        mesh.time().system(),
+        mesh,
+        IOobject::MUST_READ_IF_MODIFIED,
+        IOobject::NO_WRITE,
+        false
+    )
+);
+
+const dictionary& autocorrelationFunctionDict
+(
+    mdTransportProperitesDict.subDict("autocorrelationFunctions")
+);
+
+//- Velocity autocorrelation function
+
+Info << tab << "velocty" << endl;
+
+const dictionary& velocityACFDict
+(
+    autocorrelationFunctionDict.subDict("velocity")
+);
+
+correlationFunction<vector> vacf
+(
+    mesh,
+    velocityACFDict,
+    molecules.size()
+);
+
+bool writeVacf(Switch(velocityACFDict.lookup("writeFile")));
+
+//- Pressure autocorrelation function
+
+Info << tab << "pressure" << endl;
+
+const dictionary& pressureACFDict
+(
+    autocorrelationFunctionDict.subDict("pressure")
+);
+
+correlationFunction<vector> pacf
+(
+    mesh,
+    pressureACFDict,
+    1
+);
+
+bool writePacf(Switch(pressureACFDict.lookup("writeFile")));
+
+//- Heat flux autocorrelation function
+
+Info << tab << "heat flux" << endl;
+
+const dictionary& heatFluxACFDict
+(
+    autocorrelationFunctionDict.subDict("heatFlux")
+);
+
+correlationFunction<vector> hfacf
+(
+    mesh,
+    heatFluxACFDict,
+    1
+);
+
+bool writeHFacf(Switch(heatFluxACFDict.lookup("writeFile")));
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/createMDFields.H b/src/lagrangian/molecularDynamics/molecule/mdTools/createMDFields.H
new file mode 100644
index 0000000000..6d6617d652
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/createMDFields.H
@@ -0,0 +1,324 @@
+// Fields for data gathering
+
+List< scalarField > allSpeciesN_RU
+(
+    molecules.potential().nIds(),
+    scalarField (mesh.nCells(), 0.0)
+);
+
+List< scalarField > allSpeciesM_RU
+(
+    molecules.potential().nIds(),
+    scalarField (mesh.nCells(), 0.0)
+);
+
+List< vectorField > allSpeciesVelocitySum_RU
+(
+    molecules.potential().nIds(),
+    vectorField (mesh.nCells(), vector::zero)
+);
+
+List< scalarField > allSpeciesVelocityMagSquaredSum_RU
+(
+    molecules.potential().nIds(),
+    scalarField (mesh.nCells(), 0.0)
+);
+
+// Geometric Fields for IO
+
+Info << nl << "Creating fields." << endl;
+
+/*---------------------------------------------------------------------------*\
+    Number density
+\*---------------------------------------------------------------------------*/
+
+PtrList<volScalarField> allSpeciesRhoN
+(
+    molecules.potential().nIds()
+);
+
+forAll(allSpeciesRhoN, rN)
+{
+    Info<< "    Creating number density field for "
+        << molecules.potential().idList()[rN] << endl;
+
+    allSpeciesRhoN.set
+    (
+        rN,
+        new volScalarField
+        (
+            IOobject
+            (
+                "rhoN_" + molecules.potential().idList()[rN],
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimless/dimVolume,
+            "zeroGradient"
+        )
+    );
+    allSpeciesRhoN[rN].internalField() = scalarField (mesh.nCells(), 0.0);
+    allSpeciesRhoN[rN].correctBoundaryConditions();
+}
+
+Info << "    Creating total number density field" << endl;
+
+volScalarField totalRhoN
+(
+    IOobject
+    (
+        "rhoN_total",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimless/dimVolume,
+    "zeroGradient"
+);
+totalRhoN.internalField() = scalarField (mesh.nCells(), 0.0);
+totalRhoN.correctBoundaryConditions();
+
+/*---------------------------------------------------------------------------*\
+    Mass density
+\*---------------------------------------------------------------------------*/
+
+PtrList<volScalarField> allSpeciesRhoM
+(
+    molecules.potential().nIds()
+);
+
+forAll(allSpeciesRhoM, rM)
+{
+    Info<< "    Creating mass density field for "
+        << molecules.potential().idList()[rM] << endl;
+
+    allSpeciesRhoM.set
+    (
+        rM,
+        new volScalarField
+        (
+            IOobject
+            (
+                "rhoM_" + molecules.potential().idList()[rM],
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimMass/dimVolume,
+            "zeroGradient"
+        )
+    );
+    allSpeciesRhoM[rM].internalField() = scalarField (mesh.nCells(), 0.0);
+    allSpeciesRhoM[rM].correctBoundaryConditions();
+}
+
+Info << "    Creating total mass density field" << endl;
+
+volScalarField totalRhoM
+(
+    IOobject
+    (
+        "rhoM_total",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimMass/dimVolume,
+    "zeroGradient"
+);
+totalRhoM.internalField() = scalarField (mesh.nCells(), 0.0);
+totalRhoM.correctBoundaryConditions();
+
+/*---------------------------------------------------------------------------*\
+    Bulk velocity
+\*---------------------------------------------------------------------------*/
+
+PtrList<volVectorField> allSpeciesVelocity
+(
+    molecules.potential().nIds()
+);
+
+forAll(allSpeciesVelocity, v)
+{
+    Info<< "    Creating velocity field for "
+        << molecules.potential().idList()[v] << endl;
+
+    allSpeciesVelocity.set
+    (
+        v,
+        new volVectorField
+        (
+            IOobject
+            (
+                "velocity_" + molecules.potential().idList()[v],
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimVelocity,
+            "zeroGradient"
+        )
+    );
+    allSpeciesVelocity[v].internalField() =
+        vectorField (mesh.nCells(), vector::zero);
+    allSpeciesVelocity[v].correctBoundaryConditions();
+}
+
+Info << "    Creating total velocity field" << endl;
+
+// volVectorField totalVelocity
+// (
+//     IOobject
+//     (
+//         "velocity_total",
+//         runTime.timeName(),
+//         mesh,
+//         IOobject::NO_READ,
+//         IOobject::AUTO_WRITE
+//     ),
+//     mesh,
+//     dimVelocity,
+//     "zeroGradient"
+// );
+// totalVelocity.internalField() = vectorField (mesh.nCells(), vector::zero);
+// totalVelocity.correctBoundaryConditions();
+
+
+volVectorField totalVelocity
+
+(
+    IOobject
+    (
+
+    "velocity_total",
+    runTime.timeName(),
+    mesh,
+    IOobject::NO_READ,
+    IOobject::AUTO_WRITE
+
+    ),
+    mesh,
+    dimensionedVector("zero", dimVelocity, vector::zero)
+);
+
+/*---------------------------------------------------------------------------*\
+    Kinetic temperature
+\*---------------------------------------------------------------------------*/
+
+PtrList<volScalarField> allSpeciesTemperature
+(
+    molecules.potential().nIds()
+);
+
+forAll(allSpeciesTemperature, t)
+{
+    Info<< "    Creating temperature field for "
+        << molecules.potential().idList()[t] << endl;
+
+    allSpeciesTemperature.set
+    (
+        t,
+        new volScalarField
+        (
+            IOobject
+            (
+                "temperature_" + molecules.potential().idList()[t],
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimTemperature,
+            "zeroGradient"
+        )
+    );
+    allSpeciesTemperature[t].internalField() = scalarField (mesh.nCells(), 0.0);
+    allSpeciesTemperature[t].correctBoundaryConditions();
+}
+
+Info << "    Creating total temperature field" << endl;
+
+volScalarField totalTemperature
+(
+    IOobject
+    (
+        "temperature_total",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimTemperature,
+    "zeroGradient"
+);
+totalTemperature.internalField() = scalarField (mesh.nCells(), 0.0);
+totalTemperature.correctBoundaryConditions();
+
+/*---------------------------------------------------------------------------*\
+    Mean kinetic energy
+\*---------------------------------------------------------------------------*/
+
+
+PtrList<volScalarField> allSpeciesMeanKE
+(
+    molecules.potential().nIds()
+);
+
+forAll(allSpeciesMeanKE, mKE)
+{
+    Info<< "    Creating mean kinetic energy field for "
+        << molecules.potential().idList()[mKE] << endl;
+
+    allSpeciesMeanKE.set
+    (
+        mKE,
+        new volScalarField
+        (
+            IOobject
+            (
+                "meanKE_" + molecules.potential().idList()[mKE],
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimensionSet(1, 2, -2, 0, 0, 0, 0),
+            "zeroGradient"
+        )
+    );
+    allSpeciesMeanKE[mKE].internalField() = scalarField (mesh.nCells(), 0.0);
+    allSpeciesMeanKE[mKE].correctBoundaryConditions();
+}
+
+Info << "    Creating total mean kinetic energy field" << endl;
+
+volScalarField totalMeanKE
+(
+    IOobject
+    (
+        "meanKE_total",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionSet(1, 2, -2, 0, 0, 0, 0),
+    "zeroGradient"
+);
+totalMeanKE.internalField() = scalarField (mesh.nCells(), 0.0);
+totalMeanKE.correctBoundaryConditions();
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/createRefUnits.H b/src/lagrangian/molecularDynamics/molecule/mdTools/createRefUnits.H
new file mode 100644
index 0000000000..c1e558b221
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/createRefUnits.H
@@ -0,0 +1,22 @@
+reducedUnits refUnits;
+
+IOobject reducedUnitsDictIOobject
+(
+    "reducedUnitsDict",
+    runTime.system(),
+    mesh,
+    IOobject::READ_IF_PRESENT,
+    IOobject::NO_WRITE
+);
+
+if (reducedUnitsDictIOobject.headerOk())
+{
+    Info<< nl
+        << "Reading reference quantities from reducedUnitsDict file." << endl;
+
+    IOdictionary reducedUnitsDict(reducedUnitsDictIOobject);
+
+    refUnits.setRefValues(reducedUnitsDict);
+}
+
+Info << refUnits << endl;
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/md.H b/src/lagrangian/molecularDynamics/molecule/mdTools/md.H
new file mode 100644
index 0000000000..fc830da843
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/md.H
@@ -0,0 +1,8 @@
+#ifndef md_H
+#define md_H
+    #include "potential.H"
+    #include "moleculeCloud.H"
+    #include "correlationFunction.H"
+    #include "distribution.H"
+    #include "reducedUnits.H"
+#endif
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/meanMomentumEnergyAndNMols.H b/src/lagrangian/molecularDynamics/molecule/mdTools/meanMomentumEnergyAndNMols.H
new file mode 100644
index 0000000000..77ecb7a657
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/meanMomentumEnergyAndNMols.H
@@ -0,0 +1,180 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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/>.
+
+Global
+    meanMomentumEnergyAndNMols.H
+
+Description
+    Calculates and prints the mean momentum and energy in the system
+    and the number of molecules.
+
+\*---------------------------------------------------------------------------*/
+
+
+vector singleStepTotalLinearMomentum(vector::zero);
+
+vector singleStepTotalAngularMomentum(vector::zero);
+
+scalar singleStepMaxVelocityMag = 0.0;
+
+scalar singleStepTotalMass = 0.0;
+
+scalar singleStepTotalLinearKE = 0.0;
+
+scalar singleStepTotalAngularKE = 0.0;
+
+scalar singleStepTotalPE = 0.0;
+
+scalar singleStepTotalrDotf = 0.0;
+
+//vector singleStepCentreOfMass(vector::zero);
+
+label singleStepNMols = molecules.size();
+
+label singleStepDOFs = 0;
+
+{
+    forAllConstIter(IDLList<molecule>, molecules, mol)
+    {
+        const label molId = mol().id();
+
+        scalar molMass(molecules.constProps(molId).mass());
+
+        singleStepTotalMass += molMass;
+
+        //singleStepCentreOfMass += mol().position()*molMass;
+    }
+
+    // if (singleStepNMols)
+    // {
+    //     singleStepCentreOfMass /= singleStepTotalMass;
+    // }
+
+    forAllConstIter(IDLList<molecule>, molecules, mol)
+    {
+        const label molId = mol().id();
+
+        const molecule::constantProperties cP(molecules.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();
+
+        singleStepTotalLinearMomentum += molV * molMass;
+
+        singleStepTotalAngularMomentum += molPiGlobal;
+        //+((mol().position() - singleStepCentreOfMass) ^ (molV * molMass));
+
+        if (mag(molV) > singleStepMaxVelocityMag)
+        {
+            singleStepMaxVelocityMag = mag(molV);
+        }
+
+        singleStepTotalLinearKE += 0.5*molMass*magSqr(molV);
+
+        singleStepTotalAngularKE += 0.5*(molOmega & molMoI & molOmega);
+
+        singleStepTotalPE += mol().potentialEnergy();
+
+        singleStepTotalrDotf += tr(mol().rf());
+
+        singleStepDOFs += cP.degreesOfFreedom();
+    }
+}
+
+if (Pstream::parRun())
+{
+    reduce(singleStepTotalLinearMomentum, sumOp<vector>());
+
+    reduce(singleStepTotalAngularMomentum, sumOp<vector>());
+
+    reduce(singleStepMaxVelocityMag, maxOp<scalar>());
+
+    reduce(singleStepTotalMass, sumOp<scalar>());
+
+    reduce(singleStepTotalLinearKE, sumOp<scalar>());
+
+    reduce(singleStepTotalAngularKE, sumOp<scalar>());
+
+    reduce(singleStepTotalPE, sumOp<scalar>());
+
+    reduce(singleStepTotalrDotf, sumOp<scalar>());
+
+    reduce(singleStepNMols, sumOp<label>());
+
+    reduce(singleStepDOFs, sumOp<label>());
+}
+
+if (singleStepNMols)
+{
+    Info<< "Number of molecules in system = "
+        << singleStepNMols << nl
+        << "Overall number density = "
+        << singleStepNMols/meshVolume << nl
+        << "Overall mass density = "
+        << singleStepTotalMass/meshVolume << nl
+        << "Average linear momentum per molecule = "
+        << singleStepTotalLinearMomentum/singleStepNMols << ' '
+        << mag(singleStepTotalLinearMomentum)/singleStepNMols << nl
+        << "Average angular momentum per molecule = "
+        << singleStepTotalAngularMomentum << ' '
+        << mag(singleStepTotalAngularMomentum)/singleStepNMols << nl
+        << "Maximum |velocity| = "
+        << singleStepMaxVelocityMag << nl
+        << "Average linear KE per molecule = "
+        << singleStepTotalLinearKE/singleStepNMols << nl
+        << "Average angular KE per molecule = "
+        << singleStepTotalAngularKE/singleStepNMols << nl
+        << "Average PE per molecule = "
+        << singleStepTotalPE/singleStepNMols << nl
+        << "Average TE per molecule = "
+        <<
+        (
+            singleStepTotalLinearKE
+          + singleStepTotalAngularKE
+          + singleStepTotalPE
+        )
+        /singleStepNMols
+        << endl;
+
+        // Info<< singleStepNMols << " "
+        //     << singleStepTotalMomentum/singleStepTotalMass << " "
+        //     << singleStepMaxVelocityMag << " "
+        //     << singleStepTotalKE/singleStepNMols << " "
+        //     << singleStepTotalPE/singleStepNMols << " "
+        //     << (singleStepTotalKE + singleStepTotalPE)
+        //        /singleStepNMols << endl;
+}
+else
+{
+    Info<< "No molecules in system" << endl;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/resetMDFields.H b/src/lagrangian/molecularDynamics/molecule/mdTools/resetMDFields.H
new file mode 100644
index 0000000000..af0197b5fe
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/resetMDFields.H
@@ -0,0 +1,26 @@
+if (runTime.outputTime())
+{
+    allSpeciesN_RU = List< scalarField >
+    (
+        molecules.potential().nIds(),
+        scalarField(mesh.nCells(), 0.0)
+    );
+
+    allSpeciesM_RU = List< scalarField >
+    (
+        molecules.potential().nIds(),
+        scalarField(mesh.nCells(), 0.0)
+    );
+
+    allSpeciesVelocitySum_RU = List< vectorField >
+    (
+        molecules.potential().nIds(),
+        vectorField(mesh.nCells(), vector::zero)
+    );
+
+    allSpeciesVelocityMagSquaredSum_RU = List< scalarField >
+    (
+        molecules.potential().nIds(),
+        scalarField(mesh.nCells(), 0.0)
+    );
+}
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H
new file mode 100644
index 0000000000..247402bb9f
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H
@@ -0,0 +1,111 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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/>.
+
+Global
+    temperatureAndPressure.H
+
+Description
+    Accumulates values for temperature and pressure measurement, and
+    calculates and outputs the average values at output times.
+    Requires temperatureAndPressureVariables.H to be declared before the
+    timeloop.
+
+\*---------------------------------------------------------------------------*/
+
+accumulatedTotalLinearMomentum += singleStepTotalLinearMomentum;
+
+accumulatedTotalMass += singleStepTotalMass;
+
+accumulatedTotalLinearKE += singleStepTotalLinearKE;
+
+accumulatedTotalAngularKE += singleStepTotalAngularKE;
+
+accumulatedTotalPE += singleStepTotalPE;
+
+accumulatedTotalrDotfSum += singleStepTotalrDotf;
+
+accumulatedNMols += singleStepNMols;
+
+accumulatedDOFs += singleStepDOFs;
+
+if (runTime.outputTime())
+{
+    if (accumulatedNMols)
+    {
+        Info<< "calculating averages" << endl;
+
+        averageTemperature =
+        (
+            2.0/(physicoChemical::k.value()*accumulatedDOFs)
+            *
+            (
+                accumulatedTotalLinearKE + accumulatedTotalAngularKE
+                -
+                0.5*magSqr(accumulatedTotalLinearMomentum)/accumulatedTotalMass
+            )
+        );
+
+        averagePressure =
+        (
+            (
+                (accumulatedNMols/nAveragingSteps)
+               *physicoChemical::k.value()*averageTemperature
+              + accumulatedTotalrDotfSum/(6.0*nAveragingSteps)
+            )
+            /
+            meshVolume
+        );
+
+        Info<< "----------------------------------------" << nl
+            << "Averaged properties" << nl
+            << "Average |velocity| = "
+            << mag(accumulatedTotalLinearMomentum)/accumulatedTotalMass << nl
+            << "Average temperature = " << averageTemperature << nl
+            << "Average pressure = " << averagePressure << nl
+            << "----------------------------------------" << endl;
+    }
+    else
+    {
+        Info<< "Not averaging temperature and pressure: "
+            << "no molecules in system" << endl;
+    }
+
+    accumulatedTotalLinearMomentum = vector::zero;
+
+    accumulatedTotalMass = 0.0;
+
+    accumulatedTotalLinearKE = 0.0;
+
+    accumulatedTotalAngularKE = 0.0;
+
+    accumulatedTotalPE = 0.0;
+
+    accumulatedTotalrDotfSum = 0.0;
+
+    accumulatedNMols = 0;
+
+    accumulatedDOFs = 0;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressureVariables.H b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressureVariables.H
new file mode 100644
index 0000000000..3c671ef11c
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressureVariables.H
@@ -0,0 +1,62 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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/>.
+
+Global
+    temperatureAndPressureVariables.H
+
+Description
+    Provides accumulation variables for temperatureAndPressure.H
+
+\*---------------------------------------------------------------------------*/
+
+vector accumulatedTotalLinearMomentum(vector::zero);
+
+scalar accumulatedTotalMass = 0.0;
+
+scalar accumulatedTotalAngularKE = 0.0;
+
+scalar accumulatedTotalLinearKE = 0.0;
+
+scalar accumulatedTotalPE = 0.0;
+
+scalar accumulatedTotalrDotfSum = 0.0;
+
+label accumulatedNMols = 0;
+
+label accumulatedDOFs = 0;
+
+scalar averageTemperature = 0.0;
+
+scalar averagePressure = 0.0;
+
+const scalarField& cellVols = mesh.cellVolumes();
+
+scalar meshVolume = sum(cellVols);
+
+if (Pstream::parRun())
+{
+    reduce(meshVolume, sumOp<scalar>());
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureEquilibration.H b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureEquilibration.H
new file mode 100644
index 0000000000..a4a9133749
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureEquilibration.H
@@ -0,0 +1,42 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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/>.
+
+Global
+    temperatureEquilibration.H
+
+Description
+    Applies temperature control to molecules
+
+\*---------------------------------------------------------------------------*/
+
+if (runTime.outputTime())
+{
+    molecules.applyConstraintsAndThermostats
+    (
+        targetTemperature,
+        averageTemperature
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
new file mode 100644
index 0000000000..1163d3f184
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
@@ -0,0 +1,296 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "moleculeCloud.H"
+#include "molecule.H"
+#include "Random.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::tensor Foam::molecule::rotationTensorX(scalar phi) const
+{
+    return tensor
+    (
+        1, 0, 0,
+        0, Foam::cos(phi), -Foam::sin(phi),
+        0, Foam::sin(phi), Foam::cos(phi)
+    );
+}
+
+
+Foam::tensor Foam::molecule::rotationTensorY(scalar phi) const
+{
+    return tensor
+    (
+        Foam::cos(phi), 0, Foam::sin(phi),
+        0, 1, 0,
+        -Foam::sin(phi), 0, Foam::cos(phi)
+    );
+}
+
+
+Foam::tensor Foam::molecule::rotationTensorZ(scalar phi) const
+{
+    return tensor
+    (
+        Foam::cos(phi), -Foam::sin(phi), 0,
+        Foam::sin(phi), Foam::cos(phi), 0,
+        0, 0, 1
+    );
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::molecule::move(molecule::trackingData& td, const scalar trackTime)
+{
+    td.switchProcessor = false;
+    td.keepParticle = true;
+
+    const constantProperties& constProps(td.cloud().constProps(id_));
+
+    if (td.part() == 0)
+    {
+        // First leapfrog velocity adjust part, required before tracking+force
+        // part
+
+        v_ += 0.5*trackTime*a_;
+
+        pi_ += 0.5*trackTime*tau_;
+    }
+    else if (td.part() == 1)
+    {
+        // 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;
+        }
+    }
+    else if (td.part() == 2)
+    {
+        // Leapfrog orientation adjustment, carried out before force calculation
+        // but after tracking stage, i.e. rotation carried once linear motion
+        // complete.
+
+        if (!constProps.pointMolecule())
+        {
+            const diagTensor& momentOfInertia(constProps.momentOfInertia());
+
+            tensor R;
+
+            if (!constProps.linearMolecule())
+            {
+                R = rotationTensorX(0.5*trackTime*pi_.x()/momentOfInertia.xx());
+                pi_ = pi_ & R;
+                Q_ = Q_ & R;
+            }
+
+            R = rotationTensorY(0.5*trackTime*pi_.y()/momentOfInertia.yy());
+            pi_ = pi_ & R;
+            Q_ = Q_ & R;
+
+            R = rotationTensorZ(trackTime*pi_.z()/momentOfInertia.zz());
+            pi_ = pi_ & R;
+            Q_ = Q_ & R;
+
+            R = rotationTensorY(0.5*trackTime*pi_.y()/momentOfInertia.yy());
+            pi_ = pi_ & R;
+            Q_ = Q_ & R;
+
+            if (!constProps.linearMolecule())
+            {
+                R = rotationTensorX(0.5*trackTime*pi_.x()/momentOfInertia.xx());
+                pi_ = pi_ & R;
+                Q_ = Q_ & R;
+            }
+        }
+
+        setSitePositions(constProps);
+    }
+    else if (td.part() == 3)
+    {
+        // Second leapfrog velocity adjust part, required after tracking+force
+        // part
+
+        scalar m = constProps.mass();
+
+        a_ = vector::zero;
+
+        tau_ = vector::zero;
+
+        forAll(siteForces_, s)
+        {
+            const vector& f = siteForces_[s];
+
+            a_ += f/m;
+
+            tau_ += (constProps.siteReferencePositions()[s] ^ (Q_.T() & f));
+        }
+
+        v_ += 0.5*trackTime*a_;
+
+        pi_ += 0.5*trackTime*tau_;
+
+        if (constProps.pointMolecule())
+        {
+            tau_ = vector::zero;
+
+            pi_ = vector::zero;
+        }
+
+        if (constProps.linearMolecule())
+        {
+            tau_.x() = 0.0;
+
+            pi_.x() = 0.0;
+        }
+    }
+    else
+    {
+        FatalErrorIn("molecule::move(trackingData&, const scalar)") << nl
+            << td.part() << " is an invalid part of the integration method."
+            << abort(FatalError);
+    }
+
+    return td.keepParticle;
+}
+
+
+void Foam::molecule::transformProperties(const tensor& T)
+{
+    particle::transformProperties(T);
+
+    Q_ = T & Q_;
+
+    v_ = transform(T, v_);
+
+    a_ = transform(T, a_);
+
+    pi_ = Q_.T() & transform(T, Q_ & pi_);
+
+    tau_ = Q_.T() & transform(T, Q_ & tau_);
+
+    rf_ = transform(T, rf_);
+
+    sitePositions_ = position_ + (T & (sitePositions_ - position_));
+
+    siteForces_ = T & siteForces_;
+}
+
+
+void Foam::molecule::transformProperties(const vector& separation)
+{
+    particle::transformProperties(separation);
+
+    if (special_ == SPECIAL_TETHERED)
+    {
+        specialPosition_ += separation;
+    }
+
+    sitePositions_ = sitePositions_ + separation;
+}
+
+
+void Foam::molecule::setSitePositions(const constantProperties& constProps)
+{
+    sitePositions_ = position_ + (Q_ & constProps.siteReferencePositions());
+}
+
+
+void Foam::molecule::setSiteSizes(label size)
+{
+    sitePositions_.setSize(size);
+
+    siteForces_.setSize(size);
+}
+
+
+bool Foam::molecule::hitPatch
+(
+    const polyPatch&,
+    trackingData&,
+    const label,
+    const scalar,
+    const tetIndices&
+)
+{
+    return false;
+}
+
+
+void Foam::molecule::hitProcessorPatch
+(
+    const processorPolyPatch&,
+    trackingData& td
+)
+{
+    td.switchProcessor = true;
+}
+
+
+void Foam::molecule::hitWallPatch
+(
+    const wallPolyPatch& wpp,
+    trackingData& td,
+    const tetIndices& tetIs
+)
+{
+    // Use of the normal from tetIs is not required as
+    // hasWallImpactDistance for a moleculeCloud is false.
+    vector nw = normal();
+    nw /= mag(nw);
+
+    scalar vn = v_ & nw;
+
+    // Specular reflection
+    if (vn > 0)
+    {
+        v_ -= 2*vn*nw;
+    }
+}
+
+
+void Foam::molecule::hitPatch
+(
+    const polyPatch&,
+    trackingData& td
+)
+{
+    td.keepParticle = false;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
new file mode 100644
index 0000000000..54bdb51b9f
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
@@ -0,0 +1,392 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::molecule
+
+Description
+    Foam::molecule
+
+SourceFiles
+    moleculeI.H
+    molecule.C
+    moleculeIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef molecule_H
+#define molecule_H
+
+#include "particle.H"
+#include "IOstream.H"
+#include "autoPtr.H"
+#include "diagTensor.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Class forward declarations
+class moleculeCloud;
+
+/*---------------------------------------------------------------------------*\
+                           Class molecule Declaration
+\*---------------------------------------------------------------------------*/
+
+class molecule
+:
+    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 molecule constant properties
+    class constantProperties
+    {
+        // Private data
+
+            Field<vector> siteReferencePositions_;
+
+            List<scalar> siteMasses_;
+
+            List<scalar> siteCharges_;
+
+            List<label> siteIds_;
+
+            List<bool> pairPotentialSites_;
+
+            List<bool> electrostaticSites_;
+
+            diagTensor momentOfInertia_;
+
+            scalar mass_;
+
+
+        // Private Member Functions
+
+            void checkSiteListSizes() const;
+
+            void setInteracionSiteBools
+            (
+                const List<word>& siteIds,
+                const List<word>& pairPotSiteIds
+            );
+
+            bool linearMoleculeTest() const;
+
+
+    public:
+
+        inline constantProperties();
+
+        //- Construct from dictionary
+        inline constantProperties(const dictionary& dict);
+
+        // Member functions
+
+        inline const Field<vector>& siteReferencePositions() const;
+
+        inline const List<scalar>& siteMasses() const;
+
+        inline const List<scalar>& siteCharges() const;
+
+        inline const List<label>& siteIds() const;
+
+        inline List<label>& siteIds();
+
+        inline const List<bool>& pairPotentialSites() const;
+
+        inline bool pairPotentialSite(label sId) const;
+
+        inline const List<bool>& electrostaticSites() const;
+
+        inline bool electrostaticSite(label sId) 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;
+    };
+
+
+    //- Class used to pass tracking data to the trackToFace function
+    class trackingData
+    :
+        public particle::TrackingData<moleculeCloud>
+    {
+        // label specifying which part of the integration algorithm is taking
+        label part_;
+
+
+    public:
+
+        // Constructors
+
+            trackingData(moleculeCloud& cloud, label part)
+            :
+                particle::TrackingData<moleculeCloud>(cloud),
+                part_(part)
+            {}
+
+        // Member functions
+
+            inline label part() const
+            {
+                return part_;
+            }
+    };
+
+
+private:
+
+    // Private data
+
+        tensor Q_;
+
+        vector v_;
+
+        vector a_;
+
+        vector pi_;
+
+        vector tau_;
+
+        vector specialPosition_;
+
+        scalar potentialEnergy_;
+
+        // - r_ij f_ij, stress dyad
+        tensor rf_;
+
+        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:
+
+    friend class Cloud<molecule>;
+
+    // Constructors
+
+        //- Construct from components
+        inline molecule
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label cellI,
+            const label tetFaceI,
+            const label tetPtI,
+            const tensor& Q,
+            const vector& v,
+            const vector& a,
+            const vector& pi,
+            const vector& tau,
+            const vector& specialPosition,
+            const constantProperties& constProps,
+            const label special,
+            const label id
+        );
+
+        //- Construct from Istream
+        molecule
+        (
+            const polyMesh& mesh,
+            Istream& is,
+            bool readFields = true
+        );
+
+        //- Construct and return a clone
+        autoPtr<particle> clone() const
+        {
+            return autoPtr<particle>(new molecule(*this));
+        }
+
+        //- Factory class to read-construct particles used for
+        //  parallel transfer
+        class iNew
+        {
+            const polyMesh& mesh_;
+
+        public:
+
+            iNew(const polyMesh& mesh)
+            :
+                mesh_(mesh)
+            {}
+
+            autoPtr<molecule> operator()(Istream& is) const
+            {
+                return autoPtr<molecule>(new molecule(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 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;
+
+
+    // 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
+
+        static void readFields(Cloud<molecule>& mC);
+
+        static void writeFields(const Cloud<molecule>& mC);
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<<(Ostream&, const molecule&);
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "moleculeI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
new file mode 100644
index 0000000000..451ed5750c
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
@@ -0,0 +1,606 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::molecule::constantProperties::constantProperties()
+:
+    siteReferencePositions_(Field<vector>(0)),
+    siteMasses_(List<scalar>(0)),
+    siteCharges_(List<scalar>(0)),
+    siteIds_(List<label>(0)),
+    pairPotentialSites_(List<bool>(false)),
+    electrostaticSites_(List<bool>(false)),
+    momentOfInertia_(diagTensor(0, 0, 0)),
+    mass_(0)
+{}
+
+
+inline Foam::molecule::constantProperties::constantProperties
+(
+    const dictionary& dict
+)
+:
+    siteReferencePositions_(dict.lookup("siteReferencePositions")),
+    siteMasses_(dict.lookup("siteMasses")),
+    siteCharges_(dict.lookup("siteCharges")),
+    siteIds_(List<word>(dict.lookup("siteIds")).size(), -1),
+    pairPotentialSites_(),
+    electrostaticSites_(),
+    momentOfInertia_(),
+    mass_()
+{
+    checkSiteListSizes();
+
+    setInteracionSiteBools
+    (
+        List<word>(dict.lookup("siteIds")),
+        List<word>(dict.lookup("pairPotentialSiteIds"))
+    );
+
+    mass_ = sum(siteMasses_);
+
+    vector centreOfMass(vector::zero);
+
+    // Calculate the centre of mass of the body and subtract it from each
+    // position
+
+    forAll(siteReferencePositions_, i)
+    {
+        centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
+    }
+
+    centreOfMass /= mass_;
+
+    siteReferencePositions_ -= centreOfMass;
+
+    if (siteIds_.size() == 1)
+    {
+        // Single site molecule - no rotational motion.
+
+        siteReferencePositions_[0] = vector::zero;
+
+        momentOfInertia_ = diagTensor(-1, -1, -1);
+    }
+    else if (linearMoleculeTest())
+    {
+        // Linear molecule.
+
+        Info<< nl << "Linear molecule." << endl;
+
+        vector dir = siteReferencePositions_[1] - siteReferencePositions_[0];
+
+        dir /= mag(dir);
+
+        tensor Q = rotationTensor(dir, vector(1,0,0));
+
+        siteReferencePositions_ = (Q & siteReferencePositions_);
+
+        // The rotation was around the centre of mass but remove any
+        // components that have crept in due to floating point errors
+
+        centreOfMass = vector::zero;
+
+        forAll(siteReferencePositions_, i)
+        {
+            centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
+        }
+
+        centreOfMass /= mass_;
+
+        siteReferencePositions_ -= centreOfMass;
+
+        diagTensor momOfInertia = diagTensor::zero;
+
+        forAll(siteReferencePositions_, i)
+        {
+            const vector& p(siteReferencePositions_[i]);
+
+            momOfInertia +=
+                siteMasses_[i]*diagTensor(0, p.x()*p.x(), p.x()*p.x());
+        }
+
+        momentOfInertia_ = diagTensor
+         (
+             -1,
+             momOfInertia.yy(),
+             momOfInertia.zz()
+         );
+    }
+    else
+    {
+        // Fully 6DOF molecule
+
+        // Calculate the inertia tensor in the current orientation
+
+        tensor momOfInertia(tensor::zero);
+
+        forAll(siteReferencePositions_, i)
+        {
+            const vector& p(siteReferencePositions_[i]);
+
+            momOfInertia += siteMasses_[i]*tensor
+            (
+                p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(),
+                -p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(),
+                -p.z()*p.x(), -p.z()*p.y(), p.x()*p.x() + p.y()*p.y()
+            );
+        }
+
+        if (eigenValues(momOfInertia).x() < VSMALL)
+        {
+            FatalErrorIn("molecule::constantProperties::constantProperties")
+                << "An eigenvalue of the inertia tensor is zero, the molecule "
+                << dict.name() << " is not a valid 6DOF shape."
+                << nl << abort(FatalError);
+        }
+
+        // Normalise the inertia tensor magnitude to avoid SMALL numbers in the
+        // components causing problems
+
+        momOfInertia /= eigenValues(momOfInertia).x();
+
+        tensor e = eigenVectors(momOfInertia);
+
+        // Calculate the transformation between the principle axes to the
+        // global axes
+
+        tensor Q =
+            vector(1,0,0)*e.x() + vector(0,1,0)*e.y() + vector(0,0,1)*e.z();
+
+        // Transform the site positions
+
+        siteReferencePositions_ = (Q & siteReferencePositions_);
+
+        // Recalculating the moment of inertia with the new site positions
+
+        // The rotation was around the centre of mass but remove any
+        // components that have crept in due to floating point errors
+
+        centreOfMass = vector::zero;
+
+        forAll(siteReferencePositions_, i)
+        {
+            centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
+        }
+
+        centreOfMass /= mass_;
+
+        siteReferencePositions_ -= centreOfMass;
+
+        // Calculate the moment of inertia in the principle component
+        // reference frame
+
+        momOfInertia = tensor::zero;
+
+        forAll(siteReferencePositions_, i)
+        {
+            const vector& p(siteReferencePositions_[i]);
+
+            momOfInertia += siteMasses_[i]*tensor
+            (
+                p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(),
+                -p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(),
+                -p.z()*p.x(), -p.z()*p.y(), p.x()*p.x() + p.y()*p.y()
+            );
+        }
+
+        momentOfInertia_ = diagTensor
+        (
+            momOfInertia.xx(),
+            momOfInertia.yy(),
+            momOfInertia.zz()
+        );
+    }
+}
+
+
+inline Foam::molecule::molecule
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label cellI,
+    const label tetFaceI,
+    const label tetPtI,
+    const tensor& Q,
+    const vector& v,
+    const vector& a,
+    const vector& pi,
+    const vector& tau,
+    const vector& specialPosition,
+    const constantProperties& constProps,
+    const label special,
+    const label id
+
+)
+:
+    particle(mesh, position, cellI, tetFaceI, tetPtI),
+    Q_(Q),
+    v_(v),
+    a_(a),
+    pi_(pi),
+    tau_(tau),
+    specialPosition_(specialPosition),
+    potentialEnergy_(0.0),
+    rf_(tensor::zero),
+    special_(special),
+    id_(id),
+    siteForces_(constProps.nSites(), vector::zero),
+    sitePositions_(constProps.nSites())
+{
+    setSitePositions(constProps);
+}
+
+
+// * * * * * * * constantProperties Private Member Functions * * * * * * * * //
+
+inline void Foam::molecule::constantProperties::checkSiteListSizes() const
+{
+    if
+    (
+        siteIds_.size() != siteReferencePositions_.size()
+     || siteIds_.size() != siteCharges_.size()
+    )
+    {
+        FatalErrorIn("molecule::constantProperties::checkSiteListSizes")
+            << "Sizes of site id, charge and "
+            << "referencePositions are not the same. "
+            << nl << abort(FatalError);
+    }
+}
+
+
+inline void Foam::molecule::constantProperties::setInteracionSiteBools
+(
+    const List<word>& siteIds,
+    const List<word>& pairPotSiteIds
+)
+{
+    pairPotentialSites_.setSize(siteIds_.size());
+
+    electrostaticSites_.setSize(siteIds_.size());
+
+    forAll(siteIds_, i)
+    {
+        const word& id(siteIds[i]);
+
+        pairPotentialSites_[i] = (findIndex(pairPotSiteIds, id) > -1);
+
+        electrostaticSites_[i] = (mag(siteCharges_[i]) > VSMALL);
+    }
+}
+
+
+inline bool Foam::molecule::constantProperties::linearMoleculeTest() const
+{
+    if (siteIds_.size() == 2)
+    {
+        return true;
+    }
+
+    vector refDir = siteReferencePositions_[1] - siteReferencePositions_[0];
+
+    refDir /= mag(refDir);
+
+    for
+    (
+        label i = 2;
+        i < siteReferencePositions_.size();
+        i++
+    )
+    {
+        vector dir = siteReferencePositions_[i] - siteReferencePositions_[i-1];
+
+        dir /= mag(dir);
+
+        if (mag(refDir & dir) < 1 - SMALL)
+        {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+
+// * * * * * * * constantProperties Member Functions * * * * * * * * * * * * //
+
+inline const Foam::Field<Foam::vector>&
+Foam::molecule::constantProperties::siteReferencePositions() const
+{
+    return siteReferencePositions_;
+}
+
+
+inline const Foam::List<Foam::scalar>&
+Foam::molecule::constantProperties::siteMasses() const
+{
+    return siteMasses_;
+}
+
+
+inline const Foam::List<Foam::scalar>&
+Foam::molecule::constantProperties::siteCharges() const
+{
+    return siteCharges_;
+}
+
+
+inline const Foam::List<Foam::label>&
+Foam::molecule::constantProperties::siteIds() const
+{
+    return siteIds_;
+}
+
+
+inline Foam::List<Foam::label>&
+Foam::molecule::constantProperties::siteIds()
+{
+    return siteIds_;
+}
+
+
+inline const Foam::List<bool>&
+Foam::molecule::constantProperties::pairPotentialSites() const
+{
+    return pairPotentialSites_;
+}
+
+
+inline bool Foam::molecule::constantProperties::pairPotentialSite
+(
+    label sId
+) const
+{
+    label s = findIndex(siteIds_, sId);
+
+    if (s == -1)
+    {
+        FatalErrorIn("moleculeI.H") << nl
+            << sId << " site not found."
+            << nl << abort(FatalError);
+    }
+
+    return pairPotentialSites_[s];
+}
+
+
+inline const Foam::List<bool>&
+Foam::molecule::constantProperties::electrostaticSites() const
+{
+    return electrostaticSites_;
+}
+
+
+inline bool Foam::molecule::constantProperties::electrostaticSite
+(
+    label sId
+) const
+{
+    label s = findIndex(siteIds_, sId);
+
+    if (s == -1)
+    {
+        FatalErrorIn
+        (
+            "molecule::constantProperties::electrostaticSite(label)"
+        )   << sId << " site not found."
+            << nl << abort(FatalError);
+    }
+
+    return electrostaticSites_[s];
+}
+
+
+inline const Foam::diagTensor&
+Foam::molecule::constantProperties::momentOfInertia() const
+{
+    return momentOfInertia_;
+}
+
+
+inline bool Foam::molecule::constantProperties::linearMolecule() const
+{
+    return ((momentOfInertia_.xx() < 0) && (momentOfInertia_.yy() > 0));
+}
+
+
+inline bool Foam::molecule::constantProperties::pointMolecule() const
+{
+    return (momentOfInertia_.zz() < 0);
+}
+
+
+inline Foam::label Foam::molecule::constantProperties::degreesOfFreedom() const
+{
+    if (linearMolecule())
+    {
+        return 5;
+    }
+    else if (pointMolecule())
+    {
+        return 3;
+    }
+    else
+    {
+        return 6;
+    }
+}
+
+
+inline Foam::scalar Foam::molecule::constantProperties::mass() const
+{
+    return mass_;
+}
+
+
+inline Foam::label Foam::molecule::constantProperties::nSites() const
+{
+    return siteIds_.size();
+
+}
+
+
+// * * * * * * * * * * * * molecule Member Functions * * * * * * * * * * * * //
+
+inline const Foam::tensor& Foam::molecule::Q() const
+{
+    return Q_;
+}
+
+
+inline Foam::tensor& Foam::molecule::Q()
+{
+    return Q_;
+}
+
+
+inline const Foam::vector& Foam::molecule::v() const
+{
+    return v_;
+}
+
+
+inline Foam::vector& Foam::molecule::v()
+{
+    return v_;
+}
+
+
+inline const Foam::vector& Foam::molecule::a() const
+{
+    return a_;
+}
+
+
+inline Foam::vector& Foam::molecule::a()
+{
+    return a_;
+}
+
+
+inline const Foam::vector& Foam::molecule::pi() const
+{
+    return pi_;
+}
+
+
+inline Foam::vector& Foam::molecule::pi()
+{
+    return pi_;
+}
+
+
+inline const Foam::vector& Foam::molecule::tau() const
+{
+    return tau_;
+}
+
+
+inline Foam::vector& Foam::molecule::tau()
+{
+    return tau_;
+}
+
+
+inline const Foam::List<Foam::vector>& Foam::molecule::siteForces() const
+{
+    return siteForces_;
+}
+
+
+inline Foam::List<Foam::vector>& Foam::molecule::siteForces()
+{
+    return siteForces_;
+}
+
+
+inline const Foam::List<Foam::vector>& Foam::molecule::sitePositions() const
+{
+    return sitePositions_;
+}
+
+
+inline Foam::List<Foam::vector>& Foam::molecule::sitePositions()
+{
+    return sitePositions_;
+}
+
+
+inline const Foam::vector& Foam::molecule::specialPosition() const
+{
+    return specialPosition_;
+}
+
+
+inline Foam::vector& Foam::molecule::specialPosition()
+{
+    return specialPosition_;
+}
+
+
+inline Foam::scalar Foam::molecule::potentialEnergy() const
+{
+    return potentialEnergy_;
+}
+
+
+inline Foam::scalar& Foam::molecule::potentialEnergy()
+{
+    return potentialEnergy_;
+}
+
+
+inline const Foam::tensor& Foam::molecule::rf() const
+{
+    return rf_;
+}
+
+
+inline Foam::tensor& Foam::molecule::rf()
+{
+    return rf_;
+}
+
+
+inline Foam::label Foam::molecule::special() const
+{
+    return special_;
+}
+
+
+inline bool Foam::molecule::tethered() const
+{
+    return special_ == SPECIAL_TETHERED;
+}
+
+
+inline Foam::label Foam::molecule::id() const
+{
+    return id_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C
new file mode 100644
index 0000000000..77685e6c0d
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C
@@ -0,0 +1,312 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "molecule.H"
+#include "IOstreams.H"
+#include "moleculeCloud.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::molecule::molecule
+(
+    const polyMesh& mesh,
+    Istream& is,
+    bool readFields
+)
+:
+    particle(mesh, is, readFields),
+    Q_(tensor::zero),
+    v_(vector::zero),
+    a_(vector::zero),
+    pi_(vector::zero),
+    tau_(vector::zero),
+    specialPosition_(vector::zero),
+    potentialEnergy_(0.0),
+    rf_(tensor::zero),
+    special_(0),
+    id_(0),
+    siteForces_(0),
+    sitePositions_(0)
+{
+    if (readFields)
+    {
+        if (is.format() == IOstream::ASCII)
+        {
+            is  >> Q_;
+            is  >> v_;
+            is  >> a_;
+            is  >> pi_;
+            is  >> tau_;
+            is  >> siteForces_;
+            potentialEnergy_ = readScalar(is);
+            is  >> rf_;
+            special_ = readLabel(is);
+            id_ = readLabel(is);
+            is  >> sitePositions_;
+            is  >> specialPosition_;
+        }
+        else
+        {
+            is.read
+            (
+                reinterpret_cast<char*>(&Q_),
+                sizeof(Q_)
+              + sizeof(v_)
+              + sizeof(a_)
+              + sizeof(pi_)
+              + sizeof(tau_)
+              + sizeof(specialPosition_)
+              + sizeof(potentialEnergy_)
+              + sizeof(rf_)
+              + sizeof(special_)
+              + sizeof(id_)
+            );
+
+            is  >> siteForces_ >> sitePositions_;
+        }
+    }
+
+    // Check state of Istream
+    is.check
+    (
+        "Foam::molecule::molecule"
+        "(const Cloud<molecule>& cloud, Foam::Istream&), bool"
+    );
+}
+
+
+void Foam::molecule::readFields(Cloud<molecule>& mC)
+{
+    if (!mC.size())
+    {
+        return;
+    }
+
+    particle::readFields(mC);
+
+    IOField<tensor> Q(mC.fieldIOobject("Q", IOobject::MUST_READ));
+    mC.checkFieldIOobject(mC, Q);
+
+    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> pi(mC.fieldIOobject("pi", IOobject::MUST_READ));
+    mC.checkFieldIOobject(mC, pi);
+
+    IOField<vector> tau(mC.fieldIOobject("tau", IOobject::MUST_READ));
+    mC.checkFieldIOobject(mC, tau);
+
+    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(moleculeCloud, mC, iter)
+    {
+        molecule& mol = iter();
+
+        mol.Q_ = Q[i];
+        mol.v_ = v[i];
+        mol.a_ = a[i];
+        mol.pi_ = pi[i];
+        mol.tau_ = tau[i];
+        mol.specialPosition_ = specialPosition[i];
+        mol.special_ = special[i];
+        mol.id_ = id[i];
+        i++;
+    }
+}
+
+
+void Foam::molecule::writeFields(const Cloud<molecule>& mC)
+{
+    particle::writeFields(mC);
+
+    label np = mC.size();
+
+    IOField<tensor> Q(mC.fieldIOobject("Q", IOobject::NO_READ), np);
+    IOField<vector> v(mC.fieldIOobject("v", IOobject::NO_READ), np);
+    IOField<vector> a(mC.fieldIOobject("a", IOobject::NO_READ), np);
+    IOField<vector> pi(mC.fieldIOobject("pi", IOobject::NO_READ), np);
+    IOField<vector> tau(mC.fieldIOobject("tau", 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);
+
+    // Post processing fields
+
+    IOField<vector> piGlobal
+    (
+        mC.fieldIOobject("piGlobal", IOobject::NO_READ),
+        np
+    );
+
+    IOField<vector> tauGlobal
+    (
+        mC.fieldIOobject("tauGlobal", IOobject::NO_READ),
+        np
+    );
+
+    IOField<vector> orientation1
+    (
+        mC.fieldIOobject("orientation1", IOobject::NO_READ),
+        np
+    );
+
+    IOField<vector> orientation2
+    (
+        mC.fieldIOobject("orientation2", IOobject::NO_READ),
+        np
+    );
+
+    IOField<vector> orientation3
+    (
+        mC.fieldIOobject("orientation3", IOobject::NO_READ),
+        np
+    );
+
+    label i = 0;
+    forAllConstIter(moleculeCloud, mC, iter)
+    {
+        const molecule& mol = iter();
+
+        Q[i] = mol.Q_;
+        v[i] = mol.v_;
+        a[i] = mol.a_;
+        pi[i] = mol.pi_;
+        tau[i] = mol.tau_;
+        specialPosition[i] = mol.specialPosition_;
+        special[i] = mol.special_;
+        id[i] = mol.id_;
+
+        piGlobal[i] = mol.Q_ & mol.pi_;
+        tauGlobal[i] = mol.Q_ & mol.tau_;
+
+        orientation1[i] = mol.Q_ & vector(1,0,0);
+        orientation2[i] = mol.Q_ & vector(0,1,0);
+        orientation3[i] = mol.Q_ & vector(0,0,1);
+
+        i++;
+    }
+
+    Q.write();
+    v.write();
+    a.write();
+    pi.write();
+    tau.write();
+    specialPosition.write();
+    special.write();
+    id.write();
+
+    piGlobal.write();
+    tauGlobal.write();
+
+    orientation1.write();
+    orientation2.write();
+    orientation3.write();
+
+    Info<< "writeFields " << mC.name() << endl;
+
+    if (isA<moleculeCloud>(mC))
+    {
+        const moleculeCloud& m = dynamic_cast<const moleculeCloud&>(mC);
+
+        m.writeXYZ
+        (
+            m.mesh().time().timePath()/cloud::prefix/"moleculeCloud.xmol"
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const molecule& 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.Q_
+            << token::SPACE << mol.v_
+            << token::SPACE << mol.a_
+            << token::SPACE << mol.pi_
+            << token::SPACE << mol.tau_
+            << 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.Q_),
+            sizeof(mol.Q_)
+          + sizeof(mol.v_)
+          + sizeof(mol.a_)
+          + sizeof(mol.pi_)
+          + sizeof(mol.tau_)
+          + 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::molecule&)"
+    );
+
+    return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
new file mode 100644
index 0000000000..289735c170
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
@@ -0,0 +1,1250 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "moleculeCloud.H"
+#include "fvMesh.H"
+#include "mathematicalConstants.H"
+
+using namespace Foam::constant::mathematical;
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebug(Cloud<molecule>, 0);
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::moleculeCloud::buildConstProps()
+{
+    Info<< nl << "Reading moleculeProperties dictionary." << endl;
+
+    const List<word>& idList(pot_.idList());
+
+    constPropList_.setSize(idList.size());
+
+    const List<word>& siteIdList(pot_.siteIdList());
+
+    IOdictionary moleculePropertiesDict
+    (
+        IOobject
+        (
+            "moleculeProperties",
+            mesh_.time().constant(),
+            mesh_,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE,
+            false
+        )
+    );
+
+    forAll(idList, i)
+    {
+        const word& id = idList[i];
+        const dictionary& molDict = moleculePropertiesDict.subDict(id);
+
+        List<word> siteIdNames = molDict.lookup("siteIds");
+
+        List<label> siteIds(siteIdNames.size());
+
+        forAll(siteIdNames, sI)
+        {
+            const word& siteId = siteIdNames[sI];
+
+            siteIds[sI] = findIndex(siteIdList, siteId);
+
+            if (siteIds[sI] == -1)
+            {
+                FatalErrorIn("moleculeCloud::buildConstProps()")
+                    << siteId << " site not found."
+                    << nl << abort(FatalError);
+            }
+        }
+
+        molecule::constantProperties& constProp = constPropList_[i];
+
+        constProp = molecule::constantProperties(molDict);
+
+        constProp.siteIds() = siteIds;
+    }
+}
+
+
+void Foam::moleculeCloud::setSiteSizesAndPositions()
+{
+    forAllIter(moleculeCloud, *this, mol)
+    {
+        const molecule::constantProperties& cP = constProps(mol().id());
+
+        mol().setSiteSizes(cP.nSites());
+
+        mol().setSitePositions(cP);
+    }
+}
+
+
+void Foam::moleculeCloud::buildCellOccupancy()
+{
+    forAll(cellOccupancy_, cO)
+    {
+        cellOccupancy_[cO].clear();
+    }
+
+    forAllIter(moleculeCloud, *this, mol)
+    {
+        cellOccupancy_[mol().cell()].append(&mol());
+    }
+
+    forAll(cellOccupancy_, cO)
+    {
+        cellOccupancy_[cO].shrink();
+    }
+}
+
+
+void Foam::moleculeCloud::calculatePairForce()
+{
+    PstreamBuffers pBufs(Pstream::nonBlocking);
+
+    // Start sending referred data
+    il_.sendReferredData(cellOccupancy(), pBufs);
+
+    molecule* molI = NULL;
+    molecule* molJ = NULL;
+
+    {
+        // Real-Real interactions
+
+        const labelListList& dil = il_.dil();
+
+        forAll(dil, d)
+        {
+            forAll(cellOccupancy_[d],cellIMols)
+            {
+                molI = cellOccupancy_[d][cellIMols];
+
+                forAll(dil[d], interactingCells)
+                {
+                    List<molecule*> cellJ =
+                        cellOccupancy_[dil[d][interactingCells]];
+
+                    forAll(cellJ, cellJMols)
+                    {
+                        molJ = cellJ[cellJMols];
+
+                        evaluatePair(*molI, *molJ);
+                    }
+                }
+
+                forAll(cellOccupancy_[d], cellIOtherMols)
+                {
+                    molJ = cellOccupancy_[d][cellIOtherMols];
+
+                    if (molJ > molI)
+                    {
+                        evaluatePair(*molI, *molJ);
+                    }
+                }
+            }
+        }
+    }
+
+    // Receive referred data
+    il_.receiveReferredData(pBufs);
+
+    {
+        // Real-Referred interactions
+
+        const labelListList& ril = il_.ril();
+
+        List<IDLList<molecule> >& referredMols = il_.referredParticles();
+
+        forAll(ril, r)
+        {
+            const List<label>& realCells = ril[r];
+
+            IDLList<molecule>& refMols = referredMols[r];
+
+            forAllIter
+            (
+                IDLList<molecule>,
+                refMols,
+                refMol
+            )
+            {
+                forAll(realCells, rC)
+                {
+                    List<molecule*> cellI = cellOccupancy_[realCells[rC]];
+
+                    forAll(cellI, cellIMols)
+                    {
+                        molI = cellI[cellIMols];
+
+                        evaluatePair(*molI, refMol());
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::moleculeCloud::calculateTetherForce()
+{
+    const tetherPotentialList& tetherPot(pot_.tetherPotentials());
+
+    forAllIter(moleculeCloud, *this, mol)
+    {
+        if (mol().tethered())
+        {
+            vector rIT = mol().position() - mol().specialPosition();
+
+            label idI = mol().id();
+
+            scalar massI = constProps(idI).mass();
+
+            vector fIT = tetherPot.force(idI, rIT);
+
+            mol().a() += fIT/massI;
+
+            mol().potentialEnergy() += tetherPot.energy(idI, rIT);
+
+            // What to do here?
+            // mol().rf() += rIT*fIT;
+        }
+    }
+}
+
+
+void Foam::moleculeCloud::calculateExternalForce()
+{
+    forAllIter(moleculeCloud, *this, mol)
+    {
+        mol().a() += pot_.gravity();
+    }
+}
+
+
+void Foam::moleculeCloud::removeHighEnergyOverlaps()
+{
+    Info<< nl << "Removing high energy overlaps, limit = "
+        << pot_.potentialEnergyLimit()
+        << nl << "Removal order:";
+
+    forAll(pot_.removalOrder(), rO)
+    {
+        Info<< ' ' << pot_.idList()[pot_.removalOrder()[rO]];
+    }
+
+    Info<< nl ;
+
+    label initialSize = this->size();
+
+    buildCellOccupancy();
+
+    // Real-Real interaction
+
+    molecule* molI = NULL;
+    molecule* molJ = NULL;
+
+    {
+        DynamicList<molecule*> molsToDelete;
+
+        const labelListList& dil(il_.dil());
+
+        forAll(dil, d)
+        {
+            forAll(cellOccupancy_[d],cellIMols)
+            {
+                molI = cellOccupancy_[d][cellIMols];
+
+                forAll(dil[d], interactingCells)
+                {
+                    List<molecule*> cellJ =
+                        cellOccupancy_[dil[d][interactingCells]];
+
+                    forAll(cellJ, cellJMols)
+                    {
+                        molJ = cellJ[cellJMols];
+
+                        if (evaluatePotentialLimit(*molI, *molJ))
+                        {
+                            label idI = molI->id();
+
+                            label idJ = molJ->id();
+
+                            if
+                            (
+                                idI == idJ
+                             || findIndex(pot_.removalOrder(), idJ)
+                              < findIndex(pot_.removalOrder(), idI)
+                            )
+                            {
+                                if (findIndex(molsToDelete, molJ) == -1)
+                                {
+                                    molsToDelete.append(molJ);
+                                }
+                            }
+                            else if (findIndex(molsToDelete, molI) == -1)
+                            {
+                                molsToDelete.append(molI);
+                            }
+                        }
+                    }
+                }
+            }
+
+            forAll(cellOccupancy_[d], cellIOtherMols)
+            {
+                molJ = cellOccupancy_[d][cellIOtherMols];
+
+                if (molJ > molI)
+                {
+                    if (evaluatePotentialLimit(*molI, *molJ))
+                    {
+                        label idI = molI->id();
+
+                        label idJ = molJ->id();
+
+                        if
+                        (
+                            idI == idJ
+                         || findIndex(pot_.removalOrder(), idJ)
+                          < findIndex(pot_.removalOrder(), idI)
+                        )
+                        {
+                            if (findIndex(molsToDelete, molJ) == -1)
+                            {
+                                molsToDelete.append(molJ);
+                            }
+                        }
+                        else if (findIndex(molsToDelete, molI) == -1)
+                        {
+                            molsToDelete.append(molI);
+                        }
+                    }
+                }
+            }
+        }
+
+        forAll(molsToDelete, mTD)
+        {
+            deleteParticle(*(molsToDelete[mTD]));
+        }
+    }
+
+    buildCellOccupancy();
+
+    PstreamBuffers pBufs(Pstream::nonBlocking);
+
+    // Start sending referred data
+    il_.sendReferredData(cellOccupancy(), pBufs);
+
+        // Receive referred data
+    il_.receiveReferredData(pBufs);
+
+    // Real-Referred interaction
+
+    {
+        DynamicList<molecule*> molsToDelete;
+
+        const labelListList& ril(il_.ril());
+
+        List<IDLList<molecule> >& referredMols = il_.referredParticles();
+
+        forAll(ril, r)
+        {
+            IDLList<molecule>& refMols = referredMols[r];
+
+            forAllIter
+            (
+                IDLList<molecule>,
+                refMols,
+                refMol
+            )
+            {
+                molJ = &refMol();
+
+                const List<label>& realCells = ril[r];
+
+                forAll(realCells, rC)
+                {
+                    label cellI = realCells[rC];
+
+                    List<molecule*> cellIMols = cellOccupancy_[cellI];
+
+                    forAll(cellIMols, cIM)
+                    {
+                        molI = cellIMols[cIM];
+
+                        if (evaluatePotentialLimit(*molI, *molJ))
+                        {
+                            label idI = molI->id();
+
+                            label idJ = molJ->id();
+
+                            if
+                            (
+                                findIndex(pot_.removalOrder(), idI)
+                              < findIndex(pot_.removalOrder(), idJ)
+                            )
+                            {
+                                if (findIndex(molsToDelete, molI) == -1)
+                                {
+                                    molsToDelete.append(molI);
+                                }
+                            }
+                            else if
+                            (
+                                findIndex(pot_.removalOrder(), idI)
+                             == findIndex(pot_.removalOrder(), idJ)
+                            )
+                            {
+                                // Remove one of the molecules
+                                // arbitrarily, assuring that a
+                                // consistent decision is made for
+                                // both real-referred pairs.
+
+                                if (molI->origId() > molJ->origId())
+                                {
+                                    if (findIndex(molsToDelete, molI) == -1)
+                                    {
+                                        molsToDelete.append(molI);
+                                    }
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+
+        forAll(molsToDelete, mTD)
+        {
+            deleteParticle(*(molsToDelete[mTD]));
+        }
+    }
+
+    buildCellOccupancy();
+
+    // Start sending referred data
+    il_.sendReferredData(cellOccupancy(), pBufs);
+
+    // Receive referred data
+    il_.receiveReferredData(pBufs);
+
+    label molsRemoved = initialSize - this->size();
+
+    if (Pstream::parRun())
+    {
+        reduce(molsRemoved, sumOp<label>());
+    }
+
+    Info<< tab << molsRemoved << " molecules removed" << endl;
+}
+
+
+void Foam::moleculeCloud::initialiseMolecules
+(
+    const IOdictionary& mdInitialiseDict
+)
+{
+    Info<< nl
+        << "Initialising molecules in each zone specified in mdInitialiseDict."
+        << endl;
+
+    const cellZoneMesh& cellZones = mesh_.cellZones();
+
+    if (!cellZones.size())
+    {
+        FatalErrorIn("void Foam::moleculeCloud::initialiseMolecules")
+            << "No cellZones found in the mesh."
+            << abort(FatalError);
+    }
+
+    forAll(cellZones, z)
+    {
+        const cellZone& zone(cellZones[z]);
+
+        if (zone.size())
+        {
+            if (!mdInitialiseDict.found(zone.name()))
+            {
+                Info<< "No specification subDictionary for zone "
+                    << zone.name() << " found, skipping." << endl;
+            }
+            else
+            {
+                const dictionary& zoneDict =
+                    mdInitialiseDict.subDict(zone.name());
+
+                const scalar temperature
+                (
+                    readScalar(zoneDict.lookup("temperature"))
+                );
+
+                const vector bulkVelocity(zoneDict.lookup("bulkVelocity"));
+
+                List<word> latticeIds
+                (
+                    zoneDict.lookup("latticeIds")
+                );
+
+                List<vector> latticePositions
+                (
+                    zoneDict.lookup("latticePositions")
+                );
+
+                if (latticeIds.size() != latticePositions.size())
+                {
+                    FatalErrorIn("Foam::moleculeCloud::initialiseMolecules")
+                        << "latticeIds and latticePositions must be the same "
+                        << " size." << nl
+                        << abort(FatalError);
+                }
+
+                diagTensor latticeCellShape
+                (
+                    zoneDict.lookup("latticeCellShape")
+                );
+
+                scalar latticeCellScale = 0.0;
+
+                if (zoneDict.found("numberDensity"))
+                {
+                    scalar numberDensity = readScalar
+                    (
+                        zoneDict.lookup("numberDensity")
+                    );
+
+                    if (numberDensity < VSMALL)
+                    {
+                        WarningIn("moleculeCloud::initialiseMolecules")
+                            << "numberDensity too small, not filling zone "
+                            << zone.name() << endl;
+
+                        continue;
+                    }
+
+                    latticeCellScale = pow
+                    (
+                        latticeIds.size()/(det(latticeCellShape)*numberDensity),
+                        (1.0/3.0)
+                    );
+                }
+                else if (zoneDict.found("massDensity"))
+                {
+                    scalar unitCellMass = 0.0;
+
+                    forAll(latticeIds, i)
+                    {
+                        label id = findIndex(pot_.idList(), latticeIds[i]);
+
+                        const molecule::constantProperties& cP(constProps(id));
+
+                        unitCellMass += cP.mass();
+                    }
+
+                    scalar massDensity = readScalar
+                    (
+                        zoneDict.lookup("massDensity")
+                    );
+
+                    if (massDensity < VSMALL)
+                    {
+                        WarningIn("moleculeCloud::initialiseMolecules")
+                            << "massDensity too small, not filling zone "
+                            << zone.name() << endl;
+
+                        continue;
+                    }
+
+
+                    latticeCellScale = pow
+                    (
+                        unitCellMass/(det(latticeCellShape)*massDensity),
+                        (1.0/3.0)
+                    );
+                }
+                else
+                {
+                    FatalErrorIn("Foam::moleculeCloud::initialiseMolecules")
+                        << "massDensity or numberDensity not specified " << nl
+                        << abort(FatalError);
+                }
+
+                latticeCellShape *= latticeCellScale;
+
+                vector anchor(zoneDict.lookup("anchor"));
+
+                bool tethered = false;
+
+                if (zoneDict.found("tetherSiteIds"))
+                {
+                    tethered = bool
+                    (
+                        List<word>(zoneDict.lookup("tetherSiteIds")).size()
+                    );
+                }
+
+                const vector orientationAngles
+                (
+                    zoneDict.lookup("orientationAngles")
+                );
+
+                scalar phi(orientationAngles.x()*pi/180.0);
+
+                scalar theta(orientationAngles.y()*pi/180.0);
+
+                scalar psi(orientationAngles.z()*pi/180.0);
+
+                const tensor R
+                (
+                    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)
+                );
+
+                // Find the optimal anchor position.  Finding the approximate
+                // mid-point of the zone of cells and snapping to the nearest
+                // lattice location.
+
+                vector zoneMin = VGREAT*vector::one;
+
+                vector zoneMax = -VGREAT*vector::one;
+
+                forAll(zone, cell)
+                {
+                    const point cellCentre = mesh_.cellCentres()[zone[cell]];
+
+                    if (cellCentre.x() > zoneMax.x())
+                    {
+                        zoneMax.x() = cellCentre.x();
+                    }
+                    if (cellCentre.x() < zoneMin.x())
+                    {
+                        zoneMin.x() = cellCentre.x();
+                    }
+                    if (cellCentre.y() > zoneMax.y())
+                    {
+                        zoneMax.y() = cellCentre.y();
+                    }
+                    if (cellCentre.y() < zoneMin.y())
+                    {
+                        zoneMin.y() = cellCentre.y();
+                    }
+                    if (cellCentre.z() > zoneMax.z())
+                    {
+                        zoneMax.z() = cellCentre.z();
+                    }
+                    if (cellCentre.z() < zoneMin.z())
+                    {
+                        zoneMin.z() = cellCentre.z();
+                    }
+                }
+
+                point zoneMid = 0.5*(zoneMin + zoneMax);
+
+                point latticeMid = inv(latticeCellShape) & (R.T()
+                & (zoneMid - anchor));
+
+                point latticeAnchor
+                (
+                    label(latticeMid.x() + 0.5*sign(latticeMid.x())),
+                    label(latticeMid.y() + 0.5*sign(latticeMid.y())),
+                    label(latticeMid.z() + 0.5*sign(latticeMid.z()))
+                );
+
+                anchor += (R & (latticeCellShape & latticeAnchor));
+
+                // Continue trying to place molecules 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.
+
+                label n = 0;
+
+                label totalZoneMols = 0;
+
+                label molsPlacedThisIteration = 0;
+
+                while
+                (
+                    molsPlacedThisIteration != 0
+                 || totalZoneMols == 0
+                )
+                {
+                    label sizeBeforeIteration = this->size();
+
+                    bool partOfLayerInBounds = false;
+
+                    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+                    // start of placement of molecules
+                    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+                    if (n == 0)
+                    {
+                        // Special treatment is required for the first position,
+                        // i.e. iteration zero.
+
+                        labelVector unitCellLatticePosition(0,0,0);
+
+                        forAll(latticePositions, p)
+                        {
+                            label id = findIndex(pot_.idList(), latticeIds[p]);
+
+                            const vector& latticePosition =
+                                vector
+                                (
+                                    unitCellLatticePosition.x(),
+                                    unitCellLatticePosition.y(),
+                                    unitCellLatticePosition.z()
+                                )
+                              + latticePositions[p];
+
+                            point globalPosition =
+                                anchor
+                              + (R & (latticeCellShape & latticePosition));
+
+                            partOfLayerInBounds = mesh_.bounds().contains
+                            (
+                                globalPosition
+                            );
+
+                            label cell = -1;
+                            label tetFace = -1;
+                            label tetPt = -1;
+
+                            mesh_.findCellFacePt
+                            (
+                                globalPosition,
+                                cell,
+                                tetFace,
+                                tetPt
+                            );
+
+                            if (findIndex(zone, cell) != -1)
+                            {
+                                createMolecule
+                                (
+                                    globalPosition,
+                                    cell,
+                                    tetFace,
+                                    tetPt,
+                                    id,
+                                    tethered,
+                                    temperature,
+                                    bulkVelocity
+                                );
+                            }
+                        }
+                    }
+                    else
+                    {
+                        // Place top and bottom caps.
+
+                        labelVector unitCellLatticePosition(0,0,0);
+
+                        for
+                        (
+                            unitCellLatticePosition.z() = -n;
+                            unitCellLatticePosition.z() <= n;
+                            unitCellLatticePosition.z() += 2*n
+                        )
+                        {
+                            for
+                            (
+                                unitCellLatticePosition.y() = -n;
+                                unitCellLatticePosition.y() <= n;
+                                unitCellLatticePosition.y()++
+                            )
+                            {
+                                for
+                                (
+                                    unitCellLatticePosition.x() = -n;
+                                    unitCellLatticePosition.x() <= n;
+                                    unitCellLatticePosition.x()++
+                                )
+                                {
+                                    forAll(latticePositions, p)
+                                    {
+                                        label id = findIndex
+                                        (
+                                            pot_.idList(),
+                                            latticeIds[p]
+                                        );
+
+                                        const vector& latticePosition =
+                                            vector
+                                            (
+                                                unitCellLatticePosition.x(),
+                                                unitCellLatticePosition.y(),
+                                                unitCellLatticePosition.z()
+                                            )
+                                          + latticePositions[p];
+
+                                        point globalPosition =
+                                            anchor
+                                          + (
+                                                R
+                                              & (
+                                                    latticeCellShape
+                                                  & latticePosition
+                                                )
+                                            );
+
+                                        partOfLayerInBounds =
+                                            mesh_.bounds().contains
+                                            (
+                                                globalPosition
+                                            );
+
+                                        label cell = -1;
+                                        label tetFace = -1;
+                                        label tetPt = -1;
+
+                                        mesh_.findCellFacePt
+                                        (
+                                            globalPosition,
+                                            cell,
+                                            tetFace,
+                                            tetPt
+                                        );
+
+                                        if (findIndex(zone, cell) != -1)
+                                        {
+                                            createMolecule
+                                            (
+                                                globalPosition,
+                                                cell,
+                                                tetFace,
+                                                tetPt,
+                                                id,
+                                                tethered,
+                                                temperature,
+                                                bulkVelocity
+                                            );
+                                        }
+                                    }
+                                }
+                            }
+                        }
+
+                        for
+                        (
+                            unitCellLatticePosition.z() = -(n-1);
+                            unitCellLatticePosition.z() <= (n-1);
+                            unitCellLatticePosition.z()++
+                        )
+                        {
+                            for (label iR = 0; iR <= 2*n -1; iR++)
+                            {
+                                unitCellLatticePosition.x() = n;
+
+                                unitCellLatticePosition.y() = -n + (iR + 1);
+
+                                for (label iK = 0; iK < 4; iK++)
+                                {
+                                    forAll(latticePositions, p)
+                                    {
+                                        label id = findIndex
+                                        (
+                                            pot_.idList(),
+                                            latticeIds[p]
+                                        );
+
+                                        const vector& latticePosition =
+                                            vector
+                                            (
+                                                unitCellLatticePosition.x(),
+                                                unitCellLatticePosition.y(),
+                                                unitCellLatticePosition.z()
+                                            )
+                                          + latticePositions[p];
+
+                                        point globalPosition =
+                                            anchor
+                                          + (
+                                                R
+                                              & (
+                                                   latticeCellShape
+                                                 & latticePosition
+                                                )
+                                            );
+
+                                        partOfLayerInBounds =
+                                            mesh_.bounds().contains
+                                            (
+                                                globalPosition
+                                            );
+
+                                        label cell = -1;
+                                        label tetFace = -1;
+                                        label tetPt = -1;
+
+                                        mesh_.findCellFacePt
+                                        (
+                                            globalPosition,
+                                            cell,
+                                            tetFace,
+                                            tetPt
+                                        );
+
+                                        if (findIndex(zone, cell) != -1)
+                                        {
+                                            createMolecule
+                                            (
+                                                globalPosition,
+                                                cell,
+                                                tetFace,
+                                                tetPt,
+                                                id,
+                                                tethered,
+                                                temperature,
+                                                bulkVelocity
+                                            );
+                                        }
+                                    }
+
+                                    unitCellLatticePosition =
+                                        labelVector
+                                        (
+                                          - unitCellLatticePosition.y(),
+                                            unitCellLatticePosition.x(),
+                                            unitCellLatticePosition.z()
+                                        );
+                                }
+                            }
+                        }
+                    }
+
+                    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+                    // end of placement of molecules
+                    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+                    if
+                    (
+                        totalZoneMols == 0
+                     && !partOfLayerInBounds
+                    )
+                    {
+                        WarningIn("Foam::moleculeCloud::initialiseMolecules()")
+                            << "A whole layer of unit cells was placed "
+                            << "outside the bounds of the mesh, but no "
+                            << "molecules have been placed in zone '"
+                            << zone.name()
+                            << "'.  This is likely to be because the zone "
+                            << "has few cells ("
+                            << zone.size()
+                            << " in this case) and no lattice position "
+                            << "fell inside them.  "
+                            << "Aborting filling this zone."
+                            << endl;
+
+                        break;
+                    }
+
+                    molsPlacedThisIteration =
+                        this->size() - sizeBeforeIteration;
+
+                    totalZoneMols += molsPlacedThisIteration;
+
+                    n++;
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::moleculeCloud::createMolecule
+(
+    const point& position,
+    label cell,
+    label tetFace,
+    label tetPt,
+    label id,
+    bool tethered,
+    scalar temperature,
+    const vector& bulkVelocity
+)
+{
+    if (cell == -1)
+    {
+        mesh_.findCellFacePt(position, cell, tetFace, tetPt);
+    }
+
+    if (cell == -1)
+    {
+        FatalErrorIn("Foam::moleculeCloud::createMolecule")
+            << "Position specified does not correspond to a mesh cell." << nl
+            << abort(FatalError);
+    }
+
+    point specialPosition(vector::zero);
+
+    label special = 0;
+
+    if (tethered)
+    {
+        specialPosition = position;
+
+        special = molecule::SPECIAL_TETHERED;
+    }
+
+    const molecule::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)
+        );
+    }
+
+    addParticle
+    (
+        new molecule
+        (
+            mesh_,
+            position,
+            cell,
+            tetFace,
+            tetPt,
+            Q,
+            v,
+            vector::zero,
+            pi,
+            vector::zero,
+            specialPosition,
+            constProps(id),
+            special,
+            id
+        )
+    );
+}
+
+
+Foam::label Foam::moleculeCloud::nSites() const
+{
+    label n = 0;
+
+    forAllConstIter(moleculeCloud, *this, mol)
+    {
+        n += constProps(mol().id()).nSites();
+    }
+
+    return n;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::moleculeCloud::moleculeCloud
+(
+    const polyMesh& mesh,
+    const potential& pot,
+    bool readFields
+)
+:
+    Cloud<molecule>(mesh, "moleculeCloud", false),
+    mesh_(mesh),
+    pot_(pot),
+    cellOccupancy_(mesh_.nCells()),
+    il_(mesh_, pot_.pairPotentials().rCutMax(), false),
+    constPropList_(),
+    rndGen_(clock::getTime())
+{
+    if (readFields)
+    {
+        molecule::readFields(*this);
+    }
+
+    buildConstProps();
+
+    setSiteSizesAndPositions();
+
+    removeHighEnergyOverlaps();
+
+    calculateForce();
+}
+
+
+Foam::moleculeCloud::moleculeCloud
+(
+    const polyMesh& mesh,
+    const potential& pot,
+    const IOdictionary& mdInitialiseDict,
+    bool readFields
+)
+:
+    Cloud<molecule>(mesh, "moleculeCloud", false),
+    mesh_(mesh),
+    pot_(pot),
+    il_(mesh_, 0.0, false),
+    constPropList_(),
+    rndGen_(clock::getTime())
+{
+    if (readFields)
+    {
+        molecule::readFields(*this);
+    }
+
+    clear();
+
+    buildConstProps();
+
+    initialiseMolecules(mdInitialiseDict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::moleculeCloud::evolve()
+{
+    molecule::trackingData td0(*this, 0);
+    Cloud<molecule>::move(td0, mesh_.time().deltaTValue());
+
+    molecule::trackingData td1(*this, 1);
+    Cloud<molecule>::move(td1, mesh_.time().deltaTValue());
+
+    molecule::trackingData td2(*this, 2);
+    Cloud<molecule>::move(td2, mesh_.time().deltaTValue());
+
+    calculateForce();
+
+    molecule::trackingData td3(*this, 3);
+    Cloud<molecule>::move(td3, mesh_.time().deltaTValue());
+}
+
+
+void Foam::moleculeCloud::calculateForce()
+{
+    buildCellOccupancy();
+
+    // Set accumulated quantities to zero
+    forAllIter(moleculeCloud, *this, mol)
+    {
+        mol().siteForces() = vector::zero;
+
+        mol().potentialEnergy() = 0.0;
+
+        mol().rf() = tensor::zero;
+    }
+
+    calculatePairForce();
+
+    calculateTetherForce();
+
+    calculateExternalForce();
+}
+
+
+void Foam::moleculeCloud::applyConstraintsAndThermostats
+(
+    const scalar targetTemperature,
+    const scalar measuredTemperature
+)
+{
+    scalar temperatureCorrectionFactor =
+        sqrt(targetTemperature/max(VSMALL, measuredTemperature));
+
+    Info<< "----------------------------------------" << nl
+        << "Temperature equilibration" << nl
+        << "Target temperature = "
+        << targetTemperature << nl
+        << "Measured temperature = "
+        << measuredTemperature << nl
+        << "Temperature correction factor = "
+        << temperatureCorrectionFactor << nl
+        << "----------------------------------------"
+        << endl;
+
+    forAllIter(moleculeCloud, *this, mol)
+    {
+        mol().v() *= temperatureCorrectionFactor;
+
+        mol().pi() *= temperatureCorrectionFactor;
+    }
+}
+
+
+void Foam::moleculeCloud::writeXYZ(const fileName& fName) const
+{
+    OFstream os(fName);
+
+    os  << nSites() << nl << "moleculeCloud site positions in angstroms" << nl;
+
+    forAllConstIter(moleculeCloud, *this, mol)
+    {
+        const molecule::constantProperties& cP = constProps(mol().id());
+
+        forAll(mol().sitePositions(), i)
+        {
+            const point& sP = mol().sitePositions()[i];
+
+            os  << pot_.siteIdList()[cP.siteIds()[i]]
+                << ' ' << sP.x()*1e10
+                << ' ' << sP.y()*1e10
+                << ' ' << sP.z()*1e10
+                << nl;
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
new file mode 100644
index 0000000000..7cf3359a35
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
@@ -0,0 +1,219 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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
+
+SourceFiles
+    moleculeCloudI.H
+    moleculeCloud.C
+
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef moleculeCloud_H
+#define moleculeCloud_H
+
+#include "Cloud.H"
+#include "molecule.H"
+#include "IOdictionary.H"
+#include "potential.H"
+#include "InteractionLists.H"
+#include "labelVector.H"
+#include "Random.H"
+#include "fileName.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class moleculeCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+class moleculeCloud
+:
+    public Cloud<molecule>
+{
+
+private:
+
+    // Private data
+
+        const polyMesh& mesh_;
+
+        const potential& pot_;
+
+        List<DynamicList<molecule*> > cellOccupancy_;
+
+        InteractionLists<molecule> il_;
+
+        List<molecule::constantProperties> constPropList_;
+
+        Random rndGen_;
+
+
+    // Private Member Functions
+
+        void buildConstProps();
+
+        void setSiteSizesAndPositions();
+
+        //- Determine which molecules are in which cells
+        void buildCellOccupancy();
+
+        void calculatePairForce();
+
+        inline void evaluatePair
+        (
+            molecule& molI,
+            molecule& molJ
+        );
+
+        inline bool evaluatePotentialLimit
+        (
+            molecule& molI,
+            molecule& molJ
+        ) const;
+
+        void calculateTetherForce();
+
+        void calculateExternalForce();
+
+        void removeHighEnergyOverlaps();
+
+        void initialiseMolecules
+        (
+            const IOdictionary& mdInitialiseDict
+        );
+
+        void createMolecule
+        (
+            const point& position,
+            label cell,
+            label tetFace,
+            label tetPt,
+            label id,
+            bool tethered,
+            scalar temperature,
+            const vector& bulkVelocity
+        );
+
+        label nSites() const;
+
+        inline vector equipartitionLinearVelocity
+        (
+            scalar temperature,
+            scalar mass
+        );
+
+        inline vector equipartitionAngularMomentum
+        (
+            scalar temperature,
+            const molecule::constantProperties& cP
+        );
+
+        //- Disallow default bitwise copy construct
+        moleculeCloud(const moleculeCloud&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const moleculeCloud&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct given mesh and potential references
+        moleculeCloud
+        (
+            const polyMesh& mesh,
+            const potential& pot,
+            bool readFields = true
+        );
+
+        //- Construct given mesh, potential and mdInitialiseDict
+        moleculeCloud
+        (
+            const polyMesh& mesh,
+            const potential& pot,
+            const IOdictionary& mdInitialiseDict,
+            bool readFields = true
+        );
+
+
+    // Member Functions
+
+        //- Evolve the molecules (move, calculate forces, control state etc)
+        void evolve();
+
+        void calculateForce();
+
+        void applyConstraintsAndThermostats
+        (
+            const scalar targetTemperature,
+            const scalar measuredTemperature
+        );
+
+
+        // Access
+
+            inline const polyMesh& mesh() const;
+
+            inline const potential& pot() const;
+
+            inline const List<DynamicList<molecule*> >& cellOccupancy() const;
+
+            inline const InteractionLists<molecule>& il() const;
+
+            inline const List<molecule::constantProperties> constProps() const;
+
+            inline const molecule::constantProperties&
+                constProps(label id) const;
+
+            inline Random& rndGen();
+
+
+    // Member Operators
+
+        //- Write molecule sites in XYZ format
+        void writeXYZ(const fileName& fName) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "moleculeCloudI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
new file mode 100644
index 0000000000..dbeee755a5
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
@@ -0,0 +1,390 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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 "constants.H"
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+inline void Foam::moleculeCloud::evaluatePair
+(
+    molecule& molI,
+    molecule& molJ
+)
+{
+    const pairPotentialList& pairPot = pot_.pairPotentials();
+
+    const pairPotential& electrostatic = pairPot.electrostatic();
+
+    label idI = molI.id();
+
+    label idJ = molJ.id();
+
+    const molecule::constantProperties& constPropI(constProps(idI));
+
+    const molecule::constantProperties& constPropJ(constProps(idJ));
+
+    List<label> siteIdsI = constPropI.siteIds();
+
+    List<label> siteIdsJ = constPropJ.siteIds();
+
+    List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
+
+    List<bool> electrostaticSitesI = constPropI.electrostaticSites();
+
+    List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
+
+    List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
+
+    forAll(siteIdsI, sI)
+    {
+        label idsI(siteIdsI[sI]);
+
+        forAll(siteIdsJ, sJ)
+        {
+            label idsJ(siteIdsJ[sJ]);
+
+            if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
+            {
+                vector rsIsJ =
+                    molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
+
+                scalar rsIsJMagSq = magSqr(rsIsJ);
+
+                if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
+                {
+                    scalar rsIsJMag = mag(rsIsJ);
+
+                    vector fsIsJ =
+                        (rsIsJ/rsIsJMag)
+                       *pairPot.force(idsI, idsJ, rsIsJMag);
+
+                    molI.siteForces()[sI] += fsIsJ;
+
+                    molJ.siteForces()[sJ] += -fsIsJ;
+
+                    scalar potentialEnergy
+                    (
+                        pairPot.energy(idsI, idsJ, rsIsJMag)
+                    );
+
+                    molI.potentialEnergy() += 0.5*potentialEnergy;
+
+                    molJ.potentialEnergy() += 0.5*potentialEnergy;
+
+                    vector rIJ = molI.position() - molJ.position();
+
+                    tensor virialContribution =
+                        (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
+
+                    molI.rf() += virialContribution;
+
+                    molJ.rf() += virialContribution;
+                }
+            }
+
+            if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
+            {
+                vector rsIsJ =
+                molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
+
+                scalar rsIsJMagSq = magSqr(rsIsJ);
+
+                if (rsIsJMagSq <= electrostatic.rCutSqr())
+                {
+                    scalar rsIsJMag = mag(rsIsJ);
+
+                    scalar chargeI = constPropI.siteCharges()[sI];
+
+                    scalar chargeJ = constPropJ.siteCharges()[sJ];
+
+                    vector fsIsJ =
+                        (rsIsJ/rsIsJMag)
+                       *chargeI*chargeJ*electrostatic.force(rsIsJMag);
+
+                    molI.siteForces()[sI] += fsIsJ;
+
+                    molJ.siteForces()[sJ] += -fsIsJ;
+
+                    scalar potentialEnergy =
+                        chargeI*chargeJ
+                       *electrostatic.energy(rsIsJMag);
+
+                    molI.potentialEnergy() += 0.5*potentialEnergy;
+
+                    molJ.potentialEnergy() += 0.5*potentialEnergy;
+
+                    vector rIJ = molI.position() - molJ.position();
+
+                    tensor virialContribution =
+                        (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
+
+                    molI.rf() += virialContribution;
+
+                    molJ.rf() += virialContribution;
+                }
+            }
+        }
+    }
+}
+
+
+inline bool Foam::moleculeCloud::evaluatePotentialLimit
+(
+    molecule& molI,
+    molecule& molJ
+) const
+{
+    const pairPotentialList& pairPot = pot_.pairPotentials();
+
+    const pairPotential& electrostatic = pairPot.electrostatic();
+
+    label idI = molI.id();
+
+    label idJ = molJ.id();
+
+    const molecule::constantProperties& constPropI(constProps(idI));
+
+    const molecule::constantProperties& constPropJ(constProps(idJ));
+
+    List<label> siteIdsI = constPropI.siteIds();
+
+    List<label> siteIdsJ = constPropJ.siteIds();
+
+    List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
+
+    List<bool> electrostaticSitesI = constPropI.electrostaticSites();
+
+    List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
+
+    List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
+
+    forAll(siteIdsI, sI)
+    {
+        label idsI(siteIdsI[sI]);
+
+        forAll(siteIdsJ, sJ)
+        {
+            label idsJ(siteIdsJ[sJ]);
+
+            if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
+            {
+                vector rsIsJ =
+                    molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
+
+                scalar rsIsJMagSq = magSqr(rsIsJ);
+
+                if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
+                {
+                    scalar rsIsJMag = mag(rsIsJ);
+
+                    // Guard against pairPot.energy being evaluated
+                    // if rIJMag < SMALL. A floating point exception will
+                    // happen otherwise.
+
+                    if (rsIsJMag < SMALL)
+                    {
+                        WarningIn("moleculeCloud::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 molecules"
+                            << " twice. Removing one of the molecules."
+                            << endl;
+
+                        return true;
+                    }
+
+                    // Guard against pairPot.energy being evaluated if rIJMag <
+                    // rMin.  A tabulation lookup error will occur otherwise.
+
+                    if (rsIsJMag < pairPot.rMin(idsI, idsJ))
+                    {
+                        return true;
+                    }
+
+                    if
+                    (
+                        mag(pairPot.energy(idsI, idsJ, rsIsJMag))
+                      > pot_.potentialEnergyLimit()
+                    )
+                    {
+                        return true;
+                    };
+                }
+            }
+
+            if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
+            {
+                vector rsIsJ =
+                    molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
+
+                scalar rsIsJMagSq = magSqr(rsIsJ);
+
+                if (pairPot.rCutMaxSqr(rsIsJMagSq))
+                {
+                    scalar rsIsJMag = mag(rsIsJ);
+
+                    // Guard against pairPot.energy being evaluated
+                    // if rIJMag < SMALL. A floating point exception will
+                    // happen otherwise.
+
+                    if (rsIsJMag < SMALL)
+                    {
+                        WarningIn("moleculeCloud::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 molecules"
+                            << " twice. Removing one of the molecules."
+                            << endl;
+
+                        return true;
+                    }
+
+                    if (rsIsJMag < electrostatic.rMin())
+                    {
+                        return true;
+                    }
+
+                    scalar chargeI = constPropI.siteCharges()[sI];
+
+                    scalar chargeJ = constPropJ.siteCharges()[sJ];
+
+                    if
+                    (
+                        mag(chargeI*chargeJ*electrostatic.energy(rsIsJMag))
+                      > pot_.potentialEnergyLimit()
+                    )
+                    {
+                        return true;
+                    };
+                }
+            }
+        }
+    }
+
+    return false;
+}
+
+
+inline Foam::vector Foam::moleculeCloud::equipartitionLinearVelocity
+(
+    scalar temperature,
+    scalar mass
+)
+{
+    return sqrt(physicoChemical::k.value()*temperature/mass)*vector
+    (
+        rndGen_.GaussNormal(),
+        rndGen_.GaussNormal(),
+        rndGen_.GaussNormal()
+    );
+}
+
+
+inline Foam::vector Foam::moleculeCloud::equipartitionAngularMomentum
+(
+    scalar temperature,
+    const molecule::constantProperties& cP
+)
+{
+    scalar sqrtKbT = sqrt(physicoChemical::k.value()*temperature);
+
+    if (cP.linearMolecule())
+    {
+        return sqrtKbT*vector
+        (
+            0.0,
+            sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal(),
+            sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal()
+        );
+    }
+    else
+    {
+        return sqrtKbT*vector
+        (
+            sqrt(cP.momentOfInertia().xx())*rndGen_.GaussNormal(),
+            sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal(),
+            sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal()
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::polyMesh& Foam::moleculeCloud::mesh() const
+{
+    return mesh_;
+}
+
+
+inline const Foam::potential& Foam::moleculeCloud::pot() const
+{
+    return pot_;
+}
+
+
+inline const Foam::List<Foam::DynamicList<Foam::molecule*> >&
+    Foam::moleculeCloud::cellOccupancy() const
+{
+    return cellOccupancy_;
+}
+
+
+inline const Foam::InteractionLists<Foam::molecule>&
+    Foam::moleculeCloud::il() const
+{
+    return il_;
+}
+
+
+inline const Foam::List<Foam::molecule::constantProperties>
+    Foam::moleculeCloud::constProps() const
+{
+    return constPropList_;
+}
+
+
+inline const Foam::molecule::constantProperties&
+    Foam::moleculeCloud::constProps(label id) const
+{
+    return constPropList_[id];
+}
+
+
+inline Foam::Random& Foam::moleculeCloud::rndGen()
+{
+    return rndGen_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnits.C b/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnits.C
new file mode 100644
index 0000000000..cdcb1d4a65
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnits.C
@@ -0,0 +1,164 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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 "reducedUnits.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+const  Foam::scalar Foam::reducedUnits::kb = 1.3806504e-23;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::reducedUnits::calcRefValues()
+{
+    if
+    (
+        refTime_ < VSMALL
+     || refLength_ < VSMALL
+     || refMass_ < VSMALL
+    )
+    {
+        FatalErrorIn("Foam::reducedUnits::calcRefValues() ")
+            << "One of more referencence values too small for floating point "
+            << "calculation: "
+            << "refTime_ = " << refTime_
+            << ", refLength = " << refTemp_
+            << ", refMass = " << refMass_
+            << nl << abort(FatalError);
+    }
+
+    refEnergy_ = refLength_*refLength_*refMass_/(refTime_*refTime_);
+
+    refTemp_ = refEnergy_ / kb;
+
+    refForce_ = refEnergy_/refLength_;
+
+    refVelocity_ = Foam::sqrt(refEnergy_/refMass_);
+
+    refVolume_ = Foam::pow(refLength_,3.0);
+
+    refPressure_ = refEnergy_/refVolume_;
+
+    refMassDensity_ = refMass_/refVolume_;
+
+    refNumberDensity_ = 1.0/refVolume_;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::reducedUnits::reducedUnits()
+:
+    refLength_(1e-9),
+    refTime_(1e-12),
+    refMass_(1.660538782e-27)
+{
+    calcRefValues();
+}
+
+
+Foam::reducedUnits::reducedUnits
+(
+    scalar refLength,
+    scalar refTime,
+    scalar refMass
+)
+:
+    refLength_(refLength),
+    refTime_(refTime),
+    refMass_(refMass)
+{
+    calcRefValues();
+}
+
+
+Foam::reducedUnits::reducedUnits(const IOdictionary& reducedUnitsDict)
+:
+    refLength_(),
+    refTime_(),
+    refMass_()
+{
+    setRefValues(reducedUnitsDict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::reducedUnits::~reducedUnits()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::reducedUnits::setRefValues
+(
+    scalar refLength,
+    scalar refTime,
+    scalar refMass
+)
+{
+    refLength_ = refLength;
+
+    refTime_ = refTime;
+
+    refMass_ = refMass;
+
+    calcRefValues();
+}
+
+
+void Foam::reducedUnits::setRefValues
+(
+    const IOdictionary& reducedUnitsDict
+)
+{
+    refLength_ = readScalar(reducedUnitsDict.lookup("refLength"));
+
+    refTime_ = readScalar(reducedUnitsDict.lookup("refTime"));
+
+    refMass_  = readScalar(reducedUnitsDict.lookup("refMass"));
+
+    calcRefValues();
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+void Foam::reducedUnits::operator=(const reducedUnits& rhs)
+{
+    // Check for assignment to self
+    if (this == &rhs)
+    {
+        FatalErrorIn
+        (
+            "Foam::reducedUnits::operator=(const Foam::reducedUnits&)"
+        )   << "Attempted assignment to self"
+            << abort(FatalError);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnits.H b/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnits.H
new file mode 100644
index 0000000000..9445ff64da
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnits.H
@@ -0,0 +1,182 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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::reducedUnits
+
+Description
+
+SourceFiles
+    reducedUnitsI.H
+    reducedUnits.C
+    reducedUnitsIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef reducedUnits_H
+#define reducedUnits_H
+
+#include "scalar.H"
+#include "IOdictionary.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class reducedUnits Declaration
+\*---------------------------------------------------------------------------*/
+
+class reducedUnits
+{
+    // Private data
+
+        // Reduced units
+
+            // Fundamental values
+
+                scalar refLength_;
+
+                scalar refTime_;
+
+                scalar refMass_;
+
+            // Derived values
+
+                scalar refEnergy_;
+
+                scalar refTemp_;
+
+                scalar refForce_;
+
+                scalar refVelocity_;
+
+                scalar refVolume_;
+
+                scalar refPressure_;
+
+               scalar refMassDensity_;
+
+                scalar refNumberDensity_;
+
+
+    // Private Member Functions
+
+        void calcRefValues();
+
+        //- Disallow default bitwise copy construct
+        reducedUnits(const reducedUnits&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const reducedUnits&);
+
+
+public:
+
+    // Static data members
+
+        //- Static data someStaticData
+        static const scalar kb;
+
+
+    // Constructors
+
+        //- Construct with no argument, uses default values:
+        // length  = 1nm
+        // mass = 1.660538782e−27kg (unified atomic mass unit)
+        // temperature = 1K (therefore, energy = 1*kb)
+        reducedUnits();
+
+        //- Construct from components
+        reducedUnits
+        (
+            scalar refLength,
+            scalar refTime,
+            scalar refMass
+        );
+
+        //- Construct from dictionary
+        reducedUnits(const IOdictionary& reducedUnitsDict);
+
+
+    //- Destructor
+    ~reducedUnits();
+
+
+    // Member Functions
+
+        void setRefValues
+        (
+            scalar refLength,
+            scalar refTime,
+            scalar refMass
+        );
+
+        void setRefValues(const IOdictionary& reducedUnitsDict);
+
+
+        // Access
+
+            inline scalar refLength() const;
+
+            inline scalar refTime() const;
+
+            inline scalar refMass() const;
+
+            inline scalar refTemp() const;
+
+            inline scalar refEnergy() const;
+
+            inline scalar refForce() const;
+
+            inline scalar refVelocity() const;
+
+            inline scalar refVolume() const;
+
+            inline scalar refPressure() const;
+
+            inline scalar refMassDensity() const;
+
+            inline scalar refNumberDensity() const;
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<<(Ostream&, const reducedUnits&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "reducedUnitsI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnitsI.H b/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnitsI.H
new file mode 100644
index 0000000000..ec5fa23329
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnitsI.H
@@ -0,0 +1,94 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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  * * * * * * * * * * * * * //
+
+inline Foam::scalar Foam::reducedUnits::refLength() const
+{
+    return refLength_;
+}
+
+
+inline Foam::scalar Foam::reducedUnits::refTime() const
+{
+    return refTime_;
+}
+
+
+inline Foam::scalar Foam::reducedUnits::refMass() const
+{
+    return refMass_;
+}
+
+
+inline Foam::scalar Foam::reducedUnits::refTemp() const
+{
+    return refTemp_;
+}
+
+
+inline Foam::scalar Foam::reducedUnits::refEnergy() const
+{
+    return refEnergy_;
+}
+
+
+inline Foam::scalar Foam::reducedUnits::refForce() const
+{
+    return refForce_;
+}
+
+
+inline Foam::scalar Foam::reducedUnits::refVelocity() const
+{
+    return refVelocity_;
+}
+
+
+inline Foam::scalar Foam::reducedUnits::refVolume() const
+{
+    return refVolume_;
+}
+
+
+inline Foam::scalar Foam::reducedUnits::refPressure() const
+{
+    return refPressure_;
+}
+
+
+inline Foam::scalar Foam::reducedUnits::refMassDensity() const
+{
+    return refMassDensity_;
+}
+
+
+inline Foam::scalar Foam::reducedUnits::refNumberDensity() const
+{
+    return refNumberDensity_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnitsIO.C b/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnitsIO.C
new file mode 100644
index 0000000000..6c7bd8f6f5
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnitsIO.C
@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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 "reducedUnits.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const reducedUnits& rU)
+{
+    os  << nl << "Defined: " << nl
+        << tab << "refLength = " << rU.refLength() << " m" << nl
+        << tab << "refTime = " << rU.refTime() << " s" << nl
+        << tab << "refMass = " << rU.refMass() << " kg" << nl
+        << tab << "Boltzmann constant, kb = " << reducedUnits::kb << " J/K"
+        << nl << "Calculated: " << nl
+        << tab << "refEnergy = " << rU.refEnergy() << " J" << nl
+        << tab << "refTemp = " << rU.refTemp() << " K" << nl
+        << tab << "refForce = " << rU.refForce() << " N" << nl
+        << tab << "refVelocity = " << rU.refVelocity() << " m/s" << nl
+        << tab << "refVolume = " << rU.refVolume() << " m^3" << nl
+        << tab << "refPressure = " << rU.refPressure() << " N/m^2" << nl
+        << tab << "refMassDensity = " << rU.refMassDensity() << " kg/m^3" << nl
+        << tab << "refNumberDensity = " << rU.refNumberDensity() << " m^-3"
+        << endl;
+
+    // Check state of Ostream
+    os.check
+    (
+        "Foam::Ostream& Foam::operator<<(Foam::Ostream&, "
+        "const Foam::reducedUnits&)"
+    );
+
+    return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/electrostaticPotential/electrostaticPotential.C b/src/lagrangian/molecularDynamics/potential/electrostaticPotential/electrostaticPotential.C
index 6589a0ab74..bbf72146e1 100644
--- a/src/lagrangian/molecularDynamics/potential/electrostaticPotential/electrostaticPotential.C
+++ b/src/lagrangian/molecularDynamics/potential/electrostaticPotential/electrostaticPotential.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/electrostaticPotential/electrostaticPotential.H b/src/lagrangian/molecularDynamics/potential/electrostaticPotential/electrostaticPotential.H
index 6df5c1d184..126a42d350 100644
--- a/src/lagrangian/molecularDynamics/potential/electrostaticPotential/electrostaticPotential.H
+++ b/src/lagrangian/molecularDynamics/potential/electrostaticPotential/electrostaticPotential.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunction.C b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunction.C
index 394b82953c..83777d3a94 100644
--- a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunction.C
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunction.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunction.H b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunction.H
index 3fc0746363..29f8cc833b 100644
--- a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunction.H
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunction.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunctionNew.C b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunctionNew.C
index acd844e24a..9c9b669ef5 100644
--- a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunctionNew.C
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunctionNew.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.C b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.C
index cc0af81625..3fb9f9226b 100644
--- a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.C
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.H b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.H
index d9d620db73..c77c07a793 100644
--- a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.H
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.C b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.C
index 14debf4963..057f32f686 100644
--- a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.C
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.H b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.H
index 7347c092a8..c266084c97 100644
--- a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.H
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.C b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.C
index b2f32591f5..319d18392b 100644
--- a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.C
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.H b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.H
index d8ad9e0cd0..e49b0fce3f 100644
--- a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.H
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shiftedForce/shiftedForce.C b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shiftedForce/shiftedForce.C
index 343372d8ec..a4bd4a5291 100644
--- a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shiftedForce/shiftedForce.C
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shiftedForce/shiftedForce.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shiftedForce/shiftedForce.H b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shiftedForce/shiftedForce.H
index 052bdce21c..8fcc2be772 100644
--- a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shiftedForce/shiftedForce.H
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shiftedForce/shiftedForce.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/sigmoid/sigmoid.C b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/sigmoid/sigmoid.C
index 43cb75aa99..c0f006e46c 100644
--- a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/sigmoid/sigmoid.C
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/sigmoid/sigmoid.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/sigmoid/sigmoid.H b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/sigmoid/sigmoid.H
index 15d64e7fdf..83868e491e 100644
--- a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/sigmoid/sigmoid.H
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/sigmoid/sigmoid.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.C b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.C
index a18006f34d..2a8bf6bdbe 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.C
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.H b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.H
index 4b25ca4a93..e89736e035 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.H
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotentialI.H b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotentialI.H
index 9d4bf8c4cd..1b1db7c8ca 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotentialI.H
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotentialI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotentialIO.C b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotentialIO.C
index a3876dc2be..3d74e6c310 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotentialIO.C
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotentialIO.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotentialNew.C b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotentialNew.C
index 1703a06018..981eb655ff 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotentialNew.C
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotentialNew.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/azizChen/azizChen.C b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/azizChen/azizChen.C
index f93e3b3756..c5b40128ce 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/azizChen/azizChen.C
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/azizChen/azizChen.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/coulomb/coulomb.C b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/coulomb/coulomb.C
index cd59601ee6..c3effe0639 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/coulomb/coulomb.C
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/coulomb/coulomb.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/coulomb/coulomb.H b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/coulomb/coulomb.H
index bca4726188..08eb5c5596 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/coulomb/coulomb.H
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/coulomb/coulomb.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/dampedCoulomb/dampedCoulomb.C b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/dampedCoulomb/dampedCoulomb.C
index f37910aec2..82d073fc2e 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/dampedCoulomb/dampedCoulomb.C
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/dampedCoulomb/dampedCoulomb.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/dampedCoulomb/dampedCoulomb.H b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/dampedCoulomb/dampedCoulomb.H
index 88fbdc8902..a6e66c1aa2 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/dampedCoulomb/dampedCoulomb.H
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/dampedCoulomb/dampedCoulomb.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/exponentialRepulsion/exponentialRepulsion.C b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/exponentialRepulsion/exponentialRepulsion.C
index bd2eb772e0..5342b0cfad 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/exponentialRepulsion/exponentialRepulsion.C
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/exponentialRepulsion/exponentialRepulsion.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/exponentialRepulsion/exponentialRepulsion.H b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/exponentialRepulsion/exponentialRepulsion.H
index 80b6cec984..fc0c9b0d5c 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/exponentialRepulsion/exponentialRepulsion.H
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/exponentialRepulsion/exponentialRepulsion.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/lennardJones/lennardJones.C b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/lennardJones/lennardJones.C
index 2acb6db9c7..866b8e5bf1 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/lennardJones/lennardJones.C
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/lennardJones/lennardJones.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/lennardJones/lennardJones.H b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/lennardJones/lennardJones.H
index 3df8a25e01..03c85d2184 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/lennardJones/lennardJones.H
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/lennardJones/lennardJones.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/maitlandSmith/maitlandSmith.C b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/maitlandSmith/maitlandSmith.C
index f6f109fbfd..51cc941b77 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/maitlandSmith/maitlandSmith.C
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/maitlandSmith/maitlandSmith.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/noInteraction/noInteraction.C b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/noInteraction/noInteraction.C
index 42eec34d85..1994393d23 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/noInteraction/noInteraction.C
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/noInteraction/noInteraction.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/noInteraction/noInteraction.H b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/noInteraction/noInteraction.H
index ce8a6b3ac9..162a08c2db 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/noInteraction/noInteraction.H
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/noInteraction/noInteraction.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.C b/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.C
index d5dcc184ea..4b4a0ad57c 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.C
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.H b/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.H
index 09c9f0ee4f..db2f6e965c 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.H
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialListI.H b/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialListI.H
index 91cfb897a1..8657e9cfc8 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialListI.H
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialListI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/potential/potential.C b/src/lagrangian/molecularDynamics/potential/potential/potential.C
index 76a624a652..b34b6e38cc 100644
--- a/src/lagrangian/molecularDynamics/potential/potential/potential.C
+++ b/src/lagrangian/molecularDynamics/potential/potential/potential.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -93,10 +93,7 @@ void Foam::potential::setSiteIdList(const dictionary& moleculePropertiesDict)
 }
 
 
-void Foam::potential::potential::readPotentialDict
-(
-    const dictionary& moleculePropertiesDict
-)
+void Foam::potential::potential::readPotentialDict()
 {
     Info<< nl <<  "Reading potential dictionary:" << endl;
 
@@ -114,28 +111,31 @@ void Foam::potential::potential::readPotentialDict
 
     idList_ = List<word>(idListDict.lookup("idList"));
 
-    setSiteIdList(moleculePropertiesDict);
+    setSiteIdList
+    (
+        IOdictionary
+        (
+            IOobject
+            (
+                "moleculeProperties",
+                mesh_.time().constant(),
+                mesh_,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE,
+                false
+            )
+        )
+    );
 
     List<word> pairPotentialSiteIdList
     (
         SubList<word>(siteIdList_, nPairPotIds_)
     );
 
-    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;
+    Info<< nl << "Unique site ids found: " << siteIdList_
+        << nl << "Site Ids requiring a pair potential: "
+        << pairPotentialSiteIdList
+        << endl;
 
     List<word> tetherSiteIdList(0);
 
@@ -229,25 +229,37 @@ void Foam::potential::potential::readPotentialDict
 
     if (potentialDict.found("external"))
     {
-        Info<< nl << "Reading external forces: ";
+        Info<< nl << "Reading external forces:" << endl;
 
         const dictionary& externalDict = potentialDict.subDict("external");
 
         // gravity
         externalDict.readIfPresent("gravity", gravity_);
-
-        Info<< "gravity = " << gravity_ << nl << endl;
     }
+
+    Info<< nl << tab << "gravity = " << gravity_ << endl;
 }
 
 
 void Foam::potential::potential::readMdInitialiseDict
 (
-    const dictionary& mdInitialiseDict,
-    const dictionary& moleculePropertiesDict,
-    dictionary& idListDict
+    const IOdictionary& mdInitialiseDict,
+    IOdictionary& idListDict
 )
 {
+    IOdictionary moleculePropertiesDict
+    (
+        IOobject
+        (
+            "moleculeProperties",
+            mesh_.time().constant(),
+            mesh_,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE,
+            false
+        )
+    );
+
     DynamicList<word> idList;
 
     DynamicList<word> tetherSiteIdList;
@@ -274,12 +286,11 @@ void Foam::potential::potential::readMdInitialiseDict
                 (
                     "potential::readMdInitialiseDict"
                     "("
-                        "const dictionary&, "
-                        "dictionary&"
+                        "const IOdictionary&, "
+                        "IOdictionary&"
                     ")"
                 )   << "Molecule type " << id
-                    << " not found in " << moleculePropertiesDict.name()
-                    << " dictionary." << nl
+                    << " not found in moleculeProperties dictionary." << nl
                     << abort(FatalError);
             }
 
@@ -330,8 +341,8 @@ void Foam::potential::potential::readMdInitialiseDict
                 (
                     "potential::readMdInitialiseDict"
                     "("
-                        "const dictionary&, "
-                        "dictionary&"
+                        "const IOdictionary&, "
+                        "IOdictionary&"
                     ")"
                 )   << "Tether id  " << tetherSiteId
                     << " not found as a site of any molecule in zone." << nl
@@ -353,34 +364,24 @@ void Foam::potential::potential::readMdInitialiseDict
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::potential::potential
-(
-    const polyMesh& mesh,
-    const dictionary& moleculePropertiesDict
-)
+Foam::potential::potential(const polyMesh& mesh)
 :
     mesh_(mesh)
 {
-    readPotentialDict(moleculePropertiesDict);
+    readPotentialDict();
 }
 
 
 Foam::potential::potential
 (
     const polyMesh& mesh,
-    const dictionary& mdInitialiseDict,
-    const dictionary& moleculePropertiesDict,
-    dictionary& idListDict
+    const IOdictionary& mdInitialiseDict,
+    IOdictionary& idListDict
 )
 :
     mesh_(mesh)
 {
-    readMdInitialiseDict
-    (
-        mdInitialiseDict,
-        moleculePropertiesDict,
-        idListDict
-    );
+    readMdInitialiseDict(mdInitialiseDict, idListDict);
 }
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
diff --git a/src/lagrangian/molecularDynamics/potential/potential/potential.H b/src/lagrangian/molecularDynamics/potential/potential/potential.H
index 02e08c65d0..2aaebf25a5 100644
--- a/src/lagrangian/molecularDynamics/potential/potential/potential.H
+++ b/src/lagrangian/molecularDynamics/potential/potential/potential.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -55,49 +55,35 @@ 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(const dictionary& moleculePropertiesDict);
+        void readPotentialDict();
 
-        //-
         void readMdInitialiseDict
         (
-            const dictionary& mdInitialiseDict,
-            const dictionary& moleculePropertiesDict,
-            dictionary& idListDict
+            const IOdictionary& mdInitialiseDict,
+            IOdictionary& idListDict
         );
 
         //- Disallow default bitwise copy construct
@@ -112,19 +98,14 @@ public:
     // Constructors
 
         //- Construct from mesh reference
-        potential
-        (
-            const polyMesh& mesh,
-            const dictionary& moleculePropertiesDict
-        );
+        potential(const polyMesh& mesh);
 
         //- Construct from mdInitialiseDict
         potential
         (
             const polyMesh& mesh,
-            const dictionary& mdInitialiseDict,
-            const dictionary& moleculePropertiesDict,
-            dictionary& idListDict
+            const IOdictionary& mdInitialiseDict,
+            IOdictionary& idListDict
         );
 
 
@@ -136,31 +117,22 @@ 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/lagrangian/molecularDynamics/potential/potential/potentialI.H b/src/lagrangian/molecularDynamics/potential/potential/potentialI.H
index 79fadb00d0..6f2a0f2d06 100644
--- a/src/lagrangian/molecularDynamics/potential/potential/potentialI.H
+++ b/src/lagrangian/molecularDynamics/potential/potential/potentialI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.C b/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.C
index 9594152e64..8222274f20 100644
--- a/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.C
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.H
index a415fce717..78ab0aa927 100644
--- a/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.H
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotentialNew.C b/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotentialNew.C
index 16ed7674af..35cdef11d7 100644
--- a/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotentialNew.C
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotentialNew.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.C b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.C
index 0b5670a7e5..bad1af9a9b 100644
--- a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.C
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.H
index 12d4d44bf4..4e5b78de09 100644
--- a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.H
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.C b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.C
index ab75d048b1..66a2e7feef 100644
--- a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.C
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.H
index ed881333bd..52bb066758 100644
--- a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.H
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.C b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.C
index e56670b93f..ece3aff2aa 100644
--- a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.C
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.H
index 30d96e6691..ec4a8c8c12 100644
--- a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.H
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.C b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.C
index 016a183cf7..4a25f36f97 100644
--- a/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.C
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.H
index c26b3d8673..d862824f77 100644
--- a/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.H
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialListI.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialListI.H
index d541fed24a..e80c1a1dc4 100644
--- a/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialListI.H
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialListI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License