diff --git a/applications/solvers/molecularDynamics/mdFoam/mdFoam.C b/applications/solvers/molecularDynamics/mdFoam/mdFoam.C
index a299d558fd..99a9e90d7d 100644
--- a/applications/solvers/molecularDynamics/mdFoam/mdFoam.C
+++ b/applications/solvers/molecularDynamics/mdFoam/mdFoam.C
@@ -43,10 +43,6 @@ int main(int argc, char *argv[])
 
     potential pot(mesh);
 
-    const pairPotentialList& pairPot(pot.pairPotentials());
-
-    Info<< pairPot.energy(0, 0, 0.45e-9) << endl;
-
     moleculeCloud molecules(mesh, pot);
 
     Info << "\nStarting time loop\n" << endl;
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
index 8fb1f5ebab..89e8483e32 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
@@ -86,6 +86,10 @@ public:
 
             List<label> siteIds_;
 
+            List<bool> pairPotentialSites_;
+
+            List<bool> electrostaticSites_;
+
             diagTensor momentOfInertia_;
 
             scalar mass_;
@@ -94,21 +98,17 @@ public:
 
             void checkSiteListSizes() const;
 
+            void setInteracionSiteBools
+            (
+                const List<word>& siteIds,
+                const List<word>& pairPotSiteIds
+            );
+
     public:
 
         inline constantProperties();
 
-        //- Constructor
-        inline constantProperties
-        (
-            const List<vector>& siteReferencePositions,
-            const List<scalar>& siteCharges,
-            const List<label>& siteIds,
-            const diagTensor& momentOfInertia,
-            scalar mass
-        );
-
-        //- Constructor
+        //- Construct from dictionary
         inline constantProperties(const dictionary& dict);
 
         // Member functions
@@ -120,6 +120,10 @@ public:
         inline const List<label>& siteIds() const;
         inline List<label>& siteIds();
 
+        inline const List<bool>& pairPotentialSites() const;
+
+        inline const List<bool>& electrostaticSites() const;
+
         inline const diagTensor& momentOfInertia() const;
 
         inline scalar mass() const;
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
index 1b65f447e9..b1ee4d5814 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
@@ -31,30 +31,13 @@ inline Foam::molecule::constantProperties::constantProperties()
     siteReferencePositions_(List<vector>(0)),
     siteCharges_(List<scalar>(0)),
     siteIds_(List<label>(0)),
+    pairPotentialSites_(List<bool>(false)),
+    electrostaticSites_(List<bool>(false)),
     momentOfInertia_(diagTensor(0, 0, 0)),
     mass_(0)
 {}
 
 
-inline Foam::molecule::constantProperties::constantProperties
-(
-    const List<vector>& siteReferencePositions,
-    const List<scalar>& siteCharges,
-    const List<label>& siteIds,
-    const diagTensor& momentOfInertia,
-    scalar mass
-)
-:
-    siteReferencePositions_(siteReferencePositions),
-    siteCharges_(siteCharges),
-    siteIds_(siteIds),
-    momentOfInertia_(momentOfInertia),
-    mass_(mass)
-{
-    checkSiteListSizes();
-}
-
-
 inline Foam::molecule::constantProperties::constantProperties
 (
     const dictionary& dict
@@ -63,10 +46,18 @@ inline Foam::molecule::constantProperties::constantProperties
     siteReferencePositions_(dict.lookup("siteReferencePositions")),
     siteCharges_(dict.lookup("siteCharges")),
     siteIds_(List<word>(dict.lookup("siteIds")).size(), -1),
+    pairPotentialSites_(),
+    electrostaticSites_(),
     momentOfInertia_(dict.lookup("momentOfInertia")),
     mass_(readScalar(dict.lookup("mass")))
 {
     checkSiteListSizes();
+
+    setInteracionSiteBools
+    (
+        List<word>(dict.lookup("siteIds")),
+        List<word>(dict.lookup("pairPotentialSiteIds"))
+    );
 }
 
 
@@ -124,38 +115,75 @@ inline void Foam::molecule::constantProperties::checkSiteListSizes() const
 }
 
 
+inline void Foam::molecule::constantProperties::setInteracionSiteBools
+(
+    const List<word>& siteIds,
+    const List<word>& pairPotSiteIds
+)
+{
+    pairPotentialSites_.setSize(siteIds_.size());
+
+    electrostaticSites_.setSize(siteIds_.size());
+
+    forAll(siteIds_, i)
+    {
+        const word& id(siteIds[i]);
+
+        pairPotentialSites_[i] = (findIndex(pairPotSiteIds, id) > -1);
+
+        electrostaticSites_[i] = (mag(siteCharges_[i]) > VSMALL);
+    }
+
+    Info<< pairPotentialSites_ << nl << electrostaticSites_ << endl;
+}
+
+
 // * * * * * * * constantProperties Member Functions * * * * * * * * * * * * //
 
 inline const Foam::List<Foam::vector>&
