Modified id reading and mapping to make mols specify which sites they need pair potentials for and only look for them in the potentialDict. Ids requiring pair potentials are stored first in the siteIdList to make the look-up of them in pairPotential work as before.

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;

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;

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

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

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);
             }
         }
 
-    allSiteIdNames_.transfer(allSiteIdNames.shrink());
+        List<word> pairPotSiteIds = molDict.lookup("pairPotentialSiteIds");
 
-    Info<< nl << "Unique site ids found: " << allSiteIdNames_ << endl;
-
-    List<word> tetherIdList(0);
-
-    if (idListDict.found("tetherIdList"))
+        forAll(pairPotSiteIds, sI)
         {
-        tetherIdList = List<word>(idListDict.lookup("tetherIdList"));
+            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);
+            }
+        }
+    }
+
+    label nPairPotIds_ = pairPotentialSiteIdList.size();
+
+    forAll(siteIdList, aSIN)
+    {
+        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();
 }

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;

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

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

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