From 69ab38bc8fd9ee0730eda36fd508fe3d9ad44acd Mon Sep 17 00:00:00 2001
From: graham <graham.macpherson@strath.ac.uk>
Date: Tue, 7 Oct 2008 20:14:10 +0100
Subject: [PATCH] Reading individual sites from moleculeProperties dictionary
 entries, WIP.

---
 .../solvers/molecularDynamics/mdFoam/mdFoam.C |  4 +-
 .../mdTransportProperitesFoam.C               |  4 +-
 .../molecule/molecule/molecule.H              | 17 ++--
 .../molecule/molecule/moleculeI.H             | 58 +++++++++++---
 .../molecule/moleculeCloud/moleculeCloud.C    | 53 +++++++++++++
 .../molecule/moleculeCloud/moleculeCloud.H    |  8 +-
 .../molecule/moleculeCloud/moleculeCloudI.H   |  2 +-
 .../potential/potential/potential.C           | 79 +++++++++++++++----
 8 files changed, 187 insertions(+), 38 deletions(-)

diff --git a/applications/solvers/molecularDynamics/mdFoam/mdFoam.C b/applications/solvers/molecularDynamics/mdFoam/mdFoam.C
index 37c5cbff9a..99a9e90d7d 100644
--- a/applications/solvers/molecularDynamics/mdFoam/mdFoam.C
+++ b/applications/solvers/molecularDynamics/mdFoam/mdFoam.C
@@ -31,7 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-//#include "md.H"
+#include "md.H"
 #include "potential.H"
 
 int main(int argc, char *argv[])
@@ -43,6 +43,8 @@ int main(int argc, char *argv[])
 
     potential pot(mesh);
 
+    moleculeCloud molecules(mesh, pot);
+
     Info << "\nStarting time loop\n" << endl;
 
     while (runTime.run())
diff --git a/applications/solvers/molecularDynamics/mdTransportProperitesFoam/mdTransportProperitesFoam.C b/applications/solvers/molecularDynamics/mdTransportProperitesFoam/mdTransportProperitesFoam.C
index 61c94d81d4..a65ae3474d 100755
--- a/applications/solvers/molecularDynamics/mdTransportProperitesFoam/mdTransportProperitesFoam.C
+++ b/applications/solvers/molecularDynamics/mdTransportProperitesFoam/mdTransportProperitesFoam.C
@@ -53,7 +53,7 @@ int main(int argc, char *argv[])
 #   include "createAutoCorrelationFunctions.H"
 
     label nAveragingSteps = 0;
-    
+
     Info << "\nStarting time loop\n" << endl;
 
     while (runTime.run())
@@ -87,6 +87,6 @@ int main(int argc, char *argv[])
 #   include "calculateTransportProperties.H"
 
     Info << "End\n" << endl;
-            
+
     return(0);
 }
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
index 95a519ef59..8fb1f5ebab 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
@@ -80,16 +80,19 @@ public:
 
         // Private data
 
-        List<vector> siteReferencePositions_;
+            List<vector> siteReferencePositions_;
 
-        List<scalar> siteCharges_;
+            List<scalar> siteCharges_;
 
-        List<label> siteIds_;
+            List<label> siteIds_;
 
-        diagTensor momentOfInertia_;
+            diagTensor momentOfInertia_;
 
-        scalar mass_;
+            scalar mass_;
 
+        // Private Member Functions
+
+            void checkSiteListSizes() const;
 
     public:
 
@@ -105,6 +108,9 @@ public:
             scalar mass
         );
 
+        //- Constructor
+        inline constantProperties(const dictionary& dict);
+
         // Member functions
 
         inline const List<vector>& siteReferencePositions() const;
@@ -112,6 +118,7 @@ public:
         inline const List<scalar>& siteCharges() const;
 
         inline const List<label>& siteIds() const;
+        inline List<label>& siteIds();
 
         inline const diagTensor& momentOfInertia() const;
 
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
index b6f3438aca..8da0c8b676 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H
@@ -33,9 +33,7 @@ inline Foam::molecule::constantProperties::constantProperties()
     siteIds_(List<label>(0)),
     momentOfInertia_(diagTensor(0, 0, 0)),
     mass_(0)
-{
-
-}
+{}
 
 
 inline Foam::molecule::constantProperties::constantProperties
