ENH: Templated MolecularCloud.

applications/solvers/discreteMethods/molecularDynamics/mdFoam/createFields.H

@@ -12,6 +12,23 @@
         mesh
     );
 
-    potential pot(mesh);
+    word cloudName("polyatomicCloud");
 
-    polyatomicCloud molecules(mesh, pot);
+    potential pot
+    (
+        mesh,
+        IOdictionary
+        (
+            IOobject
+            (
+                cloudName +  "Properties",
+                mesh.time().constant(),
+                mesh,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE,
+                false
+            )
+        )
+    );
+
+    polyatomicCloud molecules(cloudName, mesh, pot);

applications/utilities/preProcessing/mdInitialise/Make/options

@@ -3,7 +3,6 @@ EXE_INC = \
     -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
 
@@ -13,6 +12,5 @@ EXE_LIBS = \
     -lfiniteVolume \
     -llagrangian \
     -lmolecule \
-    -lpotential \
+    -lpotential
-    -lmolecularMeasurements
 

applications/utilities/preProcessing/mdInitialise/mdInitialise.C

@@ -62,11 +62,28 @@ int main(int argc, char *argv[])
         )
     );
 
-    potential pot(mesh, mdInitialiseDict, idListDict);
+    word cloudName("polyatomicCloud");
 
-    polyatomicCloud molecules(mesh, pot, mdInitialiseDict);
+    potential pot
+    (
+        mesh,
+        mdInitialiseDict,
+        IOdictionary
+        (
+            IOobject
+            (
+                cloudName +  "Properties",
+                mesh.time().constant(),
+                mesh,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE,
+                false
+            )
+        ),
+        idListDict
+    );
 
-    Info<< "Cloud is called " << molecules.name() << endl;
+    polyatomicCloud molecules(cloudName, mesh, pot, mdInitialiseDict);
 
     label totalMolecules = molecules.size();
 

src/lagrangian/molecularDynamics/molecule/Make/files

@@ -1,9 +1,9 @@
-constPropSite/constPropSite.C
+clouds/baseClasses/moleculeCloud/moleculeCloud.C
 
-polyatomic/polyatomic/polyatomic.C
+molecules/constPropSite/constPropSite.C
-polyatomic/polyatomic/polyatomicIO.C
 
-polyatomic/polyatomicCloud/polyatomicCloud.C
+molecules/polyatomic/polyatomic.C
+molecules/polyatomic/polyatomicIO.C
 
 /* controllers/basic/controllers/controllers.C
 controllers/basic/stateController/stateController.C

src/lagrangian/molecularDynamics/molecule/Make/options

@@ -2,13 +2,11 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
-    -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude
-    -I$(LIB_SRC)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
     -llagrangian \
-    -lpotential \
+    -lpotential
-    -lmolecularMeasurements
 

src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.C

@@ -23,36 +23,29 @@ License
 
 \*----------------------------------------------------------------------------*/
 
