diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.H b/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.H
index ebffa912da..2140078e26 100644
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.H
+++ b/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.H
@@ -185,9 +185,11 @@ public:
 
             inline const referredCellList& referredInteractionList() const;
 
-            inline const List<sendingReferralList>& cellSendingReferralLists() const;
+            inline const List<sendingReferralList>&
+                cellSendingReferralLists() const;
 
-            inline const List<receivingReferralList>& cellReceivingReferralLists() const;
+            inline const List<receivingReferralList>&
+                cellReceivingReferralLists() const;
 
             inline label nInteractingProcs() const;
 };
diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionListsI.H b/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionListsI.H
index a3205d5aa1..e2077a52cd 100644
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionListsI.H
+++ b/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionListsI.H
@@ -40,7 +40,8 @@ const Foam::directInteractionList& Foam::interactionLists::dil() const
 }
 
 
-inline const Foam::referredCellList& Foam::interactionLists::referredInteractionList() const
+inline const Foam::referredCellList&
+    Foam::interactionLists::referredInteractionList() const
 {
     return ril_;
 }
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
index 9c2b34ef39..181e9184f3 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
@@ -85,10 +85,11 @@ bool Foam::molecule::move(molecule::trackData& td)
 
     if (td.part() == 0)
     {
-        // Leapfrog velocity adjust part, required before and after
-        // tracking+force part
+        // First leapfrog velocity adjust part, required before tracking+force part
 
-        U_ += 0.5*deltaT*A_;
+        v_ += 0.5*deltaT*a_;
+
+        omega1 += 0.5*deltaT*alpha_;
     }
     else if (td.part() == 1)
     {
@@ -99,7 +100,7 @@ bool Foam::molecule::move(molecule::trackData& td)
             // set the lagrangian time-step
             scalar dt = min(dtMax, tEnd);
 
-            dt *= trackToFace(position() + dt*U_, td);
+            dt *= trackToFace(position() + dt*v_, td);
 
             tEnd -= dt;
             stepFraction() = 1.0 - tEnd/deltaT;
@@ -115,6 +116,33 @@ bool Foam::molecule::move(molecule::trackData& td)
 
         setSitePositions(td.molCloud().constProps(id_));
     }
