rewriting molConfig, renamed to mdInitiaise. Driving all creation of molecules from within moleculeCloud.

applications/solvers/molecularDynamics/mdFoam/mdFoam.C

@@ -32,7 +32,6 @@ Description
 
 #include "fvCFD.H"
 #include "md.H"
-#include "potential.H"
 
 int main(int argc, char *argv[])
 {

applications/utilities/preProcessing/mdInitialise/Make/files

@@ -0,0 +1,3 @@
+mdInitialise.C
+
+EXE = $(FOAM_APPBIN)/mdInitialise

applications/utilities/preProcessing/mdInitialise/Make/options

@@ -0,0 +1,18 @@
+EXE_INC = \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools \
+    -ldynamicMesh \
+    -lfiniteVolume \
+    -llagrangian \
+    -lmolecule \
+    -lpotential \
+    -lmolecularMeasurements
+

applications/utilities/preProcessing/mdInitialise/mdInitialise.C

@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+ =========                   |
+ \\      /   F ield          | OpenFOAM: The Open Source CFD Toolbox
+  \\    /    O peration      |
+   \\  /     A nd            | Copyright (C) 1991-2008 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 "md.H"
+#include "fvCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+
+    IOdictionary mdInitialiseDict
+    (
+        IOobject
+        (
+            "mdInitialiseDict",
+            runTime.system(),
+            runTime,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        )
+    );
+
+    IOdictionary idListDict
+    (
+        IOobject
+        (
+            "idList",
+            mesh.time().constant(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        )
+    );
+
+    potential pot(mesh, mdInitialiseDict, idListDict);
+
+    moleculeCloud molecules(mesh, pot, mdInitialiseDict);
+
+    label totalMolecules = molecules.size();
+
+    if (Pstream::parRun())
+    {
+        reduce(totalMolecules, sumOp<label>());
+    }
+
+    Info<< nl << "Total number of molecules added: " << totalMolecules
+        << nl << endl;
+
+    IOstream::defaultPrecision(12);
+
+    if (!mesh.write())
+    {
+        FatalErrorIn(args.executable())
+            << "Failed writing moleculeCloud."
+            << nl << exit(FatalError);
+    }
+
+    Info<< nl << "ClockTime = " << runTime.elapsedClockTime() << " s"
+        << nl << endl;
+
+    Info << nl << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunctionI.H

@@ -25,43 +25,43 @@ License
 \*---------------------------------------------------------------------------*/
 
 template<class Type>
-inline const Field< Field<Type> >& Foam::correlationFunction<Type>::
-tZeroBuffers() const
+inline const Foam::Field< Foam::Field<Type> >&
+Foam::correlationFunction<Type>::tZeroBuffers() const
 {
     return tZeroBuffers_;
 }
 
 
 template<class Type>
-inline scalar Foam::correlationFunction<Type>::duration() const
+inline Foam::scalar Foam::correlationFunction<Type>::duration() const
 {
     return duration_;
 }
 
 
 template<class Type>
-inline scalar Foam::correlationFunction<Type>::sampleInterval() const
+inline Foam::scalar Foam::correlationFunction<Type>::sampleInterval() const
 {
     return sampleInterval_;
 }
 
 
 template<class Type>
-inline scalar Foam::correlationFunction<Type>::averagingInterval() const
+inline Foam::scalar Foam::correlationFunction<Type>::averagingInterval() const
 {
     return averagingInterval_;
 }
 
 
 template<class Type>
-inline label Foam::correlationFunction<Type>::sampleSteps() const
+inline Foam::label Foam::correlationFunction<Type>::sampleSteps() const
 {
     return sampleSteps_;
 }
 
 
 template<class Type>
-inline label Foam::correlationFunction<Type>::measurandFieldSize() const
+inline Foam::label Foam::correlationFunction<Type>::measurandFieldSize() const
 {
     return tZeroBuffers_[0].size();
 }

src/lagrangian/molecularDynamics/molecule/mdTools/md.H

@@ -1,6 +1,7 @@
 #ifndef md_H
 #define md_H
 
+#   include "potential.H"
 #   include "moleculeCloud.H"
 #   include "correlationFunction.H"
 #   include "distribution.H"

src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C

@@ -460,6 +460,53 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
 }
 
 
+void Foam::moleculeCloud::initialiseMolecules
+(
+    const IOdictionary& mdInitialiseDict
+)
+{
+    Info<< nl
+        << "Initialising molecules in each zone specified in mdInitialiseDict."
+        << endl;
+
+    const cellZoneMesh& cellZoneI = mesh_.cellZones();
+
+    if (!cellZoneI.size())
+    {
+        FatalErrorIn("void Foam::moleculeCloud::initialiseMolecules")
+            << "No cellZones found in the mesh."
+            << abort(FatalError);
+    }
+
+    forAll (cellZoneI, cZ)
+    {
+        if (cellZoneI[cZ].size())
+        {
+            if (!mdInitialiseDict.found(cellZoneI[cZ].name()))
+            {
+                Info<< "No specification subDictionary for zone "
+                    << cellZoneI[cZ].name() << " found.  Not filling." << endl;
+            }
+            else
+            {
+                label n = 0;
+
+                label totalZoneMols = 0;
+
+                label molsPlacedThisIteration = 0;
+
+                // Continue trying to place molecules as long as at
+                // least one molecule is placed in each iteration.
+                // The "|| totalZoneMols == 0" condition means that the
+                // algorithm will continue if the origin is outside the
+                // zone.
+
+
+            }
+        }
+    }
+}
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::moleculeCloud::moleculeCloud
@@ -479,77 +526,6 @@ Foam::moleculeCloud::moleculeCloud
 
     buildConstProps();
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    const Cloud<molecule>& cloud = *this;
-
-    if(Pstream::master())
-    {
-        vector position1 = vector(-2.4e-9,0,0.3e-9);
-
-        addParticle
-        (
-            new molecule
-            (
-                cloud,
-                position1,
-                mesh_.findCell(position1),
-                tensor(1, 0, 0, 0, 1, 0, 0, 0, 1),
-                vector(-405, 5, 5),
-                vector::zero,
-                vector(1e-36, -1.2e-35, -3.7e-35),
-                vector::zero,
-                vector::zero,
-                constProps(0),
-                0,
-                0
-            )
-        );
-
-        vector position2 = vector(-2.35e-9,0,-0.4e-9);
-
-        addParticle
-        (
-            new molecule
-            (
-                cloud,
-                position2,
-                mesh_.findCell(position2),
-                tensor(0, 0, 1, 0, 1, 0, -1, 0, 0),
-                vector(396,-3.7,-5),
-                vector::zero,
-                vector(-2.1e-35, 1.4e-35, 2e-36),
-                vector::zero,
-                vector::zero,
-                constProps(1),
-                0,
-                1
-            )
-        );
-
-        vector position3 = vector(-2e-9,0.52e-9,0);
-
-        addParticle
-        (
-            new molecule
-            (
-                cloud,
-                position3,
-                mesh_.findCell(position2),
-                tensor(-1, 0, 0, 0, 1, 0, 0, 0, -1),
-                vector(403.2, 1, -2),
-                vector::zero,
-                vector(4e-36, 1e-35, -1.3e-35),
-                vector::zero,
-                vector::zero,
-                constProps(1),
-                0,
-                1
-            )
-        );
-    }
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-
     buildCellOccupancy();
 
     removeHighEnergyOverlaps();
@@ -558,29 +534,27 @@ Foam::moleculeCloud::moleculeCloud
 }
 
 
