From c85269854f1490e160e176c6ee81d961874954ce Mon Sep 17 00:00:00 2001
From: graham <graham.macpherson@strath.ac.uk>
Date: Tue, 18 Nov 2008 17:51:53 +0000
Subject: [PATCH] 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.

---
 .../solvers/molecularDynamics/mdFoam/mdFoam.C |  4 -
 .../molecule/molecule/molecule.H              | 26 ++++---
 .../molecule/molecule/moleculeI.H             | 76 +++++++++++++------
 .../molecule/moleculeCloud/moleculeCloud.C    |  6 +-
 .../potential/potential/potential.C           | 69 +++++++++++++----
 .../potential/potential/potential.H           | 10 ++-
 .../potential/potential/potentialI.H          | 16 ++--
 .../tetherPotentialList/tetherPotentialList.C | 28 +++----
 .../tetherPotentialList/tetherPotentialList.H | 14 ++--
 .../tetherPotentialListI.H                    |  4 +-
 10 files changed, 167 insertions(+), 86 deletions(-)

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