-#include "polyatomicCloud.H"
+#include "MoleculeCloud.H"
 #include "fvMesh.H"
 #include "mathematicalConstants.H"
 
 using namespace Foam::constant::mathematical;
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTemplateTypeNameAndDebug(Cloud<polyatomic>, 0);
-}
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::polyatomicCloud::buildConstProps()
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::buildConstProps()
 {
-    Info<< nl << "Reading mdProperties dictionary." << endl;
-
     const List<word>& idList(pot_.idList());
 
     constPropList_.setSize(idList.size());
 
     const List<word>& siteIdList(pot_.siteIdList());
 
-    IOdictionary mdPropertiesDict
+    IOdictionary molPropertiesDict
     (
         IOobject
         (
-            "mdProperties",
+            this->name() + "Properties",
             mesh_.time().constant(),
             mesh_,
             IOobject::MUST_READ_IF_MODIFIED,
@@ -61,11 +54,13 @@ void Foam::polyatomicCloud::buildConstProps()
         )
     );
 
+    Info<< nl << "Reading " << molPropertiesDict.name() << endl;
+
     forAll(idList, i)
     {
         const word& id = idList[i];
 
-        const dictionary& molDict(mdPropertiesDict.subDict(id));
+        const dictionary& molDict(molPropertiesDict.subDict(id));
 
         List<word> siteIdNames = molDict.lookup("siteIds");
 
@@ -81,17 +76,14 @@ void Foam::polyatomicCloud::buildConstProps()
             {
                 FatalErrorIn
                 (
-                    "Foam::polyatomic::constantProperties::constantProperties"
+                    "void Foam::MoleculeCloud<MoleculeType>::buildConstProps()"
-                    "("
-                    "const dictionary& dict"
-                    ")"
                 )
                     << siteId << " site not found."
                     << nl << abort(FatalError);
             }
         }
 
-        constPropList_[i] = polyatomic::constantProperties
+        constPropList_[i] = typename MoleculeType::constantProperties
         (
             molDict,
             siteIds
@@ -100,11 +92,15 @@ void Foam::polyatomicCloud::buildConstProps()
 }
 
 
-void Foam::polyatomicCloud::setSiteSizesAndPositions()
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::setSiteSizesAndPositions()
 {
-    forAllIter(polyatomicCloud, *this, mol)
+    forAllIter(typename MoleculeCloud<MoleculeType>, *this, mol)
     {
-        const polyatomic::constantProperties& cP = constProps(mol().id());
+        const typename MoleculeType::constantProperties& cP
+        (
+            constProps(mol().id())
+        );
 
         mol().setSiteSizes(cP.nSites());
 
@@ -113,14 +109,15 @@ void Foam::polyatomicCloud::setSiteSizesAndPositions()
 }
 
 
-void Foam::polyatomicCloud::buildCellOccupancy()
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::buildCellOccupancy()
 {
     forAll(cellOccupancy_, cO)
     {
         cellOccupancy_[cO].clear();
     }
 
-    forAllIter(polyatomicCloud, *this, mol)
+    forAllIter(typename MoleculeCloud<MoleculeType>, *this, mol)
     {
         cellOccupancy_[mol().cell()].append(&mol());
     }
@@ -132,15 +129,16 @@ void Foam::polyatomicCloud::buildCellOccupancy()
 }
 
 
-void Foam::polyatomicCloud::calculatePairForce()
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::calculatePairForce()
 {
     PstreamBuffers pBufs(Pstream::nonBlocking);
 
     // Start sending referred data
     il_.sendReferredData(cellOccupancy(), pBufs);
 
-    polyatomic* molI = NULL;
+    MoleculeType* molI = NULL;
-    polyatomic* molJ = NULL;
+    MoleculeType* molJ = NULL;
 
     {
         // Real-Real interactions
@@ -155,7 +153,7 @@ void Foam::polyatomicCloud::calculatePairForce()
 
                 forAll(dil[d], interactingCells)
                 {
-                    List<polyatomic*> cellJ =
+                    List<MoleculeType*> cellJ =
                         cellOccupancy_[dil[d][interactingCells]];
 
                     forAll(cellJ, cellJMols)
@@ -187,24 +185,24 @@ void Foam::polyatomicCloud::calculatePairForce()
 
         const labelListList& ril = il_.ril();
 
-        List<IDLList<polyatomic> >& referredMols = il_.referredParticles();
+        List<IDLList<MoleculeType> >& referredMols = il_.referredParticles();
 
         forAll(ril, r)
         {
             const List<label>& realCells = ril[r];
 
-            IDLList<polyatomic>& refMols = referredMols[r];
+            IDLList<MoleculeType>& refMols = referredMols[r];
 
             forAllIter
             (
-                IDLList<polyatomic>,
+                typename IDLList<MoleculeType>,
                 refMols,
                 refMol
             )
             {
                 forAll(realCells, rC)
                 {
-                    List<polyatomic*> cellI = cellOccupancy_[realCells[rC]];
+                    List<MoleculeType*> cellI = cellOccupancy_[realCells[rC]];
 
                     forAll(cellI, cellIMols)
                     {
@@ -219,11 +217,12 @@ void Foam::polyatomicCloud::calculatePairForce()
 }
 
 
-void Foam::polyatomicCloud::calculateTetherForce()
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::calculateTetherForce()
 {
     const tetherPotentialList& tetherPot(pot_.tetherPotentials());
 
-    forAllIter(polyatomicCloud, *this, mol)
+    forAllIter(typename MoleculeCloud<MoleculeType>, *this, mol)
     {
         if (mol().tethered())
         {
@@ -246,16 +245,18 @@ void Foam::polyatomicCloud::calculateTetherForce()
 }
 
 
-void Foam::polyatomicCloud::calculateExternalForce()
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::calculateExternalForce()
 {
-    forAllIter(polyatomicCloud, *this, mol)
+    forAllIter(typename MoleculeCloud<MoleculeType>, *this, mol)
     {
         mol().a() += pot_.gravity();
     }
 }
 
 
-void Foam::polyatomicCloud::removeHighEnergyOverlaps()
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::removeHighEnergyOverlaps()
 {
     Info<< nl << "Removing high energy overlaps, limit = "
         << pot_.potentialEnergyLimit()
@@ -274,11 +275,11 @@ void Foam::polyatomicCloud::removeHighEnergyOverlaps()
 
     // Real-Real interaction
 
-    polyatomic* molI = NULL;
+    MoleculeType* molI = NULL;
-    polyatomic* molJ = NULL;
+    MoleculeType* molJ = NULL;
 
     {
-        DynamicList<polyatomic*> molsToDelete;
+        DynamicList<MoleculeType*> molsToDelete;
 
         const labelListList& dil(il_.dil());
 
@@ -290,7 +291,7 @@ void Foam::polyatomicCloud::removeHighEnergyOverlaps()
 
                 forAll(dil[d], interactingCells)
                 {
-                    List<polyatomic*> cellJ =
+                    List<MoleculeType*> cellJ =
                         cellOccupancy_[dil[d][interactingCells]];
 
                     forAll(cellJ, cellJMols)
@@ -376,19 +377,19 @@ void Foam::polyatomicCloud::removeHighEnergyOverlaps()
     // Real-Referred interaction
 
     {
-        DynamicList<polyatomic*> molsToDelete;
+        DynamicList<MoleculeType*> molsToDelete;
 
         const labelListList& ril(il_.ril());
 
-        List<IDLList<polyatomic> >& referredMols = il_.referredParticles();
+        List<IDLList<MoleculeType> >& referredMols = il_.referredParticles();
 
         forAll(ril, r)
         {
-            IDLList<polyatomic>& refMols = referredMols[r];
+            IDLList<MoleculeType>& refMols = referredMols[r];
 
             forAllIter
             (
-                IDLList<polyatomic>,
+                typename IDLList<MoleculeType>,
                 refMols,
                 refMol
             )
@@ -401,7 +402,7 @@ void Foam::polyatomicCloud::removeHighEnergyOverlaps()
                 {
                     label cellI = realCells[rC];
 
-                    List<polyatomic*> cellIMols = cellOccupancy_[cellI];
+                    List<MoleculeType*> cellIMols = cellOccupancy_[cellI];
 
                     forAll(cellIMols, cIM)
                     {
@@ -430,7 +431,7 @@ void Foam::polyatomicCloud::removeHighEnergyOverlaps()
                              == findIndex(pot_.removalOrder(), idJ)
                             )
                             {
-                                // Remove one of the polyatomics
+                                // Remove one of the molecules
                                 // arbitrarily, assuring that a
                                 // consistent decision is made for
                                 // both real-referred pairs.
@@ -470,17 +471,18 @@ void Foam::polyatomicCloud::removeHighEnergyOverlaps()
         reduce(molsRemoved, sumOp<label>());
     }
 
-    Info<< tab << molsRemoved << " polyatomics removed" << endl;
+    Info<< tab << molsRemoved << " molecules removed" << endl;
 }
 
 
-void Foam::polyatomicCloud::initialiseMolecules
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::initialiseMolecules
 (
     const IOdictionary& mdInitialiseDict
 )
 {
     Info<< nl
-        << "Initialising polyatomics in each zone specified in "
+        << "Initialising molecules in each zone specified in "
         << mdInitialiseDict.name()
         << endl;
 
@@ -488,7 +490,10 @@ void Foam::polyatomicCloud::initialiseMolecules
 
     if (!cellZones.size())
     {
-        FatalErrorIn("void Foam::polyatomicCloud::initialiseMolecules")
+        FatalErrorIn
+        (
+            "void Foam::MoleculeCloud<MoleculeType>::initialiseMolecules"
+        )
             << "No cellZones found in the mesh."
             << abort(FatalError);
     }
@@ -528,7 +533,10 @@ void Foam::polyatomicCloud::initialiseMolecules
 
                 if (latticeIds.size() != latticePositions.size())
                 {
-                    FatalErrorIn("Foam::polyatomicCloud::initialiseMolecules")
+                    FatalErrorIn
+                    (
+                        "Foam::MoleculeCloud<MoleculeType>::initialiseMolecules"
+                    )
                         << "latticeIds and latticePositions must be the same "
                         << " size." << nl
                         << abort(FatalError);
@@ -550,7 +558,10 @@ void Foam::polyatomicCloud::initialiseMolecules
 
                     if (numberDensity < VSMALL)
                     {
-                        WarningIn("polyatomicCloud::initialiseMolecules")
+                        WarningIn
+                        (
+                            "MoleculeCloud<MoleculeType>::initialiseMolecules"
+                        )
                             << "numberDensity too small, not filling zone "
                             << zone.name() << endl;
 
@@ -571,7 +582,7 @@ void Foam::polyatomicCloud::initialiseMolecules
                     {
                         label id = findIndex(pot_.idList(), latticeIds[i]);
 
-                        const polyatomic::constantProperties& cP
+                        const typename MoleculeType::constantProperties& cP
                         (
                             constProps(id)
                         );
@@ -586,7 +597,10 @@ void Foam::polyatomicCloud::initialiseMolecules
 
                     if (massDensity < VSMALL)
                     {
-                        WarningIn("polyatomicCloud::initialiseMolecules")
+                        WarningIn
+                        (
+                            "MoleculeCloud<MoleculeType>::initialiseMolecules"
+                        )
                             << "massDensity too small, not filling zone "
                             << zone.name() << endl;
 
@@ -602,7 +616,10 @@ void Foam::polyatomicCloud::initialiseMolecules
                 }
                 else
                 {
-                    FatalErrorIn("Foam::polyatomicCloud::initialiseMolecules")
+                    FatalErrorIn
+                    (
+                        "Foam::MoleculeCloud<MoleculeType>::initialiseMolecules"
+                    )
                         << "massDensity or numberDensity not specified " << nl
                         << abort(FatalError);
                 }
@@ -697,8 +714,8 @@ void Foam::polyatomicCloud::initialiseMolecules
 
                 anchor += (R & (latticeCellShape & latticeAnchor));
 
-                // Continue trying to place polyatomics as long as at
+                // Continue trying to place molecule as long as at
-                // least one polyatomic is placed in each iteration.
+                // 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.
@@ -719,10 +736,6 @@ void Foam::polyatomicCloud::initialiseMolecules
 
                     bool partOfLayerInBounds = false;
 
-                    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-                    // start of placement of polyatomics
-                    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
                     if (n == 0)
                     {
                         // Special treatment is required for the first position,
@@ -959,10 +972,6 @@ void Foam::polyatomicCloud::initialiseMolecules
                         }
                     }
 
-                    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-                    // end of placement of polyatomics
-                    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
                     if
                     (
                         totalZoneMols == 0
@@ -971,11 +980,12 @@ void Foam::polyatomicCloud::initialiseMolecules
                     {
                         WarningIn
                         (
-                            "Foam::polyatomicCloud::initialiseMolecules()"
+                            "Foam::MoleculeCloud<MoleculeType>::"
+                            "initialiseMolecules()"
                         )
                             << "A whole layer of unit cells was placed "
                             << "outside the bounds of the mesh, but no "
-                            << "polyatomics have been placed in zone '"
+                            << "molecules have been placed in zone '"
                             << zone.name()
                             << "'.  This is likely to be because the zone "
                             << "has few cells ("
@@ -1001,7 +1011,8 @@ void Foam::polyatomicCloud::initialiseMolecules
 }
 
 
-void Foam::polyatomicCloud::createMolecule
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::createMolecule
 (
     const point& position,
     label cell,
@@ -1020,7 +1031,7 @@ void Foam::polyatomicCloud::createMolecule
 
     if (cell == -1)
     {
-        FatalErrorIn("Foam::polyatomicCloud::createMolecule")
+        FatalErrorIn("Foam::MoleculeCloud<MoleculeType>::createMolecule")
             << "Position specified does not correspond to a mesh cell." << nl
             << abort(FatalError);
     }
@@ -1033,10 +1044,10 @@ void Foam::polyatomicCloud::createMolecule
     {
         specialPosition = position;
 
-        special = polyatomic::SPECIAL_TETHERED;
+        special = MoleculeType::SPECIAL_TETHERED;
     }
 
-    const polyatomic::constantProperties& cP(constProps(id));
+    const typename MoleculeType::constantProperties& cP(constProps(id));
 
     vector v = equipartitionLinearVelocity(temperature, cP.mass());
 
@@ -1072,7 +1083,7 @@ void Foam::polyatomicCloud::createMolecule
 
     addParticle
     (
-        new polyatomic
+        new MoleculeType
         (
             mesh_,
             position,
@@ -1093,11 +1104,12 @@ void Foam::polyatomicCloud::createMolecule
 }
 
 
-Foam::label Foam::polyatomicCloud::nSites() const
+template<class MoleculeType>
+Foam::label Foam::MoleculeCloud<MoleculeType>::nSites() const
 {
     label n = 0;
 
-    forAllConstIter(polyatomicCloud, *this, mol)
+    forAllConstIter(typename MoleculeCloud<MoleculeType>, *this, mol)
     {
         n += constProps(mol().id()).nSites();
     }
@@ -1108,14 +1120,17 @@ Foam::label Foam::polyatomicCloud::nSites() const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::polyatomicCloud::polyatomicCloud
+template<class MoleculeType>
+Foam::MoleculeCloud<MoleculeType>::MoleculeCloud
 (
+    const word& cloudName,
     const polyMesh& mesh,
     const potential& pot,
     bool readFields
 )
 :
-    Cloud<polyatomic>(mesh, "polyatomicCloud", false),
+    Cloud<MoleculeType>(mesh, cloudName, false),
+    moleculeCloud(),
     mesh_(mesh),
     pot_(pot),
     cellOccupancy_(mesh_.nCells()),
@@ -1125,7 +1140,7 @@ Foam::polyatomicCloud::polyatomicCloud
 {
     if (readFields)
     {
-        polyatomic::readFields(*this);
+        MoleculeType::readFields(*this);
     }
 
     buildConstProps();
@@ -1138,15 +1153,18 @@ Foam::polyatomicCloud::polyatomicCloud
 }
 
 
-Foam::polyatomicCloud::polyatomicCloud
+template<class MoleculeType>
+Foam::MoleculeCloud<MoleculeType>::MoleculeCloud
 (
+    const word& cloudName,
     const polyMesh& mesh,
     const potential& pot,
     const IOdictionary& mdInitialiseDict,
     bool readFields
 )
 :
-    Cloud<polyatomic>(mesh, "polyatomicCloud", false),
+    Cloud<MoleculeType>(mesh, cloudName, false),
+    moleculeCloud(),
     mesh_(mesh),
     pot_(pot),
     il_(mesh_, 0.0, false),
@@ -1155,10 +1173,10 @@ Foam::polyatomicCloud::polyatomicCloud
 {
     if (readFields)
     {
-        polyatomic::readFields(*this);
+        MoleculeType::readFields(*this);
     }
 
-    clear();
+    this->clear();
 
     buildConstProps();
 
@@ -1168,32 +1186,34 @@ Foam::polyatomicCloud::polyatomicCloud
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::polyatomicCloud::evolve()
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::evolve()
 {
-    polyatomic::trackingData td0(*this, 0);
+    typename MoleculeType::trackingData td0(*this, 0);
-    Cloud<polyatomic>::move(td0, mesh_.time().deltaTValue());
+    Cloud<MoleculeType>::move(td0, mesh_.time().deltaTValue());
 
-    polyatomic::trackingData td1(*this, 1);
+    typename MoleculeType::trackingData td1(*this, 1);
-    Cloud<polyatomic>::move(td1, mesh_.time().deltaTValue());
+    Cloud<MoleculeType>::move(td1, mesh_.time().deltaTValue());
 
-    polyatomic::trackingData td2(*this, 2);
+    typename MoleculeType::trackingData td2(*this, 2);
-    Cloud<polyatomic>::move(td2, mesh_.time().deltaTValue());
+    Cloud<MoleculeType>::move(td2, mesh_.time().deltaTValue());
 
     calculateForce();
 
-    polyatomic::trackingData td3(*this, 3);
+    typename MoleculeType::trackingData td3(*this, 3);
-    Cloud<polyatomic>::move(td3, mesh_.time().deltaTValue());
+    Cloud<MoleculeType>::move(td3, mesh_.time().deltaTValue());
 
     info();
 }
 
 
-void Foam::polyatomicCloud::calculateForce()
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::calculateForce()
 {
     buildCellOccupancy();
 
     // Set accumulated quantities to zero
-    forAllIter(polyatomicCloud, *this, mol)
+    forAllIter(typename MoleculeCloud<MoleculeType>, *this, mol)
     {
         mol().siteForces() = vector::zero;
 
@@ -1210,36 +1230,8 @@ void Foam::polyatomicCloud::calculateForce()
 }
 
 
-void Foam::polyatomicCloud::applyConstraintsAndThermostats
+template<class MoleculeType>
-(
+void Foam::MoleculeCloud<MoleculeType>::info() const
-    const scalar targetTemperature,
-    const scalar measuredTemperature
-)
-{
-    scalar temperatureCorrectionFactor =
-        sqrt(targetTemperature/max(VSMALL, measuredTemperature));
-
-    Info<< "----------------------------------------" << nl
-        << "Temperature equilibration" << nl
-        << "Target temperature = "
-        << targetTemperature << nl
-        << "Measured temperature = "
-        << measuredTemperature << nl
-        << "Temperature correction factor = "
-        << temperatureCorrectionFactor << nl
-        << "----------------------------------------"
-        << endl;
-
-    forAllIter(polyatomicCloud, *this, mol)
-    {
-        mol().v() *= temperatureCorrectionFactor;
-
-        mol().pi() *= temperatureCorrectionFactor;
-    }
-}
-
-
-void Foam::polyatomicCloud::info() const
 {
     // Calculates and prints the mean momentum and energy in the system
     // and the number of molecules.
@@ -1267,7 +1259,7 @@ void Foam::polyatomicCloud::info() const
     label dofs = 0;
 
     {
-        forAllConstIter(polyatomicCloud, *this, mol)
+        forAllConstIter(typename MoleculeCloud<MoleculeType>, *this, mol)
         {
             const label molId = mol().id();
 
@@ -1283,11 +1275,14 @@ void Foam::polyatomicCloud::info() const
         //     CentreOfMass /= totalMass;
         // }
 
-        forAllConstIter(polyatomicCloud, *this, mol)
+        forAllConstIter(typename MoleculeCloud<MoleculeType>, *this, mol)
         {
             const label molId = mol().id();
 
-            const polyatomic::constantProperties cP(this->constProps(molId));
+            const typename MoleculeType::constantProperties cP
+            (
+                this->constProps(molId)
+            );
 
             scalar molMass(cP.mass());
 
@@ -1340,27 +1335,27 @@ void Foam::polyatomicCloud::info() const
 
     if (nMols)
     {
-        Info<< "Number of molecules in " << this->name() << " = "
+        Info<< nl << "Number of molecules in " << this->name() << " = "
             << nMols << nl
-            << "Overall number density = "
+            << "    Overall number density = "
             << nMols/meshVolume << nl
-            << "Overall mass density = "
+            << "    Overall mass density = "
             << totalMass/meshVolume << nl
-            << "Average linear momentum per molecule = "
+            << "    Average linear momentum per molecule = "
             << totalLinearMomentum/nMols << ' '
             << mag(totalLinearMomentum)/nMols << nl
-            << "Average angular momentum per molecule = "
+            << "    Average angular momentum per molecule = "
             << totalAngularMomentum << ' '
             << mag(totalAngularMomentum)/nMols << nl
-            << "maximum |velocity| = "
+            << "    maximum |velocity| = "
             << maxVelocityMag << nl
-            << "Average linear KE per molecule = "
+            << "    Average linear KE per molecule = "
             << totalLinearKE/nMols << nl
-            << "Average angular KE per molecule = "
+            << "    Average angular KE per molecule = "
             << totalAngularKE/nMols << nl
-            << "Average PE per molecule = "
+            << "    Average PE per molecule = "
             << totalPE/nMols << nl
-            << "Average TE per molecule = "
+            << "    Average TE per molecule = "
             <<
             (
                   totalLinearKE
@@ -1368,25 +1363,29 @@ void Foam::polyatomicCloud::info() const
                 + totalPE
             )
             /nMols
-            << endl;
+            << nl << endl;
     }
     else
     {
-        Info<< "No molecules in " << this->name() << endl;
+        Info<< nl << "No molecules in " << this->name() << nl << endl;
     }
 }
 
 
-void Foam::polyatomicCloud::writeXYZ(const fileName& fName) const
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::writeXYZ(const fileName& fName) const
 {
     OFstream os(fName);
 
     os  << nSites() << nl
-        << "polyatomicCloud site positions in angstroms" << nl;
+        << "MoleculeCloud<MoleculeType> site positions in angstroms" << nl;
 
-    forAllConstIter(polyatomicCloud, *this, mol)
+    forAllConstIter(typename MoleculeCloud<MoleculeType>, *this, mol)
     {
-        const polyatomic::constantProperties& cP = constProps(mol().id());
+        const typename MoleculeType::constantProperties& cP
+        (
+            constProps(mol().id())
+        );
 
         forAll(mol().sitePositions(), i)
         {

src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.H

@@ -22,22 +22,22 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::polyatomicCloud
+    Foam::MoleculeCloud
 
 Description
 
 SourceFiles
-    polyatomicCloudI.H
+    MoleculeCloudI.H
-    polyatomicCloud.C
+    MoleculeCloud.C
 
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef polyatomicCloud_H
+#ifndef MoleculeCloud_H
-#define polyatomicCloud_H
+#define MoleculeCloud_H
 
 #include "Cloud.H"
-#include "polyatomic.H"
+#include "moleculeCloud.H"
 #include "IOdictionary.H"
 #include "potential.H"
 #include "InteractionLists.H"
@@ -51,65 +51,83 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                       Class polyatomicCloud Declaration
+                       Class MoleculeCloud Declaration
 \*---------------------------------------------------------------------------*/
 
-class polyatomicCloud
+template<class MoleculeType>
+class MoleculeCloud
 :
-    public Cloud<polyatomic>
+    public Cloud<MoleculeType>,
+    public moleculeCloud
 {
 
 private:
 
     // Private data
 
+        //-
         const polyMesh& mesh_;
 
+        //-
         const potential& pot_;
 
-        List<DynamicList<polyatomic*> > cellOccupancy_;
+        //-
+        List<DynamicList<MoleculeType*> > cellOccupancy_;
 
-        InteractionLists<polyatomic> il_;
+        //-
+        InteractionLists<MoleculeType> il_;
 
-        List<polyatomic::constantProperties> constPropList_;
+        //-
+        List<typename MoleculeType::constantProperties> constPropList_;
 
+        //-
         Random rndGen_;
 
 
     // Private Member Functions
 
+        //-
         void buildConstProps();
 
+        //-
         void setSiteSizesAndPositions();
 
-        //- Determine which polyatomics are in which cells
+        //- Determine which molecules are in which cells
         void buildCellOccupancy();
 
+        //-
         void calculatePairForce();
 
+        //-
         inline void evaluatePair
         (
-            polyatomic& molI,
+            MoleculeType& molI,
-            polyatomic& molJ
+            MoleculeType& molJ
         );
 
+        //-
         inline bool evaluatePotentialLimit
         (
-            polyatomic& molI,
+            MoleculeType& molI,
-            polyatomic& molJ
+            MoleculeType& molJ
         ) const;
 
+        //-
         void calculateTetherForce();
 
+        //-
         void calculateExternalForce();
 
+        //-
         void removeHighEnergyOverlaps();
 
+        //-
         void initialiseMolecules
         (
             const IOdictionary& mdInitialiseDict
         );
 
+        //-
         void createMolecule
         (
             const point& position,
@@ -122,25 +140,28 @@ private:
             const vector& bulkVelocity
         );
 
+        //-
         label nSites() const;
 
+        //-
         inline vector equipartitionLinearVelocity
         (
             scalar temperature,
             scalar mass
         );
 
+        //-
         inline vector equipartitionAngularMomentum
         (
             scalar temperature,
-            const polyatomic::constantProperties& cP
+            const typename MoleculeType::constantProperties& cP
         );
 
         //- Disallow default bitwise copy construct
-        polyatomicCloud(const polyatomicCloud&);
+        MoleculeCloud(const MoleculeCloud&);
 
         //- Disallow default bitwise assignment
-        void operator=(const polyatomicCloud&);
+        void operator=(const MoleculeCloud&);
 
 
 public:
@@ -148,16 +169,18 @@ public:
     // Constructors
 
         //- Construct given mesh and potential references
-        polyatomicCloud
+        MoleculeCloud
         (
+            const word& cloudName,
             const polyMesh& mesh,
             const potential& pot,
             bool readFields = true
         );
 
         //- Construct given mesh, potential and mdInitialiseDict
-        polyatomicCloud
+        MoleculeCloud
         (
+            const word& cloudName,
             const polyMesh& mesh,
             const potential& pot,
             const IOdictionary& mdInitialiseDict,
@@ -167,43 +190,46 @@ public:
 
     // Member Functions
 
-        //- Evolve the polyatomics (move, calculate forces, control state etc)
+        //- Evolve the molecules (move, calculate forces, control state etc)
         void evolve();
 
+        //-
         void calculateForce();
 
-        void applyConstraintsAndThermostats
-        (
-            const scalar targetTemperature,
-            const scalar measuredTemperature
-        );
-
         //- Print cloud information
         void info() const;
 
 
         // Access
 
+            //-
             inline const polyMesh& mesh() const;
 
+            //-
             inline const potential& pot() const;
 
-            inline const List<DynamicList<polyatomic*> >& cellOccupancy() const;
+            //-
+            inline const List<DynamicList<MoleculeType*> >&
+            cellOccupancy() const;
 
-            inline const InteractionLists<polyatomic>& il() const;
+            //-
+            inline const InteractionLists<MoleculeType>& il() const;
 
-            inline const List<polyatomic::constantProperties>
+            //-
+            inline const List<typename MoleculeType::constantProperties>
             constProps() const;
 
-            inline const polyatomic::constantProperties&
+            //-
+            inline const typename MoleculeType::constantProperties&
             constProps(label id) const;
 
+            //-
             inline Random& rndGen();
 
 
     // Member Operators
 
-        //- Write polyatomic sites in XYZ format
+        //- Write molecule sites in XYZ format
         void writeXYZ(const fileName& fName) const;
 };
 
@@ -214,7 +240,13 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#include "polyatomicCloudI.H"
+#include "MoleculeCloudI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "MoleculeCloud.C"
+#endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloudI.H

@@ -29,10 +29,11 @@ using namespace Foam::constant;
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-inline void Foam::polyatomicCloud::evaluatePair
+template<class MoleculeType>
+void Foam::MoleculeCloud<MoleculeType>::evaluatePair
 (
-    polyatomic& molI,
+    MoleculeType& molI,
-    polyatomic& molJ
+    MoleculeType& molJ
 )
 {
     const pairPotentialList& pairPot = pot_.pairPotentials();
@@ -43,9 +44,15 @@ inline void Foam::polyatomicCloud::evaluatePair
 
     label idJ = molJ.id();
 
-    const polyatomic::constantProperties& constPropI(constProps(idI));
+    const typename MoleculeType::constantProperties& constPropI
+    (
+        constProps(idI)
+    );
 
-    const polyatomic::constantProperties& constPropJ(constProps(idJ));
+    const typename MoleculeType::constantProperties& constPropJ
+    (
+        constProps(idJ)
+    );
 
     forAll(constPropI.pairPotSites(), pI)
     {
@@ -149,10 +156,11 @@ inline void Foam::polyatomicCloud::evaluatePair
 }
 
 
-inline bool Foam::polyatomicCloud::evaluatePotentialLimit
+template<class MoleculeType>
+bool Foam::MoleculeCloud<MoleculeType>::evaluatePotentialLimit
 (
-    polyatomic& molI,
+    MoleculeType& molI,
-    polyatomic& molJ
+    MoleculeType& molJ
 ) const
 {
     const pairPotentialList& pairPot = pot_.pairPotentials();
@@ -163,9 +171,15 @@ inline bool Foam::polyatomicCloud::evaluatePotentialLimit
 
     label idJ = molJ.id();
 
-    const polyatomic::constantProperties& constPropI(constProps(idI));
+    const typename MoleculeType::constantProperties& constPropI
+    (
+        constProps(idI)
+    );
 
-    const polyatomic::constantProperties& constPropJ(constProps(idJ));
+    const typename MoleculeType::constantProperties& constPropJ
+    (
+        constProps(idJ)
+    );
 
     forAll(constPropI.pairPotSites(), pI)
     {
@@ -194,14 +208,18 @@ inline bool Foam::polyatomicCloud::evaluatePotentialLimit
 
                 if (rsIsJMag < SMALL)
                 {
-                    WarningIn("polyatomicCloud::removeHighEnergyOverlaps()")
+                    WarningIn
+                    (
+                        "MoleculeCloud<MoleculeType>::"
+                        "removeHighEnergyOverlaps()"
+                    )
                         << "Molecule site pair closer than "
                         << SMALL
                         << ": mag separation = " << rsIsJMag
                         << ". These may have been placed on top of each"
                         << " other by a rounding error in mdInitialise in"
-                        << " parallel or a block filled with polyatomics"
+                        << " parallel or a block filled with moleculess"
-                        << " twice. Removing one of the polyatomics."
+                        << " twice. Removing one of the molecules."
                         << endl;
 
                     return true;
@@ -250,14 +268,18 @@ inline bool Foam::polyatomicCloud::evaluatePotentialLimit
 
                 if (rsIsJMag < SMALL)
                 {
-                    WarningIn("polyatomicCloud::removeHighEnergyOverlaps()")
+                    WarningIn
+                    (
+                        "MoleculeCloud<MoleculeType>::"
+                        "removeHighEnergyOverlaps()"
+                    )
                         << "Molecule site pair closer than "
                         << SMALL
                         << ": mag separation = " << rsIsJMag
                         << ". These may have been placed on top of each"
                         << " other by a rounding error in mdInitialise in"
-                        << " parallel or a block filled with polyatomics"
+                        << " parallel or a block filled with molecules"
-                        << " twice. Removing one of the polyatomics."
+                        << " twice. Removing one of the molecules."
                         << endl;
 
                     return true;
@@ -288,7 +310,9 @@ inline bool Foam::polyatomicCloud::evaluatePotentialLimit
 }
 
 
-inline Foam::vector Foam::polyatomicCloud::equipartitionLinearVelocity
+template<class MoleculeType>
+Foam::vector
+Foam::MoleculeCloud<MoleculeType>::equipartitionLinearVelocity
 (
     scalar temperature,
     scalar mass
@@ -303,10 +327,12 @@ inline Foam::vector Foam::polyatomicCloud::equipartitionLinearVelocity
 }
 
 
-inline Foam::vector Foam::polyatomicCloud::equipartitionAngularMomentum
+template<class MoleculeType>
+Foam::vector
+Foam::MoleculeCloud<MoleculeType>::equipartitionAngularMomentum
 (
     scalar temperature,
-    const polyatomic::constantProperties& cP
+    const typename MoleculeType::constantProperties& cP
 )
 {
     scalar sqrtKbT = sqrt(physicoChemical::k.value()*temperature);
@@ -334,47 +360,54 @@ inline Foam::vector Foam::polyatomicCloud::equipartitionAngularMomentum
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline const Foam::polyMesh& Foam::polyatomicCloud::mesh() const
+template<class MoleculeType>
+const Foam::polyMesh& Foam::MoleculeCloud<MoleculeType>::mesh() const
 {
     return mesh_;
 }
 
 
-inline const Foam::potential& Foam::polyatomicCloud::pot() const
+template<class MoleculeType>
+const Foam::potential& Foam::MoleculeCloud<MoleculeType>::pot() const
 {
     return pot_;
 }
 
 
-inline const Foam::List<Foam::DynamicList<Foam::polyatomic*> >&
+template<class MoleculeType>
-    Foam::polyatomicCloud::cellOccupancy() const
+const Foam::List<Foam::DynamicList<MoleculeType*> >&
+Foam::MoleculeCloud<MoleculeType>::cellOccupancy() const
 {
     return cellOccupancy_;
 }
 
 
-inline const Foam::InteractionLists<Foam::polyatomic>&
+template<class MoleculeType>
-    Foam::polyatomicCloud::il() const
+const Foam::InteractionLists<MoleculeType>&
+Foam::MoleculeCloud<MoleculeType>::il() const
 {
     return il_;
 }
 
 
-inline const Foam::List<Foam::polyatomic::constantProperties>
+template<class MoleculeType>
-    Foam::polyatomicCloud::constProps() const
+const Foam::List<typename MoleculeType::constantProperties>
+Foam::MoleculeCloud<MoleculeType>::constProps() const
 {
     return constPropList_;
 }
 
 
-inline const Foam::polyatomic::constantProperties&
+template<class MoleculeType>
-    Foam::polyatomicCloud::constProps(label id) const
+const typename MoleculeType::constantProperties&
+Foam::MoleculeCloud<MoleculeType>::constProps(label id) const
 {
     return constPropList_[id];
 }
 
 
-inline Foam::Random& Foam::polyatomicCloud::rndGen()
+template<class MoleculeType>
+Foam::Random& Foam::MoleculeCloud<MoleculeType>::rndGen()
 {
     return rndGen_;
 }

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

@@ -0,0 +1,48 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "moleculeCloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(moleculeCloud, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::moleculeCloud::moleculeCloud()
+{}
+
+
+// * * * * * * * * * * * * * * * *  Destructors  * * * * * * * * * * * * * * //
+
+Foam::moleculeCloud::~moleculeCloud()
+{}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,83 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::moleculeCloud
+
+Description
+    Virtual abstract base class for templated moleculeCloud
+
+SourceFiles
+    moleculeCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef moleculeCloud_H
+#define moleculeCloud_H
+
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class moleculeCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+class moleculeCloud
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        moleculeCloud(const moleculeCloud&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const moleculeCloud&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("moleculeCloud");
+
+    // Constructors
+
+        //- Null constructor
+        moleculeCloud();
+
+    //- Destructor
+    virtual ~moleculeCloud();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/clouds/derived/polyatomicCloud/polyatomicCloud.H

@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::polyatomicCloud
+
+Description
+    Cloud class to simulate polyatomic molecules
+
+SourceFiles
+    polyatomicCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef polyatomicCloud_H
+#define polyatomicCloud_H
+
+#include "MoleculeCloud.H"
+#include "polyatomic.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef MoleculeCloud<polyatomic> polyatomicCloud;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.C

@@ -23,11 +23,19 @@ License
 
 \*----------------------------------------------------------------------------*/
 
-#include "polyatomicCloud.H"
 #include "polyatomic.H"
 #include "Random.H"
 #include "Time.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(polyatomic, 0);
+    defineTemplateTypeNameAndDebug(Cloud<polyatomic>, 0);
+}
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 Foam::tensor Foam::polyatomic::rotationTensorX(scalar phi) const
@@ -285,7 +293,7 @@ void Foam::polyatomic::hitWallPatch
 )
 {
     // Use of the normal from tetIs is not required as
-    // hasWallImpactDistance for a polyatomicCloud is false.
+    // hasWallImpactDistance is false.
     vector nw = normal();
     nw /= mag(nw);
 

src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.H

@@ -42,14 +42,15 @@ SourceFiles
 #include "autoPtr.H"
 #include "diagTensor.H"
 #include "constPropSite.H"
+#include "MoleculeCloud.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-// Class forward declarations
+template<class ParcelType>
-class polyatomicCloud;
+class DsmcParcel;
 
 /*---------------------------------------------------------------------------*\
                            Class polyatomic Declaration
@@ -91,22 +92,19 @@ public:
             //- Moment of intertia (in principal axis configiration)
             diagTensor momentOfInertia_;
 
+            //-
             scalar mass_;
 
 
         // Private Member Functions
 
-            void setInteracionSiteBools
+            //-
-            (
-                const List<word>& siteIds,
-                const List<word>& pairPotSiteIds
-            );
-
             bool linearMoleculeTest() const;
 
 
     public:
 
+        //-
         inline constantProperties();
 
         //- Construct from dictionary
@@ -118,22 +116,31 @@ public:
 
         // Member functions
 
+        //-
         inline const List<constPropSite>& sites() const;
 
+        //-
         inline const List<label>& pairPotSites() const;
 
+        //-
         inline const List<label>& electrostaticSites() const;
 
+        //-
         inline const diagTensor& momentOfInertia() const;
 
+        //-
         inline bool linearMolecule() const;
 
+        //-
         inline bool pointMolecule() const;
 
+        //-
         inline label degreesOfFreedom() const;
 
+        //-
         inline scalar mass() const;
 
+        //-
         inline label nSites() const;
     };
 
@@ -141,7 +148,7 @@ public:
     //- Class used to pass tracking data to the trackToFace function
     class trackingData
     :
-        public particle::TrackingData<polyatomicCloud>
+        public particle::TrackingData<MoleculeCloud<polyatomic> >
     {
         // label specifying which part of the integration algorithm is taking
         label part_;
@@ -151,9 +158,9 @@ public:
 
         // Constructors
 
-            trackingData(polyatomicCloud& cloud, label part)
+            trackingData(MoleculeCloud<polyatomic>& cloud, label part)
             :
-                particle::TrackingData<polyatomicCloud>(cloud),
+                particle::TrackingData<MoleculeCloud<polyatomic> >(cloud),
                 part_(part)
             {}
 
@@ -176,10 +183,10 @@ private:
         //  bodyLocalVector = Q_.T() & globalVector
         tensor Q_;
 
-        // Linear velocity of polyatomic
+        //- Linear velocity of polyatomic
         vector v_;
 
-        // Total linear acceleration of polyatomic
+        //- Total linear acceleration of polyatomic
         vector a_;
 
         //- Angular momentum of polyatomic, in body local reference frame
@@ -188,8 +195,10 @@ private:
         //- Total torque on polyatomic, in body local reference frame
         vector tau_;
 
+        //-
         vector specialPosition_;
 
+        //-
         scalar potentialEnergy_;
 
         // - r_ij f_ij, stress dyad
@@ -198,31 +207,42 @@ private:
 //         // - r_ij outer product f_ij: virial contribution
 //         tensor rDotf_;
 
+        //-
         label special_;
 
+        //-
         label id_;
 
+        //-
         List<vector> siteForces_;
 
+        //-
         List<vector> sitePositions_;
 
 
     // Private Member Functions
 
+
+        //-
         tensor rotationTensorX(scalar deltaT) const;
 
+        //-
         tensor rotationTensorY(scalar deltaT) const;
 
+        //-
         tensor rotationTensorZ(scalar deltaT) const;
 
 
 public:
 
+    //- Runtime type information
+    TypeName("polyatomic");
+
     friend class Cloud<polyatomic>;
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from all components
         inline polyatomic
         (
             const polyMesh& mesh,
@@ -279,53 +299,91 @@ public:
 
         // Tracking
 
+            //-
             bool move(trackingData&, const scalar trackTime);
 
+            //-
             virtual void transformProperties(const tensor& T);
 
+            //-
             virtual void transformProperties(const vector& separation);
 
+            //-
             void setSitePositions(const constantProperties& constProps);
 
+            //-
             void setSiteSizes(label size);
 
 
         // Access
 
+            //-
             inline const tensor& Q() const;
+
+            //-
             inline tensor& Q();
 
+            //-
             inline const vector& v() const;
+
+            //-
             inline vector& v();
 
+            //-
             inline const vector& a() const;
+
+            //-
             inline vector& a();
 
+            //-
             inline const vector& pi() const;
+
+            //-
             inline vector& pi();
 
+            //-
             inline const vector& tau() const;
+
+            //-
             inline vector& tau();
 
+            //-
             inline const List<vector>& siteForces() const;
+
+            //-
             inline List<vector>& siteForces();
 
+            //-
             inline const List<vector>& sitePositions() const;
+
+            //-
             inline List<vector>& sitePositions();
 
+            //-
             inline const vector& specialPosition() const;
+
+            //-
             inline vector& specialPosition();
 
+            //-
             inline scalar potentialEnergy() const;
+
+            //-
             inline scalar& potentialEnergy();
 
+            //-
             inline const tensor& rf() const;
+
+            //-
             inline tensor& rf();
 
+            //-
             inline label special() const;
 
+            //-
             inline bool tethered() const;
 
+            //-
             inline label id() const;
 
 
@@ -367,8 +425,10 @@ public:
 
     // I-O
 
+        //- Read
         static void readFields(Cloud<polyatomic>& mC);
 
+        //- Write
         static void writeFields(const Cloud<polyatomic>& mC);
 
 

src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomicIO.C

@@ -25,7 +25,7 @@ License
 
 #include "polyatomic.H"
 #include "IOstreams.H"
-#include "polyatomicCloud.H"
+#include "Cloud.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -92,7 +92,11 @@ Foam::polyatomic::polyatomic
     is.check
     (
         "Foam::polyatomic::polyatomic"
-        "(const Cloud<polyatomic>& cloud, Foam::Istream&), bool"
+        "("
+            "const polyMesh& mesh,"
+            "Istream& is,"
+            "bool readFields"
+        ")"
     );
 }
 
@@ -134,7 +138,8 @@ void Foam::polyatomic::readFields(Cloud<polyatomic>& mC)
     mC.checkFieldIOobject(mC, id);
 
     label i = 0;
-    forAllIter(polyatomicCloud, mC, iter)
+
+    forAllIter(typename Cloud<polyatomic>, mC, iter)
     {
         polyatomic& mol = iter();
 
@@ -203,7 +208,7 @@ void Foam::polyatomic::writeFields(const Cloud<polyatomic>& mC)
     );
 
     label i = 0;
-    forAllConstIter(polyatomicCloud, mC, iter)
+    forAllConstIter(typename Cloud<polyatomic>, mC, iter)
     {
         const polyatomic& mol = iter();
 
@@ -241,18 +246,6 @@ void Foam::polyatomic::writeFields(const Cloud<polyatomic>& mC)
     orientation1.write();
     orientation2.write();
     orientation3.write();
-
-    Info<< "writeFields " << mC.name() << endl;
-
-    if (isA<polyatomicCloud>(mC))
-    {
-        const polyatomicCloud& m = dynamic_cast<const polyatomicCloud&>(mC);
-
-        m.writeXYZ
-        (
-            m.mesh().time().timePath()/cloud::prefix/"polyatomicCloud.xmol"
-        );
-    }
 }
 
 

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

@@ -27,7 +27,7 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::potential::setSiteIdList(const dictionary& mdPropertiesDict)
+void Foam::potential::setSiteIdList(const dictionary& moleculePropertiesDict)
 {
     DynamicList<word> siteIdList;
     DynamicList<word> pairPotentialSiteIdList;
@@ -36,14 +36,14 @@ void Foam::potential::setSiteIdList(const dictionary& mdPropertiesDict)
     {
         const word& id(idList_[i]);
 
-        if (!mdPropertiesDict.found(id))
+        if (!moleculePropertiesDict.found(id))
         {
             FatalErrorIn("potential::setSiteIdList(const dictionary&)")
                 << id << " molecule subDict not found"
                 << nl << abort(FatalError);
         }
 
-        const dictionary& molDict(mdPropertiesDict.subDict(id));
+        const dictionary& molDict(moleculePropertiesDict.subDict(id));
 
         List<word> siteIdNames = molDict.lookup("siteIds");
 
@@ -93,7 +93,10 @@ void Foam::potential::setSiteIdList(const dictionary& mdPropertiesDict)
 }
 
 
-void Foam::potential::potential::readPotentialDict()
+void Foam::potential::potential::readPotentialDict
+(
+    const dictionary& moleculePropertiesDict
+)
 {
     Info<< nl <<  "Reading potential dictionary:" << endl;
 
@@ -111,31 +114,28 @@ void Foam::potential::potential::readPotentialDict()
 
     idList_ = List<word>(idListDict.lookup("idList"));
 
-    setSiteIdList
+    setSiteIdList(moleculePropertiesDict);
-    (
-        IOdictionary
-        (
-            IOobject
-            (
-                "mdProperties",
-                mesh_.time().constant(),
-                mesh_,
-                IOobject::MUST_READ_IF_MODIFIED,
-                IOobject::NO_WRITE,
-                false
-            )
-        )
-    );
 
     List<word> pairPotentialSiteIdList
     (
         SubList<word>(siteIdList_, nPairPotIds_)
     );
 
-    Info<< nl << "Unique site ids found: " << siteIdList_
+    Info<< nl << "Unique site ids found:";
-        << nl << "Site Ids requiring a pair potential: "
+
-        << pairPotentialSiteIdList
+    forAll(siteIdList_, i)
-        << endl;
+    {
+        Info<< " " << siteIdList_[i];
+    }
+
+    Info << nl << "Site Ids requiring a pair potential:";
+
+    forAll(pairPotentialSiteIdList, i)
+    {
+        Info<< " " << pairPotentialSiteIdList[i];
+    }
+
+    Info<< nl;
 
     List<word> tetherSiteIdList(0);
 
@@ -229,37 +229,25 @@ void Foam::potential::potential::readPotentialDict()
 
     if (potentialDict.found("external"))
     {
-        Info<< nl << "Reading external forces:" << endl;
+        Info<< nl << "Reading external forces: ";
 
         const dictionary& externalDict = potentialDict.subDict("external");
 
         // gravity
         externalDict.readIfPresent("gravity", gravity_);
-    }
 
-    Info<< nl << tab << "gravity = " << gravity_ << endl;
+        Info<< "gravity = " << gravity_ << nl << endl;
+    }
 }
 
 
 void Foam::potential::potential::readMdInitialiseDict
 (
     const IOdictionary& mdInitialiseDict,
+    const dictionary& moleculePropertiesDict,
     IOdictionary& idListDict
 )
 {
-    IOdictionary mdPropertiesDict
-    (
-        IOobject
-        (
-            "mdProperties",
-            mesh_.time().constant(),
-            mesh_,
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE,
-            false
-        )
-    );
-
     DynamicList<word> idList;
 
     DynamicList<word> tetherSiteIdList;
@@ -280,7 +268,7 @@ void Foam::potential::potential::readMdInitialiseDict
         {
             const word& id = latticeIds[i];
 
-            if (!mdPropertiesDict.found(id))
+            if (!moleculePropertiesDict.found(id))
             {
                 FatalErrorIn
                 (
@@ -322,7 +310,7 @@ void Foam::potential::potential::readMdInitialiseDict
 
                 List<word> siteIds
                 (
-                    mdPropertiesDict.subDict(id).lookup("siteIds")
+                    moleculePropertiesDict.subDict(id).lookup("siteIds")
                 );
 
                 if (findIndex(siteIds, tetherSiteId) != -1)
@@ -359,16 +347,20 @@ void Foam::potential::potential::readMdInitialiseDict
 
     idListDict.add("tetherSiteIdList", tetherSiteIdList);
 
-    setSiteIdList(mdPropertiesDict);
+    setSiteIdList(moleculePropertiesDict);
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::potential::potential(const polyMesh& mesh)
+Foam::potential::potential
+(
+    const polyMesh& mesh,
+    const dictionary& moleculePropertiesDict
+)
 :
     mesh_(mesh)
 {
-    readPotentialDict();
+    readPotentialDict(moleculePropertiesDict);
 }
 
 
@@ -376,12 +368,18 @@ Foam::potential::potential
 (
     const polyMesh& mesh,
     const IOdictionary& mdInitialiseDict,
+    const dictionary& moleculePropertiesDict,
     IOdictionary& idListDict
 )
 :
     mesh_(mesh)
 {
-    readMdInitialiseDict(mdInitialiseDict, idListDict);
+    readMdInitialiseDict
+    (
+        mdInitialiseDict,
+        moleculePropertiesDict,
+        idListDict
+    );
 }
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

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

@@ -55,34 +55,48 @@ class potential
 {
     // Private data
 
+        //-
         const polyMesh& mesh_;
 
+        //-
         List<word> idList_;
 
+        //-
         List<word> siteIdList_;
 
+        //-
         label nPairPotIds_;
 
+        //-
         scalar potentialEnergyLimit_;
 
+        //-
         labelList removalOrder_;
 
+        //-
         pairPotentialList pairPotentials_;
 
+        //-
         tetherPotentialList tetherPotentials_;
 
+        //-
         vector gravity_;
 
 
     // Private Member Functions
 
-        void setSiteIdList(const dictionary& mdPropertiesDict);
 
-        void readPotentialDict();
+        //-
+        void setSiteIdList(const dictionary& moleculePropertiesDict);
 
+        //-
+        void readPotentialDict(const dictionary& moleculePropertiesDict);
+
+        //-
         void readMdInitialiseDict
         (
             const IOdictionary& mdInitialiseDict,
+            const dictionary& moleculePropertiesDict,
             IOdictionary& idListDict
         );
 
@@ -98,13 +112,18 @@ public:
     // Constructors
 
         //- Construct from mesh reference
-        potential(const polyMesh& mesh);
+        potential
+        (
+            const polyMesh& mesh,
+            const dictionary& moleculePropertiesDict
+        );
 
         //- Construct from mdInitialiseDict
         potential
         (
             const polyMesh& mesh,
             const IOdictionary& mdInitialiseDict,
+            const dictionary& moleculePropertiesDict,
             IOdictionary& idListDict
         );
 
@@ -117,22 +136,31 @@ public:
 
         // Access
 
+            //-
             inline label nIds() const;
 
+            //-
             inline const List<word>& idList() const;
 
+            //-
             inline const List<word>& siteIdList() const;
 
+            //-
             inline scalar potentialEnergyLimit() const;
 
+            //-
             inline label nPairPotentials() const;
 
+            //-
             inline const labelList& removalOrder() const;
 
+            //-
             inline const pairPotentialList& pairPotentials() const;
 
+            //-
             inline const tetherPotentialList& tetherPotentials() const;
 
+            //-
             inline const vector& gravity() const;
 };