-// Foam::moleculeCloud::moleculeCloud
-// (
-//     const polyMesh& mesh,
-//     label nMol,
-//     const labelField& id,
-//     const scalarField& mass,
-//     const vectorField& positions,
-//     const labelField& cells,
-//     const vectorField& U,
-//     const vectorField& A,
-//     const labelField& tethered,
-//     const vectorField& tetherPositions
-// )
-// :
-//     Cloud<molecule>(mesh, "moleculeCloud", false),
-//     mesh_(mesh),
-//     referredInteractionList_(*this)
-// {
+Foam::moleculeCloud::moleculeCloud
+(
+    const polyMesh& mesh,
+    const potential& pot,
+    const IOdictionary& mdInitialiseDict
+)
+    :
+    Cloud<molecule>(mesh, "moleculeCloud", false),
+    mesh_(mesh),
+    pot_(pot),
+    il_(mesh_),
+    constPropList_()
+{
+    buildConstProps();
+
+    initialiseMolecules(mdInitialiseDict);
+}
 //     // Do I need to read the fields if I'm just about to clear them?
 //     molecule::readFields(*this);
 
 //     clear();
-
 //     // This clear() is here for the moment to stop existing files
 //     // being appended to, this would be better accomplished by getting
 //     // mesh.removeFiles(mesh.instance()); (or equivalent) to work.