@@ -53,16 +51,24 @@ inline Foam::molecule::constantProperties::constantProperties
     momentOfInertia_(momentOfInertia),
     mass_(mass)
 {
-    if
-    (
-        siteIds_.size() != siteReferencePositions_.size()
-     || siteIds_.size() != siteCharges_.size()
-    )
-    {
-        FatalErrorIn("molecule::constantProperties::constantProperties")
-            << "Sizes of site id, charge and referencePositions are not the same. "
-            << abort(FatalError);
-    }
+    checkSiteListSizes();
+}
+
+
+inline Foam::molecule::constantProperties::constantProperties
+(
+    const dictionary& dict
+)
+:
+    siteReferencePositions_(dict.lookup("siteReferencePositions")),
+    siteCharges_(dict.lookup("siteCharges")),
+    siteIds_(List<word>(dict.lookup("siteIds")).size(), -1),
+    momentOfInertia_(dict.lookup("momentOfInertia")),
+    mass_(readScalar(dict.lookup("mass")))
+{
+    checkSiteListSizes();
+
+    Info<< nl << "siteIds not set yet." << endl;
 }
 
 
@@ -101,6 +107,25 @@ inline Foam::molecule::molecule
 }
 
 
+// * * * constantProperties Private Member Functions * * * * * * * * * * * * //
+
+
+inline void Foam::molecule::constantProperties::checkSiteListSizes() const
+{
+    if
+    (
+        siteIds_.size() != siteReferencePositions_.size()
+     || siteIds_.size() != siteCharges_.size()
+    )
+    {
+        FatalErrorIn("molecule::constantProperties::checkSiteListSizes")
+            << "Sizes of site id, charge and "
+            << "referencePositions are not the same. "
+            << nl << abort(FatalError);
+    }
+}
+
+
 // * * * * * * * constantProperties Member Functions * * * * * * * * * * * * //
 
 inline const Foam::List<Foam::vector>&
@@ -124,6 +149,13 @@ inline const Foam::List<Foam::label>&
 }
 
 