-    Foam::molecule::constantProperties::siteReferencePositions() const
+Foam::molecule::constantProperties::siteReferencePositions() const
 {
     return siteReferencePositions_;
 }
 
 
 inline const Foam::List<Foam::scalar>&
-    Foam::molecule::constantProperties::siteCharges() const
+Foam::molecule::constantProperties::siteCharges() const
 {
     return siteCharges_;
 }
 
 
 inline const Foam::List<Foam::label>&
-    Foam::molecule::constantProperties::siteIds() const
+Foam::molecule::constantProperties::siteIds() const
 {
     return siteIds_;
 }
 
 
 inline Foam::List<Foam::label>&
-    Foam::molecule::constantProperties::siteIds()
+Foam::molecule::constantProperties::siteIds()
 {
     return siteIds_;
 }
 
 
+inline const Foam::List<bool>&
+Foam::molecule::constantProperties::pairPotentialSites() const
+{
+    return pairPotentialSites_;
+}
+
+
+inline const Foam::List<bool>&
+Foam::molecule::constantProperties::electrostaticSites() const
+{
+    return electrostaticSites_;
+}
+
+
 inline const Foam::diagTensor&
-    Foam::molecule::constantProperties::momentOfInertia() const
+Foam::molecule::constantProperties::momentOfInertia() const
 {
     return momentOfInertia_;
 }
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
index b4a411ea14..3c95d80ca7 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
@@ -52,7 +52,7 @@ void Foam::moleculeCloud::buildConstProps()
 
     constPropList_.setSize(idList.size());
 
-    const List<word>& allSiteIdNames(pot_.allSiteIdNames());
+    const List<word>& siteIdList(pot_.siteIdList());
 
     IOdictionary moleculePropertiesDict
     (
@@ -81,7 +81,7 @@ void Foam::moleculeCloud::buildConstProps()
         {
             const word& siteId = siteIdNames[sI];
 
-            siteIds[sI] = findIndex(allSiteIdNames, siteId);
+            siteIds[sI] = findIndex(siteIdList, siteId);
 
             if (siteIds[sI] == -1)
             {
@@ -96,6 +96,8 @@ void Foam::moleculeCloud::buildConstProps()
         constProp = molecule::constantProperties(molDict);
 
         constProp.siteIds() = siteIds;
+
+        Info<< "sites " << constProp.siteIds() << endl;
     }
 }
 
diff --git a/src/lagrangian/molecularDynamics/potential/potential/potential.C b/src/lagrangian/molecularDynamics/potential/potential/potential.C
index c045a5fa48..0767e25812 100644
--- a/src/lagrangian/molecularDynamics/potential/potential/potential.C
+++ b/src/lagrangian/molecularDynamics/potential/potential/potential.C
@@ -59,7 +59,9 @@ void Foam::potential::potential::readPotentialDict()
         )
     );
 
-    DynamicList<word> allSiteIdNames;
+    DynamicList<word> siteIdList;
+
+    DynamicList<word> pairPotentialSiteIdList;
 
     forAll(idList_, i)
     {
@@ -80,22 +82,58 @@ void Foam::potential::potential::readPotentialDict()
         {
             const word& siteId = siteIdNames[sI];
 
-            if(findIndex(allSiteIdNames, siteId) == -1)
+            if(findIndex(siteIdList, siteId) == -1)
             {
-                allSiteIdNames.append(siteId);
+                siteIdList.append(siteId);
+            }
+        }
+
+        List<word> pairPotSiteIds = molDict.lookup("pairPotentialSiteIds");
+
+        forAll(pairPotSiteIds, sI)
+        {
+            const word& siteId = pairPotSiteIds[sI];
+
+            if(findIndex(siteIdNames, siteId) == -1)
+            {
+                FatalErrorIn("potential.C") << nl
+                    << siteId << " in pairPotentialSiteIds is not in siteIds: "
+                    << siteIdNames << nl << abort(FatalError);
+            }
+
+            if(findIndex(pairPotentialSiteIdList, siteId) == -1)
+            {
+                pairPotentialSiteIdList.append(siteId);
             }
         }
     }
 
-    allSiteIdNames_.transfer(allSiteIdNames.shrink());
+    label nPairPotIds_ = pairPotentialSiteIdList.size();
 
-    Info<< nl << "Unique site ids found: " << allSiteIdNames_ << endl;
-
-    List<word> tetherIdList(0);
-
-    if (idListDict.found("tetherIdList"))
+    forAll(siteIdList, aSIN)
     {
-        tetherIdList = List<word>(idListDict.lookup("tetherIdList"));
+        const word& siteId = siteIdList[aSIN];
+
+        if(findIndex(pairPotentialSiteIdList, siteId) == -1)
+        {
+            pairPotentialSiteIdList.append(siteId);
+        }
+    }
+
+    siteIdList_.transfer(pairPotentialSiteIdList.shrink());
+
+    pairPotentialSiteIdList = SubList<word>(siteIdList_, nPairPotIds_);
+
+    Info<< nl << "Unique site ids found: " << siteIdList_
+        << nl << "Site Ids requiring a pair potential: "
+        << pairPotentialSiteIdList
+        << endl;
+
+    List<word> tetherSiteIdList(0);
+
+    if (idListDict.found("tetherSiteIdList"))
+    {
+        tetherSiteIdList = List<word>(idListDict.lookup("tetherSiteIdList"));
     }
 
     IOdictionary potentialDict
@@ -149,7 +187,7 @@ void Foam::potential::potential::readPotentialDict()
 
     pairPotentials_.buildPotentials
     (
-        allSiteIdNames_,
+        pairPotentialSiteIdList,
         pairDict,
         mesh_
     );
@@ -157,7 +195,7 @@ void Foam::potential::potential::readPotentialDict()
     // *************************************************************************
     // Tether potentials
 
-    if (tetherIdList.size())
+    if (tetherSiteIdList.size())
     {
         if (!potentialDict.found("tether"))
         {
@@ -170,9 +208,9 @@ void Foam::potential::potential::readPotentialDict()
 
         tetherPotentials_.buildPotentials
         (
-            allSiteIdNames_,
+            siteIdList_,
             tetherDict,
-            tetherIdList
+            tetherSiteIdList
         );
     }
 
@@ -205,7 +243,8 @@ void Foam::potential::potential::readPotentialDict()
 
 Foam::potential::potential(const polyMesh& mesh)
 :
-    mesh_(mesh)
+    mesh_(mesh),
+    electrostaticPotential_()
 {
     readPotentialDict();
 }
diff --git a/src/lagrangian/molecularDynamics/potential/potential/potential.H b/src/lagrangian/molecularDynamics/potential/potential/potential.H
index fc7430bf0c..feae0f48af 100644
--- a/src/lagrangian/molecularDynamics/potential/potential/potential.H
+++ b/src/lagrangian/molecularDynamics/potential/potential/potential.H
@@ -60,7 +60,9 @@ class potential
 
         List<word> idList_;
 
-        List<word> allSiteIdNames_;
+        List<word> siteIdList_;
+
+        label nPairPotIds_;
 
         scalar potentialEnergyLimit_;
 
@@ -68,6 +70,8 @@ class potential
 
         pairPotentialList pairPotentials_;
 
+        electrostaticPotential electrostaticPotential_;
+
         tetherPotentialList tetherPotentials_;
 
         vector gravity_;
@@ -103,7 +107,7 @@ public:
 
             inline const List<word>& idList() const;
 
-            inline const List<word>& allSiteIdNames() const;
+            inline const List<word>& siteIdList() const;
 
             inline scalar potentialEnergyLimit() const;;
 
@@ -113,6 +117,8 @@ public:
 
             inline const pairPotentialList& pairPotentials() const;
 
+            inline const electrostaticPotential& electrostatic() const;
+
             inline const tetherPotentialList& tetherPotentials() const;
 
             inline const vector& gravity() const;
diff --git a/src/lagrangian/molecularDynamics/potential/potential/potentialI.H b/src/lagrangian/molecularDynamics/potential/potential/potentialI.H
index 9e86a3a9a0..9500a56f9c 100644
--- a/src/lagrangian/molecularDynamics/potential/potential/potentialI.H
+++ b/src/lagrangian/molecularDynamics/potential/potential/potentialI.H
@@ -40,9 +40,9 @@ inline const Foam::List<Foam::word>& Foam::potential::idList() const
 }
 
 
-inline const Foam::List<Foam::word>& Foam::potential::allSiteIdNames() const
+inline const Foam::List<Foam::word>& Foam::potential::siteIdList() const
 {
-    return allSiteIdNames_;
+    return siteIdList_;
 }
 
 
@@ -64,15 +64,21 @@ inline const Foam::labelList& Foam::potential::removalOrder() const
 }
 
 
-inline const Foam::pairPotentialList&
-    Foam::potential::pairPotentials() const
+inline const Foam::pairPotentialList& Foam::potential::pairPotentials() const
 {
     return pairPotentials_;
 }
 
 
