diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C b/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C
index f2ddf13741..c3ae69847d 100644
--- a/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C
+++ b/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C
@@ -39,6 +39,20 @@ int main(int argc, char *argv[])
 #   include "createTime.H"
 #   include "createMesh.H"
 
+    Info<< "\nReading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
     potential pot(mesh);
 
     moleculeCloud molecules(mesh, pot);
diff --git a/applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C b/applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C
index 3aa2532c3e..ade0895a2a 100644
--- a/applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C
+++ b/applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C
@@ -39,6 +39,20 @@ int main(int argc, char *argv[])
 #   include "createTime.H"
 #   include "createMesh.H"
 
+    Info<< "\nReading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
     potential pot(mesh);
 
     moleculeCloud molecules(mesh, pot);
diff --git a/src/lagrangian/molecularDynamics/molecule/Make/files b/src/lagrangian/molecularDynamics/molecule/Make/files
index a8a36e5967..9086ba7ffa 100644
--- a/src/lagrangian/molecularDynamics/molecule/Make/files
+++ b/src/lagrangian/molecularDynamics/molecule/Make/files
@@ -1,24 +1,9 @@
-interactionLists = interactionLists
-referredMolecule = $(interactionLists)/referredMolecule
-referredCellList = $(interactionLists)/referredCellList
-referredCell = $(interactionLists)/referredCell
-directInteractionList = $(interactionLists)/directInteractionList
+reducedUnits/reducedUnits.C
+reducedUnits/reducedUnitsIO.C
 
-$(referredCellList)/referredCellList.C
-$(referredCell)/referredCell.C
-$(referredMolecule)/referredMolecule.C
-$(directInteractionList)/directInteractionList.C
-$(interactionLists)/interactionLists.C
+molecule/molecule.C
+molecule/moleculeIO.C
 
-reducedUnits = reducedUnits
-$(reducedUnits)/reducedUnits.C
-$(reducedUnits)/reducedUnitsIO.C
-
-molecule = molecule
-$(molecule)/molecule.C
-$(molecule)/moleculeIO.C
-
-moleculeCloud = moleculeCloud
-$(moleculeCloud)/moleculeCloud.C
+moleculeCloud/moleculeCloud.C
 
 LIB = $(FOAM_LIBBIN)/libmolecule
diff --git a/src/lagrangian/molecularDynamics/molecule/Make/options b/src/lagrangian/molecularDynamics/molecule/Make/options
index 9354036e6f..1bdb57a0d0 100644
--- a/src/lagrangian/molecularDynamics/molecule/Make/options
+++ b/src/lagrangian/molecularDynamics/molecule/Make/options
@@ -1,11 +1,13 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
     -I$(LIB_SRC)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lmeshTools \
     -llagrangian \
     -lpotential \
     -lmolecularMeasurements
diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionList.C b/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionList.C
deleted file mode 100644
index da9511daf5..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionList.C
+++ /dev/null
@@ -1,346 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "directInteractionList.H"
-#include "interactionLists.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::directInteractionList::buildDirectInteractionList
-(
-    bool pointPointListBuild
-)
-{
-    Info<< nl << "Building list of direct interaction neighbours" << endl;
-
-    const polyMesh& mesh(il_.mesh());
-
-    List<DynamicList<label> > directInteractionList(mesh.nCells());
-
-    if (pointPointListBuild)
-    {
-        Info<< tab << "Point-Point direct interaction list build." << endl;
-
-        label pointJIndex;
-
-        forAll(mesh.points(), pointIIndex)
-        {
-            for
-            (
-                pointJIndex = pointIIndex;
-                pointJIndex != mesh.points().size();
-                ++pointJIndex
-            )
-            {
-                if (il_.testPointPointDistance(pointIIndex, pointJIndex))
-                {
-                    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);
-                                }
-                            }
-                        }
-                    }
-                }
-            }
-        }
-    }
-    else
-    {
-        Info<< tab << "Point-Face, Edge-Edge direct interaction list build."
-            << endl;
-
-        forAll(mesh.points(), p)
-        {
-            forAll(mesh.faces(), f)
-            {
-                if (il_.testPointFaceDistance(p, f))
-                {
-                    const labelList& pCells(mesh.pointCells()[p]);
-
-                    const label cellO(mesh.faceOwner()[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
-
-                            const label cellN(mesh.faceNeighbour()[f]);
-
-                            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 (il_.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);
-                                }
-                            }
-                        }
-                    }
-                }
-            }
-        }
-    }
-
-    forAll(directInteractionList, transDIL)
-    {
-        (*this)[transDIL].transfer
-        (
-            directInteractionList[transDIL].shrink()
-        );
-    }
-
-    // sorting DILs
-
-    forAll((*this), dIL)
-    {
-        sort((*this)[dIL]);
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::directInteractionList::directInteractionList
-(
-    const interactionLists& il,
-    bool pointPointListBuild
-)
-:
-    labelListList(il.mesh().nCells()),
-    il_(il)
-{
-    if ((*this).size() > 1)
-    {
-        buildDirectInteractionList(pointPointListBuild);
-    }
-    else if ((*this).size() == 1)
-    {
-        Info<< nl
-            << "Single cell mesh, no direct interaction lists required."
-            << endl;
-
-        (*this)[0].setSize(0);
-    }
-}
-
-
-Foam::directInteractionList::directInteractionList
-(
-    const interactionLists& il
-)
-:
-    labelListList(il.mesh().nCells()),
-    il_(il)
-{
-    Info<< "Read directInteractionList from disk not implemented" << endl;
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::directInteractionList::~directInteractionList()
-{}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionList.H b/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionList.H
deleted file mode 100644
index 1c7c6272ee..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionList.H
+++ /dev/null
@@ -1,124 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::directInteractionList
-
-Description
-
-SourceFiles
-    directInteractionListI.H
-    directInteractionList.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef directInteractionList_H
-#define directInteractionList_H
-
-#include "polyMesh.H"
-#include "List.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-class interactionLists;
-
-/*---------------------------------------------------------------------------*\
-                           Class directInteractionList Declaration
-\*---------------------------------------------------------------------------*/
-
-class directInteractionList
-:
-    public labelListList
-{
-    // Private data
-
-        const interactionLists& il_;
-
-
-    // Private Member Functions
-
-        void buildDirectInteractionList
-        (
-            bool pointPointListBuild
-        );
-
-        //- Disallow default bitwise copy construct
-        directInteractionList(const directInteractionList&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const directInteractionList&);
-
-
-public:
-
-    // Constructors
-
-        //- Construct lists by searching the mesh
-        directInteractionList
-        (
-            const interactionLists& il,
-            bool pointPointListBuild
-        );
-
-        //- Construct from file
-        directInteractionList
-        (
-            const interactionLists& il
-        );
-
-
-    //- Destructor
-    ~directInteractionList();
-
-
-    // Member Functions
-
-        // Access
-
-        inline const interactionLists& il() const;
-
-
-    // IOstream Operators
-
-        friend Istream& operator>>(Istream&, directInteractionList&);
-        friend Ostream& operator<<(Ostream&, const directInteractionList&);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "directInteractionListI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionListI.H b/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionListI.H
deleted file mode 100644
index 41b30a8ab1..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionListI.H
+++ /dev/null
@@ -1,34 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-inline const Foam::interactionLists& Foam::directInteractionList::il() const
-{
-    return il_;
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.C b/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.C
deleted file mode 100644
index c9e7a6aa4e..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.C
+++ /dev/null
@@ -1,691 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "interactionLists.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-Foam::scalar Foam::interactionLists::transTol = 1e-12;
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::interactionLists::buildCellReferralLists()
-{
-    Info<< nl << "Determining molecule referring schedule" << endl;
-
-    const referredCellList& refIntL(ril());
-
-    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]
-            );
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-
-Foam::interactionLists::interactionLists
-(
-    const polyMesh& mesh,
-    scalar rCutMaxSqr,
-    bool pointPointListBuild
-)
-:
-    mesh_(mesh),
-    rCutMaxSqr_(rCutMaxSqr),
-    dil_(*this, pointPointListBuild),
-    ril_(*this, pointPointListBuild),
-    cellSendingReferralLists_(),
-    cellReceivingReferralLists_()
-{
-    buildCellReferralLists();
-}
-
-
-Foam::interactionLists::interactionLists(const polyMesh& mesh)
-:
-    mesh_(mesh),
-    dil_(*this),
-    ril_(*this)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::interactionLists::~interactionLists()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::interactionLists::testPointPointDistance
-(
-    const label ptI,
-    const label ptJ
-) const
-{
-    return (magSqr(mesh_.points()[ptI] - mesh_.points()[ptJ]) <= rCutMaxSqr_);
-}
-
-
-bool Foam::interactionLists::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::interactionLists::testPointFaceDistance
-(
-    const label p,
-    const label faceNo
-) const
-{
-    const vector& pointPosition(mesh_.points()[p]);
-
-    return testPointFaceDistance(pointPosition, faceNo);
-}
-
-
-bool Foam::interactionLists::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::interactionLists::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::interactionLists::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::interactionLists::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 (magSqr(perpDist) > rCutMaxSqr_)
-    {
-        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
-    (
-        "interactionLists::testPointFaceDistance"
-        "("
-            "const vector&, "
-            "const labelList&, "
-            "const vectorList&, "
-            "const vector&, "
-            "const vector&"
-        ") const"
-    )   << "point " << p << " to face " << faceToTest
-        << " comparison did not find a nearest point"
-        << " to be inside or outside face."
-        << abort(FatalError);
-
-    return false;
-}
-
-
-bool Foam::interactionLists::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;
-}
-
-
-const Foam::labelList Foam::interactionLists::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);
-                    }
-                }
-            }
-        }
-    }
-
-    return realCellsFoundInRange.shrink();
-}
-
-
-const Foam::labelList Foam::interactionLists::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();
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.H b/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.H
deleted file mode 100644
index 67e053b79e..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.H
+++ /dev/null
@@ -1,214 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::interactionLists
-
-Description
-
-SourceFiles
-    interactionListsI.H
-    interactionLists.C
-    interactionListsIO.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef interactionLists_H
-#define interactionLists_H
-
-#include "polyMesh.H"
-#include "vector2D.H"
-#include "directInteractionList.H"
-#include "referredCell.H"
-#include "referredCellList.H"
-#include "sendingReferralList.H"
-#include "receivingReferralList.H"
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                       Class interactionLists Declaration
-\*---------------------------------------------------------------------------*/
-
-class interactionLists
-{
-    // Private data
-
-        const polyMesh& mesh_;
-
-        scalar rCutMaxSqr_;
-
-        directInteractionList dil_;
-
-        referredCellList ril_;
-
-        List<sendingReferralList> cellSendingReferralLists_;
-
-        List<receivingReferralList> cellReceivingReferralLists_;
-
-
-    // Private Member Functions
-
-        //- Build referralLists which define how to send information
-        // to referredCells to source cells
-        void buildCellReferralLists();
-
-        //- Disallow default bitwise copy construct
-        interactionLists(const interactionLists&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const interactionLists&);
-
-
-public:
-
-    // Static data members
-
-        //- Tolerance for checking that faces on a patch segment
-        static scalar transTol;
-
-
-    // Constructors
-
-        //- Construct and create all information from the mesh
-        interactionLists
-        (
-            const polyMesh& mesh,
-            scalar rCutMaxSqr,
-            bool pointPointListBuild = false
-        );
-
-        //- Construct from file
-        interactionLists(const polyMesh& mesh);
-
-
-    //- Destructor
-    ~interactionLists();
-
-
-    // Member Functions
-
-            bool testPointPointDistance
-            (
-                const label ptI,
-                const label ptJ
-            ) const;
-
-            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;
-
-
-        // Access
-
-            inline const polyMesh& mesh() const;
-
-            inline const directInteractionList& dil() const;
-
-            inline const referredCellList& ril() const;
-            inline referredCellList& ril();
-
-            inline const List<sendingReferralList>&
-                cellSendingReferralLists() const;
-
-            inline const List<receivingReferralList>&
-                cellReceivingReferralLists() const;
-
-            inline label nInteractingProcs() const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "interactionListsI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionListsI.H b/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionListsI.H
deleted file mode 100644
index 382d4c44ec..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionListsI.H
+++ /dev/null
@@ -1,77 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-const Foam::polyMesh& Foam::interactionLists::mesh() const
-{
-    return mesh_;
-}
-
-
-const Foam::directInteractionList& Foam::interactionLists::dil() const
-{
-    return dil_;
-}
-
-
-inline const Foam::referredCellList& Foam::interactionLists::ril() const
-{
-    return ril_;
-}
-
-
-inline Foam::referredCellList& Foam::interactionLists::ril()
-{
-    return ril_;
-}
-
-
-inline const Foam::List<Foam::sendingReferralList>&
-    Foam::interactionLists::cellSendingReferralLists() const
-{
-    return cellSendingReferralLists_;
-}
-
-
-inline const Foam::List<Foam::receivingReferralList>&
-    Foam::interactionLists::cellReceivingReferralLists() const
-{
-    return cellReceivingReferralLists_;
-}
-
-
-inline Foam::label Foam::interactionLists::nInteractingProcs() const
-{
-    return cellReceivingReferralLists_.size();
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCell/referredCell.C b/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCell/referredCell.C
deleted file mode 100644
index bf0a4f7021..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCell/referredCell.C
+++ /dev/null
@@ -1,448 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*----------------------------------------------------------------------------*/
-
-#include "referredCell.H"
-#include "interactionLists.H"
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void referredCell::setConstructionData
-(
-    const polyMesh& mesh,
-    const label sourceCell
-)
-{
-    // Points
-
-    const labelList& points = mesh.cellPoints()[sourceCell];
-
-    vectorList sourceCellVertices(points.size());
-
-    forAll(sourceCellVertices, sCV)
-    {
-        sourceCellVertices[sCV] = mesh.points()[points[sCV]];
-    }
-
-    vertexPositions_ = referPositions(sourceCellVertices);
-
-
-    // Edges
-
-    const labelList& edges = mesh.cellEdges()[sourceCell];
-
-    edgeList sourceCellEdges(edges.size());
-
-    forAll(sourceCellEdges, sCE)
-    {
-        sourceCellEdges[sCE] = mesh.edges()[edges[sCE]];
-    }
-
-    locallyMapEdgeList(points, sourceCellEdges);
-
-
-    // Faces
-
-    labelList faces(mesh.cells()[sourceCell]);
-
-    vectorList sourceCellFaceCentres(faces.size());
-
-    vectorList sourceCellFaceAreas(faces.size());
-
-    labelListList sourceCellFaces(faces.size());
-
-    forAll(faces, f)
-    {
-        sourceCellFaces[f] = mesh.faces()[faces[f]];
-
-        sourceCellFaceCentres[f] = mesh.faceCentres()[faces[f]];
-
-        sourceCellFaceAreas[f] = mesh.faceAreas()[faces[f]];
-    }
-
-    locallyMapFaceList(points, sourceCellFaces);
-
-    faceCentres_ = referPositions(sourceCellFaceCentres);
-
-    faceAreas_ = rotateVectors(sourceCellFaceAreas);
-}
-
-
-void referredCell::locallyMapEdgeList
-(
-    const labelList& points,
-    const edgeList& sourceCellEdges
-)
-{
-    edges_.setSize(sourceCellEdges.size());
-
-    forAll(sourceCellEdges, sCE)
-    {
-        const edge& e(sourceCellEdges[sCE]);
-
-        edges_[sCE].start() = findIndex(points, e.start());
-
-        edges_[sCE].end() = findIndex(points, e.end());
-
-        if
-        (
-            edges_[sCE].start() == -1
-         || edges_[sCE].end() == -1
-        )
-        {
-            FatalErrorIn("Foam::referredCell::locallyMapEdgeList")
-                << "edgeList and points labelList for "
-                << "referred cell do not match: "
-                << nl << "points: " << points
-                << nl << "egdes: " << sourceCellEdges
-                << abort(FatalError);
-        }
-    }
-}
-
-
-void referredCell::locallyMapFaceList
-(
-    const labelList& points,
-    const labelListList& sourceCellFaces
-)
-{
-    faces_.setSize(sourceCellFaces.size());
-
-    forAll(sourceCellFaces, sCF)
-    {
-        const labelList& sourceCellFace(sourceCellFaces[sCF]);
-
-        labelList& localFace(faces_[sCF]);
-
-        localFace.setSize(sourceCellFace.size());
-
-        forAll(sourceCellFace, p)
-        {
-            localFace[p] = findIndex(points, sourceCellFace[p]);
-
-            if (localFace[p] == -1)
-            {
-                FatalErrorIn("Foam::referredCell::locallyMapEdgeList")
-                    << "edgeList and points labelList for "
-                    << "referred cell do not match: "
-                    << nl << "points: " << points
-                    << nl << "faces: " << sourceCellFaces
-                    << abort(FatalError);
-            }
-        }
-    }
-}
-
-
-vector referredCell::referPosition(const vector& positionToRefer)
-{
-    return offset_ + (rotation_ & positionToRefer);
-}
-
-
-vectorList referredCell::referPositions(const vectorList& positionsToRefer)
-{
-    return offset_ + (rotation_ & positionsToRefer);
-}
-
-
-vector referredCell::rotateVector(const vector& vectorToRotate)
-{
-    return rotation_ & vectorToRotate;
-}
-
-
-vectorList referredCell::rotateVectors(const vectorList& vectorsToRotate)
-{
-    return rotation_ & vectorsToRotate;
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-referredCell::referredCell()
-:
-    DynamicList<referredMolecule>(),
-    sourceProc_(-1),
-    sourceCell_(-1),
-    vertexPositions_(),
-    offset_(vector::zero),
-    rotation_(I)
-{}
-
-
-referredCell::referredCell
-(
-    const polyMesh& mesh,
-    const label sourceProc,
-    const label sourceCell,
-    const vector& offset,
-    const tensor& rotation
-)
-:
-    DynamicList<referredMolecule>(),
-    sourceProc_(sourceProc),
-    sourceCell_(sourceCell),
-    offset_(offset),
-    rotation_(rotation)
-{
-    setConstructionData(mesh, sourceCell);
-}
-
-
-referredCell::referredCell
-(
-    const label sourceProc,
-    const label sourceCell,
-    const vectorList& vertexPositions,
-    const edgeList& localEdges,
-    const labelListList& localFaces,
-    const vectorList& faceCentres,
-    const vectorList& faceAreas,
-    const vector& offset,
-    const tensor& rotation
-)
-:
-    DynamicList<referredMolecule>(),
-    sourceProc_(sourceProc),
-    sourceCell_(sourceCell),
-    edges_(localEdges),
-    faces_(localFaces),
-    offset_(offset),
-    rotation_(rotation)
-{
-    // Supplied vertexPositions, faceCentres and faceAreas are of the
-    // "original" cell, and need to be transformed to the referred
-    // locations on construction
-
-    vertexPositions_ = referPositions(vertexPositions);
-
-    faceCentres_ = referPositions(faceCentres);
-
-    faceAreas_ = rotateVectors(faceAreas);
-}
-
-
-referredCell::referredCell
-(
-    const polyMesh& mesh,
-    const label sourceProc,
-    const label sourceCell,
-    const vector& cS,
-    const vector& cD,
-    const vector& nS,
-    const vector& nD
-)
-:
-    DynamicList<referredMolecule>(),
-    sourceProc_(sourceProc),
-    sourceCell_(sourceCell)
-{
-    // It is assumed that the vectors originating from the faces being referred
-    // here are correct periodic faces - i.e. they have the same area etc.
-
-    vector nA = -nS/mag(nS);
-    vector nB = nD/mag(nD);
-
-    rotation_ = rotationTensor(nA, nB);
-
-    offset_ = cD - (rotation_ & cS);
-
-    // Allow sourceCell = -1 to create a dummy referredCell
-    // to obtain the transformation
-
-    if (sourceCell >= 0)
-    {
-        setConstructionData(mesh, sourceCell);
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-referredCell::~referredCell()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-referredCell referredCell::reRefer
-(
-    const vector& cS,
-    const vector& cD,
-    const vector& nS,
-    const vector& nD
-)
-{
-    vector nA = -nS/mag(nS);
-    vector nB = nD/mag(nD);
-
-    tensor newRotation = rotationTensor(nA, nB);
-
-    vector newOffset = cD - (newRotation & cS);
-
-    tensor reReferredRotation = newRotation & rotation_;
-
-    vector reReferredOffset = newOffset + (newRotation & offset_);
-
-    return referredCell
-    (
-        sourceProc_,
-        sourceCell_,
-        rotation_.T() & (vertexPositions_ - offset_),
-        edges_,
-        faces_,
-        rotation_.T() & (faceCentres_ - offset_),
-        rotation_.T() & (faceAreas_),
-        reReferredOffset,
-        reReferredRotation
-    );
-}
-
-
-vector referredCell::referPosition(const vector& positionToRefer) const
-{
-    return offset_ + (rotation_ & positionToRefer);
-}
-
-
-vectorList referredCell::referPosition
-(
-    const vectorList& positionsToRefer
-) const
-{
-    return offset_ + (rotation_ & positionsToRefer);
-}
-
-
-vector referredCell::rotateVector(const vector& vectorToRotate) const
-{
-    return rotation_ & vectorToRotate;
-}
-
-
-vectorList referredCell::rotateVectors(const vectorList& vectorsToRotate) const
-{
-    return rotation_ & vectorsToRotate;
-}
-
-
-void referredCell::referInMols(const List<referredMolecule>& incomingMols)
-{
-    clear();
-
-    forAll(incomingMols, iM)
-    {
-        append
-        (
-            referredMolecule
-            (
-                incomingMols[iM].id(),
-                referPosition
-                (
-                    incomingMols[iM].position()
-                ),
-                referPosition
-                (
-                    incomingMols[iM].sitePositions()
-                )
-            )
-        );
-    }
-
-    shrink();
-}
-
-
-bool referredCell::duplicate(const referredCell& refCellDupl) const
-{
-    return
-    (
-        sourceProc_ == refCellDupl.sourceProc()
-     && sourceCell_ == refCellDupl.sourceCell()
-     && mag(offset_ - refCellDupl.offset()) < interactionLists::transTol
-    );
-}
-
-
-bool referredCell::duplicate(const label procNo,const label nCells) const
-{
-    return
-    (
-        sourceProc_ == procNo
-     && sourceCell_ < nCells
-     && mag(offset_) < interactionLists::transTol
-    );
-}
-
-
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
-
-Istream& operator>>(Istream& is, referredCell& rC)
-{
-
-    is  >> rC.sourceProc_
-        >> rC.sourceCell_
-        >> rC.vertexPositions_
-        >> rC.edges_
-        >> rC.faces_
-        >> rC.faceCentres_
-        >> rC.faceAreas_
-        >> rC.offset_
-        >> rC.rotation_;
-
-    is.check("Istream& operator<<(Istream& f, const referredCell& rC");
-
-    return is;
-}
-
-
-Ostream& operator<<(Ostream& os, const referredCell& rC)
-{
-
-    os  << rC.sourceProc()
-        << token::SPACE << rC.sourceCell()
-        << token::SPACE << rC.vertexPositions()
-        << token::SPACE << rC.edges()
-        << token::SPACE << rC.faces()
-        << token::SPACE << rC.faceCentres()
-        << token::SPACE << rC.faceAreas()
-        << token::SPACE << rC.offset()
-        << token::SPACE << rC.rotation();
-
-    os.check("Ostream& operator<<(Ostream& f, const referredCell& rC");
-
-    return os;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-}  // End namespace Foam
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCell/referredCell.H b/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCell/referredCell.H
deleted file mode 100644
index a954c13f8a..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCell/referredCell.H
+++ /dev/null
@@ -1,275 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::referredCell
-
-Description
-
-SourceFiles
-    referredCellI.H
-    referredCell.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef referredCell_H
-#define referredCell_H
-
-#include "vector.H"
-#include "vectorList.H"
-#include "tensor.H"
-#include "transform.H"
-#include "DynamicList.H"
-#include "labelList.H"
-#include "edgeList.H"
-#include "polyMesh.H"
-
-#include "referredMolecule.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                          Class referredCell Declaration
-\*---------------------------------------------------------------------------*/
-
-class referredCell
-:
-    public DynamicList<referredMolecule>
-{
-    // Private data
-
-        label sourceProc_;
-
-        label sourceCell_;
-
-        //- Referred vertex positions
-        vectorList vertexPositions_;
-
-        edgeList edges_;
-
-        labelListList faces_;
-
-        vectorList faceCentres_;
-
-        vectorList faceAreas_;
-
-        labelList realCellsForInteraction_;
-
-        vector offset_;
-
-        tensor rotation_;
-
-
-    // Private Member Functions
-
-        void setConstructionData
-        (
-            const polyMesh& mesh,
-            const label sourceCell
-        );
-
-        void locallyMapEdgeList
-        (
-            const labelList& points,
-            const edgeList& sourceCellEdges
-        );
-
-        void locallyMapFaceList
-        (
-            const labelList& points,
-            const labelListList& sourceCellFaces
-        );
-
-        vector referPosition(const vector& positionToRefer);
-
-        vectorList referPositions(const vectorList& positionsToRefer);
-
-        vector rotateVector(const vector& vectorToRotate);
-
-        vectorList rotateVectors(const vectorList& vectorsToRotate);
-
-
-public:
-
-    // Constructors
-
-        //- Construct null
-
-        referredCell();
-
-        //- Construct from components with external edge information
-        referredCell
-        (
-            const polyMesh& mesh,
-            const label sourceProc,
-            const label sourceCell,
-            const vector& offset,
-            const tensor& rotation
-        );
-
-        //- Construct from components with existing local edge information
-        referredCell
-        (
-            const label sourceProc,
-            const label sourceCell,
-            const vectorList& vertexPositions,
-            const edgeList& localEdges,
-            const labelListList& localFaces,
-            const vectorList& faceCentres,
-            const vectorList& faceAreas,
-            const vector& offset,
-            const tensor& rotation
-        );
-
-        //- Construct from pair of face centers (c) and plain
-        //  face normals (n) (no need to make unit vectors or
-        //  reverse one direction)
-        //  Order of vectors important (S = source, D = Destination).
-        //  External edge information.
-
-        referredCell
-        (
-            const polyMesh& mesh,
-            const label sourceProc,
-            const label sourceCell,
-            const vector& cS,
-            const vector& cD,
-            const vector& nS,
-            const vector& nD
-        );
-
-
-    //- Destructor
-    virtual ~referredCell();
-
-
-    // Member Functions
-
-        //- Take this referredCell object that has already had its transform
-        //  calculated and refer it on again, retaining same source info.
-        referredCell reRefer
-        (
-            const vector& cS,
-            const vector& cD,
-            const vector& nS,
-            const vector& nD
-        );
-
-        //- Use internal transformatation values to transform the given
-        //  postion to its new location.
-        vector referPosition(const vector& positionToRefer) const;
-
-        //- Use internal transformatation values to transform the given
-        //  list of postions to their new locations.
-        vectorList referPosition(const vectorList& positionsToRefer) const;
-
-        //- Use internal transformatation values to rotate the given vector
-        vector rotateVector(const vector& vectorToRotate) const;
-
-        //- Use internal transformatation values to rotate the given
-        // list of vectors
-        vectorList rotateVectors(const vectorList& vectorsToRotate) const;
-
-        //- referInMols clears, the stored list of referred mols and takes
-        //  in a list of referred mols coming from a source processor,
-        //  transformingtheir positions
-        void referInMols(const List<referredMolecule>& incomingMols);
-
-        //- duplicate() function to test whether a referred or real cell
-        //  supplied by arguement is a duplicate of this referredCell.
-        //  Can be used bi-directionally - i.e. can be called on an existing
-        //  referred cell with a proposed referredCell as argument,
-        //  or vice versa.  Can only be called by a proposed referredCell with
-        //  a real cell index as arguement to test to see if the proposed
-        //  referredCell is a duplicate.
-        //  A duplicate cell is defined as one which has the same source
-        //  processor, source cell, and an equal offset. Real cells have zero
-        //  offset by definition.
-        bool duplicate(const referredCell& refCellDupl) const;
-
-        bool duplicate(const label procNo, const label nCells) const;
-
-
-        // Access
-
-            inline label sourceProc() const;
-
-            inline label sourceCell() const;
-
-            inline const vector& offset() const;
-
-            inline const tensor& rotation() const;
-
-            inline const vectorList& vertexPositions() const;
-
-            inline const edgeList& edges() const;
-
-            inline const labelListList& faces() const;
-
-            inline const vectorList& faceCentres() const;
-
-            inline const vectorList& faceAreas() const;
-
-            inline labelList& realCells();
-
-            inline const labelList& realCellsForInteraction() const;
-
-
-    // Friend Operators
-
-        inline friend bool operator==
-        (
-            const referredCell& a,
-            const referredCell& b
-        );
-
-        inline friend bool operator!=
-        (
-            const referredCell& a,
-            const referredCell& b
-        );
-
-
-    // IOstream Operators
-
-        friend Istream& operator>>(Istream&, referredCell&);
-        friend Ostream& operator<<(Ostream&, const referredCell&);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "referredCellI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCell/referredCellI.H b/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCell/referredCellI.H
deleted file mode 100644
index 0170068121..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCell/referredCellI.H
+++ /dev/null
@@ -1,126 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-inline label referredCell::sourceProc() const
-{
-    return sourceProc_;
-}
-
-
-inline label referredCell::sourceCell() const
-{
-    return sourceCell_;
-}
-
-
-inline const vector& referredCell::offset() const
-{
-    return offset_;
-}
-
-
-inline const tensor& referredCell::rotation() const
-{
-    return rotation_;
-}
-
-
-inline const vectorList& referredCell::vertexPositions() const
-{
-    return vertexPositions_;
-}
-
-
-inline const edgeList& referredCell::edges() const
-{
-    return edges_;
-}
-
-
-inline const labelListList& referredCell::faces() const
-{
-    return faces_;
-}
-
-
-inline const vectorList& referredCell::faceCentres() const
-{
-    return faceCentres_;
-}
-
-
-inline const vectorList& referredCell::faceAreas() const
-{
-    return faceAreas_;
-}
-
-
-inline labelList& referredCell::realCells()
-{
-    return realCellsForInteraction_;
-}
-
-
-inline const labelList& referredCell::realCellsForInteraction() const
-{
-    return realCellsForInteraction_;
-}
-
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
-inline bool operator==
-(
-    const referredCell& a,
-    const referredCell& b
-)
-{
-    return const_cast<referredCell&>(a).duplicate
-    (
-        const_cast<const referredCell&>(b)
-    );
-}
-
-
-inline bool operator!=
-(
-    const referredCell& a,
-    const referredCell& b
-)
-{
-    return !(a == b);
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCellList/referredCellList.C b/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCellList/referredCellList.C
deleted file mode 100644
index 32a6325898..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCellList/referredCellList.C
+++ /dev/null
@@ -1,1583 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*----------------------------------------------------------------------------*/
-
-#include "referredCellList.H"
-#include "interactionLists.H"
-#include "polyBoundaryMeshEntries.H"
-#include "PstreamCombineReduceOps.H"
-#include "Time.H"
-#include "globalMeshData.H"
-#include "processorPolyPatch.H"
-#include "cyclicPolyPatch.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::referredCellList::buildReferredCellList
-(
-    bool pointPointListBuild
-)
-{
-    Info<< nl << "Building list of referred interaction neighbours" << endl;
-
-    const polyMesh& mesh(il_.mesh());
-
-    DynamicList<referredCell> referredInteractionList;
-
-    // realCellsWithinRCutMaxOfAnyReferringPatch
-    DynamicList<label> rCellsWRRP;
-
-    // realFacesWithinRCutMaxOfAnyReferringPatch
-    DynamicList<label> rFacesWRRP;
-
-    // realEdgesWithinRCutMaxOfAnyReferringPatch
-    DynamicList<label> rEdgesWRRP;
-
-    // realPointsWithinRCutMaxOfAnyReferringPatch
-    DynamicList<label> rPointsWRRP;
-
-    labelListList processorPatchSegmentMapping
-    (
-        mesh.globalData().processorPatches().size()
-    );
-
-    List<vectorList> allNeighbourFaceCentres
-    (
-        mesh.globalData().processorPatches().size()
-    );
-
-    List<vectorList> allNeighbourFaceAreas
-    (
-        mesh.globalData().processorPatches().size()
-    );
-
-    label nUndecomposedPatches = 0;
-
-    if (Pstream::parRun())
-    {
-        dictionary patchDictionary;
-
-        DynamicList<word> patchNames;
-
-        Time undecomposedTime
-        (
-            Time::controlDictName,
-            mesh.time().rootPath(),
-            mesh.time().caseName().path()
-        );
-
-        IOobject undecomposedBoundaryHeader
-        (
-            "boundary",
-            undecomposedTime.constant(),
-            polyMesh::meshSubDir,
-            undecomposedTime,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false
-        );
-
-        if (undecomposedBoundaryHeader.headerOk())
-        {
-            polyBoundaryMeshEntries undecomposedPatchEntries
-            (
-                undecomposedBoundaryHeader
-            );
-
-            forAll(undecomposedPatchEntries, patchi)
-            {
-                patchNames.append
-                (
-                    undecomposedPatchEntries[patchi].keyword()
-                );
-
-                patchDictionary.add
-                (
-                    undecomposedPatchEntries[patchi]
-                );
-            }
-        }
-        else
-        {
-            FatalErrorIn
-            (
-                "referredCellList::buildReferredCellList(bool)"
-            )   << "unable to read undecomposed boundary file from "
-                << "constant/polyMesh" << nl
-                << abort(FatalError);
-        }
-
-        labelIOList faceProcAddressing
-        (
-            IOobject
-            (
-                "faceProcAddressing",
-                mesh.time().constant(),
-                polyMesh::meshSubDir,
-                mesh,
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE,
-                false
-            )
-        );
-
-        labelList procPatches(mesh.globalData().processorPatches());
-
-        nUndecomposedPatches = patchNames.size();
-
-        // processorPatchSegmentMapping works by mapping the original patch and
-        // half that a face on a processor patch was on before decomposition.
-        // This creates a patch segment for each half of each original (cyclic)
-        // patch which can be assessed separately.  There are n =
-        // patchNames.size() original patches, k = 0 to n-1.  The mapping is:
-        // value = 0: originally an internal face value = k, was originally on
-        // the on the patch k-1, 1st half value = -k, was originally on the on
-        // the patch k-1, 2nd half
-
-        forAll(procPatches,pP)
-        {
-            const processorPolyPatch& patch = refCast<const processorPolyPatch>
-            (
-                mesh.boundaryMesh()[procPatches[pP]]
-            );
-
-            labelList& procPatchSegMap = processorPatchSegmentMapping[pP];
-
-            procPatchSegMap.setSize(patch.size());
-
-            forAll(patch, pI)
-            {
-                label decomposedMeshFace = patch.start() + pI;
-
-                label faceProcAdd = faceProcAddressing[decomposedMeshFace];
-
-                label globalFace = abs(faceProcAdd)-1;
-
-                label minStart = -1;
-
-                forAll(patchNames, pN)
-                {
-                    if (patchDictionary.found(patchNames[pN]))
-                    {
-                        const dictionary& patchDict =
-                        patchDictionary.subDict(patchNames[pN]);
-
-                        label startFace
-                        (
-                            readLabel
-                            (
-                                patchDict.lookup("startFace")
-                            )
-                        );
-
-                        label nFaces(readLabel(patchDict.lookup("nFaces")));
-
-                        if (minStart < 0 || startFace < minStart)
-                        {
-                            minStart = startFace;
-                        }
-
-                        if
-                        (
-                            globalFace >= startFace
-                         && globalFace < startFace + nFaces/2
-                        )
-                        {
-                            procPatchSegMap[pI] = pN + 1;
-                        }
-                        else if
-                        (
-                            globalFace >= startFace + nFaces/2
-                         && globalFace < startFace + nFaces
-                        )
-                        {
-                            procPatchSegMap[pI] = -(pN + 1);
-                        }
-                    }
-                }
-
-                if (globalFace < minStart)
-                {
-                    procPatchSegMap[pI] = 0;
-                }
-            }
-        }
-
-        forAll(procPatches,pP)
-        {
-            const processorPolyPatch& patch = refCast<const processorPolyPatch>
-            (
-                mesh.boundaryMesh()[procPatches[pP]]
-            );
-
-            {
-                OPstream toNeighbProc
-                (
-                    Pstream::blocking,
-                    patch.neighbProcNo()
-                );
-
-                toNeighbProc << patch.faceCentres() << patch.faceAreas();
-            }
-        }
-
-        forAll(procPatches,pP)
-        {
-            const processorPolyPatch& patch = refCast<const processorPolyPatch>
-            (
-                mesh.boundaryMesh()[procPatches[pP]]
-            );
-
-            vectorList& neighbFaceCentres = allNeighbourFaceCentres[pP];
-
-            neighbFaceCentres.setSize(patch.size());
-
-            vectorList& neighbFaceAreas = allNeighbourFaceAreas[pP];
-
-            neighbFaceAreas.setSize(patch.size());
-
-            {
-                IPstream fromNeighbProc
-                (
-                    Pstream::blocking,
-                    patch.neighbProcNo()
-                );
-
-                fromNeighbProc >> neighbFaceCentres >> neighbFaceAreas;
-            }
-        }
-
-        // *************************************************************
-        // Tests that all processor patch segments from different
-        // original patches prior to decomposition produce the same
-        // transform. Check before 1st iteration.
-        // *************************************************************
-
-        forAll(procPatches,pP)
-        {
-            const processorPolyPatch& patch = refCast<const processorPolyPatch>
-            (
-                mesh.boundaryMesh()[procPatches[pP]]
-            );
-
-            const vectorList& neighbFaceCentres = allNeighbourFaceCentres[pP];
-
-            const vectorList& neighbFaceAreas = allNeighbourFaceAreas[pP];
-
-            label nUP;
-
-            for
-            (
-                nUP = -nUndecomposedPatches;
-                nUP <= nUndecomposedPatches;
-                nUP++
-            )
-            {
-                DynamicList<vector> refOff;
-
-                DynamicList<tensor> refTrans;
-
-                forAll(patch, faceI)
-                {
-                    if (processorPatchSegmentMapping[pP][faceI] == nUP)
-                    {
-                        referredCell testRefCell
-                        (
-                            mesh,
-                            -1,
-                            -1,
-                            patch.faceCentres()[faceI],
-                            neighbFaceCentres[faceI],
-                            patch.faceNormals()[faceI],
-                            neighbFaceAreas[faceI]
-                            /(mag(neighbFaceAreas[faceI]) + VSMALL)
-                        );
-
-                        refOff.append(testRefCell.offset());
-
-                        refTrans.append(testRefCell.rotation());
-                    }
-                }
-
-                refOff.shrink();
-
-                refTrans.shrink();
-
-                if (refOff.size())
-                {
-                    if
-                    (
-                        sum(mag(refOff-refOff[0]))/refOff.size()
-                            > interactionLists::transTol
-                     || sum(mag(refTrans-refTrans[0]))/refTrans.size()
-                            > interactionLists::transTol
-                    )
-                    {
-                        FatalErrorIn
-                        (
-                            "referredCellList::buildReferredCellList(bool)"
-                        )   << "Face pairs on patch " << patch.name()
-                            << ", segment " << patchNames[nUP]
-                            << " do not give the same referring "
-                            << " transformations to within tolerance of "
-                            << interactionLists::transTol << nl
-                            << " Referring offsets:" << refOff << nl
-                            << " Average sum of mag difference: "
-                            << sum(mag(refOff-refOff[0]))/refOff.size() << nl
-                            << " Referring transforms:" << refTrans << nl
-                            << " Average sum of mag difference: "
-                            << sum(mag(refTrans-refTrans[0]))/refTrans.size()
-                            << nl << abort(FatalError);
-                    }
-                }
-            }
-        }
-    }
-
-    label cellsReferredThisIteration = 1;
-
-    label iterationNo = 0;
-
-    while (cellsReferredThisIteration)
-    {
-        label refIntListStartSize = referredInteractionList.size();
-
-        forAll(mesh.boundaryMesh(), patchI)
-        {
-            // Treat local cyclics on each processor before processor
-            // boundaries.  Separate treatment allows the serial version to run
-            // transparently.
-
-            if (isA<cyclicPolyPatch>(mesh.boundaryMesh()[patchI]))
-            {
-                const cyclicPolyPatch& patch = refCast<const cyclicPolyPatch>
-                (
-                    mesh.boundaryMesh()[patchI]
-                );
-
-                if (patch.size())
-                {
-                    if (iterationNo == 0)
-                    {
-                        // Tests that all combinations of face pairs produce the
-                        // same transform.  Only check on 1st iteration.
-
-                        label faceL;
-                        // A face in the 1st half of the patch
-
-                        label faceM;
-                        // Face corresponding to faceL in the 2nd half of the
-                        // patch. Assumes correct face ordering.
-
-                        vectorList refOff(patch.size()/2);
-
-                        List<tensor> refTrans(patch.size()/2);
-
-                        for
-                        (
-                            faceL = 0, faceM = patch.size()/2;
-                            faceL < patch.size()/2;
-                            faceL++, faceM++
-                        )
-                        {
-                            referredCell testRefCell
-                            (
-                                mesh,
-                                -1,
-                                -1,
-                                patch.faceCentres()[faceL],
-                                patch.faceCentres()[faceM],
-                                patch.faceNormals()[faceL],
-                                patch.faceNormals()[faceM]
-                            );
-
-                            refOff[faceL] = testRefCell.offset();
-
-                            refTrans[faceL] = testRefCell.rotation();
-                        }
-
-                        if
-                        (
-                            sum(mag(refOff - refOff[0]))/(patch.size()/2)
-                                > interactionLists::transTol
-                         || sum(mag(refTrans - refTrans[0]))/(patch.size()/2)
-                                > interactionLists::transTol
-                        )
-                        {
-                            FatalErrorIn
-                            (
-                                "referredCellList::buildReferredCellList(bool)"
-                            )   << "Face pairs on patch " << patch.name()
-                                << " do not give the same referring "
-                                << " transformations to within tolerance of "
-                                << interactionLists::transTol << nl
-                                << " Referring offsets:" << refOff << nl
-                                << " Average sum of mag difference: "
-                                << sum(mag(refOff-refOff[0]))/refOff.size()
-                                << nl
-                                << " Referring transforms:" << refTrans << nl
-                                << " Average sum of mag difference: "
-                                << sum(mag(refTrans-refTrans[0]))
-                                   /refTrans.size()
-                                << nl << abort(FatalError);
-                        }
-                    }
-
-                    // *********************************************************
-                    // 1st half of face list - 1st side of boundary
-                    // *********************************************************
-
-                    label faceI;
-
-                    DynamicList<label> meshFacesOnThisSegment;
-
-                    for (faceI = 0; faceI < patch.size()/2; faceI++)
-                    {
-                        // unable to use the normal accessors for the polyPatch
-                        // because points on separate halves need used
-                        // separately
-
-                        meshFacesOnThisSegment.append(faceI + patch.start());
-                    }
-
-                    meshFacesOnThisSegment.shrink();
-
-                    DynamicList<label> meshEdgesOnThisSegment;
-
-                    DynamicList<label> meshPointsOnThisSegment;
-
-                    forAll(meshFacesOnThisSegment, mFOTS)
-                    {
-                        const label segFace = meshFacesOnThisSegment[mFOTS];
-
-                        const labelList& faceEdges = mesh.faceEdges()[segFace];
-
-                        forAll(faceEdges, fE)
-                        {
-                            const label faceEdge(faceEdges[fE]);
-
-                            if
-                            (
-                                findIndex
-                                (
-                                    meshEdgesOnThisSegment,
-                                    faceEdge
-                                ) == -1
-                            )
-                            {
-                                meshEdgesOnThisSegment.append(faceEdge);
-                            }
-                        }
-
-                        const face& facePoints(mesh.faces()[segFace]);
-
-                        forAll(facePoints, fP)
-                        {
-                            const label facePoint(facePoints[fP]);
-
-                            if
-                            (
-                                findIndex
-                                (
-                                    meshPointsOnThisSegment,
-                                    facePoint
-                                )
-                             ==
-                                -1
-                            )
-                            {
-                                meshPointsOnThisSegment.append(facePoint);
-                            }
-                        }
-                    }
-
-                    meshEdgesOnThisSegment.shrink();
-
-                    meshPointsOnThisSegment.shrink();
-
-                    if (iterationNo == 0)
-                    {
-                        // Assessing real cells in range is only required on
-                        // the 1st iteration because they do not change from
-                        // iteration to iteration.
-
-                        labelList realCellsFoundInRange
-                        (
-                            il_.realCellsInRangeOfSegment
-                            (
-                                meshFacesOnThisSegment,
-                                meshEdgesOnThisSegment,
-                                meshPointsOnThisSegment
-                            )
-                        );
-
-                        forAll(realCellsFoundInRange,cFIR)
-                        {
-                            const label realCell = realCellsFoundInRange[cFIR];
-
-                            referredCell cellToRefer
-                            (
-                                mesh,
-                                Pstream::myProcNo(),
-                                realCell,
-                                patch.faceCentres()[0],
-                                patch.faceCentres()[patch.size()/2],
-                                patch.faceNormals()[0],
-                                patch.faceNormals()[patch.size()/2]
-                            );
-
-                            // Test all existing referred and real cells to
-                            // check duplicates are not being made or cells
-                            // aren't being referred back onto themselves
-
-                            bool addCellToRefer = true;
-
-                            // Check if cellToRefer is an existing referred cell
-
-                            forAll(referredInteractionList, rIL)
-                            {
-                                if
-                                (
-                                    cellToRefer.duplicate
-                                    (
-                                        referredInteractionList[rIL]
-                                    )
-                                )
-                                {
-                                    addCellToRefer = false;
-
-                                    break;
-                                }
-                            }
-
-                            // Check for cellToRefer being referred back
-                            // ontop of a real cell
-
-                            if
-                            (
-                                cellToRefer.duplicate
-                                (
-                                    Pstream::myProcNo(),
-                                    mesh.nCells()
-                                )
-                            )
-                            {
-                                addCellToRefer = false;
-                            }
-
-                            if (addCellToRefer)
-                            {
-                                referredInteractionList.append(cellToRefer);
-                            }
-
-                            // add real cells found in range of cyclic patch
-                            // to whole mesh list
-
-                            if (findIndex (rCellsWRRP, realCell) == -1)
-                            {
-                                rCellsWRRP.append(realCell);
-                            }
-                        }
-                    }
-
-                    referredInteractionList.shrink();
-
-                    labelList referredCellsFoundInRange
-                    (
-                        il_.referredCellsInRangeOfSegment
-                        (
-                            referredInteractionList,
-                            meshFacesOnThisSegment,
-                            meshEdgesOnThisSegment,
-                            meshPointsOnThisSegment
-                        )
-                    );
-
-                    forAll(referredCellsFoundInRange,cFIR)
-                    {
-                        referredCell& existingRefCell =
-                            referredInteractionList
-                            [
-                                referredCellsFoundInRange[cFIR]
-                            ];
-
-                        referredCell cellToReRefer =
-                            existingRefCell.reRefer
-                            (
-                                patch.faceCentres()[0],
-                                patch.faceCentres()[patch.size()/2],
-                                patch.faceNormals()[0],
-                                patch.faceNormals()[patch.size()/2]
-                            );
-
-                        // Test all existing referred and real cells to check
-                        // duplicates are not being made or cells aren't being
-                        // referred back onto themselves
-
-                        bool addCellToReRefer = true;
-
-                        // Check if cellToRefer is an existing referred cell
-
-                        forAll(referredInteractionList, rIL)
-                        {
-                            if
-                            (
-                                cellToReRefer.duplicate
-                                (
-                                    referredInteractionList[rIL]
-                                )
-                            )
-                            {
-                                addCellToReRefer = false;
-
-                                break;
-                            }
-                        }
-
-                        // Check for cellToRefer being referred back
-                        // ontop of a real cell
-
-                        if
-                        (
-                            cellToReRefer.duplicate
-                            (
-                                Pstream::myProcNo(),
-                                mesh.nCells()
-                            )
-                        )
-                        {
-                            addCellToReRefer = false;
-                        }
-
-                        if (addCellToReRefer)
-                        {
-                            referredInteractionList.append(cellToReRefer);
-                        }
-                    }
-
-                    // *********************************************************
-                    // 2nd half of face list - 2nd side of boundary
-                    // *********************************************************
-
-                    meshFacesOnThisSegment.clear();
-
-                    for (faceI = patch.size()/2; faceI < patch.size(); faceI++)
-                    {
-                        // unable to use the normal accessors for the polyPatch
-                        // because points on separate halves need used
-                        // separately
-
-                        meshFacesOnThisSegment.append(faceI + patch.start());
-                    }
-
-                    meshFacesOnThisSegment.shrink();
-
-                    meshEdgesOnThisSegment.clear();
-
-                    meshPointsOnThisSegment.clear();
-
-                    forAll(meshFacesOnThisSegment, mFOTS)
-                    {
-                        const label segFace = meshFacesOnThisSegment[mFOTS];
-
-                        const labelList& faceEdges = mesh.faceEdges()[segFace];
-
-                        forAll(faceEdges, fE)
-                        {
-                            const label faceEdge(faceEdges[fE]);
-
-                            if
-                            (
-                                findIndex
-                                (
-                                    meshEdgesOnThisSegment,
-                                    faceEdge
-                                )
-                             ==
-                                -1
-                            )
-                            {
-                                meshEdgesOnThisSegment.append(faceEdge);
-                            }
-                        }
-
-                        const face& facePoints(mesh.faces()[segFace]);
-
-                        forAll(facePoints, fP)
-                        {
-                            const label facePoint(facePoints[fP]);
-
-                            if
-                            (
-                                findIndex
-                                (
-                                    meshPointsOnThisSegment,
-                                    facePoint
-                                )
-                             ==
-                                -1
-                            )
-                            {
-                                meshPointsOnThisSegment.append(facePoint);
-                            }
-                        }
-                    }
-
-                    meshEdgesOnThisSegment.shrink();
-
-                    meshPointsOnThisSegment.shrink();
-
-                    if (iterationNo == 0)
-                    {
-                        // Assessing real cells in range is only required on
-                        // the 1st iteration because they do not change from
-                        // iteration to iteration.
-
-                        labelList realCellsFoundInRange
-                        (
-                            il_.realCellsInRangeOfSegment
-                            (
-                                meshFacesOnThisSegment,
-                                meshEdgesOnThisSegment,
-                                meshPointsOnThisSegment
-                            )
-                        );
-
-                        forAll(realCellsFoundInRange,cFIR)
-                        {
-                            const label realCell = realCellsFoundInRange[cFIR];
-
-                            referredCell cellToRefer
-                            (
-                                mesh,
-                                Pstream::myProcNo(),
-                                realCell,
-                                patch.faceCentres()[patch.size()/2],
-                                patch.faceCentres()[0],
-                                patch.faceNormals()[patch.size()/2],
-                                patch.faceNormals()[0]
-                            );
-
-                            // Test all existing referred and real cells to
-                            // check duplicates are not being made or cells
-                            // aren't being referred back onto themselves
-
-                            bool addCellToRefer = true;
-
-                            // Check if cellToRefer is an existing referred cell
-
-                            forAll(referredInteractionList, rIL)
-                            {
-                                if
-                                (
-                                    cellToRefer.duplicate
-                                    (
-                                        referredInteractionList[rIL]
-                                    )
-                                )
-                                {
-                                    addCellToRefer = false;
-
-                                    break;
-                                }
-                            }
-
-                            // Check for cellToRefer being referred back
-                            // ontop of a real cell
-
-                            if
-                            (
-                                cellToRefer.duplicate
-                                (
-                                    Pstream::myProcNo(),
-                                    mesh.nCells()
-                                )
-                            )
-                            {
-                                addCellToRefer = false;
-                            }
-
-                            if (addCellToRefer)
-                            {
-                                referredInteractionList.append(cellToRefer);
-                            }
-
-                            // add real cells found in range of cyclic patch
-                            // to whole mesh list
-
-                            if (findIndex (rCellsWRRP, realCell) == -1)
-                            {
-                                rCellsWRRP.append(realCell);
-                            }
-                        }
-                    }
-
-                    referredInteractionList.shrink();
-
-                    referredCellsFoundInRange =
-                        il_.referredCellsInRangeOfSegment
-                        (
-                            referredInteractionList,
-                            meshFacesOnThisSegment,
-                            meshEdgesOnThisSegment,
-                            meshPointsOnThisSegment
-                        );
-
-                    forAll(referredCellsFoundInRange,cFIR)
-                    {
-                        referredCell& existingRefCell =
-                            referredInteractionList
-                            [
-                                referredCellsFoundInRange[cFIR]
-                            ];
-
-                        referredCell cellToReRefer =
-                            existingRefCell.reRefer
-                            (
-                                patch.faceCentres()[patch.size()/2],
-                                patch.faceCentres()[0],
-                                patch.faceNormals()[patch.size()/2],
-                                patch.faceNormals()[0]
-                            );
-
-                        // Test all existing referred and real cells to check
-                        // duplicates are not being made or cells aren't being
-                        // referred back onto themselves
-
-                        bool addCellToReRefer = true;
-
-                        // Check if cellToRefer is an existing referred cell
-
-                        forAll(referredInteractionList, rIL)
-                        {
-                            if
-                            (
-                                cellToReRefer.duplicate
-                                (
-                                    referredInteractionList[rIL]
-                                )
-                            )
-                            {
-                                addCellToReRefer = false;
-
-                                break;
-                            }
-                        }
-
-                        // Check for cellToRefer being referred back
-                        // ontop of a real cell
-
-                        if
-                        (
-                            cellToReRefer.duplicate
-                            (
-                                Pstream::myProcNo(),
-                                mesh.nCells()
-                            )
-                        )
-                        {
-                            addCellToReRefer = false;
-                        }
-
-                        if (addCellToReRefer)
-                        {
-                            referredInteractionList.append(cellToReRefer);
-                        }
-                    }
-                }
-            }
-        }
-
-        if (Pstream::parRun())
-        {
-            labelList procPatches(mesh.globalData().processorPatches());
-
-            forAll(procPatches,pP)
-            {
-                const processorPolyPatch& patch =
-                    refCast<const processorPolyPatch>
-                    (
-                        mesh.boundaryMesh()[procPatches[pP]]
-                    );
-
-                DynamicList<referredCell> referredCellsToTransfer;
-
-                const vectorList& neighbFaceCentres =
-                    allNeighbourFaceCentres[pP];
-
-                const vectorList& neighbFaceAreas = allNeighbourFaceAreas[pP];
-
-                label nUP;
-
-                for
-                (
-                    nUP = -nUndecomposedPatches;
-                    nUP <= nUndecomposedPatches;
-                    nUP++
-                )
-                {
-                    // faceT is used to specify one face on this patch segment
-                    // that will be used to calculate the transformation values.
-                    // All faces are guaranteed to produce the same transform
-                    // because of the checks carried out at the start of the
-                    // function.  Setting to -1 until the 1st face on this
-                    // segment is found.
-
-                    label faceT = -1;
-
-                    DynamicList<label> meshFacesOnThisSegment;
-
-                    forAll(patch, faceI)
-                    {
-                        if (processorPatchSegmentMapping[pP][faceI] == nUP)
-                        {
-                            if (faceT == -1)
-                            {
-                                faceT = faceI;
-                            }
-
-                            meshFacesOnThisSegment.append
-                            (
-                                faceI + patch.start()
-                            );
-                        }
-                    }
-
-                    meshFacesOnThisSegment.shrink();
-
-                    DynamicList<label> meshEdgesOnThisSegment;
-
-                    DynamicList<label> meshPointsOnThisSegment;
-
-                    forAll(meshFacesOnThisSegment, mFOTS)
-                    {
-                        const label segFace = meshFacesOnThisSegment[mFOTS];
-
-                        const labelList& faceEdges = mesh.faceEdges()[segFace];
-
-                        forAll(faceEdges, fE)
-                        {
-                            const label faceEdge(faceEdges[fE]);
-
-                            if
-                            (
-                                findIndex
-                                (
-                                    meshEdgesOnThisSegment,
-                                    faceEdge
-                                )
-                             ==
-                                -1
-                            )
-                            {
-                                meshEdgesOnThisSegment.append(faceEdge);
-                            }
-                        }
-
-                        const face& facePoints(mesh.faces()[segFace]);
-
-                        forAll(facePoints, fP)
-                        {
-                            const label facePoint(facePoints[fP]);
-
-                            if
-                            (
-                                findIndex
-                                (
-                                    meshPointsOnThisSegment,
-                                    facePoint
-                                )
-                             ==
-                                -1
-                            )
-                            {
-                                meshPointsOnThisSegment.append(facePoint);
-                            }
-                        }
-                    }
-
-                    meshEdgesOnThisSegment.shrink();
-
-                    meshPointsOnThisSegment.shrink();
-
-                    if (meshFacesOnThisSegment.size())
-                    {
-                        if (faceT == -1)
-                        {
-                            FatalErrorIn
-                            (
-                                "referredCellList::buildReferredCellList(bool)"
-                            )   << "faceT == -1 encountered but "
-                                << meshFacesOnThisSegment.size()
-                                << " faces found on patch segment."
-                                << abort(FatalError);
-                        }
-
-                        if (iterationNo == 0)
-                        {
-                            // Assessing real cells in range is only required on
-                            // the 1st iteration because they do not change from
-                            // iteration to iteration.
-
-                            labelList realCellsFoundInRange
-                            (
-                                il_.realCellsInRangeOfSegment
-                                (
-                                    meshFacesOnThisSegment,
-                                    meshEdgesOnThisSegment,
-                                    meshPointsOnThisSegment
-                                )
-                            );
-
-                            forAll(realCellsFoundInRange,cFIR)
-                            {
-                                const label realCell =
-                                    realCellsFoundInRange[cFIR];
-
-                                referredCell cellToRefer
-                                (
-                                    mesh,
-                                    Pstream::myProcNo(),
-                                    realCell,
-                                    patch.faceCentres()[faceT],
-                                    neighbFaceCentres[faceT],
-                                    patch.faceNormals()[faceT],
-                                    neighbFaceAreas[faceT]
-                                    /(mag(neighbFaceAreas[faceT]) + VSMALL)
-                                );
-
-                                referredCellsToTransfer.append(cellToRefer);
-
-                                // add real cells found in range of processor
-                                // patch to whole mesh list
-
-                                if (findIndex (rCellsWRRP, realCell) == -1)
-                                {
-                                    rCellsWRRP.append(realCell);
-                                }
-                            }
-                        }
-
-                        referredInteractionList.shrink();
-
-                        labelList referredCellsFoundInRange
-                        (
-                            il_.referredCellsInRangeOfSegment
-                            (
-                                referredInteractionList,
-                                meshFacesOnThisSegment,
-                                meshEdgesOnThisSegment,
-                                meshPointsOnThisSegment
-                            )
-                        );
-
-                        forAll(referredCellsFoundInRange,cFIR)
-                        {
-                            referredCell& existingRefCell =
-                                referredInteractionList
-                                [
-                                    referredCellsFoundInRange[cFIR]
-                                ];
-
-                            referredCell cellToReRefer =
-                                existingRefCell.reRefer
-                                (
-                                    patch.faceCentres()[faceT],
-                                    neighbFaceCentres[faceT],
-                                    patch.faceNormals()[faceT],
-                                    neighbFaceAreas[faceT]
-                                    /(mag(neighbFaceAreas[faceT]) + VSMALL)
-                                );
-
-                            referredCellsToTransfer.append(cellToReRefer);
-                        }
-                    }
-                }
-
-                referredCellsToTransfer.shrink();
-
-                // Send these cells to the neighbouring processor.
-
-                {
-                    OPstream toNeighbProc
-                    (
-                        Pstream::blocking,
-                        patch.neighbProcNo()
-                    );
-
-                    toNeighbProc << referredCellsToTransfer;
-                }
-            }
-
-            forAll(procPatches,pP)
-            {
-                const processorPolyPatch& patch =
-                refCast<const processorPolyPatch>
-                (
-                    mesh.boundaryMesh()[procPatches[pP]]
-                );
-
-                // Receive the cells from neighbour
-
-                List<referredCell> referredCellsFromNeighbour(patch.size());
-
-                {
-                    IPstream fromNeighbProc
-                    (
-                        Pstream::blocking,
-                        patch.neighbProcNo()
-                    );
-
-                    fromNeighbProc >> referredCellsFromNeighbour;
-                }
-
-                // Check to see if they are duplicates, if not append
-                // them to the referredInteractionList
-
-                forAll(referredCellsFromNeighbour,rCFN)
-                {
-                    referredCell& cellToRefer =
-                    referredCellsFromNeighbour[rCFN];
-
-                    // Test all existing referred and real cells to check
-                    // duplicates are not being made or cells aren't being
-                    // referred back onto themselves
-
-                    bool addCellToRefer = true;
-
-                    // Check if cellToRefer is an existing referred cell
-
-                    forAll(referredInteractionList, rIL)
-                    {
-                        if (cellToRefer.duplicate(referredInteractionList[rIL]))
-                        {
-                            addCellToRefer = false;
-
-                            break;
-                        }
-                    }
-
-                    // Check for cellToRefer being referred back ontop of a real
-                    // cell
-
-                    if
-                    (
-                        cellToRefer.duplicate
-                        (
-                            Pstream::myProcNo(),
-                            mesh.nCells()
-                        )
-                    )
-                    {
-                        addCellToRefer = false;
-                    }
-
-                    if (addCellToRefer)
-                    {
-                        referredInteractionList.append(cellToRefer);
-                    }
-                }
-            }
-        }
-
-        if (iterationNo == 0)
-        {
-            // record all real cells in range of any referring patch (cyclic or
-            // processor) on the first iteration when the real cells are
-            // evaluated.
-
-            rCellsWRRP.shrink();
-
-            // construct {faces, edges, points}WithinRCutMaxOfAnyReferringPatch
-
-            forAll(rCellsWRRP, rCWR)
-            {
-                const label realCell(rCellsWRRP[rCWR]);
-
-                const labelList& rCFaces
-                (
-                    mesh.cells()[realCell]
-                );
-
-                forAll(rCFaces, rCF)
-                {
-                    const label f(rCFaces[rCF]);
-
-                    if (findIndex(rFacesWRRP,f) == -1)
-                    {
-                        rFacesWRRP.append(f);
-                    }
-                }
-
-                const labelList& rCEdges
-                (
-                    mesh.cellEdges()[realCell]
-                );
-
-                forAll(rCEdges, rCE)
-                {
-                    const label e(rCEdges[rCE]);
-
-                    if (findIndex(rEdgesWRRP,e) == -1)
-                    {
-                        rEdgesWRRP.append(e);
-                    }
-                }
-
-                const labelList& rCPoints
-                (
-                    mesh.cellPoints()[realCell]
-                );
-
-                forAll(rCPoints, rCP)
-                {
-                    const label p(rCPoints[rCP]);
-
-                    if (findIndex(rPointsWRRP,p) == -1)
-                    {
-                        rPointsWRRP.append(p);
-                    }
-                }
-            }
-
-            rFacesWRRP.shrink();
-
-            rEdgesWRRP.shrink();
-
-            rPointsWRRP.shrink();
-        }
-
-        iterationNo++;
-
-        cellsReferredThisIteration =
-            referredInteractionList.size() - refIntListStartSize;
-
-        reduce(cellsReferredThisIteration, sumOp<label>());
-
-        Info<< tab << "Cells added this iteration: "
-            << cellsReferredThisIteration << endl;
-    }
-
-    referredInteractionList.shrink();
-
-    //     Info<< "referredInteractionList.size() = "
-    //         << referredInteractionList.size() << endl;
-
-    //     forAll(referredInteractionList,rIL)
-    //     {
-    //         Info<< referredInteractionList[rIL];
-    //     }
-
-    (*this).setSize
-    (
-        referredInteractionList.size()
-    );
-
-    forAll(referredInteractionList, rIL)
-    {
-        (*this)[rIL] = referredInteractionList[rIL];
-    }
-
-    Info<< nl << "Finding real cells in range of referred cells" << endl;
-
-    forAll(*this, rC)
-    {
-        const polyMesh& mesh(il_.mesh());
-
-        referredCell& refCell = (*this)[rC];
-
-        DynamicList<label> realCellsFoundInRange;
-
-        const vectorList& refCellPoints = refCell.vertexPositions();
-
-        forAll(rFacesWRRP, rCF)
-        {
-            const label f(rFacesWRRP[rCF]);
-
-            if (il_.testPointFaceDistance(refCellPoints,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(rPointsWRRP, rCP)
-        {
-            const label p(rPointsWRRP[rCP]);
-
-            if (il_.testPointFaceDistance(p,refCell))
-            {
-                const labelList& pCells(mesh.pointCells()[p]);
-
-                forAll(pCells, pC)
-                {
-                    const label cellI(pCells[pC]);
-
-                    if (findIndex(realCellsFoundInRange, cellI) == -1)
-                    {
-                        realCellsFoundInRange.append(cellI);
-                    }
-                }
-            }
-        }
-
-
-        const edgeList& refCellEdges = refCell.edges();
-
-        forAll(rEdgesWRRP, rCE)
-        {
-            const label edgeIIndex(rEdgesWRRP[rCE]);
-
-            const edge& eI(mesh.edges()[edgeIIndex]);
-
-            forAll(refCellEdges, rCE)
-            {
-                const edge& eJ(refCellEdges[rCE]);
-
-                if
-                (
-                    il_.testEdgeEdgeDistance
-                    (
-                        eI,
-                        refCellPoints[eJ.start()],
-                        refCellPoints[eJ.end()]
-                    )
-                )
-                {
-                    const labelList& eICells(mesh.edgeCells()[edgeIIndex]);
-
-                    forAll(eICells, eIC)
-                    {
-                        const label cellI(eICells[eIC]);
-
-                        if (findIndex(realCellsFoundInRange, cellI) == -1)
-                        {
-                            realCellsFoundInRange.append(cellI);
-                        }
-                    }
-                }
-            }
-        }
-
-//         scalar rCutMaxSqr = molCloud_.rCutMax()*molCloud_.rCutMax();
-//
-//         forAll(molCloud_.mesh().points(), pointIIndex)
-//         {
-//             const point& ptI
-//             (
-//                 molCloud_.mesh().points()[pointIIndex]
-//             );
-//
-//             forAll(refCellPoints, rCP)
-//             {
-//                 if (magSqr(ptI - refCellPoints[rCP]) <= rCutMaxSqr)
-//                 {
-//                     const labelList& ptICells
-//                     (
-//                         molCloud_.mesh().pointCells()[pointIIndex]
-//                     );
-//
-//                     forAll(ptICells, pIC)
-//                     {
-//                         const label cellI(ptICells[pIC]);
-//
-//                         if (findIndex(realCellsFoundInRange, cellI) == -1)
-//                         {
-//                             realCellsFoundInRange.append(cellI);
-//                         }
-//                     }
-//                 }
-//             }
-//         }
-
-        refCell.realCells() = realCellsFoundInRange.shrink();
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::referredCellList::referredCellList
-(
-    interactionLists& il,
-    bool pointPointListBuild
-)
-:
-    List<referredCell>(),
-    il_(il)
-{
-    buildReferredCellList(pointPointListBuild);
-}
-
-
-Foam::referredCellList::referredCellList(interactionLists& il)
-:
-    List<referredCell>(),
-    il_(il)
-{
-    Info<< "Read referredCellList from disk not implemented" << endl;
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::referredCellList::~referredCellList()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::referredCellList::referMolecules
-(
-    const List<DynamicList<molecule*> >& cellOccupancy
-)
-{
-    // Create referred molecules for sending using cell occupancy and
-    // cellSendingReferralLists
-
-    forAll(il_.cellSendingReferralLists(), cSRL)
-    {
-        const sendingReferralList& sRL
-        (
-            il_.cellSendingReferralLists()[cSRL]
-        );
-
-        List<DynamicList<referredMolecule> > molsToReferOut(sRL.size());
-
-        forAll(sRL, sRLI)
-        {
-            List<molecule*> realMols = cellOccupancy[sRL[sRLI]];
-
-            forAll(realMols, rM)
-            {
-                molecule* mol = realMols[rM];
-
-                molsToReferOut[sRLI].append
-                (
-                    referredMolecule
-                    (
-                        mol->id(),
-                        mol->position(),
-                        mol->sitePositions()
-                    )
-                );
-            }
-
-            molsToReferOut[sRLI].shrink();
-        }
-
-        // Send lists of referred molecules to other processors
-
-        if (sRL.destinationProc() != Pstream::myProcNo())
-        {
-            if (Pstream::parRun())
-            {
-                OPstream toInteractingProc
-                (
-                    Pstream::blocking,
-                    sRL.destinationProc()
-                );
-
-                toInteractingProc << molsToReferOut;
-            }
-        }
-        else
-        {
-            // Refer molecules directly for referred cells on the same
-            // processor.
-
-            const receivingReferralList& rRL
-            (
-                il_.cellReceivingReferralLists()[cSRL]
-            );
-
-            forAll(rRL, rRLI)
-            {
-                forAll(rRL[rRLI], rC)
-                {
-                    referredCell& refCellToRefMolsTo = (*this)[rRL[rRLI][rC]];
-
-                    refCellToRefMolsTo.referInMols(molsToReferOut[rRLI]);
-                }
-            }
-        }
-    }
-
-    // Receive referred molecule lists to and distribute to referredCells
-    // according tocellReceivingReferralLists, referredCells deal with the
-    // transformations themselves
-
-    forAll(il_.cellReceivingReferralLists(), cRRL)
-    {
-        const receivingReferralList& rRL
-        (
-            il_.cellReceivingReferralLists()[cRRL]
-        );
-
-        List<List<referredMolecule> > molsToReferIn(rRL.size());
-
-        if (rRL.sourceProc() != Pstream::myProcNo())
-        {
-            if (Pstream::parRun())
-            {
-                IPstream fromInteractingProc
-                (
-                    Pstream::blocking,
-                    rRL.sourceProc()
-                );
-
-                fromInteractingProc >> molsToReferIn;
-            }
-
-            forAll(rRL, rRLI)
-            {
-                forAll(rRL[rRLI], rC)
-                {
-                    referredCell& refCellToRefMolsTo = (*this)[rRL[rRLI][rC]];
-
-                    refCellToRefMolsTo.referInMols(molsToReferIn[rRLI]);
-                }
-            }
-        }
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCellList/referredCellList.H b/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCellList/referredCellList.H
deleted file mode 100644
index 5f6f8284c7..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCellList/referredCellList.H
+++ /dev/null
@@ -1,109 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::referredCellList
-
-Description
-
-SourceFiles
-    referredCellListI.H
-    referredCellList.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef referredCellList_H
-#define referredCellList_H
-
-#include "referredCell.H"
-#include "molecule.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-class interactionLists;
-
-/*---------------------------------------------------------------------------*\
-                      Class referredCellList Declaration
-\*---------------------------------------------------------------------------*/
-
-class referredCellList
-:
-    public List<referredCell>
-{
-    // Private data
-
-        const interactionLists& il_;
-
-
-    // Private Member Functions
-
-        void buildReferredCellList
-        (
-            bool pointPointListBuild
-        );
-
-
-public:
-
-    // Constructors
-
-        //- Construct lists by searching the mesh
-        referredCellList
-        (
-            interactionLists& il,
-            bool pointPointListBuild
-        );
-
-        //- Construct from file
-        referredCellList (interactionLists& il);
-
-
-    //- Destructor
-    ~referredCellList();
-
-
-    // Member Functions
-
-        void referMolecules(const List<DynamicList<molecule*> >& cellOccupancy);
-
-        inline const interactionLists& il() const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "referredCellListI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCellList/referredCellListI.H b/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCellList/referredCellListI.H
deleted file mode 100644
index 624667100f..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCellList/referredCellListI.H
+++ /dev/null
@@ -1,34 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-inline const Foam::interactionLists& Foam::referredCellList::il() const
-{
-    return il_;
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredMolecule/referredMolecule.C b/src/lagrangian/molecularDynamics/molecule/interactionLists/referredMolecule/referredMolecule.C
deleted file mode 100644
index de4428c3b2..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredMolecule/referredMolecule.C
+++ /dev/null
@@ -1,85 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*----------------------------------------------------------------------------*/
-
-#include "referredMolecule.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::referredMolecule::referredMolecule()
-{}
-
-
-Foam::referredMolecule::referredMolecule
-(
-    const label id,
-    const vector& position,
-    const List<vector>& sitePositions
-)
-:
-    id_(id),
-    position_(position),
-    sitePositions_(sitePositions)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::referredMolecule::~referredMolecule()
-{}
-
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
-Foam::Istream& Foam::operator>>
-(
-    Istream& is,
-    referredMolecule& rM
-)
-{
-    is >> rM.id_ >> rM.position_ >> rM.sitePositions_;
-
-    is.check("Istream& operator<<(Istream& f, const referredMolecule& sRL");
-
-    return is;
-}
-
-
-Foam::Ostream& Foam::operator<<
-(
-    Ostream& os,
-    const referredMolecule& rM
-)
-{
-    os  << rM.id()
-        << token::SPACE << rM.position()
-        << token::SPACE << rM.sitePositions();
-
-    os.check("Ostream& operator<<(Ostream& f, const referredMolecule& rM");
-
-    return os;
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredMolecule/referredMolecule.H b/src/lagrangian/molecularDynamics/molecule/interactionLists/referredMolecule/referredMolecule.H
deleted file mode 100644
index c55c3ba263..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredMolecule/referredMolecule.H
+++ /dev/null
@@ -1,135 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::referredMolecule
-
-Description
-
-SourceFiles
-    referredMoleculeI.H
-    referredMolecule.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef referredMolecule_H
-#define referredMolecule_H
-
-#include "vector.H"
-#include "List.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                      Class referredMolecule Declaration
-\*---------------------------------------------------------------------------*/
-
-class referredMolecule
-{
-    // Private data
-
-        label id_;
-
-        vector position_;
-
-        List<vector> sitePositions_;
-
-
-public:
-
-    // Constructors
-
-        //- Construct null
-        referredMolecule();
-
-        //- Construct from components
-        referredMolecule
-        (
-            const label id,
-            const vector& position,
-            const List<vector>& sitePositions
-        );
-
-
-    //- Destructor
-    virtual ~referredMolecule();
-
-
-    // Member Functions
-
-        // Access
-
-            inline label id() const;
-
-            inline const vector& position() const;
-
-            inline const List<vector>& sitePositions() const;
-
-
-    // Friend Operators
-
-        inline friend bool operator==
-        (
-            const referredMolecule& a,
-            const referredMolecule& b
-        );
-
-        inline friend bool operator!=
-        (
-            const referredMolecule& a,
-            const referredMolecule& b
-        );
-
-
-    // IOstream Operators
-
-        friend Istream& operator>>
-        (
-            Istream&,
-            referredMolecule&
-        );
-
-        friend Ostream& operator<<
-        (
-            Ostream&,
-            const referredMolecule&
-        );
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "referredMoleculeI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredMolecule/referredMoleculeI.H b/src/lagrangian/molecularDynamics/molecule/interactionLists/referredMolecule/referredMoleculeI.H
deleted file mode 100644
index 317b8cca30..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredMolecule/referredMoleculeI.H
+++ /dev/null
@@ -1,74 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-inline Foam::label Foam::referredMolecule::id() const
-{
-    return id_;
-}
-
-
-inline const Foam::vector& Foam::referredMolecule::position() const
-{
-    return position_;
-}
-
-
-inline const Foam::List<Foam::vector>&
-Foam::referredMolecule::sitePositions() const
-{
-    return sitePositions_;
-}
-
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
-inline bool Foam::operator==
-(
-    const Foam::referredMolecule& a,
-    const Foam::referredMolecule& b
-)
-{
-    return
-    (
-        a.id() == b.id()
-     && a.position() == b.position()
-    );
-}
-
-
-inline bool Foam::operator!=
-(
-    const Foam::referredMolecule& a,
-    const Foam::referredMolecule& b
-)
-{
-    return !(a == b);
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
index 9dca0cbc40..2af297c6ac 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
@@ -209,8 +209,20 @@ bool Foam::molecule::move(molecule::trackData& td)
 
 void Foam::molecule::transformProperties(const tensor& T)
 {
+    Particle<molecule>::transformProperties(T);
+
     Q_ = T & Q_;
 
+    v_ = transform(T, v_);
+
+    a_ = transform(T, a_);
+
+    pi_ = Q_.T() & transform(T, Q_ & pi_);
+
+    tau_ = Q_.T() & transform(T, Q_ & tau_);
+
+    rf_ = transform(T, rf_);
+
     sitePositions_ = position_ + (T & (sitePositions_ - position_));
 
     siteForces_ = T & siteForces_;
@@ -219,10 +231,14 @@ void Foam::molecule::transformProperties(const tensor& T)
 
 void Foam::molecule::transformProperties(const vector& separation)
 {
+    Particle<molecule>::transformProperties(separation);
+
     if (special_ == SPECIAL_TETHERED)
     {
         specialPosition_ += separation;
     }
+
+    sitePositions_ = sitePositions_ + separation;
 }
 
 
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
index 8671356ffd..448cb0955b 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
@@ -182,13 +182,6 @@ private:
 
     // Private data
 
-        //- Be careful with the ordering of data.
-        //  It has an impact on binary transfer:
-        //    -# Put the largest data members 1st
-        //    -# Pair up labels,
-        //    -# Don't go scalar-label, scalar-label, because in 64bit mode,
-        //       the labels will be padded by 4bytes.
-
         tensor Q_;
 
         vector v_;
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C
index 7c8884ae31..e7ac44a972 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C
@@ -242,11 +242,17 @@ void Foam::molecule::writeFields(const Cloud<molecule>& mC)
     orientation2.write();
     orientation3.write();
 
-    const moleculeCloud& m = dynamic_cast<const moleculeCloud&>(mC);
-    m.writeXYZ
-    (
-        m.mesh().time().timePath() + "/lagrangian" + "/moleculeCloud.xmol"
-    );
+    Info<< "writeFields " << mC.name() << endl;
+
+    if (isA<moleculeCloud>(mC))
+    {
+        const moleculeCloud& m = dynamic_cast<const moleculeCloud&>(mC);
+
+        m.writeXYZ
+        (
+            m.mesh().time().timePath()/cloud::prefix/"moleculeCloud.xmol"
+        );
+    }
 }
 
 
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
index 58b6b0603e..55f24f100f 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
@@ -123,23 +123,23 @@ void Foam::moleculeCloud::buildCellOccupancy()
     {
         cellOccupancy_[cO].shrink();
     }
-
-    il_.ril().referMolecules(cellOccupancy());
 }
 
 
 void Foam::moleculeCloud::calculatePairForce()
 {
-    iterator mol(this->begin());
+    PstreamBuffers pBufs(Pstream::nonBlocking);
+
+    // Start sending referred data
+    il_.sendReferredData(cellOccupancy(), pBufs);
+
+    molecule* molI = NULL;
+    molecule* molJ = NULL;
 
     {
         // Real-Real interactions
 
-        molecule* molI = &mol();
-
-        molecule* molJ = &mol();
-
-        const directInteractionList& dil(il_.dil());
+        const labelListList& dil = il_.dil();
 
         forAll(dil, d)
         {
@@ -149,54 +149,62 @@ void Foam::moleculeCloud::calculatePairForce()
 
                 forAll(dil[d], interactingCells)
                 {
-                    List< molecule* > cellJ =
-                    cellOccupancy_[dil[d][interactingCells]];
+                    List<molecule*> cellJ =
+                        cellOccupancy_[dil[d][interactingCells]];
 
                     forAll(cellJ, cellJMols)
                     {
                         molJ = cellJ[cellJMols];
 
-                        evaluatePair(molI, molJ);
+                        evaluatePair(*molI, *molJ);
                     }
                 }
 
-                forAll(cellOccupancy_[d],cellIOtherMols)
+                forAll(cellOccupancy_[d], cellIOtherMols)
                 {
                     molJ = cellOccupancy_[d][cellIOtherMols];
 
                     if (molJ > molI)
                     {
-                        evaluatePair(molI, molJ);
+                        evaluatePair(*molI, *molJ);
                     }
                 }
             }
         }
     }
 
+    // Receive referred data
+    il_.receiveReferredData(pBufs);
+
     {
         // Real-Referred interactions
 
-        molecule* molK = &mol();
+        const labelListList& ril = il_.ril();
 
-        referredCellList& ril(il_.ril());
+        List<IDLList<molecule> >& referredMols = il_.referredParticles();
 
         forAll(ril, r)
         {
-            const List<label>& realCells = ril[r].realCellsForInteraction();
+            const List<label>& realCells = ril[r];
 
-            forAll(ril[r], refMols)
+            IDLList<molecule>& refMols = referredMols[r];
+
+            forAllIter
+            (
+                IDLList<molecule>,
+                refMols,
+                refMol
+            )
             {
-                referredMolecule* molL(&(ril[r][refMols]));
-
                 forAll(realCells, rC)
                 {
-                    List<molecule*> cellK = cellOccupancy_[realCells[rC]];
+                    List<molecule*> cellI = cellOccupancy_[realCells[rC]];
 
-                    forAll(cellK, cellKMols)
+                    forAll(cellI, cellIMols)
                     {
-                        molK = cellK[cellKMols];
+                        molI = cellI[cellIMols];
 
-                        evaluatePair(molK, molL);
+                        evaluatePair(*molI, refMol());
                     }
                 }
             }
@@ -259,18 +267,14 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
     buildCellOccupancy();
 
     // Real-Real interaction
-    iterator mol(this->begin());
+
+    molecule* molI = NULL;
+    molecule* molJ = NULL;
 
     {
-        mol = this->begin();
-
-        molecule* molI = &mol();
-
-        molecule* molJ = &mol();
-
         DynamicList<molecule*> molsToDelete;
 
-        const directInteractionList& dil(il_.dil());
+        const labelListList& dil(il_.dil());
 
         forAll(dil, d)
         {
@@ -280,14 +284,14 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
 
                 forAll(dil[d], interactingCells)
                 {
-                    List< molecule* > cellJ =
-                    cellOccupancy_[dil[d][interactingCells]];
+                    List<molecule*> cellJ =
+                        cellOccupancy_[dil[d][interactingCells]];
 
                     forAll(cellJ, cellJMols)
                     {
                         molJ = cellJ[cellJMols];
 
-                        if (evaluatePotentialLimit(molI, molJ))
+                        if (evaluatePotentialLimit(*molI, *molJ))
                         {
                             label idI = molI->id();
 
@@ -297,7 +301,7 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
                             (
                                 idI == idJ
                              || findIndex(pot_.removalOrder(), idJ)
-                                < findIndex(pot_.removalOrder(), idI)
+                              < findIndex(pot_.removalOrder(), idI)
                             )
                             {
                                 if (findIndex(molsToDelete, molJ) == -1)
@@ -314,13 +318,13 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
                 }
             }
 
-            forAll(cellOccupancy_[d],cellIOtherMols)
+            forAll(cellOccupancy_[d], cellIOtherMols)
             {
                 molJ = cellOccupancy_[d][cellIOtherMols];
 
                 if (molJ > molI)
                 {
-                    if (evaluatePotentialLimit(molI, molJ))
+                    if (evaluatePotentialLimit(*molI, *molJ))
                     {
                         label idI = molI->id();
 
@@ -330,7 +334,7 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
                         (
                             idI == idJ
                          || findIndex(pot_.removalOrder(), idJ)
-                            < findIndex(pot_.removalOrder(), idI)
+                          < findIndex(pot_.removalOrder(), idI)
                         )
                         {
                             if (findIndex(molsToDelete, molJ) == -1)
@@ -355,79 +359,81 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
 
     buildCellOccupancy();
 
+    PstreamBuffers pBufs(Pstream::nonBlocking);
+
+    // Start sending referred data
+    il_.sendReferredData(cellOccupancy(), pBufs);
+
+        // Receive referred data
+    il_.receiveReferredData(pBufs);
+
     // Real-Referred interaction
 
     {
-        molecule* molK = &mol();
-
         DynamicList<molecule*> molsToDelete;
 
-        referredCellList& ril(il_.ril());
+        const labelListList& ril(il_.ril());
+
+        List<IDLList<molecule> >& referredMols = il_.referredParticles();
 
         forAll(ril, r)
         {
-            referredCell& refCell = ril[r];
+            IDLList<molecule>& refMols = referredMols[r];
 
-            forAll(refCell, refMols)
+            forAllIter
+            (
+                IDLList<molecule>,
+                refMols,
+                refMol
+            )
             {
-                referredMolecule* molL = &(refCell[refMols]);
+                molJ = &refMol();
 
-                List <label> realCells = refCell.realCellsForInteraction();
+                const List<label>& realCells = ril[r];
 
                 forAll(realCells, rC)
                 {
-                    label cellK = realCells[rC];
+                    label cellI = realCells[rC];
 
-                    List<molecule*> cellKMols = cellOccupancy_[cellK];
+                    List<molecule*> cellIMols = cellOccupancy_[cellI];
 
-                    forAll(cellKMols, cKM)
+                    forAll(cellIMols, cIM)
                     {
-                        molK = cellKMols[cKM];
+                        molI = cellIMols[cIM];
 
-                        if (evaluatePotentialLimit(molK, molL))
+                        if (evaluatePotentialLimit(*molI, *molJ))
                         {
-                            label idK = molK->id();
+                            label idI = molI->id();
 
-                            label idL = molL->id();
+                            label idJ = molJ->id();
 
                             if
                             (
-                                findIndex(pot_.removalOrder(), idK)
-                                    < findIndex(pot_.removalOrder(), idL)
+                                findIndex(pot_.removalOrder(), idI)
+                              < findIndex(pot_.removalOrder(), idJ)
                             )
                             {
-                                if (findIndex(molsToDelete, molK) == -1)
+                                if (findIndex(molsToDelete, molI) == -1)
                                 {
-                                    molsToDelete.append(molK);
+                                    molsToDelete.append(molI);
                                 }
                             }
-
                             else if
                             (
-                                findIndex(pot_.removalOrder(), idK)
-                                    == findIndex(pot_.removalOrder(), idL)
+                                findIndex(pot_.removalOrder(), idI)
+                             == findIndex(pot_.removalOrder(), idJ)
                             )
                             {
-                                if
-                                (
-                                    Pstream::myProcNo() == refCell.sourceProc()
-                                 && cellK <= refCell.sourceCell()
-                                )
-                                {
-                                    if (findIndex(molsToDelete, molK) == -1)
-                                    {
-                                        molsToDelete.append(molK);
-                                    }
-                                }
+                                // Remove one of the molecules
+                                // arbitrarily, assuring that a
+                                // consistent decision is made for
+                                // both real-referred pairs.
 
-                                else if
-                                (
-                                    Pstream::myProcNo() < refCell.sourceProc()
-                                )
+                                if (molI->origId() > molJ->origId())
                                 {
-                                    if (findIndex(molsToDelete, molK) == -1)
+                                    if (findIndex(molsToDelete, molI) == -1)
                                     {
-                                        molsToDelete.append(molK);
+                                        molsToDelete.append(molI);
                                     }
                                 }
                             }
@@ -445,6 +451,12 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
 
     buildCellOccupancy();
 
+    // Start sending referred data
+    il_.sendReferredData(cellOccupancy(), pBufs);
+
+    // Receive referred data
+    il_.receiveReferredData(pBufs);
+
     label molsRemoved = initialSize - this->size();
 
     if (Pstream::parRun())
@@ -713,19 +725,22 @@ void Foam::moleculeCloud::initialiseMolecules
                             label id = findIndex(pot_.idList(), latticeIds[p]);
 
                             const vector& latticePosition =
-                            vector
-                            (
-                                unitCellLatticePosition.x(),
-                                unitCellLatticePosition.y(),
-                                unitCellLatticePosition.z()
-                            )
-                            + latticePositions[p];
+                                vector
+                                (
+                                    unitCellLatticePosition.x(),
+                                    unitCellLatticePosition.y(),
+                                    unitCellLatticePosition.z()
+                                )
+                              + latticePositions[p];
 
                             point globalPosition =
-                            anchor + (R & (latticeCellShape & latticePosition));
+                                anchor
+                              + (R & (latticeCellShape & latticePosition));
 
-                            partOfLayerInBounds =
-                            mesh_.bounds().contains(globalPosition);
+                            partOfLayerInBounds = mesh_.bounds().contains
+                            (
+                                globalPosition
+                            );
 
                             label cell = mesh_.findCell(globalPosition);
 
@@ -779,20 +794,29 @@ void Foam::moleculeCloud::initialiseMolecules
                                         );
 
                                         const vector& latticePosition =
-                                        vector
-                                        (
-                                            unitCellLatticePosition.x(),
-                                            unitCellLatticePosition.y(),
-                                            unitCellLatticePosition.z()
-                                        )
-                                        + latticePositions[p];
+                                            vector
+                                            (
+                                                unitCellLatticePosition.x(),
+                                                unitCellLatticePosition.y(),
+                                                unitCellLatticePosition.z()
+                                            )
+                                          + latticePositions[p];
 
-                                        point globalPosition = anchor
-                                        + (R
-                                        &(latticeCellShape & latticePosition));
+                                        point globalPosition =
+                                            anchor
+                                          + (
+                                                R
+                                              & (
+                                                    latticeCellShape
+                                                  & latticePosition
+                                                )
+                                            );
 
                                         partOfLayerInBounds =
-                                        mesh_.bounds().contains(globalPosition);
+                                            mesh_.bounds().contains
+                                            (
+                                                globalPosition
+                                            );
 
                                         label cell = mesh_.findCell
                                         (
@@ -840,20 +864,29 @@ void Foam::moleculeCloud::initialiseMolecules
                                         );
 
                                         const vector& latticePosition =
-                                        vector
-                                        (
-                                            unitCellLatticePosition.x(),
-                                            unitCellLatticePosition.y(),
-                                            unitCellLatticePosition.z()
-                                        )
-                                        + latticePositions[p];
+                                            vector
+                                            (
+                                                unitCellLatticePosition.x(),
+                                                unitCellLatticePosition.y(),
+                                                unitCellLatticePosition.z()
+                                            )
+                                          + latticePositions[p];
 
-                                        point globalPosition = anchor
-                                        + (R
-                                        &(latticeCellShape & latticePosition));
+                                        point globalPosition =
+                                            anchor
+                                          + (
+                                                R
+                                              & (
+                                                   latticeCellShape
+                                                 & latticePosition
+                                                )
+                                            );
 
                                         partOfLayerInBounds =
-                                        mesh_.bounds().contains(globalPosition);
+                                            mesh_.bounds().contains
+                                            (
+                                                globalPosition
+                                            );
 
                                         label cell = mesh_.findCell
                                         (
@@ -875,12 +908,12 @@ void Foam::moleculeCloud::initialiseMolecules
                                     }
 
                                     unitCellLatticePosition =
-                                    labelVector
-                                    (
-                                        - unitCellLatticePosition.y(),
-                                        unitCellLatticePosition.x(),
-                                        unitCellLatticePosition.z()
-                                    );
+                                        labelVector
+                                        (
+                                          - unitCellLatticePosition.y(),
+                                            unitCellLatticePosition.x(),
+                                            unitCellLatticePosition.z()
+                                        );
                                 }
                             }
                         }
@@ -913,7 +946,7 @@ void Foam::moleculeCloud::initialiseMolecules
                     }
 
                     molsPlacedThisIteration =
-                    this->size() - sizeBeforeIteration;
+                        this->size() - sizeBeforeIteration;
 
                     totalZoneMols += molsPlacedThisIteration;
 
@@ -1041,7 +1074,7 @@ Foam::moleculeCloud::moleculeCloud
     mesh_(mesh),
     pot_(pot),
     cellOccupancy_(mesh_.nCells()),
-    il_(mesh_, pot_.pairPotentials().rCutMaxSqr(), false),
+    il_(mesh_, pot_.pairPotentials().rCutMax(), true),
     constPropList_(),
     rndGen_(clock::getTime())
 {
@@ -1071,7 +1104,7 @@ Foam::moleculeCloud::moleculeCloud
     Cloud<molecule>(mesh, "moleculeCloud", false),
     mesh_(mesh),
     pot_(pot),
-    il_(mesh_),
+    il_(mesh_, 0.0, false),
     constPropList_(),
     rndGen_(clock::getTime())
 {
@@ -1137,7 +1170,7 @@ void Foam::moleculeCloud::applyConstraintsAndThermostats
 )
 {
     scalar temperatureCorrectionFactor =
-        sqrt(targetTemperature/measuredTemperature);
+        sqrt(targetTemperature/max(VSMALL, measuredTemperature));
 
     Info<< "----------------------------------------" << nl
         << "Temperature equilibration" << nl
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
index 9ddc1302c5..b92f4e5e7b 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
@@ -40,7 +40,7 @@ SourceFiles
 #include "molecule.H"
 #include "IOdictionary.H"
 #include "potential.H"
-#include "interactionLists.H"
+#include "InteractionLists.H"
 #include "labelVector.H"
 #include "Random.H"
 #include "fileName.H"
@@ -69,7 +69,7 @@ private:
 
         List<DynamicList<molecule*> > cellOccupancy_;
 
-        interactionLists il_;
+        InteractionLists<molecule> il_;
 
         List<molecule::constantProperties> constPropList_;
 
@@ -89,26 +89,14 @@ private:
 
         inline void evaluatePair
         (
-            molecule* molI,
-            molecule* molJ
-        );
-
-        inline void evaluatePair
-        (
-            molecule* molReal,
-            referredMolecule* molRef
+            molecule& molI,
+            molecule& molJ
         );
 
         inline bool evaluatePotentialLimit
         (
-            molecule* molI,
-            molecule* molJ
-        ) const;
-
-        inline bool evaluatePotentialLimit
-        (
-            molecule* molReal,
-            referredMolecule* molRef
+            molecule& molI,
+            molecule& molJ
         ) const;
 
         void calculateTetherForce();
@@ -197,7 +185,7 @@ public:
 
             inline const List<DynamicList<molecule*> >& cellOccupancy() const;
 
-            inline const interactionLists& il() const;
+            inline const InteractionLists<molecule>& il() const;
 
             inline const List<molecule::constantProperties> constProps() const;
 
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
index 1eb02fad7f..dbeee755a5 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
@@ -31,17 +31,17 @@ using namespace Foam::constant;
 
 inline void Foam::moleculeCloud::evaluatePair
 (
-    molecule* molI,
-    molecule* molJ
+    molecule& molI,
+    molecule& molJ
 )
 {
     const pairPotentialList& pairPot = pot_.pairPotentials();
 
     const pairPotential& electrostatic = pairPot.electrostatic();
 
-    label idI = molI->id();
+    label idI = molI.id();
 
-    label idJ = molJ->id();
+    label idJ = molJ.id();
 
     const molecule::constantProperties& constPropI(constProps(idI));
 
@@ -70,7 +70,7 @@ inline void Foam::moleculeCloud::evaluatePair
             if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
             {
                 vector rsIsJ =
-                    molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
+                    molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
 
                 scalar rsIsJMagSq = magSqr(rsIsJ);
 
@@ -82,38 +82,34 @@ inline void Foam::moleculeCloud::evaluatePair
                         (rsIsJ/rsIsJMag)
                        *pairPot.force(idsI, idsJ, rsIsJMag);
 
-                    molI->siteForces()[sI] += fsIsJ;
+                    molI.siteForces()[sI] += fsIsJ;
 
-                    molJ->siteForces()[sJ] += -fsIsJ;
+                    molJ.siteForces()[sJ] += -fsIsJ;
 
                     scalar potentialEnergy
                     (
                         pairPot.energy(idsI, idsJ, rsIsJMag)
                     );
 
-                    molI->potentialEnergy() += 0.5*potentialEnergy;
+                    molI.potentialEnergy() += 0.5*potentialEnergy;
 
-                    molJ->potentialEnergy() += 0.5*potentialEnergy;
+                    molJ.potentialEnergy() += 0.5*potentialEnergy;
 
-                    vector rIJ = molI->position() - molJ->position();
+                    vector rIJ = molI.position() - molJ.position();
 
                     tensor virialContribution =
                         (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
 
-                    molI->rf() += virialContribution;
+                    molI.rf() += virialContribution;
 
-                    molJ->rf() += virialContribution;
-
-                    // molI->rf() += rsIsJ * fsIsJ;
-
-                    // molJ->rf() += rsIsJ * fsIsJ;
+                    molJ.rf() += virialContribution;
                 }
             }
 
             if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
             {
                 vector rsIsJ =
-                molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
+                molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
 
                 scalar rsIsJMagSq = magSqr(rsIsJ);
 
@@ -129,149 +125,26 @@ inline void Foam::moleculeCloud::evaluatePair
                         (rsIsJ/rsIsJMag)
                        *chargeI*chargeJ*electrostatic.force(rsIsJMag);
 
-                    molI->siteForces()[sI] += fsIsJ;
+                    molI.siteForces()[sI] += fsIsJ;
 
-                    molJ->siteForces()[sJ] += -fsIsJ;
+                    molJ.siteForces()[sJ] += -fsIsJ;
 
                     scalar potentialEnergy =
                         chargeI*chargeJ
                        *electrostatic.energy(rsIsJMag);
 
-                    molI->potentialEnergy() += 0.5*potentialEnergy;
+                    molI.potentialEnergy() += 0.5*potentialEnergy;
 
-                    molJ->potentialEnergy() += 0.5*potentialEnergy;
+                    molJ.potentialEnergy() += 0.5*potentialEnergy;
 
-                    vector rIJ = molI->position() - molJ->position();
+                    vector rIJ = molI.position() - molJ.position();
 
                     tensor virialContribution =
                         (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
 
-                    molI->rf() += virialContribution;
+                    molI.rf() += virialContribution;
 
-                    molJ->rf() += virialContribution;
-
-                    // molI->rf() += rsIsJ * fsIsJ;
-
-                    // molJ->rf() += rsIsJ * fsIsJ;
-                }
-            }
-        }
-    }
-}
-
-
-inline void Foam::moleculeCloud::evaluatePair
-(
-    molecule* molReal,
-    referredMolecule* molRef
-)
-{
-    const pairPotentialList& pairPot = pot_.pairPotentials();
-
-    const pairPotential& electrostatic = pairPot.electrostatic();
-
-    label idReal = molReal->id();
-
-    label idRef = molRef->id();
-
-    const molecule::constantProperties& constPropReal(constProps(idReal));
-
-    const molecule::constantProperties& constPropRef(constProps(idRef));
-
-    List<label> siteIdsReal = constPropReal.siteIds();
-
-    List<label> siteIdsRef = constPropRef.siteIds();
-
-    List<bool> pairPotentialSitesReal = constPropReal.pairPotentialSites();
-
-    List<bool> electrostaticSitesReal = constPropReal.electrostaticSites();
-
-    List<bool> pairPotentialSitesRef = constPropRef.pairPotentialSites();
-
-    List<bool> electrostaticSitesRef = constPropRef.electrostaticSites();
-
-    forAll(siteIdsReal, sReal)
-    {
-        label idsReal(siteIdsReal[sReal]);
-
-        forAll(siteIdsRef, sRef)
-        {
-            label idsRef(siteIdsRef[sRef]);
-
-            if (pairPotentialSitesReal[sReal] && pairPotentialSitesRef[sRef])
-            {
-                vector rsRealsRef =
-                    molReal->sitePositions()[sReal]
-                  - molRef->sitePositions()[sRef];
-
-                scalar rsRealsRefMagSq = magSqr(rsRealsRef);
-
-                if (pairPot.rCutSqr(idsReal, idsRef, rsRealsRefMagSq))
-                {
-                    scalar rsRealsRefMag = mag(rsRealsRef);
-
-                    vector fsRealsRef =
-                        (rsRealsRef/rsRealsRefMag)
-                       *pairPot.force(idsReal, idsRef, rsRealsRefMag);
-
-                    molReal->siteForces()[sReal] += fsRealsRef;
-
-                    scalar potentialEnergy
-                    (
-                        pairPot.energy(idsReal, idsRef, rsRealsRefMag)
-                    );
-
-                    molReal->potentialEnergy() += 0.5*potentialEnergy;
-
-                    vector rRealRef = molReal->position() - molRef->position();
-
-                    molReal->rf() +=
-                        (rsRealsRef*fsRealsRef)
-                       *(rsRealsRef & rRealRef)
-                       /rsRealsRefMagSq;
-
-                    // molReal->rf() += rsRealsRef * fsRealsRef;
-
-                }
-            }
-
-            if (electrostaticSitesReal[sReal] && electrostaticSitesRef[sRef])
-            {
-                vector rsRealsRef =
-                    molReal->sitePositions()[sReal]
-                  - molRef->sitePositions()[sRef];
-
-                scalar rsRealsRefMagSq = magSqr(rsRealsRef);
-
-                if (rsRealsRefMagSq <= electrostatic.rCutSqr())
-                {
-                    scalar rsRealsRefMag = mag(rsRealsRef);
-
-                    scalar chargeReal = constPropReal.siteCharges()[sReal];
-
-                    scalar chargeRef = constPropRef.siteCharges()[sRef];
-
-                    vector fsRealsRef =
-                        (rsRealsRef/rsRealsRefMag)
-                       *chargeReal*chargeRef
-                       *electrostatic.force(rsRealsRefMag);
-
-                    molReal->siteForces()[sReal] += fsRealsRef;
-
-                    scalar potentialEnergy =
-                        chargeReal*chargeRef
-                       *electrostatic.energy(rsRealsRefMag);
-
-                    molReal->potentialEnergy() += 0.5*potentialEnergy;
-
-                    vector rRealRef = molReal->position() - molRef->position();
-
-                    molReal->rf() +=
-                        (rsRealsRef*fsRealsRef)
-                       *(rsRealsRef & rRealRef)
-                       /rsRealsRefMagSq;
-
-                    // molReal->rf() += rsRealsRef * fsRealsRef;
+                    molJ.rf() += virialContribution;
                 }
             }
         }
@@ -281,17 +154,17 @@ inline void Foam::moleculeCloud::evaluatePair
 
 inline bool Foam::moleculeCloud::evaluatePotentialLimit
 (
-    molecule* molI,
-    molecule* molJ
+    molecule& molI,
+    molecule& molJ
 ) const
 {
     const pairPotentialList& pairPot = pot_.pairPotentials();
 
     const pairPotential& electrostatic = pairPot.electrostatic();
 
-    label idI = molI->id();
+    label idI = molI.id();
 
-    label idJ = molJ->id();
+    label idJ = molJ.id();
 
     const molecule::constantProperties& constPropI(constProps(idI));
 
@@ -320,7 +193,7 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
             if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
             {
                 vector rsIsJ =
-                    molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
+                    molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
 
                 scalar rsIsJMagSq = magSqr(rsIsJ);
 
@@ -369,7 +242,7 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
             if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
             {
                 vector rsIsJ =
-                    molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
+                    molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
 
                 scalar rsIsJMagSq = magSqr(rsIsJ);
 
@@ -422,157 +295,6 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
 }
 
 
-inline bool Foam::moleculeCloud::evaluatePotentialLimit
-(
-    molecule* molReal,
-    referredMolecule* molRef
-) const
-{
-    const pairPotentialList& pairPot = pot_.pairPotentials();
-
-    const pairPotential& electrostatic = pairPot.electrostatic();
-
-    label idReal = molReal->id();
-
-    label idRef = molRef->id();
-
-    const molecule::constantProperties& constPropReal(constProps(idReal));
-
-    const molecule::constantProperties& constPropRef(constProps(idRef));
-
-    List<label> siteIdsReal = constPropReal.siteIds();
-
-    List<label> siteIdsRef = constPropRef.siteIds();
-
-    List<bool> pairPotentialSitesReal = constPropReal.pairPotentialSites();
-
-    List<bool> electrostaticSitesReal = constPropReal.electrostaticSites();
-
-    List<bool> pairPotentialSitesRef = constPropRef.pairPotentialSites();
-
-    List<bool> electrostaticSitesRef = constPropRef.electrostaticSites();
-
-    forAll(siteIdsReal, sReal)
-    {
-        label idsReal(siteIdsReal[sReal]);
-
-        forAll(siteIdsRef, sRef)
-        {
-            label idsRef(siteIdsRef[sRef]);
-
-            if (pairPotentialSitesReal[sReal] && pairPotentialSitesRef[sRef])
-            {
-                vector rsRealsRef =
-                    molReal->sitePositions()[sReal]
-                  - molRef->sitePositions()[sRef];
-
-                scalar rsRealsRefMagSq = magSqr(rsRealsRef);
-
-                if (pairPot.rCutSqr(idsReal, idsRef, rsRealsRefMagSq))
-                {
-                    scalar rsRealsRefMag = mag(rsRealsRef);
-
-                    // Guard against pairPot.energy being evaluated
-                    // if rRealRefMag < SMALL. A floating point exception will
-                    // happen otherwise.
-
-                    if (rsRealsRefMag < SMALL)
-                    {
-                        WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
-                            << "Molecule site pair closer than "
-                            << SMALL
-                            << ": mag separation = " << rsRealsRefMag
-                            << ". These may have been placed on top of each"
-                            << " other by a rounding error in mdInitialise in"
-                            << " parallel or a block filled with molecules"
-                            << " twice. Removing one of the molecules."
-                            << endl;
-
-                        return true;
-                    }
-
-                    // Guard against pairPot.energy being evaluated if
-                    // rRealRefMag < rMin.  A tabulation lookup error will occur
-                    // otherwise.
-
-                    if (rsRealsRefMag < pairPot.rMin(idsReal, idsRef))
-                    {
-                        return true;
-                    }
-
-                    if
-                    (
-                        mag(pairPot.energy(idsReal, idsRef, rsRealsRefMag))
-                      > pot_.potentialEnergyLimit()
-                    )
-                    {
-                        return true;
-                    };
-                }
-            }
-
-            if (electrostaticSitesReal[sReal] && electrostaticSitesRef[sRef])
-            {
-                vector rsRealsRef =
-                    molReal->sitePositions()[sReal]
-                  - molRef->sitePositions()[sRef];
-
-                scalar rsRealsRefMagSq = magSqr(rsRealsRef);
-
-                if (pairPot.rCutMaxSqr(rsRealsRefMagSq))
-                {
-                    scalar rsRealsRefMag = mag(rsRealsRef);
-
-                    // Guard against pairPot.energy being evaluated
-                    // if rRealRefMag < SMALL. A floating point exception will
-                    // happen otherwise.
-
-                    if (rsRealsRefMag < SMALL)
-                    {
-                        WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
-                            << "Molecule site pair closer than "
-                            << SMALL
-                            << ": mag separation = " << rsRealsRefMag
-                            << ". These may have been placed on top of each"
-                            << " other by a rounding error in mdInitialise in"
-                            << " parallel or a block filled with molecules"
-                            << " twice. Removing one of the molecules."
-                            << endl;
-
-                        return true;
-                    }
-
-                    if (rsRealsRefMag < electrostatic.rMin())
-                    {
-                        return true;
-                    }
-
-                    scalar chargeReal = constPropReal.siteCharges()[sReal];
-
-                    scalar chargeRef = constPropRef.siteCharges()[sRef];
-
-                    if
-                    (
-                        mag
-                        (
-                            chargeReal
-                           *chargeRef
-                           *electrostatic.energy(rsRealsRefMag)
-                        )
-                      > pot_.potentialEnergyLimit()
-                    )
-                    {
-                        return true;
-                    };
-                }
-            }
-        }
-    }
-
-    return false;
-}
-
-
 inline Foam::vector Foam::moleculeCloud::equipartitionLinearVelocity
 (
     scalar temperature,
@@ -638,7 +360,7 @@ inline const Foam::List<Foam::DynamicList<Foam::molecule*> >&
 }
 
 
-inline const Foam::interactionLists&
+inline const Foam::InteractionLists<Foam::molecule>&
     Foam::moleculeCloud::il() const
 {
     return il_;
diff --git a/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCubeArgon/0/U b/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCubeArgon/0/U
new file mode 100644
index 0000000000..5e3b0435c4
--- /dev/null
+++ b/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCubeArgon/0/U
@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    periodicX
+    {
+        type            cyclic;
+    }
+    periodicY
+    {
+        type            cyclic;
+    }
+    periodicZ
+    {
+        type            cyclic;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCubeWater/0/U b/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCubeWater/0/U
new file mode 100644
index 0000000000..5e3b0435c4
--- /dev/null
+++ b/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCubeWater/0/U
@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    periodicX
+    {
+        type            cyclic;
+    }
+    periodicY
+    {
+        type            cyclic;
+    }
+    periodicZ
+    {
+        type            cyclic;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/discreteMethods/molecularDynamics/mdFoam/nanoNozzle/0/U b/tutorials/discreteMethods/molecularDynamics/mdFoam/nanoNozzle/0/U
new file mode 100644
index 0000000000..b9dc8a00f4
--- /dev/null
+++ b/tutorials/discreteMethods/molecularDynamics/mdFoam/nanoNozzle/0/U
@@ -0,0 +1,55 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    sectionAEnd
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    sectionCEnd
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    front
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    back
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    top
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    bottom
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //