Merge branch 'master' of /home/noisy2/OpenFOAM/OpenFOAM-dev/

applications/solvers/molecularDynamics/gnemdFoam/Make/files

@@ -0,0 +1,3 @@
+gnemdFoam.C
+
+EXE = $(FOAM_APPBIN)/gnemdFoam

applications/solvers/molecularDynamics/gnemdFoam/Make/options

@@ -0,0 +1,11 @@
+EXE_INC = \
+    -I$(LIB_SRC)/lagrangian/molecule/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools \
+    -lfiniteVolume \
+    -llagrangian \
+    -lmolecule

applications/solvers/molecularDynamics/gnemdFoam/gnemdFoam.C

@@ -0,0 +1,90 @@
+/*---------------------------------------------------------------------------*\
+ =========                   |
+ \\      /   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
+
+Application
+    gnemdFOAM
+
+Description
+    MD for Fluid Mechanics and hybridising with a continuum solver.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "md.H"
+
+int main(int argc, char *argv[])
+{
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+
+    moleculeCloud molecules(mesh);
+
+#   include "createMDFields.H"
+
+    molecules.removeHighEnergyOverlaps();
+
+#   include "temperatureAndPressureVariables.H"
+
+    label nAveragingSteps = 0;
+
+    Info << "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        runTime++;
+
+        nAveragingSteps++;
+
+        Info << "Time = " << runTime.timeName() << endl;
+
+        molecules.integrateEquationsOfMotion();
+
+#       include "meanMomentumEnergyAndNMols.H"
+
+#       include "temperatureAndPressure.H"
+
+#       include "calculateMDFields.H"
+
+#       include "averageMDFields.H"
+
+        runTime.write();
+
+#       include "resetMDFields.H"
+
+        if (runTime.outputTime())
+        {
+            nAveragingSteps = 0;
+        }
+
+        Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info << "End\n" << endl;
+
+    return(0);
+}

applications/solvers/molecularDynamics/mdEquilibrationFoam/Make/files

@@ -0,0 +1,3 @@
+mdEquilibrationFoam.C
+
+EXE = $(FOAM_APPBIN)/mdEquilibrationFoam

applications/solvers/molecularDynamics/mdEquilibrationFoam/Make/options

@@ -0,0 +1,11 @@
+EXE_INC = \
+    -I$(LIB_SRC)/lagrangian/molecule/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools \
+    -lfiniteVolume \
+    -llagrangian \
+    -lmolecule

applications/solvers/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C

@@ -0,0 +1,86 @@
+/*---------------------------------------------------------------------------*\
+ =========                   |
+ \\      /   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
+
+Application
+    mdEquilibrationFOAM
+
+Description
+    Equilibrates and/or preconditions MD systems
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "md.H"
+
+int main(int argc, char *argv[])
+{
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+
+    moleculeCloud molecules(mesh);
+
+    molecules.removeHighEnergyOverlaps();
+
+#   include "temperatureAndPressureVariables.H"
+
+#   include "readmdEquilibrationDict.H"
+
+    label nAveragingSteps = 0;
+
+    Info << "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        runTime++;
+
+        nAveragingSteps++;
+
+        Info << "Time = " << runTime.timeName() << endl;
+
+        molecules.integrateEquationsOfMotion();
+
+#       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);
+}

applications/solvers/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H

@@ -0,0 +1,18 @@
+Info<< "Reading MD Equilibration Dictionary" << nl << endl;
+
+IOdictionary mdEquilibrationDict
+(
+    IOobject
+    (
+        "mdEquilibrationDict",
+        runTime.system(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+scalar targetTemperature = readScalar
+(
+    mdEquilibrationDict.lookup("equilibrationTargetTemperature")
+);

applications/utilities/postProcessing/field/components/components.C

@@ -35,44 +35,7 @@ Description
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-template <class Type>
+#include "writeComponentFields.C"
-void writeComponents
-(
-    const IOobject& header,
-    const fvMesh& mesh,
-    bool& processed
-)
-{
-    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
-
-    if (header.headerClassName() == fieldType::typeName)
-    {
-        Info<< "    Reading " << header.name() << endl;
-        fieldType field(header, mesh);
-
-        for (direction i=0; i<Type::nComponents; i++)
-        {
-            Info<< "    Calculating " << header.name()
-                << Type::componentNames[i] << endl;
-
-            volScalarField componentField
-            (
-                IOobject
-                (
-                    header.name() + word(Type::componentNames[i]),
-                    mesh.time().timeName(),
-                    mesh,
-                    IOobject::NO_READ
-                ),
-                field.component(i)
-            );
-            componentField.write();
-        }
-
-        processed = true;
-    }
-}
-
 
 int main(int argc, char *argv[])
 {
@@ -115,10 +78,15 @@ int main(int argc, char *argv[])
             mesh.readUpdate();
 
             bool processed = false;
-            writeComponents<vector>(fieldHeader, mesh, processed);
+            writeComponentFields<vector>(fieldHeader, mesh, processed);
-            writeComponents<sphericalTensor>(fieldHeader, mesh, processed);
+            writeComponentFields<sphericalTensor>
-            writeComponents<symmTensor>(fieldHeader, mesh, processed);
+            (
-            writeComponents<tensor>(fieldHeader, mesh, processed);
+                fieldHeader,
+                mesh,
+                processed
+            );
+            writeComponentFields<symmTensor>(fieldHeader, mesh, processed);
+            writeComponentFields<tensor>(fieldHeader, mesh, processed);
             if (!processed)
             {
                 FatalError

applications/utilities/postProcessing/field/components/writeComponentFields.C

@@ -0,0 +1,77 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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
+
+Application
+    components
+
+Description
+    Writes scalar fields corresponding to each component of the supplied
+    field (name) for each time.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template <class Type>
+void writeComponentFields
+(
+    const IOobject& header,
+    const fvMesh& mesh,
+    bool& processed
+)
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    if (header.headerClassName() == fieldType::typeName)
+    {
+        Info<< "    Reading " << header.name() << endl;
+        fieldType field(header, mesh);
+
+        for (direction i=0; i<Type::nComponents; i++)
+        {
+            Info<< "    Calculating " << header.name()
+                << Type::componentNames[i] << endl;
+
+            volScalarField componentField
+            (
+                IOobject
+                (
+                    header.name() + word(Type::componentNames[i]),
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ
+                ),
+                field.component(i)
+            );
+            componentField.write();
+        }
+
+        processed = true;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/utilities/postProcessing/field/mag/mag.C

@@ -34,39 +34,7 @@ Description
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-template<class Type>
+#include "writeMagField.C"
-void writeMagField
-(
-    const IOobject& header,
-    const fvMesh& mesh,
-    bool& processed
-)
-{
-    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
-
-    if (header.headerClassName() == fieldType::typeName)
-    {
-        Info<< "    Reading " << header.name() << endl;
-        fieldType field(header, mesh);
-
-        Info<< "    Calculating mag" << header.name() << endl;
-        volScalarField magField
-        (
-            IOobject
-            (
-                "mag" + header.name(),
-                mesh.time().timeName(),
-                mesh,
-                IOobject::NO_READ
-            ),
-            mag(field)
-        );
-        magField.write();
-
-        processed = true;
-    }
-}
-
 
 int main(int argc, char *argv[])
 {

applications/utilities/postProcessing/field/mag/writeMagField.C

@@ -0,0 +1,71 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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
+
+Application
+    mag
+
+Description
+    Calculates and writes the magnitude of a field for each time
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+void writeMagField
+(
+    const IOobject& header,
+    const fvMesh& mesh,
+    bool& processed
+)
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    if (header.headerClassName() == fieldType::typeName)
+    {
+        Info<< "    Reading " << header.name() << endl;
+        fieldType field(header, mesh);
+
+        Info<< "    Calculating mag" << header.name() << endl;
+        volScalarField magField
+        (
+            IOobject
+            (
+                "mag" + header.name(),
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ
+            ),
+            mag(field)
+        );
+        magField.write();
+
+        processed = true;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/utilities/postProcessing/field/magGrad/magGrad.C

@@ -26,23 +26,26 @@ Application
     magGrad
 
 Description
-    Calculates and writes the scalar magnitude of a scalar or vector field
+    Calculates and writes the magnitude of the gradient of a field for each
-    at each time
+    time
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-//  Main program:
+
+#include "writeMagGradField.C"
 
 int main(int argc, char *argv[])
 {
-    argList::validArgs.append("field");
+    argList::validArgs.append("fieldName");
 
 #   include "addTimeOptions.H"
 #   include "setRootCase.H"
 
+    word fieldName(args.additionalArgs()[0]);
+
 #   include "createTime.H"
 
     // Get times list
@@ -51,8 +54,6 @@ int main(int argc, char *argv[])
     // set startTime and endTime depending on -time and -latestTime options
 #   include "checkTimeOptions.H"
 
-    const word fieldName(args.additionalArgs()[0]);
-
     runTime.setTime(Times[startTime], startTime);
 
 #   include "createMesh.H"
@@ -63,7 +64,7 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << endl;
 
-        IOobject header
+        IOobject fieldHeader
         (
             fieldName,
             runTime.timeName(),
@@ -71,56 +72,33 @@ int main(int argc, char *argv[])
             IOobject::MUST_READ
         );
 
-        // Check U exists
+        // Check field "fieldName" exists
-        if (header.headerOk())
+        if (fieldHeader.headerOk())
         {
             mesh.readUpdate();
 
-            if (header.headerClassName() == volVectorField::typeName)
+            bool processed = false;
+            writeMagGradField<scalar>(fieldHeader, mesh, processed);
+            writeMagGradField<vector>(fieldHeader, mesh, processed);
+            if (!processed)
             {
-                Info<< "    Reading " << fieldName << endl;
+                FatalError
-                volVectorField U(header, mesh);
+                    << "Unable to process " << fieldName << nl
-
+                    << "No call to magGrad for fields of type "
-                volScalarField magGrad
+                    << fieldHeader.headerClassName() << nl << nl
-                (
+                    << exit(FatalError);
-                    IOobject
-                    (
-                        "magGrad" + fieldName,
-                        runTime.timeName(),
-                        mesh,
-                        IOobject::NO_READ
-                    ),
-                    mag(fvc::grad(U))
-                );
-                Info<< "    Calculating " << magGrad.name() << endl;
-                magGrad.write();
-            }
-            else if (header.headerClassName() == volScalarField::typeName)
-            {
-                Info<< "    Reading " << fieldName << endl;
-                volScalarField U(header, mesh);
-
-                volScalarField magGrad
-                (
-                    IOobject
-                    (
-                        "magGrad" + fieldName,
-                        runTime.timeName(),
-                        mesh,
-                        IOobject::NO_READ
-                    ),
-                    mag(fvc::grad(U))
-                );
-                Info<< "    Calculating " << magGrad.name() << endl;
-                magGrad.write();
             }
         }
         else
         {
             Info<< "    No " << fieldName << endl;
         }
+
+        Info<< endl;
     }
 
+    Info<< "End\n" << endl;
+
     return(0);
 }
 

applications/utilities/postProcessing/field/magGrad/writeMagGradField.C

@@ -0,0 +1,72 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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
+
+Application
+    magGrad
+
+Description
+    Calculates and writes the magnitude of the gradient of a field for each
+    time
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+void writeMagGradField
+(
+    const IOobject& header,
+    const fvMesh& mesh,
+    bool& processed
+)
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    if (header.headerClassName() == fieldType::typeName)
+    {
+        Info<< "    Reading " << header.name() << endl;
+        fieldType field(header, mesh);
+
+        Info<< "    Calculating magGrad" << header.name() << endl;
+        volScalarField magGradField
+        (
+            IOobject
+            (
+                "magGrad" + header.name(),
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ
+            ),
+            mag(fvc::grad(field))
+        );
+        magGradField.write();
+
+        processed = true;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/utilities/postProcessing/field/magSqr/magSqr.C

@@ -34,39 +34,7 @@ Description
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-template<class Type>
+#include "writeMagSqrField.C"
-void writeMagSqrField
-(
-    const IOobject& header,
-    const fvMesh& mesh,
-    bool& processed
-)
-{
-    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
-
-    if (header.headerClassName() == fieldType::typeName)
-    {
-        Info<< "    Reading " << header.name() << endl;
-        fieldType field(header, mesh);
-
-        Info<< "    Calculating magSqr" << header.name() << endl;
-        volScalarField magSqrField
-        (
-            IOobject
-            (
-                "magSqr" + header.name(),
-                mesh.time().timeName(),
-                mesh,
-                IOobject::NO_READ
-            ),
-            magSqr(field)
-        );
-        magSqrField.write();
-
-        processed = true;
-    }
-}
-
 
 int main(int argc, char *argv[])
 {

applications/utilities/postProcessing/field/magSqr/writeMagSqrField.C

@@ -0,0 +1,72 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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
+
+Application
+    mag
+
+Description
+    Calculates and writes the magnitude-squared of a field for each time
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+void writeMagSqrField
+(
+    const IOobject& header,
+    const fvMesh& mesh,
+    bool& processed
+)
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    if (header.headerClassName() == fieldType::typeName)
+    {
+        Info<< "    Reading " << header.name() << endl;
+        fieldType field(header, mesh);
+
+        Info<< "    Calculating magSqr" << header.name() << endl;
+        volScalarField magSqrField
+        (
+            IOobject
+            (
+                "magSqr" + header.name(),
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ
+            ),
+            magSqr(field)
+        );
+        magSqrField.write();
+
+        processed = true;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+

applications/utilities/preProcessing/molConfig/Make/files

@@ -0,0 +1,8 @@
+latticeStructures = latticeStructures
+velocityDistributions = velocityDistributions
+
+createMolecules.C
+molConfig.C
+genMolConfig.C
+
+EXE = $(FOAM_APPBIN)/molConfig

applications/utilities/preProcessing/molConfig/Make/options

@@ -0,0 +1,15 @@
+EXE_INC = \
+    -I$(latticeStructures) \
+    -I$(velocityDistributions) \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecule/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools \
+    -ldynamicMesh \
+    -lfiniteVolume \
+    -llagrangian \
+    -lmolecule

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;

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;
+}
+
+
+// ************************************************************************* //

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);
+}

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"
+}

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;
+}
+
+
+// ************************************************************************* //

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()
+                    );
+            }
+        }
+    }
+}

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()
+                    );
+            }
+        }
+    }
+}

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);
+            }
+        }
+    }
+}
+

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()
+{}
+
+
+// ************************************************************************* //

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
+
+// ************************************************************************* //

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
+
+// ************************************************************************* //

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;
+

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;
+

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;
+

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);
+}
+
+

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);
+}

src/lagrangian/molecule/Make/files

@@ -0,0 +1,57 @@
+correlationFunction = correlationFunction
+
+distribution = distribution
+
+molecule = molecule
+moleculeCloud = moleculeCloud
+
+reducedUnits = reducedUnits
+
+referredMolecule = referredMolecule
+referredCellList = referredCellList
+referredCell = referredCell
+referralLists = referralLists
+
+potentials = potentials
+pairPotential = $(potentials)/pairPotential
+tetherPotential = $(potentials)/tetherPotential
+
+$(distribution)/distribution.C
+
+$(reducedUnits)/reducedUnits.C
+$(reducedUnits)/reducedUnitsIO.C
+
+$(molecule)/molecule.C
+$(molecule)/moleculeIO.C
+
+$(moleculeCloud)/moleculeCloud.C
+$(moleculeCloud)/moleculeCloudBuildCellOccupancy.C
+$(moleculeCloud)/moleculeCloudBuildCellInteractionLists.C
+$(moleculeCloud)/moleculeCloudBuildCellReferralLists.C
+
+$(moleculeCloud)/moleculeCloudTestEdgeEdgeDistance.C
+$(moleculeCloud)/moleculeCloudTestPointFaceDistance.C
+$(moleculeCloud)/moleculeCloudRealCellsInRangeOfSegment.C
+$(moleculeCloud)/moleculeCloudReferredCellsInRangeOfSegment.C
+
+$(moleculeCloud)/moleculeCloudCalculateForce.C
+$(moleculeCloud)/moleculeCloudCalculatePairForce.C
+$(moleculeCloud)/moleculeCloudCalculateTetherForce.C
+$(moleculeCloud)/moleculeCloudCalculateExternalForce.C
+$(moleculeCloud)/moleculeCloudIntegrateEquationsOfMotion.C
+$(moleculeCloud)/moleculeCloudRemoveHighEnergyOverlaps.C
+$(moleculeCloud)/moleculeCloudApplyConstraintsAndThermostats.C
+
+$(pairPotential)/basic/pairPotential.C
+$(pairPotential)/basic/pairPotentialList.C
+$(tetherPotential)/tetherPotential.C
+$(tetherPotential)/tetherPotentialList.C
+
+$(referralLists)/receivingReferralList.C
+$(referralLists)/sendingReferralList.C
+$(referredCellList)/referredCellList.C
+$(referredCell)/referredCell.C
+$(referredMolecule)/referredMolecule.C
+
+LIB = $(FOAM_LIBBIN)/libmolecule
+

src/lagrangian/molecule/Make/options

@@ -0,0 +1,8 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -llagrangian
+

src/lagrangian/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulator.C

@@ -0,0 +1,238 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "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"
+
+// ************************************************************************* //

src/lagrangian/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulator.H

@@ -0,0 +1,178 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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::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 const 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
+
+// ************************************************************************* //

src/lagrangian/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulatorI.H

@@ -0,0 +1,80 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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
+{
+
+// * * * * * * * * * * * * * 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 const 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
+
+// ************************************************************************* //

src/lagrangian/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulatorIO.C

@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source bAD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenbAD 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 "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;
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/correlationFunction/correlationFunction.C

@@ -0,0 +1,224 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "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().deltaT().value());
+
+    sampleInterval_ = sampleSteps_*mesh_.time().deltaT().value();
+
+    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"
+
+// ************************************************************************* //

src/lagrangian/molecule/correlationFunction/correlationFunction.H

@@ -0,0 +1,180 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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::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 const scalar duration() const;
+
+        inline const scalar sampleInterval() const;
+
+        inline const scalar averagingInterval() const;
+
+        inline const 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
+
+// ************************************************************************* //

src/lagrangian/molecule/correlationFunction/correlationFunctionI.H

@@ -0,0 +1,70 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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
+
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+inline const Field< Field<Type> >& Foam::correlationFunction<Type>::
+tZeroBuffers() const
+{
+    return tZeroBuffers_;
+}
+
+
+template<class Type>
+inline const scalar Foam::correlationFunction<Type>::duration() const
+{
+    return duration_;
+}
+
+
+template<class Type>
+inline const scalar Foam::correlationFunction<Type>::sampleInterval() const
+{
+    return sampleInterval_;
+}
+
+
+template<class Type>
+inline const scalar Foam::correlationFunction<Type>::averagingInterval() const
+{
+    return averagingInterval_;
+}
+
+
+template<class Type>
+inline const label Foam::correlationFunction<Type>::sampleSteps() const
+{
+    return sampleSteps_;
+}
+
+
+template<class Type>
+inline label Foam::correlationFunction<Type>::measurandFieldSize() const
+{
+    return tZeroBuffers_[0].size();
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/correlationFunction/correlationFunctionIO.C

@@ -0,0 +1,72 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "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;
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/distribution/distribution.C

@@ -0,0 +1,449 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "distribution.H"
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+distribution::distribution()
+:
+    Map<label>(),
+    binWidth_(1)
+{}
+
+
+distribution::distribution(const scalar binWidth)
+:
+    Map<label>(),
+    binWidth_(binWidth)
+{}
+
+
+distribution::distribution(const distribution& d)
+:
+    Map<label>(static_cast< Map<label> >(d)),
+    binWidth_(d.binWidth())
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+distribution::~distribution()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+label 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;
+}
+
+
+scalar distribution::approxTotalEntries() const
+{
+    scalar sumOfEntries = 0;
+
+    forAllConstIter(Map<label>, *this, iter)
+    {
+        sumOfEntries += scalar(iter());
+    }
+
+    return sumOfEntries;
+}
+
+
+scalar 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;
+}
+
+
+scalar 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 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 distribution::add(const label valueToAdd)
+{
+    add(scalar(valueToAdd));
+}
+
+
+void distribution::insertMissingKeys()
+{
+    iterator iter(this->begin());
+
+    List<label> keys = toc();
+
+    sort(keys);
+
+    label k;
+
+    if (keys.size() > 0)
+    {
+        for (k = keys[0]; k < keys[keys.size()-1]; k++)
+        {
+            iter = find(k);
+
+            if (iter == this->end())
+            {
+                this->insert(k,0);
+            }
+        }
+    }
+}
+
+
+List< Pair<scalar> > 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_;
+    }
+
+    return normDist;
+}
+
+
+List< Pair<scalar> > distribution::normalisedMinusMean()
+{
+    return normalisedShifted(mean());
+}
+
+
+List< Pair<scalar> > distribution::normalisedShifted(const 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;
+
+//     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++;
+    }
+
+//     checkNormalisation = 0;
+
+//     forAll (newDist, nD)
+//     {
+//         checkNormalisation += newDist[nD].second()*binWidth_;
+//     }
+
+//     Info << "Shifted normalisation = " << checkNormalisation << endl;
+
+    return newDist;
+}
+
+
+List< Pair<scalar> > 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 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  * * * * * * * * * * * * * //
+
+Ostream& 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;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/molecule/distribution/distribution.H

@@ -0,0 +1,133 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::distribution
+
+Description
+
+SourceFiles
+    distributionI.H
+    distribution.C
+    distributionIO.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:
+
+    // Constructors
+
+        //- Construct null
+        distribution();
+
+        //- Construct from binWidth
+        distribution(const scalar binWidth);
+
+        //- Construct as copy
+        distribution(const distribution&);
+
+
+    // Destructor
+
+        ~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(const scalar shiftValue);
+
+        List< Pair<scalar> > raw();
+
+        // Access
+
+            inline const scalar binWidth() const;
+
+
+    // Member Operators
+
+        void operator=(const distribution&);
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<<(Ostream&, const distribution&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "distributionI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/molecule/distribution/distributionI.H

@@ -0,0 +1,42 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 scalar distribution::binWidth() const
+{
+    return binWidth_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/molecule/distribution/distributionIO.C

@@ -0,0 +1,35 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "distribution.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// construct from Istream
+
+
+// ************************************************************************* //

src/lagrangian/molecule/mdTools/averageMDFields.H

@@ -0,0 +1,227 @@
+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;
+        }
+    }
+}
+

src/lagrangian/molecule/mdTools/calculateAutoCorrelationFunctions.H

@@ -0,0 +1,83 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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
+
+\*---------------------------------------------------------------------------*/
+
+if (mesh.time().timeIndex() % vacf.sampleSteps() == 0)
+{
+    IDLList<molecule>::iterator mol(molecules.begin());
+
+    Field<vector> uVals(molecules.size());
+
+    label uV = 0;
+
+    for
+    (
+        mol = molecules.begin();
+        mol != molecules.end();
+        ++mol, uV++
+    )
+    {
+        uVals[uV] = mol().U();
+    }
+
+    vacf.calculateCorrelationFunction(uVals);
+}
+
+if (mesh.time().timeIndex() % pacf.sampleSteps() == 0)
+{
+    IDLList<molecule>::iterator mol(molecules.begin());
+
+    vector p = vector::zero;
+
+    for
+    (
+        mol = molecules.begin();
+        mol != molecules.end();
+        ++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)
+{
+//     hFacf.calculateCorrelationFunction();
+}

src/lagrangian/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;
+    }
+}

