diff --git a/applications/solvers/molecularDynamics/gnemdFoam/Make/options b/applications/solvers/molecularDynamics/gnemdFoam/Make/options
index 80f977345684ed70cfcff7d09a17b0dc1a7770d0..89431d062560428f33fd5f5986f6b937270edee1 100755
--- a/applications/solvers/molecularDynamics/gnemdFoam/Make/options
+++ b/applications/solvers/molecularDynamics/gnemdFoam/Make/options
@@ -1,5 +1,6 @@
 EXE_INC = \
-    -I$(LIB_SRC)/lagrangian/molecule/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
@@ -8,4 +9,5 @@ EXE_LIBS = \
     -lmeshTools \
     -lfiniteVolume \
     -llagrangian \
-    -lmolecule
+    -lmolecule \
+    -lpotential
diff --git a/applications/solvers/molecularDynamics/mdEquilibrationFoam/Make/options b/applications/solvers/molecularDynamics/mdEquilibrationFoam/Make/options
index 80f977345684ed70cfcff7d09a17b0dc1a7770d0..89431d062560428f33fd5f5986f6b937270edee1 100755
--- a/applications/solvers/molecularDynamics/mdEquilibrationFoam/Make/options
+++ b/applications/solvers/molecularDynamics/mdEquilibrationFoam/Make/options
@@ -1,5 +1,6 @@
 EXE_INC = \
-    -I$(LIB_SRC)/lagrangian/molecule/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
@@ -8,4 +9,5 @@ EXE_LIBS = \
     -lmeshTools \
     -lfiniteVolume \
     -llagrangian \
-    -lmolecule
+    -lmolecule \
+    -lpotential
diff --git a/applications/utilities/preProcessing/molConfig/Make/options b/applications/utilities/preProcessing/molConfig/Make/options
index 196a4d66ddba259f0d7547f781a9c925db63b0eb..aab9a2ca4f7e9119ebfa389b5b14b74911317f1a 100755
--- a/applications/utilities/preProcessing/molConfig/Make/options
+++ b/applications/utilities/preProcessing/molConfig/Make/options
@@ -3,7 +3,8 @@ EXE_INC = \
     -I$(velocityDistributions) \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/lagrangian/molecule/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
@@ -12,4 +13,5 @@ EXE_LIBS = \
     -ldynamicMesh \
     -lfiniteVolume \
     -llagrangian \
-    -lmolecule
+    -lmolecule \
+    -lpotential
diff --git a/src/lagrangian/Allwmake b/src/lagrangian/Allwmake
index 567f7a10defa3445265239cb9018bddfeed28790..088a4a3b5db2884a8cb70fbe7e8912b0e29d35cb 100755
--- a/src/lagrangian/Allwmake
+++ b/src/lagrangian/Allwmake
@@ -5,4 +5,7 @@ wmake libso basic
 wmake libso solidParticle
 wmake libso intermediate
 wmake libso dieselSpray
-wmake libso molecule
+wmake libso molecularDynamics/potential
+wmake libso molecularDynamics/molecule
+
+ 
diff --git a/src/lagrangian/molecule/Make/files b/src/lagrangian/molecularDynamics/molecule/Make/files
similarity index 84%
rename from src/lagrangian/molecule/Make/files
rename to src/lagrangian/molecularDynamics/molecule/Make/files
index d8063d097437f7ca2dfa6b8e20a74cafa47723c9..a490edafe1ae49ff8c38791440fe3300e214078a 100755
--- a/src/lagrangian/molecule/Make/files
+++ b/src/lagrangian/molecularDynamics/molecule/Make/files
@@ -12,9 +12,7 @@ referredCellList = referredCellList
 referredCell = referredCell
 referralLists = referralLists
 
-potentials = potentials
-pairPotential = $(potentials)/pairPotential
-tetherPotential = $(potentials)/tetherPotential
+tetherPotential = tetherPotential
 
 $(distribution)/distribution.C
 
@@ -42,11 +40,6 @@ $(moleculeCloud)/moleculeCloudIntegrateEquationsOfMotion.C
 $(moleculeCloud)/moleculeCloudRemoveHighEnergyOverlaps.C
 $(moleculeCloud)/moleculeCloudApplyConstraintsAndThermostats.C
 
-$(pairPotential)/basic/pairPotential.C
-$(pairPotential)/basic/pairPotentialList.C
-$(tetherPotential)/tetherPotential.C
-$(tetherPotential)/tetherPotentialList.C
-
 $(referralLists)/receivingReferralList.C
 $(referralLists)/sendingReferralList.C
 $(referredCellList)/referredCellList.C
diff --git a/src/lagrangian/molecule/Make/options b/src/lagrangian/molecularDynamics/molecule/Make/options
similarity index 57%
rename from src/lagrangian/molecule/Make/options
rename to src/lagrangian/molecularDynamics/molecule/Make/options
index a0fb316eecd4a5e63b9f261787d2ecf29ff957b9..cc91fd4bbef065a768be4f1abd66bf51346d8ef4 100755
--- a/src/lagrangian/molecule/Make/options
+++ b/src/lagrangian/molecularDynamics/molecule/Make/options
@@ -1,8 +1,10 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
-    -llagrangian
+    -llagrangian \
+    -lpotential
 
diff --git a/src/lagrangian/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulator.C b/src/lagrangian/molecularDynamics/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulator.C
similarity index 100%
rename from src/lagrangian/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulator.C
rename to src/lagrangian/molecularDynamics/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulator.C
diff --git a/src/lagrangian/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulator.H b/src/lagrangian/molecularDynamics/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulator.H
similarity index 98%
rename from src/lagrangian/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulator.H
rename to src/lagrangian/molecularDynamics/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulator.H
index a42f75681d72dfa26d70464355a7526c1c388187..27bd52771ccf64271138184a77a7ec06f980d8d3 100755
--- a/src/lagrangian/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulator.H
+++ b/src/lagrangian/molecularDynamics/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulator.H
@@ -123,7 +123,7 @@ public:
 
         // Access
 
-            inline const label averagesTaken() const;
+            inline label averagesTaken() const;
 
             inline label nBuffers() const;
 
diff --git a/src/lagrangian/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulatorI.H b/src/lagrangian/molecularDynamics/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulatorI.H
similarity index 97%
rename from src/lagrangian/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulatorI.H
rename to src/lagrangian/molecularDynamics/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulatorI.H
index 7f27130f2bc25214aeac54ac6244558f543ed449..cb3bcc617e478bdbd46037cbdbf8afba4bf00a9d 100755
--- a/src/lagrangian/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulatorI.H
+++ b/src/lagrangian/molecularDynamics/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulatorI.H
@@ -46,7 +46,7 @@ inline const Field<Type>& bufferedAccumulator<Type>::accumulationBuffer() const
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-inline const label bufferedAccumulator<Type>::averagesTaken() const
+inline label bufferedAccumulator<Type>::averagesTaken() const
 {
     return averagesTaken_;
 }
diff --git a/src/lagrangian/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulatorIO.C b/src/lagrangian/molecularDynamics/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulatorIO.C
similarity index 100%
rename from src/lagrangian/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulatorIO.C
rename to src/lagrangian/molecularDynamics/molecule/correlationFunction/bufferedAccumulator/bufferedAccumulatorIO.C
diff --git a/src/lagrangian/molecule/correlationFunction/correlationFunction.C b/src/lagrangian/molecularDynamics/molecule/correlationFunction/correlationFunction.C
similarity index 100%
rename from src/lagrangian/molecule/correlationFunction/correlationFunction.C
rename to src/lagrangian/molecularDynamics/molecule/correlationFunction/correlationFunction.C
diff --git a/src/lagrangian/molecule/correlationFunction/correlationFunction.H b/src/lagrangian/molecularDynamics/molecule/correlationFunction/correlationFunction.H
similarity index 95%
rename from src/lagrangian/molecule/correlationFunction/correlationFunction.H
rename to src/lagrangian/molecularDynamics/molecule/correlationFunction/correlationFunction.H
index f5fae016cdff177f5f42c795f5be448eb1486f76..fcd023e36374f10fa089269afab0407cfd5466b3 100644
--- a/src/lagrangian/molecule/correlationFunction/correlationFunction.H
+++ b/src/lagrangian/molecularDynamics/molecule/correlationFunction/correlationFunction.H
@@ -141,13 +141,13 @@ public:
 
         inline const Field< Field<Type> >& tZeroBuffers() const;
 
-        inline const scalar duration() const;
+        inline scalar duration() const;
 
-        inline const scalar sampleInterval() const;
+        inline scalar sampleInterval() const;
 
-        inline const scalar averagingInterval() const;
+        inline scalar averagingInterval() const;
 
-        inline const label sampleSteps() const;
+        inline label sampleSteps() const;
 
         inline label measurandFieldSize() const;
 
diff --git a/src/lagrangian/molecule/correlationFunction/correlationFunctionI.H b/src/lagrangian/molecularDynamics/molecule/correlationFunction/correlationFunctionI.H
similarity index 86%
rename from src/lagrangian/molecule/correlationFunction/correlationFunctionI.H
rename to src/lagrangian/molecularDynamics/molecule/correlationFunction/correlationFunctionI.H
index 2c2a3743c01aeabfd85d3ebd991710ed90656a4d..6858f647a64478873d3b3a5d35c97033db3bd0ca 100644
--- a/src/lagrangian/molecule/correlationFunction/correlationFunctionI.H
+++ b/src/lagrangian/molecularDynamics/molecule/correlationFunction/correlationFunctionI.H
@@ -33,28 +33,28 @@ tZeroBuffers() const
 
 
 template<class Type>
-inline const scalar Foam::correlationFunction<Type>::duration() const
+inline scalar Foam::correlationFunction<Type>::duration() const
 {
     return duration_;
 }
 
 
 template<class Type>
-inline const scalar Foam::correlationFunction<Type>::sampleInterval() const
+inline scalar Foam::correlationFunction<Type>::sampleInterval() const
 {
     return sampleInterval_;
 }
 
 
 template<class Type>
-inline const scalar Foam::correlationFunction<Type>::averagingInterval() const
+inline scalar Foam::correlationFunction<Type>::averagingInterval() const
 {
     return averagingInterval_;
 }
 
 
 template<class Type>
-inline const label Foam::correlationFunction<Type>::sampleSteps() const
+inline label Foam::correlationFunction<Type>::sampleSteps() const
 {
     return sampleSteps_;
 }
diff --git a/src/lagrangian/molecule/correlationFunction/correlationFunctionIO.C b/src/lagrangian/molecularDynamics/molecule/correlationFunction/correlationFunctionIO.C
similarity index 100%
rename from src/lagrangian/molecule/correlationFunction/correlationFunctionIO.C
rename to src/lagrangian/molecularDynamics/molecule/correlationFunction/correlationFunctionIO.C
diff --git a/src/lagrangian/molecule/distribution/distribution.C b/src/lagrangian/molecularDynamics/molecule/distribution/distribution.C
similarity index 100%
rename from src/lagrangian/molecule/distribution/distribution.C
rename to src/lagrangian/molecularDynamics/molecule/distribution/distribution.C
diff --git a/src/lagrangian/molecule/distribution/distribution.H b/src/lagrangian/molecularDynamics/molecule/distribution/distribution.H
similarity index 98%
rename from src/lagrangian/molecule/distribution/distribution.H
rename to src/lagrangian/molecularDynamics/molecule/distribution/distribution.H
index 8cefd24a81a839b9fdb2e784bb662936db26d78d..e9bfdfccf82af2bc6fee244e15fbf0f9cd7e4973 100755
--- a/src/lagrangian/molecule/distribution/distribution.H
+++ b/src/lagrangian/molecularDynamics/molecule/distribution/distribution.H
@@ -104,7 +104,7 @@ public:
 
         // Access
 
-            inline const scalar binWidth() const;
+            inline scalar binWidth() const;
 
 
     // Member Operators
diff --git a/src/lagrangian/molecule/distribution/distributionI.H b/src/lagrangian/molecularDynamics/molecule/distribution/distributionI.H
similarity index 96%
rename from src/lagrangian/molecule/distribution/distributionI.H
rename to src/lagrangian/molecularDynamics/molecule/distribution/distributionI.H
index fed3118eb5b4f0ae25ed4f1fce018e8965e03273..85d5934da00d4f938dba4aadd3187762f2b5a012 100755
--- a/src/lagrangian/molecule/distribution/distributionI.H
+++ b/src/lagrangian/molecularDynamics/molecule/distribution/distributionI.H
@@ -29,7 +29,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline const scalar distribution::binWidth() const
+inline scalar distribution::binWidth() const
 {
     return binWidth_;
 }
diff --git a/src/lagrangian/molecule/distribution/distributionIO.C b/src/lagrangian/molecularDynamics/molecule/distribution/distributionIO.C
similarity index 100%
rename from src/lagrangian/molecule/distribution/distributionIO.C
rename to src/lagrangian/molecularDynamics/molecule/distribution/distributionIO.C
diff --git a/src/lagrangian/molecule/mdTools/averageMDFields.H b/src/lagrangian/molecularDynamics/molecule/mdTools/averageMDFields.H
similarity index 100%
rename from src/lagrangian/molecule/mdTools/averageMDFields.H
rename to src/lagrangian/molecularDynamics/molecule/mdTools/averageMDFields.H
diff --git a/src/lagrangian/molecule/mdTools/calculateAutoCorrelationFunctions.H b/src/lagrangian/molecularDynamics/molecule/mdTools/calculateAutoCorrelationFunctions.H
similarity index 84%
rename from src/lagrangian/molecule/mdTools/calculateAutoCorrelationFunctions.H
rename to src/lagrangian/molecularDynamics/molecule/mdTools/calculateAutoCorrelationFunctions.H
index a7c4c326bb3c4dc17aa510f095bc44054bce25ab..8f7574a07727601831c71e016ff55f246b856d49 100644
--- a/src/lagrangian/molecule/mdTools/calculateAutoCorrelationFunctions.H
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/calculateAutoCorrelationFunctions.H
@@ -79,5 +79,27 @@ if (mesh.time().timeIndex() % pacf.sampleSteps() == 0)
 
 if (mesh.time().timeIndex() % hfacf.sampleSteps() == 0)
 {
-//     hFacf.calculateCorrelationFunction();
+
+    IDLList<molecule>::iterator mol(molecules.begin());
+
+    vector s = vector::zero;
+
+    for
+    (
+        mol = molecules.begin();
+        mol != molecules.end();
+        ++mol
+    )
+    {
+        s +=
+        (
+            0.5 * mol().mass() * magSqr(mol().U())
+            +
+            mol().potentialEnergy()
+        ) * mol().U()
+        +
+        0.5 * ( mol().rf() & mol().U() );
+    }
+
+    hfacf.calculateCorrelationFunction(s);
 }
diff --git a/src/lagrangian/molecule/mdTools/calculateMDFields.H b/src/lagrangian/molecularDynamics/molecule/mdTools/calculateMDFields.H
similarity index 100%
rename from src/lagrangian/molecule/mdTools/calculateMDFields.H
rename to src/lagrangian/molecularDynamics/molecule/mdTools/calculateMDFields.H
diff --git a/src/lagrangian/molecule/mdTools/calculateTransportProperties.H b/src/lagrangian/molecularDynamics/molecule/mdTools/calculateTransportProperties.H
similarity index 81%
rename from src/lagrangian/molecule/mdTools/calculateTransportProperties.H
rename to src/lagrangian/molecularDynamics/molecule/mdTools/calculateTransportProperties.H
index 15e0dab7e34c468856704e799a537d1889169ee7..b0923c42db8453dee655139be69644e93c81932e 100644
--- a/src/lagrangian/molecule/mdTools/calculateTransportProperties.H
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/calculateTransportProperties.H
@@ -59,7 +59,24 @@ Info<< "Viscosity = "
 
 if(writeHFacf)
 {
-
+    OFstream hfacfFile
+    (
+        runTime.path()/ + "hfacf"
+    );
+    
+    if (!hfacf.writeAveraged(hfacfFile))
+    {
+        FatalErrorIn(args.executable())
+            << "Failed writing to "
+            << hfacfFile.name()
+            << abort(FatalError);
+    }
 }
 
-
+Info << "Thermal conductivity = "
+    << hfacf.integral()
+        /averageTemperature
+        /averageTemperature
+        /moleculeCloud::kb
+        / meshVolume
+    << endl;
diff --git a/src/lagrangian/molecule/mdTools/createAutoCorrelationFunctions.H b/src/lagrangian/molecularDynamics/molecule/mdTools/createAutoCorrelationFunctions.H
similarity index 100%
rename from src/lagrangian/molecule/mdTools/createAutoCorrelationFunctions.H
rename to src/lagrangian/molecularDynamics/molecule/mdTools/createAutoCorrelationFunctions.H
diff --git a/src/lagrangian/molecule/mdTools/createMDFields.H b/src/lagrangian/molecularDynamics/molecule/mdTools/createMDFields.H
similarity index 80%
rename from src/lagrangian/molecule/mdTools/createMDFields.H
rename to src/lagrangian/molecularDynamics/molecule/mdTools/createMDFields.H
index f5742fd2c8fd89afdb9428abb5b90ff3e6523967..dccd1891a7fa966df363fa44687bf654571bed3a 100755
--- a/src/lagrangian/molecule/mdTools/createMDFields.H
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/createMDFields.H
@@ -2,25 +2,25 @@
 
 List< scalarField > allSpeciesN_RU
 (
-    molecules.nIds(),
+    molecules.pairPotentials().nIds(),
     scalarField (mesh.nCells(), 0.0)
 );
 
 List< scalarField > allSpeciesM_RU
 (
-    molecules.nIds(),
+    molecules.pairPotentials().nIds(),
     scalarField (mesh.nCells(), 0.0)
 );
 
 List< vectorField > allSpeciesVelocitySum_RU
 (
-    molecules.nIds(),
+    molecules.pairPotentials().nIds(),
     vectorField (mesh.nCells(), vector::zero)
 );
 
 List< scalarField > allSpeciesVelocityMagSquaredSum_RU
 (
-    molecules.nIds(),
+    molecules.pairPotentials().nIds(),
     scalarField (mesh.nCells(), 0.0)
 );
 
