diff --git a/applications/solvers/molecularDynamics/mdFoam/mdFoam.C b/applications/solvers/molecularDynamics/mdFoam/mdFoam.C
index 265514c3ac..37c5cbff9a 100644
--- a/applications/solvers/molecularDynamics/mdFoam/mdFoam.C
+++ b/applications/solvers/molecularDynamics/mdFoam/mdFoam.C
@@ -31,8 +31,8 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-//#include "md.h"
-#include "interactionLists.H"
+//#include "md.H"
+#include "potential.H"
 
 int main(int argc, char *argv[])
 {
@@ -41,9 +41,7 @@ int main(int argc, char *argv[])
 #   include "createTime.H"
 #   include "createMesh.H"
 
-    interactionLists il(mesh, (1e-9*1e-9), false);
-
-    //    Info<< labelListList(il.dil()) << endl;
+    potential pot(mesh);
 
     Info << "\nStarting time loop\n" << endl;
 
diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCellList/referredCellList.H b/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCellList/referredCellList.H
index f04f03530f..335ccc00f2 100644
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCellList/referredCellList.H
+++ b/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCellList/referredCellList.H
@@ -37,7 +37,7 @@ SourceFiles
 #define referredCellList_H
 
 #include "referredCell.H"
-//#include "molecule.H"
+#include "molecule.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/createMDFields.H b/src/lagrangian/molecularDynamics/molecule/mdTools/createMDFields.H
index dccd1891a7..730db2d0a8 100644
--- a/src/lagrangian/molecularDynamics/molecule/mdTools/createMDFields.H
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/createMDFields.H
@@ -2,25 +2,25 @@
 
 List< scalarField > allSpeciesN_RU
 (
-    molecules.pairPotentials().nIds(),
+    molecules.potential().nIds(),
     scalarField (mesh.nCells(), 0.0)
 );
 
 List< scalarField > allSpeciesM_RU
 (
-    molecules.pairPotentials().nIds(),
+    molecules.potential().nIds(),
     scalarField (mesh.nCells(), 0.0)
 );
 
 List< vectorField > allSpeciesVelocitySum_RU
 (
-    molecules.pairPotentials().nIds(),
+    molecules.potential().nIds(),
     vectorField (mesh.nCells(), vector::zero)
 );
 
 List< scalarField > allSpeciesVelocityMagSquaredSum_RU
 (
-    molecules.pairPotentials().nIds(),
+    molecules.potential().nIds(),
     scalarField (mesh.nCells(), 0.0)
 );
 
@@ -34,13 +34,13 @@ Info << nl << "Creating fields." << endl;
 
 PtrList<volScalarField> allSpeciesRhoN
 (
-    molecules.pairPotentials().nIds()
+    molecules.potential().nIds()
 );
 
 forAll (allSpeciesRhoN, rN)
 {
     Info << "    Creating number density field for "
-        << molecules.pairPotentials().idList()[rN] << endl;
+        << molecules.potential().idList()[rN] << endl;
 
     allSpeciesRhoN.set
     (
@@ -49,7 +49,7 @@ forAll (allSpeciesRhoN, rN)
         (
             IOobject
             (
-                "rhoN_" + molecules.pairPotentials().idList()[rN],
+                "rhoN_" + molecules.potential().idList()[rN],
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,
@@ -89,13 +89,13 @@ totalRhoN.correctBoundaryConditions();
 
 PtrList<volScalarField> allSpeciesRhoM
 (
-    molecules.pairPotentials().nIds()
+    molecules.potential().nIds()
 );
 
 forAll (allSpeciesRhoM, rM)
 {
     Info << "    Creating mass density field for "
-        << molecules.pairPotentials().idList()[rM] << endl;
+        << molecules.potential().idList()[rM] << endl;
 
     allSpeciesRhoM.set
     (
@@ -104,7 +104,7 @@ forAll (allSpeciesRhoM, rM)
         (
             IOobject
             (
-                "rhoM_" + molecules.pairPotentials().idList()[rM],
+                "rhoM_" + molecules.potential().idList()[rM],
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,
@@ -144,13 +144,13 @@ totalRhoM.correctBoundaryConditions();
 
 PtrList<volVectorField> allSpeciesVelocity
 (
-    molecules.pairPotentials().nIds()
+    molecules.potential().nIds()
 );
 
 forAll (allSpeciesVelocity, v)
 {
     Info << "    Creating velocity field for "
-        << molecules.pairPotentials().idList()[v] << endl;
+        << molecules.potential().idList()[v] << endl;
 
     allSpeciesVelocity.set
     (
@@ -159,7 +159,7 @@ forAll (allSpeciesVelocity, v)
         (
             IOobject
             (
-                "velocity_" + molecules.pairPotentials().idList()[v],
+                "velocity_" + molecules.potential().idList()[v],
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,
@@ -218,13 +218,13 @@ volVectorField totalVelocity
 
 PtrList<volScalarField> allSpeciesTemperature
 (
-    molecules.pairPotentials().nIds()
+    molecules.potential().nIds()
 );
 
 forAll (allSpeciesTemperature, t)
 {
     Info << "    Creating temperature field for "
-        << molecules.pairPotentials().idList()[t] << endl;
+        << molecules.potential().idList()[t] << endl;
 
     allSpeciesTemperature.set
     (
@@ -233,7 +233,7 @@ forAll (allSpeciesTemperature, t)
         (
             IOobject
             (
-                "temperature_" + molecules.pairPotentials().idList()[t],
+                "temperature_" + molecules.potential().idList()[t],
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,
@@ -274,13 +274,13 @@ totalTemperature.correctBoundaryConditions();
 
 PtrList<volScalarField> allSpeciesMeanKE
 (
-    molecules.pairPotentials().nIds()
+    molecules.potential().nIds()
 );
 
 forAll (allSpeciesMeanKE, mKE)
 {
     Info << "    Creating mean kinetic energy field for "
-        << molecules.pairPotentials().idList()[mKE] << endl;
+        << molecules.potential().idList()[mKE] << endl;
 
     allSpeciesMeanKE.set
     (
@@ -289,7 +289,7 @@ forAll (allSpeciesMeanKE, mKE)
         (
             IOobject
             (
-                "meanKE_" + molecules.pairPotentials().idList()[mKE],
+                "meanKE_" + molecules.potential().idList()[mKE],
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/resetMDFields.H b/src/lagrangian/molecularDynamics/molecule/mdTools/resetMDFields.H
index 337647ec02..ff274a3336 100644
--- a/src/lagrangian/molecularDynamics/molecule/mdTools/resetMDFields.H
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/resetMDFields.H
@@ -2,25 +2,25 @@ if (runTime.outputTime())
 {
     allSpeciesN_RU = List< scalarField >
     (
-        molecules.pairPotentials().nIds(),
+        molecules.potential().nIds(),
         scalarField (mesh.nCells(), 0.0)
     );
 
     allSpeciesM_RU = List< scalarField >
     (
-        molecules.pairPotentials().nIds(),
+        molecules.potential().nIds(),
         scalarField (mesh.nCells(), 0.0)
     );
 
     allSpeciesVelocitySum_RU = List< vectorField >
     (
-        molecules.pairPotentials().nIds(),
+        molecules.potential().nIds(),
         vectorField (mesh.nCells(), vector::zero)
     );
 
     allSpeciesVelocityMagSquaredSum_RU = List< scalarField >
     (
-        molecules.pairPotentials().nIds(),
+        molecules.potential().nIds(),
         scalarField (mesh.nCells(), 0.0)
     );
 }
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
index 1341695907..c2f3cda083 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
@@ -226,7 +226,7 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
 
     forAll(removalOrder_, rO)
     {
-        Info << " " << pairPotentials_.idList()[removalOrder_[rO]];
+        Info << " " << pot_.idList()[removalOrder_[rO]];
     }
 
     Info << nl ;
@@ -439,10 +439,15 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::moleculeCloud::moleculeCloud(const polyMesh& mesh)
+Foam::moleculeCloud::moleculeCloud
+(
+    const polyMesh& mesh,
+    const potential& pot
+)
 :
     Cloud<molecule>(mesh, "moleculeCloud", false),
     mesh_(mesh),
+    pot_(pot),
     referredInteractionList_(*this)
     cellOccupancy_(mesh_.nCells())
 {
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
index 14cee4ddac..7dd06e926e 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
@@ -134,8 +134,12 @@ public:
 
     // Constructors
 
-        //- Construct given mesh
-        moleculeCloud(const polyMesh& mesh);
+        //- Construct given mesh and potential references
+        moleculeCloud
+        (
+            const polyMesh& mesh,
+            const potential& pot
+        );
 
         //- Construct given polyMesh and fields of components.  Intended for use
         //- by the molConfig utility
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
index 6b163bb407..7406e6130d 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
@@ -32,6 +32,8 @@ inline void Foam::moleculeCloud::evaluatePair
     molecule* molJ
 )
 {
+    const pairPotentialList& pairPot(pot_.pairPotentials());
+
     idI = molI->id();
 
     idJ = molJ->id();
@@ -40,11 +42,11 @@ inline void Foam::moleculeCloud::evaluatePair
 
     rIJMagSq = magSqr(rIJ);
 
-    if(pairPotentials_.rCutSqr(idI, idJ, rIJMagSq))
+    if(pairPot.rCutSqr(idI, idJ, rIJMagSq))
     {
         rIJMag = mag(rIJ);
 
-        fIJ = (rIJ/rIJMag)*pairPotentials_.force(idI, idJ, rIJMag);
+        fIJ = (rIJ/rIJMag)*pairPot.force(idI, idJ, rIJMag);
 
         // Acceleration increment for molI
         molI->A() += fIJ/(molI->mass());
@@ -54,7 +56,7 @@ inline void Foam::moleculeCloud::evaluatePair
 
         scalar potentialEnergy
         (
-            pairPotentials_.energy(idI, idJ, rIJMag)
+            pairPot.energy(idI, idJ, rIJMag)
         );
 
         molI->potentialEnergy() += 0.5*potentialEnergy;
@@ -73,6 +75,8 @@ inline void Foam::moleculeCloud::evaluatePair
     referredMolecule* molRef
 )
 {
+    const pairPotentialList& pairPot(pot_.pairPotentials());
+
     idReal = molReal->id();
 
     idRef = molRef->id();
@@ -81,12 +85,12 @@ inline void Foam::moleculeCloud::evaluatePair
 
     rRealRefMagSq = magSqr(rRealRef);
 
-    if (pairPotentials_.rCutSqr(idReal, idRef, rRealRefMagSq))
+    if (pairPot.rCutSqr(idReal, idRef, rRealRefMagSq))
     {
         rRealRefMag = mag(rRealRef);
 
         fRealRef = (rRealRef/rRealRefMag)
-            *pairPotentials_.force(idReal, idRef, rRealRefMag);
+            *pairPot.force(idReal, idRef, rRealRefMag);
 
         // Acceleration increment for molReal
         molReal->A() += fRealRef/(molReal->mass());
@@ -94,7 +98,7 @@ inline void Foam::moleculeCloud::evaluatePair
         // Adding a contribution of 1/2 of the potential energy
         // from this interaction
         molReal->potentialEnergy() +=
-        0.5*pairPotentials_.energy(idReal, idRef, rRealRefMag);
+            0.5*pairPot.energy(idReal, idRef, rRealRefMag);
 
         molReal->rf() += rRealRef * fRealRef;
     }
@@ -112,13 +116,15 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
         << "(label idI, label idJ, scalar rIJMagSq) const "
         << "to pairPotential class."
 
+    const pairPotentialList& pairPot(pot_.pairPotentials());
+
     bool remove = false;
 
-    if (pairPotentials_.rCutSqr(idI, idJ, rIJMagSq))
+    if (pairPot.rCutSqr(idI, idJ, rIJMagSq))
     {
         rIJMag = mag(rIJ);
 
-        // Guard against pairPotentials_.energy being evaluated
+        // Guard against pairPot.energy being evaluated
         // if rIJMag < SMALL. A floating point exception will
         // happen otherwise.
 
@@ -141,11 +147,11 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
             remove = true;
         }
 
-        // Guard against pairPotentials_.energy being
+        // Guard against pairPot.energy being
         // evaluated if rIJMag < rMin. A tabulation lookup
         // error will occur otherwise.
 
-        if (rIJMag < pairPotentials_.rMin(idI, idJ))
+        if (rIJMag < pairPot.rMin(idI, idJ))
         {
             remove = true;
         }
@@ -154,7 +160,7 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
         {
             if
             (
-                pairPotentials_.energy(idI, idJ, rIJMag)
+                pairPot.energy(idI, idJ, rIJMag)
               > potentialEnergyLimit_
             )
             {
@@ -245,47 +251,4 @@ inline const molecule::constantProperties& constProps(label id) const
 }
 
 
-// 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/Make/files b/src/lagrangian/molecularDynamics/potential/Make/files
index 2e85ac11bc..294bac6a58 100644
--- a/src/lagrangian/molecularDynamics/potential/Make/files
+++ b/src/lagrangian/molecularDynamics/potential/Make/files
@@ -1,7 +1,3 @@
-potential = potential
-
-$(potential) = potential.C
-
 pairPotential = pairPotential
 
 $(pairPotential)/pairPotentialList/pairPotentialList.C
@@ -38,5 +34,8 @@ $(tetherPotential)/derived/harmonicSpring/harmonicSpring.C
 $(tetherPotential)/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.C
 $(tetherPotential)/derived/pitchForkRing/pitchForkRing.C
 
-LIB = $(FOAM_LIBBIN)/libpotential
+potential = potential
 
+$(potential)/potential.C
+
+LIB = $(FOAM_LIBBIN)/libpotential
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.C b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.C
index 2566941f24..58b64072d6 100644
--- a/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.C
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.C
@@ -27,15 +27,14 @@ License
 #include "pairPotential.H"
 #include "energyScalingFunction.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
+    defineTypeNameAndDebug(pairPotential, 0);
+    defineRunTimeSelectionTable(pairPotential, dictionary);
+}
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-defineTypeNameAndDebug(pairPotential, 0);
-defineRunTimeSelectionTable(pairPotential, dictionary);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -230,8 +229,4 @@ bool Foam::pairPotential::read(const dictionary& pairPotentialProperties)
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/potential/potential.C b/src/lagrangian/molecularDynamics/potential/potential/potential.C
index 28dcf45472..c358d6b376 100644
--- a/src/lagrangian/molecularDynamics/potential/potential/potential.C
+++ b/src/lagrangian/molecularDynamics/potential/potential/potential.C
@@ -26,12 +26,8 @@ License
 
 #include "potential.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-
 Foam::potential::potential(const polyMesh& mesh)
 :
     mesh_(mesh)
@@ -85,7 +81,7 @@ Foam::potential::potential(const polyMesh& mesh)
 
         forAll(removalOrder_, rO)
         {
-            removalOrder_[rO] = findIndex(pairPotentials_.idList(), remOrd[rO]);
+            removalOrder_[rO] = findIndex(idList_, remOrd[rO]);
         }
     }
 
@@ -101,7 +97,7 @@ Foam::potential::potential(const polyMesh& mesh)
 
     const dictionary& pairDict = potentialDict.subDict("pair");
 
-    pairPotentials_.buildPotentials(idList, pairDict, mesh_);
+    pairPotentials_.buildPotentials(idList_, pairDict, mesh_);
 
     // ****************************************************************************
     // Tether potentials
@@ -117,7 +113,7 @@ Foam::potential::potential(const polyMesh& mesh)
 
         const dictionary& tetherDict = potentialDict.subDict("tether");
 
-        tetherPotentials_.buildPotentials(idList, tetherDict, tetherIdList);
+        tetherPotentials_.buildPotentials(idList_, tetherDict, tetherIdList);
     }
 
     // ****************************************************************************
diff --git a/src/lagrangian/molecularDynamics/potential/potential/potential.H b/src/lagrangian/molecularDynamics/potential/potential/potential.H
index e36dfc2af1..b319a76ee2 100644
--- a/src/lagrangian/molecularDynamics/potential/potential/potential.H
+++ b/src/lagrangian/molecularDynamics/potential/potential/potential.H
@@ -38,6 +38,7 @@ SourceFiles
 
 #include "polyMesh.H"
 #include "IOdictionary.H"
+#include "Time.H"
 #include "pairPotentialList.H"
 #include "tetherPotentialList.H"
 
@@ -93,7 +94,13 @@ public:
 
         // Access
 
-            inline scalar potentialEnergyLimit() const;
+            inline label nIds() const;
+
+            inline const List<word>& idList() const;
+
+            inline scalar potentialEnergyLimit() const;;
+
+            inline label nPairPotentials() const;
 
             inline const labelList& removalOrder() const;
 
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.C b/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.C
index e7cf17e32c..588ebd786b 100644
--- a/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.C
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.C
@@ -26,15 +26,14 @@ License
 
 #include "tetherPotential.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
+    defineTypeNameAndDebug(tetherPotential, 0);
+    defineRunTimeSelectionTable(tetherPotential, dictionary);
+}
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-defineTypeNameAndDebug(tetherPotential, 0);
-defineRunTimeSelectionTable(tetherPotential, dictionary);
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -60,8 +59,4 @@ bool Foam::tetherPotential::read(const dictionary& tetherPotentialProperties)
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //