diff --git a/applications/utilities/preProcessing/molConfig/Make/files b/applications/utilities/preProcessing/molConfig/Make/files
new file mode 100755
index 0000000000..6501e68426
--- /dev/null
+++ b/applications/utilities/preProcessing/molConfig/Make/files
@@ -0,0 +1,8 @@
+latticeStructures = latticeStructures
+velocityDistributions = velocityDistributions
+
+createMolecules.C
+molConfig.C
+genMolConfig.C
+
+EXE = $(FOAM_APPBIN)/molConfig
diff --git a/applications/utilities/preProcessing/molConfig/correctVelocities.H b/applications/utilities/preProcessing/molConfig/correctVelocities.H
new file mode 100644
index 0000000000..0c69f15c7c
--- /dev/null
+++ b/applications/utilities/preProcessing/molConfig/correctVelocities.H
@@ -0,0 +1,21 @@
+for (molN = totalMols; molN < totalMols + totalZoneMols; molN++)
+{
+
+    // Remove bulk momentum introduced by random numbers and add
+    // desired bulk velocity
+
+    // For systems with molecules of significantly differing masses, this may
+    // need to be an iterative process or employ a better algorithm for
+    // removing an appropriate share of the excess momentum from each molecule.
+
+    initialVelocities(molN) += bulkVelocity - momentumSum/totalZoneMols/mass;
+}
+
+// momentumSum = vector::zero;
+//
+// for (molN = totalMols; molN < totalMols + totalZoneMols; molN++)
+// {
+//     momentumSum += mass*initialVelocities(molN);
+// }
+//
+// Info << "Check momentum adjustment: " << momentumSum << endl;
diff --git a/applications/utilities/preProcessing/molConfig/createMolecules.C b/applications/utilities/preProcessing/molConfig/createMolecules.C
new file mode 100644
index 0000000000..20318f06d6
--- /dev/null
+++ b/applications/utilities/preProcessing/molConfig/createMolecules.C
@@ -0,0 +1,253 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2005 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "molConfig.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::molConfig::createMolecules()
+{
+    Info<< nl << "Creating molecules from zone specifications\n" << endl;
+
+    DynamicList<vector> initialPositions(0);
+
+    DynamicList<label> initialIds(0);
+
+    DynamicList<scalar> initialMasses(0);
+
+    DynamicList<label> initialCelli(0);
+
+    DynamicList<vector> initialVelocities(0);
+
+    DynamicList<vector> initialAccelerations(0);
+
+    DynamicList<label> initialTethered(0);
+
+    DynamicList<vector> initialTetherPositions(0);
+
+    label totalMols = 0;
+
+    label idAssign;
+
+    Random rand(clock::getTime());
+
+// * * * * * * * * Building the IdList * * * * * * * * * //
+
+//Notes: - each processor will have an identical idList_.
+//       - The order of id's inside the idList_ depends on the order
+//         of subDicts inside the molConigDict.
+
+    Info<< "Building the idList: " ;
+
+    forAll(molConfigDescription_.toc(), cZs)
+    {
+        word subDictName (molConfigDescription_.toc()[cZs]);
+
+        word iD (molConfigDescription_.subDict(subDictName).lookup("id"));
+
+        if (findIndex(idList_,iD) == -1)
+        {
+            idList_.append(iD);
+        }
+    }
+
+    forAll(idList_, i)
+    {
+        Info << " " << idList_[i];
+    }
+
+    Info << nl << endl;
+
+// * * * * * * * * Filling the Mesh * * * * * * * * * //
+
+    const cellZoneMesh& cellZoneI = mesh_.cellZones();
+
+    if (cellZoneI.size())
+    {
+        Info<< "Filling the zones with molecules." << nl << endl;
+    }
+    else
+    {
+        FatalErrorIn("void createMolecules()\n")
+            << "No cellZones found in mesh description."
+            << abort(FatalError);
+    }
+
+    forAll (cellZoneI, cZ)
+    {
+        if (cellZoneI[cZ].size())
+        {
+            if (!molConfigDescription_.found(cellZoneI[cZ].name()))
+            {
+                Info << "Zone specification subDictionary: "
+                    << cellZoneI[cZ].name() << " not found." << endl;
+            }
+            else
+            {
+                label n = 0;
+
+                label totalZoneMols = 0;
+
+                label molsPlacedThisIteration;
+
+#               include "readZoneSubDict.H"
+
+                idAssign = findIndex(idList_,id);
+
+#               include "startingPoint.H"
+
+                // 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 - it will cause an infinite loop if no molecules
+                // are ever placed by the algorithm.
+
+                if (latticeStructure != "empty")
+                {
+
+                    while
+                    (
+                        molsPlacedThisIteration != 0
+                     || totalZoneMols == 0
+                    )
+                    {
+                        molsPlacedThisIteration = 0;
+
+                        bool partOfLayerInBounds = false;
+
+#                       include "createPositions.H"
+
+                        if
+                        (
+                            totalZoneMols == 0
+                         && !partOfLayerInBounds
+                        )
+                        {
+                            WarningIn("molConfig::createMolecules()")
+                                << "A whole layer of unit cells was placed "
+                                << "outside the bounds of the mesh, but no "
+                                << "molecules have been placed in zone '"
+                                << cellZoneI[cZ].name()
+                                << "'.  This is likely to be because the zone "
+                                << "has few cells ("
+                                << cellZoneI[cZ].size()
+                                << " in this case) and no lattice position "
+                                << "fell inside them.  "
+                                << "Aborting filling this zone."
+                                << endl;
+
+                            break;
+                        }
+
+                        totalZoneMols += molsPlacedThisIteration;
+
+                        n++;
+                    }
+
+                    label molN;
+
+                    for
+                    (
+                        molN = totalMols;
+                        molN < totalMols + totalZoneMols;
+                        molN++
+                    )
+                    {
+                        initialIds.append(idAssign);
+
+                        initialMasses.append(mass);
+
+                        initialAccelerations.append(vector::zero);
+
+                        if (tethered)
+                        {
+                            initialTethered.append(1);
+
+                            initialTetherPositions.append
+                            (
+                                initialPositions[molN]
+                            );
+                        }
+
+                        else
+                        {
+                            initialTethered.append(0);
+
+                            initialTetherPositions.append(vector::zero);
+                        }
+                    }
+
+#                   include "createVelocities.H"
+
+#                   include "correctVelocities.H"
+
+                }
+
+                totalMols += totalZoneMols;
+            }
+        }
+    }
+
+    idList_.shrink();
+
+    positions_ = initialPositions;
+
+    positions_.setSize(initialPositions.size());
+
+    id_ = initialIds;
+
+    id_.setSize(initialIds.size());
+
+    mass_ = initialMasses;
+
+    mass_.setSize(initialMasses.size());
+
+    cells_ = initialCelli;
+
+    cells_.setSize(initialCelli.size());
+
+    U_ = initialVelocities;
+
+    U_.setSize(initialVelocities.size());
+
+    A_ = initialAccelerations;
+
+    A_.setSize(initialAccelerations.size());
+
+    tethered_ = initialTethered;
+
+    tethered_.setSize(initialTethered.size());
+
+    tetherPositions_ = initialTetherPositions;
+
+    tetherPositions_.setSize(initialTetherPositions.size());
+
+    nMol_ = totalMols;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/molConfig/createPositions.H b/applications/utilities/preProcessing/molConfig/createPositions.H
new file mode 100644
index 0000000000..b26486338c
--- /dev/null
+++ b/applications/utilities/preProcessing/molConfig/createPositions.H
@@ -0,0 +1,26 @@
+vector latticePosition;
+
+vector globalPosition;
+
+if (latticeStructure == "SC")
+{
+#   include "SC.H"
+}
+
+else if (latticeStructure == "FCC")
+{
+#   include "FCC.H"
+}
+
+else if (latticeStructure == "BCC")
+{
+#   include "BCC.H"
+}
+
+else
+{
+    FatalErrorIn("createPositions.H\n")
+        << "latticeStructure " << latticeStructure
+        << " not supported."
+        << abort(FatalError);
+}
diff --git a/applications/utilities/preProcessing/molConfig/createVelocities.H b/applications/utilities/preProcessing/molConfig/createVelocities.H
new file mode 100644
index 0000000000..d014ffd146
--- /dev/null
+++ b/applications/utilities/preProcessing/molConfig/createVelocities.H
@@ -0,0 +1,13 @@
+vector velocity;
+
+vector momentumSum = vector::zero;
+
+if (velocityDistribution == "uniform")
+{
+#    include "uniform.H"
+}
+
+if (velocityDistribution == "maxwellian")
+{
+#    include "maxwellian.H"
+}
diff --git a/applications/utilities/preProcessing/molConfig/genMolConfig.C b/applications/utilities/preProcessing/molConfig/genMolConfig.C
new file mode 100644
index 0000000000..91bc5c8d24
--- /dev/null
+++ b/applications/utilities/preProcessing/molConfig/genMolConfig.C
@@ -0,0 +1,129 @@
+/*---------------------------------------------------------------------------*\
+ =========                   |
+ \\      /   F ield          | OpenFOAM: The Open Source CFD Toolbox
+  \\    /    O peration      |
+   \\  /     A nd            | Copyright (C) 1991-2005 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "molConfig.H"
+#include "fvCFD.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+
+    Info<< nl << "Reading molecular configuration description dictionary"
+        << endl;
+
+    IOobject molConfigDescriptionIOobject
+    (
+        "molConfigDict",
+        runTime.system(),
+        runTime,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE,
+        false
+    );
+
+    if (!molConfigDescriptionIOobject.headerOk())
+    {
+        FatalErrorIn(args.executable())
+            << "Cannot find molConfig description file " << nl
+            << args.caseName()/runTime.system()/"molConfig"/"molConfigDict"
+            << nl << exit(FatalError);
+    }
+
+    IOdictionary molConfigDescription(molConfigDescriptionIOobject);
+
+
+    // Create molCloud, registering object with mesh
+
+    Info<< nl << "Creating molecular configuration" << endl;
+
+    molConfig molecules(molConfigDescription, mesh);
+
+    label totalMolecules = molecules.nMol();
+
+    if (Pstream::parRun())
+    {
+        reduce(totalMolecules, sumOp<label>());
+    }
+
+    Info<< nl << "Total number of molecules added: " << totalMolecules
+        << nl << endl;
+
+    moleculeCloud molCloud
+    (
+        mesh,
+        molecules.nMol(),
+        molecules.id(),
+        molecules.mass(),
+        molecules.positions(),
+        molecules.cells(),
+        molecules.U(),
+        molecules.A(),
+        molecules.tethered(),
+        molecules.tetherPositions()
+    );
+
+    IOdictionary idListDict
+    (
+        IOobject
+        (
+            "idList",
+            mesh.time().constant(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        )
+    );
+
+    idListDict.add("idList", molecules.molIdList());
+
+    IOstream::defaultPrecision(12);
+
+    Info << nl << "Writing molecular configuration" << endl;
+
+    if (!mesh.write())
+    {
+        FatalErrorIn(args.executable())
+            << "Failed writing moleculeCloud."
+            << nl << exit(FatalError);
+    }
+
+    Info<< nl << "ClockTime = " << runTime.elapsedClockTime() << " s"
+        << nl << endl;
+
+    Info << nl << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/molConfig/latticeStructures/BCC.H b/applications/utilities/preProcessing/molConfig/latticeStructures/BCC.H
new file mode 100644
index 0000000000..cd3cbdf100
--- /dev/null
+++ b/applications/utilities/preProcessing/molConfig/latticeStructures/BCC.H
@@ -0,0 +1,179 @@
+labelVector iN(0,0,0);
+
+vector gap = (vector::one)*pow((numberDensity/2.0),-(1.0/3.0));
+
+#include "origin.H"
+
+// Info<< "gap = " << gap << endl;
+
+// Special treatment is required for the first position, i.e. iteration zero.
+
+if (n == 0)
+{
+    latticePosition.x() = (iN.x()*gap.x());
+
+    latticePosition.y() = (iN.y()*gap.y());
+
+    latticePosition.z() = (iN.z()*gap.z());
+
+    // Placing 2 molecules in each unit cell, using the algorithm from
+    // D. Rapaport, The Art of Molecular Dynamics Simulation, 2nd Ed, p68
+
+    for (label iU = 0; iU < 2; iU++)
+    {
+        vector unitCellLatticePosition = latticePosition;
+
+        if (iU == 1)
+        {
+            unitCellLatticePosition += 0.5 * gap;
+        }
+
+        if (originSpecifies == "corner")
+        {
+            unitCellLatticePosition -= 0.25*gap;
+        }
+
+//         Info << nl << n << ", " << unitCellLatticePosition;
+
+        globalPosition =
+            origin + transform(latticeToGlobal,unitCellLatticePosition);
+
+        partOfLayerInBounds = mesh_.bounds().contains(globalPosition);
+
+        if
+        (
+             findIndex(mesh_.cellZones()[cZ], mesh_.findCell(globalPosition))
+           != -1
+        )
+        {
+            molsPlacedThisIteration++;
+
+            initialPositions.append(globalPosition);
+
+            initialCelli.append(mesh_.findCell(globalPosition));
+        }
+    }
+}
+else
+{
+    // Place top and bottom caps.
+
+    for (iN.z() = -n; iN.z() <= n; iN.z() += 2*n)
+    {
+        for (iN.y() = -n; iN.y() <= n; iN.y()++)
+        {
+            for (iN.x() = -n; iN.x() <= n; iN.x()++)
+            {
+                latticePosition.x() = (iN.x() * gap.x());
+
+                latticePosition.y() = (iN.y() * gap.y());
+
+                latticePosition.z() = (iN.z() * gap.z());
+
+                for (label iU = 0; iU < 2; iU++)
+                {
+                    vector unitCellLatticePosition = latticePosition;
+
+                    if (iU == 1)
+                    {
+                        unitCellLatticePosition += 0.5*gap;
+                    }
+
+                    if(originSpecifies == "corner")
+                    {
+                        unitCellLatticePosition -= 0.25*gap;
+                    }
+
+//                     Info << nl << iN << ", " << unitCellLatticePosition;
+
+                    globalPosition =
+                        origin
+                      + transform(latticeToGlobal,unitCellLatticePosition);
+
+                    partOfLayerInBounds =
+                        mesh_.bounds().contains(globalPosition);
+
+                    if
+                    (
+                        findIndex
+                        (
+                             mesh_.cellZones()[cZ],
+                             mesh_.findCell(globalPosition)
+                        )
+                     != -1)
+                    {
+                        molsPlacedThisIteration++;
+
+                        initialPositions.append(globalPosition);
+
+                        initialCelli.append(mesh_.findCell(globalPosition));
+                    }
+                }
+            }
+        }
+    }
+
+//     Placing sides
+
+    for (iN.z() = -(n-1); iN.z() <= (n-1); iN.z()++)
+    {
+        for (label iR = 0; iR <= 2*n -1; iR++)
+        {
+            latticePosition.x() = (n*gap.x());
+
+            latticePosition.y() = ((-n + (iR + 1))*gap.y());
+
+            latticePosition.z() = (iN.z() * gap.z());
+
+            for (label iK = 0; iK < 4; iK++)
+            {
+                for (label iU = 0; iU < 2; iU++)
+                {
+                    vector unitCellLatticePosition = latticePosition;
+
+                    if (iU == 1)
+                    {
+                        unitCellLatticePosition += 0.5 * gap;
+                    }
+
+                    if (originSpecifies == "corner")
+                    {
+                        unitCellLatticePosition -= 0.25*gap;
+                    }
+
+                    globalPosition =
+                        origin
+                      + transform(latticeToGlobal,unitCellLatticePosition);
+
+                    partOfLayerInBounds =
+                        mesh_.bounds().contains(globalPosition);
+
+                    if
+                    (
+                        findIndex
+                         (
+                             mesh_.cellZones()[cZ],
+                             mesh_.findCell(globalPosition)
+                         )
+                      != -1
+                    )
+                    {
+                        molsPlacedThisIteration++;
+
+                        initialPositions.append(globalPosition);
+
+                        initialCelli.append(mesh_.findCell(globalPosition));
+                    }
+                }
+
+                latticePosition =
+                    vector
+                    (
+                      - latticePosition.y(),
+                        latticePosition.x(),
+                        latticePosition.z()
+                    );
+            }
+        }
+    }
+}
diff --git a/applications/utilities/preProcessing/molConfig/latticeStructures/FCC.H b/applications/utilities/preProcessing/molConfig/latticeStructures/FCC.H
new file mode 100644
index 0000000000..313ebf96cb
--- /dev/null
+++ b/applications/utilities/preProcessing/molConfig/latticeStructures/FCC.H
@@ -0,0 +1,217 @@
+labelVector iN(0,0,0);
+
+vector gap = (vector::one)*pow((numberDensity/4.0),-(1.0/3.0));
+
+#include "origin.H"
+
+// Info<< "gap = " << gap << endl;
+
+// Special treatment is required for the first position, i.e. iteration zero.
+
+if (n == 0)
+{
+    latticePosition.x() = (iN.x() * gap.x());
+
+    latticePosition.y() = (iN.y() * gap.y());
+
+    latticePosition.z() = (iN.z() * gap.z());
+
+    // Placing 4 molecules in each unit cell, using the algorithm from
+    // D. Rapaport, The Art of Molecular Dynamics Simulation, 2nd Ed, p68
+
+    for (label iU = 0; iU < 4; iU++)
+    {
+        vector unitCellLatticePosition = latticePosition;
+
+        if (iU != 3)
+        {
+            if (iU != 0)
+            {
+                unitCellLatticePosition.x() += 0.5 * gap.x();
+            }
+
+            if (iU != 1)
+            {
+                unitCellLatticePosition.y() += 0.5 * gap.y();
+            }
+
+            if (iU != 2)
+            {
+                unitCellLatticePosition.z() += 0.5 * gap.z();
+            }
+        }
+
+        if (originSpecifies == "corner")
+        {
+            unitCellLatticePosition -= 0.25*gap;
+        }
+
+//         Info << nl << n << ", " << unitCellLatticePosition;
+
+        globalPosition =
+            origin + transform(latticeToGlobal,unitCellLatticePosition);
+
+        partOfLayerInBounds = mesh_.bounds().contains(globalPosition);
+
+        if
+        (
+            findIndex(mesh_.cellZones()[cZ], mesh_.findCell(globalPosition))
+          != -1
+        )
+        {
+            molsPlacedThisIteration++;
+
+            initialPositions.append(globalPosition);
+
+            initialCelli.append(mesh_.findCell(globalPosition));
+        }
+    }
+}
+else
+{
+    // Place top and bottom caps.
+
+    for (iN.z() = -n; iN.z() <= n; iN.z() += 2*n)
+    {
+        for (iN.y() = -n; iN.y() <= n; iN.y()++)
+        {
+            for (iN.x() = -n; iN.x() <= n; iN.x()++)
+            {
+                latticePosition.x() = (iN.x() * gap.x());
+
+                latticePosition.y() = (iN.y() * gap.y());
+
+                latticePosition.z() = (iN.z() * gap.z());
+
+                for (label iU = 0; iU < 4; iU++)
+                {
+                    vector unitCellLatticePosition = latticePosition;
+
+                    if (iU != 3)
+                    {
+                        if (iU != 0)
+                        {
+                            unitCellLatticePosition.x() += 0.5 * gap.x();
+                        }
+
+                        if (iU != 1)
+                        {
+                            unitCellLatticePosition.y() += 0.5 * gap.y();
+                        }
+
+                        if (iU != 2)
+                        {
+                            unitCellLatticePosition.z() += 0.5 * gap.z();
+                        }
+                    }
+
+                    if (originSpecifies == "corner")
+                    {
+                        unitCellLatticePosition -= 0.25*gap;
+                    }
+
+                    globalPosition =
+                        origin
+                      + transform(latticeToGlobal,unitCellLatticePosition);
+
+                    partOfLayerInBounds =
+                        mesh_.bounds().contains(globalPosition);
+
+                    if
+                    (
+                        findIndex
+                        (
+                            mesh_.cellZones()[cZ],
+                            mesh_.findCell(globalPosition)
+                        )
+                     != -1
+                    )
+                    {
+                        molsPlacedThisIteration++;
+
+                        initialPositions.append(globalPosition);
+
+                        initialCelli.append(mesh_.findCell(globalPosition));
+                    }
+                }
+            }
+        }
+    }
+
+//     Placing sides
+
+    for (iN.z() = -(n-1); iN.z() <= (n-1); iN.z()++)
+    {
+        for (label iR = 0; iR <= 2*n -1; iR++)
+        {
+            latticePosition.x() = (n * gap.x());
+
+            latticePosition.y() = ((-n + (iR + 1)) * gap.y());
+
+            latticePosition.z() = (iN.z() * gap.z());
+
+            for (label iK = 0; iK < 4; iK++)
+            {
+                for (label iU = 0; iU < 4; iU++)
+                {
+                    vector unitCellLatticePosition = latticePosition;
+
+                    if (iU != 3)
+                    {
+                        if (iU != 0)
+                        {
+                            unitCellLatticePosition.x() += 0.5 * gap.x();
+                        }
+
+                        if (iU != 1)
+                        {
+                            unitCellLatticePosition.y() += 0.5 * gap.y();
+                        }
+
+                        if (iU != 2)
+                        {
+                            unitCellLatticePosition.z() += 0.5 * gap.z();
+                        }
+                    }
+
+                    if (originSpecifies == "corner")
+                    {
+                        unitCellLatticePosition -= 0.25*gap;
+                    }
+
+                    globalPosition =
+                        origin
+                      + transform(latticeToGlobal,unitCellLatticePosition);
+
+                    partOfLayerInBounds =
+                        mesh_.bounds().contains(globalPosition);
+
+                    if
+                    (
+                        findIndex
+                        (
+                            mesh_.cellZones()[cZ],
+                            mesh_.findCell(globalPosition)
+                        )
+                     != -1
+                    )
+                    {
+                        molsPlacedThisIteration++;
+
+                        initialPositions.append(globalPosition);
+
+                        initialCelli.append(mesh_.findCell(globalPosition));
+                    }
+                }
+
+                latticePosition =
+                    vector
+                    (
+                      - latticePosition.y(),
+                        latticePosition.x(),
+                        latticePosition.z()
+                    );
+            }
+        }
+    }
+}
diff --git a/applications/utilities/preProcessing/molConfig/latticeStructures/SC.H b/applications/utilities/preProcessing/molConfig/latticeStructures/SC.H
new file mode 100644
index 0000000000..dd0ff7c00f
--- /dev/null
+++ b/applications/utilities/preProcessing/molConfig/latticeStructures/SC.H
@@ -0,0 +1,127 @@
+labelVector iN(0,0,0);
+
+vector gap = (vector::one)*pow(numberDensity, -(1.0/3.0));
+
+#include "origin.H"
+
+// Info<< "gap = " << gap << endl;
+
+// Special treatment is required for the first position, i.e. iteration zero.
+
+if (n == 0)
+{
+    latticePosition = vector::zero;
+
+    if (originSpecifies == "corner")
+    {
+        latticePosition += 0.5*gap;
+    }
+
+    globalPosition = origin + transform(latticeToGlobal,latticePosition);
+
+    partOfLayerInBounds = mesh_.bounds().contains(globalPosition);
+
+    if (findIndex(mesh_.cellZones()[cZ], mesh_.findCell(globalPosition)) != -1)
+    {
+        molsPlacedThisIteration++;
+
+        initialPositions.append(globalPosition);
+
+        initialCelli.append(mesh_.findCell(globalPosition));
+    }
+}
+else
+{
+    for (iN.z() = -n; iN.z() <= n; iN.z() += 2*n)
+    {
+        for (iN.y() = -n; iN.y() <= n; iN.y()++)
+        {
+            for (iN.x() = -n; iN.x() <= n; iN.x()++)
+            {
+                latticePosition.x() = (iN.x() * gap.x());
+
+                latticePosition.y() = (iN.y() * gap.y());
+
+                latticePosition.z() = (iN.z() * gap.z());
+
+                if (originSpecifies == "corner")
+                {
+                    latticePosition += 0.5*gap;
+                }
+
+                globalPosition =
+                    origin + transform(latticeToGlobal,latticePosition);
+
+                partOfLayerInBounds = mesh_.bounds().contains(globalPosition);
+
+                if
+                (
+                    findIndex
+                    (
+                        mesh_.cellZones()[cZ],
+                        mesh_.findCell(globalPosition)
+                    )
+                 != -1
+                )
+                {
+                    molsPlacedThisIteration++;
+
+                    initialPositions.append(globalPosition);
+
+                    initialCelli.append(mesh_.findCell(globalPosition));
+                }
+            }
+        }
+    }
+
+    tensor quarterRotate(EulerCoordinateRotation(-90, 0, 0, true).R());
+
+    iN.x() = n;
+    for (iN.z() = -(n-1); iN.z() <= (n-1); iN.z()++)
+    {
+        for (iN.y() = -(n-1); iN.y() <= n; iN.y()++)
+        {
+            latticePosition.x() = (iN.x()*gap.x());
+
+            latticePosition.y() = (iN.y()*gap.y());
+
+            latticePosition.z() = (iN.z()*gap.z());
+
+            for (label iR = 0; iR < 4; iR++)
+            {
+                vector offsetCorrectedLatticePosition = latticePosition;
+
+                if (originSpecifies == "corner")
+                {
+                    offsetCorrectedLatticePosition += 0.5*gap;
+                }
+
+                globalPosition =
+                    origin
+                  + transform(latticeToGlobal,offsetCorrectedLatticePosition);
+
+                partOfLayerInBounds = mesh_.bounds().contains(globalPosition);
+
+                if
+                (
+                    findIndex
+                    (
+                        mesh_.cellZones()[cZ],
+                        mesh_.findCell(globalPosition)
+                    )
+                 != -1
+                )
+                {
+                    molsPlacedThisIteration++;
+
+                    initialPositions.append(globalPosition);
+
+                    initialCelli.append(mesh_.findCell(globalPosition));
+                }
+
+                latticePosition = transform(quarterRotate,latticePosition);
+            }
+        }
+    }
+}
+
diff --git a/applications/utilities/preProcessing/molConfig/molConfig.C b/applications/utilities/preProcessing/molConfig/molConfig.C
new file mode 100644
index 0000000000..e837346b97
--- /dev/null
+++ b/applications/utilities/preProcessing/molConfig/molConfig.C
@@ -0,0 +1,50 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2005 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "molConfig.H"
+
+// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
+
+Foam::molConfig::molConfig
+(
+    IOdictionary& molConfigDescription,
+    const polyMesh& mesh
+)
+:
+    molConfigDescription_(molConfigDescription),
+    mesh_(mesh)
+{
+    createMolecules();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::molConfig::~molConfig()
+{}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/molConfig/molConfig.H b/applications/utilities/preProcessing/molConfig/molConfig.H
new file mode 100644
index 0000000000..a9d549da12
--- /dev/null
+++ b/applications/utilities/preProcessing/molConfig/molConfig.H
@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2005 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::molConfig
+
+Description
+
+SourceFiles
+    molConfigI.H
+    molConfig.C
+    molConfigIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef molConfig_H
+#define molConfig_H
+
+#include "labelVector.H"
+#include "scalar.H"
+#include "vector.H"
+#include "labelField.H"
+#include "scalarField.H"
+#include "vectorField.H"
+#include "IOField.H"
+#include "EulerCoordinateRotation.H"
+#include "Random.H"
+
+#include "Time.H"
+#include "IOdictionary.H"
+#include "IOstreams.H"
+#include "moleculeCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class molConfig Declaration
+\*---------------------------------------------------------------------------*/
+
+class molConfig
+{
+    // Private data
+
+        const IOdictionary& molConfigDescription_;
+
+        const polyMesh& mesh_;
+
+        DynamicList<word> idList_;
+
+        labelField id_;
+
+        scalarField mass_;
+
+        vectorField positions_;
+
+        labelField cells_;
+
+        vectorField U_;
+
+        vectorField A_;
+
+        labelField tethered_;
+
+        vectorField tetherPositions_;
+
+        label nMol_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from IOdictionary and mesh
+        molConfig(IOdictionary&, const polyMesh&);
+
+
+    // Destructor
+
+        ~molConfig();
+
+
+    // Member Functions
+
+        void createMolecules();
+
+
+        // Access
+
+            inline const List<word>& molIdList() const;
+
+            inline const labelField& id() const;
+
+            inline const scalarField& mass() const;
+
+            inline const vectorField& positions() const;
+
+            inline const labelField& cells() const;
+
+            inline const vectorField& U() const;
+
+            inline const vectorField& A() const;
+
+            inline const labelField& tethered() const;
+
+            inline const vectorField& tetherPositions() const;
+
+            inline label nMol() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "molConfigI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/molConfig/molConfigI.H b/applications/utilities/preProcessing/molConfig/molConfigI.H
new file mode 100644
index 0000000000..b85a656c03
--- /dev/null
+++ b/applications/utilities/preProcessing/molConfig/molConfigI.H
@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2005 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const List<word>& molConfig::molIdList() const
+{
+    return idList_;
+}
+
+
+inline const labelField& molConfig::id() const
+{
+    return id_;
+}
+
+
+inline const scalarField& molConfig::mass() const
+{
+    return mass_;
+}
+
+
+inline const vectorField& molConfig::positions() const
+{
+    return positions_;
+}
+
+
+inline const labelField& molConfig::cells() const
+{
+    return cells_;
+}
+
+
+inline const vectorField& molConfig::U() const
+{
+    return U_;
+}
+
+
+inline const vectorField& molConfig::A() const
+{
+    return A_;
+}
+
+
+inline const labelField& molConfig::tethered() const
+{
+    return tethered_;
+}
+
+
+inline const vectorField& molConfig::tetherPositions() const
+{
+    return tetherPositions_;
+}
+
+
+inline label molConfig::nMol() const
+{
+    return nMol_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/molConfig/origin.H b/applications/utilities/preProcessing/molConfig/origin.H
new file mode 100644
index 0000000000..24bc994d34
--- /dev/null
+++ b/applications/utilities/preProcessing/molConfig/origin.H
@@ -0,0 +1,49 @@
+// Please refer to notes
+
+// 1. Determine the unit cell dimensions: xU, yU and zU
+
+const scalar xU = gap.x();
+const scalar yU = gap.y();
+const scalar zU = gap.z();
+
+// 2. Determine the anchorPoint co-ordinates: xA, yA  and zA
+
+const scalar xA = anchorPoint.x();
+const scalar yA = anchorPoint.y();
+const scalar zA = anchorPoint.z();
+
+// 3. Determine the vector rAB from global co-ordinate system:
+
+const vector rAB((xMid - xA), (yMid - yA), (zMid - zA));
+
+// 4. Transform vector rAS into lattice co-ordinate system:
+
+const vector rASTransf = transform(latticeToGlobal.T(), rAB);
+
+// Info << "The vector rAS = " << rAS << endl;
+// Info << "The vector rAStransf = " << rAStransf << endl;
+
+// 5. Calculate the integer values: ni, nj and nk
+scalar nIscalar = rASTransf.x()/xU;
+scalar nJscalar = rASTransf.y()/yU;
+scalar nKscalar = rASTransf.z()/zU;
+
+// Info << "The nI, nJ, nK values before are: " << nIscalar <<" "<< nJscalar <<" "<< nKscalar << endl;
+
+label nI = label(nIscalar + 0.5*sign(nIscalar));
+label nJ = label(nJscalar + 0.5*sign(nJscalar));
+label nK = label(nKscalar + 0.5*sign(nKscalar));
+
+// Info << "The nI, nJ, nK values after are: " << nI <<" "<< nJ <<" "<< nK << endl;
+
+// 6. Calculate the corrected starting point, rAC (in the lattice co-ordinate system):
+const vector rAC((nI*xU), (nJ*yU), (nK*zU));
+
+// 7. Transform the corrected starting point in the global co-ordinate system, rC:
+const vector rC = anchorPoint + transform(latticeToGlobal, rAC);
+
+
+const vector& origin = rC;
+
+// Pout << "The Corrected Starting Point: " << origin << endl;
+
diff --git a/applications/utilities/preProcessing/molConfig/readZoneSubDict.H b/applications/utilities/preProcessing/molConfig/readZoneSubDict.H
new file mode 100644
index 0000000000..72705c761e
--- /dev/null
+++ b/applications/utilities/preProcessing/molConfig/readZoneSubDict.H
@@ -0,0 +1,93 @@
+
+// Info << "Zone description subDict " << cZ <<": "  << cellZoneI[cZ].name() << endl;
+
+const dictionary& subDictI =
+    molConfigDescription_.subDict(cellZoneI[cZ].name());
+
+const scalar temperature(readScalar(subDictI.lookup("temperature")));
+
+const word velocityDistribution(subDictI.lookup("velocityDistribution"));
+
+const vector bulkVelocity(subDictI.lookup("bulkVelocity"));
+
+const word id(subDictI.lookup("id"));
+
+const scalar mass(readScalar(subDictI.lookup("mass")));
+
+scalar numberDensity_read(0.0);
+
+if (subDictI.found("numberDensity"))
+{
+    numberDensity_read = readScalar(subDictI.lookup("numberDensity"));
+}
+else if (subDictI.found("massDensity"))
+{
+    numberDensity_read = readScalar(subDictI.lookup("massDensity"))/mass;
+}
+else
+{
+    FatalErrorIn("readZoneSubDict.H\n")
+        << "massDensity or numberDensity not specified " << nl
+        << abort(FatalError);
+}
+
+const scalar numberDensity(numberDensity_read);
+
+const word latticeStructure(subDictI.lookup("latticeStructure"));
+
+const vector anchorPoint(subDictI.lookup("anchor"));
+
+const word originSpecifies(subDictI.lookup("anchorSpecifies"));
+
+if
+(
+    originSpecifies != "corner"
+ && originSpecifies != "molecule"
+)
+{
+    FatalErrorIn("readZoneSubDict.H\n")
+        << "anchorSpecifies must be either 'corner' or 'molecule', found "
+        << originSpecifies << nl
+        << abort(FatalError);
+}
+
+bool tethered = false;
+
+if (subDictI.found("tethered"))
+{
+    tethered = Switch(subDictI.lookup("tethered"));
+}
+
+const vector orientationAngles(subDictI.lookup("orientationAngles"));
+
+scalar phi(orientationAngles.x()*mathematicalConstant::pi/180.0);
+scalar theta(orientationAngles.y()*mathematicalConstant::pi/180.0);
+scalar psi(orientationAngles.z()*mathematicalConstant::pi/180.0);
+
+const tensor latticeToGlobal
+(
+    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)
+);
+
+// Info << "\tcells: " << cellZoneI[cZ].size() << endl;
+// Info << "\tnumberDensity: " << numberDensity << endl;
+// Info << "\ttemperature: " << temperature << endl;
+// Info << "\tvelocityDistribution: " << velocityDistribution << endl;
+// Info << "\tbulkVelocity: " << bulkVelocity << endl;
+// Info << "\tid: " << id << endl;
+// Info << "\tmass: " << mass << endl;
+// Info << "\tlatticeStructure: " << latticeStructure << endl;
+// Info << "\tanchor: " << anchorPoint << endl;
+// Info << "\toriginSpecifies: " << originSpecifies << endl;
+// Info << "\ttethered: " << tethered << endl;
+// Info << "\torientationAngles: " << orientationAngles << endl;
+// Info << "\tlatticeToGlobal: " << latticeToGlobal << endl;
+
diff --git a/applications/utilities/preProcessing/molConfig/startingPoint.H b/applications/utilities/preProcessing/molConfig/startingPoint.H
new file mode 100644
index 0000000000..6ae01c83f1
--- /dev/null
+++ b/applications/utilities/preProcessing/molConfig/startingPoint.H
@@ -0,0 +1,97 @@
+scalar xMax = 0;
+
+scalar yMax = 0;
+
+scalar zMax = 0;
+
+scalar xMin = 0;
+
+scalar yMin = 0;
+
+scalar zMin = 0;
+
+label xMaxPtLabel = 0;
+
+label yMaxPtLabel = 0;
+
+label zMaxPtLabel = 0;
+
+label xMinPtLabel = 0;
+
+label yMinPtLabel = 0;
+
+label zMinPtLabel = 0;
+
+forAll (cellZoneI[cZ], nC)
+{
+    const labelList& cellPointsJ = mesh_.cellPoints()[cellZoneI[cZ][nC]];
+
+    forAll(cellPointsJ, nP)
+    {
+        const point& ptI = mesh_.points()[cellPointsJ[nP]];
+
+        const label& ptILabel = cellPointsJ[nP];
+
+        if (ptI.x() > xMax || nC == 0)
+        {
+            xMax = ptI.x();
+            xMaxPtLabel = ptILabel;
+        }
+        if (ptI.y() > yMax || nC == 0)
+        {
+            yMax = ptI.y();
+            yMaxPtLabel = ptILabel;
+        }
+        if (ptI.z() > zMax || nC == 0)
+        {
+            zMax = ptI.z();
+            zMaxPtLabel = ptILabel;
+        }
+        if (ptI.x() < xMin || nC == 0)
+        {
+            xMin = ptI.x();
+            xMinPtLabel = ptILabel;
+        }
+        if (ptI.y() < yMin || nC == 0)
+        {
+            yMin = ptI.y();
+            yMinPtLabel = ptILabel;
+        }
+        if (ptI.z() < zMin || nC == 0)
+        {
+            zMin = ptI.z();
+            zMinPtLabel = ptILabel;
+        }
+    }
+}
+
+// Info << "Xmax: label = " << xMaxPtLabel2 << "; vector = " <<mesh_.points()[xMaxPtLabel2]
+//      <<"; x-component = " << mesh_.points()[xMaxPtLabel2].x() << endl;
+// Info << "Ymax: label = " << yMaxPtLabel2 << "; vector = " <<mesh_.points()[yMaxPtLabel2]
+//      <<"; y-component = " << mesh_.points()[yMaxPtLabel2].y() << endl;
+// Info << "Zmax: label = " << zMaxPtLabel2 << "; vector = " <<mesh_.points()[zMaxPtLabel2]
+//      <<"; z-component = " << mesh_.points()[zMaxPtLabel2].z() << endl;
+//
+// Info << "Xmin: label = " << xMinPtLabel << "; vector = " <<mesh_.points()[xMinPtLabel]
+//      <<"; x-component = " << mesh_.points()[xMinPtLabel].x() << endl;
+// Info << "Ymin: label = " << yMinPtLabel << "; vector = " <<mesh_.points()[yMinPtLabel]
+//      <<"; y-component = " << mesh_.points()[yMinPtLabel].y() << endl;
+// Info << "Zmin: label = " << zMinPtLabel << "; vector = " <<mesh_.points()[zMinPtLabel]
+//      <<"; z-component = " << mesh_.points()[zMinPtLabel].z() << endl;
+
+scalar xMid =
+    (mesh_.points()[xMaxPtLabel].x()
+  + mesh_.points()[xMinPtLabel].x()) / 2;
+
+scalar yMid =
+    (mesh_.points()[yMaxPtLabel].y()
+  + mesh_.points()[yMinPtLabel].y()) / 2;
+
+scalar zMid =
+    (mesh_.points()[zMaxPtLabel].z()
+  + mesh_.points()[zMinPtLabel].z()) / 2;
+
+vector rS(xMid, yMid, zMid);
+
+// Info << "\t The Estimated Starting Point: " << rS << endl;
+
diff --git a/applications/utilities/preProcessing/molConfig/velocityDistributions/maxwellian.H b/applications/utilities/preProcessing/molConfig/velocityDistributions/maxwellian.H
new file mode 100644
index 0000000000..01c3ff8716
--- /dev/null
+++ b/applications/utilities/preProcessing/molConfig/velocityDistributions/maxwellian.H
@@ -0,0 +1,26 @@
+scalar velCmptMag = sqrt(moleculeCloud::kb*temperature/mass);
+
+for (molN = totalMols; molN < totalMols + totalZoneMols; molN++)
+{
+    // Assign velocity: random direction, magnitude determined by desired
+    // maxwellian distribution at temperature
+
+    // Temperature gradients could be created by specifying a gradient in the
+    // zone subDict, or by reading a field from a mesh.
+
+    // The velocities are treated on a zone-by-zone basis for the purposes of
+    // removal of bulk momentum - hence nMols becomes totalZoneMols
+
+    velocity = vector
+    (
+        velCmptMag*rand.GaussNormal(),
+        velCmptMag*rand.GaussNormal(),
+        velCmptMag*rand.GaussNormal()
+    );
+
+    momentumSum += mass*velocity;
+
+    initialVelocities.append(velocity);
+}
+
+
diff --git a/applications/utilities/preProcessing/molConfig/velocityDistributions/uniform.H b/applications/utilities/preProcessing/molConfig/velocityDistributions/uniform.H
new file mode 100644
index 0000000000..9c178af262
--- /dev/null
+++ b/applications/utilities/preProcessing/molConfig/velocityDistributions/uniform.H
@@ -0,0 +1,27 @@
+scalar initVelMag =
+    sqrt
+    (
+        3.0*(1.0 - 1.0 / totalZoneMols)
+       *moleculeCloud::kb*temperature
+       /mass
+    );
+
+for (molN = totalMols; molN < totalMols + totalZoneMols; molN++)
+{
+    // Assign velocity: random direction, magnitude determined by desired
+    // temperature
+
+    // Temperature gradients could be created by specifying a gradient in the
+    // zone subDict, or by reading a field from a mesh.
+
+    // The velocities are treated on a zone-by-zone basis for the purposes of
+    // removal of bulk momentum - hence nMols becomes totalZoneMols
+
+    velocity = (2.0*rand.vector01() - vector::one);
+
+    velocity *= initVelMag/mag(velocity);
+
+    momentumSum += mass*velocity;
+
+    initialVelocities.append(velocity);
+}
diff --git a/tutorials/mdEquilibrationFoam/periodicCube/constant/polyMesh/blockMeshDict b/tutorials/mdEquilibrationFoam/periodicCube/constant/polyMesh/blockMeshDict
new file mode 100644
index 0000000000..23a2ee9874
--- /dev/null
+++ b/tutorials/mdEquilibrationFoam/periodicCube/constant/polyMesh/blockMeshDict
@@ -0,0 +1,73 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.3                                   |
+|   \\  /    A nd           | Web:      http://www.openfoam.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+    root            "";
+    case            "";
+    instance        "";
+    local           "";
+
+    class           dictionary;
+    object          blockMeshDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 2.462491658e-9;
+
+vertices
+(
+    (-1 -1 -1)
+    (1 -1 -1)
+    (1 1 -1)
+    (-1 1 -1)
+    (-1 -1 1)
+    (1 -1 1)
+    (1 1 1)
+    (-1 1 1)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) liquid (12 12 12) simpleGrading (1 1 1)
+);
+
+patches
+(
+    cyclic
+    periodicX
+    (
+        (1 2 6 5)
+        (0 4 7 3)
+    )
+
+    cyclic
+    periodicY
+    (
+        (2 3 7 6)
+        (0 1 5 4)
+    )
+
+    cyclic
+    periodicZ
+    (
+        (0 3 2 1)
+        (4 5 6 7)
+    )
+)
+
+mergePatchPairs
+(
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/mdEquilibrationFoam/periodicCube/system/fvSchemes b/tutorials/mdEquilibrationFoam/periodicCube/system/fvSchemes
new file mode 100644
index 0000000000..8504e0fd45
--- /dev/null
+++ b/tutorials/mdEquilibrationFoam/periodicCube/system/fvSchemes
@@ -0,0 +1,67 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.3                                   |
+|   \\  /    A nd           | Web:      http://www.openfoam.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+    root            "";
+    case            "";
+    instance        "";
+    local           "";
+
+    class           dictionary;
+    object          fvSchemes;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(p)         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         none;
+    laplacian(nu,U) Gauss linear corrected;
+    laplacian(1|A(U),p) Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+    interpolate(HbyA) linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mdEquilibrationFoam/periodicCube/system/fvSolution b/tutorials/mdEquilibrationFoam/periodicCube/system/fvSolution
new file mode 100644
index 0000000000..c85c7b791b
--- /dev/null
+++ b/tutorials/mdEquilibrationFoam/periodicCube/system/fvSolution
@@ -0,0 +1,40 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.3                                   |
+|   \\  /    A nd           | Web:      http://www.openfoam.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+    root            "";
+    case            "";
+    instance        "";
+    local           "";
+
+    class           dictionary;
+    object          fvSolution;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p               ICCG 1e-06 0;
+    U               BICCG 1e-05 0;
+}
+
+PISO
+{
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+    pRefCell        0;
+    pRefValue       0;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mdEquilibrationFoam/periodicCube/system/mdEquilibrationDict b/tutorials/mdEquilibrationFoam/periodicCube/system/mdEquilibrationDict
new file mode 100644
index 0000000000..e2323961f5
--- /dev/null
+++ b/tutorials/mdEquilibrationFoam/periodicCube/system/mdEquilibrationDict
@@ -0,0 +1,27 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.4.1                                  |
+|   \\  /    A nd           | Web:      http://www.openfoam.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+    root            "";
+    case            "";
+    instance        "";
+    local           "";
+
+    class           dictionary;
+    object          mdEquilibrationDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+equilibrationTargetTemperature  300.0;
+
+// ************************************************************************* //
\ No newline at end of file
diff --git a/tutorials/mdEquilibrationFoam/periodicCube/system/molConfigDict b/tutorials/mdEquilibrationFoam/periodicCube/system/molConfigDict
new file mode 100644
index 0000000000..b01cd00edb
--- /dev/null
+++ b/tutorials/mdEquilibrationFoam/periodicCube/system/molConfigDict
@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.3                                   |
+|   \\  /    A nd           | Web:      http://www.openfoam.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+    root            "";
+    case            "";
+    instance        "";
+    local           "";
+
+    class           dictionary;
+    object          molConfigDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Euler angles, expressed in degrees as phi, theta, psi,
+// see http://mathworld.wolfram.com/EulerAngles.html
+
+liquid
+{
+    massDensity             1220.0;
+    temperature             300.0;
+    velocityDistribution    maxwellian;
+    bulkVelocity            (0.0 0.0 0.0);
+    id                      Ar;
+    mass                    6.63352033e-26;
+    latticeStructure        SC;
+    anchor                  (0.0 0.0 0.0);
+    anchorSpecifies         molecule;
+    tethered                no;
+    orientationAngles       (0 0 0);
+}
+
+// ************************************************************************* //