+    else if (td.part() == 3)
+    {
+        // Second leapfrog velocity adjust part, required after tracking+force part
+
+        const diagTensor& I;  // TAKE REFERENCE TO CONST PROP TO GET THIS
+
+        scalar m;  // TAKE REFERENCE TO CONST PROP TO GET THIS
+
+        a_ = vector::zero;
+
+        vector tau(vector::zero);
+
+        forAll(siteForces_, s)
+        {
+            const vector& f = siteForces_[s];
+
+            a_ += f/m;
+
+            tau += (sitePositions_[s] - position_) ^ f;
+        }
+
+        alpha_ = R_ & inv(I) & R_.T() & tau;
+
+        v_ += 0.5*deltaT*a_;
+
+        omega_ += 0.5*deltaT*alpha_;
+    }
     else
     {
         FatalErrorIn("molecule::move(molecule::trackData& td)") << nl
@@ -173,8 +201,8 @@ void Foam::molecule::hitWallPatch
     vector nw = wpp.faceAreas()[wpp.whichFace(face())];
     nw /= mag(nw);
 
-    scalar Un = U_ & nw;
-//     vector Ut = U_ - Un*nw;
+    scalar vn = v_ & nw;
+//     vector vt = v_ - vn*nw;
 
 //     Random rand(clock::getTime());
 
@@ -186,11 +214,11 @@ void Foam::molecule::hitWallPatch
 //     {
 //         // Diffuse reflection
 //
-//         vector tw1 = Ut/mag(Ut);
+//         vector tw1 = vt/mag(vt);
 //
 //         vector tw2 = nw ^ tw1;
 //
-//         U_ = sqrt(wallTemp/mass_)*rand.GaussNormal()*tw1
+//         V_ = sqrt(wallTemp/mass_)*rand.GaussNormal()*tw1
 //                 + sqrt(wallTemp/mass_)*rand.GaussNormal()*tw2
 //                 - mag(sqrt(wallTemp/mass_)*rand.GaussNormal())*nw;
 //     }
@@ -199,9 +227,9 @@ void Foam::molecule::hitWallPatch
 //     {
         // Specular reflection
 
-        if (Un > 0)
+        if (vn > 0)
         {
-            U_ -= 2*Un*nw;
+            v_ -= 2*vn*nw;
         }
 
 //     }
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
index d1d740206b..1341695907 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
@@ -42,6 +42,401 @@ Foam::scalar Foam::moleculeCloud::elementaryCharge = 1.602176487e-19;
 Foam::scalar Foam::moleculeCloud::vacuumPermittivity = 8.854187817e-12;
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::moleculeCloud::buildCellOccupancy()
+{
+    forAll(cellOccupancy_, cO)
+    {
+        cellOccupancy_[cO].clear();
+    }
+
+    iterator mol(this->begin());
+
+    for
+    (
+        mol = this->begin();
+        mol != this->end();
+        ++mol
+    )
+    {
+        cellOccupancy_[mol().cell()].append(&mol());
+    }
+
+    forAll(cellOccupancy_, cO)
+    {
+        cellOccupancy_[cO].shrink();
+    }
+
+    il_ril().referMolecules(cellOccupancy_);
+}
+
+
+void Foam::moleculeCloud::calculatePairForce()
+{
+    {
+        iterator mol(this->begin());
+
+        // Real-Real interactions
+        vector rIJ;
+
+        scalar rIJMag;
+
+        scalar rIJMagSq;
+
+        vector fIJ;
+
+        label idI;
+
+        label idJ;
+
+        mol = this->begin();
+
+        molecule* molI = &mol();
+
+        molecule* molJ = &mol();
+
+        forAll(directInteractionList_, dIL)
+        {
+            forAll(cellOccupancy_[dIL],cellIMols)
+            {
+                molI = cellOccupancy_[dIL][cellIMols];
+
+                forAll(directInteractionList_[dIL], interactingCells)
+                {
+                    List< molecule* > cellJ =
+                    cellOccupancy_[directInteractionList_[dIL][interactingCells]];
+
+                    forAll(cellJ, cellJMols)
+                    {
+                        molJ = cellJ[cellJMols];
+
+                        evaluatelPair(molI, molJ);
+                    }
+                }
+
+                forAll(cellOccupancy_[dIL],cellIOtherMols)
+                {
+                    molJ = cellOccupancy_[dIL][cellIOtherMols];
+
+                    if (molJ > molI)
+                    {
+                        evaluatePair(molI, molJ);
+                    }
+                }
+            }
+        }
+    }
+
+    {
+        // Real-Referred interactions
+
+        vector rKL;
+
+        scalar rKLMag;
+
+        scalar rKLMagSq;
+
+        vector fKL;
+
+        label idK;
+
+        label idL;
+
+        molecule* molK = &mol();
+
+        forAll(referredInteractionList_, rIL)
+        {
+            const List<label>& realCells =
+            referredInteractionList_[rIL].realCellsForInteraction();
+
+            forAll(referredInteractionList_[rIL], refMols)
+            {
+                referredMolecule* molL = &(referredInteractionList_[rIL][refMols]);
+
+                forAll(realCells, rC)
+                {
+                    List<molecule*> cellK = cellOccupancy_[realCells[rC]];
+
+                    forAll(cellK, cellKMols)
+                    {
+                        molK = cellK[cellKMols];
+
+                        evaluatePair(molK, molL);
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::moleculeCloud::calculateTetherForce()
+{
+    iterator mol(this->begin());
+
+    vector rIT;
+
+    vector fIT;
+
+    for (mol = this->begin(); mol != this->end(); ++mol)
+    {
+        if (mol().tethered())
+        {
+            rIT = mol().position() - mol().tetherPosition();
+
+            fIT = tetherPotentials_.force
+            (
+                mol().id(),
+                rIT
+            );
+
+            mol().A() += fIT/(mol().mass());
+
+            mol().potentialEnergy() += tetherPotentials_.energy
+            (
+                mol().id(),
+                rIT
+            );
+
+            mol().rf() += rIT*fIT;
+        }
+    }
+}
+
+
+void Foam::moleculeCloud::calculateExternalForce()
+{
+    iterator mol(this->begin());
+
+    // Info<< "Warning! Includes dissipation term!" << endl;
+
+    for (mol = this->begin(); mol != this->end(); ++mol)
+    {
+        mol().A() += gravity_;
+
+        // mol().A() += -1.0 * mol().U() /mol().mass();
+    }
+}
+
+
+void Foam::moleculeCloud::removeHighEnergyOverlaps()
+{
+    Info << nl << "Removing high energy overlaps, removal order:";
+
+    forAll(removalOrder_, rO)
+    {
+        Info << " " << pairPotentials_.idList()[removalOrder_[rO]];
+    }
+
+    Info << nl ;
+
+    label initialSize = this->size();
+
+    buildCellOccupancy();
+
+    // Real-Real interaction
+    iterator mol(this->begin());
+
+    {
+        vector rIJ;
+
+        scalar rIJMag;
+
+        scalar rIJMagSq;
+
+        label idI;
+
+        label idJ;
+
+        mol = this->begin();
+
+        molecule* molI = &mol();
+
+        molecule* molJ = &mol();
+
+        DynamicList<molecule*> molsToDelete;
+
+        forAll(directInteractionList_, dIL)
+        {
+            forAll(cellOccupancy_[dIL],cellIMols)
+            {
+                molI = cellOccupancy_[dIL][cellIMols];
+
+                forAll(directInteractionList_[dIL], interactingCells)
+                {
+                    List< molecule* > cellJ =
+                    cellOccupancy_[directInteractionList_[dIL][interactingCells]];
+
+                    forAll(cellJ, cellJMols)
+                    {
+                        molJ = cellJ[cellJMols];
+
+                        if (evaluatePotentialLimit(molI, molJ))
+                        {
+                            if
+                            (
+                                idI == idJ
+                             || findIndex(removalOrder_, idJ)
+                                    < findIndex(removalOrder_, idI)
+                            )
+                            {
+                                if (findIndex(molsToDelete, molJ) == -1)
+                                {
+                                    molsToDelete.append(molJ);
+                                }
+                            }
+                            else if (findIndex(molsToDelete, molI) == -1)
+                            {
+                                molsToDelete.append(molI);
+                            }
+                        }
+                    }
+                }
+            }
+
+            forAll(cellOccupancy_[dIL],cellIOtherMols)
+            {
+                molJ = cellOccupancy_[dIL][cellIOtherMols];
+
+                if (molJ > molI)
+                {
+                    if (evaluatePotentialLimit(molI, molJ))
+                    {
+                        if
+                        (
+                            idI == idJ
+                         || findIndex(removalOrder_, idJ)
+                                < findIndex(removalOrder_, idI)
+                        )
+                        {
+                            if (findIndex(molsToDelete, molJ) == -1)
+                            {
+                                molsToDelete.append(molJ);
+                            }
+                        }
+                        else if (findIndex(molsToDelete, molI) == -1)
+                        {
+                            molsToDelete.append(molI);
+                        }
+                    }
+                }
+            }
+        }
+
+        forAll (molsToDelete, mTD)
+        {
+            deleteParticle(*(molsToDelete[mTD]));
+        }
+    }
+
+    buildCellOccupancy();
+
+    // Real-Referred interaction
+
+    {
+        vector rKL;
+
+        scalar rKLMag;
+
+        scalar rKLMagSq;
+
+        label idK;
+
+        label idL;
+
+        molecule* molK = &mol();
+
+        DynamicList<molecule*> molsToDelete;
+
+        forAll(referredInteractionList_, rIL)
+        {
+            referredCell& refCell = referredInteractionList_[rIL];
+
+            forAll(refCell, refMols)
+            {
+                referredMolecule* molL = &(refCell[refMols]);
+
+                List <label> realCells = refCell.realCellsForInteraction();
+
+                forAll(realCells, rC)
+                {
+                    label cellK = realCells[rC];
+
+                    List<molecule*> cellKMols = cellOccupancy_[cellK];
+
+                    forAll(cellKMols, cKM)
+                    {
+                        molK = cellKMols[cKM];
+
+                        if (evaluatePotentialLimit(molK, molL))
+                        {
+                            if
+                            (
+                                findIndex(removalOrder_, idK)
+                                < findIndex(removalOrder_, idL)
+                            )
+                            {
+                                if (findIndex(molsToDelete, molK) == -1)
+                                {
+                                    molsToDelete.append(molK);
+                                }
+                            }
+
+                            else if
+                            (
+                                findIndex(removalOrder_, idK)
+                                == findIndex(removalOrder_, idL)
+                            )
+                            {
+                                if
+                                (
+                                    Pstream::myProcNo() == refCell.sourceProc()
+                                    && cellK <= refCell.sourceCell()
+                                )
+                                {
+                                    if (findIndex(molsToDelete, molK) == -1)
+                                    {
+                                        molsToDelete.append(molK);
+                                    }
+                                }
+
+                                else if
+                                (
+                                    Pstream::myProcNo() < refCell.sourceProc()
+                                )
+                                {
+                                    if (findIndex(molsToDelete, molK) == -1)
+                                    {
+                                        molsToDelete.append(molK);
+                                    }
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+
+        forAll (molsToDelete, mTD)
+        {
+            deleteParticle(*(molsToDelete[mTD]));
+        }
+    }
+
+    buildCellOccupancy();
+
+    label molsRemoved = initialSize - this->size();
+
+    if (Pstream::parRun())
+    {
+        reduce(molsRemoved, sumOp<label>());
+    }
+
+    Info << tab << molsRemoved << " molecules removed" << endl;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::moleculeCloud::moleculeCloud(const polyMesh& mesh)
@@ -115,48 +510,74 @@ Foam::moleculeCloud::moleculeCloud(const polyMesh& mesh)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::moleculeCloud::buildCellOccupancy()
-{
-    forAll(cellOccupancy_, cO)
-    {
-        cellOccupancy_[cO].clear();
-    }
-
-    iterator mol(this->begin());
-
-    for
-    (
-        mol = this->begin();
-        mol != this->end();
-        ++mol
-    )
-    {
-        cellOccupancy_[mol().cell()].append(&mol());
-    }
-
-    forAll(cellOccupancy_, cO)
-    {
-        cellOccupancy_[cO].shrink();
-    }
-
-    il_ril().referMolecules(cellOccupancy_);
-}
-
 void Foam::moleculeCloud::evolve()
 {
     molecule::trackData td0(*this, 0);
-    molecule::trackData td1(*this, 1);
-    molecule::trackData td2(*this, 2);
-
     Cloud<molecule>::move(td0);
 
+    molecule::trackData td1(*this, 1);
+    Cloud<molecule>::move(td1);
+
+    molecule::trackData td2(*this, 2);
+    Cloud<molecule>::move(td2);
+
     calculateForce();
 
-    Cloud<molecule>::move(td1);
+    molecule::trackData td2(*this, 3);
+    Cloud<molecule>::move(td3);
+}
 
-    Cloud<molecule>::move(td2);
 
-    Cloud<molecule>::move(td0);
+void Foam::moleculeCloud::calculateForce()
+{
+    buildCellOccupancy();
+
+    iterator mol(this->begin());
+
+    // Set accumulated quantities to zero
+    for (mol = this->begin(); mol != this->end(); ++mol)
+    {
+        mol().siteForces() = vector::zero;
+
+        mol().potentialEnergy() = 0.0;
+
+        mol().rf() = tensor::zero;
+    }
+
+    calculatePairForce();
+
+    calculateTetherForce();
+
+    calculateExternalForce();
+}
+
+
+void Foam::moleculeCloud::applyConstraintsAndThermostats
+(
+    const scalar targetTemperature,
+    const scalar measuredTemperature
+)
+{
+    scalar temperatureCorrectionFactor =
+        sqrt(targetTemperature/measuredTemperature);
+
+    Info<< "----------------------------------------" << nl
+        << "Temperature equilibration" << nl
+        << "Target temperature = "
+        << targetTemperature << nl
+        << "Measured temperature = "
+        << measuredTemperature << nl
+        << "Temperature correction factor = "
+        << temperatureCorrectionFactor << nl
+        << "----------------------------------------"
+        << endl;
+
+    iterator mol(this->begin());
+
+    for (mol = this->begin(); mol != this->end(); ++mol)
+    {
+        mol().u() *= temperatureCorrectionFactor;
+    }
 }
 
 
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
index 3fc7950fef..14d3d57664 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
@@ -63,13 +63,13 @@ private:
 
         const polyMesh& mesh_;
 
-            // MD solution parameters
+        // MD solution parameters
 
-            List<DynamicList<molecule*> > cellOccupancy_;
+        List<DynamicList<molecule*> > cellOccupancy_;
 
-            interactionLists il_;
+        interactionLists il_;
 
-            List<molecule::constantProperties> constPropList_;;
+        List<molecule::constantProperties> constPropList_;
 
     // Move to potential class and store reference to it ->
             scalar potentialEnergyLimit_;
@@ -83,29 +83,55 @@ private:
             vector gravity_;
     // <- Move to potential class and store reference to it
 
-        // Private Member Functions
+    // Private Member Functions
 
-            void applyConstraintsAndThermostats
-            (
-                const scalar targetTemperature,
-                const scalar measuredTemperature
-            );
+        //- Determine which molecules are in which cells
+        void buildCellOccupancy();
 
-            void calculateForce();
+        void calculatePairForce();
 
-            void calculatePairForce();
+        inline void evaluatePair
+        (
+            molecule* molI,
+            molecule* molJ
+        );
 
-            void calculateTetherForce();
+        inline void evaluatePair
+        (
+            molecule* molReal,
+            referredMolecule* molRef
+        )
 
-            void calculateExternalForce();
+        inline bool evaluatePotentialLimit
+        (
+            label idI,
+            label idJ,
+            scalar rIJMagSq
+        ) const;
 
-            void removeHighEnergyOverlaps();
+        inline bool evaluatePotentialLimit
+        (
+            molecule* molI,
+            molecule* molJ
+        ) const;
 
-            //- Disallow default bitwise copy construct
-            moleculeCloud(const moleculeCloud&);
+        inline bool evaluatePotentialLimit
+        (
+            molecule* molReal,
+            referredMolecule* molRef
+        ) const;
 
-            //- Disallow default bitwise assignment
-            void operator=(const moleculeCloud&);
+        void calculateTetherForce();
+
+        void calculateExternalForce();
+
+        void removeHighEnergyOverlaps();
+
+        //- Disallow default bitwise copy construct
+        moleculeCloud(const moleculeCloud&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const moleculeCloud&);
 
 
 public:
@@ -142,22 +168,31 @@ public:
 
     // Member Functions
 
-    //- Determine which molecules are in which cells
-        void buildCellOccupancy();
 
         //- Evolve the molecules (move, calculate forces, control state etc)
         void evolve();
 
-        //- Integrate the equations of motion using algorithm selected at
-        //  runtime from a dictionary. This will also call the function
-        //  to calculate the intermolecular forces (calculatePairForce()).
-
+        void calculateForce();
 
+        void applyConstraintsAndThermostats
+        (
+            const scalar targetTemperature,
+            const scalar measuredTemperature
+        );
 
         // Access
 
             inline const polyMesh& mesh() const;
 
+            inline const interactionLists& interactionLists() const;
+
+            inline const List<DynamicList<molecule*> >& cellOccupancy() const;
+
+            inline const List<molecule::constantProperties> constProps() const;
+
+            inline const molecule::constantProperties& constProps(label id) const;
+
+    // Move to potential class and store reference to it ->
             inline scalar potentialEnergyLimit() const;
 
             inline const labelList& removalOrder() const;
@@ -167,10 +202,7 @@ public:
             inline const tetherPotentialList& tetherPotentials() const;
 
             inline const vector& gravity() const;
-
-            inline const List<DynamicList<molecule*> >& cellOccupancy() const;
-
-
+    // <- Move to potential class and store reference to it
 
     // Member Operators
 
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
new file mode 100644
index 0000000000..7f45b64a60
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
@@ -0,0 +1,285 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+inline void Foam::moleculeCloud::evaluatePair
+(
+    molecule* molI,
+    molecule* molJ
+)
+{
+    idI = molI->id();
+
+    idJ = molJ->id();
+
+    rIJ = molI->position() - molJ->position();
+
+    rIJMagSq = magSqr(rIJ);
+
+    if(pairPotentials_.rCutSqr(idI, idJ, rIJMagSq))
+    {
+        rIJMag = mag(rIJ);
+
+        fIJ = (rIJ/rIJMag)*pairPotentials_.force(idI, idJ, rIJMag);
+
+        // Acceleration increment for molI
+        molI->A() += fIJ/(molI->mass());
+
+        // Acceleration increment for molJ
+        molJ->A() += -fIJ/(molJ->mass());
+
+        scalar potentialEnergy
+        (
+            pairPotentials_.energy(idI, idJ, rIJMag)
+        );
+
+        molI->potentialEnergy() += 0.5*potentialEnergy;
+
+        molJ->potentialEnergy() += 0.5*potentialEnergy;
+
+        molI->rf() += rIJ * fIJ;
+
+        molJ->rf() += rIJ * fIJ;
+    }
+}
+
+inline void Foam::moleculeCloud::evaluatePair
+(
+    molecule* molReal,
+    referredMolecule* molRef
+)
+{
+    idReal = molReal->id();
+
+    idRef = molRef->id();
+
+    rRealRef = molReal->position() - molRef->position();
+
+    rRealRefMagSq = magSqr(rRealRef);
+
+    if (pairPotentials_.rCutSqr(idReal, idRef, rRealRefMagSq))
+    {
+        rRealRefMag = mag(rRealRef);
+
+        fRealRef = (rRealRef/rRealRefMag)
+            *pairPotentials_.force(idReal, idRef, rRealRefMag);
+
+        // Acceleration increment for molReal
+        molReal->A() += fRealRef/(molReal->mass());
+
+        // Adding a contribution of 1/2 of the potential energy
+        // from this interaction
+        molReal->potentialEnergy() +=
+        0.5*pairPotentials_.energy(idReal, idRef, rRealRefMag);
+
+        molReal->rf() += rRealRef * fRealRef;
+    }
+}
+
+
+inline bool Foam::moleculeCloud::evaluatePotentialLimit
+(
+    label idI,
+    label idJ,
+    scalar rIJMagSq
+) const
+{
+    Info << "Move Foam::moleculeCloud::evaluatePotentialLimit "
+        << "(label idI, label idJ, scalar rIJMagSq) const "
+        << "to pairPotential class."
+
+    bool remove = false;
+
+    if (pairPotentials_.rCutSqr(idI, idJ, rIJMagSq))
+    {
+        rIJMag = mag(rIJ);
+
+        // Guard against pairPotentials_.energy being evaluated
+        // if rIJMag < SMALL. A floating point exception will
+        // happen otherwise.
+
+        if (rIJMag < SMALL)
+        {
+            WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
+            << "Real molecule pair "
+                << " idI = " << idI
+                << " at position " << molI->position()
+                << " idJ = " << idJ
+                << " at position " << molJ->position()
+                << " are closer than " << SMALL
+                << ": mag separation = " << rIJMag
+                << ". These may have been placed on top of each"
+                << " other by a rounding error in molConfig in parallel"
+                << " or a block filled with molecules twice."
+                << " Removing one of the molecules."
+                << endl;
+
+            remove = true;
+        }
+
+        // Guard against pairPotentials_.energy being
+        // evaluated if rIJMag < rMin. A tabulation lookup
+        // error will occur otherwise.
+
+        if (rIJMag < pairPotentials_.rMin(idI, idJ))
+        {
+            remove = true;
+        }
+
+        if (!remove)
+        {
+            if
+            (
+                pairPotentials_.energy(idI, idJ, rIJMag)
+              > potentialEnergyLimit_
+            )
+            {
+
+                remove = true;
+            }
+        }
+    }
+
+    return remove;
+}
+
+
+inline bool Foam::moleculeCloud::evaluatePotentialLimit
+(
+    molecule* molI,
+    molecule* molJ
+) const
+{
+    idI = molI->id();
+
+    idJ = molJ->id();
+
+    rIJ = molI->position() - molJ->position();
+
+    rIJMagSq = magSqr(rIJ);
+
+    return evaluatePotentialLimit(idI, idJ, rIJMagSq);
+}
+
+
+inline bool Foam::moleculeCloud::evaluatePotentialLimit
+(
+    molecule* molReal,
+    referredMolecule* molRef
+) const
+{
+
+    idReal = molReal->id();
+
+    idRef = molRef->id();
+
+    rRealRef = molReal->position() - molRef->position();
+
+    rRealRefMagSq = magSqr(rRealRef);
+
+    return evaluatePotentialRefimit(idReal, idRef, rRealRefMagSq);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::polyMesh& moleculeCloud::mesh() const
+{
+    return mesh_;
+}
+
+
+inline const Foam::interactionLists&
+    Foam::moleculeCloud::interactionLists() const
+{
+    return il_;
+}
+
+
+inline const Foam::List<DynamicList<molecule*> >&
+    Foam::moleculeCloud::cellOccupancy() const
+{
+    return cellOccupancy_;
+}
+
+
+inline const Foam::List<molecule::constantProperties> constProps() const
+{
+    return constPropList_;
+}
+
+
+inline const molecule::constantProperties& constProps(label id) const
+{
+    return constPropList_[id];
+}
+
+
+// Move to potential class and store reference to it ->
+
+inline Foam::scalar Foam::moleculeCloud::potentialEnergyLimit() const
+{
+    return potentialEnergyLimit_;
+}
+
+
+inline Foam::label Foam::moleculeCloud::nPairPotentials() const
+{
+    return pairPotentials_.size();
+}
+
+
+inline const Foam::labelList& Foam::moleculeCloud::removalOrder() const
+{
+    return removalOrder_;
+}
+
+
+inline const Foam::pairPotentialList&
+    Foam::moleculeCloud::pairPotentials() const
+{
+    return pairPotentials_;
+}
+
+
+inline const Foam::tetherPotentialList&
+    Foam::moleculeCloud::tetherPotentials() const
+{
+    return tetherPotentials_;
+}
+
+
+inline const Foam::vector& Foam::moleculeCloud::gravity() const
+{
+    return gravity_;
+}
+// <- Move to potential class and store reference to it
+
+
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.H b/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.H
index 4ae264d373..c56808e267 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.H
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.H
@@ -57,7 +57,7 @@ class pairPotentialList
     // Private data
 
         List<word> idList_;
-        
+
         scalar rCutMax_;
 
         scalar rCutMaxSqr_;
@@ -101,7 +101,7 @@ public:
         ~pairPotentialList();
 
     // Member Functions
-        
+
         void buildPotentials
         (
             const dictionary& idListDict,