+inline Foam::List<Foam::label>&
+    Foam::molecule::constantProperties::siteIds()
+{
+    return siteIds_;
+}
+
+
 inline const Foam::diagTensor&
     Foam::molecule::constantProperties::momentOfInertia() const
 {
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
index 42555d8c5c..02473e90a5 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
@@ -44,6 +44,57 @@ Foam::scalar Foam::moleculeCloud::vacuumPermittivity = 8.854187817e-12;
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+void Foam::moleculeCloud::buildConstProps()
+{
+    Info<< nl << "Reading moleculeProperties dictionary." << endl;
+
+    const List<word>& idList(pot_.idList());
+
+    constPropList_.setSize(idList.size());
+
+    const List<word>& allSiteIdNames(pot_.allSiteIdNames());
+
+    forAll(idList, i)
+    {
+        const word& id(idList[i]);
+
+        const dictionary& molDict(moleculePropertiesDict.subDict(id));
+
+        List<word> siteIdNames = molDict.lookup("siteIds");
+
+        List<label> siteIds(siteIdNames.size());
+
+        forAll(siteIdNames, sI)
+        {
+            const word& siteId = siteIdNames[sI];
+
+            siteIds[sI] = findIndex(allSiteIdNames, siteId);
+
+            if (siteIds[sI] == -1)
+            {
+                FatalErrorIn("moleculeCloud.C") << nl
+                    << siteId << " site not found."
+                    << nl << abort(FatalError);
+            }
+        }
+
+        molecule::constantProperties& constProp = constPropList_[i];
+
+        constProp = molecule::constantProperties(molDict);
+
+        constProp.siteIds() = siteIds;
+
+        Info<< constPropList_[i].siteReferencePositions()
+            << nl << constPropList_[i].siteCharges()
+            << nl << constPropList_[i].siteIds()
+            << nl << constPropList_[i].momentOfInertia()
+            << nl << constPropList_[i].mass()
+            << nl << constPropList_[i].nSites()
+            << endl;
+    }
+}
+
+
 void Foam::moleculeCloud::buildCellOccupancy()
 {
     forAll(cellOccupancy_, cO)
@@ -422,6 +473,8 @@ Foam::moleculeCloud::moleculeCloud
 {
     molecule::readFields(*this);
 
+    buildConstProps();
+
     buildCellOccupancy();
 
     removeHighEnergyOverlaps();
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
index ebb9f00906..3e54f65e46 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
@@ -73,6 +73,8 @@ private:
 
     // Private Member Functions
 
+        void buildConstProps();
+
         //- Determine which molecules are in which cells
         void buildCellOccupancy();
 
@@ -141,7 +143,8 @@ public:
             const potential& pot
         );
 
-        //- Construct given polyMesh and fields of components.  Intended for use
+        //- Construct given polyMesh and fields of components.
+        //- Intended for use
         //- by the molConfig utility
         // moleculeCloud
         // (
@@ -184,7 +187,8 @@ public:
 
             inline const List<molecule::constantProperties> constProps() const;
 
-            inline const molecule::constantProperties& constProps(label id) const;
+            inline const molecule::constantProperties&
+            constProps(label id) const;
 
     // Member Operators
 
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
index 76fae4db46..17145089cd 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
@@ -101,7 +101,7 @@ inline void Foam::moleculeCloud::evaluatePair
         // Acceleration increment for molReal
         molReal->a() += fRealRef/massReal;
 
-        // Adding a contribution 19f 1/2 of the potential energy
+        // Adding a contribution of 1/2 of the potential energy
         // from this interaction
         molReal->potentialEnergy() +=
             0.5*pairPot.energy(idReal, idRef, rRealRefMag);
diff --git a/src/lagrangian/molecularDynamics/potential/potential/potential.C b/src/lagrangian/molecularDynamics/potential/potential/potential.C
index c358d6b376..5593b12e37 100644
--- a/src/lagrangian/molecularDynamics/potential/potential/potential.C
+++ b/src/lagrangian/molecularDynamics/potential/potential/potential.C
@@ -31,7 +31,22 @@ License
 Foam::potential::potential(const polyMesh& mesh)
 :
     mesh_(mesh)
+{}
 
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::potential::~potential()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+void Foam::potential::potential::readPotentialDict
+(
+    const List<word>& allSiteIdNames
+)
 {
     Info<< nl <<  "Reading potential dictionary:" << endl;
 
@@ -49,6 +64,52 @@ Foam::potential::potential(const polyMesh& mesh)
 
     idList_ = List<word>(idListDict.lookup("idList"));
 
+    IOdictionary moleculePropertiesDict
+    (
+        IOobject
+        (
+            "moleculeProperties",
+            mesh_.time().constant(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        )
+    );
+
+    DynamicList<word> allSiteIdNames;
+
+    forAll(idList, i)
+    {
+        const word& id(idList[i]);
+
+        if (!moleculePropertiesDict.found(id))
+        {
+            FatalErrorIn("potential.C") << nl
+                << id << " molecule subDict not found"
+                << nl << abort(FatalError);
+        }
+
+        const dictionary& molDict(moleculePropertiesDict.subDict(id));
+
+        List<word> siteIdNames = molDict.lookup("siteIds");
+
+        forAll(siteIdNames, sI)
+        {
+            const word& siteId = siteIdNames[sI];
+
+            if(findIndex(allSiteIdNames, siteId) == -1)
+            {
+                allSiteIdNames.append(siteId);
+            }
+        }
+    }
+
+    // CREATE allSiteIdNames_ DATA MEMBER
+    allSiteIdNames_.transfer(allSiteIdNames.shrink());
+
+    Info<< nl << "Unique site ids found: " << allSiteIdNames << endl;
+
     List<word> tetherIdList(0);
 
     if (idListDict.found("tetherIdList"))
@@ -85,7 +146,7 @@ Foam::potential::potential(const polyMesh& mesh)
         }
     }
 
-    // ****************************************************************************
+    // *************************************************************************
     // Pair potentials
 
     if (!potentialDict.found("pair"))
@@ -99,7 +160,7 @@ Foam::potential::potential(const polyMesh& mesh)
 
     pairPotentials_.buildPotentials(idList_, pairDict, mesh_);
 
-    // ****************************************************************************
+    // *************************************************************************
     // Tether potentials
 
     if (tetherIdList.size())
@@ -116,7 +177,7 @@ Foam::potential::potential(const polyMesh& mesh)
         tetherPotentials_.buildPotentials(idList_, tetherDict, tetherIdList);
     }
 
-    // ****************************************************************************
+    // *************************************************************************
     // External Forces
 
     gravity_ = vector::zero;
@@ -128,7 +189,7 @@ Foam::potential::potential(const polyMesh& mesh)
 
         const dictionary& externalDict = potentialDict.subDict("external");
 
-        // ************************************************************************
+        // *********************************************************************
         // gravity
 
         if (externalDict.found("gravity"))
@@ -140,16 +201,6 @@ Foam::potential::potential(const polyMesh& mesh)
     Info << nl << tab << "gravity = " << gravity_ << endl;
 }
 
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::potential::~potential()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 // ************************************************************************* //