@@ -34,13 +34,13 @@ Info << nl << "Creating fields." << endl;
 
 PtrList<volScalarField> allSpeciesRhoN
 (
-    molecules.nIds()
+    molecules.pairPotentials().nIds()
 );
 
 forAll (allSpeciesRhoN, rN)
 {
     Info << "    Creating number density field for "
-        << molecules.idList()[rN] << endl;
+        << molecules.pairPotentials().idList()[rN] << endl;
 
     allSpeciesRhoN.set
     (
@@ -49,7 +49,7 @@ forAll (allSpeciesRhoN, rN)
         (
             IOobject
             (
-                "rhoN_" + molecules.idList()[rN],
+                "rhoN_" + molecules.pairPotentials().idList()[rN],
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,
@@ -89,13 +89,13 @@ totalRhoN.correctBoundaryConditions();
 
 PtrList<volScalarField> allSpeciesRhoM
 (
-    molecules.nIds()
+    molecules.pairPotentials().nIds()
 );
 
 forAll (allSpeciesRhoM, rM)
 {
     Info << "    Creating mass density field for "
-        << molecules.idList()[rM] << endl;
+        << molecules.pairPotentials().idList()[rM] << endl;
 
     allSpeciesRhoM.set
     (
@@ -104,7 +104,7 @@ forAll (allSpeciesRhoM, rM)
         (
             IOobject
             (
-                "rhoM_" + molecules.idList()[rM],
+                "rhoM_" + molecules.pairPotentials().idList()[rM],
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,
@@ -144,13 +144,13 @@ totalRhoM.correctBoundaryConditions();
 
 PtrList<volVectorField> allSpeciesVelocity
 (
-    molecules.nIds()
+    molecules.pairPotentials().nIds()
 );
 
 forAll (allSpeciesVelocity, v)
 {
     Info << "    Creating velocity field for "
-        << molecules.idList()[v] << endl;
+        << molecules.pairPotentials().idList()[v] << endl;
 
     allSpeciesVelocity.set
     (
@@ -159,7 +159,7 @@ forAll (allSpeciesVelocity, v)
         (
             IOobject
             (
-                "velocity_" + molecules.idList()[v],
+                "velocity_" + molecules.pairPotentials().idList()[v],
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,
@@ -177,22 +177,40 @@ forAll (allSpeciesVelocity, v)
 
 Info << "    Creating total velocity field" << endl;
 
+// volVectorField totalVelocity
+// (
+//     IOobject
+//     (
+//         "velocity_total",
+//         runTime.timeName(),
+//         mesh,
+//         IOobject::NO_READ,
+//         IOobject::AUTO_WRITE
+//     ),
+//     mesh,
+//     dimVelocity,
+//     "zeroGradient"
+// );
+// totalVelocity.internalField() = vectorField (mesh.nCells(), vector::zero);
+// totalVelocity.correctBoundaryConditions();
+
+
 volVectorField totalVelocity
+
 (
     IOobject
     (
-        "velocity_total",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
+
+    "velocity_total",
+    runTime.timeName(),
+    mesh,
+    IOobject::NO_READ,
+    IOobject::AUTO_WRITE
+
     ),
     mesh,
-    dimVelocity,
-    "zeroGradient"
+    dimensionedVector("zero", dimVelocity, vector::zero)
 );
-totalVelocity.internalField() = vectorField (mesh.nCells(), vector::zero);
-totalVelocity.correctBoundaryConditions();
 
 /*---------------------------------------------------------------------------*\
     Kinetic temperature
@@ -200,13 +218,13 @@ totalVelocity.correctBoundaryConditions();
 
 PtrList<volScalarField> allSpeciesTemperature
 (
-    molecules.nIds()
+    molecules.pairPotentials().nIds()
 );
 
 forAll (allSpeciesTemperature, t)
 {
     Info << "    Creating temperature field for "
-        << molecules.idList()[t] << endl;
+        << molecules.pairPotentials().idList()[t] << endl;
 
     allSpeciesTemperature.set
     (
@@ -215,7 +233,7 @@ forAll (allSpeciesTemperature, t)
         (
             IOobject
             (
-                "temperature_" + molecules.idList()[t],
+                "temperature_" + molecules.pairPotentials().idList()[t],
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,
@@ -256,13 +274,13 @@ totalTemperature.correctBoundaryConditions();
 
 PtrList<volScalarField> allSpeciesMeanKE
 (
-    molecules.nIds()
+    molecules.pairPotentials().nIds()
 );
 
 forAll (allSpeciesMeanKE, mKE)
 {
     Info << "    Creating mean kinetic energy field for "
-        << molecules.idList()[mKE] << endl;
+        << molecules.pairPotentials().idList()[mKE] << endl;
 
     allSpeciesMeanKE.set
     (
@@ -271,7 +289,7 @@ forAll (allSpeciesMeanKE, mKE)
         (
             IOobject
             (
-                "meanKE_" + molecules.idList()[mKE],
+                "meanKE_" + molecules.pairPotentials().idList()[mKE],
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,
diff --git a/src/lagrangian/molecule/mdTools/createRefUnits.H b/src/lagrangian/molecularDynamics/molecule/mdTools/createRefUnits.H
similarity index 100%
rename from src/lagrangian/molecule/mdTools/createRefUnits.H
rename to src/lagrangian/molecularDynamics/molecule/mdTools/createRefUnits.H
diff --git a/src/lagrangian/molecule/mdTools/md.H b/src/lagrangian/molecularDynamics/molecule/mdTools/md.H
similarity index 100%
rename from src/lagrangian/molecule/mdTools/md.H
rename to src/lagrangian/molecularDynamics/molecule/mdTools/md.H
diff --git a/src/lagrangian/molecule/mdTools/meanMomentumEnergyAndNMols.H b/src/lagrangian/molecularDynamics/molecule/mdTools/meanMomentumEnergyAndNMols.H
similarity index 100%
rename from src/lagrangian/molecule/mdTools/meanMomentumEnergyAndNMols.H
rename to src/lagrangian/molecularDynamics/molecule/mdTools/meanMomentumEnergyAndNMols.H
diff --git a/src/lagrangian/molecule/mdTools/resetMDFields.H b/src/lagrangian/molecularDynamics/molecule/mdTools/resetMDFields.H
similarity index 72%
rename from src/lagrangian/molecule/mdTools/resetMDFields.H
rename to src/lagrangian/molecularDynamics/molecule/mdTools/resetMDFields.H
index e5fdbbd8a46b4947abf691e179ca205c3ab6f03a..337647ec02207b6c8662142244664e4c9350d447 100755
--- a/src/lagrangian/molecule/mdTools/resetMDFields.H
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/resetMDFields.H
@@ -2,25 +2,25 @@ if (runTime.outputTime())
 {
     allSpeciesN_RU = List< scalarField >
     (
-        molecules.nIds(),
+        molecules.pairPotentials().nIds(),
         scalarField (mesh.nCells(), 0.0)
     );
 
     allSpeciesM_RU = List< scalarField >
     (
-        molecules.nIds(),
+        molecules.pairPotentials().nIds(),
         scalarField (mesh.nCells(), 0.0)
     );
 
     allSpeciesVelocitySum_RU = List< vectorField >
     (
-        molecules.nIds(),
+        molecules.pairPotentials().nIds(),
         vectorField (mesh.nCells(), vector::zero)
     );
 
     allSpeciesVelocityMagSquaredSum_RU = List< scalarField >
     (
-        molecules.nIds(),
+        molecules.pairPotentials().nIds(),
         scalarField (mesh.nCells(), 0.0)
     );
 }
diff --git a/src/lagrangian/molecule/mdTools/temperatureAndPressure.H b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H
similarity index 97%
rename from src/lagrangian/molecule/mdTools/temperatureAndPressure.H
rename to src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H
index 0563700cb73d1008a4ce10857143e907abb3b29e..0ae8b5d1eac1f93a8f1b478a0cbd602f8df89739 100755
--- a/src/lagrangian/molecule/mdTools/temperatureAndPressure.H
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H
@@ -86,6 +86,8 @@ if (runTime.outputTime())
                 << " m/s" << nl
             << "Average temperature = "
                 << averageTemperature << " K" << nl
+            << "accumulatedTotalrDotfSum = "
+                << accumulatedTotalrDotfSum << nl
             << "Average pressure = "
                 << averagePressure << " N/m^2" << nl
             << "----------------------------------------" << endl;
diff --git a/src/lagrangian/molecule/mdTools/temperatureAndPressureVariables.H b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressureVariables.H
similarity index 100%
rename from src/lagrangian/molecule/mdTools/temperatureAndPressureVariables.H
rename to src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressureVariables.H
diff --git a/src/lagrangian/molecule/mdTools/temperatureEquilibration.H b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureEquilibration.H
similarity index 100%
rename from src/lagrangian/molecule/mdTools/temperatureEquilibration.H
rename to src/lagrangian/molecularDynamics/molecule/mdTools/temperatureEquilibration.H
diff --git a/src/lagrangian/molecule/molecule/molecule.C b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
similarity index 100%
rename from src/lagrangian/molecule/molecule/molecule.C
rename to src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
diff --git a/src/lagrangian/molecule/molecule/molecule.H b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
similarity index 96%
rename from src/lagrangian/molecule/molecule/molecule.H
rename to src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
index 68cff646b1a2718020ea785d3ec68b5d976b1d10..cbc0a7a134618552a44adeee06d9b00e77e1da4a 100755
--- a/src/lagrangian/molecule/molecule/molecule.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
@@ -23,10 +23,9 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::molecule
+    molecule
 
 Description
-    Foam::molecule
 
 SourceFiles
     moleculeI.H
@@ -166,10 +165,10 @@ public:
         // Access
 
             //- Return id
-            inline const label id() const;
+            inline label id() const;
 
             //- Return mass
-            inline const scalar mass() const;
+            inline scalar mass() const;
 
             //- Return velocity
             inline const vector& U() const;
@@ -180,7 +179,7 @@ public:
             inline vector& A();
 
             //- Return potential energy
-            inline const scalar potentialEnergy() const;
+            inline scalar potentialEnergy() const;
             inline scalar& potentialEnergy();
 
             //- Return stress contribution
@@ -188,10 +187,11 @@ public:
             inline tensor& rf();
 
             //- Return tethered
-            inline const label tethered() const;
+            inline label tethered() const;
 
             //- Return tetherPosition
             inline const vector& tetherPosition() const;
+            inline vector& tetherPosition();
 
 
         //- Tracking
diff --git a/src/lagrangian/molecule/molecule/moleculeI.H b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
similarity index 92%
rename from src/lagrangian/molecule/molecule/moleculeI.H
rename to src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
index 6dd50e337e6123d21ee54068aef697e0b6e51815..c9f1b9eb73925fd1d227fa8d071a76fe26346b02 100755
--- a/src/lagrangian/molecule/molecule/moleculeI.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
@@ -68,13 +68,13 @@ inline molecule::trackData::trackData
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline const label molecule::id() const
+inline label molecule::id() const
 {
     return id_;
 }
 
 
-inline const scalar molecule::mass() const
+inline scalar molecule::mass() const
 {
     return mass_;
 }
@@ -104,7 +104,7 @@ inline vector& molecule::A()
 }
 
 
-inline const scalar molecule::potentialEnergy() const
+inline scalar molecule::potentialEnergy() const
 {
     return potentialEnergy_;
 }
@@ -128,7 +128,7 @@ inline tensor& molecule::rf()
 }
 
 
-inline const label molecule::tethered() const
+inline label molecule::tethered() const
 {
     return tethered_;
 }
@@ -140,6 +140,12 @@ inline const vector& molecule::tetherPosition() const
 }
 
 
+inline vector& molecule::tetherPosition()
+{
+    return tetherPosition_;
+}
+
+
 inline moleculeCloud& molecule::trackData::molCloud()
 {
     return molCloud_;
diff --git a/src/lagrangian/molecule/molecule/moleculeIO.C b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C
similarity index 100%
rename from src/lagrangian/molecule/molecule/moleculeIO.C
rename to src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloud.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
similarity index 100%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloud.C
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloud.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
similarity index 94%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloud.H
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
index 5175fdd0d80b4a8f298847ace66b2c5cddcaf2e5..600c874346f0f3c288ebd8d06af414f70648901e 100755
--- a/src/lagrangian/molecule/moleculeCloud/moleculeCloud.H
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
@@ -97,18 +97,6 @@ private:
 
             scalar potentialEnergyLimit_;
 
-            scalar guardRadius_;
-
-            scalar rCutMax_;
-
-            //- storing rCutMaxSqr in class as well as rCutMax to
-            //- avoid needing to calculate it often.
-            //- Possibilty of inconsistency if tinkered with.
-
-            scalar rCutMaxSqr_;
-
-            List<word> idList_;
-
             labelList removalOrder_;
 
             labelListList directInteractionList_;
@@ -189,20 +177,10 @@ public:
 
             inline scalar potentialEnergyLimit() const;
 
-            inline scalar guardRadius() const;
-
-            inline const scalar rCutMax() const;
-
-            inline const scalar rCutMaxSqr() const;
-
-            inline const List<word>& idList() const;
-
             inline const labelList& removalOrder() const;
 
             inline label nPairPotentials() const;
 
-            inline label nIds() const;
-
             inline const labelListList& directInteractionList() const;
 
             inline const referredCellList& referredInteractionList() const;
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudApplyConstraintsAndThermostats.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudApplyConstraintsAndThermostats.C
similarity index 100%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudApplyConstraintsAndThermostats.C
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudApplyConstraintsAndThermostats.C
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudBuildCellInteractionLists.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudBuildCellInteractionLists.C
similarity index 100%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudBuildCellInteractionLists.C
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudBuildCellInteractionLists.C
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudBuildCellOccupancy.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudBuildCellOccupancy.C
similarity index 89%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudBuildCellOccupancy.C
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudBuildCellOccupancy.C
index f02be61f7efc30adfca2faee72eea507cd81c983..89ad6cd3e6db527294daf736c1e1f9b46b4f1ba9 100755
--- a/src/lagrangian/molecule/moleculeCloud/moleculeCloudBuildCellOccupancy.C
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudBuildCellOccupancy.C
@@ -22,13 +22,21 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
+Class
+    moleculeCloud
+
+Description
+
 \*----------------------------------------------------------------------------*/
 
 #include "moleculeCloud.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-void Foam::moleculeCloud::buildCellOccupancy()
+namespace Foam
+{
+
+void moleculeCloud::buildCellOccupancy()
 {
     forAll(cellOccupancy_, cO)
     {
@@ -37,7 +45,12 @@ void Foam::moleculeCloud::buildCellOccupancy()
 
     iterator mol(this->begin());
 
-    for (mol = this->begin(); mol != this->end(); ++mol)
+    for
+    (
+        mol = this->begin();
+        mol != this->end();
+        ++mol
+    )
     {
         cellOccupancy_[mol().cell()].append(&mol());
     }
@@ -50,5 +63,6 @@ void Foam::moleculeCloud::buildCellOccupancy()
     referredInteractionList_.referMolecules();
 }
 
+} // End namespace Foam
 
 // ************************************************************************* //
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudBuildCellReferralLists.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudBuildCellReferralLists.C
similarity index 100%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudBuildCellReferralLists.C
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudBuildCellReferralLists.C
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudBuildDirectInteractionList.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudBuildDirectInteractionList.H
similarity index 100%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudBuildDirectInteractionList.H
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudBuildDirectInteractionList.H
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudBuildReferredInteractionList.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudBuildReferredInteractionList.H
similarity index 100%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudBuildReferredInteractionList.H
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudBuildReferredInteractionList.H
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculateExternalForce.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculateExternalForce.C
similarity index 100%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculateExternalForce.C
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculateExternalForce.C
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculateForce.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculateForce.C
similarity index 100%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculateForce.C
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculateForce.C
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculatePairForce.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculatePairForce.C
similarity index 100%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculatePairForce.C
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculatePairForce.C
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculatePairForceRealCells.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculatePairForceRealCells.H
similarity index 100%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculatePairForceRealCells.H
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculatePairForceRealCells.H
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculatePairForceRealCellsCalculationStep.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculatePairForceRealCellsCalculationStep.H
similarity index 92%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculatePairForceRealCellsCalculationStep.H
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculatePairForceRealCellsCalculationStep.H
index 49e1e2d781b2462e80d44c6133be77acd9c5cdf9..9f2080d02e555d3922a9a7bd3d03af188473ca28 100755
--- a/src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculatePairForceRealCellsCalculationStep.H
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculatePairForceRealCellsCalculationStep.H
@@ -6,7 +6,7 @@ rIJ = molI->position() - molJ->position();
 
 rIJMagSq = magSqr(rIJ);
 
-if (pairPotentials_.rCutSqr(idI, idJ, rIJMagSq))
+if(pairPotentials_.rCutSqr(idI, idJ, rIJMagSq))
 {
     rIJMag = mag(rIJ);
 
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculatePairForceReferredCells.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculatePairForceReferredCells.H
similarity index 100%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculatePairForceReferredCells.H
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculatePairForceReferredCells.H
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculateTetherForce.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculateTetherForce.C
similarity index 100%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudCalculateTetherForce.C
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculateTetherForce.C
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCodeSnippets.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCodeSnippets.H
new file mode 100755
index 0000000000000000000000000000000000000000..a66b409b450acb56cc0df28b36e6b955d76d646c
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCodeSnippets.H
@@ -0,0 +1,193 @@
+// Parallel coding to access boundary information to build up interaction cell info
+
+// See preservePatchTypes for how to read the boundary file.
+
+// Read faceProcAddressing, as per reconstructPar, to get hold of the original,
+// undecomposed face label from a face on a processor mesh.  See email from Mattijs:
+
+// > Is it a case of reading the faceProcAddressing file, in the same way as
+// > something like reconstructPar?
+//         Correct.
+// 
+//         Note that faceProcAddressing is a bit weird since it also includes which side 
+// of an internal face we have. If I remember correctly:
+// 
+//         faceI == 0 illegal
+//         faceI > 0 we have the original owner of  faceI-1 i.e. we have the face in the 
+//         original order.
+//         faceI < 0 we have the original neighbour of -faceI-1 so the face is flipped.
+
+// Use the same functionality as
+// label polyBoundaryMesh::whichPatch(const label faceIndex) const
+// To determine which patch a face was on originally.
+
+if (Pstream::parRun())
+{
+//     if (Pstream::myProcNo() == Pstream::masterNo())
+// //     {
+//         dictionary patchDictionary;
+// 
+//         DynamicList<word> patchNames;
+// 
+//         {
+//             IOobject undecomposedBoundaryHeader
+//             (
+//                 "undecomposedBoundary",
+//                 mesh_.time().constant(),
+//                 polyMesh::meshSubDir,
+//                 mesh_,
+//                 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
+//                 (
+//                     "moleculeCloudBuildCellInteractionLists.C\n"
+//                 )
+//                     << "undecomposedBoundary file not found in "
+//                         "constant/polyMesh"
+//                     << abort(FatalError);
+//             }
+//         }
+// 
+//         labelIOList faceProcAddressing
+//         (
+//             IOobject
+//             (
+//                 "faceProcAddressing",
+//                 mesh_.time().constant(),
+//                 polyMesh::meshSubDir,
+//                 mesh_,
+//                 IOobject::MUST_READ,
+//                 IOobject::NO_WRITE,
+//                 false
+//             )
+//         );
+
+        labelList procPatches(mesh_.globalData().processorPatches());
+
+        forAll(procPatches,pP)
+        {
+            const processorPolyPatch& patch =
+                refCast<const processorPolyPatch>
+            (
+                mesh_.boundaryMesh()[procPatches[pP]]
+            );
+// 
+//             Pout << nl << "name: " << patch.name() << nl
+//                     << "start: " << patch.start() << nl
+//                     << "size: " << patch.size() << nl
+//                     << "separated: " << Switch(patch.separated()) << nl
+//                     << "parallel: " << Switch(patch.parallel()) << nl << endl;
+// 
+//             forAll (patch, pI)
+//             {
+//                 label decomposedMeshFace = patch.start() + pI;
+// 
+//                 label faceProcAdd = faceProcAddressing[decomposedMeshFace];
+// 
+//                 label globalFace = abs(faceProcAdd)-1;
+// 
+//                 Pout << "Patch index: " << pI
+//                         << " " << patch[pI]
+//                         << " Mesh index: " << decomposedMeshFace
+//                         << " faceProcAdd: " << faceProcAdd
+//                         << " globalFace:" << globalFace;
+// 
+//                 label minStart = -1;
+// 
+//                 // Scanning the dictionary each time is a very ugly way of
+//                 // finding out what patch a face originally belonged to, but
+//                 // it proves the concept.  Read the patch info a container
+//                 // class and have a neat way of tell which patch a face is from
+//                 // embedded in that.  Split each processor face down into
+//                 // separate lists for each different originiating patch.
+// 
+//                 forAll(patchNames, patchi)
+//                 {
+//                     if (patchDictionary.found(patchNames[patchi]))
+//                     {
+//                         const dictionary& patchDict = 
+//                                 patchDictionary.subDict(patchNames[patchi]);
+// 
+//                         word faceName(patchNames[patchi]);
+//                         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
+//                         )
+//                         {
+//                             Pout << " original patch: " << faceName  << endl;
+//                         }
+//                     }
+//                 }
+// 
+//                 if (globalFace < minStart)
+//                 {
+//                     Pout << " originally an internal face" << endl;
+//                 }
+//             }
+// 
+            if (patch.separated())
+            {
+                 Pout << patch.separation();
+            }
+
+            if (!patch.parallel())
+            {
+                 Pout << patch.forwardT();
+            }
+        }
+//     }
+//     else
+//     {
+// 
+//     }
+
+        // Get coords of my shared points
+//     vector sharedPoints(vector::one*(Pstream::myProcNo()+1));
+//     label testRedLab(Pstream::myProcNo()+1);
+    
+//     Pout << testRedLab << endl;
+
+        // Append from all processors
+//     combineReduce(sharedPoints, plusEqOp<vector>());
+//     reduce(testRedLab, plusOp<label>());
+
+//     Pout << testRedLab << endl;
+}
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudI.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
similarity index 90%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudI.H
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
index 2d74b7b4b87b38388d0b44ab919721be63e3d70d..215f59d30938ddc1c5e6c0208aeee1b6642c1293 100755
--- a/src/lagrangian/molecule/moleculeCloud/moleculeCloudI.H
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
@@ -50,42 +50,12 @@ inline scalar moleculeCloud::potentialEnergyLimit() const
 }
 
 
-inline scalar moleculeCloud::guardRadius() const
-{
-    return guardRadius_;
-}
-
-
-inline const scalar moleculeCloud::rCutMax() const
-{
-    return rCutMax_;
-}
-
-
-inline const scalar moleculeCloud::rCutMaxSqr() const
-{
-    return rCutMaxSqr_;
-}
-
-
-inline const List<word>& moleculeCloud::idList() const
-{
-    return idList_;
-}
-
-
 inline label moleculeCloud::nPairPotentials() const
 {
     return pairPotentials_.size();
 }
 
 
-inline label moleculeCloud::nIds() const
-{
-    return idList_.size();
-}
-
-
 inline const labelList& moleculeCloud::removalOrder() const
 {
     return removalOrder_;
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudIntegrateEquationsOfMotion.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudIntegrateEquationsOfMotion.C
similarity index 100%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudIntegrateEquationsOfMotion.C
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudIntegrateEquationsOfMotion.C
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudReadMDParameters.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudReadMDParameters.H
new file mode 100755
index 0000000000000000000000000000000000000000..b74f9c5abf44570bec494bc046599e29304c544c
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudReadMDParameters.H
@@ -0,0 +1,177 @@
+Info<< nl <<  "Reading MD solution parameters:" << endl;
+
+IOdictionary mdSolution
+(
+    IOobject
+    (
+        "mdSolution",
+        mesh_.time().system(),
+        mesh_,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+integrationMethod_ = integrationMethodNames_.read
+(
+    mdSolution.lookup("integrationMethod")
+);
+
+potentialEnergyLimit_ = readScalar
+(
+    mdSolution.lookup("potentialEnergyLimit")
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+IOdictionary potentialDict
+(
+    IOobject
+    (
+        "potentialDict",
+        mesh_.time().system(),
+        mesh_,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+IOdictionary idListDict
+(
+    IOobject
+    (
+        "idList",
+        mesh_.time().constant(),
+        mesh_,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+// ****************************************************************************
+// Pair potentials
+
+if (!potentialDict.found("pair"))
+{
+    FatalErrorIn("moleculeCloudReadPotentials.H") << nl
+            << "pair potential specification subDict not found"
+            << abort(FatalError);
+}
+
+const dictionary& pairDict = potentialDict.subDict("pair");
+
+pairPotentials_.buildPotentials(idListDict, pairDict, mesh_);
+
+if (potentialDict.found("removalOrder"))
+{
+    List<word> remOrd = potentialDict.lookup("removalOrder");
+
+    removalOrder_.setSize(remOrd.size());
+
+    forAll(removalOrder_, rO)
+    {
+        removalOrder_[rO] = findIndex(pairPotentials_.idList(), remOrd[rO]);
+    }
+}
+
+// ****************************************************************************
+// Tether potentials
+
+iterator mol(this->begin());
+
+DynamicList<label> tetherIds;
+
+for
+(
+ mol = this->begin();
+ mol != this->end();
+ ++mol
+)
+{
+    if (mol().tethered())
+    {
+        if (findIndex(tetherIds, mol().id()) == -1)
+        {
+            tetherIds.append
+            (
+                mol().id()
+            );
+        }
+    }
+}
+
+if (Pstream::parRun())
+{
+    List< labelList > allTetherIds(Pstream::nProcs());
+
+    allTetherIds[Pstream::myProcNo()] = tetherIds;
+
+    Pstream::gatherList(allTetherIds);
+
+    if (Pstream::master())
+    {
+        DynamicList<label> globalTetherIds;
+
+        forAll(allTetherIds, procN)
+        {
+            const labelList& procNTetherIds = allTetherIds[procN];
+
+            forAll(procNTetherIds, id)
+            {
+                if (findIndex(globalTetherIds, procNTetherIds[id]) == -1)
+                {
+                    globalTetherIds.append
+                    (
+                        procNTetherIds[id]
+                    );
+                }
+            }
+        }
+
+        globalTetherIds.shrink();
+
+        tetherIds = globalTetherIds;
+    }
+
+    Pstream::scatter(tetherIds);
+}
+
+tetherIds.shrink();
+
+if (tetherIds.size())
+{
+    if (!potentialDict.found("tether"))
+    {
+        FatalErrorIn("moleculeCloudReadPotentials.H") << nl
+                << "tether potential specification subDict not found"
+                << abort(FatalError);
+    }
+
+    const dictionary& tetherDict = potentialDict.subDict("tether");
+
+    tetherPotentials_.buildPotentials(idListDict, tetherDict, tetherIds);
+}
+
+// ****************************************************************************
+// External Forces
+
+gravity_ = vector::zero;
+
+if (potentialDict.found("external"))
+{
+
+    Info << nl << "Reading external forces:" << endl;
+
+    const dictionary& externalDict = potentialDict.subDict("external");
+
+    // ************************************************************************
+    // gravity
+
+    if (externalDict.found("gravity"))
+    {
+        gravity_ = externalDict.lookup("gravity");
+    }
+}
+
+Info << nl << tab << "gravity = " << gravity_ << endl;
+
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudRealCellsInRangeOfSegment.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudRealCellsInRangeOfSegment.C
similarity index 100%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudRealCellsInRangeOfSegment.C
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudRealCellsInRangeOfSegment.C
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudReferredCellsInRangeOfSegment.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudReferredCellsInRangeOfSegment.C
similarity index 100%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudReferredCellsInRangeOfSegment.C
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudReferredCellsInRangeOfSegment.C
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlaps.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlaps.C
similarity index 96%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlaps.C
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlaps.C
index ac4385574fa69fa65a1f2116fff5debe6bbd1a76..e12ea9dbcea3a364f140d24318abc4eddb5e8740 100755
--- a/src/lagrangian/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlaps.C
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlaps.C
@@ -34,7 +34,7 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
 
     forAll(removalOrder_, rO)
     {
-        Info << " " << idList_[removalOrder_[rO]];
+        Info << " " << pairPotentials_.idList()[removalOrder_[rO]];
     }
 
     Info << nl ;
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlapsRealCells.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlapsRealCells.H
similarity index 100%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlapsRealCells.H
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlapsRealCells.H
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlapsRealCellsCalculationStep.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlapsRealCellsCalculationStep.H
similarity index 87%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlapsRealCellsCalculationStep.H
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlapsRealCellsCalculationStep.H
index bc36a600958e02e58c2b33d4b4e828d71063ffa5..04f9ffd8baec9d7f34a3537197b575af40147e95 100755
--- a/src/lagrangian/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlapsRealCellsCalculationStep.H
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlapsRealCellsCalculationStep.H
@@ -12,6 +12,10 @@ if (pairPotentials_.rCutSqr(idI, idJ, rIJMagSq))
 
     bool remove = false;
 
+    // Guard against pairPotentials_.energy being evaluated
+    // if rIJMag < SMALL. A floating point exception will
+    // happen otherwise.
+
     if (rIJMag < SMALL)
     {
         WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
@@ -30,10 +34,14 @@ if (pairPotentials_.rCutSqr(idI, idJ, rIJMagSq))
 
         remove = true;
     }
-
+    
     // Guard against pairPotentials_.energy being evaluated
-    // if rIJMag < SMALL. A floating point exception will
-    // happen otherwise.
+    // if rIJMag < rMin. A tabulation lookup error will occur otherwise.
+
+    if (rIJMag < pairPotentials_.rMin(idI, idJ))
+    {
+        remove = true;
+    }
 
     if (!remove)
     {
@@ -42,6 +50,7 @@ if (pairPotentials_.rCutSqr(idI, idJ, rIJMagSq))
             pairPotentials_.energy(idI, idJ, rIJMag) > potentialEnergyLimit_
         )
         {
+           
             remove = true;
         }
     }
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlapsReferredCells.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlapsReferredCells.H
similarity index 93%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlapsReferredCells.H
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlapsReferredCells.H
index 5b6b76135f2066dd765012858fe1ef12c4a21236..c0914b1a11c7e2f3a4c5d7fb9badf237ec483c1f 100755
--- a/src/lagrangian/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlapsReferredCells.H
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudRemoveHighEnergyOverlapsReferredCells.H
@@ -47,6 +47,10 @@
 
                         bool remove = false;
 
+                        // Guard against pairPotentials_.energy being evaluated
+                        // if rKLMag < SMALL. A floating point exception will
+                        // happen otherwise.
+
                         if (rKLMag < SMALL)
                         {
                             WarningIn
@@ -72,8 +76,12 @@
                         }
 
                         // Guard against pairPotentials_.energy being evaluated
-                        // if rKLMag < SMALL. A floating point exception will
-                        // happen otherwise.
+                        // if rIJMag < rMin. A tubulation lookup error will occur otherwise.
+
+                        if (rKLMag < pairPotentials_.rMin(idK, idL))
+                        {
+                            remove = true;
+                        }
 
                         if (!remove)
                         {
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudTestEdgeEdgeDistance.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudTestEdgeEdgeDistance.C
similarity index 97%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudTestEdgeEdgeDistance.C
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudTestEdgeEdgeDistance.C
index 285447974cf216b8ee508747e96918e9deb6200d..1e8fd26630cb02589d2fa740049321040f9d2c56 100755
--- a/src/lagrangian/molecule/moleculeCloud/moleculeCloudTestEdgeEdgeDistance.C
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudTestEdgeEdgeDistance.C
@@ -75,7 +75,7 @@ bool Foam::moleculeCloud::testEdgeEdgeDistance
          && s <= 1
          && t >= 0
          && t <= 1
-         && magSqr(eIs + a*s - eJs - b*t) <= rCutMaxSqr()
+         && magSqr(eIs + a*s - eJs - b*t) <= pairPotentials_.rCutMaxSqr()
         );
     }
 
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudTestPointFaceDistance.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudTestPointFaceDistance.C
similarity index 93%
rename from src/lagrangian/molecule/moleculeCloud/moleculeCloudTestPointFaceDistance.C
rename to src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudTestPointFaceDistance.C
index e01be864769728a17f9727db4a3fd4c194439e21..03a30873a750e704f7a4c2f437392181e17309a7 100755
--- a/src/lagrangian/molecule/moleculeCloud/moleculeCloudTestPointFaceDistance.C
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudTestPointFaceDistance.C
@@ -128,14 +128,14 @@ bool Foam::moleculeCloud::testPointFaceDistance
 
     scalar perpDist((p - faceC) & faceN);
 
-    if (mag(perpDist) > rCutMax())
+    if (mag(perpDist) > pairPotentials_.rCutMax())
     {
         return false;
     }
 
     vector pointOnPlane = (p - faceN * perpDist);
 
-    if (magSqr(faceC - pointOnPlane) < rCutMaxSqr()*1e-8)
+    if (magSqr(faceC - pointOnPlane) < pairPotentials_.rCutMaxSqr()*1e-8)
     {
         // If pointOnPlane is very close to the face centre
         // then defining the local axes will be inaccurate
@@ -143,7 +143,7 @@ bool Foam::moleculeCloud::testPointFaceDistance
         // inside the face, so return true if the points
         // are in range to be safe
 
-        return (magSqr(pointOnPlane - p) <= rCutMaxSqr());
+        return (magSqr(pointOnPlane - p) <= pairPotentials_.rCutMaxSqr());
     }
 
     vector xAxis = (faceC - pointOnPlane)/mag(faceC - pointOnPlane);
@@ -158,7 +158,7 @@ bool Foam::moleculeCloud::testPointFaceDistance
     {
         const vector& V(points[faceToTest[fTT]]);
 
-        if (magSqr(V-p) <= rCutMaxSqr())
+        if (magSqr(V-p) <= pairPotentials_.rCutMaxSqr())
         {
             return true;
         }
@@ -206,7 +206,7 @@ bool Foam::moleculeCloud::testPointFaceDistance
     if (la_valid < 0)
     {
         // perpendicular point inside face, nearest point is pointOnPlane;
-        return (magSqr(pointOnPlane-p) <= rCutMaxSqr());
+        return (magSqr(pointOnPlane-p) <= pairPotentials_.rCutMaxSqr());
     }
     else
     {
@@ -215,7 +215,7 @@ bool Foam::moleculeCloud::testPointFaceDistance
         return
         (
             magSqr(pointOnPlane + la_valid*(faceC - pointOnPlane) - p)
-         <= rCutMaxSqr()
+         <= pairPotentials_.rCutMaxSqr()
         );
     }
 
diff --git a/src/lagrangian/molecule/reducedUnits/reducedUnits.C b/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnits.C
similarity index 100%
rename from src/lagrangian/molecule/reducedUnits/reducedUnits.C
rename to src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnits.C
diff --git a/src/lagrangian/molecule/reducedUnits/reducedUnits.H b/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnits.H
similarity index 98%
rename from src/lagrangian/molecule/reducedUnits/reducedUnits.H
rename to src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnits.H
index 7adc1051fd414ac75a66e64c79bacada4a5ba045..6a1379ecac206c022e7fd53f6eeae863bb12e25a 100644
--- a/src/lagrangian/molecule/reducedUnits/reducedUnits.H
+++ b/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnits.H
@@ -105,7 +105,7 @@ public:
 
         //- Construct with no argument, uses default values:
         // length  = 1nm
-        // mass = 1.660538782e−27kg (unified atomic mass unit)
+        // mass = 1.660538782e−27kg (unified atomic mass unit)
         // temperature = 1K (therefore, energy = 1*kb)
         reducedUnits();
 
diff --git a/src/lagrangian/molecule/reducedUnits/reducedUnitsI.H b/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnitsI.H
similarity index 100%
rename from src/lagrangian/molecule/reducedUnits/reducedUnitsI.H
rename to src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnitsI.H
diff --git a/src/lagrangian/molecule/reducedUnits/reducedUnitsIO.C b/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnitsIO.C
similarity index 100%
rename from src/lagrangian/molecule/reducedUnits/reducedUnitsIO.C
rename to src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnitsIO.C
diff --git a/src/lagrangian/molecule/referralLists/receivingReferralList.C b/src/lagrangian/molecularDynamics/molecule/referralLists/receivingReferralList.C
similarity index 100%
rename from src/lagrangian/molecule/referralLists/receivingReferralList.C
rename to src/lagrangian/molecularDynamics/molecule/referralLists/receivingReferralList.C
diff --git a/src/lagrangian/molecule/referralLists/receivingReferralList.H b/src/lagrangian/molecularDynamics/molecule/referralLists/receivingReferralList.H
similarity index 100%
rename from src/lagrangian/molecule/referralLists/receivingReferralList.H
rename to src/lagrangian/molecularDynamics/molecule/referralLists/receivingReferralList.H
diff --git a/src/lagrangian/molecule/referralLists/receivingReferralListI.H b/src/lagrangian/molecularDynamics/molecule/referralLists/receivingReferralListI.H
similarity index 100%
rename from src/lagrangian/molecule/referralLists/receivingReferralListI.H
rename to src/lagrangian/molecularDynamics/molecule/referralLists/receivingReferralListI.H
diff --git a/src/lagrangian/molecule/referralLists/sendingReferralList.C b/src/lagrangian/molecularDynamics/molecule/referralLists/sendingReferralList.C
similarity index 100%
rename from src/lagrangian/molecule/referralLists/sendingReferralList.C
rename to src/lagrangian/molecularDynamics/molecule/referralLists/sendingReferralList.C
diff --git a/src/lagrangian/molecule/referralLists/sendingReferralList.H b/src/lagrangian/molecularDynamics/molecule/referralLists/sendingReferralList.H
similarity index 100%
rename from src/lagrangian/molecule/referralLists/sendingReferralList.H
rename to src/lagrangian/molecularDynamics/molecule/referralLists/sendingReferralList.H
diff --git a/src/lagrangian/molecule/referralLists/sendingReferralListI.H b/src/lagrangian/molecularDynamics/molecule/referralLists/sendingReferralListI.H
similarity index 100%
rename from src/lagrangian/molecule/referralLists/sendingReferralListI.H
rename to src/lagrangian/molecularDynamics/molecule/referralLists/sendingReferralListI.H
diff --git a/src/lagrangian/molecule/referredCell/referredCell.C b/src/lagrangian/molecularDynamics/molecule/referredCell/referredCell.C
similarity index 100%
rename from src/lagrangian/molecule/referredCell/referredCell.C
rename to src/lagrangian/molecularDynamics/molecule/referredCell/referredCell.C
diff --git a/src/lagrangian/molecule/referredCell/referredCell.H b/src/lagrangian/molecularDynamics/molecule/referredCell/referredCell.H
similarity index 100%
rename from src/lagrangian/molecule/referredCell/referredCell.H
rename to src/lagrangian/molecularDynamics/molecule/referredCell/referredCell.H
diff --git a/src/lagrangian/molecule/referredCell/referredCellI.H b/src/lagrangian/molecularDynamics/molecule/referredCell/referredCellI.H
similarity index 100%
rename from src/lagrangian/molecule/referredCell/referredCellI.H
rename to src/lagrangian/molecularDynamics/molecule/referredCell/referredCellI.H
diff --git a/src/lagrangian/molecule/referredCellList/referredCellList.C b/src/lagrangian/molecularDynamics/molecule/referredCellList/referredCellList.C
similarity index 100%
rename from src/lagrangian/molecule/referredCellList/referredCellList.C
rename to src/lagrangian/molecularDynamics/molecule/referredCellList/referredCellList.C
diff --git a/src/lagrangian/molecule/referredCellList/referredCellList.H b/src/lagrangian/molecularDynamics/molecule/referredCellList/referredCellList.H
similarity index 100%
rename from src/lagrangian/molecule/referredCellList/referredCellList.H
rename to src/lagrangian/molecularDynamics/molecule/referredCellList/referredCellList.H
diff --git a/src/lagrangian/molecule/referredCellList/referredCellListI.H b/src/lagrangian/molecularDynamics/molecule/referredCellList/referredCellListI.H
similarity index 100%
rename from src/lagrangian/molecule/referredCellList/referredCellListI.H
rename to src/lagrangian/molecularDynamics/molecule/referredCellList/referredCellListI.H
diff --git a/src/lagrangian/molecule/referredMolecule/referredMolecule.C b/src/lagrangian/molecularDynamics/molecule/referredMolecule/referredMolecule.C
similarity index 100%
rename from src/lagrangian/molecule/referredMolecule/referredMolecule.C
rename to src/lagrangian/molecularDynamics/molecule/referredMolecule/referredMolecule.C
diff --git a/src/lagrangian/molecule/referredMolecule/referredMolecule.H b/src/lagrangian/molecularDynamics/molecule/referredMolecule/referredMolecule.H
similarity index 100%
rename from src/lagrangian/molecule/referredMolecule/referredMolecule.H
rename to src/lagrangian/molecularDynamics/molecule/referredMolecule/referredMolecule.H
diff --git a/src/lagrangian/molecule/referredMolecule/referredMoleculeI.H b/src/lagrangian/molecularDynamics/molecule/referredMolecule/referredMoleculeI.H
similarity index 100%
rename from src/lagrangian/molecule/referredMolecule/referredMoleculeI.H
rename to src/lagrangian/molecularDynamics/molecule/referredMolecule/referredMoleculeI.H
diff --git a/src/lagrangian/molecule/potentials/tetherPotential/tetherPotential.C b/src/lagrangian/molecularDynamics/molecule/tetherPotential/tetherPotential.C
similarity index 100%
rename from src/lagrangian/molecule/potentials/tetherPotential/tetherPotential.C
rename to src/lagrangian/molecularDynamics/molecule/tetherPotential/tetherPotential.C
diff --git a/src/lagrangian/molecule/potentials/tetherPotential/tetherPotential.H b/src/lagrangian/molecularDynamics/molecule/tetherPotential/tetherPotential.H
similarity index 100%
rename from src/lagrangian/molecule/potentials/tetherPotential/tetherPotential.H
rename to src/lagrangian/molecularDynamics/molecule/tetherPotential/tetherPotential.H
diff --git a/src/lagrangian/molecule/potentials/tetherPotential/tetherPotentialI.H b/src/lagrangian/molecularDynamics/molecule/tetherPotential/tetherPotentialI.H
similarity index 100%
rename from src/lagrangian/molecule/potentials/tetherPotential/tetherPotentialI.H
rename to src/lagrangian/molecularDynamics/molecule/tetherPotential/tetherPotentialI.H
diff --git a/src/lagrangian/molecule/potentials/tetherPotential/tetherPotentialList.C b/src/lagrangian/molecularDynamics/molecule/tetherPotential/tetherPotentialList.C
similarity index 100%
rename from src/lagrangian/molecule/potentials/tetherPotential/tetherPotentialList.C
rename to src/lagrangian/molecularDynamics/molecule/tetherPotential/tetherPotentialList.C
diff --git a/src/lagrangian/molecule/potentials/tetherPotential/tetherPotentialList.H b/src/lagrangian/molecularDynamics/molecule/tetherPotential/tetherPotentialList.H
similarity index 100%
rename from src/lagrangian/molecule/potentials/tetherPotential/tetherPotentialList.H
rename to src/lagrangian/molecularDynamics/molecule/tetherPotential/tetherPotentialList.H
diff --git a/src/lagrangian/molecule/potentials/tetherPotential/tetherPotentialListI.H b/src/lagrangian/molecularDynamics/molecule/tetherPotential/tetherPotentialListI.H
similarity index 100%
rename from src/lagrangian/molecule/potentials/tetherPotential/tetherPotentialListI.H
rename to src/lagrangian/molecularDynamics/molecule/tetherPotential/tetherPotentialListI.H
diff --git a/src/lagrangian/molecularDynamics/potential/Make/files b/src/lagrangian/molecularDynamics/potential/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..43cf0bb73728b1799b5b1c92f911b36045bf614b
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/Make/files
@@ -0,0 +1,35 @@
+pairPotential = pairPotential
+
+$(pairPotential)/pairPotentialList/pairPotentialList.C
+
+$(pairPotential)/basic/pairPotential.C
+$(pairPotential)/basic/pairPotentialIO.C
+$(pairPotential)/basic/newPairPotential.C
+
+$(pairPotential)/derived/lennardJones/lennardJones.C
+$(pairPotential)/derived/maitlandSmith/maitlandSmith.C
+$(pairPotential)/derived/azizChen/azizChen.C
+
+energyScalingFunction = energyScalingFunction
+
+$(energyScalingFunction)/basic/energyScalingFunction.C
+$(energyScalingFunction)/basic/newEnergyScalingFunction.C
+
+$(energyScalingFunction)/derived/shifted/shifted.C
+$(energyScalingFunction)/derived/shiftedForce/shiftedForce.C
+$(energyScalingFunction)/derived/noScaling/noScaling.C
+$(energyScalingFunction)/derived/sigmoid/sigmoid.C
+$(energyScalingFunction)/derived/doubleSigmoid/doubleSigmoid.C
+
+tetherPotential = tetherPotential
+
+$(tetherPotential)/tetherPotentialList/tetherPotentialList.C
+
+$(tetherPotential)/basic/tetherPotential.C
+$(tetherPotential)/basic/newTetherPotential.C
+
+$(tetherPotential)/derived/harmonicSpring/harmonicSpring.C
+$(tetherPotential)/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.C
+
+LIB = $(FOAM_LIBBIN)/libpotential
+
diff --git a/src/lagrangian/molecularDynamics/potential/Make/options b/src/lagrangian/molecularDynamics/potential/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..a412632ebb87b7a78b4466b6d94fedc04f71922d
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/Make/options
@@ -0,0 +1,6 @@
+EXE_INC = \
+    -I.. \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+LIB_LIBS = \
+    -lfiniteVolume
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunction.C b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunction.C
new file mode 100644
index 0000000000000000000000000000000000000000..1ebd4b4005ddd9cf408978a5092e15281ecff87e
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunction.C
@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    energyScalingFunction
+
+\*---------------------------------------------------------------------------*/
+
+#include "energyScalingFunction.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(energyScalingFunction, 0);
+defineRunTimeSelectionTable(energyScalingFunction, dictionary);
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::energyScalingFunction::energyScalingFunction
+(
+    const word& name,
+    const dictionary& energyScalingFunctionProperties,
+    const pairPotential& pairPot
+)
+:
+    name_(name),
+    energyScalingFunctionProperties_(energyScalingFunctionProperties),
+    pairPot_(pairPot)
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+
+bool Foam::energyScalingFunction::read
+(
+    const dictionary& energyScalingFunctionProperties
+)
+{
+    energyScalingFunctionProperties_ = energyScalingFunctionProperties;
+
+    return true;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunction.H b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunction.H
new file mode 100644
index 0000000000000000000000000000000000000000..e2f38b2a83f1d3ea9a3149b767305bc914af4b47
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunction.H
@@ -0,0 +1,149 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    energyScalingFunction
+
+Description
+
+SourceFiles
+    energyScalingFunction.C
+    newEnergyScalingFunction.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef energyScalingFunction_H
+#define energyScalingFunction_H
+
+#include "IOdictionary.H"
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "autoPtr.H"
+#include "pairPotential.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                Class energyScalingFunction Declaration
+\*---------------------------------------------------------------------------*/
+
+class energyScalingFunction
+{
+
+protected:
+
+    // Protected data
+
+        word name_;
+        
+        dictionary energyScalingFunctionProperties_;
+
+        const pairPotential& pairPot_;
+
+    // Private Member Functions
+
+        //- Disallow copy construct
+        energyScalingFunction(const energyScalingFunction&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const energyScalingFunction&);
+
+public:
+
+    //- Runtime type information
+    TypeName("energyScalingFunction");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            energyScalingFunction,
+            dictionary,
+            (
+                const word& name,
+                const dictionary& energyScalingFunctionProperties,
+                const pairPotential& pairPot
+            ),
+            (name, energyScalingFunctionProperties, pairPot)
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected viscosity model
+        static autoPtr<energyScalingFunction> New
+        (
+            const word& name,
+            const dictionary& energyScalingFunctionProperties,
+            const pairPotential& pairPot
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        energyScalingFunction
+        (
+            const word& name,
+            const dictionary& energyScalingFunctionProperties,
+            const pairPotential& pairPot
+        );
+
+
+    // Destructor
+
+        virtual ~energyScalingFunction()
+        {}
+
+
+    // Member Functions
+
+        virtual void scaleEnergy(scalar& e, const scalar r) const = 0;
+
+        const dictionary& energyScalingFunctionProperties() const
+        {
+            return energyScalingFunctionProperties_;
+        }
+
+        //- Read energyScalingFunction dictionary
+        virtual bool read(const dictionary& energyScalingFunctionProperties) = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/newEnergyScalingFunction.C b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/newEnergyScalingFunction.C
new file mode 100644
index 0000000000000000000000000000000000000000..bb8153a38346700ddb411c519eaddb3c8ac0ac3c
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/newEnergyScalingFunction.C
@@ -0,0 +1,79 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    energyScalingFunction
+
+\*---------------------------------------------------------------------------*/
+
+#include "energyScalingFunction.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+autoPtr<energyScalingFunction> energyScalingFunction::New
+(
+    const word& name,
+    const dictionary& energyScalingFunctionProperties,
+    const pairPotential& pairPot
+)
+{
+    word energyScalingFunctionTypeName
+    (
+        energyScalingFunctionProperties.lookup("energyScalingFunction")
+    );
+
+    Info<< "Selecting energy scaling function "
+        << energyScalingFunctionTypeName << " for "
+        << name << " potential energy." << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(energyScalingFunctionTypeName);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "energyScalingFunction::New()"
+        )   << "Unknown energyScalingFunction type "
+            << energyScalingFunctionTypeName << endl << endl
+            << "Valid  energyScalingFunctions are : " << endl
+            << dictionaryConstructorTablePtr_->toc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<energyScalingFunction>
+        (cstrIter()(name, energyScalingFunctionProperties, pairPot));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.C b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.C
new file mode 100644
index 0000000000000000000000000000000000000000..63fabd6ad9d0e561a97576cda5aa860289e242cb
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.C
@@ -0,0 +1,104 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "doubleSigmoid.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace energyScalingFunctions
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(doubleSigmoid, 0);
+
+addToRunTimeSelectionTable
+(
+    energyScalingFunction,
+    doubleSigmoid,
+    dictionary
+);
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+scalar doubleSigmoid::sigmoidScale
+    (
+        const scalar r,
+        const scalar shift,
+        const scalar scale
+    ) const
+{
+    return 1.0 / (1.0 + exp( scale * (r - shift)));
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+doubleSigmoid::doubleSigmoid
+(
+    const word& name,
+    const dictionary& energyScalingFunctionProperties,
+    const pairPotential& pairPot
+)
+:
+    energyScalingFunction(name, energyScalingFunctionProperties, pairPot),
+    doubleSigmoidCoeffs_(energyScalingFunctionProperties.subDict(typeName + "Coeffs")),
+    shift1_(readScalar(doubleSigmoidCoeffs_.lookup("shift1"))),
+    scale1_(readScalar(doubleSigmoidCoeffs_.lookup("scale1"))),
+    shift2_(readScalar(doubleSigmoidCoeffs_.lookup("shift2"))),
+    scale2_(readScalar(doubleSigmoidCoeffs_.lookup("scale2")))
+{}
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void doubleSigmoid::scaleEnergy(scalar& e, const scalar r) const
+{
+    e *= sigmoidScale(r, shift1_, scale1_) * sigmoidScale(r, shift2_, scale2_);
+}
+
+bool doubleSigmoid::read(const dictionary& energyScalingFunctionProperties)
+{
+    energyScalingFunction::read(energyScalingFunctionProperties);
+
+    doubleSigmoidCoeffs_ = energyScalingFunctionProperties.subDict(typeName + "Coeffs");
+
+    doubleSigmoidCoeffs_.lookup("shift1") >> shift1_;
+    doubleSigmoidCoeffs_.lookup("scale1") >> scale1_;
+    doubleSigmoidCoeffs_.lookup("shift2") >> shift2_;
+    doubleSigmoidCoeffs_.lookup("scale2") >> scale2_;
+    
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace energyScalingFunctions
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.H b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.H
new file mode 100644
index 0000000000000000000000000000000000000000..38dba03a923a178b91667ad83adc03612bba290e
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.H
@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    doubleSigmoid
+
+Description
+     
+
+SourceFiles
+    doubleSigmoid.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef doubleSigmoid_H
+#define doubleSigmoid_H
+
+#include "energyScalingFunction.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace energyScalingFunctions
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class doubleSigmoid Declaration
+\*---------------------------------------------------------------------------*/
+
+class doubleSigmoid
+:
+    public energyScalingFunction
+{
+    // Private data
+        
+        dictionary doubleSigmoidCoeffs_;
+
+        scalar shift1_;
+        scalar scale1_;
+        scalar shift2_;
+        scalar scale2_;
+        
+    // Private Member Functions
+
+        scalar sigmoidScale
+        (
+            const scalar r,
+            const scalar shift,
+            const scalar scale
+        ) const;
+
+public:
+
+    //- Runtime type information
+    TypeName("doubleSigmoid");
+
+
+    // Constructors
+
+        //- Construct from components
+        doubleSigmoid
+        (
+            const word& name,
+            const dictionary& energyScalingFunctionProperties,
+            const pairPotential& pairPot
+        );
+
+
+    // Destructor
+
+        ~doubleSigmoid()
+        {}
+
+    // Member Functions
+
+        void scaleEnergy(scalar& e, const scalar r) const;
+
+        //- Read transportProperties dictionary
+        bool read(const dictionary& energyScalingFunctionProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace energyScalingFunctions
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.C b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.C
new file mode 100644
index 0000000000000000000000000000000000000000..08daf130cbcc23d0b4770260808e73cbf6cb8f16
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.C
@@ -0,0 +1,81 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "noScaling.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace energyScalingFunctions
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(noScaling, 0);
+
+addToRunTimeSelectionTable
+(
+    energyScalingFunction,
+    noScaling,
+    dictionary
+);
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+noScaling::noScaling
+(
+    const word& name,
+    const dictionary& energyScalingFunctionProperties,
+    const pairPotential& pairPot
+)
+:
+    energyScalingFunction(name, energyScalingFunctionProperties, pairPot)
+{}
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void noScaling::scaleEnergy(scalar& e, const scalar r) const
+{}
+
+bool noScaling::read(const dictionary& energyScalingFunctionProperties)
+{
+    energyScalingFunction::read(energyScalingFunctionProperties);
+    
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace energyScalingFunctions
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.H b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.H
new file mode 100644
index 0000000000000000000000000000000000000000..01ee2c2f90d5b037a001254ba77b19b782294f17
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.H
@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    noScaling
+
+Description
+     
+
+SourceFiles
+    noScaling.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef noScaling_H
+#define noScaling_H
+
+#include "energyScalingFunction.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace energyScalingFunctions
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class noScaling Declaration
+\*---------------------------------------------------------------------------*/
+
+class noScaling
+:
+    public energyScalingFunction
+{
+    // Private data
+
+public:
+
+    //- Runtime type information
+    TypeName("noScaling");
+
+
+    // Constructors
+
+        //- Construct from components
+        noScaling
+        (
+            const word& name,
+            const dictionary& energyScalingFunctionProperties,
+            const pairPotential& pairPot
+        );
+
+
+    // Destructor
+
+        ~noScaling()
+        {}
+
+    // Member Functions
+
+        void scaleEnergy(scalar& e, const scalar r) const;
+
+        //- Read transportProperties dictionary
+        bool read(const dictionary& energyScalingFunctionProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace energyScalingFunctions
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.C b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.C
new file mode 100644
index 0000000000000000000000000000000000000000..28882fb1188b349bb679b41e1121a297a598aeac
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.C
@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "shifted.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace energyScalingFunctions
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(shifted, 0);
+
+addToRunTimeSelectionTable
+(
+    energyScalingFunction,
+    shifted,
+    dictionary
+);
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+shifted::shifted
+(
+    const word& name,
+    const dictionary& energyScalingFunctionProperties,
+    const pairPotential& pairPot
+)
+:
+    energyScalingFunction(name, energyScalingFunctionProperties, pairPot),
+    e_at_rCut_(pairPot.unscaledEnergy(pairPot.rCut()))
+{}
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void shifted::scaleEnergy(scalar& e, const scalar r) const
+{
+    e -= e_at_rCut_;
+}
+
+bool shifted::read(const dictionary& energyScalingFunctionProperties)
+{
+    energyScalingFunction::read(energyScalingFunctionProperties);
+    
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace energyScalingFunctions
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.H b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.H
new file mode 100644
index 0000000000000000000000000000000000000000..4d7d77f22e0a62c66b6669adcd290b79c6c7ee22
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.H
@@ -0,0 +1,100 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    shifted
+
+Description
+     
+
+SourceFiles
+    shifted.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef shifted_H
+#define shifted_H
+
+#include "energyScalingFunction.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace energyScalingFunctions
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class shifted Declaration
+\*---------------------------------------------------------------------------*/
+
+class shifted
+:
+    public energyScalingFunction
+{
+    // Private data
+        
+        scalar e_at_rCut_;
+
+public:
+
+    //- Runtime type information
+    TypeName("shifted");
+
+
+    // Constructors
+
+        //- Construct from components
+        shifted
+        (
+            const word& name,
+            const dictionary& energyScalingFunctionProperties,
+            const pairPotential& pairPot
+        );
+
+
+    // Destructor
+
+        ~shifted()
+        {}
+
+    // Member Functions
+
+        void scaleEnergy(scalar& e, const scalar r) const;
+
+        //- Read transportProperties dictionary
+        bool read(const dictionary& energyScalingFunctionProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace energyScalingFunctions
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shiftedForce/shiftedForce.C b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shiftedForce/shiftedForce.C
new file mode 100644
index 0000000000000000000000000000000000000000..c619c3fc0676d18dc9abbf3fd238cb5cb2af2038
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shiftedForce/shiftedForce.C
@@ -0,0 +1,86 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "shiftedForce.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace energyScalingFunctions
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(shiftedForce, 0);
+
+addToRunTimeSelectionTable
+(
+    energyScalingFunction,
+    shiftedForce,
+    dictionary
+);
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+shiftedForce::shiftedForce
+(
+    const word& name,
+    const dictionary& energyScalingFunctionProperties,
+    const pairPotential& pairPot
+)
+:
+    energyScalingFunction(name, energyScalingFunctionProperties, pairPot),
+    rCut_(pairPot.rCut()),
+    e_at_rCut_(pairPot.unscaledEnergy(rCut_)),
+    de_dr_at_rCut_(pairPot.energyDerivative(rCut_, false))
+{}
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void shiftedForce::scaleEnergy(scalar& e, const scalar r) const
+{
+    e -= ( e_at_rCut_ + de_dr_at_rCut_ * (r - rCut_) );
+}
+
+bool shiftedForce::read(const dictionary& energyScalingFunctionProperties)
+{
+    energyScalingFunction::read(energyScalingFunctionProperties);
+    
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace energyScalingFunctions
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shiftedForce/shiftedForce.H b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shiftedForce/shiftedForce.H
new file mode 100644
index 0000000000000000000000000000000000000000..5900b2a869cca43bbba128b31b347efb9c4497f7
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shiftedForce/shiftedForce.H
@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    shiftedForce
+
+Description
+     
+
+SourceFiles
+    shiftedForce.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef shiftedForce_H
+#define shiftedForce_H
+
+#include "energyScalingFunction.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace energyScalingFunctions
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class shiftedForce Declaration
+\*---------------------------------------------------------------------------*/
+
+class shiftedForce
+:
+    public energyScalingFunction
+{
+    // Private data
+        
+        scalar rCut_;
+
+        scalar e_at_rCut_;
+
+        scalar de_dr_at_rCut_;
+
+    // Private Member Functions
+
+public:
+
+    //- Runtime type information
+    TypeName("shiftedForce");
+
+
+    // Constructors
+
+        //- Construct from components
+        shiftedForce
+        (
+            const word& name,
+            const dictionary& energyScalingFunctionProperties,
+            const pairPotential& pairPot
+        );
+
+
+    // Destructor
+
+        ~shiftedForce()
+        {}
+
+    // Member Functions
+        
+        void scaleEnergy(scalar& e, const scalar r) const;
+
+        //- Read transportProperties dictionary
+        bool read(const dictionary& energyScalingFunctionProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace energyScalingFunctions
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/sigmoid/sigmoid.C b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/sigmoid/sigmoid.C
new file mode 100644
index 0000000000000000000000000000000000000000..5523cff338cd0a2c0fe42f5f9f6785d11d98fa33
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/sigmoid/sigmoid.C
@@ -0,0 +1,100 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "sigmoid.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace energyScalingFunctions
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(sigmoid, 0);
+
+addToRunTimeSelectionTable
+(
+    energyScalingFunction,
+    sigmoid,
+    dictionary
+);
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+scalar sigmoid::sigmoidScale
+    (
+        const scalar r,
+        const scalar shift,
+        const scalar scale
+    ) const
+{
+    return 1.0 / (1.0 + exp( scale * (r - shift)));
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+sigmoid::sigmoid
+(
+    const word& name,
+    const dictionary& energyScalingFunctionProperties,
+    const pairPotential& pairPot
+)
+:
+    energyScalingFunction(name, energyScalingFunctionProperties, pairPot),
+    sigmoidCoeffs_(energyScalingFunctionProperties.subDict(typeName + "Coeffs")),
+    shift_(readScalar(sigmoidCoeffs_.lookup("shift"))),
+    scale_(readScalar(sigmoidCoeffs_.lookup("scale")))
+{}
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void sigmoid::scaleEnergy(scalar& e, const scalar r) const
+{
+    e *= sigmoidScale(r, shift_, scale_);
+}
+
+bool sigmoid::read(const dictionary& energyScalingFunctionProperties)
+{
+    energyScalingFunction::read(energyScalingFunctionProperties);
+
+    sigmoidCoeffs_ = energyScalingFunctionProperties.subDict(typeName + "Coeffs");
+
+    sigmoidCoeffs_.lookup("shift") >> shift_;
+    sigmoidCoeffs_.lookup("scale") >> shift_;
+    
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace energyScalingFunctions
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/sigmoid/sigmoid.H b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/sigmoid/sigmoid.H
new file mode 100644
index 0000000000000000000000000000000000000000..0ceb56b0dfa7e81cbbebfc2e6ce20678e80463f8
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/sigmoid/sigmoid.H
@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    sigmoid
+
+Description
+     
+
+SourceFiles
+    sigmoid.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sigmoid_H
+#define sigmoid_H
+
+#include "energyScalingFunction.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace energyScalingFunctions
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class sigmoid Declaration
+\*---------------------------------------------------------------------------*/
+
+class sigmoid
+:
+    public energyScalingFunction
+{
+    // Private data
+        
+        dictionary sigmoidCoeffs_;
+
+        scalar shift_;
+        scalar scale_;
+        
+    // Private Member Functions
+
+        scalar sigmoidScale
+        (
+            const scalar r,
+            const scalar shift,
+            const scalar scale
+        ) const;
+
+public:
+
+    //- Runtime type information
+    TypeName("sigmoid");
+
+
+    // Constructors
+
+        //- Construct from components
+        sigmoid
+        (
+            const word& name,
+            const dictionary& energyScalingFunctionProperties,
+            const pairPotential& pairPot
+        );
+
+
+    // Destructor
+
+        ~sigmoid()
+        {}
+
+    // Member Functions
+
+        void scaleEnergy(scalar& e, const scalar r) const;
+
+        //- Read transportProperties dictionary
+        bool read(const dictionary& energyScalingFunctionProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace energyScalingFunctions
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/basic/newPairPotential.C b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/newPairPotential.C
new file mode 100644
index 0000000000000000000000000000000000000000..f10b0aa84d46f1bd4f565a6ca79de0e7c6d0cec8
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/newPairPotential.C
@@ -0,0 +1,75 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    pairPotential
+
+\*---------------------------------------------------------------------------*/
+
+#include "pairPotential.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+autoPtr<pairPotential> pairPotential::New
+(
+    const word& name,
+    const dictionary& pairPotentialProperties
+)
+{
+    word pairPotentialTypeName(pairPotentialProperties.lookup("pairPotential"));
+
+    Info<< nl << "Selecting intermolecular pair potential "
+        << pairPotentialTypeName << " for "
+        << name << " interaction." << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(pairPotentialTypeName);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "pairPotential::New()"
+        )   << "Unknown pairPotential type "
+            << pairPotentialTypeName << endl << endl
+            << "Valid  pairPotentials are : " << endl
+            << dictionaryConstructorTablePtr_->toc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<pairPotential>
+        (cstrIter()(name, pairPotentialProperties));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.C b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.C
new file mode 100644
index 0000000000000000000000000000000000000000..86e6ae386fbeb6f0072afe5c9f784cb47fc96447
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.C
@@ -0,0 +1,238 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    pairPotential
+
+\*---------------------------------------------------------------------------*/
+
+#include "pairPotential.H"
+#include "energyScalingFunction.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(pairPotential, 0);
+defineRunTimeSelectionTable(pairPotential, dictionary);
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+void Foam::pairPotential::scaleEnergy(scalar& e, const scalar r) const
+{
+    if (!esfPtr_)
+    {
+        esfPtr_ = energyScalingFunction::New
+        (
+            name_, pairPotentialProperties_, *this
+        ).ptr();
+    }
+    
+    esfPtr_->scaleEnergy(e, r);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::pairPotential::pairPotential
+(
+    const word& name,
+    const dictionary& pairPotentialProperties
+)
+:
+    name_(name),
+    pairPotentialProperties_(pairPotentialProperties),
+    rCut_(readScalar(pairPotentialProperties_.lookup("rCut"))),
+    rCutSqr_(rCut_*rCut_),
+    rMin_(readScalar(pairPotentialProperties_.lookup("rMin"))),
+    dr_(readScalar(pairPotentialProperties_.lookup("dr"))),
+    forceLookup_(0),
+    energyLookup_(0),
+    esfPtr_(NULL),
+    writeTables_(Switch(pairPotentialProperties_.lookup("writeTables")))
+{}
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::pairPotential::setLookupTables()
+{
+    label N = label((rCut_ - rMin_)/dr_) + 1;
+
+    forceLookup_.setSize(N);
+
+    energyLookup_.setSize(N);
+
+    forAll(forceLookup_, k)
+    {
+        energyLookup_[k] = scaledEnergy(k*dr_ + rMin_);
+
+        forceLookup_[k] = -energyDerivative( (k * dr_ + rMin_), true);
+    }
+}
+
+
+Foam::scalar Foam::pairPotential::forceLookup(const scalar r) const
+{
+    scalar k_rIJ = (r - rMin_)/dr_;
+
+    label k(k_rIJ);
+
+    if (k < 0)
+    {
+        FatalErrorIn("pairPotential.C") << nl
+        << "r less than rMin" << nl
+        << abort(FatalError);
+    }
+
+    scalar f =
+        (k_rIJ - k)*forceLookup_[k+1]
+      + (k + 1 - k_rIJ)*forceLookup_[k];
+
+    return f;
+}
+
+
+Foam::List< Foam::Pair< Foam::scalar > >
+Foam::pairPotential::forceTable() const
+{
+    List< Pair<scalar> > forceTab(forceLookup_.size());
+
+    forAll(forceLookup_,k)
+    {
+        forceTab[k].first() = rMin_ + k*dr_;
+
+        forceTab[k].second() = forceLookup_[k];
+    }
+
+    return forceTab;
+}
+
+
+Foam::scalar Foam::pairPotential::energyLookup(const scalar r) const
+{
+    scalar k_rIJ = (r - rMin_)/dr_;
+
+    label k(k_rIJ);
+
+    if (k < 0)
+    {
+        FatalErrorIn("pairPotential.C") << nl
+        << "r less than rMin" << nl
+        << abort(FatalError);
+    }
+
+    scalar e =
+        (k_rIJ - k)*energyLookup_[k+1]
+      + (k + 1 - k_rIJ)*energyLookup_[k];
+
+    return e;
+}
+
+
+Foam::List< Foam::Pair< Foam::scalar > >
+    Foam::pairPotential::energyTable() const
+{
+    List< Pair<scalar> > energyTab(energyLookup_.size());
+
+    forAll(energyLookup_,k)
+    {
+        energyTab[k].first() = rMin_ + k*dr_;
+
+        energyTab[k].second() = energyLookup_[k];
+    }
+
+    return energyTab;
+}
+
+
+Foam::scalar Foam::pairPotential::scaledEnergy(const scalar r) const
+{
+    scalar e = unscaledEnergy(r);
+    
+    scaleEnergy(e, r);
+
+    return e;
+}
+
+
+Foam::scalar Foam::pairPotential::energyDerivative
+(
+    const scalar r,
+    const bool scaledEnergyDerivative
+) const
+{
+    // Local quadratic fit to energy: E = a0 + a1*r + a2*r^2
+    // Differentiate to give f = -dE/dr = -a1 - 2*a2*r
+
+    scalar ra = r - dr_;
+    scalar rf = r;
+    scalar rb = r + dr_;
+    
+    scalar Ea, Ef, Eb;
+
+    if (scaledEnergyDerivative)
+    {
+        Ea = scaledEnergy(ra);
+        Ef = scaledEnergy(rf);
+        Eb = scaledEnergy(rb);
+    }
+    else
+    {
+        Ea = unscaledEnergy(ra);
+        Ef = unscaledEnergy(rf);
+        Eb = unscaledEnergy(rb);
+    }
+    
+    scalar denominator = (ra - rf)*(ra - rb)*(rf - rb);
+    
+    scalar a1 =
+    (
+        rb*rb*(Ea - Ef) + ra*ra*(Ef - Eb) + rf*rf*(Eb - Ea)
+    ) / denominator;
+
+    scalar a2 = 
+    (
+        rb*(Ef - Ea) + rf*(Ea - Eb) + ra*(Eb - Ef)
+    ) / denominator;
+
+    return a1 + 2.0 * a2 * r;
+}
+
+
+bool Foam::pairPotential::read(const dictionary& pairPotentialProperties)
+{
+    pairPotentialProperties_ = pairPotentialProperties;
+
+    return true;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecule/potentials/pairPotential/basic/pairPotential.H b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.H
old mode 100755
new mode 100644
similarity index 50%
rename from src/lagrangian/molecule/potentials/pairPotential/basic/pairPotential.H
rename to src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.H
index 7427f193315dfbdea463534da0ff8fcf295572e5..a303750ba55b058b52e12858a414666cf093b066
--- a/src/lagrangian/molecule/potentials/pairPotential/basic/pairPotential.H
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2005 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,161 +23,160 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::pairPotential
+    pairPotential
 
 Description
 
-    At the moment this is hard coded to be a shifted force
-    "maitlandSmith" potential.
-    In future use templated classes virtual functions,
-    function pointers and all kinds of good stuff to make
-    this a generic pair force,
-
 SourceFiles
-    pairPotentialI.H
     pairPotential.C
+    newPairPotential.C
 
 \*---------------------------------------------------------------------------*/
 
 #ifndef pairPotential_H
 #define pairPotential_H
 
-#include "vector.H"
-#include "dictionary.H"
+#include "IOdictionary.H"
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "autoPtr.H"
 #include "List.H"
 #include "Pair.H"
+#include "Switch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+// Forward declaration of classes
+class energyScalingFunction;
+
 /*---------------------------------------------------------------------------*\
-                           Class pairPotential Declaration
+                Class pairPotential Declaration
 \*---------------------------------------------------------------------------*/
 
 class pairPotential
 {
-    // Private data
-
-        List<scalar> forceLookup_;
-
-        List<scalar> energyLookup_;
-
-        scalar m_;
 
-        scalar gamma_;
+protected:
 
-        scalar rm_;
+    // Protected data
 
-        scalar epsilon_;
+        word name_;
+        dictionary pairPotentialProperties_;
 
         scalar rCut_;
-
         scalar rCutSqr_;
 
-        scalar u_at_rCut_;
-
-        scalar du_by_dr_at_rCut_;
-
         scalar rMin_;
-
         scalar dr_;
 
+        List<scalar> forceLookup_;
+        List<scalar> energyLookup_;
 
-    // Private Member Functions
+        mutable energyScalingFunction* esfPtr_;
 
-        void setLookupTables();
+        bool writeTables_;
 
-        void setConstants();
+    // Private Member Functions
+
+        void scaleEnergy(scalar& e, const scalar r) const;
+        
+        //- Disallow copy construct
+        pairPotential(const pairPotential&);
 
+        //- Disallow default bitwise assignment
+        void operator=(const pairPotential&);
 
 public:
 
-    // Constructors
+    //- Runtime type information
+    TypeName("pairPotential");
 
-        //- Construct null
-        pairPotential();
 
-        //- Construct from dictionary
-        pairPotential(const dictionary& pPDict);
+    // Declare run-time constructor selection table
 
-        //- Construct from components
-        pairPotential
+        declareRunTimeSelectionTable
         (
-            const scalar m,
-            const scalar gamma,
-            const scalar rm,
-            const scalar epsilon,
-            const scalar rCut,
-            const scalar rMin,
-            const scalar dr
+            autoPtr,
+            pairPotential,
+            dictionary,
+            (
+                const word& name,
+                const dictionary& pairPotentialProperties
+            ),
+            (name, pairPotentialProperties)
         );
 
 
-    // Destructor
-
-        virtual ~pairPotential();
-
-
-    // Member Functions
+    // Selectors
 
-        // Access
-
-            inline const scalar m() const;
+        //- Return a reference to the selected viscosity model
+        static autoPtr<pairPotential> New
+        (
+            const word& name,
+            const dictionary& pairPotentialProperties
+        );
 
-            inline const scalar gamma() const ;
 
-            inline const scalar rm() const;
+    // Constructors
 
-            inline const scalar epsilon() const;
+        //- Construct from components
+        pairPotential
+        (
+            const word& name,
+            const dictionary& pairPotentialProperties
+        );
 
-            inline const scalar rCut() const;
 
-            inline const scalar rCutSqr() const;
+    // Destructor
 
-            inline const scalar rMin() const;
+        virtual ~pairPotential()
+        {}
 
-            inline const scalar dr() const;
 
+    // Member Functions
 
-        // Write
+        void setLookupTables();
 
-            virtual void write(Ostream&) const;
+        inline scalar rMin() const;
 
-            scalar n(const scalar rIJMag) const;
+        inline scalar dr() const;
 
-            scalar force(const scalar rIJMag) const;
+        inline scalar rCut() const;
 
-            scalar forceLookup(const scalar rIJMag) const;
+        inline scalar rCutSqr() const;
 
-            List<Pair<scalar> > forceTable() const;
+        scalar energyLookup (const scalar r) const;
 
-            scalar energy(const scalar rIJMag) const;
+        scalar forceLookup (const scalar r) const;
 
-            scalar energyLookup(const scalar rIJMag) const;
+        List< Pair <scalar> > energyTable() const;
 
-            List<Pair<scalar> > energyTable() const;
+        List< Pair <scalar> > forceTable() const;
 
+        inline bool writeTables() const;
 
-    // Friend Operators
+        virtual scalar unscaledEnergy (const scalar r) const = 0;
 
-        inline friend bool operator==
-        (
-            const pairPotential& a,
-            const pairPotential& b
-        );
+        scalar scaledEnergy(const scalar r) const;
 
-        inline friend bool operator!=
+        scalar energyDerivative
         (
-            const pairPotential& a,
-            const pairPotential& b
-        );
+            const scalar r,
+            const bool scaledEnergyDerivative = true
+        ) const;
 
+        const dictionary& pairPotentialProperties() const
+        {
+            return pairPotentialProperties_;
+        }
 
-    // IOstream Operators
+        bool writeEnergyAndForceTables(Ostream& os) const;
 
-        friend Ostream& operator<<(Ostream&, const pairPotential&);
+        //- Read pairPotential dictionary
+        virtual bool read(const dictionary& pairPotentialProperties) = 0;
 };
 
 
diff --git a/src/lagrangian/molecule/potentials/pairPotential/derived/lennardJonesI.H b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotentialI.H
similarity index 76%
rename from src/lagrangian/molecule/potentials/pairPotential/derived/lennardJonesI.H
rename to src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotentialI.H
index 4e3a1a2f27c63584a5e6e0cdd8e512a4660b94d0..67dbbfd253e70a7aab8b718714c1e9f3a565027e 100755
--- a/src/lagrangian/molecule/potentials/pairPotential/derived/lennardJonesI.H
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotentialI.H
@@ -24,46 +24,35 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline const word& pairPotential::pairPotentialName() const
-{
-    return pairPotentialName_;
-}
-
-
-inline const word& pairPotential::pairPotentialType() const
+inline Foam::scalar Foam::pairPotential::rMin() const
 {
-    return pairPotentialType_;
+    return rMin_;
 }
 
 
-inline scalar pairPotential::sigma() const
+inline Foam::scalar Foam::pairPotential::dr() const
 {
-    return sigma_;
+    return dr_;
 }
 
 
-inline scalar pairPotential::epsilon() const
+inline Foam::scalar Foam::pairPotential::rCut() const
 {
-    return epsilon_;
+    return rCut_;
 }
 
 
-inline scalar pairPotential::rCut() const
+inline Foam::scalar Foam::pairPotential::rCutSqr() const
 {
-    return rCut_;
+    return rCutSqr_;
 }
 
 
-inline scalar pairPotential::rCutSqr() const
+inline bool Foam::pairPotential::writeTables() const
 {
-    return rCutSqr_;
+    return writeTables_;
 }
 
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotentialIO.C b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotentialIO.C
new file mode 100644
index 0000000000000000000000000000000000000000..79ab55e737b4ff454e7dc9ab08b669c635910646
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotentialIO.C
@@ -0,0 +1,53 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "pairPotential.H"
+#include "IOstreams.H"
+
+bool Foam::pairPotential::writeEnergyAndForceTables(Ostream& os) const
+{
+    Info<< "Writing energy and force tables to file for potential "
+        << name_ << endl;
+
+    List< Pair <scalar> > eTab(energyTable());
+
+    List< Pair <scalar> > fTab(forceTable());
+
+    forAll(eTab, e)
+    {
+        os<< eTab[e].first()
+          << token::SPACE
+          << eTab[e].second()
+          << token::SPACE
+          << fTab[e].second()
+          << nl;
+    }
+
+    return os.good();
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/azizChen/azizChen.C b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/azizChen/azizChen.C
new file mode 100644
index 0000000000000000000000000000000000000000..b9cbbb1e5488c3a7cc083d6b6dfaf81b3b21f786
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/azizChen/azizChen.C
@@ -0,0 +1,125 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "azizChen.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace pairPotentials
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(azizChen, 0);
+
+addToRunTimeSelectionTable
+(
+    pairPotential,
+    azizChen,
+    dictionary
+);
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+azizChen::azizChen
+(
+    const word& name,
+    const dictionary& azizChen
+)
+:
+    pairPotential(name, azizChen),
+    azizChenCoeffs_(azizChen.subDict(typeName + "Coeffs")),
+    epsilon_(readScalar(azizChenCoeffs_.lookup("epsilon"))),
+    rm_(readScalar(azizChenCoeffs_.lookup("rm"))),
+    A_(readScalar(azizChenCoeffs_.lookup("A"))),
+    alpha_(readScalar(azizChenCoeffs_.lookup("alpha"))),
+    C6_(readScalar(azizChenCoeffs_.lookup("C6"))),
+    C8_(readScalar(azizChenCoeffs_.lookup("C8"))),
+    C10_(readScalar(azizChenCoeffs_.lookup("C10"))),
+    D_(readScalar(azizChenCoeffs_.lookup("D"))),
+    gamma_(readScalar(azizChenCoeffs_.lookup("gamma")))
+{
+    setLookupTables();
+}
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+scalar azizChen::unscaledEnergy(const scalar r) const
+{
+    scalar x = r/rm_;
+
+    scalar F = 1.0;
+
+    if (x < D_)
+    {
+        F = exp(-pow(((D_ / x) - 1.0),2));
+    }
+
+    return epsilon_ *
+    (
+        A_ * Foam::pow(x, gamma_) * exp(-alpha_ * x)
+      - (
+            (C6_/ Foam::pow(x, 6))
+          + (C8_/ Foam::pow(x, 8))
+          + (C10_/ Foam::pow(x, 10))
+        )
+      * F
+    );
+}
+
+
+bool azizChen::read(const dictionary& azizChen)
+{
+    pairPotential::read(azizChen);
+
+    azizChenCoeffs_ = azizChen.subDict(typeName + "Coeffs");
+
+    azizChenCoeffs_.lookup("epsilon") >> epsilon_;
+    azizChenCoeffs_.lookup("rm") >> rm_;
+    azizChenCoeffs_.lookup("A") >> A_;
+    azizChenCoeffs_.lookup("alpha") >> alpha_;
+    azizChenCoeffs_.lookup("C6") >> C6_;
+    azizChenCoeffs_.lookup("C8") >> C8_;
+    azizChenCoeffs_.lookup("C10") >> C10_;
+    azizChenCoeffs_.lookup("D") >> D_;
+    azizChenCoeffs_.lookup("gamma") >> gamma_;
+    
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace pairPotentials
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/azizChen/azizChen.H b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/azizChen/azizChen.H
new file mode 100644
index 0000000000000000000000000000000000000000..d0513b277fced26acd6a12367d4ec8d89bf917dd
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/azizChen/azizChen.H
@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    azizChen
+
+Description
+
+    From:
+
+    @article{MA_Aziz_Chen,
+    author = {R. A. Aziz and H. H. Chen},
+    collaboration = {},
+    title = {An accurate intermolecular potential for argon},
+    publisher = {AIP},
+    year = {1977},
+    journal = {The Journal of Chemical Physics},
+    volume = {67},
+    number = {12},
+    pages = {5719-5726},
+    url = {http://link.aip.org/link/?JCP/67/5719/1},
+    doi = {10.1063/1.434827}
+    }
+
+SourceFiles
+    azizChen.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef azizChen_H
+#define azizChen_H
+
+#include "pairPotential.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace pairPotentials
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class azizChen Declaration
+\*---------------------------------------------------------------------------*/
+
+class azizChen
+:
+    public pairPotential
+{
+    // Private data
+
+        dictionary azizChenCoeffs_;
+
+        scalar epsilon_;
+        scalar rm_;
+        scalar A_;
+        scalar alpha_;
+        scalar C6_;
+        scalar C8_;
+        scalar C10_;
+        scalar D_;
+        scalar gamma_;
+
+public:
+
+    //- Runtime type information
+    TypeName("azizChen");
+
+
+    // Constructors
+
+        //- Construct from components
+        azizChen
+        (
+            const word& name,
+            const dictionary& pairPotentialProperties
+        );
+
+
+    // Destructor
+
+        ~azizChen()
+        {}
+
+
+    // Member Functions
+
+        scalar unscaledEnergy(const scalar r) const;
+
+        //- Read transportProperties dictionary
+        bool read(const dictionary& pairPotentialProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace pairPotentials
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/lennardJones/lennardJones.C b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/lennardJones/lennardJones.C
new file mode 100644
index 0000000000000000000000000000000000000000..58968a6d290f423bf41db29cf332cb5a69fa464d
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/lennardJones/lennardJones.C
@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "lennardJones.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace pairPotentials
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(lennardJones, 0);
+
+addToRunTimeSelectionTable
+(
+    pairPotential,
+    lennardJones,
+    dictionary
+);
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+lennardJones::lennardJones
+(
+    const word& name,
+    const dictionary& pairPotentialProperties
+)
+:
+    pairPotential(name, pairPotentialProperties),
+    lennardJonesCoeffs_(pairPotentialProperties.subDict(typeName + "Coeffs")),
+    sigma_(readScalar(lennardJonesCoeffs_.lookup("sigma"))),
+    epsilon_(readScalar(lennardJonesCoeffs_.lookup("epsilon")))
+{
+    setLookupTables();
+}
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+scalar lennardJones::unscaledEnergy(const scalar r) const
+{
+    // (rIJ/sigma)^-2
+    scalar ir2 = (sigma_/r)*(sigma_/r);
+
+    // (rIJ/sigma)^-6
+    scalar ir6 = ir2*ir2*ir2;
+
+    return 4.0 * epsilon_ * (ir6 * (ir6 - 1.0));
+}
+
+
+bool lennardJones::read(const dictionary& pairPotentialProperties)
+{
+    pairPotential::read(pairPotentialProperties);
+
+    lennardJonesCoeffs_ = pairPotentialProperties.subDict(typeName + "Coeffs");
+
+    lennardJonesCoeffs_.lookup("sigma") >> sigma_;
+    lennardJonesCoeffs_.lookup("epsilon") >> epsilon_;
+    
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace pairPotentials
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/lennardJones/lennardJones.H b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/lennardJones/lennardJones.H
new file mode 100644
index 0000000000000000000000000000000000000000..f879549ef85d0d98cfcb345130a3d195d9af2972
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/lennardJones/lennardJones.H
@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    lennardJones
+
+Description
+     
+
+SourceFiles
+    lennardJones.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef lennardJones_H
+#define lennardJones_H
+
+#include "pairPotential.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace pairPotentials
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class lennardJones Declaration
+\*---------------------------------------------------------------------------*/
+
+class lennardJones
+:
+    public pairPotential
+{
+    // Private data
+
+        dictionary lennardJonesCoeffs_;
+
+        scalar sigma_;
+        scalar epsilon_;
+
+public:
+
+    //- Runtime type information
+    TypeName("lennardJones");
+
+
+    // Constructors
+
+        //- Construct from components
+        lennardJones
+        (
+            const word& name,
+            const dictionary& pairPotentialProperties
+        );
+
+
+    // Destructor
+
+        ~lennardJones()
+        {}
+
+
+    // Member Functions
+
+        scalar unscaledEnergy(const scalar r) const;
+
+        //- Read transportProperties dictionary
+        bool read(const dictionary& pairPotentialProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace pairPotentials
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/maitlandSmith/maitlandSmith.C b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/maitlandSmith/maitlandSmith.C
new file mode 100644
index 0000000000000000000000000000000000000000..4915f2a1fd03dc36a4294a5e665e0baf4d05f1ad
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/maitlandSmith/maitlandSmith.C
@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "maitlandSmith.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace pairPotentials
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(maitlandSmith, 0);
+
+addToRunTimeSelectionTable
+(
+    pairPotential,
+    maitlandSmith,
+    dictionary
+);
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+maitlandSmith::maitlandSmith
+(
+    const word& name,
+    const dictionary& maitlandSmith
+)
+:
+    pairPotential(name, maitlandSmith),
+    maitlandSmithCoeffs_(maitlandSmith.subDict(typeName + "Coeffs")),
+    m_(readScalar(maitlandSmithCoeffs_.lookup("m"))),
+    gamma_(readScalar(maitlandSmithCoeffs_.lookup("gamma"))),
+    rm_(readScalar(maitlandSmithCoeffs_.lookup("rm"))),
+    epsilon_(readScalar(maitlandSmithCoeffs_.lookup("epsilon")))
+{
+    setLookupTables();
+}
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+scalar maitlandSmith::unscaledEnergy(const scalar r) const
+{
+    scalar nr = (m_ + gamma_*(r/rm_ - 1.0));
+
+    return epsilon_
+       *(
+            (6.0 / (nr - 6.0))*Foam::pow(r/rm_, -nr)
+          - (nr / (nr - 6.0))*Foam::pow(r/rm_, -6)
+        );
+}
+
+
+bool maitlandSmith::read(const dictionary& maitlandSmith)
+{
+    pairPotential::read(maitlandSmith);
+
+    maitlandSmithCoeffs_ = maitlandSmith.subDict(typeName + "Coeffs");
+
+    maitlandSmithCoeffs_.lookup("m") >> m_;
+    maitlandSmithCoeffs_.lookup("gamma") >> gamma_;
+    maitlandSmithCoeffs_.lookup("rm") >> rm_;
+    maitlandSmithCoeffs_.lookup("epsilon") >> epsilon_;
+    
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace pairPotentials
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/maitlandSmith/maitlandSmith.H b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/maitlandSmith/maitlandSmith.H
new file mode 100644
index 0000000000000000000000000000000000000000..64747a369a7045fe71253d3a5510643ddcc9d098
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/maitlandSmith/maitlandSmith.H
@@ -0,0 +1,127 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    maitlandSmith
+
+Description
+
+    From:
+
+    @ARTICLE{MA_Maitland_Smith,
+    author = {{Maitland}, G.~C. and {Smith}, E.~B.},
+    title = {A simplified representation of intermolecular potential energy},
+    journal = {Chemical Physics Letters},
+    year = 1973,
+    month = oct,
+    volume = 22,
+    pages = {443-446},
+    adsurl = {http://adsabs.harvard.edu/abs/1973CPL....22..443M},
+    adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+    }
+
+    Parameters for other monoatomics from:
+    @BOOK{MD_Maitland_Rigby_Smith_Wakeham,
+    AUTHOR =       {Geoffrey C. Maitland and Maurice Rigby and E. Brian Smith and William A. Wakeham},
+    TITLE =        {Intermolecular Forces: Their Origin and Determination},
+    PUBLISHER =    {Oxford University Press},
+    YEAR =         {1981}
+    }
+
+SourceFiles
+    maitlandSmith.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef maitlandSmith_H
+#define maitlandSmith_H
+
+#include "pairPotential.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace pairPotentials
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class maitlandSmith Declaration
+\*---------------------------------------------------------------------------*/
+
+class maitlandSmith
+:
+    public pairPotential
+{
+    // Private data
+
+        dictionary maitlandSmithCoeffs_;
+
+        scalar m_;
+        scalar gamma_;
+        scalar rm_;
+        scalar epsilon_;
+
+public:
+
+    //- Runtime type information
+    TypeName("maitlandSmith");
+
+
+    // Constructors
+
+        //- Construct from components
+        maitlandSmith
+        (
+            const word& name,
+            const dictionary& pairPotentialProperties
+        );
+
+
+    // Destructor
+
+        ~maitlandSmith()
+        {}
+
+
+    // Member Functions
+
+        scalar unscaledEnergy(const scalar r) const;
+
+        //- Read transportProperties dictionary
+        bool read(const dictionary& pairPotentialProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace pairPotentials
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.C b/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.C
new file mode 100755
index 0000000000000000000000000000000000000000..44207d13ba094b89a6cbd833d8a41c934c244047
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.C
@@ -0,0 +1,284 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2005 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "pairPotentialList.H"
+#include "OFstream.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::pairPotentialList::readPairPotentialDict
+(
+    const dictionary& pairPotentialDict,
+    const polyMesh& mesh
+)
+{
+    Info<< nl << "Building pair potentials." << endl;   
+
+    rCutMax_ = 0.0;
+
+    for (label a = 0; a < nIds(); ++a)
+    {
+        word idA = idList_[a];
+    
+        for (label b = a; b < nIds(); ++b)
+        {
+            word idB = idList_[b];
+
+            word pairPotentialName;
+
+            if (a==b)
+            {
+                if(pairPotentialDict.found(idA+"-"+idB))
+                {
+                    pairPotentialName = idA+"-"+idB;
+                }
+                else
+                {
+                    FatalErrorIn("pairPotentialList::buildPotentials") << nl
+                            << "Pair pairPotential specification subDict "
+                            << idA+"-"+idB << " not found"
+                            << abort(FatalError);
+                }
+            }
+            else
+            {
+                if(pairPotentialDict.found(idA+"-"+idB))
+                {
+                    pairPotentialName = idA+"-"+idB;
+                }
+    
+                else if(pairPotentialDict.found(idB+"-"+idA))
+                {
+                    pairPotentialName = idB+"-"+idA;
+                }
+    
+                else
+                {
+                    FatalErrorIn("pairPotentialList::buildPotentials") << nl
+                            << "Pair pairPotential specification subDict "
+                            << idA+"-"+idB << " or " << idB+"-"+idA << " not found"
+                            << abort(FatalError);
+                }
+    
+                if
+                (
+                    pairPotentialDict.found(idA+"-"+idB)
+                 && pairPotentialDict.found(idB+"-"+idA)
+                )
+                {
+                    FatalErrorIn("pairPotentialList::buildPotentials") << nl
+                            << "Pair pairPotential specification subDict "
+                            << idA+"-"+idB << " and "
+                            << idB+"-"+idA << " found - multiple definition"
+                            << abort(FatalError);
+                }
+            }
+
+            (*this).set
+            (
+                pairPotentialIndex(a,b),
+                pairPotential::New
+                (
+                    pairPotentialName,
+                    pairPotentialDict.subDict(pairPotentialName)
+                )
+            );
+
+            if ((*this)[pairPotentialIndex(a,b)].rCut() > rCutMax_)
+            {
+                rCutMax_ = (*this)[pairPotentialIndex(a,b)].rCut();
+            }
+
+//             Info<< pairPotentialName << ": "
+//                 << (*this)[pairPotentialIndex(a,b)].pairPotentialProperties()
+//                 << endl;
+
+            if ((*this)[pairPotentialIndex(a,b)].writeTables())
+            {
+                OFstream ppTabFile(mesh.time().path()/pairPotentialName);
+            
+                if 
+                (
+                    !(*this)[pairPotentialIndex(a,b)].writeEnergyAndForceTables
+                    (
+                        ppTabFile
+                    )
+                )
+                {
+                    FatalErrorIn("pairPotentialList::readPairPotentialDict")
+                        << "Failed writing to "
+                        << ppTabFile.name()
+                        << abort(FatalError);
+                }
+            }
+        }
+    }
+
+    rCutMaxSqr_ = rCutMax_ * rCutMax_;
+
+    Info << nl << "rCutMax = " << rCutMax_ << endl;
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::pairPotentialList::pairPotentialList()
+:
+    PtrList<pairPotential>(),
+    idList_()
+{}
+
+Foam::pairPotentialList::pairPotentialList
+(
+    const dictionary& idListDict,
+    const dictionary& pairPotentialDict,
+    const polyMesh& mesh
+)
+:
+    PtrList<pairPotential>(),
+    idList_()
+{
+    buildPotentials(idListDict, pairPotentialDict, mesh);
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::pairPotentialList::~pairPotentialList()
+{}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::pairPotentialList::buildPotentials
+(
+    const dictionary& idListDict,
+    const dictionary& pairPotentialDict,
+    const polyMesh& mesh
+)
+{
+    idList_ = List<word>(idListDict.lookup("idList"));
+
+    setSize(((idList_.size() * (idList_.size() + 1))/2));
+
+    readPairPotentialDict(pairPotentialDict, mesh);
+}
+
+
+const Foam::pairPotential& Foam::pairPotentialList::pairPotentialFunction
+(
+    const label a,
+    const label b
+) const
+{
+    return (*this)[pairPotentialIndex (a,b)];
+}
+
+
+bool Foam::pairPotentialList::rCutSqr
+(
+    const label a,
+    const label b,
+    const scalar rIJMagSqr
+) const
+{
+    if(rIJMagSqr <= rCutSqr (a,b))
+    {
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+Foam::scalar Foam::pairPotentialList::rMin
+(
+    const label a,
+    const label b
+) const
+{
+    return (*this)[pairPotentialIndex (a,b)].rMin();
+}
+
+
+Foam::scalar Foam::pairPotentialList::dr
+(
+    const label a,
+    const label b
+) const
+{
+    return (*this)[pairPotentialIndex (a,b)].dr();
+}
+
+
+Foam::scalar Foam::pairPotentialList::rCutSqr
+(
+    const label a,
+    const label b
+) const
+{
+    return (*this)[pairPotentialIndex (a,b)].rCutSqr();
+}
+
+
+Foam::scalar Foam::pairPotentialList::rCut
+(
+    const label a,
+    const label b
+) const
+{
+    return (*this)[pairPotentialIndex (a,b)].rCut();
+}
+
+
+Foam::scalar Foam::pairPotentialList::force
+(
+    const label a,
+    const label b,
+    const scalar rIJMag
+) const
+{
+    scalar f = (*this)[pairPotentialIndex (a,b)].forceLookup(rIJMag);
+
+    return f;
+}
+
+
+Foam::scalar Foam::pairPotentialList::energy
+(
+    const label a,
+    const label b,
+    const scalar rIJMag
+) const
+{
+    scalar e = (*this)[pairPotentialIndex (a,b)].energyLookup(rIJMag);
+
+    return e;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecule/potentials/pairPotential/basic/pairPotentialList.H b/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.H
similarity index 76%
rename from src/lagrangian/molecule/potentials/pairPotential/basic/pairPotentialList.H
rename to src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.H
index 26cbc0ca642383fa322c82b847b4239bbcf785a2..92526cf4792c40cd8a71062a59300b722bd1b74a 100755
--- a/src/lagrangian/molecule/potentials/pairPotential/basic/pairPotentialList.H
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.H
@@ -36,7 +36,10 @@ SourceFiles
 #ifndef pairPotentialList_H
 #define pairPotentialList_H
 
+#include "PtrList.H"
+#include "word.H"
 #include "pairPotential.H"
+#include "polyMesh.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,12 +52,15 @@ namespace Foam
 
 class pairPotentialList
 :
-    public List<pairPotential>
+    public PtrList<pairPotential>
 {
     // Private data
 
-        label nIds_;
+        List<word> idList_;
+        
+        scalar rCutMax_;
 
+        scalar rCutMaxSqr_;
 
     // Private Member Functions
 
@@ -64,53 +70,54 @@ class pairPotentialList
             const label b
         ) const;
 
+        void readPairPotentialDict
+        (
+            const dictionary& pairPotentialDict,
+            const polyMesh& mesh
+        );
+
         //- Disallow default bitwise assignment
         void operator=(const pairPotentialList&);
 
         //- Disallow default bitwise copy construct
         pairPotentialList(const pairPotentialList&);
 
-
 public:
 
     // Constructors
 
-        //- Construct null
         pairPotentialList();
 
-        //- Construct from number of Ids
+        //- Construct from idList and potental dictionaries
         pairPotentialList
         (
-            const label nIds
+            const dictionary& idListDict,
+            const dictionary& pairPotentialDict,
+            const polyMesh& mesh
         );
 
-        //- Construct from components
-        pairPotentialList
-        (
-            const List<pairPotential>& pairPotentials,
-            const label nIds
-        );
-
-
     // Destructor
 
         ~pairPotentialList();
 
-
     // Member Functions
+        
+        void buildPotentials
+        (
+            const dictionary& idListDict,
+            const dictionary& pairPotentialDict,
+            const polyMesh& mesh
+        );
 
         // Access
 
+            inline const List<word>& idList() const;
+
             inline label nIds() const;
 
-            void setSizeAndNIds (const label);
+            inline scalar rCutMax() const;
 
-            void addPotential
-            (
-                const label a,
-                const label b,
-                const pairPotential& pot
-            );
+            inline scalar rCutMaxSqr() const;
 
             const pairPotential& pairPotentialFunction
             (
@@ -126,17 +133,13 @@ public:
                 const scalar rIJMagSqr
             ) const;
 
-            const scalar rCutSqr
-            (
-                const label a,
-                const label b
-            ) const;
+            scalar rMin (const label a, const label b) const;
 
-            const scalar rCut
-            (
-                const label a,
-                const label b
-            ) const;
+            scalar dr (const label a, const label b) const;
+
+            scalar rCutSqr (const label a, const label b) const;
+
+            scalar rCut (const label a, const label b) const;
 
             scalar force
             (
diff --git a/src/lagrangian/molecule/potentials/pairPotential/basic/pairPotentialI.H b/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialListI.H
similarity index 54%
rename from src/lagrangian/molecule/potentials/pairPotential/basic/pairPotentialI.H
rename to src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialListI.H
index 71424d740316e928c57603a9ddbbbf1aa3a718fd..ed02ca3b7cffb24febaf4e9c7cd94d99b5256347 100755
--- a/src/lagrangian/molecule/potentials/pairPotential/basic/pairPotentialI.H
+++ b/src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialListI.H
@@ -24,79 +24,61 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-inline const Foam::scalar Foam::pairPotential::m() const
-{
-    return m_;
-}
-
-
-inline const Foam::scalar Foam::pairPotential::gamma() const
-{
-    return gamma_;
-}
-
-
-inline const Foam::scalar Foam::pairPotential::rm() const
-{
-    return rm_;
-}
-
-
-inline const Foam::scalar Foam::pairPotential::epsilon() const
-{
-    return epsilon_;
-}
-
-
-inline const Foam::scalar Foam::pairPotential::rCut() const
-{
-    return rCut_;
-}
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-
-inline const Foam::scalar Foam::pairPotential::rCutSqr() const
+inline Foam::label Foam::pairPotentialList::pairPotentialIndex
+(
+    const label a,
+    const label b
+) const
 {
-    return rCutSqr_;
+    label index;
+
+    if (a < b)
+    {
+        index = a*(2 * nIds() - a - 1) / 2 + b;
+    }
+
+    else
+    {
+        index = b*(2 * nIds() - b - 1) / 2 + a;
+    }
+
+    if (index > size() - 1)
+    {
+        FatalErrorIn("Foam::pairPotentialList::pairPotentialIndex ")
+            << "Attempting to access a pairPotential with too high an index."
+            << nl << "a = " << a << ", b = " << b << ", index = " << index
+            << nl << "max index = " << size() - 1
+            << nl << abort(FatalError);
+    }
+
+    return index;
 }
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline const Foam::scalar Foam::pairPotential::rMin() const
+inline Foam::label Foam::pairPotentialList::nIds() const
 {
-    return rMin_;
+    return idList_.size();
 }
 
 
-inline const Foam::scalar Foam::pairPotential::dr() const
+inline const Foam::List<Foam::word>& Foam::pairPotentialList::idList() const
 {
-    return dr_;
+    return idList_;
 }
 
 
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
-inline bool Foam::operator==
-(
-    const pairPotential& a,
-    const pairPotential& b
-)
+inline Foam::scalar Foam::pairPotentialList::rCutMax() const
 {
-    return
-        (a.m() == b.m())
-     && (a.gamma() == b.gamma())
-     && (a.rm() == b.rm())
-     && (a.epsilon() == b.epsilon())
-     && (a.rCut() == b.rCut());
+    return rCutMax_;
 }
 
 
-inline bool Foam::operator!=
-(
-    const pairPotential& a,
-    const pairPotential& b
-)
+inline Foam::scalar Foam::pairPotentialList::rCutMaxSqr() const
 {
-    return !(a == b);
+    return rCutMaxSqr_;
 }
 
-
 // ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/newTetherPotential.C b/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/newTetherPotential.C
new file mode 100644
index 0000000000000000000000000000000000000000..d9aaf566fca94cbf8ef45afce4d318c28d382cca
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/newTetherPotential.C
@@ -0,0 +1,75 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    tetherPotential
+
+\*---------------------------------------------------------------------------*/
+
+#include "tetherPotential.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+autoPtr<tetherPotential> tetherPotential::New
+(
+    const word& name,
+    const dictionary& tetherPotentialProperties
+)
+{
+    word tetherPotentialTypeName(tetherPotentialProperties.lookup("tetherPotential"));
+
+    Info<< nl << "Selecting tether potential "
+        << tetherPotentialTypeName << " for "
+        << name << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(tetherPotentialTypeName);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "tetherPotential::New()"
+        )   << "Unknown tetherPotential type "
+            << tetherPotentialTypeName << endl << endl
+            << "Valid  tetherPotentials are : " << endl
+            << dictionaryConstructorTablePtr_->toc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<tetherPotential>
+        (cstrIter()(name, tetherPotentialProperties));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.C b/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.C
new file mode 100644
index 0000000000000000000000000000000000000000..de06d6607c3d83ec4a37616d283dca110f7bc516
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.C
@@ -0,0 +1,71 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    tetherPotential
+
+\*---------------------------------------------------------------------------*/
+
+#include "tetherPotential.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(tetherPotential, 0);
+defineRunTimeSelectionTable(tetherPotential, dictionary);
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::tetherPotential::tetherPotential
+(
+    const word& name,
+    const dictionary& tetherPotentialProperties
+)
+:
+    name_(name),
+    tetherPotentialProperties_(tetherPotentialProperties)
+{}
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+
+bool Foam::tetherPotential::read(const dictionary& tetherPotentialProperties)
+{
+    tetherPotentialProperties_ = tetherPotentialProperties;
+
+    return true;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.H
new file mode 100644
index 0000000000000000000000000000000000000000..4fb5bb8e734e6e4900d78cd50edc1e091045d36f
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.H
@@ -0,0 +1,143 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    tetherPotential
+
+Description
+
+SourceFiles
+    tetherPotential.C
+    newTetherPotential.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef tetherPotential_H
+#define tetherPotential_H
+
+#include "IOdictionary.H"
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "autoPtr.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                Class tetherPotential Declaration
+\*---------------------------------------------------------------------------*/
+
+class tetherPotential
+{
+
+protected:
+
+    // Protected data
+
+        word name_;
+        dictionary tetherPotentialProperties_;
+
+
+    // Private Member Functions
+        
+        //- Disallow copy construct
+        tetherPotential(const tetherPotential&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const tetherPotential&);
+
+public:
+
+    //- Runtime type information
+    TypeName("tetherPotential");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            tetherPotential,
+            dictionary,
+            (
+                const word& name,
+                const dictionary& tetherPotentialProperties
+            ),
+            (name, tetherPotentialProperties)
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected viscosity model
+        static autoPtr<tetherPotential> New
+        (
+            const word& name,
+            const dictionary& tetherPotentialProperties
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        tetherPotential
+        (
+            const word& name,
+            const dictionary& tetherPotentialProperties
+        );
+
+
+    // Destructor
+
+        virtual ~tetherPotential()
+        {}
+
+
+    // Member Functions
+
+        virtual scalar energy (const scalar r) const = 0;
+
+        virtual scalar force (const scalar r) const = 0;
+
+        const dictionary& tetherPotentialProperties() const
+        {
+            return tetherPotentialProperties_;
+        }
+
+        //- Read tetherPotential dictionary
+        virtual bool read(const dictionary& tetherPotentialProperties) = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.C b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.C
new file mode 100644
index 0000000000000000000000000000000000000000..4e5248f96c1e8c60da6cd7628002d499261d604a
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.C
@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "harmonicSpring.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace tetherPotentials
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(harmonicSpring, 0);
+
+addToRunTimeSelectionTable
+(
+    tetherPotential,
+    harmonicSpring,
+    dictionary
+);
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+harmonicSpring::harmonicSpring
+(
+    const word& name,
+    const dictionary& tetherPotentialProperties
+)
+:
+    tetherPotential(name, tetherPotentialProperties),
+    harmonicSpringCoeffs_(tetherPotentialProperties.subDict(typeName + "Coeffs")),
+    springConstant_(readScalar(harmonicSpringCoeffs_.lookup("springConstant")))
+{}
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+scalar harmonicSpring::energy(const scalar r) const
+{
+    return 0.5 * springConstant_ * r * r;
+}
+
+scalar harmonicSpring::force(const scalar r) const
+{
+    return -springConstant_ * r;
+}
+
+bool harmonicSpring::read(const dictionary& tetherPotentialProperties)
+{
+    tetherPotential::read(tetherPotentialProperties);
+
+    harmonicSpringCoeffs_ = tetherPotentialProperties.subDict(typeName + "Coeffs");
+
+    harmonicSpringCoeffs_.lookup("springConstant") >> springConstant_;
+    
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace tetherPotentials
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.H
new file mode 100644
index 0000000000000000000000000000000000000000..92f9908c0abf2f51910a414128bc7049cedb7bdb
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.H
@@ -0,0 +1,104 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    harmonicSpring
+
+Description
+     
+
+SourceFiles
+    harmonicSpring.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef harmonicSpring_H
+#define harmonicSpring_H
+
+#include "tetherPotential.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace tetherPotentials
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class harmonicSpring Declaration
+\*---------------------------------------------------------------------------*/
+
+class harmonicSpring
+:
+    public tetherPotential
+{
+    // Private data
+
+        dictionary harmonicSpringCoeffs_;
+
+        scalar springConstant_;
+
+public:
+
+    //- Runtime type information
+    TypeName("harmonicSpring");
+
+
+    // Constructors
+
+        //- Construct from components
+        harmonicSpring
+        (
+            const word& name,
+            const dictionary& tetherPotentialProperties
+        );
+
+
+    // Destructor
+
+        ~harmonicSpring()
+        {}
+
+
+    // Member Functions
+
+        scalar energy(const scalar r) const;
+
+        scalar force(const scalar r) const;
+
+        //- Read transportProperties dictionary
+        bool read(const dictionary& tetherPotentialProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace tetherPotentials
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.C b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.C
new file mode 100644
index 0000000000000000000000000000000000000000..7baa584b30da4941d574fd282876fdb22d8340e3
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.C
@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "restrainedHarmonicSpring.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace tetherPotentials
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(restrainedHarmonicSpring, 0);
+
+addToRunTimeSelectionTable
+(
+    tetherPotential,
+    restrainedHarmonicSpring,
+    dictionary
+);
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+restrainedHarmonicSpring::restrainedHarmonicSpring
+(
+    const word& name,
+    const dictionary& tetherPotentialProperties
+)
+:
+    tetherPotential(name, tetherPotentialProperties),
+    restrainedHarmonicSpringCoeffs_
+    (
+        tetherPotentialProperties.subDict(typeName + "Coeffs")
+    ),
+    springConstant_
+    (
+        readScalar(restrainedHarmonicSpringCoeffs_.lookup("springConstant"))
+    ),
+    rR_
+    (
+        readScalar(restrainedHarmonicSpringCoeffs_.lookup("rR"))
+    )
+{}
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+scalar restrainedHarmonicSpring::energy(const scalar r) const
+{
+    if (r < rR_)
+    {
+        return 0.5 * springConstant_ * r * r;
+    }
+    else
+    {
+        return 0.5 * springConstant_ * rR_ * rR_
+            + springConstant_ * rR_ * (r - rR_);
+    }
+}
+
+scalar restrainedHarmonicSpring::force(const scalar r) const
+{
+    if (r < rR_)
+    {
+        return -springConstant_ * r;
+    }
+    else
+    {
+        return -springConstant_ * rR_;
+    }
+    
+}
+
+bool restrainedHarmonicSpring::read(const dictionary& tetherPotentialProperties)
+{
+    tetherPotential::read(tetherPotentialProperties);
+
+    restrainedHarmonicSpringCoeffs_ =
+        tetherPotentialProperties.subDict(typeName + "Coeffs");
+
+    restrainedHarmonicSpringCoeffs_.lookup("springConstant") >> springConstant_;
+    restrainedHarmonicSpringCoeffs_.lookup("rR") >> rR_;
+    
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace tetherPotentials
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.H
new file mode 100644
index 0000000000000000000000000000000000000000..021fdd1f6ca2e98d862bc9cd162483038040af37
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.H
@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    restrainedHarmonicSpring
+
+Description
+     
+
+SourceFiles
+    restrainedHarmonicSpring.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef restrainedHarmonicSpring_H
+#define restrainedHarmonicSpring_H
+
+#include "tetherPotential.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace tetherPotentials
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class restrainedHarmonicSpring Declaration
+\*---------------------------------------------------------------------------*/
+
+class restrainedHarmonicSpring
+:
+    public tetherPotential
+{
+    // Private data
+
+        dictionary restrainedHarmonicSpringCoeffs_;
+
+        scalar springConstant_;
+
+        scalar rR_;
+
+public:
+
+    //- Runtime type information
+    TypeName("restrainedHarmonicSpring");
+
+
+    // Constructors
+
+        //- Construct from components
+        restrainedHarmonicSpring
+        (
+            const word& name,
+            const dictionary& tetherPotentialProperties
+        );
+
+
+    // Destructor
+
+        ~restrainedHarmonicSpring()
+        {}
+
+
+    // Member Functions
+
+        scalar energy(const scalar r) const;
+
+        scalar force(const scalar r) const;
+
+        //- Read transportProperties dictionary
+        bool read(const dictionary& tetherPotentialProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace tetherPotentials
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.C b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.C
new file mode 100755
index 0000000000000000000000000000000000000000..8bb11b72cd20397ba8ddac98b9666a269fd5a02a
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.C
@@ -0,0 +1,166 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2005 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "tetherPotentialList.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::tetherPotentialList::readTetherPotentialDict
+(
+    const dictionary& tetherPotentialDict,  
+    const List<word>& idList,
+    const List<label>& tetherIds
+)
+{
+    if (!tetherIds.size())
+    {
+        Info<< nl << "No tethered molecules found." << endl;
+
+        idMap_.setSize(0);
+    }
+    else
+    {
+        Info<< nl << "Building tether potentials." << endl;
+
+        idMap_ = List<label>(idList.size(), -1);
+
+        label tetherMapIndex = 0;
+
+        forAll(tetherIds, t)
+        {
+            const label tetherId = tetherIds[t];
+
+            word tetherPotentialName = idList[tetherId];
+
+            if (!tetherPotentialDict.found(tetherPotentialName))
+            {
+                FatalErrorIn("tetherPotentialList::readTetherPotentialDict")
+                        << nl << "tether potential specification subDict "
+                        << tetherPotentialName << " not found"
+                        << abort(FatalError);
+            }
+
+            (*this).set
+            (
+                tetherMapIndex,
+                tetherPotential::New
+                (
+                    tetherPotentialName,
+                    tetherPotentialDict.subDict(tetherPotentialName)
+                )
+            );
+
+            idMap_[tetherId] = tetherMapIndex;
+
+            tetherMapIndex++;
+
+//             Info<< tetherPotentialName << ": "
+//                 << (*this)[tetherPotentialIndex(tetherId)].tetherPotentialProperties()
+//                 << endl;
+
+        }
+    }
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::tetherPotentialList::tetherPotentialList()
+:
+    PtrList<tetherPotential>(),
+    idMap_()
+{}
+
+
+Foam::tetherPotentialList::tetherPotentialList
+(
+    const dictionary& idListDict,
+    const dictionary& tetherPotentialDict,
+    const List<label>& tetherIds
+)
+:
+    PtrList<tetherPotential>(),
+    idMap_()
+{
+    buildPotentials(idListDict, tetherPotentialDict, tetherIds);
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::tetherPotentialList::~tetherPotentialList()
+{}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::tetherPotentialList::buildPotentials
+(
+    const dictionary& idListDict,
+    const dictionary& tetherPotentialDict,
+    const List<label>& tetherIds
+)
+{
+    setSize(tetherIds.size());
+
+    const List<word>& idList = List<word>(idListDict.lookup("idList"));
+
+    readTetherPotentialDict(tetherPotentialDict, idList, tetherIds);
+}
+
+
+const Foam::tetherPotential& Foam::tetherPotentialList::tetherPotentialFunction
+(
+    const label a
+) const
+{
+    return (*this)[tetherPotentialIndex (a)];
+}
+
+
+Foam::scalar Foam::tetherPotentialList::force
+(
+    const label a,
+    const scalar rITMag
+) const
+{
+    scalar f = (*this)[tetherPotentialIndex (a)].force(rITMag);
+
+    return f;
+}
+
+
+Foam::scalar Foam::tetherPotentialList::energy
+(
+    const label a,
+    const scalar rITMag
+) const
+{
+    scalar e = (*this)[tetherPotentialIndex (a)].energy(rITMag);
+
+    return e;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// ************************************************************************* //
diff --git a/src/lagrangian/molecule/potentials/pairPotential/derived/lennardJones.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.H
similarity index 52%
rename from src/lagrangian/molecule/potentials/pairPotential/derived/lennardJones.H
rename to src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.H
index 57dea3ab725b00c805cf34f9518dad9b57d83516..515a9d07108917e1c2ad7533aef21f6f720a3122 100755
--- a/src/lagrangian/molecule/potentials/pairPotential/derived/lennardJones.H
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.H
@@ -23,113 +23,106 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::pairPotential
+    Foam::tetherPotentialList
 
 Description
 
 SourceFiles
-    pairPotentialI.H
-    pairPotential.C
+    tetherPotentialList.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef pairPotential_H
-#define pairPotential_H
+#ifndef tetherPotentialList_H
+#define tetherPotentialList_H
 
-#include "vector.H"
+#include "PtrList.H"
+#include "word.H"
+#include "tetherPotential.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-
 /*---------------------------------------------------------------------------*\
-                        Class pairPotential Declaration
+                           Class tetherPotentialList Declaration
 \*---------------------------------------------------------------------------*/
 
-class pairPotential
+class tetherPotentialList
+:
+    public PtrList<tetherPotential>
 {
     // Private data
 
-        // At the moment this is hard coded to be a "shiftedLennardJones"
-        // potential.  In future use templated classes virtual functions,
-        // function pointers and all kinds of good stuff to make this a
-        // generic pair force,
-
-        // Not strictly necessary, but useful to keep track of what's what.
-        word pairPotentialName_;
-
-        // To be used when runTime selection of potential function is working.
-        word pairPotentialType_;
+        List<label> idMap_;
 
-        scalar sigma_;
+    // Private Member Functions
 
-        scalar epsilon_;
-
-        scalar rCut_;
-
-        scalar rCutSqr_;
+        inline label tetherPotentialIndex
+        (
+            const label a
+        ) const;
 
-        // rCutShiftedForcePotentialConstTerm
-        scalar rCutSFPotConst_;
+        void readTetherPotentialDict
+        (
+            const dictionary& tetherPotentialDict,
+            const List<word>& idList,
+            const List<label>& tetherIds
+        );
 
-        // rCutShiftedForceConstTerm
-        scalar rCutSFConst_;
+        //- Disallow default bitwise assignment
+        void operator=(const tetherPotentialList&);
 
+        //- Disallow default bitwise copy construct
+        tetherPotentialList(const tetherPotentialList&);
 
 public:
 
     // Constructors
 
-        //- Construct null
-        pairPotential();
+        tetherPotentialList();
 
-        //- Construct from components
-        pairPotential
+        //- Construct from idList and potental dictionaries
+        tetherPotentialList
         (
-            const word& pairPotentialName,
-            const word& pairPotentialType,
-            const scalar sigma,
-            const scalar epsilon,
-            const scalar rCut
+            const dictionary& idListDict,
+            const dictionary& tetherPotentialDict,
+            const List<label>& tetherIds
         );
 
-
     // Destructor
 
-        virtual ~pairPotential();
-
+        ~tetherPotentialList();
 
     // Member Functions
+        
+        void buildPotentials
+        (
+            const dictionary& idListDict,
+            const dictionary& tetherPotentialDict,
+            const List<label>& tetherIds
+        );
 
         // Access
 
-            inline const word& pairPotentialName() const;
-
-            inline const word& pairPotentialType() const;
-
-            inline scalar sigma() const;
-
-            inline scalar epsilon() const;
-
-            inline scalar rCut() const;
-
-            inline scalar rCutSqr() const;
-
-
-        // Write
-
-            virtual void write(Ostream&) const;
-
-            scalar force(const scalar rIJMag) const;
-
-            scalar energy(const scalar rIJMag) const;
+            inline const List<word>& idMap() const;
 
+            const tetherPotential& tetherPotentialFunction
+            (
+                const label a
+            ) const;
 
-    // IOstream Operators
+            scalar force
+            (
+                const label a,
+                const scalar rIJMag
+            ) const;
 
-        friend Ostream& operator<<(Ostream&, const pairPotential&);
+            scalar energy
+            (
+                const label a,
+                const scalar rIJMag
+            ) const;
 };
 
 
@@ -139,7 +132,7 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#include "pairPotentialI.H"
+#include "tetherPotentialListI.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/lagrangian/molecule/potentials/pairPotential/basic/pairPotentialListI.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialListI.H
similarity index 77%
rename from src/lagrangian/molecule/potentials/pairPotential/basic/pairPotentialListI.H
rename to src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialListI.H
index d58fa545df485ccc2984f7aa03fa586155a02372..fb75f16db71ec3cbcdb34e51d94e6ded79947451 100755
--- a/src/lagrangian/molecule/potentials/pairPotential/basic/pairPotentialListI.H
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialListI.H
@@ -26,34 +26,28 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-inline Foam::label Foam::pairPotentialList::pairPotentialIndex
+inline Foam::label Foam::tetherPotentialList::tetherPotentialIndex
 (
-    const label a,
-    const label b
+    const label a
 ) const
 {
-    label index;
+    label index = idMap_[a];
 
-    if (a < b)
+    if (index == -1 || a >= idMap_.size())
     {
-        index = a*(2*nIds_ - a - 1)/2 + b;
+        FatalErrorIn
+        (
+            "Foam::tetherPotentialList::tetherPotentialIndex(const label a)"
+        )
+            << "Attempting to access an undefined tetherPotential."
+            << abort(FatalError);
+
+        return -1;
     }
-
     else
     {
-        index = b*(2*nIds_ - b - 1)/2 + a;
+        return index;
     }
-
-    return index;
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-inline Foam::label Foam::pairPotentialList::nIds() const
-{
-    return nIds_;
 }
 
-
 // ************************************************************************* //
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudBuildCellsFoundInRange.H b/src/lagrangian/molecule/moleculeCloud/moleculeCloudBuildCellsFoundInRange.H
deleted file mode 100755
index 16c58efea3af289a1d8b518ece6b14225cdf8683..0000000000000000000000000000000000000000
--- a/src/lagrangian/molecule/moleculeCloud/moleculeCloudBuildCellsFoundInRange.H
+++ /dev/null
@@ -1,52 +0,0 @@
-forAll (pointsOnThisSegment, pS)
-{
-    const point& ptP = mesh_.points()[pointsOnThisSegment[pS]];
-
-    // Assessing real cells in range is only required on the 1st iteration
-    // because they do not change from iteration to iteration.
-
-    if (iterationNo == 0)
-    {
-        forAll (mesh_.points(), pointIIndex)
-        {
-
-            // Compare separation of ptP to all other points in the mesh,
-            // add unique reference to cell with any point within rCutMax_
-            // to realCellsFoundInRange.
-
-            const point& ptI(mesh_.points()[pointIIndex]);
-
-            if (magSqr(ptP - ptI) <= rCutMaxSqr)
-            {
-                const labelList& ptICells(mesh_.pointCells()[pointIIndex]);
-
-                forAll(ptICells, pIC)
-                {
-                    const label cellI(ptICells[pIC]);
-
-                    if(findIndex(realCellsFoundInRange, cellI) == -1)
-                    {
-                        realCellsFoundInRange.append(cellI);
-                    }
-                }
-            }
-        }
-    }
-
-    forAll(referredInteractionList, rIL)
-    {
-        const vectorList& refCellPoints =
-            referredInteractionList[rIL].vertexPositions();
-
-        forAll(refCellPoints, rCP)
-        {
-            if (magSqr(ptP - refCellPoints[rCP]) <= rCutMaxSqr)
-            {
-                if(findIndex(referredCellsFoundInRange, rIL) == -1)
-                {
-                    referredCellsFoundInRange.append(rIL);
-                }
-            }
-        }
-    }
-}
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudReadMDParameters/moleculeCloudReadIdList.H b/src/lagrangian/molecule/moleculeCloud/moleculeCloudReadMDParameters/moleculeCloudReadIdList.H
deleted file mode 100755
index 8ddc55ed45dcbfe4e545976f849d44bf0d81dcd0..0000000000000000000000000000000000000000
--- a/src/lagrangian/molecule/moleculeCloud/moleculeCloudReadMDParameters/moleculeCloudReadIdList.H
+++ /dev/null
@@ -1,14 +0,0 @@
-IOdictionary idListDict
-(
-    IOobject
-    (
-        "idList",
-        mesh_.time().constant(),
-        mesh_,
-        IOobject::MUST_READ,
-        IOobject::NO_WRITE
-    )
-);
-
-idList_ = List<word>(idListDict.lookup("idList"));
-
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudReadMDParameters/moleculeCloudReadMDParameters.H b/src/lagrangian/molecule/moleculeCloud/moleculeCloudReadMDParameters/moleculeCloudReadMDParameters.H
deleted file mode 100755
index f6d090bfb8c382da5d1b2d75da03e3bda06e1c53..0000000000000000000000000000000000000000
--- a/src/lagrangian/molecule/moleculeCloud/moleculeCloudReadMDParameters/moleculeCloudReadMDParameters.H
+++ /dev/null
@@ -1,6 +0,0 @@
-
-#    include "moleculeCloudReadMDSolution.H"
-
-#    include "moleculeCloudReadIdList.H"
-
-#    include "moleculeCloudReadPotentials.H"
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudReadMDParameters/moleculeCloudReadMDSolution.H b/src/lagrangian/molecule/moleculeCloud/moleculeCloudReadMDParameters/moleculeCloudReadMDSolution.H
deleted file mode 100755
index 07c6b30ff04bf4cac42342ad3e882a5c66499c15..0000000000000000000000000000000000000000
--- a/src/lagrangian/molecule/moleculeCloud/moleculeCloudReadMDParameters/moleculeCloudReadMDSolution.H
+++ /dev/null
@@ -1,27 +0,0 @@
-Info<< nl <<  "Reading MD solution parameters:" << endl;
-
-IOdictionary mdSolution
-(
-    IOobject
-    (
-        "mdSolution",
-        mesh_.time().system(),
-        mesh_,
-        IOobject::MUST_READ,
-        IOobject::NO_WRITE
-    )
-);
-
-integrationMethod_ = integrationMethodNames_.read
-(
-    mdSolution.lookup("integrationMethod")
-);
-Info<< "integrationMethod =\t\t\t"
-    << integrationMethodNames_[integrationMethod_] << endl;
-
-potentialEnergyLimit_ = readScalar(mdSolution.lookup("potentialEnergyLimit"));
-Info<< "potentialEnergyLimit =\t\t\t" << potentialEnergyLimit_ << endl;
-
-guardRadius_ = readScalar(mdSolution.lookup("guardRadius"));
-Info<< "guardRadius =\t\t\t\t" << guardRadius_ << endl;
-
diff --git a/src/lagrangian/molecule/moleculeCloud/moleculeCloudReadMDParameters/moleculeCloudReadPotentials.H b/src/lagrangian/molecule/moleculeCloud/moleculeCloudReadMDParameters/moleculeCloudReadPotentials.H
deleted file mode 100755
index b1bdb2045cce399db459b4b3eab29ede64bebacb..0000000000000000000000000000000000000000
--- a/src/lagrangian/molecule/moleculeCloud/moleculeCloudReadMDParameters/moleculeCloudReadPotentials.H
+++ /dev/null
@@ -1,327 +0,0 @@
-IOdictionary potentialDict
-(
-    IOobject
-    (
-        "potentialDict",
-        mesh_.time().system(),
-        mesh_,
-        IOobject::MUST_READ,
-        IOobject::NO_WRITE
-    )
-);
-
-// ****************************************************************************
-// Pair potentials
-
-Info << nl << "Creating pair potentials:" << nl << endl;
-
-if (!potentialDict.found("pair"))
-{
-    FatalErrorIn("moleculeCloudReadPotentials.H") << nl
-            << "pair potential specification subDict not found"
-            << abort(FatalError);
-}
-
-const dictionary& pairDict = potentialDict.subDict("pair");
-
-pairPotentials_.setSizeAndNIds(nIds());
-
-label a;
-
-label b;
-
-rCutMax_ = 0.0;
-
-for
-(
-    a = 0;
-    a < idList_.size();
-    ++a
-)
-{
-    word idA = idList_[a];
-
-    for
-    (
-        b = a;
-        b < idList_.size();
-        ++b
-    )
-    {
-        word idB = idList_[b];
-
-        word pairPotentialName;
-
-        if (a==b)
-        {
-
-
-            if(pairDict.found(idA+"-"+idB))
-            {
-                pairPotentialName = idA+"-"+idB;
-            }
-
-            else
-            {
-                FatalErrorIn("moleculeCloudReadPotentials.H") << nl
-                        << "Pair pairPotential specification subDict "
-                        << idA+"-"+idB << " not found"
-                        << abort(FatalError);
-            }
-        }
-
-        else
-        {
-            if(pairDict.found(idA+"-"+idB))
-            {
-                pairPotentialName = idA+"-"+idB;
-            }
-
-            else if(pairDict.found(idB+"-"+idA))
-            {
-                pairPotentialName = idB+"-"+idA;
-            }
-
-            else
-            {
-                FatalErrorIn("moleculeCloudReadPotentials.H") << nl
-                        << "Pair pairPotential specification subDict "
-                        << idA+"-"+idB << " or " << idB+"-"+idA << " not found"
-                        << abort(FatalError);
-            }
-
-            if
-            (
-                pairDict.found(idA+"-"+idB)
-                && pairDict.found(idB+"-"+idA)
-            )
-            {
-                FatalErrorIn("moleculeCloudReadPotentials.H") << nl
-                        << "Pair pairPotential specification subDict "
-                        << idA+"-"+idB << " and "
-                        << idB+"-"+idA << " found - multiple definition"
-                        << abort(FatalError);
-            }
-        }
-
-        word pairPotentialType
-        (
-            pairDict.subDict(pairPotentialName).lookup("potentialType")
-        );
-
-        if
-        (
-            pairPotentialType == "maitlandSmithTabulated"
-        )
-        {
-            scalar rCut = readScalar
-            (
-                pairDict.subDict(pairPotentialName).lookup("rCut")
-            );
-
-            if(rCut > rCutMax_)
-            {
-                rCutMax_ = rCut;
-            }
-
-            pairPotentials_.addPotential
-            (
-                a,
-                b,
-                pairPotential
-                (
-                    pairDict.subDict(pairPotentialName)
-                )
-            );
-        }
-        else
-        {
-            FatalErrorIn("moleculeCloudReadPotentials.H") << nl
-                    << "pairPotentialType "<< pairPotentialType << " unknown"
-                    << abort(FatalError);
-        }
-
-        Info << pairPotentialName << ": "
-            << pairPotentials_.pairPotentialFunction(a,b) << endl;
-    }
-}
-
-rCutMax_ += guardRadius_;
-
-rCutMaxSqr_ = rCutMax_ * rCutMax_;
-
-Info << nl << "rCutMax + guardRadius = " << rCutMax_ << endl;
-
-if (potentialDict.found("removalOrder"))
-{
-    List<word> remOrd = potentialDict.lookup("removalOrder");
-
-    removalOrder_.setSize(remOrd.size());
-
-    forAll(removalOrder_, rO)
-    {
-        removalOrder_[rO] = findIndex(idList_, remOrd[rO]);
-    }
-}
-
-// ****************************************************************************
-// Tether potentials
-
-iterator mol(this->begin());
-
-DynamicList<label> tetherIds;
-
-for
-(
- mol = this->begin();
- mol != this->end();
- ++mol
-)
-{
-    if (mol().tethered())
-    {
-        if (findIndex(tetherIds, mol().id()) == -1)
-        {
-            tetherIds.append
-            (
-                mol().id()
-            );
-        }
-    }
-}
-
-if (Pstream::parRun())
-{
-    List< labelList > allTetherIds(Pstream::nProcs());
-
-    allTetherIds[Pstream::myProcNo()] = tetherIds;
-
-    Pstream::gatherList(allTetherIds);
-
-    if (Pstream::master())
-    {
-        DynamicList<label> globalTetherIds;
-
-        forAll(allTetherIds, procN)
-        {
-            const labelList& procNTetherIds = allTetherIds[procN];
-
-            forAll(procNTetherIds, id)
-            {
-                if (findIndex(globalTetherIds, procNTetherIds[id]) == -1)
-                {
-                    globalTetherIds.append
-                    (
-                        procNTetherIds[id]
-                    );
-                }
-            }
-        }
-
-        globalTetherIds.shrink();
-
-        tetherIds = globalTetherIds;
-    }
-
-    Pstream::scatter(tetherIds);
-}
-
-tetherIds.shrink();
-
-if (tetherIds.size())
-{
-    Info << nl << "Creating tether potentials:" << endl;
-
-    if (!potentialDict.found("tether"))
-    {
-        FatalErrorIn("moleculeCloudReadPotentials.H") << nl
-                << "tether potential specification subDict not found"
-                << abort(FatalError);
-    }
-
-    const dictionary& tetherDict = potentialDict.subDict("tether");
-
-    tetherPotentials_.setSizeAndTetherIds(tetherIds);
-
-    forAll (tetherIds, tid)
-    {
-        word tetherPotentialName = idList_[tetherIds[tid]];
-
-        if (!tetherDict.found(tetherPotentialName))
-        {
-            FatalErrorIn("moleculeCloudReadPotentials.H") << nl
-                    << "tether potential specification subDict "
-                    << tetherPotentialName << " not found"
-                    << abort(FatalError);
-        }
-
-        Info << nl << tetherPotentialName << endl;
-
-        word tetherPotentialType =
-            tetherDict.subDict(tetherPotentialName).lookup("potentialType");
-
-        if
-        (
-            tetherPotentialType == "harmonicSpring"
-        )
-        {
-
-            Info << "\tpotentialType = " << tetherPotentialType << endl;
-
-            scalar springConstant
-            (
-                readScalar
-                (
-                    tetherDict.subDict(tetherPotentialName).lookup
-                    (
-                        "springConstant"
-                    )
-                )
-            );
-
-            Info << "\tspringConstant = " << springConstant << endl;
-
-            tetherPotentials_.addPotential
-            (
-                tetherIds[tid],
-                tetherPotential
-                (
-                    tetherPotentialName,
-                    tetherPotentialType,
-                    springConstant
-                )
-            );
-        }
-
-        else
-        {
-            FatalErrorIn("moleculeCloudReadPotentials.H") << nl
-                    << "tetherPotentialType "
-                    << tetherPotentialType << " unknown"
-                    << abort(FatalError);
-        }
-    }
-}
-
-// ****************************************************************************
-// External Forces
-
-gravity_ = vector::zero;
-
-if (potentialDict.found("external"))
-{
-
-    Info << nl << "Reading external forces:" << endl;
-
-    const dictionary& externalDict = potentialDict.subDict("external");
-
-    // ************************************************************************
-    // gravity
-
-    if (externalDict.found("gravity"))
-    {
-        gravity_ = externalDict.lookup("gravity");
-    }
-}
-
-Info << nl << tab << "gravity = " << gravity_ << endl;
diff --git a/src/lagrangian/molecule/potentials/pairPotential/basic/pairPotential.C b/src/lagrangian/molecule/potentials/pairPotential/basic/pairPotential.C
deleted file mode 100755
index f7c96414e4282553fe26bfe9dfd39318035ddcde..0000000000000000000000000000000000000000
--- a/src/lagrangian/molecule/potentials/pairPotential/basic/pairPotential.C
+++ /dev/null
@@ -1,308 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2005 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*----------------------------------------------------------------------------*/
-
-#include "pairPotential.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::pairPotential::setLookupTables()
-{
-    label N = label((rCut_ - rMin_)/dr_) + 1;
-
-    forceLookup_.setSize(N);
-
-    energyLookup_.setSize(N);
-
-    forAll(forceLookup_, k)
-    {
-        forceLookup_[k] = force(k*dr_ + rMin_);
-
-        energyLookup_[k] = energy(k*dr_ + rMin_);
-    }
-
-    forceLookup_[N-1] = 0.0;
-    energyLookup_[N-1] = 0.0;
-}
-
-
-void Foam::pairPotential::setConstants()
-{
-    scalar nr = n(rCut_);
-
-    u_at_rCut_ =
-        epsilon_
-       *(
-            (6.0/(nr - 6.0))*Foam::pow( (rCut_/rm_), -nr)
-          - (nr/(nr - 6.0))*Foam::pow( (rCut_/rm_), -6)
-        );
-
-    du_by_dr_at_rCut_ =
-        -6.0 * epsilon_ * gamma_
-       *(
-            Foam::pow( (rCut_/rm_),-nr)
-           *(
-                (rCut_/rm_)
-               *(1.0 / (nr - 6.0) + log(rCut_ / rm_) + 1.0)
-              + (m_/gamma_)
-              - 1.0
-            )
-          - Foam::pow((rCut_/rm_),-6.0)
-            *(rCut_ / ( (nr - 6.0) * rm_) + (nr/gamma_))
-        )
-       /(rCut_*(nr - 6.0));
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::pairPotential::pairPotential()
-{}
-
-
-Foam::pairPotential::pairPotential(const dictionary& pPDict)
-:
-    forceLookup_(),
-    energyLookup_(),
-    m_(),
-    gamma_(),
-    rm_(),
-    epsilon_(),
-    rCut_(),
-    rCutSqr_(),
-    u_at_rCut_(),
-    du_by_dr_at_rCut_(),
-    rMin_(),
-    dr_()
-{
-
-    m_ = readScalar(pPDict.lookup("m"));
-
-    gamma_ = readScalar(pPDict.lookup("gamma"));
-
-    rm_ = readScalar(pPDict.lookup("rm"));
-
-    epsilon_ = readScalar(pPDict.lookup("epsilon"));
-
-    rCut_ = readScalar(pPDict.lookup("rCut"));
-
-    rCutSqr_ = rCut_ * rCut_;
-
-    rMin_ = readScalar(pPDict.lookup("rMin"));
-
-    dr_ = readScalar(pPDict.lookup("dr"));
-
-    setConstants();
-
-    setLookupTables();
-}
-
-
-Foam::pairPotential::pairPotential
-(
-    const scalar m,
-    const scalar gamma,
-    const scalar rm,
-    const scalar epsilon,
-    const scalar rCut,
-    const scalar rMin,
-    const scalar dr
-)
-:
-    forceLookup_(),
-    energyLookup_(),
-    m_(m),
-    gamma_(gamma),
-    rm_(rm),
-    epsilon_(epsilon),
-    rCut_(rCut),
-    rCutSqr_(rCut*rCut),
-    u_at_rCut_(),
-    du_by_dr_at_rCut_(),
-    rMin_(rMin),
-    dr_(dr)
-{
-    setConstants();
-
-    setLookupTables();
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::pairPotential::~pairPotential()
-{}
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::pairPotential::write(Ostream& os) const
-{
-    os<< "Shifted force Maitland Smith pair potential."
-        << nl << tab << "m = " << m_
-        << nl << tab << "gamma = " << gamma_
-        << nl << tab << "rm = " << rm_
-        << nl << tab << "epsilon = " << epsilon_
-        << nl << tab << "rCut = " << rCut_
-        << nl << tab << "rMin = " << rMin_
-        << nl << tab << "dr = " << dr_
-        << endl;
-}
-
-
-Foam::scalar Foam::pairPotential::n(const scalar rIJMag) const
-{
-    return (m_ + gamma_*(rIJMag/rm_ - 1.0));
-}
-
-
-Foam::scalar Foam::pairPotential::force(const scalar rIJMag) const
-{
-    scalar nr(n(rIJMag));
-
-    scalar r = rIJMag/rm_;
-
-    scalar rN = Foam::pow(r,-nr);
-
-    scalar f =
-      - epsilon_
-       *(
-           - 6.0/(nr - 6.0)/(nr - 6.0)*(gamma_/rm_)*rN
-           + 6.0/(nr - 6.0)*rN
-             *((gamma_ - m_)/rIJMag - (gamma_/rm_)*(log(r) + 1.0))
-           + 6.0/(nr - 6.0)
-             *(1.0/(nr - 6.0)*(gamma_/rm_) + nr/rIJMag)*Foam::pow(r, -6)
-        )
-      + du_by_dr_at_rCut_;
-
-    return f;
-}
-
-
-Foam::scalar Foam::pairPotential::forceLookup(const scalar rIJMag) const
-{
-    scalar k_rIJ = (rIJMag - rMin_)/dr_;
-
-    label k(k_rIJ);
-
-    if (k < 0)
-    {
-        FatalErrorIn("pairPotential.C") << nl
-        << "rIJMag less than rMin" << nl
-        << abort(FatalError);
-    }
-
-    scalar f =
-        (k_rIJ - k)*forceLookup_[k+1]
-      + (k + 1 - k_rIJ)*forceLookup_[k];
-
-    return f;
-}
-
-
-Foam::List< Foam::Pair< Foam::scalar > >
-Foam::pairPotential::forceTable() const
-{
-    List< Pair<scalar> > forceTab(forceLookup_.size());
-
-    forAll(forceLookup_,k)
-    {
-        forceTab[k].first() = rMin_ + k*dr_;
-
-        forceTab[k].second() = forceLookup_[k];
-    }
-
-    return forceTab;
-}
-
-
-Foam::scalar Foam::pairPotential::energy(const scalar rIJMag) const
-{
-    scalar nr = n(rIJMag);
-
-    scalar e =
-        epsilon_
-       *(
-            (6.0 / (nr - 6.0))*Foam::pow( rIJMag/rm_, -nr)
-          - (nr / (nr - 6.0))*Foam::pow( rIJMag/rm_, -6)
-        )
-      - u_at_rCut_
-      - (rIJMag - rCut_)
-       *du_by_dr_at_rCut_;
-
-    return e;
-}
-
-
-Foam::scalar Foam::pairPotential::energyLookup(const scalar rIJMag) const
-{
-    scalar k_rIJ = (rIJMag - rMin_)/dr_;
-
-    label k(k_rIJ);
-
-    if (k < 0)
-    {
-        FatalErrorIn("pairPotential.C") << nl
-        << "rIJMag less than rMin" << nl
-        << abort(FatalError);
-    }
-
-    scalar e =
-        (k_rIJ - k)*energyLookup_[k+1]
-      + (k + 1 - k_rIJ)*energyLookup_[k];
-
-    return e;
-}
-
-
-Foam::List< Foam::Pair< Foam::scalar > >
-    Foam::pairPotential::energyTable() const
-{
-    List< Pair<scalar> > energyTab(energyLookup_.size());
-
-    forAll(energyLookup_,k)
-    {
-        energyTab[k].first() = rMin_ + k*dr_;
-
-        energyTab[k].second() = energyLookup_[k];
-    }
-
-    return energyTab;
-}
-
-
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
-
-Foam::Ostream& Foam::operator<<(Ostream& os, const pairPotential& pot)
-{
-    pot.write(os);
-    os.check
-    (
-        "Foam::Ostream& Foam::pperator<<(Ostream& f, const pairPotential& pot"
-    );
-    return os;
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecule/potentials/pairPotential/basic/pairPotentialList.C b/src/lagrangian/molecule/potentials/pairPotential/basic/pairPotentialList.C
deleted file mode 100755
index 7728490ee0eda57b4bc4502c4b6eab5eb87eefa0..0000000000000000000000000000000000000000
--- a/src/lagrangian/molecule/potentials/pairPotential/basic/pairPotentialList.C
+++ /dev/null
@@ -1,180 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2005 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*----------------------------------------------------------------------------*/
-
-#include "pairPotentialList.H"
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::pairPotentialList::pairPotentialList()
-:
-    List<pairPotential> ()
-{}
-
-
-Foam::pairPotentialList::pairPotentialList
-(
-    const label nIds
-)
-:
-    List<pairPotential> ((nIds * (nIds + 1))/2),
-    nIds_(nIds)
-{}
-
-
-Foam::pairPotentialList::pairPotentialList
-(
-    const List<pairPotential>& pairPotentials,
-    const label nIds
-)
-:
-    List<pairPotential> (pairPotentials),
-    nIds_(nIds)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-pairPotentialList::~pairPotentialList()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void pairPotentialList::setSizeAndNIds (const label nIds)
-{
-    nIds_ = nIds;
-    setSize((nIds * (nIds + 1))/2);
-}
-
-
-void pairPotentialList::addPotential
-(
-    const label a,
-    const label b,
-    const pairPotential& pot
-)
-{
-    if (pairPotentialIndex (a,b) > size()-1)
-    {
-        FatalErrorIn
-        (
-            "Foam::pairPotentialList::addPotential "
-            "(const label a, const label b, const pairPotential& pot)"
-        )<< "Attempting to add a pairPotential with too high an index"
-         << nl
-         << "Check if the pairPotentialList has been constructed "
-         << "with the number of ids to expect"
-         << nl << abort(FatalError);
-    }
-    else
-    {
-        (*this)[pairPotentialIndex (a,b)] = pot;
-    }
-}
-
-
-const pairPotential& pairPotentialList::pairPotentialFunction
-(
-    const label a,
-    const label b
-) const
-{
-    return (*this)[pairPotentialIndex (a,b)];
-}
-
-
-bool pairPotentialList::rCutSqr
-(
-    const label a,
-    const label b,
-    const scalar rIJMagSqr
-) const
-{
-    if(rIJMagSqr <= rCutSqr (a,b))
-    {
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-const scalar pairPotentialList::rCutSqr
-(
-    const label a,
-    const label b
-) const
-{
-    return (*this)[pairPotentialIndex (a,b)].rCutSqr();
-}
-
-
-const scalar pairPotentialList::rCut
-(
-    const label a,
-    const label b
-) const
-{
-    return (*this)[pairPotentialIndex (a,b)].rCut();
-}
-
-
-scalar pairPotentialList::force
-(
-    const label a,
-    const label b,
-    const scalar rIJMag
-) const
-{
-    scalar f = (*this)[pairPotentialIndex (a,b)].forceLookup(rIJMag);
-
-    return f;
-}
-
-
-scalar pairPotentialList::energy
-(
-    const label a,
-    const label b,
-    const scalar rIJMag
-) const
-{
-    scalar e = (*this)[pairPotentialIndex (a,b)].energyLookup(rIJMag);
-
-    return e;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecule/potentials/pairPotential/derived/lennardJones.C b/src/lagrangian/molecule/potentials/pairPotential/derived/lennardJones.C
deleted file mode 100755
index 1eab1aa40a1fb753f089f7257de7f4f3387bcac5..0000000000000000000000000000000000000000
--- a/src/lagrangian/molecule/potentials/pairPotential/derived/lennardJones.C
+++ /dev/null
@@ -1,132 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2005 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*----------------------------------------------------------------------------*/
-
-#include "pairPotential.H"
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-pairPotential::pairPotential()
-{}
-
-
-// Construct from components
-pairPotential::pairPotential
-(
-    const word& pairPotentialName,
-    const word& pairPotentialType,
-    const scalar sigma,
-    const scalar epsilon,
-    const scalar rCut
-)
-:
-    pairPotentialName_(pairPotentialName),
-    pairPotentialType_(pairPotentialType),
-    sigma_(sigma),
-    epsilon_(epsilon),
-    rCut_(rCut),
-    rCutSqr_(rCut*rCut)
-{
-    // rCutShiftedForcePotentialConstTerm
-    rCutSFPotConst_ = (pow((rCut_/sigma_),-12) - pow((rCut_/sigma_),-6));
-
-    // rCutShiftedForceConstTerm
-    rCutSFConst_ =
-        (pow((rCut_/sigma_),-14) - 0.5*pow((rCut_/sigma_),-8))*rCut_;
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-pairPotential::~pairPotential()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void pairPotential::write(Ostream& os) const
-{
-    os << pairPotentialName_
-       << nl << pairPotentialType_
-       << nl << sigma_
-       << nl << epsilon_
-       << nl << rCut_
-       << endl;
-}
-
-
-scalar pairPotential::force(const scalar rIJMag) const
-{
-    // (rIJ/sigma)^-2
-    scalar rIJoSMagSqInv = (sigma_/rIJMag)*(sigma_/rIJMag);
-
-    // (rIJ/sigma)^-6
-    scalar rIJoSMagSqInvCu = rIJoSMagSqInv*rIJoSMagSqInv*rIJoSMagSqInv;
-
-    scalar f = 48.0*epsilon_/(sigma_*sigma_)
-               *(rIJoSMagSqInvCu*(rIJoSMagSqInvCu - 0.5)*rIJoSMagSqInv
-             - rCutSFConst_/rIJMag);
-
-    return f;
-}
-
-
-scalar pairPotential::energy(const scalar rIJMag) const
-{
-    // (rIJ/sigma)^-2
-    scalar rIJoSMagSqInv = (sigma_/rIJMag)*(sigma_/rIJMag);
-
-    // (rIJ/sigma)^-6
-    scalar rIJoSMagSqInvCu = rIJoSMagSqInv*rIJoSMagSqInv*rIJoSMagSqInv;
-
-    scalar e = 4*epsilon_
-               *(
-                    rIJoSMagSqInvCu*(rIJoSMagSqInvCu - 1.0)
-                  - rCutSFPotConst_
-                  + 12*(rIJMag - rCut_)/(sigma_*sigma_)*rCutSFConst_
-                );
-
-    return e;
-}
-
-
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
-
-Ostream& operator<<(Ostream& os, const pairPotential& pot)
-{
-    pot.write(os);
-    os.check("Ostream& operator<<(Ostream& f, const pairPotential& pot");
-    return os;
-}
-
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //
diff --git a/tutorials/mdEquilibrationFoam/Allclean b/tutorials/mdEquilibrationFoam/Allclean
new file mode 100755
index 0000000000000000000000000000000000000000..3d688adfb3e4b393c27afd30a49c00e5fcc7ef7c
--- /dev/null
+++ b/tutorials/mdEquilibrationFoam/Allclean
@@ -0,0 +1,13 @@
+#!/bin/sh
+
+currDir=`pwd`
+application=`basename $currDir`
+cases="periodicCube"
+
+tutorialPath=`dirname $0`/..
+. $tutorialPath/CleanFunctions
+
+for case in $cases
+do
+    cleanCase $case
+done
diff --git a/tutorials/mdEquilibrationFoam/Allrun b/tutorials/mdEquilibrationFoam/Allrun
new file mode 100755
index 0000000000000000000000000000000000000000..eefbf1aded0c0a5d5123f285ce3e6a1d6d64a8f3
--- /dev/null
+++ b/tutorials/mdEquilibrationFoam/Allrun
@@ -0,0 +1,19 @@
+#!/bin/sh
+
+currDir=`pwd`
+application=`basename $currDir`
+cases="periodicCube"
+
+tutorialPath=`dirname $0`/..
+. $tutorialPath/RunFunctions
+. $tutorialPath/CleanFunctions
+
+
+for case in $cases
+do
+    runApplication blockMesh $case
+
+    runApplication molConfig $case
+
+    runApplication $application $case
+done
diff --git a/tutorials/mdEquilibrationFoam/periodicCube/Allrun b/tutorials/mdEquilibrationFoam/periodicCube/Allrun
deleted file mode 100755
index 2121179f1cf87b8138d978cfb95c245e29edbff1..0000000000000000000000000000000000000000
--- a/tutorials/mdEquilibrationFoam/periodicCube/Allrun
+++ /dev/null
@@ -1,10 +0,0 @@
-#!/bin/sh
-# Source tutorial run functions
-. $WM_PROJECT_DIR/bin/tools/RunFunctions
-
-# Get application name from directory
-application=`basename $PWD`
-
-runApplication blockMesh
-runApplication molConfig
-runApplication $application
diff --git a/tutorials/mdEquilibrationFoam/periodicCube/system/mdSolution b/tutorials/mdEquilibrationFoam/periodicCube/system/mdSolution
index 681963175de629fcb9ebab3a3a3ebd0b5893e716..c72bb260facf92c9cacb0cc02e4f7679d0e35926 100644
--- a/tutorials/mdEquilibrationFoam/periodicCube/system/mdSolution
+++ b/tutorials/mdEquilibrationFoam/periodicCube/system/mdSolution
@@ -26,6 +26,4 @@ integrationMethod               verletLeapfrog;
 
 potentialEnergyLimit            5.256e-20;
 
-guardRadius                     0;
-
 // ************************************************************************* //
\ No newline at end of file
diff --git a/tutorials/mdEquilibrationFoam/periodicCube/system/potentialDict b/tutorials/mdEquilibrationFoam/periodicCube/system/potentialDict
index 78fa1c77465f435999e8df2b11c6388c64fa3b6b..f4d27d298b0f82c26d5dd86fc1be3a153f1fa708 100644
--- a/tutorials/mdEquilibrationFoam/periodicCube/system/potentialDict
+++ b/tutorials/mdEquilibrationFoam/periodicCube/system/potentialDict
@@ -22,7 +22,7 @@ FoamFile
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Subdictionaries specifying types of intermoleular potential.
+// Subdictionaries specifying types of intermolecular potential.
 // Sub-sub dictionaries specify the potentials themselves.
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -52,14 +52,24 @@ pair
 {
     Ar-Ar
     {
-        potentialType   maitlandSmithTabulated;
-        m               13.0;
-        gamma           7.5;
-        rm              0.3756e-9;
-        epsilon         1.990108438e-21;
-        rCut            1.0e-9;
+        pairPotential   azizChen;
+        rCut            1.2e-9;
         rMin            0.15e-9;
         dr              2e-14;
+        azizChenCoeffs
+        {
+            epsilon     1.97742255e-21;
+            rm          3.759e-10;
+            A           0.9502720e7;
+            alpha       16.345655;
+            C6          1.0914254;
+            C8          0.6002595;
+            C10         0.3700113;
+            D           1.4;
+            gamma       2.0;
+        }
+        energyScalingFunction   noScaling;
+        writeTables     yes;
     }
 }
 
@@ -70,8 +80,12 @@ tether
 {
     Ar
     {
-        potentialType   harmonicSpring;
-        springConstant  0.0277;
+        tetherPotential restrainedHarmonicSpring;
+        restrainedHarmonicSpringCoeffs
+        {
+            springConstant  0.0277;
+            rR              1.2e-9;
+        }
     }
 }