src/lagrangian/molecule/mdTools/calculateTransportProperties.H

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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
+
+\*---------------------------------------------------------------------------*/
+
+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)
+{
+
+}
+
+

src/lagrangian/molecule/mdTools/createAutoCorrelationFunctions.H

@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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
+
+\*---------------------------------------------------------------------------*/
+
+Info << nl << "Creating autocorrelation functions." << endl;
+
+IOdictionary mdTransportProperitesDict
+(
+    IOobject
+    (
+        "mdTransportProperitesDict",
+        mesh.time().system(),
+        mesh,
+        IOobject::MUST_READ,
+        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")));

src/lagrangian/molecule/mdTools/createMDFields.H

@@ -0,0 +1,306 @@
+// Fields for data gathering
+
+List< scalarField > allSpeciesN_RU
+(
+    molecules.nIds(),
+    scalarField (mesh.nCells(), 0.0)
+);
+
+List< scalarField > allSpeciesM_RU
+(
+    molecules.nIds(),
+    scalarField (mesh.nCells(), 0.0)
+);
+
+List< vectorField > allSpeciesVelocitySum_RU
+(
+    molecules.nIds(),
+    vectorField (mesh.nCells(), vector::zero)
+);
+
+List< scalarField > allSpeciesVelocityMagSquaredSum_RU
+(
+    molecules.nIds(),
+    scalarField (mesh.nCells(), 0.0)
+);
+
+// Geometric Fields for IO
+
+Info << nl << "Creating fields." << endl;
+
+/*---------------------------------------------------------------------------*\
+    Number density
+\*---------------------------------------------------------------------------*/
+
+PtrList<volScalarField> allSpeciesRhoN
+(
+    molecules.nIds()
+);
+
+forAll (allSpeciesRhoN, rN)
+{
+    Info << "    Creating number density field for "
+        << molecules.idList()[rN] << endl;
+
+    allSpeciesRhoN.set
+    (
+        rN,
+        new volScalarField
+        (
+            IOobject
+            (
+                "rhoN_" + molecules.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.nIds()
+);
+
+forAll (allSpeciesRhoM, rM)
+{
+    Info << "    Creating mass density field for "
+        << molecules.idList()[rM] << endl;
+
+    allSpeciesRhoM.set
+    (
+        rM,
+        new volScalarField
+        (
+            IOobject
+            (
+                "rhoM_" + molecules.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.nIds()
+);
+
+forAll (allSpeciesVelocity, v)
+{
+    Info << "    Creating velocity field for "
+        << molecules.idList()[v] << endl;
+
+    allSpeciesVelocity.set
+    (
+        v,
+        new volVectorField
+        (
+            IOobject
+            (
+                "velocity_" + molecules.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();
+
+/*---------------------------------------------------------------------------*\
+    Kinetic temperature
+\*---------------------------------------------------------------------------*/
+
+PtrList<volScalarField> allSpeciesTemperature
+(
+    molecules.nIds()
+);
+
+forAll (allSpeciesTemperature, t)
+{
+    Info << "    Creating temperature field for "
+        << molecules.idList()[t] << endl;
+
+    allSpeciesTemperature.set
+    (
+        t,
+        new volScalarField
+        (
+            IOobject
+            (
+                "temperature_" + molecules.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.nIds()
+);
+
+forAll (allSpeciesMeanKE, mKE)
+{
+    Info << "    Creating mean kinetic energy field for "
+        << molecules.idList()[mKE] << endl;
+
+    allSpeciesMeanKE.set
+    (
+        mKE,
+        new volScalarField
+        (
+            IOobject
+            (
+                "meanKE_" + molecules.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();

src/lagrangian/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;

src/lagrangian/molecule/mdTools/md.H

@@ -0,0 +1,11 @@
+#ifndef md_H
+#define md_H
+
+#   include "moleculeCloud.H"
+
+#   include "correlationFunction.H"
+#   include "distribution.H"
+#   include "reducedUnits.H"
+
+#endif
+

src/lagrangian/molecule/mdTools/meanMomentumEnergyAndNMols.H

@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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
+
+Global
+    meanMomentumEnergyAndNMols.H
+
+Description
+    Calculates and prints the mean momentum and energy in the system
+    and the number of molecules.
+
+\*---------------------------------------------------------------------------*/
+
+
+vector singleStepTotalMomentum(vector::zero);
+
+scalar singleStepMaxVelocityMag = 0.0;
+
+scalar singleStepTotalMass = 0.0;
+
+scalar singleStepTotalKE = 0.0;
+
+scalar singleStepTotalPE = 0.0;
+
+scalar singleStepTotalrDotf = 0.0;
+
+{
+    IDLList<molecule>::iterator mol(molecules.begin());
+
+    for
+    (
+        mol = molecules.begin();
+        mol != molecules.end();
+        ++mol
+    )
+    {
+        const scalar molM(mol().mass());
+
+        const vector& molU(mol().U());
+
+        singleStepTotalMomentum += molU * molM;
+
+        singleStepTotalMass += molM;
+
+        if(mag(molU) > singleStepMaxVelocityMag)
+        {
+            singleStepMaxVelocityMag = mag(molU);
+        }
+
+        singleStepTotalKE += 0.5*molM*magSqr(molU);
+
+        singleStepTotalPE += mol().potentialEnergy();
+
+        singleStepTotalrDotf += tr(mol().rf());
+    }
+}
+
+label singleStepNMols = molecules.size();
+
+if (Pstream::parRun())
+{
+    reduce(singleStepTotalMomentum, sumOp<vector>());
+
+    reduce(singleStepMaxVelocityMag, maxOp<scalar>());
+
+    reduce(singleStepTotalMass, sumOp<scalar>());
+
+    reduce(singleStepTotalKE, sumOp<scalar>());
+
+    reduce(singleStepTotalPE, sumOp<scalar>());
+
+    reduce(singleStepTotalrDotf, sumOp<scalar>());
+
+    reduce(singleStepNMols, sumOp<label>());
+}
+
+if (singleStepNMols)
+{
+    Info << "Number of mols in system = "
+            << singleStepNMols << nl
+        << "Overall number density = "
+            << singleStepNMols/meshVolume << " m^-3" << nl
+        << "Overall mass density = "
+            << singleStepTotalMass/meshVolume << " kg/m^3" << nl
+        << "Average velocity per mol = "
+            << singleStepTotalMomentum/singleStepTotalMass << " m/s" << nl
+        << "Maximum |velocity| = "
+            << singleStepMaxVelocityMag << " m/s" << nl
+        << "Average KE per mol = "
+            << singleStepTotalKE/singleStepNMols << " J" << nl
+        << "Average PE per mol = "
+            << singleStepTotalPE/singleStepNMols << " J" << nl
+        << "Average TE per mol = "
+            << (singleStepTotalKE + singleStepTotalPE)/singleStepNMols << " J"
+        << endl;
+
+//     Info << singleStepNMols << " "
+// //     << singleStepTotalMomentum/singleStepTotalMass << " "
+//     << singleStepMaxVelocityMag << " "
+//     << singleStepTotalKE/singleStepNMols << " "
+//     << singleStepTotalPE/singleStepNMols << " "
+//     << (singleStepTotalKE + singleStepTotalPE)/singleStepNMols << endl;
+}
+else
+{
+    Info << "No molecules in system" << endl;
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/mdTools/resetMDFields.H

@@ -0,0 +1,26 @@
+if (runTime.outputTime())
+{
+    allSpeciesN_RU = List< scalarField >
+    (
+        molecules.nIds(),
+        scalarField (mesh.nCells(), 0.0)
+    );
+
+    allSpeciesM_RU = List< scalarField >
+    (
+        molecules.nIds(),
+        scalarField (mesh.nCells(), 0.0)
+    );
+
+    allSpeciesVelocitySum_RU = List< vectorField >
+    (
+        molecules.nIds(),
+        vectorField (mesh.nCells(), vector::zero)
+    );
+
+    allSpeciesVelocityMagSquaredSum_RU = List< scalarField >
+    (
+        molecules.nIds(),
+        scalarField (mesh.nCells(), 0.0)
+    );
+}

src/lagrangian/molecule/mdTools/temperatureAndPressure.H

@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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
+
+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.
+
+\*---------------------------------------------------------------------------*/
+
+accumulatedTotalMomentum += singleStepTotalMomentum;
+
+accumulatedTotalMass += singleStepTotalMass;
+
+accumulatedTotalKE += singleStepTotalKE;
+
+accumulatedTotalPE += singleStepTotalPE;
+
+accumulatedTotalrDotfSum += singleStepTotalrDotf;
+
+accumulatedNMols += singleStepNMols;
+
+if (runTime.outputTime())
+{
+    // calculate averages
+
+    if (accumulatedNMols)
+    {
+        Info << "calculating averages" << endl;
+
+        averageTemperature =
+        (
+            2.0/(3.0 * moleculeCloud::kb * accumulatedNMols)
+            *
+            (
+                accumulatedTotalKE
+                -
+                0.5*magSqr(accumulatedTotalMomentum)/accumulatedTotalMass
+            )
+        );
+
+        averagePressure =
+        (
+            (
+                (accumulatedNMols/nAveragingSteps)
+                *
+                moleculeCloud::kb * averageTemperature
+                +
+                accumulatedTotalrDotfSum/(6.0 * nAveragingSteps)
+            )
+            /
+            meshVolume
+        );
+
+        // output values
+
+        Info << "----------------------------------------" << nl
+            << "Averaged properties" << nl
+            << "Average |velocity| = "
+                << mag(accumulatedTotalMomentum)/accumulatedTotalMass
+                << " m/s" << nl
+            << "Average temperature = "
+                << averageTemperature << " K" << nl
+            << "Average pressure = "
+                << averagePressure << " N/m^2" << nl
+            << "----------------------------------------" << endl;
+    }
+    else
+    {
+        Info << "Not averaging temperature and pressure: "
+            << "no molecules in system" << endl;
+    }
+
+    // reset counters
+
+    accumulatedTotalMomentum = vector::zero;
+
+    accumulatedTotalMass = 0.0;
+
+    accumulatedTotalKE = 0.0;
+
+    accumulatedTotalPE = 0.0;
+
+    accumulatedTotalrDotfSum = 0.0;
+
+    accumulatedNMols = 0;
+}
+
+// ************************************************************************* //

src/lagrangian/molecule/mdTools/temperatureAndPressureVariables.H

@@ -0,0 +1,59 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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
+
+Global
+    temperatureAndPressureVariables.H
+
+Description
+    Provides accumulation variables for temperatureAndPressure.H
+
+\*---------------------------------------------------------------------------*/
+
+vector accumulatedTotalMomentum(vector::zero);
+
+scalar accumulatedTotalMass = 0.0;
+
+scalar accumulatedTotalKE = 0.0;
+
+scalar accumulatedTotalPE = 0.0;
+
+scalar accumulatedTotalrDotfSum = 0.0;
+
+label accumulatedNMols = 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>());
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/mdTools/temperatureEquilibration.H

@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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
+
+Global
+    temperatureEquilibration.H
+
+Description
+    Applies temperature control to molecules
+
+\*---------------------------------------------------------------------------*/
+
+if (runTime.outputTime())
+{
+    molecules.applyConstraintsAndThermostats
+    (
+        targetTemperature,
+        averageTemperature
+    );
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/molecule/molecule.C

@@ -0,0 +1,214 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "moleculeCloud.H"
+#include "Random.H"
+#include "Time.H"
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+bool molecule::move(molecule::trackData& td)
+{
+    td.switchProcessor = false;
+    td.keepParticle = true;
+
+    scalar deltaT = cloud().pMesh().time().deltaT().value();
+    scalar tEnd = (1.0 - stepFraction())*deltaT;
+    scalar dtMax = tEnd;
+
+    moleculeCloud::integrationMethods method
+            = td.molCloud().integrationMethod();
+
+    if (method == moleculeCloud::imVerletLeapfrog)
+    {
+        if (td.part() == 1)  // Leapfrog 1st Part
+        {
+            if (stepFraction() < VSMALL)
+            {
+                U_ += 0.5*deltaT*A_;
+            }
+
+            while (td.keepParticle && !td.switchProcessor && tEnd > SMALL)
+            {
+                // set the lagrangian time-step
+                scalar dt = min(dtMax, tEnd);
+
+                dt *= trackToFace(position() + dt*U_, td);
+
+                tEnd -= dt;
+                stepFraction() = 1.0 - tEnd/deltaT;
+            }
+        }
+        else if (td.part() == 2)  // Leapfrog 2nd Part
+        {
+            U_ += 0.5*deltaT*A_;
+        }
+        else
+        {
+            FatalErrorIn("molecule::move(molecule::trackData& td)") << nl
+                << td.part()
+                << " is an invalid part of integration method: "
+                << method << nl
+                << abort(FatalError);
+        }
+    }
+    else if (method == moleculeCloud::imPredictorCorrector)
+    {
+        if (td.part() == 1) // Predictor Part
+        {
+
+        }
+        else if (td.part() == 2) // Corrector Part
+        {
+
+        }
+        else
+        {
+            FatalErrorIn("molecule::move(molecule::trackData& td)") << nl
+                << td.part() << " is an invalid part of integration method: "
+                << method
+                << abort(FatalError);
+        }
+    }
+    else
+    {
+        FatalErrorIn("molecule::move(molecule::trackData& td)") << nl
+                << "Unknown integration method: "
+                << method
+                << abort(FatalError);
+    }
+
+    return td.keepParticle;
+}
+
+
+void molecule::transformProperties(const tensor& T)
+{}
+
+
+void molecule::transformProperties(const vector& separation)
+{
+    if (tethered_)
+    {
+        tetherPosition_ += separation;
+    }
+}
+
+
+void molecule::hitProcessorPatch
+(
+    const processorPolyPatch&,
+    molecule::trackData& td
+)
+{
+    td.switchProcessor = true;
+}
+
+
+void molecule::hitProcessorPatch
+(
+    const processorPolyPatch&,
+    int&
+)
+{}
+
+
+void molecule::hitWallPatch
+(
+    const wallPolyPatch& wpp,
+    molecule::trackData& td
+)
+{
+    vector nw = wpp.faceAreas()[wpp.whichFace(face())];
+    nw /= mag(nw);
+
+    scalar Un = U_ & nw;
+//     vector Ut = U_ - Un*nw;
+
+//     Random rand(clock::getTime());
+
+//     scalar tmac = 0.8;
+
+//     scalar wallTemp = 2.5;
+
+//     if (rand.scalar01() < tmac)
+//     {
+//         // Diffuse reflection
+//
+//         vector tw1 = Ut/mag(Ut);
+//
+//         vector tw2 = nw ^ tw1;
+//
+//         U_ = sqrt(wallTemp/mass_)*rand.GaussNormal()*tw1
+//                 + sqrt(wallTemp/mass_)*rand.GaussNormal()*tw2
+//                 - mag(sqrt(wallTemp/mass_)*rand.GaussNormal())*nw;
+//     }
+
+//     else
+//     {
+        // Specular reflection
+
+        if (Un > 0)
+        {
+            U_ -= 2*Un*nw;
+        }
+//     }
+
+}
+
+
+void molecule::hitWallPatch
+(
+    const wallPolyPatch&,
+    int&
+)
+{}
+
+
+void molecule::hitPatch
+(
+    const polyPatch&,
+    molecule::trackData& td
+)
+{
+    td.keepParticle = false;
+}
+
+
+void molecule::hitPatch
+(
+    const polyPatch&,
+    int&
+)
+{}
+
+} // End namespace Foam
+
+
+// ************************************************************************* //

src/lagrangian/molecule/molecule/molecule.H

@@ -0,0 +1,270 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    molecule
+
+Description
+
+SourceFiles
+    moleculeI.H
+    molecule.C
+    moleculeIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef molecule_H
+#define molecule_H
+
+#include "Particle.H"
+#include "IOstream.H"
+#include "autoPtr.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Class forward declarations
+class moleculeCloud;
+
+/*---------------------------------------------------------------------------*\
+                           Class molecule Declaration
+\*---------------------------------------------------------------------------*/
+
+class molecule
+:
+    public Particle<molecule>
+{
+    // Private data
+
+        //- Be careful with the ordering of data.  It has an impact on binary
+        //  transfer:
+        //  1) Put the largest data members 1st
+        //  2) Pair up labels,
+        //  3) Don't go scalar-label, scalar-label, becasue in 64bit mode,
+        //     the labels will be padded by 4bytes.
+
+        // - mass of molecule
+        scalar mass_;
+
+        // - Velocity of molecule
+        vector U_;
+
+        // - Acceleration of molecule
+        vector A_;
+
+        // - Tether position
+        vector tetherPosition_;
+
+        // - Potential energy that this molecules posseses
+        scalar potentialEnergy_;
+
+        // - r_ij f_ij, stress dyad
+        tensor rf_;
+
+        // - Is the molecule tethered?
+        label tethered_;
+
+        // - id (type) of molecule
+        label id_;
+
+
+public:
+
+    friend class Cloud<molecule>;
+
+    //- Class used to pass tracking data to the trackToFace function
+    class trackData
+    :
+        public Particle<molecule>::trackData
+    {
+        //- Reference to the cloud containing this particle
+        moleculeCloud& molCloud_;
+
+        // label specifying which part of the integration algorithm is taking
+        // place (i.e. leapfrog 1 or leapfrog 2. Predictor or Corrector)
+        label part_;
+
+
+    public:
+
+        // Constructors
+
+            inline trackData
+            (
+                moleculeCloud& molCloud,
+                label part
+            );
+
+        // Member functions
+
+            inline moleculeCloud& molCloud();
+
+            inline label part() const;
+    };
+
+    // Constructors
+
+        //- Construct from components
+        inline molecule
+        (
+            const Cloud<molecule>& c,
+            const vector& position,
+            const label celli,
+            const scalar mass,
+            const vector& U,
+            const vector& A,
+            const vector& tetherPosition,
+            const label tethered,
+            const label id
+        );
+
+        //- Construct from Istream
+        molecule
+        (
+            const Cloud<molecule>& c,
+            Istream& is,
+            bool readFields = true
+        );
+
+        //- Construct and return a clone
+        autoPtr<molecule> clone() const
+        {
+            return autoPtr<molecule>(new molecule(*this));
+        }
+
+
+    // Member Functions
+
+        void transformProperties(const tensor& T);
+
+        void transformProperties(const vector& separation);
+
+        // Access
+
+            //- Return id
+            inline const label id() const;
+
+            //- Return mass
+            inline const scalar mass() const;
+
+            //- Return velocity
+            inline const vector& U() const;
+            inline vector& U();
+
+            //- Return acceleration
+            inline const vector& A() const;
+            inline vector& A();
+
+            //- Return potential energy
+            inline const scalar potentialEnergy() const;
+            inline scalar& potentialEnergy();
+
+            //- Return stress contribution
+            inline const tensor& rf() const;
+            inline tensor& rf();
+
+            //- Return tethered
+            inline const label tethered() const;
+
+            //- Return tetherPosition
+            inline const vector& tetherPosition() const;
+
+
+        //- Tracking
+        bool move(trackData&);
+
+
+    // Member Operators
+
+        //- Overridable function to handle the particle hitting a
+        //  processorPatch
+        void hitProcessorPatch
+        (
+            const processorPolyPatch&,
+            molecule::trackData& td
+        );
+
+        //- Overridable function to handle the particle hitting a
+        //  processorPatch without trackData
+        void hitProcessorPatch
+        (
+            const processorPolyPatch&,
+            int&
+        );
+
+        //- Overridable function to handle the particle hitting a wallPatch
+        void hitWallPatch
+        (
+            const wallPolyPatch&,
+            molecule::trackData& td
+        );
+
+        //- Overridable function to handle the particle hitting a wallPatch
+        //- without trackData
+        void hitWallPatch
+        (
+            const wallPolyPatch&,
+            int&
+        );
+
+        //- Overridable function to handle the particle hitting a polyPatch
+        void hitPatch
+        (
+            const polyPatch&,
+            molecule::trackData& td
+        );
+
+        //- Overridable function to handle the particle hitting a polyPatch
+        //- without trackData
+        void hitPatch
+        (
+            const polyPatch&,
+            int&
+        );
+
+        static void readFields(moleculeCloud& mC);
+
+        static void writeFields(const moleculeCloud& mC);
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<<(Ostream&, const molecule&);
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "moleculeI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/molecule/molecule/moleculeI.H

@@ -0,0 +1,159 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+inline molecule::molecule
+(
+    const Cloud<molecule>& c,
+    const vector& position,
+    const label celli,
+    const scalar mass,
+    const vector& U,
+    const vector& A,
+    const vector& tetherPosition,
+    const label tethered,
+    const label id
+)
+:
+    Particle<molecule>(c, position, celli),
+    mass_(mass),
+    U_(U),
+    A_(A),
+    tetherPosition_(tetherPosition),
+    potentialEnergy_(0.0),
+    rf_(tensor::zero),
+    tethered_(tethered),
+    id_(id)
+{}
+
+
+inline molecule::trackData::trackData
+(
+    moleculeCloud& molCloud,
+    label part
+)
+:
+    Particle<molecule>::trackData(refCast<Cloud<molecule> >(molCloud)),
+    molCloud_(molCloud),
+    part_(part)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const label molecule::id() const
+{
+    return id_;
+}
+
+
+inline const scalar molecule::mass() const
+{
+    return mass_;
+}
+
+
+inline const vector& molecule::U() const
+{
+    return U_;
+}
+
+
+inline vector& molecule::U()
+{
+    return U_;
+}
+
+
+inline const vector& molecule::A() const
+{
+    return A_;
+}
+
+
+inline vector& molecule::A()
+{
+    return A_;
+}
+
+
+inline const scalar molecule::potentialEnergy() const
+{
+    return potentialEnergy_;
+}
+
+
+inline scalar& molecule::potentialEnergy()
+{
+    return potentialEnergy_;
+}
+
+
+inline const tensor& molecule::rf() const
+{
+    return rf_;
+}
+
+
+inline tensor& molecule::rf()
+{
+    return rf_;
+}
+
+
+inline const label molecule::tethered() const
+{
+    return tethered_;
+}
+
+
+inline const vector& molecule::tetherPosition() const
+{
+    return tetherPosition_;
+}
+
+
+inline moleculeCloud& molecule::trackData::molCloud()
+{
+    return molCloud_;
+}
+
+
+inline label molecule::trackData::part() const
+{
+    return part_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/molecule/molecule/moleculeIO.C

@@ -0,0 +1,208 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "molecule.H"
+#include "IOstreams.H"
+
+#include "moleculeCloud.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::molecule::molecule
+(
+    const Cloud<molecule>& cloud,
+    Istream& is,
+    bool readFields
+)
+:
+    Particle<molecule>(cloud, is)
+{
+    if (readFields)
+    {
+        if (is.format() == IOstream::ASCII)
+        {
+            id_ = readLabel(is);
+            mass_ = readScalar(is);
+            is >> U_;
+            is >> A_;
+            is >> potentialEnergy_;
+            is >> rf_;
+            is >> tethered_;
+            is >> tetherPosition_;
+        }
+        else
+        {
+            is.read
+            (
+                reinterpret_cast<char*>(&mass_),
+                sizeof(mass_)
+                + sizeof(U_)
+                + sizeof(A_)
+                + sizeof(tetherPosition_)
+                + sizeof(potentialEnergy_)
+                + sizeof(rf_)
+                + sizeof(tethered_)
+                + sizeof(id_)
+            );
+        }
+    }
+
+    // Check state of Istream
+    is.check("Foam::molecule::molecule(Foam::Istream&)");
+}
+
+
+namespace Foam
+{
+
+void molecule::readFields(moleculeCloud& mC)
+{
+    if (!mC.size())
+    {
+        return;
+    }
+
+    IOField<label> id(mC.fieldIOobject("id"));
+    mC.checkFieldIOobject(mC, id);
+
+    IOField<scalar> mass(mC.fieldIOobject("mass"));
+    mC.checkFieldIOobject(mC, mass);
+
+    IOField<vector> U(mC.fieldIOobject("U"));
+    mC.checkFieldIOobject(mC, U);
+
+    IOField<vector> A(mC.fieldIOobject("A"));
+    mC.checkFieldIOobject(mC, A);
+
+    IOField<label> tethered(mC.fieldIOobject("tethered"));
+    mC.checkFieldIOobject(mC, tethered);
+
+    IOField<vector> tetherPositions(mC.fieldIOobject("tetherPositions"));
+    mC.checkFieldIOobject(mC, tetherPositions);
+
+    label i = 0;
+    forAllIter(moleculeCloud, mC, iter)
+    {
+        molecule& mol = iter();
+
+        mol.id_ = id[i];
+        mol.mass_ = mass[i];
+        mol.U_ = U[i];
+        mol.A_ = A[i];
+        mol.potentialEnergy_ = 0.0;
+        mol.rf_ = tensor::zero;
+        mol.tethered_ = tethered[i];
+        mol.tetherPosition_ = tetherPositions[i];
+        i++;
+    }
+}
+
+
+void molecule::writeFields(const moleculeCloud& mC)
+{
+    Particle<molecule>::writeFields(mC);
+
+    label np =  mC.size();
+
+    IOField<label> id(mC.fieldIOobject("id"), np);
+    IOField<scalar> mass(mC.fieldIOobject("mass"), np);
+    IOField<vector> U(mC.fieldIOobject("U"), np);
+    IOField<vector> A(mC.fieldIOobject("A"), np);
+    IOField<label> tethered(mC.fieldIOobject("tethered"), np);
+    IOField<vector> tetherPositions(mC.fieldIOobject("tetherPositions"), np);
+
+    label i = 0;
+    forAllConstIter(moleculeCloud, mC, iter)
+    {
+        const molecule& mol = iter();
+
+        id[i] = mol.id_;
+        mass[i] = mol.mass_;
+        U[i] = mol.U_;
+        A[i] = mol.A_;
+        tethered[i] = mol.tethered_;
+        tetherPositions[i] = mol.tetherPosition_;
+        i++;
+    }
+
+    id.write();
+    mass.write();
+    U.write();
+    A.write();
+    tethered.write();
+    tetherPositions.write();
+}
+
+};  // end of namespace Foam
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const molecule& mol)
+{
+    if (os.format() == IOstream::ASCII)
+    {
+        os  << mol.id_
+                << token::SPACE << mol.mass_
+                << token::SPACE << static_cast<const Particle<molecule>&>(mol)
+                << token::SPACE << mol.face()
+                << token::SPACE << mol.stepFraction()
+                << token::SPACE << mol.U_
+                << token::SPACE << mol.A_
+                << token::SPACE << mol.potentialEnergy_
+                << token::SPACE << mol.rf_
+                << token::SPACE << mol.tethered_
+                << token::SPACE << mol.tetherPosition_;
+    }
+    else
+    {
+        os  << static_cast<const Particle<molecule>&>(mol);
+        os.write
+        (
+            reinterpret_cast<const char*>(&mol.mass_),
+            sizeof(mol.mass_)
+            + sizeof(mol.U_)
+            + sizeof(mol.A_)
+            + sizeof(mol.tetherPosition_)
+            + sizeof(mol.potentialEnergy_)
+            + sizeof(mol.rf_)
+            + sizeof(mol.tethered_)
+            + sizeof(mol.id_)
+        );
+    }
+
+    // Check state of Ostream
+    os.check
+    (
+        "Foam::Ostream& Foam::operator<<"
+        "(Foam::Ostream&, const Foam::molecule&)"
+    );
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/moleculeCloud/moleculeCloud.C

@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "moleculeCloud.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineParticleTypeNameAndDebug(molecule, 0);
+    defineTemplateTypeNameAndDebug(Cloud<molecule>, 0);
+};
+
+
+template<>
+const char* Foam::NamedEnum<Foam::moleculeCloud::
+integrationMethods, 2>::names[] =
+{
+    "verletLeapfrog",
+    "predictorCorrector"
+};
+
+const Foam::NamedEnum<Foam::moleculeCloud::integrationMethods, 2>
+        Foam::moleculeCloud::integrationMethodNames_;
+
+
+Foam::scalar Foam::moleculeCloud::transTol = 1e-12;
+
+Foam::scalar Foam::moleculeCloud::kb = 1.380650277e-23;
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::moleculeCloud::moleculeCloud(const polyMesh& mesh)
+:
+    Cloud<molecule>(mesh, "moleculeCloud", false),
+    mesh_(mesh),
+    referredInteractionList_(*this)
+{
+    molecule::readFields(*this);
+
+#   include "moleculeCloudReadMDParameters.H"
+
+    buildCellInteractionLists();
+
+    buildCellReferralLists();
+
+    buildCellOccupancy();
+}
+
+
+Foam::moleculeCloud::moleculeCloud
+(
+    const polyMesh& mesh,
+    label nMol,
+    const labelField& id,
+    const scalarField& mass,
+    const vectorField& positions,
+    const labelField& cells,
+    const vectorField& U,
+    const vectorField& A,
+    const labelField& tethered,
+    const vectorField& tetherPositions
+)
+:
+    Cloud<molecule>(mesh, "moleculeCloud", false),
+    mesh_(mesh),
+    referredInteractionList_(*this)
+{
+    molecule::readFields(*this);
+
+    clear();
+
+    // This clear ()is here for the moment to stop existing files
+    // being appended to, this would be better accomplished by getting
+    // mesh.removeFiles(mesh.instance()); (or equivalent) to work.
+
+    int i;
+
+    const Cloud<molecule>& cloud = *this;
+
+    for (i=0; i<nMol; i++)
+    {
+        addParticle
+        (
+            new molecule
+            (
+                cloud,
+                positions[i],
+                cells[i],
+                mass[i],
+                U[i],
+                A[i],
+                tetherPositions[i],
+                tethered[i],
+                id[i]
+            )
+        );
+    }
+};
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::moleculeCloud::writeFields() const
+{
+    molecule::writeFields(*this);
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/moleculeCloud/moleculeCloud.H

@@ -0,0 +1,350 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::moleculeCloud
+
+Description
+
+SourceFiles
+    moleculeCloudApplyConstraintsAndThermostats.C
+    moleculeCloudBuildCellInteractionLists.C
+    moleculeCloudBuildCellOccupancy.C
+    moleculeCloudBuildCellReferralLists.C
+    moleculeCloud.C
+    moleculeCloudCalculateAndAccumulateProperties.C
+    moleculeCloudCalculateExternalForce.C
+    moleculeCloudCalculateForce.C
+    moleculeCloudCalculatePairForce.C
+    moleculeCloudCalculateTetherForce.C
+    moleculeCloudIntegrateEquationsOfMotion.C
+    moleculeCloudRemoveHighEnergyOverlaps.C
+
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef moleculeCloud_H
+#define moleculeCloud_H
+
+#include "Cloud.H"
+#include "molecule.H"
+#include "IOdictionary.H"
+#include "vector2D.H"
+
+#include "pairPotentialList.H"
+#include "tetherPotentialList.H"
+
+#include "receivingReferralList.H"
+#include "sendingReferralList.H"
+#include "referredCellList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class moleculeCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+class moleculeCloud
+:
+    public Cloud<molecule>
+{
+
+public:
+
+    enum integrationMethods
+    {
+        imVerletLeapfrog,
+        imPredictorCorrector
+    };
+
+
+private:
+
+    // Private data
+
+        const polyMesh& mesh_;
+
+            // MD solution parameters
+
+            static const NamedEnum<integrationMethods, 2>
+                integrationMethodNames_;
+
+            integrationMethods integrationMethod_;
+
+            scalar potentialEnergyLimit_;
+
+            scalar guardRadius_;
+
+            scalar rCutMax_;
+
+            //- storing rCutMaxSqr in class as well as rCutMax to
+            //- avoid needing to calculate it often.
+            //- Possibilty of inconsistency if tinkered with.
+
+            scalar rCutMaxSqr_;
+
+            List<word> idList_;
+
+            labelList removalOrder_;
+
+            labelListList directInteractionList_;
+
+            referredCellList referredInteractionList_;
+
+            labelList realCellsWithinRCutMaxOfAnyReferringPatch_;
+
+            labelList realFacesWithinRCutMaxOfAnyReferringPatch_;
+
+            labelList realEdgesWithinRCutMaxOfAnyReferringPatch_;
+
+            labelList realPointsWithinRCutMaxOfAnyReferringPatch_;
+
+            List<sendingReferralList> cellSendingReferralLists_;
+
+            List<receivingReferralList> cellReceivingReferralLists_;
+
+            pairPotentialList pairPotentials_;
+
+            tetherPotentialList tetherPotentials_;
+
+            vector gravity_;
+
+            List< DynamicList<molecule*> > cellOccupancy_;
+
+
+        // Private Member Functions
+
+            //- Disallow default bitwise copy construct
+            moleculeCloud(const moleculeCloud&);
+
+            //- Disallow default bitwise assignment
+            void operator=(const moleculeCloud&);
+
+
+public:
+
+    // Static data members
+
+        //- Tolerance for checking that faces on a patch segment
+        static scalar transTol;
+
+        static scalar kb;
+
+
+    // Constructors
+
+        //- Construct given mesh
+        moleculeCloud(const polyMesh&);
+
+        //- Construct given polyMesh and fields of position, cell, mass,
+        //- id, U ands A.  Intended for use by the molConfig utility
+        moleculeCloud
+        (
+            const polyMesh& mesh,
+            label nMol,
+            const labelField& id,
+            const scalarField& mass,
+            const vectorField& positions,
+            const labelField& cells,
+            const vectorField& U,
+            const vectorField& A,
+            const labelField& tethered,
+            const vectorField& tetherPositions
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            inline const polyMesh& mesh() const;
+
+        // MD solution parameters
+
+            inline const integrationMethods& integrationMethod() const;
+
+            inline scalar potentialEnergyLimit() const;
+
+            inline scalar guardRadius() const;
+
+            inline const scalar rCutMax() const;
+
+            inline const scalar rCutMaxSqr() const;
+
+            inline const List<word>& idList() const;
+
+            inline const labelList& removalOrder() const;
+
+            inline label nPairPotentials() const;
+
+            inline label nIds() const;
+
+            inline const labelListList& directInteractionList() const;
+
+            inline const referredCellList& referredInteractionList() const;
+
+            inline const labelList&
+                realCellsWithinRCutMaxOfAnyReferringPatch() const;
+
+            inline const labelList&
+                realFacesWithinRCutMaxOfAnyReferringPatch() const;
+
+            inline const labelList&
+                realEdgesWithinRCutMaxOfAnyReferringPatch() const;
+
+            inline const labelList&
+                realPointsWithinRCutMaxOfAnyReferringPatch() const;
+
+            inline const List<sendingReferralList>&
+                cellSendingReferralLists() const;
+
+            inline const List<receivingReferralList>&
+                cellReceivingReferralLists() const;
+
+            inline label nInteractingProcs() const;
+
+            inline const pairPotentialList& pairPotentials() const;
+
+            inline const tetherPotentialList& tetherPotentials() const;
+
+            inline const vector& gravity() const;
+
+            inline const List< DynamicList<molecule*> >& cellOccupancy() const;
+
+            void buildCellInteractionLists();
+
+            //- Build referralLists which define how to send information
+            // to referredCells to source cells
+            void buildCellReferralLists();
+
+            bool testPointFaceDistance
+            (
+                const label p,
+                const label faceNo
+            ) const;
+
+            bool testPointFaceDistance
+            (
+                const label p,
+                const referredCell& refCell
+            ) const;
+
+            bool testPointFaceDistance
+            (
+                const vectorList& pointsToTest,
+                const label faceNo
+            ) const;
+
+            bool testPointFaceDistance
+            (
+                const vector& p,
+                const label faceNo
+            ) const;
+
+            bool testPointFaceDistance
+            (
+                const vector& p,
+                const labelList& faceToTest,
+                const vectorList& points,
+                const vector& faceC,
+                const vector& faceA
+            ) const;
+
+            bool testEdgeEdgeDistance
+            (
+                const edge& eI,
+                const edge& eJ
+            ) const;
+
+            bool testEdgeEdgeDistance
+            (
+                const edge& eI,
+                const vector& eJs,
+                const vector& eJe
+            ) const;
+
+            const labelList realCellsInRangeOfSegment
+            (
+                const labelList& segmentFaces,
+                const labelList& segmentEdges,
+                const labelList& segmentPoints
+            ) const;
+
+            const labelList referredCellsInRangeOfSegment
+            (
+                const List<referredCell>& referredInteractionList,
+                const labelList& segmentFaces,
+                const labelList& segmentEdges,
+                const labelList& segmentPoints
+            ) const;
+
+            //- Determine which molecules are in which cells
+            void buildCellOccupancy();
+
+            //- Integrate the equations of motion using algorithm selected at
+            //  runtime from a dictionary. This will also call the function
+            //  to calculate the intermolecular forces (calculatePairForce()).
+
+            void integrateEquationsOfMotion();
+
+            void applyConstraintsAndThermostats
+            (
+                const scalar targetTemperature,
+                const scalar measuredTemperature
+            );
+
+            void calculateForce();
+
+            void calculatePairForce();
+
+            void calculateTetherForce();
+
+            void calculateExternalForce();
+
+            void removeHighEnergyOverlaps();
+
+
+    // Member Operators
+
+        //- Write fields
+        void writeFields() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "moleculeCloudI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/molecule/moleculeCloud/moleculeCloudApplyConstraintsAndThermostats.C

@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "moleculeCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::moleculeCloud::applyConstraintsAndThermostats
+(
+    const scalar targetTemperature,
+    const scalar measuredTemperature
+)
+{
+    scalar temperatureCorrectionFactor =
+        sqrt(targetTemperature/measuredTemperature);
+
+    Info<< "----------------------------------------" << nl
+        << "Temperature equilibration" << nl
+        << "Target temperature = "
+        << targetTemperature << nl
+        << "Measured temperature = "
+        << measuredTemperature << nl
+        << "Temperature correction factor ="
+        << temperatureCorrectionFactor << nl
+        << "----------------------------------------"
+        << endl;
+
+    iterator mol(this->begin());
+
+    for (mol = this->begin(); mol != this->end(); ++mol)
+    {
+        mol().U() *= temperatureCorrectionFactor;
+    }
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/moleculeCloud/moleculeCloudBuildCellInteractionLists.C

@@ -0,0 +1,80 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "moleculeCloud.H"
+#include "polyBoundaryMeshEntries.H"
+#include "PstreamCombineReduceOps.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::moleculeCloud::buildCellInteractionLists()
+{
+// #       include "moleculeCloudCodeSnippets.H"
+
+    List<DynamicList<label> > directInteractionList(mesh_.nCells());
+
+#   include "moleculeCloudBuildDirectInteractionList.H"
+
+    directInteractionList_.setSize(mesh_.nCells());
+
+    forAll(directInteractionList, transDIL)
+    {
+        directInteractionList_[transDIL].transfer
+        (
+            directInteractionList[transDIL].shrink()
+        );
+    }
+
+//     for sorting DILs
+
+//     forAll(directInteractionList_, dIL)
+//     {
+//         sort(directInteractionList_[dIL]);
+//     }
+
+    cellOccupancy_.setSize(mesh_.nCells());
+
+    DynamicList<referredCell> referredInteractionList;
+
+#   include "moleculeCloudBuildReferredInteractionList.H"
+
+//#   include "moleculeCloudBuildReferredInteractionList_old.H"
+
+    referredInteractionList_.setSize
+    (
+        referredInteractionList.size()
+    );
+
+    forAll(referredInteractionList, rIL)
+    {
+        referredInteractionList_[rIL] = referredInteractionList[rIL];
+    }
+
+    referredInteractionList_.setRealCellsInRange();
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/moleculeCloud/moleculeCloudBuildCellOccupancy.C

@@ -0,0 +1,54 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "moleculeCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::moleculeCloud::buildCellOccupancy()
+{
+    forAll(cellOccupancy_, cO)
+    {
+        cellOccupancy_[cO].clear();
+    }
+
+    iterator mol(this->begin());
+
+    for (mol = this->begin(); mol != this->end(); ++mol)
+    {
+        cellOccupancy_[mol().cell()].append(&mol());
+    }
+
+    forAll(cellOccupancy_, cO)
+    {
+        cellOccupancy_[cO].shrink();
+    }
+
+    referredInteractionList_.referMolecules();
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/moleculeCloud/moleculeCloudBuildCellReferralLists.C

@@ -0,0 +1,194 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "moleculeCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::moleculeCloud::buildCellReferralLists()
+{
+    Info<< nl << "Determining molecule referring schedule" << endl;
+
+    const referredCellList& refIntL(referredInteractionList());
+
+    DynamicList<label> referralProcs;
+
+    // Run through all referredCells to build list of interacting processors
+
+    forAll(refIntL, rIL)
+    {
+        const referredCell& rC(refIntL[rIL]);
+
+        if (findIndex(referralProcs, rC.sourceProc()) == -1)
+        {
+            referralProcs.append(rC.sourceProc());
+        }
+    }
+
+    referralProcs.shrink();
+
+//     Pout << "referralProcs: " << nl << referralProcs << endl;
+
+    List<DynamicList<label> > cellSendingReferralLists(referralProcs.size());
+
+    List<DynamicList<DynamicList<label> > >
+        cellReceivingReferralLists(referralProcs.size());
+
+    // Run through all referredCells again building up send and receive info
+
+    forAll(refIntL, rIL)
+    {
+        const referredCell& rC(refIntL[rIL]);
+
+        label rPI = findIndex(referralProcs, rC.sourceProc());
+
+        DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
+
+        DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
+
+        label existingSource = findIndex(sRL, rC.sourceCell());
+
+        // Check to see if this source cell has already been allocated to
+        // come to this processor.  If not, add the source cell to the sending
+        // list and add the current referred cell to the receiving list.
+
+        // It shouldn't be possible for the sending and receiving lists to be
+        // different lengths, because their append operations happen at the
+        // same time.
+
+        if (existingSource == -1)
+        {
+            sRL.append(rC.sourceCell());
+
+            rRL.append
+            (
+                DynamicList<label> (labelList(1,rIL))
+            );
+        }
+        else
+        {
+            rRL[existingSource].append(rIL);
+
+            rRL[existingSource].shrink();
+        }
+    }
+
+    forAll(referralProcs, rPI)
+    {
+        DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
+
+        DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
+
+        sRL.shrink();
+
+        rRL.shrink();
+    }
+
+    // It is assumed that cell exchange is reciprocal, if proc A has cells to
+    // send to proc B, then proc B must have some to send to proc A.
+
+    cellReceivingReferralLists_.setSize(referralProcs.size());
+
+    cellSendingReferralLists_.setSize(referralProcs.size());
+
+    forAll(referralProcs, rPI)
+    {
+        DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
+
+        labelListList translLL(rRL.size());
+
+        forAll(rRL, rRLI)
+        {
+            translLL[rRLI] = rRL[rRLI];
+        }
+
+        cellReceivingReferralLists_[rPI] = receivingReferralList
+        (
+            referralProcs[rPI],
+            translLL
+        );
+    }
+
+    // Send sendingReferralLists to each interacting processor.
+
+    forAll(referralProcs, rPI)
+    {
+
+        DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
+
+        if (referralProcs[rPI] != Pstream::myProcNo())
+        {
+            if (Pstream::parRun())
+            {
+                OPstream toInteractingProc
+                (
+                    Pstream::blocking,
+                    referralProcs[rPI]
+                );
+
+                toInteractingProc << sendingReferralList
+                (
+                    Pstream::myProcNo(),
+                    sRL
+                );
+            }
+        }
+    }
+
+    // Receive sendingReferralLists from each interacting processor.
+
+    forAll(referralProcs, rPI)
+    {
+        if (referralProcs[rPI] != Pstream::myProcNo())
+        {
+            if (Pstream::parRun())
+            {
+                IPstream fromInteractingProc
+                (
+                    Pstream::blocking,
+                    referralProcs[rPI]
+                );
+
+                fromInteractingProc >> cellSendingReferralLists_[rPI];
+            }
+        }
+        else
+        {
+            cellSendingReferralLists_[rPI] = sendingReferralList
+            (
+                Pstream::myProcNo(),
+                cellSendingReferralLists[rPI]
+            );
+        }
+    }
+
+//     Pout << "receiving list: " << nl << cellReceivingReferralLists_ << endl;
+
+//     Pout << "sending list: " << nl << cellSendingReferralLists_ << endl;
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/moleculeCloud/moleculeCloudBuildCellsFoundInRange.H

@@ -0,0 +1,52 @@
+forAll (pointsOnThisSegment, pS)
+{
+    const point& ptP = mesh_.points()[pointsOnThisSegment[pS]];
+
+    // Assessing real cells in range is only required on the 1st iteration
+    // because they do not change from iteration to iteration.
+
+    if (iterationNo == 0)
+    {
+        forAll (mesh_.points(), pointIIndex)
+        {
+
+            // Compare separation of ptP to all other points in the mesh,
+            // add unique reference to cell with any point within rCutMax_
+            // to realCellsFoundInRange.
+
+            const point& ptI(mesh_.points()[pointIIndex]);
+
+            if (magSqr(ptP - ptI) <= rCutMaxSqr)
+            {
+                const labelList& ptICells(mesh_.pointCells()[pointIIndex]);
+
+                forAll(ptICells, pIC)
+                {
+                    const label cellI(ptICells[pIC]);
+
+                    if(findIndex(realCellsFoundInRange, cellI) == -1)
+                    {
+                        realCellsFoundInRange.append(cellI);
+                    }
+                }
+            }
+        }
+    }
+
+    forAll(referredInteractionList, rIL)
+    {
+        const vectorList& refCellPoints =
+            referredInteractionList[rIL].vertexPositions();
+
+        forAll(refCellPoints, rCP)
+        {
+            if (magSqr(ptP - refCellPoints[rCP]) <= rCutMaxSqr)
+            {
+                if(findIndex(referredCellsFoundInRange, rIL) == -1)
+                {
+                    referredCellsFoundInRange.append(rIL);
+                }
+            }
+        }
+    }
+}

src/lagrangian/molecule/moleculeCloud/moleculeCloudBuildDirectInteractionList.H

@@ -0,0 +1,188 @@
+Info << nl << "Building list of direct interaction neighbours" << endl;
+
+forAll (mesh_.points(), p)
+{
+    forAll(mesh_.faces(), f)
+    {
+        if(testPointFaceDistance(p, f))
+        {
+            const labelList& pCells(mesh_.pointCells()[p]);
+
+            const label cellO(mesh_.faceOwner()[f]);
+
+            const label cellN(mesh_.faceNeighbour()[f]);
+
+            forAll(pCells, pC)
+            {
+                const label cellI(pCells[pC]);
+
+                // cells are not added to their own DIL
+
+                if (cellO > cellI)
+                {
+                    if (findIndex(directInteractionList[cellI], cellO) == -1)
+                    {
+                        directInteractionList[cellI].append(cellO);
+                    }
+                }
+
+                if (cellI > cellO)
+                {
+                    if (findIndex(directInteractionList[cellO], cellI) == -1)
+                    {
+                        directInteractionList[cellO].append(cellI);
+                    }
+                }
+
+                if (mesh_.isInternalFace(f))
+                {
+                    // boundary faces will not have neighbour information
+
+                    if (cellN > cellI)
+                    {
+                        if
+                        (
+                             findIndex(directInteractionList[cellI], cellN)
+                          == -1
+                        )
+                        {
+                            directInteractionList[cellI].append(cellN);
+                        }
+                    }
+
+                    if (cellI > cellN)
+                    {
+                        if
+                        (
+                            findIndex(directInteractionList[cellN], cellI)
+                         == -1
+                        )
+                        {
+                            directInteractionList[cellN].append(cellI);
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+label edgeJIndex;
+
+forAll (mesh_.edges(), edgeIIndex)
+{
+    const edge& eI(mesh_.edges()[edgeIIndex]);
+
+    for
+    (
+        edgeJIndex = edgeIIndex + 1;
+        edgeJIndex != mesh_.edges().size();
+        ++edgeJIndex
+    )
+    {
+        const edge& eJ(mesh_.edges()[edgeJIndex]);
+
+        if (testEdgeEdgeDistance(eI, eJ))
+        {
+            const labelList& eICells(mesh_.edgeCells()[edgeIIndex]);
+
+            const labelList& eJCells(mesh_.edgeCells()[edgeJIndex]);
+
+            forAll(eICells, eIC)
+            {
+                const label cellI(eICells[eIC]);
+
+                forAll(eJCells, eJC)
+                {
+                    const label cellJ(eJCells[eJC]);
+
+                    if (cellJ > cellI)
+                    {
+                        if
+                        (
+                            findIndex(directInteractionList[cellI], cellJ)
+                         == -1
+                        )
+                        {
+                            directInteractionList[cellI].append(cellJ);
+                        }
+                    }
+
+                    if (cellI > cellJ)
+                    {
+                        if
+                        (
+                            findIndex(directInteractionList[cellJ], cellI)
+                         == -1
+                        )
+                        {
+                            directInteractionList[cellJ].append(cellI);
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+// label pointJIndex;
+//
+// forAll (mesh_.points(), pointIIndex)
+// {
+//     const point& ptI
+//     (
+//         mesh_.points()[pointIIndex]
+//     );
+//
+//     for
+//     (
+//         pointJIndex = pointIIndex;
+//         pointJIndex != mesh_.points().size();
+//         ++pointJIndex
+//     )
+//     {
+//         const point& ptJ
+//         (
+//             mesh_.points()[pointJIndex]
+//         );
+//
+//         if (magSqr(ptI - ptJ) <= rCutMaxSqr)
+//         {
+//             const labelList& ptICells
+//             (
+//                 mesh_.pointCells()[pointIIndex]
+//             );
+//
+//             const labelList& ptJCells
+//             (
+//                 mesh_.pointCells()[pointJIndex]
+//             );
+//
+//             forAll(ptICells, pIC)
+//             {
+//                 const label cellI(ptICells[pIC]);
+//
+//                 forAll(ptJCells, pJC)
+//                 {
+//                     const label cellJ(ptJCells[pJC]);
+//
+//                     if (cellJ > cellI)
+//                     {
+//                         if(findIndex(directInteractionList[cellI], cellJ) == -1)
+//                         {
+//                             directInteractionList[cellI].append(cellJ);
+//                         }
+//                     }
+//
+//                     if (cellI > cellJ)
+//                     {
+//                         if(findIndex(directInteractionList[cellJ], cellI) == -1)
+//                         {
+//                             directInteractionList[cellJ].append(cellI);
+//                         }
+//                     }
+//                 }
+//             }
+//         }
+//     }
+// }

src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculateExternalForce.C

@@ -0,0 +1,42 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "moleculeCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::moleculeCloud::calculateExternalForce()
+{
+    iterator mol(this->begin());
+
+    for (mol = this->begin(); mol != this->end(); ++mol)
+    {
+        mol().A() += gravity_;
+    }
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculateForce.C

@@ -0,0 +1,59 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "moleculeCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::moleculeCloud::calculateForce()
+{
+    buildCellOccupancy();
+
+    iterator mol(this->begin());
+
+    // Set all accumulated quantities to zero
+    for (mol = this->begin(); mol != this->end(); ++mol)
+    {
+        mol().A() = vector::zero;
+
+        mol().potentialEnergy() = 0.0;
+
+        mol().rf() = tensor::zero;
+    }
+
+    calculatePairForce();
+
+    calculateTetherForce();
+
+    calculateExternalForce();
+}
+
+
+// ************************************************************************* //
+
+
+
+

src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculatePairForce.C

@@ -0,0 +1,41 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "moleculeCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::moleculeCloud::calculatePairForce()
+{
+    iterator mol(this->begin());
+
+#   include "moleculeCloudCalculatePairForceRealCells.H"
+
+#   include "moleculeCloudCalculatePairForceReferredCells.H"
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculatePairForceRealCells.H

@@ -0,0 +1,49 @@
+
+vector rIJ;
+
+scalar rIJMag;
+
+scalar rIJMagSq;
+
+vector fIJ;
+
+label idI;
+
+label idJ;
+
+mol = this->begin();
+
+molecule* molI = &mol();
+
+molecule* molJ = &mol();
+
+forAll(directInteractionList_, dIL)
+{
+    forAll(cellOccupancy_[dIL],cellIMols)
+    {
+        molI = cellOccupancy_[dIL][cellIMols];
+
+        forAll(directInteractionList_[dIL], interactingCells)
+        {
+            List< molecule* > cellJ =
+                cellOccupancy_[directInteractionList_[dIL][interactingCells]];
+
+            forAll(cellJ, cellJMols)
+            {
+                molJ = cellJ[cellJMols];
+
+#               include "moleculeCloudCalculatePairForceRealCellsCalculationStep.H"
+            }
+        }
+
+        forAll(cellOccupancy_[dIL],cellIOtherMols)
+        {
+            molJ = cellOccupancy_[dIL][cellIOtherMols];
+
+            if (molJ > molI)
+            {
+#                include "moleculeCloudCalculatePairForceRealCellsCalculationStep.H"
+            }
+        }
+    }
+}

src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculatePairForceRealCellsCalculationStep.H

@@ -0,0 +1,33 @@
+idI = molI->id();
+
+idJ = molJ->id();
+
+rIJ = molI->position() - molJ->position();
+
+rIJMagSq = magSqr(rIJ);
+
+if (pairPotentials_.rCutSqr(idI, idJ, rIJMagSq))
+{
+    rIJMag = mag(rIJ);
+
+    fIJ = (rIJ/rIJMag)*pairPotentials_.force(idI, idJ, rIJMag);
+
+    // Acceleration increment for molI
+    molI->A() += fIJ/(molI->mass());
+
+    // Acceleration increment for molJ
+    molJ->A() += -fIJ/(molJ->mass());
+
+    scalar potentialEnergy
+    (
+        pairPotentials_.energy(idI, idJ, rIJMag)
+    );
+
+    molI->potentialEnergy() += 0.5*potentialEnergy;
+
+    molJ->potentialEnergy() += 0.5*potentialEnergy;
+
+    molI->rf() += rIJ * fIJ;
+
+    molJ->rf() += rIJ * fIJ;
+}

src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculatePairForceReferredCells.H

@@ -0,0 +1,59 @@
+vector rKL;
+
+scalar rKLMag;
+
+scalar rKLMagSq;
+
+vector fKL;
+
+label idK;
+
+label idL;
+
+molecule* molK = &mol();
+
+forAll(referredInteractionList_, rIL)
+{
+    const List<label>& realCells =
+        referredInteractionList_[rIL].realCellsForInteraction();
+
+    forAll(referredInteractionList_[rIL], refMols)
+    {
+        referredMolecule* molL = &(referredInteractionList_[rIL][refMols]);
+
+        forAll(realCells, rC)
+        {
+            List<molecule*> cellK = cellOccupancy_[realCells[rC]];
+
+            forAll(cellK, cellKMols)
+            {
+                molK = cellK[cellKMols];
+
+                idK = molK->id();
+
+                idL = molL->id();
+
+                rKL = molK->position() - molL->position();
+
+                rKLMagSq = magSqr(rKL);
+
+                if (pairPotentials_.rCutSqr(idK, idL, rKLMagSq))
+                {
+                    rKLMag = mag(rKL);
+
+                    fKL = (rKL/rKLMag)*pairPotentials_.force(idK, idL, rKLMag);
+
+                    // Acceleration increment for molK
+                    molK->A() += fKL/(molK->mass());
+
+                    // Adding a contribution of 1/2 of the potential energy
+                    // from this interaction
+                    molK->potentialEnergy() +=
+                        0.5*pairPotentials_.energy(idK, idL, rKLMag);
+
+                    molK->rf() += rKL * fKL;
+                }
+            }
+        }
+    }
+}

src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculateTetherForce.C

@@ -0,0 +1,69 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "moleculeCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::moleculeCloud::calculateTetherForce()
+{
+    iterator mol(this->begin());
+
+    vector rIT;
+
+    scalar rITMag;
+
+    vector fIT;
+
+    for (mol = this->begin(); mol != this->end(); ++mol)
+    {
+        if (mol().tethered())
+        {
+            rIT = mol().position() - mol().tetherPosition();
+
+            rITMag = mag(rIT);
+
+            fIT = (rIT/rITMag) * tetherPotentials_.force
+            (
+                mol().id(),
+                rITMag
+            );
+
+            mol().A() += fIT/(mol().mass());
+
+            mol().potentialEnergy() += tetherPotentials_.energy
+            (
+                mol().id(),
+                rITMag
+            );
+
+            mol().rf() += rIT*fIT;
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/moleculeCloud/moleculeCloudI.H

@@ -0,0 +1,184 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 polyMesh& moleculeCloud::mesh() const
+{
+    return mesh_;
+}
+
+
+// MD solution parameters
+
+inline const moleculeCloud::integrationMethods&
+    moleculeCloud::integrationMethod() const
+{
+    return integrationMethod_;
+}
+
+
+inline scalar moleculeCloud::potentialEnergyLimit() const
+{
+    return potentialEnergyLimit_;
+}
+
+
+inline scalar moleculeCloud::guardRadius() const
+{
+    return guardRadius_;
+}
+
+
+inline const scalar moleculeCloud::rCutMax() const
+{
+    return rCutMax_;
+}
+
+
+inline const scalar moleculeCloud::rCutMaxSqr() const
+{
+    return rCutMaxSqr_;
+}
+
+
+inline const List<word>& moleculeCloud::idList() const
+{
+    return idList_;
+}
+
+
+inline label moleculeCloud::nPairPotentials() const
+{
+    return pairPotentials_.size();
+}
+
+
+inline label moleculeCloud::nIds() const
+{
+    return idList_.size();
+}
+
+
+inline const labelList& moleculeCloud::removalOrder() const
+{
+    return removalOrder_;
+}
+
+
+inline const labelListList& moleculeCloud::directInteractionList() const
+{
+    return directInteractionList_;
+}
+
+
+inline const referredCellList& moleculeCloud::referredInteractionList() const
+{
+    return referredInteractionList_;
+}
+
+
+inline const labelList&
+    moleculeCloud::realCellsWithinRCutMaxOfAnyReferringPatch() const
+{
+    return realCellsWithinRCutMaxOfAnyReferringPatch_;
+}
+
+
+inline const labelList&
+    moleculeCloud::realFacesWithinRCutMaxOfAnyReferringPatch() const
+{
+    return realFacesWithinRCutMaxOfAnyReferringPatch_;
+}
+
+
+inline const labelList&
+    moleculeCloud::realEdgesWithinRCutMaxOfAnyReferringPatch() const
+{
+    return realEdgesWithinRCutMaxOfAnyReferringPatch_;
+}
+
+
+inline const labelList&
+    moleculeCloud::realPointsWithinRCutMaxOfAnyReferringPatch() const
+{
+    return realPointsWithinRCutMaxOfAnyReferringPatch_;
+}
+
+
+inline const List<sendingReferralList>&
+    moleculeCloud::cellSendingReferralLists() const
+{
+    return cellSendingReferralLists_;
+}
+
+
+inline const List<receivingReferralList>&
+    moleculeCloud::cellReceivingReferralLists() const
+{
+    return cellReceivingReferralLists_;
+}
+
+
+inline label moleculeCloud::nInteractingProcs() const
+{
+    return cellReceivingReferralLists_.size();
+}
+
+
+inline const pairPotentialList& moleculeCloud::pairPotentials() const
+{
+    return pairPotentials_;
+}
+
+
+inline const tetherPotentialList& moleculeCloud::tetherPotentials() const
+{
+    return tetherPotentials_;
+}
+
+
+inline const vector& moleculeCloud::gravity() const
+{
+    return gravity_;
+}
+
+
+inline const List< DynamicList<molecule*> >&
+    moleculeCloud::cellOccupancy() const
+{
+    return cellOccupancy_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+}  // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/molecule/moleculeCloud/moleculeCloudIntegrateEquationsOfMotion.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 "moleculeCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::moleculeCloud::integrateEquationsOfMotion()
+{
+    // Supply trackData to the move method to tell it whether this is the
+    // leapfrog stage 1 or stage 2
+    // or whether it is the
+    // predictor or corrector step.
+
+    molecule::trackData td1(*this, 1);
+
+    Cloud<molecule>::move(td1);
+
+    calculateForce();
+
+    molecule::trackData td2(*this, 2);
+
+    Cloud<molecule>::move(td2);
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/moleculeCloud/moleculeCloudReadMDParameters/moleculeCloudReadIdList.H

@@ -0,0 +1,14 @@
+IOdictionary idListDict
+(
+    IOobject
+    (
+        "idList",
+        mesh_.time().constant(),
+        mesh_,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+idList_ = List<word>(idListDict.lookup("idList"));
+

src/lagrangian/molecule/moleculeCloud/moleculeCloudReadMDParameters/moleculeCloudReadMDParameters.H

@@ -0,0 +1,6 @@
+
+#    include "moleculeCloudReadMDSolution.H"
+
+#    include "moleculeCloudReadIdList.H"
+
+#    include "moleculeCloudReadPotentials.H"

src/lagrangian/molecule/moleculeCloud/moleculeCloudReadMDParameters/moleculeCloudReadMDSolution.H

@@ -0,0 +1,27 @@
+Info<< nl <<  "Reading MD solution parameters:" << endl;
+
+IOdictionary mdSolution
+(
+    IOobject
+    (
+        "mdSolution",
+        mesh_.time().system(),
+        mesh_,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+integrationMethod_ = integrationMethodNames_.read
+(
+    mdSolution.lookup("integrationMethod")
+);
+Info<< "integrationMethod =\t\t\t"
+    << integrationMethodNames_[integrationMethod_] << endl;
+
+potentialEnergyLimit_ = readScalar(mdSolution.lookup("potentialEnergyLimit"));
+Info<< "potentialEnergyLimit =\t\t\t" << potentialEnergyLimit_ << endl;
+
+guardRadius_ = readScalar(mdSolution.lookup("guardRadius"));
+Info<< "guardRadius =\t\t\t\t" << guardRadius_ << endl;
+

src/lagrangian/molecule/moleculeCloud/moleculeCloudReadMDParameters/moleculeCloudReadPotentials.H

@@ -0,0 +1,327 @@
+IOdictionary potentialDict
+(
+    IOobject
+    (
+        "potentialDict",
+        mesh_.time().system(),
+        mesh_,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+// ****************************************************************************
+// Pair potentials
+
+Info << nl << "Creating pair potentials:" << nl << endl;
+
+if (!potentialDict.found("pair"))
+{
+    FatalErrorIn("moleculeCloudReadPotentials.H") << nl
+            << "pair potential specification subDict not found"
+            << abort(FatalError);
+}
+
+const dictionary& pairDict = potentialDict.subDict("pair");
+
+pairPotentials_.setSizeAndNIds(nIds());
+
+label a;
+
+label b;
+
+rCutMax_ = 0.0;
+
+for
+(
+    a = 0;
+    a < idList_.size();
+    ++a
+)
+{
+    word idA = idList_[a];
+
+    for
+    (
+        b = a;
+        b < idList_.size();
+        ++b
+    )
+    {
+        word idB = idList_[b];
+
+        word pairPotentialName;
+
+        if (a==b)
+        {
+
+
+            if(pairDict.found(idA+"-"+idB))
+            {
+                pairPotentialName = idA+"-"+idB;
+            }
+
+            else
+            {
+                FatalErrorIn("moleculeCloudReadPotentials.H") << nl
+                        << "Pair pairPotential specification subDict "
+                        << idA+"-"+idB << " not found"
+                        << abort(FatalError);
+            }
+        }
+
+        else
+        {
+            if(pairDict.found(idA+"-"+idB))
+            {
+                pairPotentialName = idA+"-"+idB;
+            }
+
+            else if(pairDict.found(idB+"-"+idA))
+            {
+                pairPotentialName = idB+"-"+idA;
+            }
+
+            else
+            {
+                FatalErrorIn("moleculeCloudReadPotentials.H") << nl
+                        << "Pair pairPotential specification subDict "
+                        << idA+"-"+idB << " or " << idB+"-"+idA << " not found"
+                        << abort(FatalError);
+            }
+
+            if
+            (
+                pairDict.found(idA+"-"+idB)
+                && pairDict.found(idB+"-"+idA)
+            )
+            {
+                FatalErrorIn("moleculeCloudReadPotentials.H") << nl
+                        << "Pair pairPotential specification subDict "
+                        << idA+"-"+idB << " and "
+                        << idB+"-"+idA << " found - multiple definition"
+                        << abort(FatalError);
+            }
+        }
+
+        word pairPotentialType
+        (
+            pairDict.subDict(pairPotentialName).lookup("potentialType")
+        );
+
+        if
+        (
+            pairPotentialType == "maitlandSmithTabulated"
+        )
+        {
+            scalar rCut = readScalar
+            (
+                pairDict.subDict(pairPotentialName).lookup("rCut")
+            );
+
+            if(rCut > rCutMax_)
+            {
+                rCutMax_ = rCut;
+            }
+
+            pairPotentials_.addPotential
+            (
+                a,
+                b,
+                pairPotential
+                (
+                    pairDict.subDict(pairPotentialName)
+                )
+            );
+        }
+        else
+        {
+            FatalErrorIn("moleculeCloudReadPotentials.H") << nl
+                    << "pairPotentialType "<< pairPotentialType << " unknown"
+                    << abort(FatalError);
+        }
+
+        Info << pairPotentialName << ": "
+            << pairPotentials_.pairPotentialFunction(a,b) << endl;
+    }
+}
+
+rCutMax_ += guardRadius_;
+
+rCutMaxSqr_ = rCutMax_ * rCutMax_;
+
+Info << nl << "rCutMax + guardRadius = " << rCutMax_ << endl;
+
+if (potentialDict.found("removalOrder"))
+{
+    List<word> remOrd = potentialDict.lookup("removalOrder");
+
+    removalOrder_.setSize(remOrd.size());
+
+    forAll(removalOrder_, rO)
+    {
+        removalOrder_[rO] = findIndex(idList_, remOrd[rO]);
+    }
+}
+
+// ****************************************************************************
+// Tether potentials
+
+iterator mol(this->begin());
+
+DynamicList<label> tetherIds;
+
+for
+(
+ mol = this->begin();
+ mol != this->end();
+ ++mol
+)
+{
+    if (mol().tethered())
+    {
+        if (findIndex(tetherIds, mol().id()) == -1)
+        {
+            tetherIds.append
+            (
+                mol().id()
+            );
+        }
+    }
+}
+
+if (Pstream::parRun())
+{
+    List< labelList > allTetherIds(Pstream::nProcs());
+
+    allTetherIds[Pstream::myProcNo()] = tetherIds;
+
+    Pstream::gatherList(allTetherIds);
+
+    if (Pstream::master())
+    {
+        DynamicList<label> globalTetherIds;
+
+        forAll(allTetherIds, procN)
+        {
+            const labelList& procNTetherIds = allTetherIds[procN];
+
+            forAll(procNTetherIds, id)
+            {
+                if (findIndex(globalTetherIds, procNTetherIds[id]) == -1)
+                {
+                    globalTetherIds.append
+                    (
+                        procNTetherIds[id]
+                    );
+                }
+            }
+        }
+
+        globalTetherIds.shrink();
+
+        tetherIds = globalTetherIds;
+    }
+
+    Pstream::scatter(tetherIds);
+}
+
+tetherIds.shrink();
+
+if (tetherIds.size())
+{
+    Info << nl << "Creating tether potentials:" << endl;
+
+    if (!potentialDict.found("tether"))
+    {
+        FatalErrorIn("moleculeCloudReadPotentials.H") << nl
+                << "tether potential specification subDict not found"
+                << abort(FatalError);
+    }
+
+    const dictionary& tetherDict = potentialDict.subDict("tether");
+
+    tetherPotentials_.setSizeAndTetherIds(tetherIds);
+
+    forAll (tetherIds, tid)
+    {
+        word tetherPotentialName = idList_[tetherIds[tid]];
+
+        if (!tetherDict.found(tetherPotentialName))
+        {
+            FatalErrorIn("moleculeCloudReadPotentials.H") << nl
+                    << "tether potential specification subDict "
+                    << tetherPotentialName << " not found"
+                    << abort(FatalError);
+        }
+
+        Info << nl << tetherPotentialName << endl;
+
+        word tetherPotentialType =
+            tetherDict.subDict(tetherPotentialName).lookup("potentialType");
+
+        if
+        (
+            tetherPotentialType == "harmonicSpring"
+        )
+        {
+
+            Info << "\tpotentialType = " << tetherPotentialType << endl;
+
+            scalar springConstant
+            (
+                readScalar
+                (
+                    tetherDict.subDict(tetherPotentialName).lookup
+                    (
+                        "springConstant"
+                    )
+                )
+            );
+
+            Info << "\tspringConstant = " << springConstant << endl;
+
+            tetherPotentials_.addPotential
+            (
+                tetherIds[tid],
+                tetherPotential
+                (
+                    tetherPotentialName,
+                    tetherPotentialType,
+                    springConstant
+                )
+            );
+        }
+
+        else
+        {
+            FatalErrorIn("moleculeCloudReadPotentials.H") << nl
+                    << "tetherPotentialType "
+                    << tetherPotentialType << " unknown"
+                    << abort(FatalError);
+        }
+    }
+}
+
+// ****************************************************************************
+// External Forces
+
+gravity_ = vector::zero;
+
+if (potentialDict.found("external"))
+{
+
+    Info << nl << "Reading external forces:" << endl;
+
+    const dictionary& externalDict = potentialDict.subDict("external");
+
+    // ************************************************************************
+    // gravity
+
+    if (externalDict.found("gravity"))
+    {
+        gravity_ = externalDict.lookup("gravity");
+    }
+}
+
+Info << nl << tab << "gravity = " << gravity_ << endl;

src/lagrangian/molecule/moleculeCloud/moleculeCloudRealCellsInRangeOfSegment.C

@@ -0,0 +1,152 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "moleculeCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+const Foam::labelList Foam::moleculeCloud::realCellsInRangeOfSegment
+(
+    const labelList& segmentFaces,
+    const labelList& segmentEdges,
+    const labelList& segmentPoints
+) const
+{
+    DynamicList<label> realCellsFoundInRange;
+
+    forAll(segmentFaces, sF)
+    {
+        const label f = segmentFaces[sF];
+
+        forAll (mesh_.points(), p)
+        {
+            if (testPointFaceDistance(p, f))
+            {
+                const labelList& pCells(mesh_.pointCells()[p]);
+
+                forAll(pCells, pC)
+                {
+                    const label cellI(pCells[pC]);
+
+                    if (findIndex(realCellsFoundInRange, cellI) == -1)
+                    {
+                        realCellsFoundInRange.append(cellI);
+                    }
+                }
+            }
+        }
+    }
+
+    forAll(segmentPoints, sP)
+    {
+        const label p = segmentPoints[sP];
+
+        forAll(mesh_.faces(), f)
+        {
+            if (testPointFaceDistance(p, f))
+            {
+                const label cellO(mesh_.faceOwner()[f]);
+
+                if (findIndex(realCellsFoundInRange, cellO) == -1)
+                {
+                    realCellsFoundInRange.append(cellO);
+                }
+
+                if (mesh_.isInternalFace(f))
+                {
+                    // boundary faces will not have neighbour information
+
+                    const label cellN(mesh_.faceNeighbour()[f]);
+
+                    if (findIndex(realCellsFoundInRange, cellN) == -1)
+                    {
+                        realCellsFoundInRange.append(cellN);
+                    }
+                }
+            }
+        }
+    }
+
+    forAll(segmentEdges, sE)
+    {
+        const edge& eJ(mesh_.edges()[segmentEdges[sE]]);
+
+        forAll (mesh_.edges(), edgeIIndex)
+        {
+            const edge& eI(mesh_.edges()[edgeIIndex]);
+
+            if (testEdgeEdgeDistance(eI, eJ))
+            {
+                const labelList& eICells(mesh_.edgeCells()[edgeIIndex]);
+
+                forAll(eICells, eIC)
+                {
+                    const label cellI(eICells[eIC]);
+
+                    if (findIndex(realCellsFoundInRange, cellI) == -1)
+                    {
+                        realCellsFoundInRange.append(cellI);
+                    }
+                }
+            }
+        }
+    }
+
+//     forAll (points, p)
+//     {
+//         const point& ptI = mesh_.points()[points[p]];
+//
+//         forAll(mesh_.faces(), f)
+//         {
+//             if (testPointFaceDistance(ptI, f))
+//             {
+//                 const label cellO(mesh_.faceOwner()[f]);
+//
+//                 if (findIndex(realCellsFoundInRange, cellO) == -1)
+//                 {
+//                     realCellsFoundInRange.append(cellO);
+//                 }
+//
+//                 if (mesh_.isInternalFace(f))
+//                 {
+//                     // boundary faces will not have neighbour information
+//
+//                     const label cellN(mesh_.faceNeighbour()[f]);
+//
+//                     if (findIndex(realCellsFoundInRange, cellN) == -1)
+//                     {
+//                         realCellsFoundInRange.append(cellN);
+//                     }
+//                 }
+//             }
+//         }
+//     }
+//
+    return realCellsFoundInRange.shrink();
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/moleculeCloud/moleculeCloudReferredCellsInRangeOfSegment.C

@@ -0,0 +1,117 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "moleculeCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+const Foam::labelList Foam::moleculeCloud::referredCellsInRangeOfSegment
+(
+    const List<referredCell>& referredInteractionList,
+    const labelList& segmentFaces,
+    const labelList& segmentEdges,
+    const labelList& segmentPoints
+) const
+{
+    DynamicList<label> referredCellsFoundInRange;
+
+    forAll(segmentFaces, sF)
+    {
+        const label f = segmentFaces[sF];
+
+        forAll(referredInteractionList, rIL)
+        {
+            const vectorList& refCellPoints
+                = referredInteractionList[rIL].vertexPositions();
+
+            if (testPointFaceDistance(refCellPoints, f))
+            {
+                if (findIndex(referredCellsFoundInRange, rIL) == -1)
+                {
+                    referredCellsFoundInRange.append(rIL);
+                }
+            }
+        }
+    }
+
+    forAll(segmentPoints, sP)
+    {
+        const label p = segmentPoints[sP];
+
+        forAll(referredInteractionList, rIL)
+        {
+            const referredCell& refCell(referredInteractionList[rIL]);
+
+            if (testPointFaceDistance(p, refCell))
+            {
+                if (findIndex(referredCellsFoundInRange, rIL) == -1)
+                {
+                    referredCellsFoundInRange.append(rIL);
+                }
+            }
+        }
+    }
+
+    forAll(segmentEdges, sE)
+    {
+        const edge& eI(mesh_.edges()[segmentEdges[sE]]);
+
+        forAll(referredInteractionList, rIL)
+        {
+            const vectorList& refCellPoints
+                = referredInteractionList[rIL].vertexPositions();
+
+            const edgeList& refCellEdges
+                = referredInteractionList[rIL].edges();
+
+            forAll(refCellEdges, rCE)
+            {
+                const edge& eJ(refCellEdges[rCE]);
+
+                if
+                (
+                    testEdgeEdgeDistance
+                    (
+                        eI,
+                        refCellPoints[eJ.start()],
+                        refCellPoints[eJ.end()]
+                    )
+                )
+                {
+                    if(findIndex(referredCellsFoundInRange, rIL) == -1)
+                    {
+                        referredCellsFoundInRange.append(rIL);
+                    }
+                }
+            }
+        }
+    }
+
+    return referredCellsFoundInRange.shrink();
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlaps.C

@@ -0,0 +1,67 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "moleculeCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::moleculeCloud::removeHighEnergyOverlaps()
+{
+    Info << nl << "Removing high energy overlaps, removal order:";
+
+    forAll(removalOrder_, rO)
+    {
+        Info << " " << idList_[removalOrder_[rO]];
+    }
+
+    Info << nl ;
+
+    label initialSize = this->size();
+
+    buildCellOccupancy();
+
+    iterator mol(this->begin());
+
+#   include "moleculeCloudRemoveHighEnergyOverlapsRealCells.H"
+
+    buildCellOccupancy();
+
+#   include "moleculeCloudRemoveHighEnergyOverlapsReferredCells.H"
+
+    buildCellOccupancy();
+
+    label molsRemoved = initialSize - this->size();
+
+    if (Pstream::parRun())
+    {
+        reduce(molsRemoved, sumOp<label>());
+    }
+
+    Info << tab << molsRemoved << " molecules removed" << endl;
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlapsRealCells.H

@@ -0,0 +1,55 @@
+{ 
+    vector rIJ;
+
+    scalar rIJMag;
+
+    scalar rIJMagSq;
+
+    label idI;
+
+    label idJ;
+
+    mol = this->begin();
+
+    molecule* molI = &mol();
+
+    molecule* molJ = &mol();
+
+    DynamicList<molecule*> molsToDelete;
+
+    forAll(directInteractionList_, dIL)
+    {
+        forAll(cellOccupancy_[dIL],cellIMols)
+        {
+            molI = cellOccupancy_[dIL][cellIMols];
+
+            forAll(directInteractionList_[dIL], interactingCells)
+            {
+                List< molecule* > cellJ =
+                    cellOccupancy_[directInteractionList_[dIL][interactingCells]];
+
+                forAll(cellJ, cellJMols)
+                {
+                    molJ = cellJ[cellJMols];
+
+#                   include "moleculeCloudRemoveHighEnergyOverlapsRealCellsCalculationStep.H"
+                }
+            }
+
+            forAll(cellOccupancy_[dIL],cellIOtherMols)
+            {
+                molJ = cellOccupancy_[dIL][cellIOtherMols];
+
+                if (molJ > molI)
+                {
+#                   include "moleculeCloudRemoveHighEnergyOverlapsRealCellsCalculationStep.H"
+                }
+            }
+        }
+    }
+
+    forAll (molsToDelete, mTD)
+    {
+        deleteParticle(*(molsToDelete[mTD]));
+    }
+}

src/lagrangian/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlapsRealCellsCalculationStep.H

@@ -0,0 +1,70 @@
+idI = molI->id();
+
+idJ = molJ->id();
+
+rIJ = molI->position() - molJ->position();
+
+rIJMagSq = magSqr(rIJ);
+
+if (pairPotentials_.rCutSqr(idI, idJ, rIJMagSq))
+{
+    rIJMag = mag(rIJ);
+
+    bool remove = false;
+
+    if (rIJMag < SMALL)
+    {
+        WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
+            << "Real molecule pair "
+            << " idI = " << idI
+            << " at position " << molI->position()
+            << " idJ = " << idJ
+            << " at position " << molJ->position()
+            << " are closer than " << SMALL
+            << ": mag separation = " << rIJMag
+            << ". These may have been placed on top of each"
+            << " other by a rounding error in molConfig in parallel"
+            << " or a block filled with molecules twice."
+            << " Removing one of the molecules."
+            << endl;
+
+        remove = true;
+    }
+
+    // Guard against pairPotentials_.energy being evaluated
+    // if rIJMag < SMALL. A floating point exception will
+    // happen otherwise.
+
+    if (!remove)
+    {
+        if
+        (
+            pairPotentials_.energy(idI, idJ, rIJMag) > potentialEnergyLimit_
+        )
+        {
+            remove = true;
+        }
+    }
+
+    if (remove)
+    {
+        if
+        (
+            idI == idJ
+         || findIndex(removalOrder_, idJ) < findIndex(removalOrder_, idI)
+        )
+        {
+            if (findIndex(molsToDelete, molJ) == -1)
+            {
+                molsToDelete.append(molJ);
+            }
+        }
+        else
+        {
+            if (findIndex(molsToDelete, molI) == -1)
+            {
+                molsToDelete.append(molI);
+            }
+        }
+    }
+}

src/lagrangian/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlapsReferredCells.H

@@ -0,0 +1,144 @@
+{
+    vector rKL;
+
+    scalar rKLMag;
+
+    scalar rKLMagSq;
+
+    label idK;
+
+    label idL;
+
+    molecule* molK = &mol();
+
+    DynamicList<molecule*> molsToDelete;
+
+    forAll(referredInteractionList_, rIL)
+    {
+        referredCell& refCell = referredInteractionList_[rIL];
+
+        forAll(refCell, refMols)
+        {
+            referredMolecule* molL = &(refCell[refMols]);
+
+            List <label> realCells = refCell.realCellsForInteraction();
+
+            forAll(realCells, rC)
+            {
+                label cellK = realCells[rC];
+
+                List<molecule*> cellKMols = cellOccupancy_[cellK];
+
+                forAll(cellKMols, cKM)
+                {
+                    molK = cellKMols[cKM];
+
+                    idK = molK->id();
+
+                    idL = molL->id();
+
+                    rKL = molK->position() - molL->position();
+
+                    rKLMagSq = magSqr(rKL);
+
+                    if (pairPotentials_.rCutSqr(idK, idL, rKLMagSq))
+                    {
+                        rKLMag = mag(rKL);
+
+                        bool remove = false;
+
+                        if (rKLMag < SMALL)
+                        {
+                            WarningIn
+                            (
+                                "moleculeCloud::removeHighEnergyOverlaps()"
+                            )
+                                << "Real-referred molecule pair "
+                                << " idK = " << idK << " (real)"
+                                << " at position " << molK->position()
+                                << " idL = " << idL << " (referred)"
+                                << " at position " << molL->position()
+                                << " are closer than " << SMALL
+                                << ": mag separation = " << rKLMag
+                                << ". These may have been placed on top of each"
+                                << " other by a rounding error in molConfig in"
+                                << " parallel, a block filled with molecules"
+                                << " twice, or a lattice ending very close"
+                                << " to the edge of the mesh."
+                                << " Removing one of the molecules."
+                                << endl;
+
+                                remove = true;
+                        }
+
+                        // Guard against pairPotentials_.energy being evaluated
+                        // if rKLMag < SMALL. A floating point exception will
+                        // happen otherwise.
+
+                        if (!remove)
+                        {
+                            if
+                            (
+                                pairPotentials_.energy(idK, idL, rKLMag)
+                              > potentialEnergyLimit_
+                            )
+                            {
+                                remove = true;
+                            }
+                        }
+
+                        if (remove)
+                        {
+                            if
+                            (
+                                findIndex(removalOrder_, idK)
+                              < findIndex(removalOrder_, idL)
+                            )
+                            {
+                                if (findIndex(molsToDelete, molK) == -1)
+                                {
+                                    molsToDelete.append(molK);
+                                }
+                            }
+
+                            else if
+                            (
+                                findIndex(removalOrder_, idK)
+                             == findIndex(removalOrder_, idL)
+                            )
+                            {
+                                if
+                                (
+                                    Pstream::myProcNo() == refCell.sourceProc()
+                                 && cellK <= refCell.sourceCell()
+                                )
+                                {
+                                    if (findIndex(molsToDelete, molK) == -1)
+                                    {
+                                        molsToDelete.append(molK);
+                                    }
+                                }
+
+                                else if
+                                (
+                                    Pstream::myProcNo() < refCell.sourceProc()
+                                )
+                                {
+                                    if (findIndex(molsToDelete, molK) == -1)
+                                    {
+                                        molsToDelete.append(molK);
+                                    }
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    forAll (molsToDelete, mTD)
+    {
+        deleteParticle(*(molsToDelete[mTD]));
+    }
+}

src/lagrangian/molecule/moleculeCloud/moleculeCloudTestEdgeEdgeDistance.C

@@ -0,0 +1,86 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "moleculeCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+bool Foam::moleculeCloud::testEdgeEdgeDistance
+(
+    const edge& eI,
+    const edge& eJ
+) const
+{
+    const vector& eJs(mesh_.points()[eJ.start()]);
+    const vector& eJe(mesh_.points()[eJ.end()]);
+
+    return testEdgeEdgeDistance(eI, eJs, eJe);
+}
+
+bool Foam::moleculeCloud::testEdgeEdgeDistance
+(
+    const edge& eI,
+    const vector& eJs,
+    const vector& eJe
+) const
+{
+    vector a(eI.vec(mesh_.points()));
+    vector b(eJe - eJs);
+
+    const vector& eIs(mesh_.points()[eI.start()]);
+
+    vector c(eJs - eIs);
+
+    vector crossab = a ^ b;
+    scalar magCrossSqr = magSqr(crossab);
+
+    if (magCrossSqr > VSMALL)
+    {
+        // If the edges are parallel then a point-face
+        // search will pick them up
+
+        scalar s = ((c ^ b) & crossab)/magCrossSqr;
+        scalar t = ((c ^ a) & crossab)/magCrossSqr;
+
+        // Check for end points outside of range 0..1
+        // If the closest point is outside this range
+        // a point-face search will have found it.
+
+        return
+        (
+            s >= 0
+         && s <= 1
+         && t >= 0
+         && t <= 1
+         && magSqr(eIs + a*s - eJs - b*t) <= rCutMaxSqr()
+        );
+    }
+
+    return false;
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/moleculeCloud/moleculeCloudTestPointFaceDistance.C

@@ -0,0 +1,235 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "moleculeCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+bool Foam::moleculeCloud::testPointFaceDistance
+(
+    const label p,
+    const label faceNo
+) const
+{
+    const vector& pointPosition(mesh_.points()[p]);
+
+    return testPointFaceDistance(pointPosition, faceNo);
+}
+
+bool Foam::moleculeCloud::testPointFaceDistance
+(
+    const label p,
+    const referredCell& refCell
+) const
+{
+    const vector& pointPosition(mesh_.points()[p]);
+
+    forAll (refCell.faces(), rCF)
+    {
+        if
+        (
+            testPointFaceDistance
+            (
+                pointPosition,
+                refCell.faces()[rCF],
+                refCell.vertexPositions(),
+                refCell.faceCentres()[rCF],
+                refCell.faceAreas()[rCF]
+            )
+        )
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+bool Foam::moleculeCloud::testPointFaceDistance
+(
+    const vectorList& pointsToTest,
+    const label faceNo
+) const
+{
+    forAll(pointsToTest, pTT)
+    {
+        const vector& p(pointsToTest[pTT]);
+
+        // if any point in the list is in range of the face
+        // then the rest do not need to be tested and
+        // true can be returned
+
+        if (testPointFaceDistance(p, faceNo))
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+bool Foam::moleculeCloud::testPointFaceDistance
+(
+    const vector& p,
+    const label faceNo
+) const
+{
+    const face& faceToTest(mesh_.faces()[faceNo]);
+
+    const vector& faceC(mesh_.faceCentres()[faceNo]);
+
+    const vector& faceA(mesh_.faceAreas()[faceNo]);
+
+    const vectorList& points(mesh_.points());
+
+    return testPointFaceDistance
+    (
+        p,
+        faceToTest,
+        points,
+        faceC,
+        faceA
+    );
+}
+
+bool Foam::moleculeCloud::testPointFaceDistance
+(
+    const vector& p,
+    const labelList& faceToTest,
+    const vectorList& points,
+    const vector& faceC,
+    const vector& faceA
+) const
+{
+    vector faceN(faceA/mag(faceA));
+
+    scalar perpDist((p - faceC) & faceN);
+
+    if (mag(perpDist) > rCutMax())
+    {
+        return false;
+    }
+
+    vector pointOnPlane = (p - faceN * perpDist);
+
+    if (magSqr(faceC - pointOnPlane) < rCutMaxSqr()*1e-8)
+    {
+        // If pointOnPlane is very close to the face centre
+        // then defining the local axes will be inaccurate
+        // and it is very likely that pointOnPlane will be
+        // inside the face, so return true if the points
+        // are in range to be safe
+
+        return (magSqr(pointOnPlane - p) <= rCutMaxSqr());
+    }
+
+    vector xAxis = (faceC - pointOnPlane)/mag(faceC - pointOnPlane);
+
+    vector yAxis =
+        ((faceC - pointOnPlane) ^ faceN)
+       /mag((faceC - pointOnPlane) ^ faceN);
+
+    List<vector2D> local2DVertices(faceToTest.size());
+
+    forAll(faceToTest, fTT)
+    {
+        const vector& V(points[faceToTest[fTT]]);
+
+        if (magSqr(V-p) <= rCutMaxSqr())
+        {
+            return true;
+        }
+
+        local2DVertices[fTT] = vector2D
+        (
+            ((V - pointOnPlane) & xAxis),
+            ((V - pointOnPlane) & yAxis)
+        );
+    }
+
+    scalar localFaceCx((faceC - pointOnPlane) & xAxis);
+
+    scalar la_valid = -1;
+
+    forAll(local2DVertices, fV)
+    {
+        const vector2D& va(local2DVertices[fV]);
+
+        const vector2D& vb
+        (
+            local2DVertices[(fV + 1) % local2DVertices.size()]
+        );
+
+        if (mag(vb.y()-va.y()) > SMALL)
+        {
+            scalar la =
+                (
+                    va.x() - va.y()*((vb.x() - va.x())/(vb.y() - va.y()))
+                )
+               /localFaceCx;
+
+            scalar lv = -va.y()/(vb.y() - va.y());
+
+
+            if (la >= 0 && la <= 1 && lv >= 0 && lv <= 1)
+            {
+                la_valid = la;
+
+                break;
+            }
+        }
+    }
+
+    if (la_valid < 0)
+    {
+        // perpendicular point inside face, nearest point is pointOnPlane;
+        return (magSqr(pointOnPlane-p) <= rCutMaxSqr());
+    }
+    else
+    {
+        // perpendicular point outside face, nearest point is
+        // on edge that generated la_valid;
+        return
+        (
+            magSqr(pointOnPlane + la_valid*(faceC - pointOnPlane) - p)
+         <= rCutMaxSqr()
+        );
+    }
+
+    // if the algorithm hasn't returned anything by now then something has
+    // gone wrong.
+
+    FatalErrorIn("moleculeCloudTestPointAndFaceDistance.C") << nl
+        << "point " << p << " to face " << faceToTest
+        << " comparison did not find a nearest point"
+        << " to be inside or outside face."
+        << abort(FatalError);
+
+    return false;
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/potentials/pairPotential/basic/pairPotential.C

@@ -0,0 +1,308 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "pairPotential.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::pairPotential::setLookupTables()
+{
+    label N = label((rCut_ - rMin_)/dr_) + 1;
+
+    forceLookup_.setSize(N);
+
+    energyLookup_.setSize(N);
+
+    forAll(forceLookup_, k)
+    {
+        forceLookup_[k] = force(k*dr_ + rMin_);
+
+        energyLookup_[k] = energy(k*dr_ + rMin_);
+    }
+
+    forceLookup_[N-1] = 0.0;
+    energyLookup_[N-1] = 0.0;
+}
+
+
+void Foam::pairPotential::setConstants()
+{
+    scalar nr = n(rCut_);
+
+    u_at_rCut_ =
+        epsilon_
+       *(
+            (6.0/(nr - 6.0))*Foam::pow( (rCut_/rm_), -nr)
+          - (nr/(nr - 6.0))*Foam::pow( (rCut_/rm_), -6)
+        );
+
+    du_by_dr_at_rCut_ =
+        -6.0 * epsilon_ * gamma_
+       *(
+            Foam::pow( (rCut_/rm_),-nr)
+           *(
+                (rCut_/rm_)
+               *(1.0 / (nr - 6.0) + log(rCut_ / rm_) + 1.0)
+              + (m_/gamma_)
+              - 1.0
+            )
+          - Foam::pow((rCut_/rm_),-6.0)
+            *(rCut_ / ( (nr - 6.0) * rm_) + (nr/gamma_))
+        )
+       /(rCut_*(nr - 6.0));
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::pairPotential::pairPotential()
+{}
+
+
+Foam::pairPotential::pairPotential(const dictionary& pPDict)
+:
+    forceLookup_(),
+    energyLookup_(),
+    m_(),
+    gamma_(),
+    rm_(),
+    epsilon_(),
+    rCut_(),
+    rCutSqr_(),
+    u_at_rCut_(),
+    du_by_dr_at_rCut_(),
+    rMin_(),
+    dr_()
+{
+
+    m_ = readScalar(pPDict.lookup("m"));
+
+    gamma_ = readScalar(pPDict.lookup("gamma"));
+
+    rm_ = readScalar(pPDict.lookup("rm"));
+
+    epsilon_ = readScalar(pPDict.lookup("epsilon"));
+
+    rCut_ = readScalar(pPDict.lookup("rCut"));
+
+    rCutSqr_ = rCut_ * rCut_;
+
+    rMin_ = readScalar(pPDict.lookup("rMin"));
+
+    dr_ = readScalar(pPDict.lookup("dr"));
+
+    setConstants();
+
+    setLookupTables();
+}
+
+
+Foam::pairPotential::pairPotential
+(
+    const scalar m,
+    const scalar gamma,
+    const scalar rm,
+    const scalar epsilon,
+    const scalar rCut,
+    const scalar rMin,
+    const scalar dr
+)
+:
+    forceLookup_(),
+    energyLookup_(),
+    m_(m),
+    gamma_(gamma),
+    rm_(rm),
+    epsilon_(epsilon),
+    rCut_(rCut),
+    rCutSqr_(rCut*rCut),
+    u_at_rCut_(),
+    du_by_dr_at_rCut_(),
+    rMin_(rMin),
+    dr_(dr)
+{
+    setConstants();
+
+    setLookupTables();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::pairPotential::~pairPotential()
+{}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::pairPotential::write(Ostream& os) const
+{
+    os<< "Shifted force Maitland Smith pair potential."
+        << nl << tab << "m = " << m_
+        << nl << tab << "gamma = " << gamma_
+        << nl << tab << "rm = " << rm_
+        << nl << tab << "epsilon = " << epsilon_
+        << nl << tab << "rCut = " << rCut_
+        << nl << tab << "rMin = " << rMin_
+        << nl << tab << "dr = " << dr_
+        << endl;
+}
+
+
+Foam::scalar Foam::pairPotential::n(const scalar rIJMag) const
+{
+    return (m_ + gamma_*(rIJMag/rm_ - 1.0));
+}
+
+
+Foam::scalar Foam::pairPotential::force(const scalar rIJMag) const
+{
+    scalar nr(n(rIJMag));
+
+    scalar r = rIJMag/rm_;
+
+    scalar rN = Foam::pow(r,-nr);
+
+    scalar f =
+      - epsilon_
+       *(
+           - 6.0/(nr - 6.0)/(nr - 6.0)*(gamma_/rm_)*rN
+           + 6.0/(nr - 6.0)*rN
+             *((gamma_ - m_)/rIJMag - (gamma_/rm_)*(log(r) + 1.0))
+           + 6.0/(nr - 6.0)
+             *(1.0/(nr - 6.0)*(gamma_/rm_) + nr/rIJMag)*Foam::pow(r, -6)
+        )
+      + du_by_dr_at_rCut_;
+
+    return f;
+}
+
+
+Foam::scalar Foam::pairPotential::forceLookup(const scalar rIJMag) const
+{
+    scalar k_rIJ = (rIJMag - rMin_)/dr_;
+
+    label k(k_rIJ);
+
+    if (k < 0)
+    {
+        FatalErrorIn("pairPotential.C") << nl
+        << "rIJMag less than rMin" << nl
+        << abort(FatalError);
+    }
+
+    scalar f =
+        (k_rIJ - k)*forceLookup_[k+1]
+      + (k + 1 - k_rIJ)*forceLookup_[k];
+
+    return f;
+}
+
+
+Foam::List< Foam::Pair< Foam::scalar > >
+Foam::pairPotential::forceTable() const
+{
+    List< Pair<scalar> > forceTab(forceLookup_.size());
+
+    forAll(forceLookup_,k)
+    {
+        forceTab[k].first() = rMin_ + k*dr_;
+
+        forceTab[k].second() = forceLookup_[k];
+    }
+
+    return forceTab;
+}
+
+
+Foam::scalar Foam::pairPotential::energy(const scalar rIJMag) const
+{
+    scalar nr = n(rIJMag);
+
+    scalar e =
+        epsilon_
+       *(
+            (6.0 / (nr - 6.0))*Foam::pow( rIJMag/rm_, -nr)
+          - (nr / (nr - 6.0))*Foam::pow( rIJMag/rm_, -6)
+        )
+      - u_at_rCut_
+      - (rIJMag - rCut_)
+       *du_by_dr_at_rCut_;
+
+    return e;
+}
+
+
+Foam::scalar Foam::pairPotential::energyLookup(const scalar rIJMag) const
+{
+    scalar k_rIJ = (rIJMag - rMin_)/dr_;
+
+    label k(k_rIJ);
+
+    if (k < 0)
+    {
+        FatalErrorIn("pairPotential.C") << nl
+        << "rIJMag less than rMin" << nl
+        << abort(FatalError);
+    }
+
+    scalar e =
+        (k_rIJ - k)*energyLookup_[k+1]
+      + (k + 1 - k_rIJ)*energyLookup_[k];
+
+    return e;
+}
+
+
+Foam::List< Foam::Pair< Foam::scalar > >
+    Foam::pairPotential::energyTable() const
+{
+    List< Pair<scalar> > energyTab(energyLookup_.size());
+
+    forAll(energyLookup_,k)
+    {
+        energyTab[k].first() = rMin_ + k*dr_;
+
+        energyTab[k].second() = energyLookup_[k];
+    }
+
+    return energyTab;
+}
+
+
+// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const pairPotential& pot)
+{
+    pot.write(os);
+    os.check
+    (
+        "Foam::Ostream& Foam::pperator<<(Ostream& f, const pairPotential& pot"
+    );
+    return os;
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/potentials/pairPotential/basic/pairPotential.H

@@ -0,0 +1,196 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::pairPotential
+
+Description
+
+    At the moment this is hard coded to be a shifted force
+    "maitlandSmith" potential.
+    In future use templated classes virtual functions,
+    function pointers and all kinds of good stuff to make
+    this a generic pair force,
+
+SourceFiles
+    pairPotentialI.H
+    pairPotential.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pairPotential_H
+#define pairPotential_H
+
+#include "vector.H"
+#include "dictionary.H"
+#include "List.H"
+#include "Pair.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class pairPotential Declaration
+\*---------------------------------------------------------------------------*/
+
+class pairPotential
+{
+    // Private data
+
+        List<scalar> forceLookup_;
+
+        List<scalar> energyLookup_;
+
+        scalar m_;
+
+        scalar gamma_;
+
+        scalar rm_;
+
+        scalar epsilon_;
+
+        scalar rCut_;
+
+        scalar rCutSqr_;
+
+        scalar u_at_rCut_;
+
+        scalar du_by_dr_at_rCut_;
+
+        scalar rMin_;
+
+        scalar dr_;
+
+
+    // Private Member Functions
+
+        void setLookupTables();
+
+        void setConstants();
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        pairPotential();
+
+        //- Construct from dictionary
+        pairPotential(const dictionary& pPDict);
+
+        //- Construct from components
+        pairPotential
+        (
+            const scalar m,
+            const scalar gamma,
+            const scalar rm,
+            const scalar epsilon,
+            const scalar rCut,
+            const scalar rMin,
+            const scalar dr
+        );
+
+
+    // Destructor
+
+        virtual ~pairPotential();
+
+
+    // Member Functions
+
+        // Access
+
+            inline const scalar m() const;
+
+            inline const scalar gamma() const ;
+
+            inline const scalar rm() const;
+
+            inline const scalar epsilon() const;
+
+            inline const scalar rCut() const;
+
+            inline const scalar rCutSqr() const;
+
+            inline const scalar rMin() const;
+
+            inline const scalar dr() const;
+
+
+        // Write
+
+            virtual void write(Ostream&) const;
+
+            scalar n(const scalar rIJMag) const;
+
+            scalar force(const scalar rIJMag) const;
+
+            scalar forceLookup(const scalar rIJMag) const;
+
+            List<Pair<scalar> > forceTable() const;
+
+            scalar energy(const scalar rIJMag) const;
+
+            scalar energyLookup(const scalar rIJMag) const;
+
+            List<Pair<scalar> > energyTable() const;
+
+
+    // Friend Operators
+
+        inline friend bool operator==
+        (
+            const pairPotential& a,
+            const pairPotential& b
+        );
+
+        inline friend bool operator!=
+        (
+            const pairPotential& a,
+            const pairPotential& b
+        );
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<<(Ostream&, const pairPotential&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "pairPotentialI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/molecule/potentials/pairPotential/basic/pairPotentialI.H

@@ -0,0 +1,102 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+inline const Foam::scalar Foam::pairPotential::m() const
+{
+    return m_;
+}
+
+
+inline const Foam::scalar Foam::pairPotential::gamma() const
+{
+    return gamma_;
+}
+
+
+inline const Foam::scalar Foam::pairPotential::rm() const
+{
+    return rm_;
+}
+
+
+inline const Foam::scalar Foam::pairPotential::epsilon() const
+{
+    return epsilon_;
+}
+
+
+inline const Foam::scalar Foam::pairPotential::rCut() const
+{
+    return rCut_;
+}
+
+
+inline const Foam::scalar Foam::pairPotential::rCutSqr() const
+{
+    return rCutSqr_;
+}
+
+
+inline const Foam::scalar Foam::pairPotential::rMin() const
+{
+    return rMin_;
+}
+
+
+inline const Foam::scalar Foam::pairPotential::dr() const
+{
+    return dr_;
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+inline bool Foam::operator==
+(
+    const pairPotential& a,
+    const pairPotential& b
+)
+{
+    return
+        (a.m() == b.m())
+     && (a.gamma() == b.gamma())
+     && (a.rm() == b.rm())
+     && (a.epsilon() == b.epsilon())
+     && (a.rCut() == b.rCut());
+}
+
+
+inline bool Foam::operator!=
+(
+    const pairPotential& a,
+    const pairPotential& b
+)
+{
+    return !(a == b);
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecule/potentials/pairPotential/basic/pairPotentialList.C

@@ -0,0 +1,180 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "pairPotentialList.H"
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::pairPotentialList::pairPotentialList()
+:
+    List<pairPotential> ()
+{}
+
+
+Foam::pairPotentialList::pairPotentialList
+(
+    const label nIds
+)
+:
+    List<pairPotential> ((nIds * (nIds + 1))/2),
+    nIds_(nIds)
+{}
+
+
+Foam::pairPotentialList::pairPotentialList
+(
+    const List<pairPotential>& pairPotentials,
+    const label nIds
+)
+:
+    List<pairPotential> (pairPotentials),
+    nIds_(nIds)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+pairPotentialList::~pairPotentialList()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void pairPotentialList::setSizeAndNIds (const label nIds)
+{
+    nIds_ = nIds;
+    setSize((nIds * (nIds + 1))/2);
+}
+
+
+void pairPotentialList::addPotential
+(
+    const label a,
+    const label b,
+    const pairPotential& pot
+)
+{
+    if (pairPotentialIndex (a,b) > size()-1)
+    {
+        FatalErrorIn
+        (
+            "Foam::pairPotentialList::addPotential "
+            "(const label a, const label b, const pairPotential& pot)"
+        )<< "Attempting to add a pairPotential with too high an index"
+         << nl
+         << "Check if the pairPotentialList has been constructed "
+         << "with the number of ids to expect"
+         << nl << abort(FatalError);
+    }
+    else
+    {
+        (*this)[pairPotentialIndex (a,b)] = pot;
+    }
+}
+
+
+const pairPotential& pairPotentialList::pairPotentialFunction
+(
+    const label a,
+    const label b
+) const
+{
+    return (*this)[pairPotentialIndex (a,b)];
+}
+
+
+bool pairPotentialList::rCutSqr
+(
+    const label a,
+    const label b,
+    const scalar rIJMagSqr
+) const
+{
+    if(rIJMagSqr <= rCutSqr (a,b))
+    {
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+const scalar pairPotentialList::rCutSqr
+(
+    const label a,
+    const label b
+) const
+{
+    return (*this)[pairPotentialIndex (a,b)].rCutSqr();
+}
+
+
+const scalar pairPotentialList::rCut
+(
+    const label a,
+    const label b
+) const
+{
+    return (*this)[pairPotentialIndex (a,b)].rCut();
+}
+
+
+scalar pairPotentialList::force
+(
+    const label a,
+    const label b,
+    const scalar rIJMag
+) const
+{
+    scalar f = (*this)[pairPotentialIndex (a,b)].forceLookup(rIJMag);
+
+    return f;
+}
+
+
+scalar pairPotentialList::energy
+(
+    const label a,
+    const label b,
+    const scalar rIJMag
+) const
+{
+    scalar e = (*this)[pairPotentialIndex (a,b)].energyLookup(rIJMag);
+
+    return e;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/molecule/potentials/pairPotential/basic/pairPotentialList.H

@@ -0,0 +1,169 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::pairPotentialList
+
+Description
+
+SourceFiles
+    pairPotentialListI.H
+    pairPotentialList.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pairPotentialList_H
+#define pairPotentialList_H
+
+#include "pairPotential.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class pairPotentialList Declaration
+\*---------------------------------------------------------------------------*/
+
+class pairPotentialList
+:
+    public List<pairPotential>
+{
+    // Private data
+
+        label nIds_;
+
+
+    // Private Member Functions
+
+        inline label pairPotentialIndex
+        (
+            const label a,
+            const label b
+        ) const;
+
+        //- Disallow default bitwise assignment
+        void operator=(const pairPotentialList&);
+
+        //- Disallow default bitwise copy construct
+        pairPotentialList(const pairPotentialList&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        pairPotentialList();
+
+        //- Construct from number of Ids
+        pairPotentialList
+        (
+            const label nIds
+        );
+
+        //- Construct from components
+        pairPotentialList
+        (
+            const List<pairPotential>& pairPotentials,
+            const label nIds
+        );
+
+
+    // Destructor
+
+        ~pairPotentialList();
+
+
+    // Member Functions
+
+        // Access
+
+            inline label nIds() const;
+
+            void setSizeAndNIds (const label);
+
+            void addPotential
+            (
+                const label a,
+                const label b,
+                const pairPotential& pot
+            );
+
+            const pairPotential& pairPotentialFunction
+            (
+                const label a,
+                const label b
+            ) const;
+
+            // Return true if rIJ is within rCut for this pair.
+            bool rCutSqr
+            (
+                const label a,
+                const label b,
+                const scalar rIJMagSqr
+            ) const;
+
+            const scalar rCutSqr
+            (
+                const label a,
+                const label b
+            ) const;
+
+            const scalar rCut
+            (
+                const label a,
+                const label b
+            ) const;
+
+            scalar force
+            (
+                const label a,
+                const label b,
+                const scalar rIJMag
+            ) const;
+
+            scalar energy
+            (
+                const label a,
+                const label b,
+                const scalar rIJMag
+            ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "pairPotentialListI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/molecule/potentials/pairPotential/basic/pairPotentialListI.H

@@ -0,0 +1,59 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+inline Foam::label Foam::pairPotentialList::pairPotentialIndex
+(
+    const label a,
+    const label b
+) const
+{
+    label index;
+
+    if (a < b)
+    {
+        index = a*(2*nIds_ - a - 1)/2 + b;
+    }
+
+    else
+    {
+        index = b*(2*nIds_ - b - 1)/2 + a;
+    }
+
+    return index;
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline Foam::label Foam::pairPotentialList::nIds() const
+{
+    return nIds_;
+}
+
+
+// ************************************************************************* //