src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H

@@ -110,6 +110,11 @@ private:
 
         void removeHighEnergyOverlaps();
 
+        void initialiseMolecules
+        (
+            const IOdictionary& mdInitialiseDict
+        );
+
         //- Disallow default bitwise copy construct
         moleculeCloud(const moleculeCloud&);
 
@@ -136,23 +141,13 @@ public:
             const potential& pot
         );
 
-        //- Construct given polyMesh and fields of components.
-        //- Intended for use
-        //- by the molConfig utility
-        // moleculeCloud
-        // (
-        //     const polyMesh& mesh,
-        //     label nMol,
-        //     const labelField& id,
-        //     const scalarField& mass,
-        //     const vectorField& positions,
-        //     const labelField& cells,
-        //     const vectorField& U,
-        //     const vectorField& A,
-        //     const labelField& tethered,
-        //     const vectorField& tetherPositions
-        // );
-
+        //- Construct given mesh, potential and mdInitialiseDict
+        moleculeCloud
+        (
+            const polyMesh& mesh,
+            const potential& pot,
+            const IOdictionary& mdInitialiseDict
+        );
 
     // Member Functions
 

src/lagrangian/molecularDynamics/potential/potential/potential.C

@@ -28,37 +28,8 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::potential::potential::readPotentialDict()
+void Foam::potential::setSiteIdList(const IOdictionary& moleculePropertiesDict)
 {
-    Info<< nl <<  "Reading potential dictionary:" << endl;
-
-    IOdictionary idListDict
-    (
-        IOobject
-        (
-            "idList",
-            mesh_.time().constant(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-    idList_ = List<word>(idListDict.lookup("idList"));
-
-    IOdictionary moleculePropertiesDict
-    (
-        IOobject
-        (
-            "moleculeProperties",
-            mesh_.time().constant(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false
-        )
-    );
-
     DynamicList<word> siteIdList;
 
     DynamicList<word> pairPotentialSiteIdList;
@@ -108,7 +79,7 @@ void Foam::potential::potential::readPotentialDict()
         }
     }
 
-    label nPairPotIds_ = pairPotentialSiteIdList.size();
+    nPairPotIds_ = pairPotentialSiteIdList.size();
 
     forAll(siteIdList, aSIN)
     {
@@ -121,8 +92,46 @@ void Foam::potential::potential::readPotentialDict()
     }
 
     siteIdList_.transfer(pairPotentialSiteIdList.shrink());
+}
 
-    pairPotentialSiteIdList = SubList<word>(siteIdList_, nPairPotIds_);
+
+void Foam::potential::potential::readPotentialDict()
+{
+    Info<< nl <<  "Reading potential dictionary:" << endl;
+
+    IOdictionary idListDict
+    (
+        IOobject
+        (
+            "idList",
+            mesh_.time().constant(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    idList_ = List<word>(idListDict.lookup("idList"));
+
+    IOdictionary moleculePropertiesDict
+    (
+        IOobject
+        (
+            "moleculeProperties",
+            mesh_.time().constant(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        )
+    );
+
+    setSiteIdList(moleculePropertiesDict);
+
+    List<word> pairPotentialSiteIdList
+    (
+        SubList<word>(siteIdList_, nPairPotIds_)
+    );
 
     Info<< nl << "Unique site ids found: " << siteIdList_
         << nl << "Site Ids requiring a pair potential: "
@@ -239,6 +248,119 @@ void Foam::potential::potential::readPotentialDict()
 }
 
 
+void Foam::potential::potential::readMdInitialiseDict
+(
+    const IOdictionary& mdInitialiseDict,
+    IOdictionary& idListDict
+)
+{
+    IOdictionary moleculePropertiesDict
+    (
+        IOobject
+        (
+            "moleculeProperties",
+            mesh_.time().constant(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        )
+    );
+
+    DynamicList<word> idList;
+
+    DynamicList<word> tetherSiteIdList;
+
+    forAll(mdInitialiseDict.toc(), zone)
+    {
+        const dictionary& zoneDict = mdInitialiseDict.subDict
+        (
+            mdInitialiseDict.toc()[zone]
+        );
+
+        List<word> ids
+        (
+            zoneDict.lookup("latticeIds")
+        );
+
+        forAll(ids, i)
+        {
+            const word& id = ids[i];
+
+            if(!moleculePropertiesDict.found(id))
+            {
+                FatalErrorIn("potential.C") << nl
+                    << "Molecule type "
+                    << id
+                    << " not found in moleculeProperties dictionary."
+                    << nl
+                    << abort(FatalError);
+            }
+
+            if (findIndex(idList,id) == -1)
+            {
+                idList.append(id);
+            }
+        }
+
+        List<word> tetherSiteIds
+        (
+            zoneDict.lookup("tetherSiteIds")
+        );
+
+        forAll(tetherSiteIds, t)
+        {
+            const word& tetherSiteId = tetherSiteIds[t];
+
+            bool idFound = false;
+
+            forAll(ids, i)
+            {
+                if (idFound)
+                {
+                    break;
+                }
+
+                const word& id = ids[i];
+
+                List<word> siteIds
+                (
+                    moleculePropertiesDict.subDict(id).lookup("siteIds")
+                );
+
+                if (findIndex(siteIds, tetherSiteId) != -1)
+                {
+                    idFound = true;
+                }
+            }
+
+            if (idFound)
+            {
+                tetherSiteIdList.append(tetherSiteId);
+            }
+            else
+            {
+                FatalErrorIn("potential.C") << nl
+                    << "Tether id  "
+                    << tetherSiteId
+                    << " not found as a site of any molecule in zone."
+                    << nl
+                    << abort(FatalError);
+            }
+        }
+    }
+
+    idList_.transfer(idList);
+
+    tetherSiteIdList.shrink();
+
+    idListDict.add("idList", idList_);
+
+    idListDict.add("tetherSiteIdList", tetherSiteIdList);
+
+    setSiteIdList(moleculePropertiesDict);
+}
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::potential::potential(const polyMesh& mesh)
@@ -250,6 +372,19 @@ Foam::potential::potential(const polyMesh& mesh)
 }
 
 
+Foam::potential::potential
+(
+    const polyMesh& mesh,
+    const IOdictionary& mdInitialiseDict,
+    IOdictionary& idListDict
+)
+:
+    mesh_(mesh),
+    electrostaticPotential_()
+{
+    readMdInitialiseDict(mdInitialiseDict, idListDict);
+}
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::potential::~potential()

src/lagrangian/molecularDynamics/potential/potential/potential.H

@@ -78,8 +78,16 @@ class potential
 
     // Private Member Functions
 
+        void setSiteIdList(const IOdictionary& moleculePropertiesDict);
+
         void readPotentialDict();
 
+        void readMdInitialiseDict
+        (
+            const IOdictionary& mdInitialiseDict,
+            IOdictionary& idListDict
+        );
+
         //- Disallow default bitwise copy construct
         potential(const potential&);
 
@@ -94,6 +102,14 @@ public:
         //- Construct from mesh reference
         potential(const polyMesh& mesh);
 
+        //- Construct from mdInitialiseDict
+        potential
+        (
+            const polyMesh& mesh,
+            const IOdictionary& mdInitialiseDict,
+            IOdictionary& idListDict
+        );
+
     // Destructor
 
         ~potential();