+inline const Foam::electrostaticPotential&
+Foam::potential::electrostatic() const
+{
+    return electrostaticPotential_;
+}
+
+
 inline const Foam::tetherPotentialList&
-    Foam::potential::tetherPotentials() const
+Foam::potential::tetherPotentials() const
 {
     return tetherPotentials_;
 }
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.C b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.C
index 2313f99c69..d021a18cb6 100644
--- a/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.C
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.C
@@ -30,28 +30,28 @@ License
 
 void Foam::tetherPotentialList::readTetherPotentialDict
 (
-    const List<word>& idList,
+    const List<word>& siteIdList,
     const dictionary& tetherPotentialDict,
-    const List<word>& tetherIdList
+    const List<word>& tetherSiteIdList
 )
 {
 
     Info<< nl << "Building tether potentials." << endl;
 
-    idMap_ = List<label>(idList.size(), -1);
+    idMap_ = List<label>(siteIdList.size(), -1);
 
     label tetherMapIndex = 0;
 
-    forAll(tetherIdList, t)
+    forAll(tetherSiteIdList, t)
     {
-        word tetherPotentialName = tetherIdList[t];
+        word tetherPotentialName = tetherSiteIdList[t];
 
-        label tetherId = findIndex(idList, tetherPotentialName);
+        label tetherId = findIndex(siteIdList, tetherPotentialName);
 
         if (tetherId == -1)
         {
             FatalErrorIn("tetherPotentialList::readTetherPotentialDict")
-                << nl << "No matching entry found in idList for tether name "
+                << nl << "No matching entry found in siteIdList for tether name "
                 << tetherPotentialName
                 << abort(FatalError);
         }
@@ -93,15 +93,15 @@ Foam::tetherPotentialList::tetherPotentialList()
 
 Foam::tetherPotentialList::tetherPotentialList
 (
-    const List<word>& idList,
+    const List<word>& siteIdList,
     const dictionary& tetherPotentialDict,
-    const List<word>& tetherIdList
+    const List<word>& tetherSiteIdList
 )
 :
     PtrList<tetherPotential>(),
     idMap_()
 {
-    buildPotentials(idList, tetherPotentialDict, tetherIdList);
+    buildPotentials(siteIdList, tetherPotentialDict, tetherSiteIdList);
 }
 
 
@@ -115,14 +115,14 @@ Foam::tetherPotentialList::~tetherPotentialList()
 
 void Foam::tetherPotentialList::buildPotentials
 (
-    const List<word>& idList,
+    const List<word>& siteIdList,
     const dictionary& tetherPotentialDict,
-    const List<word>& tetherIdList
+    const List<word>& tetherSiteIdList
 )
 {
-    setSize(tetherIdList.size());
+    setSize(tetherSiteIdList.size());
 
-    readTetherPotentialDict(idList, tetherPotentialDict, tetherIdList);
+    readTetherPotentialDict(siteIdList, tetherPotentialDict, tetherSiteIdList);
 }
 
 
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.H
index e45dbca9e0..cf6930d14b 100644
--- a/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.H
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.H
@@ -66,9 +66,9 @@ class tetherPotentialList
 
         void readTetherPotentialDict
         (
-            const List<word>& idList,
+            const List<word>& siteIdList,
             const dictionary& tetherPotentialDict,
-            const List<word>& tetherIdList
+            const List<word>& tetherSiteIdList
         );
 
         //- Disallow default bitwise assignment
@@ -84,12 +84,12 @@ public:
 
         tetherPotentialList();
 
-        //- Construct from idList and potental dictionaries
+        //- Construct from siteIdList and potental dictionaries
         tetherPotentialList
         (
-            const List<word>& idList,
+            const List<word>& siteIdList,
             const dictionary& tetherPotentialDict,
-            const List<word>& tetherIdList
+            const List<word>& tetherSiteIdList
         );
 
     // Destructor
@@ -100,9 +100,9 @@ public:
 
         void buildPotentials
         (
-            const List<word>& idList,
+            const List<word>& siteIdList,
             const dictionary& tetherPotentialDict,
-            const List<word>& tetherIdList
+            const List<word>& tetherSiteIdList
         );
 
         // Access
diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialListI.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialListI.H
index fe95e60c97..02066da30f 100644
--- a/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialListI.H
+++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialListI.H
@@ -39,8 +39,8 @@ inline Foam::label Foam::tetherPotentialList::tetherPotentialIndex
         (
             "Foam::tetherPotentialList::tetherPotentialIndex(const label a)"
         )
-            << "Attempting to access an undefined tetherPotential."
-            << abort(FatalError);
+        << "Attempting to access an undefined tetherPotential."
+        << abort(FatalError);
 
         return -1;
     }