recovered bad master merge

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

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

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

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

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

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
 

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;
 

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_;
 }

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;
 

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_;
 }

src/lagrangian/molecularDynamics/molecule/distribution/distribution.H

@@ -104,7 +104,7 @@ public:
 
         // Access
 
-            inline const scalar binWidth() const;
+            inline scalar binWidth() const;
 
 
     // Member Operators

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_;
 }

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

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;

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,

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

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;

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

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

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;

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

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

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

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

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

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 ;

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

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)
                         {

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

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

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

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
+

src/lagrangian/molecularDynamics/potential/Make/options

@@ -0,0 +1,6 @@
+EXE_INC = \
+    -I.. \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+LIB_LIBS = \
+    -lfiniteVolume

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

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

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

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

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

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

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

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

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

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

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