applying update from Graham

applications/solvers/molecularDynamics/gnemdFoam/Make/files

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

applications/solvers/molecularDynamics/mdEquilibrationFoam/Make/options

@@ -1,6 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \
     -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
@@ -10,4 +11,6 @@ EXE_LIBS = \
     -lfiniteVolume \
     -llagrangian \
     -lmolecule \
-    -lpotential
+    -lpotential \
+    -lmolecularMeasurements
+

applications/solvers/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C

@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    mdEquilibrationFOAM
+    mdEquilibrationFoam
 
 Description
     Equilibrates and/or preconditions MD systems
@@ -40,9 +40,9 @@ int main(int argc, char *argv[])
 #   include "createTime.H"
 #   include "createMesh.H"
 
-    moleculeCloud molecules(mesh);
+    potential pot(mesh);
 
-    molecules.removeHighEnergyOverlaps();
+    moleculeCloud molecules(mesh, pot);
 
 #   include "temperatureAndPressureVariables.H"
 
@@ -60,7 +60,7 @@ int main(int argc, char *argv[])
 
         Info << "Time = " << runTime.timeName() << endl;
 
-        molecules.integrateEquationsOfMotion();
+        molecules.evolve();
 
 #       include "meanMomentumEnergyAndNMols.H"
 

applications/solvers/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H

@@ -1,4 +1,4 @@
-Info<< "Reading MD Equilibration Dictionary" << nl << endl;
+Info<< nl << "Reading MD Equilibration Dictionary" << nl << endl;
 
 IOdictionary mdEquilibrationDict
 (

applications/solvers/molecularDynamics/mdFoam/Make/files

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

applications/solvers/molecularDynamics/mdFoam/Make/options

@@ -1,6 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \
     -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
@@ -10,4 +11,6 @@ EXE_LIBS = \
     -lfiniteVolume \
     -llagrangian \
     -lmolecule \
-    -lpotential
+    -lpotential \
+    -lmolecularMeasurements
+

applications/solvers/molecularDynamics/mdFoam/mdFoam.C

@@ -23,10 +23,10 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    gnemdFOAM
+    mdFoam
 
 Description
-    MD for Fluid Mechanics and hybridising with a continuum solver.
+    molecular dynamics solver for fluid dynamics
 
 \*---------------------------------------------------------------------------*/
 
@@ -40,11 +40,9 @@ int main(int argc, char *argv[])
 #   include "createTime.H"
 #   include "createMesh.H"
 
-    moleculeCloud molecules(mesh);
+    potential pot(mesh);
 
-#   include "createMDFields.H"
+    moleculeCloud molecules(mesh, pot);
-
-    molecules.removeHighEnergyOverlaps();
 
 #   include "temperatureAndPressureVariables.H"
 
@@ -60,20 +58,14 @@ int main(int argc, char *argv[])
 
         Info << "Time = " << runTime.timeName() << endl;
 
-        molecules.integrateEquationsOfMotion();
+        molecules.evolve();
 
 #       include "meanMomentumEnergyAndNMols.H"
 
 #       include "temperatureAndPressure.H"
 
-#       include "calculateMDFields.H"
-
-#       include "averageMDFields.H"
-
         runTime.write();
 
-#       include "resetMDFields.H"
-
         if (runTime.outputTime())
         {
             nAveragingSteps = 0;

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(15);
+
+    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;
+}
+
+
+// ************************************************************************* //

applications/utilities/preProcessing/mdInitialise/molConfig_old/molConfig.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/lagrangian/molecularDynamics/Allwmake

@@ -3,6 +3,7 @@ cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wmake libso potential
+wmake libso molecularMeasurements
 wmake libso molecule
 
 # ----------------------------------------------------------------- end-of-file

src/lagrangian/molecularDynamics/molecularMeasurements/Make/files

@@ -0,0 +1,3 @@
+distribution/distribution.C
+
+LIB = $(FOAM_LIBBIN)/libmolecularMeasurements

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

@@ -25,43 +25,43 @@ License
 \*---------------------------------------------------------------------------*/
 
 template<class Type>
-inline const Field< Field<Type> >& Foam::correlationFunction<Type>::
+inline const Foam::Field< Foam::Field<Type> >&
-tZeroBuffers() const
+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/Make/files

@@ -1,47 +1,27 @@
-correlationFunction = correlationFunction
+interactionLists = interactionLists
+referredMolecule = $(interactionLists)/referredMolecule
+referredCellList = $(interactionLists)/referredCellList
+referredCell = $(interactionLists)/referredCell
+referralLists = $(interactionLists)/referralLists
+directInteractionList = $(interactionLists)/directInteractionList
 
-distribution = distribution
-
-molecule = molecule
-moleculeCloud = moleculeCloud
-
-reducedUnits = reducedUnits
-
-referredMolecule = referredMolecule
-referredCellList = referredCellList
-referredCell = referredCell
-referralLists = referralLists
-$(distribution)/distribution.C
-
-$(reducedUnits)/reducedUnits.C
-$(reducedUnits)/reducedUnitsIO.C
-
-$(molecule)/molecule.C
-$(molecule)/moleculeIO.C
-
-$(moleculeCloud)/moleculeCloud.C
-$(moleculeCloud)/moleculeCloudBuildCellOccupancy.C
-$(moleculeCloud)/moleculeCloudBuildCellInteractionLists.C
-$(moleculeCloud)/moleculeCloudBuildCellReferralLists.C
-
-$(moleculeCloud)/moleculeCloudTestEdgeEdgeDistance.C
-$(moleculeCloud)/moleculeCloudTestPointFaceDistance.C
-$(moleculeCloud)/moleculeCloudRealCellsInRangeOfSegment.C
-$(moleculeCloud)/moleculeCloudReferredCellsInRangeOfSegment.C
-
-$(moleculeCloud)/moleculeCloudCalculateForce.C
-$(moleculeCloud)/moleculeCloudCalculatePairForce.C
-$(moleculeCloud)/moleculeCloudCalculateTetherForce.C
-$(moleculeCloud)/moleculeCloudCalculateExternalForce.C
-$(moleculeCloud)/moleculeCloudIntegrateEquationsOfMotion.C
-$(moleculeCloud)/moleculeCloudRemoveHighEnergyOverlaps.C
-$(moleculeCloud)/moleculeCloudApplyConstraintsAndThermostats.C
-
-$(referralLists)/receivingReferralList.C
 $(referralLists)/sendingReferralList.C
+$(referralLists)/receivingReferralList.C
 $(referredCellList)/referredCellList.C
 $(referredCell)/referredCell.C
 $(referredMolecule)/referredMolecule.C
+$(directInteractionList)/directInteractionList.C
+$(interactionLists)/interactionLists.C
+
+reducedUnits = reducedUnits
+$(reducedUnits)/reducedUnits.C
+$(reducedUnits)/reducedUnitsIO.C
+
+molecule = molecule
+$(molecule)/molecule.C
+$(molecule)/moleculeIO.C
+
+moleculeCloud = moleculeCloud
+$(moleculeCloud)/moleculeCloud.C
 
 LIB = $(FOAM_LIBBIN)/libmolecule
-

src/lagrangian/molecularDynamics/molecule/Make/options

@@ -1,10 +1,12 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/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
 
 EXE_LIBS = \
     -lfiniteVolume \
     -llagrangian \
-    -lpotential
+    -lpotential \
+	-lmolecularMeasurements
 

src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionList.C

@@ -0,0 +1,320 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "directInteractionList.H"
+#include "interactionLists.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::directInteractionList::buildDirectInteractionList
+(
+    bool pointPointListBuild
+)
+{
+    Info<< nl << "Building list of direct interaction neighbours" << endl;
+
+    const polyMesh& mesh(il_.mesh());
+
+    List<DynamicList<label> > directInteractionList(mesh.nCells());
+
+    if (pointPointListBuild)
+    {
+        Info<< tab << "Point-Point direct interaction list build." << endl;
+
+        label pointJIndex;
+
+        forAll (mesh.points(), pointIIndex)
+        {
+            for
+            (
+                pointJIndex = pointIIndex;
+                pointJIndex != mesh.points().size();
+                ++pointJIndex
+            )
+            {
+                if (il_.testPointPointDistance(pointIIndex, pointJIndex))
+                {
+                    const labelList& ptICells
+                    (
+                        mesh.pointCells()[pointIIndex]
+                    );
+
+                    const labelList& ptJCells
+                    (
+                        mesh.pointCells()[pointJIndex]
+                    );
+
+                    forAll(ptICells, pIC)
+                    {
+                        const label cellI(ptICells[pIC]);
+
+                        forAll(ptJCells, pJC)
+                        {
+                            const label cellJ(ptJCells[pJC]);
+
+                            if (cellJ > cellI)
+                            {
+                                if
+                                (
+                                    findIndex(directInteractionList[cellI],
+                                    cellJ) == -1
+                                )
+                                {
+                                    directInteractionList[cellI].append(cellJ);
+                                }
+                            }
+
+                            if (cellI > cellJ)
+                            {
+                                if
+                                (
+                                    findIndex(directInteractionList[cellJ],
+                                    cellI) == -1
+                                )
+                                {
+                                    directInteractionList[cellJ].append(cellI);
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+    else
+    {
+        Info<< tab << "Point-Face, Edge-Edge direct interaction list build."
+            << endl;
+
+        forAll (mesh.points(), p)
+        {
+            forAll(mesh.faces(), f)
+            {
+                if(il_.testPointFaceDistance(p, f))
+                {
+                    const labelList& pCells(mesh.pointCells()[p]);
+
+                    const label cellO(mesh.faceOwner()[f]);
+
+                    const label cellN(mesh.faceNeighbour()[f]);
+
+                    forAll(pCells, pC)
+                    {
+                        const label cellI(pCells[pC]);
+
+                        // cells are not added to their own DIL
+
+                        if (cellO > cellI)
+                        {
+                            if
+                            (
+                                findIndex(directInteractionList[cellI],
+                                cellO) == -1
+                            )
+                            {
+                                directInteractionList[cellI].append(cellO);
+                            }
+                        }
+
+                        if (cellI > cellO)
+                        {
+                            if
+                            (
+                                findIndex(directInteractionList[cellO],
+                                cellI) == -1
+                            )
+                            {
+                                directInteractionList[cellO].append(cellI);
+                            }
+                        }
+
+                        if (mesh.isInternalFace(f))
+                        {
+                            // boundary faces will not have neighbour
+                            // information
+
+                            if (cellN > cellI)
+                            {
+                                if
+                                (
+                                    findIndex(directInteractionList[cellI],
+                                    cellN) == -1
+                                )
+                                {
+                                    directInteractionList[cellI].append(cellN);
+                                }
+                            }
+
+                            if (cellI > cellN)
+                            {
+                                if
+                                (
+                                    findIndex(directInteractionList[cellN],
+                                    cellI) == -1
+                                )
+                                {
+                                    directInteractionList[cellN].append(cellI);
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+
+        label edgeJIndex;
+
+        forAll (mesh.edges(), edgeIIndex)
+        {
+            const edge& eI(mesh.edges()[edgeIIndex]);
+
+            for
+            (
+                edgeJIndex = edgeIIndex + 1;
+                edgeJIndex != mesh.edges().size();
+                ++edgeJIndex
+            )
+            {
+                const edge& eJ(mesh.edges()[edgeJIndex]);
+
+                if (il_.testEdgeEdgeDistance(eI, eJ))
+                {
+                    const labelList& eICells(mesh.edgeCells()[edgeIIndex]);
+
+                    const labelList& eJCells(mesh.edgeCells()[edgeJIndex]);
+
+                    forAll(eICells, eIC)
+                    {
+                        const label cellI(eICells[eIC]);
+
+                        forAll(eJCells, eJC)
+                        {
+                            const label cellJ(eJCells[eJC]);
+
+                            if (cellJ > cellI)
+                            {
+                                if
+                                (
+                                    findIndex(directInteractionList[cellI],
+                                    cellJ) == -1
+                                )
+                                {
+                                    directInteractionList[cellI].append(cellJ);
+                                }
+                            }
+
+                            if (cellI > cellJ)
+                            {
+                                if
+                                (
+                                    findIndex(directInteractionList[cellJ],
+                                    cellI) == -1
+                                )
+                                {
+                                    directInteractionList[cellJ].append(cellI);
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    forAll(directInteractionList, transDIL)
+    {
+        (*this)[transDIL].transfer
+        (
+            directInteractionList[transDIL].shrink()
+        );
+    }
+
+    // sorting DILs
+
+    forAll((*this), dIL)
+    {
+        sort((*this)[dIL]);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::directInteractionList::directInteractionList
+(
+    const interactionLists& il,
+    bool pointPointListBuild
+)
+:
+    labelListList(il.mesh().nCells()),
+    il_(il)
+{
+    if((*this).size() > 1)
+    {
+        buildDirectInteractionList(pointPointListBuild);
+    }
+    else if((*this).size() == 1)
+    {
+        Info<< nl
+            << "Single cell mesh, no direct interaction lists required."
+            << endl;
+
+        (*this)[0].setSize(0);
+    }
+}
+
+
+Foam::directInteractionList::directInteractionList
+(
+    const interactionLists& il
+)
+:
+    labelListList(il.mesh().nCells()),
+    il_(il)
+{
+    Info<< "Read directInteractionList from disk not implemented" << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::directInteractionList::~directInteractionList()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionList.H

@@ -0,0 +1,129 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Class
+    Foam::directInteractionList
+
+Description
+
+SourceFiles
+    directInteractionListI.H
+    directInteractionList.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef directInteractionList_H
+#define directInteractionList_H
+
+#include "polyMesh.H"
+#include "List.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class interactionLists;
+
+/*---------------------------------------------------------------------------*\
+                           Class directInteractionList Declaration
+\*---------------------------------------------------------------------------*/
+
+class directInteractionList
+:
+    public labelListList
+{
+    // Private data
+
+        const interactionLists& il_;
+
+    // Private Member Functions
+
+        void buildDirectInteractionList
+        (
+            bool pointPointListBuild
+        );
+
+        //- Disallow default bitwise copy construct
+        directInteractionList(const directInteractionList&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const directInteractionList&);
+
+public:
+
+    // Constructors
+
+        //- Construct lists by searching the mesh
+        directInteractionList
+        (
+            const interactionLists& il,
+            bool pointPointListBuild
+        );
+
+        //- Construct from file
+        directInteractionList
+        (
+            const interactionLists& il
+        );
+
+    // Destructor
+
+        ~directInteractionList();
+
+
+    // Member Functions
+
+        // Access
+
+        inline const interactionLists& il() const;
+
+        // Check
+
+        // Edit
+
+        // Write
+
+
+    // IOstream Operators
+
+        friend Istream& operator>>(Istream&, directInteractionList&);
+        friend Ostream& operator<<(Ostream&, const directInteractionList&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "directInteractionListI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionListI.H

@@ -22,20 +22,19 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-\*----------------------------------------------------------------------------*/
+\*---------------------------------------------------------------------------*/
-
-#include "moleculeCloud.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-void Foam::moleculeCloud::calculatePairForce()
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::interactionLists& Foam::directInteractionList::il() const
 {
-    iterator mol(this->begin());
+    return il_;
-
-#   include "moleculeCloudCalculatePairForceRealCells.H"
-
-#   include "moleculeCloudCalculatePairForceReferredCells.H"
 }
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.C

@@ -0,0 +1,682 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "interactionLists.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+Foam::scalar Foam::interactionLists::transTol = 1e-12;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::interactionLists::buildCellReferralLists()
+{
+    Info<< nl << "Determining molecule referring schedule" << endl;
+
+    const referredCellList& refIntL(ril());
+
+    DynamicList<label> referralProcs;
+
+    // Run through all referredCells to build list of interacting processors
+
+    forAll(refIntL, rIL)
+    {
+        const referredCell& rC(refIntL[rIL]);
+
+        if (findIndex(referralProcs, rC.sourceProc()) == -1)
+        {
+            referralProcs.append(rC.sourceProc());
+        }
+    }
+
+    referralProcs.shrink();
+
+//     Pout << "referralProcs: " << nl << referralProcs << endl;
+
+    List<DynamicList<label> > cellSendingReferralLists(referralProcs.size());
+
+    List<DynamicList<DynamicList<label> > >
+        cellReceivingReferralLists(referralProcs.size());
+
+    // Run through all referredCells again building up send and receive info
+
+    forAll(refIntL, rIL)
+    {
+        const referredCell& rC(refIntL[rIL]);
+
+        label rPI = findIndex(referralProcs, rC.sourceProc());
+
+        DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
+
+        DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
+
+        label existingSource = findIndex(sRL, rC.sourceCell());
+
+        // Check to see if this source cell has already been allocated to
+        // come to this processor.  If not, add the source cell to the sending
+        // list and add the current referred cell to the receiving list.
+
+        // It shouldn't be possible for the sending and receiving lists to be
+        // different lengths, because their append operations happen at the
+        // same time.
+
+        if (existingSource == -1)
+        {
+            sRL.append(rC.sourceCell());
+
+            rRL.append
+            (
+                DynamicList<label> (labelList(1,rIL))
+            );
+        }
+        else
+        {
+            rRL[existingSource].append(rIL);
+
+            rRL[existingSource].shrink();
+        }
+    }
+
+    forAll(referralProcs, rPI)
+    {
+        DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
+
+        DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
+
+        sRL.shrink();
+
+        rRL.shrink();
+    }
+
+    // It is assumed that cell exchange is reciprocal, if proc A has cells to
+    // send to proc B, then proc B must have some to send to proc A.
+
+    cellReceivingReferralLists_.setSize(referralProcs.size());
+
+    cellSendingReferralLists_.setSize(referralProcs.size());
+
+    forAll(referralProcs, rPI)
+    {
+        DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
+
+        labelListList translLL(rRL.size());
+
+        forAll(rRL, rRLI)
+        {
+            translLL[rRLI] = rRL[rRLI];
+        }
+
+        cellReceivingReferralLists_[rPI] = receivingReferralList
+        (
+            referralProcs[rPI],
+            translLL
+        );
+    }
+
+    // Send sendingReferralLists to each interacting processor.
+
+    forAll(referralProcs, rPI)
+    {
+
+        DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
+
+        if (referralProcs[rPI] != Pstream::myProcNo())
+        {
+            if (Pstream::parRun())
+            {
+                OPstream toInteractingProc
+                (
+                    Pstream::blocking,
+                    referralProcs[rPI]
+                );
+
+                toInteractingProc << sendingReferralList
+                (
+                    Pstream::myProcNo(),
+                    sRL
+                );
+            }
+        }
+    }
+
+    // Receive sendingReferralLists from each interacting processor.
+
+    forAll(referralProcs, rPI)
+    {
+        if (referralProcs[rPI] != Pstream::myProcNo())
+        {
+            if (Pstream::parRun())
+            {
+                IPstream fromInteractingProc
+                (
+                    Pstream::blocking,
+                    referralProcs[rPI]
+                );
+
+                fromInteractingProc >> cellSendingReferralLists_[rPI];
+            }
+        }
+        else
+        {
+            cellSendingReferralLists_[rPI] = sendingReferralList
+            (
+                Pstream::myProcNo(),
+                cellSendingReferralLists[rPI]
+            );
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+
+Foam::interactionLists::interactionLists
+(
+    const polyMesh& mesh,
+    scalar rCutMaxSqr,
+    bool pointPointListBuild
+)
+:
+    mesh_(mesh),
+    rCutMaxSqr_(rCutMaxSqr),
+    dil_(*this, pointPointListBuild),
+    ril_(*this, pointPointListBuild),
+    cellSendingReferralLists_(),
+    cellReceivingReferralLists_()
+{
+    buildCellReferralLists();
+}
+
+
+Foam::interactionLists::interactionLists(const polyMesh& mesh)
+:
+    mesh_(mesh),
+    dil_(*this),
+    ril_(*this)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::interactionLists::~interactionLists()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::interactionLists::testPointPointDistance
+(
+    const label ptI,
+    const label ptJ
+) const
+{
+    return (magSqr(mesh_.points()[ptI] - mesh_.points()[ptJ]) <= rCutMaxSqr_);
+}
+
+
+bool Foam::interactionLists::testEdgeEdgeDistance
+(
+    const edge& eI,
+    const edge& eJ
+) const
+{
+    const vector& eJs(mesh_.points()[eJ.start()]);
+    const vector& eJe(mesh_.points()[eJ.end()]);
+
+    return testEdgeEdgeDistance(eI, eJs, eJe);
+}
+
+
+bool Foam::interactionLists::testPointFaceDistance
+(
+    const label p,
+    const label faceNo
+) const
+{
+    const vector& pointPosition(mesh_.points()[p]);
+
+    return testPointFaceDistance(pointPosition, faceNo);
+}
+
+
+bool Foam::interactionLists::testPointFaceDistance
+(
+    const label p,
+    const referredCell& refCell
+) const
+{
+    const vector& pointPosition(mesh_.points()[p]);
+
+    forAll (refCell.faces(), rCF)
+    {
+        if
+        (
+            testPointFaceDistance
+            (
+                pointPosition,
+                refCell.faces()[rCF],
+                refCell.vertexPositions(),
+                refCell.faceCentres()[rCF],
+                refCell.faceAreas()[rCF]
+            )
+        )
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
+bool Foam::interactionLists::testPointFaceDistance
+(
+    const vectorList& pointsToTest,
+    const label faceNo
+) const
+{
+    forAll(pointsToTest, pTT)
+    {
+        const vector& p(pointsToTest[pTT]);
+
+        // if any point in the list is in range of the face
+        // then the rest do not need to be tested and
+        // true can be returned
+
+        if (testPointFaceDistance(p, faceNo))
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
+bool Foam::interactionLists::testPointFaceDistance
+(
+    const vector& p,
+    const label faceNo
+) const
+{
+    const face& faceToTest(mesh_.faces()[faceNo]);
+
+    const vector& faceC(mesh_.faceCentres()[faceNo]);
+
+    const vector& faceA(mesh_.faceAreas()[faceNo]);
+
+    const vectorList& points(mesh_.points());
+
+    return testPointFaceDistance
+    (
+        p,
+        faceToTest,
+        points,
+        faceC,
+        faceA
+    );
+}
+
+bool Foam::interactionLists::testPointFaceDistance
+(
+    const vector& p,
+    const labelList& faceToTest,
+    const vectorList& points,
+    const vector& faceC,
+    const vector& faceA
+) const
+{
+    vector faceN(faceA/mag(faceA));
+
+    scalar perpDist((p - faceC) & faceN);
+
+    if (magSqr(perpDist) > rCutMaxSqr_)
+    {
+        return false;
+    }
+
+    vector pointOnPlane = (p - faceN * perpDist);
+
+    if (magSqr(faceC - pointOnPlane) < rCutMaxSqr_*1e-8)
+    {
+        // If pointOnPlane is very close to the face centre
+        // then defining the local axes will be inaccurate
+        // and it is very likely that pointOnPlane will be
+        // inside the face, so return true if the points
+        // are in range to be safe
+
+        return (magSqr(pointOnPlane - p) <= rCutMaxSqr_);
+    }
+
+    vector xAxis = (faceC - pointOnPlane)/mag(faceC - pointOnPlane);
+
+    vector yAxis =
+        ((faceC - pointOnPlane) ^ faceN)
+       /mag((faceC - pointOnPlane) ^ faceN);
+
+    List<vector2D> local2DVertices(faceToTest.size());
+
+    forAll(faceToTest, fTT)
+    {
+        const vector& V(points[faceToTest[fTT]]);
+
+        if (magSqr(V-p) <= rCutMaxSqr_)
+        {
+            return true;
+        }
+
+        local2DVertices[fTT] = vector2D
+        (
+            ((V - pointOnPlane) & xAxis),
+            ((V - pointOnPlane) & yAxis)
+        );
+    }
+
+    scalar localFaceCx((faceC - pointOnPlane) & xAxis);
+
+    scalar la_valid = -1;
+
+    forAll(local2DVertices, fV)
+    {
+        const vector2D& va(local2DVertices[fV]);
+
+        const vector2D& vb
+        (
+            local2DVertices[(fV + 1) % local2DVertices.size()]
+        );
+
+        if (mag(vb.y()-va.y()) > SMALL)
+        {
+            scalar la =
+                (
+                    va.x() - va.y()*((vb.x() - va.x())/(vb.y() - va.y()))
+                )
+               /localFaceCx;
+
+            scalar lv = -va.y()/(vb.y() - va.y());
+
+
+            if (la >= 0 && la <= 1 && lv >= 0 && lv <= 1)
+            {
+                la_valid = la;
+
+                break;
+            }
+        }
+    }
+
+    if (la_valid < 0)
+    {
+        // perpendicular point inside face, nearest point is pointOnPlane;
+        return (magSqr(pointOnPlane-p) <= rCutMaxSqr_);
+    }
+    else
+    {
+        // perpendicular point outside face, nearest point is
+        // on edge that generated la_valid;
+        return
+        (
+            magSqr(pointOnPlane + la_valid*(faceC - pointOnPlane) - p)
+         <= rCutMaxSqr_
+        );
+    }
+
+    // if the algorithm hasn't returned anything by now then something has
+    // gone wrong.
+
+    FatalErrorIn("interactionLists.C") << nl
+        << "point " << p << " to face " << faceToTest
+        << " comparison did not find a nearest point"
+        << " to be inside or outside face."
+        << abort(FatalError);
+
+    return false;
+}
+
+
+bool Foam::interactionLists::testEdgeEdgeDistance
+(
+    const edge& eI,
+    const vector& eJs,
+    const vector& eJe
+) const
+{
+    vector a(eI.vec(mesh_.points()));
+    vector b(eJe - eJs);
+
+    const vector& eIs(mesh_.points()[eI.start()]);
+
+    vector c(eJs - eIs);
+
+    vector crossab = a ^ b;
+    scalar magCrossSqr = magSqr(crossab);
+
+    if (magCrossSqr > VSMALL)
+    {
+        // If the edges are parallel then a point-face
+        // search will pick them up
+
+        scalar s = ((c ^ b) & crossab)/magCrossSqr;
+        scalar t = ((c ^ a) & crossab)/magCrossSqr;
+
+        // Check for end points outside of range 0..1
+        // If the closest point is outside this range
+        // a point-face search will have found it.
+
+        return
+        (
+            s >= 0
+         && s <= 1
+         && t >= 0
+         && t <= 1
+         && magSqr(eIs + a*s - eJs - b*t) <= rCutMaxSqr_
+        );
+    }
+
+    return false;
+}
+
+
+const Foam::labelList Foam::interactionLists::realCellsInRangeOfSegment
+(
+    const labelList& segmentFaces,
+    const labelList& segmentEdges,
+    const labelList& segmentPoints
+) const
+{
+    DynamicList<label> realCellsFoundInRange;
+
+    forAll(segmentFaces, sF)
+    {
+        const label f = segmentFaces[sF];
+
+        forAll (mesh_.points(), p)
+        {
+            if (testPointFaceDistance(p, f))
+            {
+                const labelList& pCells(mesh_.pointCells()[p]);
+
+                forAll(pCells, pC)
+                {
+                    const label cellI(pCells[pC]);
+
+                    if (findIndex(realCellsFoundInRange, cellI) == -1)
+                    {
+                        realCellsFoundInRange.append(cellI);
+                    }
+                }
+            }
+        }
+    }
+
+    forAll(segmentPoints, sP)
+    {
+        const label p = segmentPoints[sP];
+
+        forAll(mesh_.faces(), f)
+        {
+            if (testPointFaceDistance(p, f))
+            {
+                const label cellO(mesh_.faceOwner()[f]);
+
+                if (findIndex(realCellsFoundInRange, cellO) == -1)
+                {
+                    realCellsFoundInRange.append(cellO);
+                }
+
+                if (mesh_.isInternalFace(f))
+                {
+                    // boundary faces will not have neighbour information
+
+                    const label cellN(mesh_.faceNeighbour()[f]);
+
+                    if (findIndex(realCellsFoundInRange, cellN) == -1)
+                    {
+                        realCellsFoundInRange.append(cellN);
+                    }
+                }
+            }
+        }
+    }
+
+    forAll(segmentEdges, sE)
+    {
+        const edge& eJ(mesh_.edges()[segmentEdges[sE]]);
+
+        forAll (mesh_.edges(), edgeIIndex)
+        {
+            const edge& eI(mesh_.edges()[edgeIIndex]);
+
+            if (testEdgeEdgeDistance(eI, eJ))
+            {
+                const labelList& eICells(mesh_.edgeCells()[edgeIIndex]);
+
+                forAll(eICells, eIC)
+                {
+                    const label cellI(eICells[eIC]);
+
+                    if (findIndex(realCellsFoundInRange, cellI) == -1)
+                    {
+                        realCellsFoundInRange.append(cellI);
+                    }
+                }
+            }
+        }
+    }
+
+    return realCellsFoundInRange.shrink();
+}
+
+
+const Foam::labelList Foam::interactionLists::referredCellsInRangeOfSegment
+(
+    const List<referredCell>& referredInteractionList,
+    const labelList& segmentFaces,
+    const labelList& segmentEdges,
+    const labelList& segmentPoints
+) const
+{
+    DynamicList<label> referredCellsFoundInRange;
+
+    forAll(segmentFaces, sF)
+    {
+        const label f = segmentFaces[sF];
+
+        forAll(referredInteractionList, rIL)
+        {
+            const vectorList& refCellPoints
+                = referredInteractionList[rIL].vertexPositions();
+
+            if (testPointFaceDistance(refCellPoints, f))
+            {
+                if (findIndex(referredCellsFoundInRange, rIL) == -1)
+                {
+                    referredCellsFoundInRange.append(rIL);
+                }
+            }
+        }
+    }
+
+    forAll(segmentPoints, sP)
+    {
+        const label p = segmentPoints[sP];
+
+        forAll(referredInteractionList, rIL)
+        {
+            const referredCell& refCell(referredInteractionList[rIL]);
+
+            if (testPointFaceDistance(p, refCell))
+            {
+                if (findIndex(referredCellsFoundInRange, rIL) == -1)
+                {
+                    referredCellsFoundInRange.append(rIL);
+                }
+            }
+        }
+    }
+
+    forAll(segmentEdges, sE)
+    {
+        const edge& eI(mesh_.edges()[segmentEdges[sE]]);
+
+        forAll(referredInteractionList, rIL)
+        {
+            const vectorList& refCellPoints
+                = referredInteractionList[rIL].vertexPositions();
+
+            const edgeList& refCellEdges
+                = referredInteractionList[rIL].edges();
+
+            forAll(refCellEdges, rCE)
+            {
+                const edge& eJ(refCellEdges[rCE]);
+
+                if
+                (
+                    testEdgeEdgeDistance
+                    (
+                        eI,
+                        refCellPoints[eJ.start()],
+                        refCellPoints[eJ.end()]
+                    )
+                )
+                {
+                    if(findIndex(referredCellsFoundInRange, rIL) == -1)
+                    {
+                        referredCellsFoundInRange.append(rIL);
+                    }
+                }
+            }
+        }
+    }
+
+    return referredCellsFoundInRange.shrink();
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.H

@@ -0,0 +1,211 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Class
+    Foam::interactionLists
+
+Description
+
+SourceFiles
+    interactionListsI.H
+    interactionLists.C
+    interactionListsIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef interactionLists_H
+#define interactionLists_H
+
+#include "polyMesh.H"
+#include "vector2D.H"
+#include "directInteractionList.H"
+#include "referredCell.H"
+#include "referredCellList.H"
+#include "sendingReferralList.H"
+#include "receivingReferralList.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class interactionLists Declaration
+\*---------------------------------------------------------------------------*/
+
+class interactionLists
+{
+    // Private data
+
+        const polyMesh& mesh_;
+
+        scalar rCutMaxSqr_;
+
+        directInteractionList dil_;
+
+        referredCellList ril_;
+
+        List<sendingReferralList> cellSendingReferralLists_;
+
+        List<receivingReferralList> cellReceivingReferralLists_;
+
+    // Private Member Functions
+
+        //- Build referralLists which define how to send information
+        // to referredCells to source cells
+        void buildCellReferralLists();
+
+        //- Disallow default bitwise copy construct
+        interactionLists(const interactionLists&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const interactionLists&);
+
+public:
+
+    // Static data members
+
+        //- Tolerance for checking that faces on a patch segment
+        static scalar transTol;
+
+    // Constructors
+
+        //- Construct and create all information from the mesh
+        interactionLists
+        (
+            const polyMesh& mesh,
+            scalar rCutMaxSqr,
+            bool pointPointListBuild = false
+        );
+
+        //- Construct from file
+        interactionLists(const polyMesh& mesh);
+
+    // Destructor
+
+        ~interactionLists();
+
+
+    // Member Functions
+
+            bool testPointPointDistance
+            (
+                const label ptI,
+                const label ptJ
+            ) const;
+
+            bool testPointFaceDistance
+            (
+                const label p,
+                const label faceNo
+            ) const;
+
+            bool testPointFaceDistance
+            (
+                const label p,
+                const referredCell& refCell
+            ) const;
+
+            bool testPointFaceDistance
+            (
+                const vectorList& pointsToTest,
+                const label faceNo
+            ) const;
+
+            bool testPointFaceDistance
+            (
+                const vector& p,
+                const label faceNo
+            ) const;
+
+            bool testPointFaceDistance
+            (
+                const vector& p,
+                const labelList& faceToTest,
+                const vectorList& points,
+                const vector& faceC,
+                const vector& faceA
+            ) const;
+
+            bool testEdgeEdgeDistance
+            (
+                const edge& eI,
+                const edge& eJ
+            ) const;
+
+            bool testEdgeEdgeDistance
+            (
+                const edge& eI,
+                const vector& eJs,
+                const vector& eJe
+            ) const;
+
+            const labelList realCellsInRangeOfSegment
+            (
+                const labelList& segmentFaces,
+                const labelList& segmentEdges,
+                const labelList& segmentPoints
+            ) const;
+
+            const labelList referredCellsInRangeOfSegment
+            (
+                const List<referredCell>& referredInteractionList,
+                const labelList& segmentFaces,
+                const labelList& segmentEdges,
+                const labelList& segmentPoints
+            ) const;
+
+        // Access
+
+            inline const polyMesh& mesh() const;
+
+            inline const directInteractionList& dil() const;
+
+            inline const referredCellList& ril() const;
+            inline referredCellList& ril();
+
+            inline const List<sendingReferralList>&
+                cellSendingReferralLists() const;
+
+            inline const List<receivingReferralList>&
+                cellReceivingReferralLists() const;
+
+            inline label nInteractingProcs() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "interactionListsI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/interactionListsI.H

@@ -22,46 +22,57 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-\*----------------------------------------------------------------------------*/
+\*---------------------------------------------------------------------------*/
-
-#include "moleculeCloud.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-void Foam::moleculeCloud::removeHighEnergyOverlaps()
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-{
-    Info << nl << "Removing high energy overlaps, removal order:";
 
-    forAll(removalOrder_, rO)
+const Foam::polyMesh& Foam::interactionLists::mesh() const
 {
-        Info << " " << pairPotentials_.idList()[removalOrder_[rO]];
+    return mesh_;
 }
 
-    Info << nl ;
 
-    label initialSize = this->size();
+const Foam::directInteractionList& Foam::interactionLists::dil() const
-
-    buildCellOccupancy();
-
-    iterator mol(this->begin());
-
-#   include "moleculeCloudRemoveHighEnergyOverlapsRealCells.H"
-
-    buildCellOccupancy();
-
-#   include "moleculeCloudRemoveHighEnergyOverlapsReferredCells.H"
-
-    buildCellOccupancy();
-
-    label molsRemoved = initialSize - this->size();
-
-    if (Pstream::parRun())
 {
-        reduce(molsRemoved, sumOp<label>());
+    return dil_;
 }
 
-    Info << tab << molsRemoved << " molecules removed" << endl;
+
+inline const Foam::referredCellList& Foam::interactionLists::ril() const
+{
+    return ril_;
 }
 
 
+inline Foam::referredCellList& Foam::interactionLists::ril()
+{
+    return ril_;
+}
+
+
+inline const Foam::List<Foam::sendingReferralList>&
+    Foam::interactionLists::cellSendingReferralLists() const
+{
+    return cellSendingReferralLists_;
+}
+
+
+inline const Foam::List<Foam::receivingReferralList>&
+    Foam::interactionLists::cellReceivingReferralLists() const
+{
+    return cellReceivingReferralLists_;
+}
+
+
+inline Foam::label Foam::interactionLists::nInteractingProcs() const
+{
+    return cellReceivingReferralLists_.size();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/receivingReferralList.C

@@ -26,19 +26,17 @@ License
 
 #include "receivingReferralList.H"
 
-namespace Foam
-{
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-receivingReferralList::receivingReferralList()
+Foam::receivingReferralList::receivingReferralList()
 :
     labelListList(),
     sourceProc_(-1)
 {}
 
 
-receivingReferralList::receivingReferralList
+Foam::receivingReferralList::receivingReferralList
 (
     const label sourceProc,
     const labelListList& refCellsToSendTo
@@ -49,7 +47,7 @@ receivingReferralList::receivingReferralList
 {}
 
 
-receivingReferralList::receivingReferralList
+Foam::receivingReferralList::receivingReferralList
 (
     const receivingReferralList& rL
 )
@@ -61,13 +59,13 @@ receivingReferralList::receivingReferralList
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-receivingReferralList::~receivingReferralList()
+Foam::receivingReferralList::~receivingReferralList()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-void receivingReferralList::operator=(const receivingReferralList& rhs)
+void Foam::receivingReferralList::operator=(const receivingReferralList& rhs)
 {
     // Check for assignment to self
     if (this == &rhs)
@@ -91,8 +89,8 @@ void receivingReferralList::operator=(const receivingReferralList& rhs)
 
 bool operator==
 (
-    const receivingReferralList& a,
+    const Foam::receivingReferralList& a,
-    const receivingReferralList& b
+    const Foam::receivingReferralList& b
 )
 {
     // Trivial reject: lists are different size
@@ -107,11 +105,11 @@ bool operator==
         return false;
     }
 
-    List<bool> fnd(a.size(), false);
+    Foam::List<bool> fnd(a.size(), false);
 
     forAll (b, bI)
     {
-        labelList curLList = b[bI];
+        Foam::labelList curLList = b[bI];
 
         bool found = false;
 
@@ -143,7 +141,7 @@ bool operator==
 }
 
 
-Istream& operator>>(Istream& is, receivingReferralList& rRL)
+Foam::Istream& Foam::operator>>(Istream& is, receivingReferralList& rRL)
 {
     is >> rRL.sourceProc_ >> static_cast<labelListList&>(rRL);
 
@@ -156,7 +154,7 @@ Istream& operator>>(Istream& is, receivingReferralList& rRL)
 }
 
 
-Ostream& operator<<
+Foam::Ostream& Foam::operator<<
 (
     Ostream& os,
     const receivingReferralList& rRL
@@ -175,7 +173,5 @@ Ostream& operator<<
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace Foam
-
 
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/receivingReferralListI.H

@@ -26,12 +26,9 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline label receivingReferralList::sourceProc() const
+inline Foam::label Foam::receivingReferralList::sourceProc() const
 {
     return sourceProc_;
 }
@@ -41,8 +38,8 @@ inline label receivingReferralList::sourceProc() const
 
 inline bool operator!=
 (
-    const receivingReferralList& a,
+    const Foam::receivingReferralList& a,
-    const receivingReferralList& b
+    const Foam::receivingReferralList& b
 )
 {
     return (!(a == b));
@@ -51,6 +48,5 @@ inline bool operator!=
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace Foam
 
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/sendingReferralList.C

@@ -26,19 +26,16 @@ License
 
 #include "sendingReferralList.H"
 
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-sendingReferralList::sendingReferralList()
+Foam::sendingReferralList::sendingReferralList()
 :
     labelList(),
     destinationProc_(-1)
 {}
 
 
-sendingReferralList::sendingReferralList
+Foam::sendingReferralList::sendingReferralList
 (
     const label destinationProc,
     const labelList& cellsToSend
@@ -49,7 +46,7 @@ sendingReferralList::sendingReferralList
 {}
 
 
-sendingReferralList::sendingReferralList
+Foam::sendingReferralList::sendingReferralList
 (
     const sendingReferralList& rL
 )
@@ -61,13 +58,13 @@ sendingReferralList::sendingReferralList
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-sendingReferralList::~sendingReferralList()
+Foam::sendingReferralList::~sendingReferralList()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-void sendingReferralList::operator=(const sendingReferralList& rhs)
+void Foam::sendingReferralList::operator=(const sendingReferralList& rhs)
 {
     // Check for assignment to self
     if (this == &rhs)
@@ -90,8 +87,8 @@ void sendingReferralList::operator=(const sendingReferralList& rhs)
 
 bool operator==
 (
-    const sendingReferralList& a,
+    const Foam::sendingReferralList& a,
-    const sendingReferralList& b
+    const Foam::sendingReferralList& b
 )
 {
     // Trivial reject: lists are different size
@@ -106,11 +103,11 @@ bool operator==
         return false;
     }
 
-    List<bool> fnd(a.size(), false);
+    Foam::List<bool> fnd(a.size(), false);
 
     forAll (b, bI)
     {
-        label curLabel = b[bI];
+        Foam::label curLabel = b[bI];
 
         bool found = false;
 
@@ -142,7 +139,7 @@ bool operator==
 }
 
 
-Istream& operator>>
+Foam::Istream& Foam::operator>>
 (
     Istream& is,
     sendingReferralList& sRL
@@ -156,7 +153,7 @@ Istream& operator>>
 }
 
 
-Ostream& operator<<
+Foam::Ostream& Foam::operator<<
 (
     Ostream& os,
     const sendingReferralList& rL
@@ -173,6 +170,4 @@ Ostream& operator<<
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace Foam
-
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/sendingReferralListI.H

@@ -26,12 +26,9 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline label sendingReferralList::destinationProc() const
+inline Foam::label Foam::sendingReferralList::destinationProc() const
 {
     return destinationProc_;
 }
@@ -41,8 +38,8 @@ inline label sendingReferralList::destinationProc() const
 
 inline bool operator!=
 (
-    const sendingReferralList& a,
+    const Foam::sendingReferralList& a,
-    const sendingReferralList& b
+    const Foam::sendingReferralList& b
 )
 {
     return (!(a == b));
@@ -51,6 +48,4 @@ inline bool operator!=
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace Foam
-
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/referredCell/referredCell.C

@@ -25,7 +25,7 @@ License
 \*----------------------------------------------------------------------------*/
 
 #include "referredCell.H"
-#include "moleculeCloud.H"
+#include "interactionLists.H"
 
 namespace Foam
 {
@@ -330,7 +330,7 @@ vector referredCell::referPosition(const vector& positionToRefer) const
 }
 
 
-vectorList referredCell::referPositions
+vectorList referredCell::referPosition
 (
     const vectorList& positionsToRefer
 ) const
@@ -365,10 +365,16 @@ void referredCell::referInMols(const List<referredMolecule>& incomingMols)
                 referPosition
                 (
                     incomingMols[iM].position()
+                ),
+                referPosition
+                (
+                    incomingMols[iM].sitePositions()
                 )
             )
         );
     }
+
+    shrink();
 }
 
 
@@ -378,7 +384,7 @@ bool referredCell::duplicate(const referredCell& refCellDupl) const
     (
         sourceProc_ == refCellDupl.sourceProc()
         && sourceCell_ == refCellDupl.sourceCell()
-        && mag(offset_ - refCellDupl.offset()) < moleculeCloud::transTol
+        && mag(offset_ - refCellDupl.offset()) < interactionLists::transTol
     );
 }
 
@@ -389,7 +395,7 @@ bool referredCell::duplicate(const label procNo,const label nCells) const
     (
         sourceProc_ == procNo
         && sourceCell_ < nCells
-        && mag(offset_) < moleculeCloud::transTol
+        && mag(offset_) < interactionLists::transTol
     );
 }
 

src/lagrangian/molecularDynamics/molecule/interactionLists/referredCell/referredCell.H

@@ -186,7 +186,7 @@ public:
 
         //- Use internal transformatation values to transform the given
         //  list of postions to their new locations.
-        vectorList referPositions(const vectorList& positionsToRefer) const;
+        vectorList referPosition(const vectorList& positionsToRefer) const;
 
         //- Use internal transformatation values to rotate the given vector
         vector rotateVector(const vector& vectorToRotate) const;

src/lagrangian/molecularDynamics/molecule/interactionLists/referredCellList/referredCellList.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -44,6 +44,9 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
+class interactionLists;
+
 /*---------------------------------------------------------------------------*\
                         Class referredCellList Declaration
 \*---------------------------------------------------------------------------*/
@@ -54,24 +57,30 @@ class referredCellList
 {
     // Private data
 
-        moleculeCloud& molCloud_;
+        const interactionLists& il_;
+
+        // Private Member Functions
+
+        void buildReferredCellList
+        (
+            bool pointPointListBuild
+        );
 
 
 public:
 
     // Constructors
 
-        //- Construct from moleculeCloud
+        //- Construct lists by searching the mesh
-        referredCellList(moleculeCloud& molCloud);
-
-        //- Construct from components
         referredCellList
         (
-            moleculeCloud& molCloud,
+            interactionLists& il,
-            const List<referredCell>& referredCells,
+            bool pointPointListBuild
-            const List<label>& realCells
         );
 
+        //- Construct from file
+        referredCellList (interactionLists& il);
+
 
     // Destructor
 
@@ -80,11 +89,9 @@ public:
 
     // Member Functions
 
-        void setRealCellsInRange();
+        void referMolecules(const List<DynamicList<molecule*> >& cellOccupancy);
 
-        void referMolecules();
+        inline const interactionLists& il() const;
-
-        inline const moleculeCloud& molCloud();
 };
 
 

src/lagrangian/molecularDynamics/molecule/interactionLists/referredCellList/referredCellListI.H

@@ -26,9 +26,9 @@ License
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-inline const Foam::moleculeCloud& Foam::referredCellList::molCloud()
+inline const Foam::interactionLists& Foam::referredCellList::il() const
 {
-    return molCloud_;
+    return il_;
 }
 
 

src/lagrangian/molecularDynamics/molecule/interactionLists/referredMolecule/referredMolecule.C

@@ -26,41 +26,40 @@ License
 
 #include "referredMolecule.H"
 
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-referredMolecule::referredMolecule()
+Foam::referredMolecule::referredMolecule()
 {}
 
 
-referredMolecule::referredMolecule
+Foam::referredMolecule::referredMolecule
 (
     const label id,
-    const vector& position
+    const vector& position,
+    const List<vector>& sitePositions
 )
 :
     id_(id),
-    position_(position)
+    position_(position),
+    sitePositions_(sitePositions)
 {}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-referredMolecule::~referredMolecule()
+Foam::referredMolecule::~referredMolecule()
 {}
 
 
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
-Istream& operator>>
+Foam::Istream& Foam::operator>>
 (
     Istream& is,
     referredMolecule& rM
 )
 {
-    is >> rM.id_ >> rM.position_;
+    is >> rM.id_ >> rM.position_ >> rM.sitePositions_;
 
     is.check("Istream& operator<<(Istream& f, const referredMolecule& sRL");
 
@@ -68,13 +67,15 @@ Istream& operator>>
 }
 
 
-Ostream& operator<<
+Foam::Ostream& Foam::operator<<
 (
     Ostream& os,
     const referredMolecule& rM
 )
 {
-    os << rM.id() << token::SPACE << rM.position();
+    os  << rM.id()
+        << token::SPACE << rM.position()
+        << token::SPACE << rM.sitePositions();
 
     os.check("Ostream& operator<<(Ostream& f, const referredMolecule& rM");
 
@@ -84,6 +85,4 @@ Ostream& operator<<
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace Foam
-
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/referredMolecule/referredMolecule.H

@@ -37,11 +37,13 @@ SourceFiles
 #define referredMolecule_H
 
 #include "vector.H"
+#include "List.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
+
 /*---------------------------------------------------------------------------*\
                            Class referredMolecule Declaration
 \*---------------------------------------------------------------------------*/
@@ -54,6 +56,8 @@ class referredMolecule
 
         vector position_;
 
+        List<vector> sitePositions_;
+
 
 public:
 
@@ -66,7 +70,8 @@ public:
         referredMolecule
         (
             const label id,
-            const vector& position
+            const vector& position,
+            const List<vector>& sitePositions
         );
 
 
@@ -83,6 +88,8 @@ public:
 
             inline const vector& position() const;
 
+            inline const List<vector>& sitePositions() const;
+
 
     // Friend Operators
 

src/lagrangian/molecularDynamics/molecule/interactionLists/referredMolecule/referredMoleculeI.H

@@ -24,29 +24,33 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-namespace Foam
-{
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline label referredMolecule::id() const
+inline Foam::label Foam::referredMolecule::id() const
 {
     return id_;
 }
 
 
-inline const vector& referredMolecule::position() const
+inline const Foam::vector& Foam::referredMolecule::position() const
 {
     return position_;
 }
 
 
+inline const Foam::List<Foam::vector>&
+Foam::referredMolecule::sitePositions() const
+{
+    return sitePositions_;
+}
+
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
-inline bool operator==
+inline bool Foam::operator==
 (
-    const referredMolecule& a,
+    const Foam::referredMolecule& a,
-    const referredMolecule& b
+    const Foam::referredMolecule& b
 )
 {
     return
@@ -57,10 +61,10 @@ inline bool operator==
 }
 
 
-inline bool operator!=
+inline bool Foam::operator!=
 (
-    const referredMolecule& a,
+    const Foam::referredMolecule& a,
-    const referredMolecule& b
+    const Foam::referredMolecule& b
 )
 {
     return !(a == b);
@@ -69,7 +73,4 @@ inline bool operator!=
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace Foam
-
-
 // ************************************************************************* //

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

@@ -2,25 +2,25 @@
 
 List< scalarField > allSpeciesN_RU
 (
-    molecules.pairPotentials().nIds(),
+    molecules.potential().nIds(),
     scalarField (mesh.nCells(), 0.0)
 );
 
 List< scalarField > allSpeciesM_RU
 (
-    molecules.pairPotentials().nIds(),
+    molecules.potential().nIds(),
     scalarField (mesh.nCells(), 0.0)
 );
 
 List< vectorField > allSpeciesVelocitySum_RU
 (
-    molecules.pairPotentials().nIds(),
+    molecules.potential().nIds(),
     vectorField (mesh.nCells(), vector::zero)
 );
 
 List< scalarField > allSpeciesVelocityMagSquaredSum_RU
 (
-    molecules.pairPotentials().nIds(),
+    molecules.potential().nIds(),
     scalarField (mesh.nCells(), 0.0)
 );
 
@@ -34,13 +34,13 @@ Info << nl << "Creating fields." << endl;
 
 PtrList<volScalarField> allSpeciesRhoN
 (
-    molecules.pairPotentials().nIds()
+    molecules.potential().nIds()
 );
 
 forAll (allSpeciesRhoN, rN)
 {
     Info << "    Creating number density field for "
-        << molecules.pairPotentials().idList()[rN] << endl;
+        << molecules.potential().idList()[rN] << endl;
 
     allSpeciesRhoN.set
     (
@@ -49,7 +49,7 @@ forAll (allSpeciesRhoN, rN)
         (
             IOobject
             (
-                "rhoN_" + molecules.pairPotentials().idList()[rN],
+                "rhoN_" + molecules.potential().idList()[rN],
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,
@@ -89,13 +89,13 @@ totalRhoN.correctBoundaryConditions();
 
 PtrList<volScalarField> allSpeciesRhoM
 (
-    molecules.pairPotentials().nIds()
+    molecules.potential().nIds()
 );
 
 forAll (allSpeciesRhoM, rM)
 {
     Info << "    Creating mass density field for "
-        << molecules.pairPotentials().idList()[rM] << endl;
+        << molecules.potential().idList()[rM] << endl;
 
     allSpeciesRhoM.set
     (
@@ -104,7 +104,7 @@ forAll (allSpeciesRhoM, rM)
         (
             IOobject
             (
-                "rhoM_" + molecules.pairPotentials().idList()[rM],
+                "rhoM_" + molecules.potential().idList()[rM],
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,
@@ -144,13 +144,13 @@ totalRhoM.correctBoundaryConditions();
 
 PtrList<volVectorField> allSpeciesVelocity
 (
-    molecules.pairPotentials().nIds()
+    molecules.potential().nIds()
 );
 
 forAll (allSpeciesVelocity, v)
 {
     Info << "    Creating velocity field for "
-        << molecules.pairPotentials().idList()[v] << endl;
+        << molecules.potential().idList()[v] << endl;
 
     allSpeciesVelocity.set
     (
@@ -159,7 +159,7 @@ forAll (allSpeciesVelocity, v)
         (
             IOobject
             (
-                "velocity_" + molecules.pairPotentials().idList()[v],
+                "velocity_" + molecules.potential().idList()[v],
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,
@@ -218,13 +218,13 @@ volVectorField totalVelocity
 
 PtrList<volScalarField> allSpeciesTemperature
 (
-    molecules.pairPotentials().nIds()
+    molecules.potential().nIds()
 );
 
 forAll (allSpeciesTemperature, t)
 {
     Info << "    Creating temperature field for "
-        << molecules.pairPotentials().idList()[t] << endl;
+        << molecules.potential().idList()[t] << endl;
 
     allSpeciesTemperature.set
     (
@@ -233,7 +233,7 @@ forAll (allSpeciesTemperature, t)
         (
             IOobject
             (
-                "temperature_" + molecules.pairPotentials().idList()[t],
+                "temperature_" + molecules.potential().idList()[t],
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,
@@ -274,13 +274,13 @@ totalTemperature.correctBoundaryConditions();
 
 PtrList<volScalarField> allSpeciesMeanKE
 (
-    molecules.pairPotentials().nIds()
+    molecules.potential().nIds()
 );
 
 forAll (allSpeciesMeanKE, mKE)
 {
     Info << "    Creating mean kinetic energy field for "
-        << molecules.pairPotentials().idList()[mKE] << endl;
+        << molecules.potential().idList()[mKE] << endl;
 
     allSpeciesMeanKE.set
     (
@@ -289,7 +289,7 @@ forAll (allSpeciesMeanKE, mKE)
         (
             IOobject
             (
-                "meanKE_" + molecules.pairPotentials().idList()[mKE],
+                "meanKE_" + molecules.potential().idList()[mKE],
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,

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

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

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

@@ -32,18 +32,26 @@ Description
 \*---------------------------------------------------------------------------*/
 
 
-vector singleStepTotalMomentum(vector::zero);
+vector singleStepTotalLinearMomentum(vector::zero);
+
+vector singleStepTotalAngularMomentum(vector::zero);
 
 scalar singleStepMaxVelocityMag = 0.0;
 
 scalar singleStepTotalMass = 0.0;
 
-scalar singleStepTotalKE = 0.0;
+scalar singleStepTotalLinearKE = 0.0;
+
+scalar singleStepTotalAngularKE = 0.0;
 
 scalar singleStepTotalPE = 0.0;
 
 scalar singleStepTotalrDotf = 0.0;
 
+//vector singleStepCentreOfMass(vector::zero);
+
+label singleStepNMols = molecules.size();
+
 {
     IDLList<molecule>::iterator mol(molecules.begin());
 
@@ -54,20 +62,52 @@ scalar singleStepTotalrDotf = 0.0;
         ++mol
     )
     {
-        const scalar molM(mol().mass());
+        label molId = mol().id();
 
-        const vector& molU(mol().U());
+        scalar molMass(molecules.constProps(molId).mass());
 
-        singleStepTotalMomentum += molU * molM;
+        singleStepTotalMass += molMass;
 
-        singleStepTotalMass += molM;
+        //singleStepCentreOfMass += mol().position()*molMass;
-
-        if(mag(molU) > singleStepMaxVelocityMag)
-        {
-            singleStepMaxVelocityMag = mag(molU);
     }
 
-        singleStepTotalKE += 0.5*molM*magSqr(molU);
+    // if(singleStepNMols)
+    // {
+    //     singleStepCentreOfMass /= singleStepTotalMass;
+    // }
+
+    for
+    (
+        mol = molecules.begin();
+        mol != molecules.end();
+        ++mol
+    )
+    {
+        label molId = mol().id();
+
+        scalar molMass(molecules.constProps(molId).mass());
+
+        const diagTensor& molMoI(molecules.constProps(molId).momentOfInertia());
+
+        const vector& molV(mol().v());
+
+        const vector& molOmega(inv(molMoI) & mol().pi());
+
+        vector molPiGlobal = mol().Q() & mol().pi();
+
+        singleStepTotalLinearMomentum += molV * molMass;
+
+        singleStepTotalAngularMomentum += molPiGlobal;
+        //+((mol().position() - singleStepCentreOfMass) ^ (molV * molMass));
+
+        if(mag(molV) > singleStepMaxVelocityMag)
+        {
+            singleStepMaxVelocityMag = mag(molV);
+        }
+
+        singleStepTotalLinearKE += 0.5*molMass*magSqr(molV);
+
+        singleStepTotalAngularKE += 0.5*(molOmega & molMoI & molOmega);
 
         singleStepTotalPE += mol().potentialEnergy();
 
@@ -75,17 +115,19 @@ scalar singleStepTotalrDotf = 0.0;
     }
 }
 
-label singleStepNMols = molecules.size();
-
 if (Pstream::parRun())
 {
-    reduce(singleStepTotalMomentum, sumOp<vector>());
+    reduce(singleStepTotalLinearMomentum, sumOp<vector>());
+
+    reduce(singleStepTotalAngularMomentum, sumOp<vector>());
 
     reduce(singleStepMaxVelocityMag, maxOp<scalar>());
 
     reduce(singleStepTotalMass, sumOp<scalar>());
 
-    reduce(singleStepTotalKE, sumOp<scalar>());
+    reduce(singleStepTotalLinearKE, sumOp<scalar>());
+
+    reduce(singleStepTotalAngularKE, sumOp<scalar>());
 
     reduce(singleStepTotalPE, sumOp<scalar>());
 
@@ -99,19 +141,31 @@ if (singleStepNMols)
     Info<< "Number of mols in system = "
         << singleStepNMols << nl
         << "Overall number density = "
-            << singleStepNMols/meshVolume << " m^-3" << nl
+        << singleStepNMols/meshVolume << nl
         << "Overall mass density = "
-            << singleStepTotalMass/meshVolume << " kg/m^3" << nl
+        << singleStepTotalMass/meshVolume << nl
-        << "Average velocity per mol = "
+        << "Average linear momentum per mol = "
-            << singleStepTotalMomentum/singleStepTotalMass << " m/s" << nl
+        << singleStepTotalLinearMomentum/singleStepNMols << ' '
+        << mag(singleStepTotalLinearMomentum)/singleStepNMols << nl
+        << "Average angular momentum per mol = "
+        << singleStepTotalAngularMomentum << ' '
+        << mag(singleStepTotalAngularMomentum)/singleStepNMols << nl
         << "Maximum |velocity| = "
-            << singleStepMaxVelocityMag << " m/s" << nl
+        << singleStepMaxVelocityMag << nl
-        << "Average KE per mol = "
+        << "Average linear KE per mol = "
-            << singleStepTotalKE/singleStepNMols << " J" << nl
+        << singleStepTotalLinearKE/singleStepNMols << nl
+        << "Average angular KE per mol = "
+        << singleStepTotalAngularKE/singleStepNMols << nl
         << "Average PE per mol = "
-            << singleStepTotalPE/singleStepNMols << " J" << nl
+        << singleStepTotalPE/singleStepNMols << nl
         << "Average TE per mol = "
-            << (singleStepTotalKE + singleStepTotalPE)/singleStepNMols << " J"
+        <<
+        (
+            singleStepTotalLinearKE
+            + singleStepTotalAngularKE
+            + singleStepTotalPE
+        )
+        /singleStepNMols
         << endl;
 
 //     Info << singleStepNMols << " "

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

@@ -2,25 +2,25 @@ if (runTime.outputTime())
 {
     allSpeciesN_RU = List< scalarField >
     (
-        molecules.pairPotentials().nIds(),
+        molecules.potential().nIds(),
         scalarField (mesh.nCells(), 0.0)
     );
 
     allSpeciesM_RU = List< scalarField >
     (
-        molecules.pairPotentials().nIds(),
+        molecules.potential().nIds(),
         scalarField (mesh.nCells(), 0.0)
     );
 
     allSpeciesVelocitySum_RU = List< vectorField >
     (
-        molecules.pairPotentials().nIds(),
+        molecules.potential().nIds(),
         vectorField (mesh.nCells(), vector::zero)
     );
 
     allSpeciesVelocityMagSquaredSum_RU = List< scalarField >
     (
-        molecules.pairPotentials().nIds(),
+        molecules.potential().nIds(),
         scalarField (mesh.nCells(), 0.0)
     );
 }

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

@@ -33,11 +33,13 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-accumulatedTotalMomentum += singleStepTotalMomentum;
+accumulatedTotalLinearMomentum += singleStepTotalLinearMomentum;
 
 accumulatedTotalMass += singleStepTotalMass;
 
-accumulatedTotalKE += singleStepTotalKE;
+accumulatedTotalLinearKE += singleStepTotalLinearKE;
+
+accumulatedTotalAngularKE += singleStepTotalAngularKE;
 
 accumulatedTotalPE += singleStepTotalPE;
 
@@ -55,12 +57,12 @@ if (runTime.outputTime())
 
         averageTemperature =
         (
-            2.0/(3.0 * moleculeCloud::kb * accumulatedNMols)
+            2.0/(6.0 * moleculeCloud::kb * accumulatedNMols)
             *
             (
-                accumulatedTotalKE
+                accumulatedTotalLinearKE + accumulatedTotalAngularKE
                 -
-                0.5*magSqr(accumulatedTotalMomentum)/accumulatedTotalMass
+                0.5*magSqr(accumulatedTotalLinearMomentum)/accumulatedTotalMass
             )
         );
 
@@ -82,7 +84,7 @@ if (runTime.outputTime())
         Info << "----------------------------------------" << nl
             << "Averaged properties" << nl
             << "Average |velocity| = "
-                << mag(accumulatedTotalMomentum)/accumulatedTotalMass
+                << mag(accumulatedTotalLinearMomentum)/accumulatedTotalMass
                 << " m/s" << nl
             << "Average temperature = "
                 << averageTemperature << " K" << nl
@@ -98,11 +100,13 @@ if (runTime.outputTime())
 
     // reset counters
 
-    accumulatedTotalMomentum = vector::zero;
+    accumulatedTotalLinearMomentum = vector::zero;
 
     accumulatedTotalMass = 0.0;
 
-    accumulatedTotalKE = 0.0;
+    accumulatedTotalLinearKE = 0.0;
+
+    accumulatedTotalAngularKE = 0.0;
 
     accumulatedTotalPE = 0.0;
 

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

@@ -30,11 +30,13 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-vector accumulatedTotalMomentum(vector::zero);
+vector accumulatedTotalLinearMomentum(vector::zero);
 
 scalar accumulatedTotalMass = 0.0;
 
-scalar accumulatedTotalKE = 0.0;
+scalar accumulatedTotalAngularKE = 0.0;
+
+scalar accumulatedTotalLinearKE = 0.0;
 
 scalar accumulatedTotalPE = 0.0;
 

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

@@ -25,82 +25,182 @@ License
 \*----------------------------------------------------------------------------*/
 
 #include "moleculeCloud.H"
+#include "molecule.H"
 #include "Random.H"
 #include "Time.H"
 
-namespace Foam
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::tensor Foam::molecule::rotationTensorX(scalar phi) const
 {
+    return tensor
+    (
+        1, 0, 0,
+        0, Foam::cos(phi), -Foam::sin(phi),
+        0, Foam::sin(phi), Foam::cos(phi)
+    );
+}
 
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-bool molecule::move(molecule::trackData& td)
+Foam::tensor Foam::molecule::rotationTensorY(scalar phi) const
+{
+    return tensor
+    (
+        Foam::cos(phi), 0, Foam::sin(phi),
+        0, 1, 0,
+        -Foam::sin(phi), 0, Foam::cos(phi)
+    );
+}
+
+
+Foam::tensor Foam::molecule::rotationTensorZ(scalar phi) const
+{
+    return tensor
+    (
+        Foam::cos(phi), -Foam::sin(phi), 0,
+        Foam::sin(phi), Foam::cos(phi), 0,
+        0, 0, 1
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::molecule::trackData::trackData
+(
+    moleculeCloud& molCloud,
+    label part
+)
+:
+    Particle<molecule>::trackData(molCloud),
+    molCloud_(molCloud),
+    part_(part)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+bool Foam::molecule::move(molecule::trackData& td)
 {
     td.switchProcessor = false;
     td.keepParticle = true;
 
+    const constantProperties& constProps(td.molCloud().constProps(id_));
+
     scalar deltaT = cloud().pMesh().time().deltaT().value();
+
+    if (td.part() == 0)
+    {
+        // First leapfrog velocity adjust part, required before tracking+force
+        // part
+
+        v_ += 0.5*deltaT*a_;
+
+        pi_ += 0.5*deltaT*tau_;
+    }
+    else if (td.part() == 1)
+    {
+        // Leapfrog tracking part
+
         scalar tEnd = (1.0 - stepFraction())*deltaT;
         scalar dtMax = tEnd;
 
-    moleculeCloud::integrationMethods method
-            = td.molCloud().integrationMethod();
-
-    if (method == moleculeCloud::imVerletLeapfrog)
-    {
-        if (td.part() == 1)  // Leapfrog 1st Part
-        {
-            if (stepFraction() < VSMALL)
-            {
-                U_ += 0.5*deltaT*A_;
-            }
-
         while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
         {
             // set the lagrangian time-step
             scalar dt = min(dtMax, tEnd);
 
-                dt *= trackToFace(position() + dt*U_, td);
+            dt *= trackToFace(position() + dt*v_, td);
 
             tEnd -= dt;
             stepFraction() = 1.0 - tEnd/deltaT;
         }
     }
-        else if (td.part() == 2)  // Leapfrog 2nd Part
+    else if (td.part() == 2)
     {
-            U_ += 0.5*deltaT*A_;
+        // Leapfrog orientation adjustment, carried out before force calculation
+        // but after tracking stage, i.e. rotation carried once linear motion
+        // complete.
+
+        if (!constProps.pointMolecule())
+        {
+            const diagTensor& momentOfInertia(constProps.momentOfInertia());
+
+            tensor R;
+
+            if (!constProps.linearMolecule())
+            {
+                R = rotationTensorX(0.5*deltaT*pi_.x()/momentOfInertia.xx());
+                pi_ = pi_ & R;
+                Q_ = Q_ & R;
+            }
+
+            R = rotationTensorY(0.5*deltaT*pi_.y()/momentOfInertia.yy());
+            pi_ = pi_ & R;
+            Q_ = Q_ & R;
+
+            R = rotationTensorZ(deltaT*pi_.z()/momentOfInertia.zz());
+            pi_ = pi_ & R;
+            Q_ = Q_ & R;
+
+            R = rotationTensorY(0.5*deltaT*pi_.y()/momentOfInertia.yy());
+            pi_ = pi_ & R;
+            Q_ = Q_ & R;
+
+            if (!constProps.linearMolecule())
+            {
+                R = rotationTensorX(0.5*deltaT*pi_.x()/momentOfInertia.xx());
+                pi_ = pi_ & R;
+                Q_ = Q_ & R;
+            }
+        }
+
+        setSitePositions(constProps);
+    }
+    else if (td.part() == 3)
+    {
+        // Second leapfrog velocity adjust part, required after tracking+force
+        // part
+
+        scalar m = constProps.mass();
+
+        a_ = vector::zero;
+
+        tau_ = vector::zero;
+
+        forAll(siteForces_, s)
+        {
+            const vector& f = siteForces_[s];
+
+            a_ += f/m;
+
+            tau_ += (constProps.siteReferencePositions()[s] ^ (Q_.T() & f));
+        }
+
+        v_ += 0.5*deltaT*a_;
+
+        pi_ += 0.5*deltaT*tau_;
+
+        if (constProps.pointMolecule())
+        {
+            tau_ = vector::zero;
+
+            pi_ = vector::zero;
+        }
+
+        if (constProps.linearMolecule())
+        {
+            tau_.x() = 0.0;
+
+            pi_.x() = 0.0;
+        }
     }
     else
     {
         FatalErrorIn("molecule::move(molecule::trackData& td)") << nl
             << td.part()
-                << " is an invalid part of integration method: "
+            << " is an invalid part of the integration method."
-                << method << nl
-                << abort(FatalError);
-        }
-    }
-    else if (method == moleculeCloud::imPredictorCorrector)
-    {
-        if (td.part() == 1) // Predictor Part
-        {
-
-        }
-        else if (td.part() == 2) // Corrector Part
-        {
-
-        }
-        else
-        {
-            FatalErrorIn("molecule::move(molecule::trackData& td)") << nl
-                << td.part() << " is an invalid part of integration method: "
-                << method
-                << abort(FatalError);
-        }
-    }
-    else
-    {
-        FatalErrorIn("molecule::move(molecule::trackData& td)") << nl
-                << "Unknown integration method: "
-                << method
             << abort(FatalError);
     }
 
@@ -108,20 +208,40 @@ bool molecule::move(molecule::trackData& td)
 }
 
 
-void molecule::transformProperties(const tensor& T)
+void Foam::molecule::transformProperties(const tensor& T)
-{}
+{
+    Q_ = T & Q_;
+
+    sitePositions_ = position_ + (T & (sitePositions_ - position_));
+
+    siteForces_ = T & siteForces_;
+}
 
 
-void molecule::transformProperties(const vector& separation)
+void Foam::molecule::transformProperties(const vector& separation)
 {
-    if (tethered_)
+    if (special_ == SPECIAL_TETHERED)
     {
-        tetherPosition_ += separation;
+        specialPosition_ += separation;
     }
 }
 
 
-void molecule::hitProcessorPatch
+void Foam::molecule::setSitePositions(const constantProperties& constProps)
+{
+    sitePositions_ = position_ + (Q_ & constProps.siteReferencePositions());
+}
+
+
+void Foam::molecule::setSiteSizes(label size)
+{
+    sitePositions_.setSize(size);
+
+    siteForces_.setSize(size);
+}
+
+
+void Foam::molecule::hitProcessorPatch
 (
     const processorPolyPatch&,
     molecule::trackData& td
@@ -131,7 +251,7 @@ void molecule::hitProcessorPatch
 }
 
 
-void molecule::hitProcessorPatch
+void Foam::molecule::hitProcessorPatch
 (
     const processorPolyPatch&,
     int&
@@ -139,7 +259,7 @@ void molecule::hitProcessorPatch
 {}
 
 
-void molecule::hitWallPatch
+void Foam::molecule::hitWallPatch
 (
     const wallPolyPatch& wpp,
     molecule::trackData& td
@@ -148,43 +268,17 @@ void molecule::hitWallPatch
     vector nw = wpp.faceAreas()[wpp.whichFace(face())];
     nw /= mag(nw);
 
-    scalar Un = U_ & nw;
+    scalar vn = v_ & nw;
-//     vector Ut = U_ - Un*nw;
 
-//     Random rand(clock::getTime());
-
-//     scalar tmac = 0.8;
-
-//     scalar wallTemp = 2.5;
-
-//     if (rand.scalar01() < tmac)
-//     {
-//         // Diffuse reflection
-//
-//         vector tw1 = Ut/mag(Ut);
-//
-//         vector tw2 = nw ^ tw1;
-//
-//         U_ = sqrt(wallTemp/mass_)*rand.GaussNormal()*tw1
-//                 + sqrt(wallTemp/mass_)*rand.GaussNormal()*tw2
-//                 - mag(sqrt(wallTemp/mass_)*rand.GaussNormal())*nw;
-//     }
-
-//     else
-//     {
     // Specular reflection
-
+    if (vn > 0)
-        if (Un > 0)
     {
-            U_ -= 2*Un*nw;
+        v_ -= 2*vn*nw;
     }
-
-//     }
-
 }
 
 
-void molecule::hitWallPatch
+void Foam::molecule::hitWallPatch
 (
     const wallPolyPatch&,
     int&
@@ -192,7 +286,7 @@ void molecule::hitWallPatch
 {}
 
 
-void molecule::hitPatch
+void Foam::molecule::hitPatch
 (
     const polyPatch&,
     molecule::trackData& td
@@ -202,14 +296,12 @@ void molecule::hitPatch
 }
 
 
-void molecule::hitPatch
+void Foam::molecule::hitPatch
 (
     const polyPatch&,
     int&
 )
 {}
 
-} // End namespace Foam
-
 
 // ************************************************************************* //

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

@@ -41,6 +41,7 @@ SourceFiles
 #include "Particle.H"
 #include "IOstream.H"
 #include "autoPtr.H"
+#include "diagTensor.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -58,54 +59,98 @@ class molecule
 :
     public Particle<molecule>
 {
-    // Private data
-
-        //- Be careful with the ordering of data.
-        //  It has an impact on binary transfer:
-        //    -# Put the largest data members 1st
-        //    -# Pair up labels,
-        //    -# Don't go scalar-label, scalar-label, because in 64bit mode,
-        //       the labels will be padded by 4bytes.
-
-        // - mass of molecule
-        scalar mass_;
-
-        // - Velocity of molecule
-        vector U_;
-
-        // - Acceleration of molecule
-        vector A_;
-
-        // - Tether position
-        vector tetherPosition_;
-
-        // - Potential energy that this molecules posseses
-        scalar potentialEnergy_;
-
-        // - r_ij f_ij, stress dyad
-        tensor rf_;
-
-        // - Is the molecule tethered?
-        label tethered_;
-
-        // - id (type) of molecule
-        label id_;
-
 
 public:
 
-    friend class Cloud<molecule>;
+    // Values of special that are less than zero are for built-in functionality.
+    // Values greater than zero are user specifiable/expandable (i.e. test
+    // special_ >= SPECIAL_USER)
+
+    enum specialTypes
+    {
+        SPECIAL_TETHERED = -1,
+        SPECIAL_FROZEN   = -2,
+        NOT_SPECIAL      = 0,
+        SPECIAL_USER     = 1
+    };
+
+    //- Class to hold molecule constant properties
+    class constantProperties
+    {
+
+        // Private data
+
+            Field<vector> siteReferencePositions_;
+
+            List<scalar> siteCharges_;
+
+            List<label> siteIds_;
+
+            List<bool> pairPotentialSites_;
+
+            List<bool> electrostaticSites_;
+
+            diagTensor momentOfInertia_;
+
+            scalar mass_;
+
+        // Private Member Functions
+
+            void checkSiteListSizes() const;
+
+            void setInteracionSiteBools
+            (
+                const List<word>& siteIds,
+                const List<word>& pairPotSiteIds
+            );
+
+            bool linearMoleculeTest() const;
+
+    public:
+
+        inline constantProperties();
+
+        //- Construct from dictionary
+        inline constantProperties(const dictionary& dict);
+
+        // Member functions
+
+        inline const Field<vector>& siteReferencePositions() const;
+
+        inline const List<scalar>& siteCharges() const;
+
+        inline const List<label>& siteIds() const;
+
+        inline List<label>& siteIds();
+
+        inline const List<bool>& pairPotentialSites() const;
+
+        inline bool pairPotentialSite(label sId) const;
+
+        inline const List<bool>& electrostaticSites() const;
+
+        inline bool electrostaticSite(label sId) const;
+
+        inline const diagTensor& momentOfInertia() const;
+
+        inline bool linearMolecule() const;
+
+        inline bool pointMolecule() const;
+
+        inline scalar mass() const;
+
+        inline label nSites() const;
+    };
+
 
     //- Class used to pass tracking data to the trackToFace function
     class trackData
     :
         public Particle<molecule>::trackData
     {
-        //- Reference to the cloud containing this particle
         moleculeCloud& molCloud_;
 
         // label specifying which part of the integration algorithm is taking
-        // place (i.e. leapfrog 1 or leapfrog 2. Predictor or Corrector)
         label part_;
 
 
@@ -113,7 +158,7 @@ public:
 
         // Constructors
 
-            inline trackData
+            trackData
             (
                 moleculeCloud& molCloud,
                 label part
@@ -126,6 +171,55 @@ public:
             inline label part() const;
     };
 
+
+private:
+
+    // Private data
+
+        //- Be careful with the ordering of data.
+        //  It has an impact on binary transfer:
+        //    -# Put the largest data members 1st
+        //    -# Pair up labels,
+        //    -# Don't go scalar-label, scalar-label, because in 64bit mode,
+        //       the labels will be padded by 4bytes.
+
+        tensor Q_;
+
+        vector v_;
+
+        vector a_;
+
+        vector pi_;
+
+        vector tau_;
+
+        vector specialPosition_;
+
+        scalar potentialEnergy_;
+
+        // - r_ij f_ij, stress dyad
+        tensor rf_;
+
+        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:
+
+    friend class Cloud<molecule>;
+
     // Constructors
 
         //- Construct from components
@@ -134,11 +228,14 @@ public:
             const Cloud<molecule>& c,
             const vector& position,
             const label celli,
-            const scalar mass,
+            const tensor& Q,
-            const vector& U,
+            const vector& v,
-            const vector& A,
+            const vector& a,
-            const vector& tetherPosition,
+            const vector& pi,
-            const label tethered,
+            const vector& tau,
+            const vector& specialPosition,
+            const constantProperties& constProps,
+            const label special,
             const label id
         );
 
@@ -159,45 +256,54 @@ public:
 
     // Member Functions
 
+        //- Tracking
+            bool move(trackData&);
+
             void transformProperties(const tensor& T);
 
             void transformProperties(const vector& separation);
 
+            void setSitePositions(const constantProperties& constProps);
+
+            void setSiteSizes(label size);
+
         // Access
 
-            //- Return id
+            inline const tensor& Q() const;
-            inline label id() const;
+            inline tensor& Q();
 
-            //- Return mass
+            inline const vector& v() const;
-            inline scalar mass() const;
+            inline vector& v();
 
-            //- Return velocity
+            inline const vector& a() const;
-            inline const vector& U() const;
+            inline vector& a();
-            inline vector& U();
 
-            //- Return acceleration
+            inline const vector& pi() const;
-            inline const vector& A() const;
+            inline vector& pi();
-            inline vector& A();
+
+            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();
 
-            //- Return potential energy
             inline scalar potentialEnergy() const;
             inline scalar& potentialEnergy();
 
-            //- Return stress contribution
             inline const tensor& rf() const;
             inline tensor& rf();
 
-            //- Return tethered
+            inline label special() const;
-            inline label tethered() const;
 
-            //- Return tetherPosition
+            inline bool tethered() const;
-            inline const vector& tetherPosition() const;
-            inline vector& tetherPosition();
-
-
-        //- Tracking
-        bool move(trackData&);
 
+            inline label id() const;
 
     // Member Operators
 

src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H

@@ -24,142 +24,574 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-inline molecule::molecule
+inline Foam::molecule::constantProperties::constantProperties()
+:
+    siteReferencePositions_(Field<vector>(0)),
+    siteCharges_(List<scalar>(0)),
+    siteIds_(List<label>(0)),
+    pairPotentialSites_(List<bool>(false)),
+    electrostaticSites_(List<bool>(false)),
+    momentOfInertia_(diagTensor(0, 0, 0)),
+    mass_(0)
+{}
+
+
+inline Foam::molecule::constantProperties::constantProperties
+(
+    const dictionary& dict
+)
+:
+    siteReferencePositions_(dict.lookup("siteReferencePositions")),
+    siteCharges_(dict.lookup("siteCharges")),
+    siteIds_(List<word>(dict.lookup("siteIds")).size(), -1),
+    pairPotentialSites_(),
+    electrostaticSites_(),
+    momentOfInertia_(),
+    mass_()
+{
+    checkSiteListSizes();
+
+    setInteracionSiteBools
+    (
+        List<word>(dict.lookup("siteIds")),
+        List<word>(dict.lookup("pairPotentialSiteIds"))
+    );
+
+    scalarList siteMasses(dict.lookup("siteMasses"));
+
+    mass_ = sum(siteMasses);
+
+    vector centreOfMass(vector::zero);
+
+    // Calculate the centre of mass of the body and subtract it from each
+    // position
+
+    forAll(siteReferencePositions_, i)
+    {
+        centreOfMass += siteReferencePositions_[i]*siteMasses[i];
+    }
+
+    centreOfMass /= mass_;
+
+    siteReferencePositions_ -= centreOfMass;
+
+    if(siteIds_.size() == 1)
+    {
+        // Single site molecule - no rotational motion.
+
+        siteReferencePositions_[0] = vector::zero;
+
+        momentOfInertia_ = diagTensor(-1, -1, -1);
+    }
+    else if(linearMoleculeTest())
+    {
+        // Linear molecule.
+
+        Info<< nl << "Linear molecule." << endl;
+
+        vector dir = siteReferencePositions_[1] - siteReferencePositions_[0];
+
+        dir /= mag(dir);
+
+        tensor Q = rotationTensor(dir, vector(1,0,0));
+
+        siteReferencePositions_ = (Q & siteReferencePositions_);
+
+        // The rotation was around the centre of mass but remove any
+        // components that have crept in due to floating point errors
+
+        centreOfMass = vector::zero;
+
+        forAll(siteReferencePositions_, i)
+        {
+            centreOfMass += siteReferencePositions_[i]*siteMasses[i];
+        }
+
+        centreOfMass /= mass_;
+
+        siteReferencePositions_ -= centreOfMass;
+
+        diagTensor momOfInertia = diagTensor::zero;
+
+        forAll(siteReferencePositions_, i)
+        {
+            const vector& p(siteReferencePositions_[i]);
+
+            momOfInertia += siteMasses[i]*diagTensor
+            (
+                0, p.x()*p.x(), p.x()*p.x()
+            );
+        }
+
+        momentOfInertia_ = diagTensor
+         (
+             -1,
+             momOfInertia.yy(),
+             momOfInertia.zz()
+         );
+    }
+    else
+    {
+        // Fully 6DOF molecule
+
+        // Calculate the inertia tensor in the current orientation
+
+        tensor momOfInertia(tensor::zero);
+
+        forAll(siteReferencePositions_, i)
+        {
+            const vector& p(siteReferencePositions_[i]);
+
+            momOfInertia += siteMasses[i]*tensor
+            (
+                p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(),
+                -p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(),
+                -p.z()*p.x(), -p.z()*p.y(), p.x()*p.x() + p.y()*p.y()
+            );
+        }
+
+        if (eigenValues(momOfInertia).x() < VSMALL)
+        {
+            FatalErrorIn("molecule::constantProperties::constantProperties")
+                << "An eigenvalue of the inertia tensor is zero, the molecule "
+                << dict.name()
+                << " is not a valid 6DOF shape."
+                << nl << abort(FatalError);
+        }
+
+        // Normalise the inertia tensor magnitude to avoid SMALL numbers in the
+        // components causing problems
+
+        momOfInertia /= eigenValues(momOfInertia).x();
+
+        tensor e = eigenVectors(momOfInertia);
+
+        // Calculate the transformation between the principle axes to the
+        // global axes
+
+        tensor Q =
+        vector(1,0,0)*e.x() + vector(0,1,0)*e.y() + vector(0,0,1)*e.z();
+
+        // Transform the site positions
+
+        siteReferencePositions_ = (Q & siteReferencePositions_);
+
+        // Recalculating the moment of inertia with the new site positions
+
+        // The rotation was around the centre of mass but remove any
+        // components that have crept in due to floating point errors
+
+        centreOfMass = vector::zero;
+
+        forAll(siteReferencePositions_, i)
+        {
+            centreOfMass += siteReferencePositions_[i]*siteMasses[i];
+        }
+
+        centreOfMass /= mass_;
+
+        siteReferencePositions_ -= centreOfMass;
+
+        // Calculate the moment of inertia in the principle component
+        // reference frame
+
+        momOfInertia = tensor::zero;
+
+        forAll(siteReferencePositions_, i)
+        {
+            const vector& p(siteReferencePositions_[i]);
+
+            momOfInertia += siteMasses[i]*tensor
+            (
+                p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(),
+                -p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(),
+                -p.z()*p.x(), -p.z()*p.y(), p.x()*p.x() + p.y()*p.y()
+            );
+        }
+
+        momentOfInertia_ = diagTensor
+        (
+            momOfInertia.xx(),
+            momOfInertia.yy(),
+            momOfInertia.zz()
+        );
+    }
+}
+
+
+inline Foam::molecule::molecule
 (
     const Cloud<molecule>& c,
     const vector& position,
     const label celli,
-    const scalar mass,
+    const tensor& Q,
-    const vector& U,
+    const vector& v,
-    const vector& A,
+    const vector& a,
-    const vector& tetherPosition,
+    const vector& pi,
-    const label tethered,
+    const vector& tau,
+    const vector& specialPosition,
+    const constantProperties& constProps,
+    const label special,
     const label id
+
 )
 :
     Particle<molecule>(c, position, celli),
-    mass_(mass),
+    Q_(Q),
-    U_(U),
+    v_(v),
-    A_(A),
+    a_(a),
-    tetherPosition_(tetherPosition),
+    pi_(pi),
+    tau_(tau),
+    specialPosition_(specialPosition),
     potentialEnergy_(0.0),
     rf_(tensor::zero),
-    tethered_(tethered),
+    special_(special),
-    id_(id)
+    id_(id),
-{}
+    siteForces_(constProps.nSites(), vector::zero),
-
+    sitePositions_(constProps.nSites())
-
-inline molecule::trackData::trackData
-(
-    moleculeCloud& molCloud,
-    label part
-)
-:
-    Particle<molecule>::trackData(refCast<Cloud<molecule> >(molCloud)),
-    molCloud_(molCloud),
-    part_(part)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-inline label molecule::id() const
 {
-    return id_;
+    setSitePositions(constProps);
 }
 
 
-inline scalar molecule::mass() const
+// * * * 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);
+    }
+}
+
+
+inline void Foam::molecule::constantProperties::setInteracionSiteBools
+(
+    const List<word>& siteIds,
+    const List<word>& pairPotSiteIds
+)
+{
+    pairPotentialSites_.setSize(siteIds_.size());
+
+    electrostaticSites_.setSize(siteIds_.size());
+
+    forAll(siteIds_, i)
+    {
+        const word& id(siteIds[i]);
+
+        pairPotentialSites_[i] = (findIndex(pairPotSiteIds, id) > -1);
+
+        electrostaticSites_[i] = (mag(siteCharges_[i]) > VSMALL);
+    }
+}
+
+
+inline bool Foam::molecule::constantProperties::linearMoleculeTest() const
+{
+    if (siteIds_.size() == 2)
+    {
+        return true;
+    }
+
+    vector refDir = siteReferencePositions_[1] - siteReferencePositions_[0];
+
+    refDir /= mag(refDir);
+
+    for
+    (
+        label i = 2;
+        i < siteReferencePositions_.size();
+        i++
+    )
+    {
+        vector dir = siteReferencePositions_[i] - siteReferencePositions_[i-1];
+
+        dir /= mag(dir);
+
+        if (mag(refDir & dir) < 1 - SMALL)
+        {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+
+// * * * * * * * constantProperties Member Functions * * * * * * * * * * * * //
+
+inline const Foam::Field<Foam::vector>&
+Foam::molecule::constantProperties::siteReferencePositions() const
+{
+    return siteReferencePositions_;
+}
+
+
+inline const Foam::List<Foam::scalar>&
+Foam::molecule::constantProperties::siteCharges() const
+{
+    return siteCharges_;
+}
+
+
+inline const Foam::List<Foam::label>&
+Foam::molecule::constantProperties::siteIds() const
+{
+    return siteIds_;
+}
+
+
+inline Foam::List<Foam::label>&
+Foam::molecule::constantProperties::siteIds()
+{
+    return siteIds_;
+}
+
+
+inline const Foam::List<bool>&
+Foam::molecule::constantProperties::pairPotentialSites() const
+{
+    return pairPotentialSites_;
+}
+
+
+inline bool Foam::molecule::constantProperties::pairPotentialSite
+(
+    label sId
+) const
+{
+    label s = findIndex(siteIds_, sId);
+
+    if(s == -1)
+    {
+        FatalErrorIn("moleculeI.H") << nl
+            << sId << " site not found."
+            << nl << abort(FatalError);
+    }
+
+    return pairPotentialSites_[s];
+}
+
+
+inline const Foam::List<bool>&
+Foam::molecule::constantProperties::electrostaticSites() const
+{
+    return electrostaticSites_;
+}
+
+inline bool Foam::molecule::constantProperties::electrostaticSite
+(
+    label sId
+) const
+{
+    label s = findIndex(siteIds_, sId);
+
+    if(s == -1)
+    {
+        FatalErrorIn("moleculeI.H") << nl
+            << sId << " site not found."
+            << nl << abort(FatalError);
+    }
+
+    return electrostaticSites_[s];
+}
+
+
+inline const Foam::diagTensor&
+Foam::molecule::constantProperties::momentOfInertia() const
+{
+    return momentOfInertia_;
+}
+
+
+inline bool Foam::molecule::constantProperties::linearMolecule() const
+{
+    return (momentOfInertia_.xx() < 0);
+}
+
+
+inline bool Foam::molecule::constantProperties::pointMolecule() const
+{
+    return (momentOfInertia_.zz() < 0);
+}
+
+
+inline Foam::scalar Foam::molecule::constantProperties::mass() const
 {
     return mass_;
 }
 
 
-inline const vector& molecule::U() const
+inline Foam::label Foam::molecule::constantProperties::nSites() const
 {
-    return U_;
+    return siteIds_.size();
+
 }
 
+// * * * * * * * * * * *  trackData Member Functions * * * * * * * * * * * * //
 
-inline vector& molecule::U()
+inline Foam::moleculeCloud& Foam::molecule::trackData::molCloud()
-{
-    return U_;
-}
-
-
-inline const vector& molecule::A() const
-{
-    return A_;
-}
-
-
-inline vector& molecule::A()
-{
-    return A_;
-}
-
-
-inline scalar molecule::potentialEnergy() const
-{
-    return potentialEnergy_;
-}
-
-
-inline scalar& molecule::potentialEnergy()
-{
-    return potentialEnergy_;
-}
-
-
-inline const tensor& molecule::rf() const
-{
-    return rf_;
-}
-
-
-inline tensor& molecule::rf()
-{
-    return rf_;
-}
-
-
-inline label molecule::tethered() const
-{
-    return tethered_;
-}
-
-
-inline const vector& molecule::tetherPosition() const
-{
-    return tetherPosition_;
-}
-
-
-inline vector& molecule::tetherPosition()
-{
-    return tetherPosition_;
-}
-
-
-inline moleculeCloud& molecule::trackData::molCloud()
 {
     return molCloud_;
 }
 
 
-inline label molecule::trackData::part() const
+inline Foam::label Foam::molecule::trackData::part() const
 {
     return part_;
 }
 
 
+// * * * * * * * * * * * * molecule Member Functions * * * * * * * * * * * * //
+
+inline const Foam::tensor& Foam::molecule::Q() const
+{
+    return Q_;
+}
+
+
+inline Foam::tensor& Foam::molecule::Q()
+{
+    return Q_;
+}
+
+
+inline const Foam::vector& Foam::molecule::v() const
+{
+    return v_;
+}
+
+
+inline Foam::vector& Foam::molecule::v()
+{
+    return v_;
+}
+
+
+inline const Foam::vector& Foam::molecule::a() const
+{
+    return a_;
+}
+
+
+inline Foam::vector& Foam::molecule::a()
+{
+    return a_;
+}
+
+
+inline const Foam::vector& Foam::molecule::pi() const
+{
+    return pi_;
+}
+
+
+inline Foam::vector& Foam::molecule::pi()
+{
+    return pi_;
+}
+
+
+inline const Foam::vector& Foam::molecule::tau() const
+{
+    return tau_;
+}
+
+
+inline Foam::vector& Foam::molecule::tau()
+{
+    return tau_;
+}
+
+
+inline const Foam::List<Foam::vector>& Foam::molecule::siteForces() const
+{
+    return siteForces_;
+}
+
+
+inline Foam::List<Foam::vector>& Foam::molecule::siteForces()
+{
+    return siteForces_;
+}
+
+
+inline const Foam::List<Foam::vector>& Foam::molecule::sitePositions() const
+{
+    return sitePositions_;
+}
+
+
+inline Foam::List<Foam::vector>& Foam::molecule::sitePositions()
+{
+    return sitePositions_;
+}
+
+
+inline const Foam::vector& Foam::molecule::specialPosition() const
+{
+    return specialPosition_;
+}
+
+
+inline Foam::vector& Foam::molecule::specialPosition()
+{
+    return specialPosition_;
+}
+
+
+inline Foam::scalar Foam::molecule::potentialEnergy() const
+{
+    return potentialEnergy_;
+}
+
+
+inline Foam::scalar& Foam::molecule::potentialEnergy()
+{
+    return potentialEnergy_;
+}
+
+
+inline const Foam::tensor& Foam::molecule::rf() const
+{
+    return rf_;
+}
+
+
+inline Foam::tensor& Foam::molecule::rf()
+{
+    return rf_;
+}
+
+
+inline Foam::label Foam::molecule::special() const
+{
+    return special_;
+}
+
+
+inline bool Foam::molecule::tethered() const
+{
+    return special_ == SPECIAL_TETHERED;
+}
+
+
+inline Foam::label Foam::molecule::id() const
+{
+    return id_;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace Foam
-
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C

@@ -26,7 +26,6 @@ License
 
 #include "molecule.H"
 #include "IOstreams.H"
-
 #include "moleculeCloud.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -38,137 +37,164 @@ Foam::molecule::molecule
     bool readFields
 )
 :
-    Particle<molecule>(cloud, is)
+    Particle<molecule>(cloud, is, true),
+    Q_(tensor::zero),
+    v_(vector::zero),
+    a_(vector::zero),
+    pi_(vector::zero),
+    tau_(vector::zero),
+    specialPosition_(vector::zero),
+    potentialEnergy_(0.0),
+    rf_(tensor::zero),
+    special_(0),
+    id_(0),
+    siteForces_(0),
+    sitePositions_(0)
 {
     if (readFields)
     {
         if (is.format() == IOstream::ASCII)
         {
-            id_ = readLabel(is);
+            is >> Q_;
-            mass_ = readScalar(is);
+            is >> v_;
-            is >> U_;
+            is >> a_;
-            is >> A_;
+            is >> pi_;
-            is >> potentialEnergy_;
+            is >> tau_;
+            is >> siteForces_;
+            is >> sitePositions_;
+            is >> specialPosition_;
+            potentialEnergy_ = readScalar(is);
             is >> rf_;
-            is >> tethered_;
+            special_ = readLabel(is);
-            is >> tetherPosition_;
+            id_ = readLabel(is);
         }
         else
         {
             is.read
             (
-                reinterpret_cast<char*>(&mass_),
+                reinterpret_cast<char*>(&Q_),
-                sizeof(mass_)
+                sizeof(Q_)
-                + sizeof(U_)
+                + sizeof(v_)
-                + sizeof(A_)
+                + sizeof(a_)
-                + sizeof(tetherPosition_)
+                + sizeof(pi_)
+                + sizeof(tau_)
+                + sizeof(specialPosition_)
                 + sizeof(potentialEnergy_)
                 + sizeof(rf_)
-                + sizeof(tethered_)
+                + sizeof(special_)
                 + sizeof(id_)
             );
+
+            is >> siteForces_ >> sitePositions_;
         }
     }
 
     // Check state of Istream
-    is.check("Foam::molecule::molecule(Foam::Istream&)");
+    is.check
+    (
+        "Foam::molecule::molecule"
+        "(const Cloud<molecule>& cloud, Foam::Istream&), bool"
+    );
 }
 
 
-namespace Foam
+void Foam::molecule::readFields(moleculeCloud& mC)
-{
-
-void molecule::readFields(moleculeCloud& mC)
 {
     if (!mC.size())
     {
         return;
     }
 
+    IOField<tensor> Q(mC.fieldIOobject("Q", IOobject::MUST_READ));
+    mC.checkFieldIOobject(mC, Q);
+
+    IOField<vector> v(mC.fieldIOobject("v", IOobject::MUST_READ));
+    mC.checkFieldIOobject(mC, v);
+
+    IOField<vector> a(mC.fieldIOobject("a", IOobject::MUST_READ));
+    mC.checkFieldIOobject(mC, a);
+
+    IOField<vector> pi(mC.fieldIOobject("pi", IOobject::MUST_READ));
+    mC.checkFieldIOobject(mC, pi);
+
+    IOField<vector> tau(mC.fieldIOobject("tau", IOobject::MUST_READ));
+    mC.checkFieldIOobject(mC, tau);
+
+    IOField<vector> specialPosition
+    (
+        mC.fieldIOobject("specialPosition", IOobject::MUST_READ)
+    );
+    mC.checkFieldIOobject(mC, specialPosition);
+
+    IOField<label> special(mC.fieldIOobject("special", IOobject::MUST_READ));
+    mC.checkFieldIOobject(mC, special);
+
     IOField<label> id(mC.fieldIOobject("id", IOobject::MUST_READ));
     mC.checkFieldIOobject(mC, id);
 
-    IOField<scalar> mass(mC.fieldIOobject("mass", IOobject::MUST_READ));
-    mC.checkFieldIOobject(mC, mass);
-
-    IOField<vector> U(mC.fieldIOobject("U", IOobject::MUST_READ));
-    mC.checkFieldIOobject(mC, U);
-
-    IOField<vector> A(mC.fieldIOobject("A", IOobject::MUST_READ));
-    mC.checkFieldIOobject(mC, A);
-
-    IOField<label> tethered(mC.fieldIOobject("tethered", IOobject::MUST_READ));
-    mC.checkFieldIOobject(mC, tethered);
-
-    IOField<vector> tetherPositions
-    (
-        mC.fieldIOobject("tetherPositions", IOobject::MUST_READ)
-    );
-    mC.checkFieldIOobject(mC, tetherPositions);
-
     label i = 0;
     forAllIter(moleculeCloud, mC, iter)
     {
         molecule& mol = iter();
 
+        mol.Q_ = Q[i];
+        mol.v_ = v[i];
+        mol.a_ = a[i];
+        mol.pi_ = pi[i];
+        mol.tau_ = tau[i];
+        mol.specialPosition_ = specialPosition[i];
+        mol.special_ = special[i];
         mol.id_ = id[i];
-        mol.mass_ = mass[i];
-        mol.U_ = U[i];
-        mol.A_ = A[i];
-        mol.potentialEnergy_ = 0.0;
-        mol.rf_ = tensor::zero;
-        mol.tethered_ = tethered[i];
-        mol.tetherPosition_ = tetherPositions[i];
         i++;
     }
 }
 
 
-void molecule::writeFields(const moleculeCloud& mC)
+void Foam::molecule::writeFields(const moleculeCloud& mC)
 {
     Particle<molecule>::writeFields(mC);
 
     label np = mC.size();
 
+    IOField<tensor> Q(mC.fieldIOobject("Q", IOobject::NO_READ), np);
+    IOField<vector> v(mC.fieldIOobject("v", IOobject::NO_READ), np);
+    IOField<vector> a(mC.fieldIOobject("a", IOobject::NO_READ), np);
+    IOField<vector> pi(mC.fieldIOobject("pi", IOobject::NO_READ), np);
+    IOField<vector> tau(mC.fieldIOobject("tau", IOobject::NO_READ), np);
+    IOField<vector> specialPosition
+    (
+        mC.fieldIOobject("specialPosition", IOobject::NO_READ),
+        np
+    );
+    IOField<label> special(mC.fieldIOobject("special", IOobject::NO_READ), np);
     IOField<label> id(mC.fieldIOobject("id", IOobject::NO_READ), np);
-    IOField<scalar> mass(mC.fieldIOobject("mass", IOobject::NO_READ), np);
-    IOField<vector> U(mC.fieldIOobject("U", IOobject::NO_READ), np);
-    IOField<vector> A(mC.fieldIOobject("A", IOobject::NO_READ), np);
-    IOField<label> tethered
-    (
-        mC.fieldIOobject("tethered", IOobject::NO_READ),
-        np
-    );
-    IOField<vector> tetherPositions
-    (
-        mC.fieldIOobject("tetherPositions", IOobject::NO_READ),
-        np
-    );
 
     label i = 0;
     forAllConstIter(moleculeCloud, mC, iter)
     {
         const molecule& mol = iter();
 
+        Q[i] = mol.Q_;
+        v[i] = mol.v_;
+        a[i] = mol.a_;
+        pi[i] = mol.pi_;
+        tau[i] = mol.tau_;
+        specialPosition[i] = mol.specialPosition_;
+        special[i] = mol.special_;
         id[i] = mol.id_;
-        mass[i] = mol.mass_;
-        U[i] = mol.U_;
-        A[i] = mol.A_;
-        tethered[i] = mol.tethered_;
-        tetherPositions[i] = mol.tetherPosition_;
         i++;
     }
 
+    Q.write();
+    v.write();
+    a.write();
+    pi.write();
+    tau.write();
+    specialPosition.write();
+    special.write();
     id.write();
-    mass.write();
-    U.write();
-    A.write();
-    tethered.write();
-    tetherPositions.write();
 }
 
-};  // end of namespace Foam
-
 
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
@@ -176,33 +202,40 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const molecule& mol)
 {
     if (os.format() == IOstream::ASCII)
     {
-        os  << mol.id_
+        os  << token::SPACE << static_cast<const Particle<molecule>&>(mol)
-                << token::SPACE << mol.mass_
-                << token::SPACE << static_cast<const Particle<molecule>&>(mol)
             << token::SPACE << mol.face()
             << token::SPACE << mol.stepFraction()
-                << token::SPACE << mol.U_
+            << token::SPACE << mol.Q_
-                << token::SPACE << mol.A_
+            << token::SPACE << mol.v_
+            << token::SPACE << mol.a_
+            << token::SPACE << mol.pi_
+            << token::SPACE << mol.tau_
+            << token::SPACE << mol.specialPosition_
             << token::SPACE << mol.potentialEnergy_
             << token::SPACE << mol.rf_
-                << token::SPACE << mol.tethered_
+            << token::SPACE << mol.special_
-                << token::SPACE << mol.tetherPosition_;
+            << token::SPACE << mol.id_
+            << token::SPACE << mol.siteForces_
+            << token::SPACE << mol.sitePositions_;
     }
     else
     {
         os  << static_cast<const Particle<molecule>&>(mol);
         os.write
         (
-            reinterpret_cast<const char*>(&mol.mass_),
+            reinterpret_cast<const char*>(&mol.Q_),
-            sizeof(mol.mass_)
+            sizeof(mol.Q_)
-            + sizeof(mol.U_)
+            + sizeof(mol.v_)
-            + sizeof(mol.A_)
+            + sizeof(mol.a_)
-            + sizeof(mol.tetherPosition_)
+            + sizeof(mol.pi_)
+            + sizeof(mol.tau_)
+            + sizeof(mol.specialPosition_)
             + sizeof(mol.potentialEnergy_)
             + sizeof(mol.rf_)
-            + sizeof(mol.tethered_)
+            + sizeof(mol.special_)
             + sizeof(mol.id_)
         );
+        os << mol.siteForces_ << mol.sitePositions_;
     }
 
     // Check state of Ostream

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

@@ -28,18 +28,8 @@ Class
 Description
 
 SourceFiles
-    moleculeCloudApplyConstraintsAndThermostats.C
+    moleculeCloudI.H
-    moleculeCloudBuildCellInteractionLists.C
-    moleculeCloudBuildCellOccupancy.C
-    moleculeCloudBuildCellReferralLists.C
     moleculeCloud.C
-    moleculeCloudCalculateAndAccumulateProperties.C
-    moleculeCloudCalculateExternalForce.C
-    moleculeCloudCalculateForce.C
-    moleculeCloudCalculatePairForce.C
-    moleculeCloudCalculateTetherForce.C
-    moleculeCloudIntegrateEquationsOfMotion.C
-    moleculeCloudRemoveHighEnergyOverlaps.C
 
 
 \*---------------------------------------------------------------------------*/
@@ -50,14 +40,10 @@ SourceFiles
 #include "Cloud.H"
 #include "molecule.H"
 #include "IOdictionary.H"
-#include "vector2D.H"
+#include "potential.H"
-
+#include "interactionLists.H"
-#include "pairPotentialList.H"
+#include "labelVector.H"
-#include "tetherPotentialList.H"
+#include "Random.H"
-
-#include "receivingReferralList.H"
-#include "sendingReferralList.H"
-#include "referredCellList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,59 +59,90 @@ class moleculeCloud
     public Cloud<molecule>
 {
 
-public:
-
-    enum integrationMethods
-    {
-        imVerletLeapfrog,
-        imPredictorCorrector
-    };
-
-
 private:
 
     // Private data
 
         const polyMesh& mesh_;
 
-            // MD solution parameters
+        const potential& pot_;
-
-            static const NamedEnum<integrationMethods, 2>
-                integrationMethodNames_;
-
-            integrationMethods integrationMethod_;
-
-            scalar potentialEnergyLimit_;
-
-            labelList removalOrder_;
-
-            labelListList directInteractionList_;
-
-            referredCellList referredInteractionList_;
-
-            labelList realCellsWithinRCutMaxOfAnyReferringPatch_;
-
-            labelList realFacesWithinRCutMaxOfAnyReferringPatch_;
-
-            labelList realEdgesWithinRCutMaxOfAnyReferringPatch_;
-
-            labelList realPointsWithinRCutMaxOfAnyReferringPatch_;
-
-            List<sendingReferralList> cellSendingReferralLists_;
-
-            List<receivingReferralList> cellReceivingReferralLists_;
-
-            pairPotentialList pairPotentials_;
-
-            tetherPotentialList tetherPotentials_;
-
-            vector gravity_;
 
         List<DynamicList<molecule*> > cellOccupancy_;
 
+        interactionLists il_;
+
+        List<molecule::constantProperties> constPropList_;
+
+        Random rndGen_;
 
     // Private Member Functions
 
+        void buildConstProps();
+
+        void setSiteSizesAndPositions();
+
+        //- Determine which molecules are in which cells
+        void buildCellOccupancy();
+
+        void calculatePairForce();
+
+        inline void evaluatePair
+        (
+            molecule* molI,
+            molecule* molJ
+        );
+
+        inline void evaluatePair
+        (
+            molecule* molReal,
+            referredMolecule* molRef
+        );
+
+        inline bool evaluatePotentialLimit
+        (
+            molecule* molI,
+            molecule* molJ
+        ) const;
+
+        inline bool evaluatePotentialLimit
+        (
+            molecule* molReal,
+            referredMolecule* molRef
+        ) const;
+
+        void calculateTetherForce();
+
+        void calculateExternalForce();
+
+        void removeHighEnergyOverlaps();
+
+        void initialiseMolecules
+        (
+            const IOdictionary& mdInitialiseDict
+        );
+
+        void createMolecule
+        (
+            const point& position,
+            label cell,
+            label id,
+            bool tethered,
+            scalar temperature,
+            const vector& bulkVelocity
+        );
+
+        inline vector equipartitionLinearVelocity
+        (
+            scalar temperature,
+            scalar mass
+        );
+
+        inline vector equipartitionAngularMomentum
+        (
+            scalar temperature,
+            const molecule::constantProperties& cP
+        );
+
         //- Disallow default bitwise copy construct
         moleculeCloud(const moleculeCloud&);
 
@@ -137,157 +154,36 @@ public:
 
     // Static data members
 
-        //- Tolerance for checking that faces on a patch segment
-        static scalar transTol;
-
         static scalar kb;
 
+        static scalar elementaryCharge;
+
+        static scalar vacuumPermittivity;
 
     // Constructors
 
-        //- Construct given mesh
+        //- Construct given mesh and potential references
-        moleculeCloud(const polyMesh&);
-
-        //- Construct given polyMesh and fields of position, cell, mass,
-        //- id, U ands A.  Intended for use by the molConfig utility
         moleculeCloud
         (
             const polyMesh& mesh,
-            label nMol,
+            const potential& pot
-            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
 
-        // Access
 
-            inline const polyMesh& mesh() const;
+        //- Evolve the molecules (move, calculate forces, control state etc)
+        void evolve();
 
-        // MD solution parameters
+        void calculateForce();
-
-            inline const integrationMethods& integrationMethod() const;
-
-            inline scalar potentialEnergyLimit() const;
-
-            inline const labelList& removalOrder() const;
-
-            inline label nPairPotentials() const;
-
-            inline const labelListList& directInteractionList() const;
-
-            inline const referredCellList& referredInteractionList() const;
-
-            inline const labelList&
-                realCellsWithinRCutMaxOfAnyReferringPatch() const;
-
-            inline const labelList&
-                realFacesWithinRCutMaxOfAnyReferringPatch() const;
-
-            inline const labelList&
-                realEdgesWithinRCutMaxOfAnyReferringPatch() const;
-
-            inline const labelList&
-                realPointsWithinRCutMaxOfAnyReferringPatch() const;
-
-            inline const List<sendingReferralList>&
-                cellSendingReferralLists() const;
-
-            inline const List<receivingReferralList>&
-                cellReceivingReferralLists() const;
-
-            inline label nInteractingProcs() const;
-
-            inline const pairPotentialList& pairPotentials() const;
-
-            inline const tetherPotentialList& tetherPotentials() const;
-
-            inline const vector& gravity() const;
-
-            inline const List< DynamicList<molecule*> >& cellOccupancy() const;
-
-            void buildCellInteractionLists();
-
-            //- Build referralLists which define how to send information
-            // to referredCells to source cells
-            void buildCellReferralLists();
-
-            bool testPointFaceDistance
-            (
-                const label p,
-                const label faceNo
-            ) const;
-
-            bool testPointFaceDistance
-            (
-                const label p,
-                const referredCell& refCell
-            ) const;
-
-            bool testPointFaceDistance
-            (
-                const vectorList& pointsToTest,
-                const label faceNo
-            ) const;
-
-            bool testPointFaceDistance
-            (
-                const vector& p,
-                const label faceNo
-            ) const;
-
-            bool testPointFaceDistance
-            (
-                const vector& p,
-                const labelList& faceToTest,
-                const vectorList& points,
-                const vector& faceC,
-                const vector& faceA
-            ) const;
-
-            bool testEdgeEdgeDistance
-            (
-                const edge& eI,
-                const edge& eJ
-            ) const;
-
-            bool testEdgeEdgeDistance
-            (
-                const edge& eI,
-                const vector& eJs,
-                const vector& eJe
-            ) const;
-
-            const labelList realCellsInRangeOfSegment
-            (
-                const labelList& segmentFaces,
-                const labelList& segmentEdges,
-                const labelList& segmentPoints
-            ) const;
-
-            const labelList referredCellsInRangeOfSegment
-            (
-                const List<referredCell>& referredInteractionList,
-                const labelList& segmentFaces,
-                const labelList& segmentEdges,
-                const labelList& segmentPoints
-            ) const;
-
-            //- Determine which molecules are in which cells
-            void buildCellOccupancy();
-
-            //- Integrate the equations of motion using algorithm selected at
-            //  runtime from a dictionary. This will also call the function
-            //  to calculate the intermolecular forces (calculatePairForce()).
-
-            void integrateEquationsOfMotion();
 
         void applyConstraintsAndThermostats
         (
@@ -295,20 +191,27 @@ public:
             const scalar measuredTemperature
         );
 
-            void calculateForce();
+        // Access
 
-            void calculatePairForce();
+            inline const polyMesh& mesh() const;
 
-            void calculateTetherForce();
+            inline const potential& pot() const;
 
-            void calculateExternalForce();
+            inline const List<DynamicList<molecule*> >& cellOccupancy() const;
 
-            void removeHighEnergyOverlaps();
+            inline const interactionLists& il() const;
 
+            inline const List<molecule::constantProperties> constProps() const;
+
+            inline const molecule::constantProperties&
+            constProps(label id) const;
+
+            inline Random& rndGen();
 
     // Member Operators
 
         //- Write fields
+
             void writeFields() const;
 };
 

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

@@ -1,60 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 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 "moleculeCloud.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-void Foam::moleculeCloud::applyConstraintsAndThermostats
-(
-    const scalar targetTemperature,
-    const scalar measuredTemperature
-)
-{
-    scalar temperatureCorrectionFactor =
-        sqrt(targetTemperature/measuredTemperature);
-
-    Info<< "----------------------------------------" << nl
-        << "Temperature equilibration" << nl
-        << "Target temperature = "
-        << targetTemperature << nl
-        << "Measured temperature = "
-        << measuredTemperature << nl
-        << "Temperature correction factor ="
-        << temperatureCorrectionFactor << nl
-        << "----------------------------------------"
-        << endl;
-
-    iterator mol(this->begin());
-
-    for (mol = this->begin(); mol != this->end(); ++mol)
-    {
-        mol().U() *= temperatureCorrectionFactor;
-    }
-}
-
-
-// ************************************************************************* //

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

@@ -1,194 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 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 "moleculeCloud.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-void Foam::moleculeCloud::buildCellReferralLists()
-{
-    Info<< nl << "Determining molecule referring schedule" << endl;
-
-    const referredCellList& refIntL(referredInteractionList());
-
-    DynamicList<label> referralProcs;
-
-    // Run through all referredCells to build list of interacting processors
-
-    forAll(refIntL, rIL)
-    {
-        const referredCell& rC(refIntL[rIL]);
-
-        if (findIndex(referralProcs, rC.sourceProc()) == -1)
-        {
-            referralProcs.append(rC.sourceProc());
-        }
-    }
-
-    referralProcs.shrink();
-
-//     Pout << "referralProcs: " << nl << referralProcs << endl;
-
-    List<DynamicList<label> > cellSendingReferralLists(referralProcs.size());
-
-    List<DynamicList<DynamicList<label> > >
-        cellReceivingReferralLists(referralProcs.size());
-
-    // Run through all referredCells again building up send and receive info
-
-    forAll(refIntL, rIL)
-    {
-        const referredCell& rC(refIntL[rIL]);
-
-        label rPI = findIndex(referralProcs, rC.sourceProc());
-
-        DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
-
-        DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
-
-        label existingSource = findIndex(sRL, rC.sourceCell());
-
-        // Check to see if this source cell has already been allocated to
-        // come to this processor.  If not, add the source cell to the sending
-        // list and add the current referred cell to the receiving list.
-
-        // It shouldn't be possible for the sending and receiving lists to be
-        // different lengths, because their append operations happen at the
-        // same time.
-
-        if (existingSource == -1)
-        {
-            sRL.append(rC.sourceCell());
-
-            rRL.append
-            (
-                DynamicList<label> (labelList(1,rIL))
-            );
-        }
-        else
-        {
-            rRL[existingSource].append(rIL);
-
-            rRL[existingSource].shrink();
-        }
-    }
-
-    forAll(referralProcs, rPI)
-    {
-        DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
-
-        DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
-
-        sRL.shrink();
-
-        rRL.shrink();
-    }
-
-    // It is assumed that cell exchange is reciprocal, if proc A has cells to
-    // send to proc B, then proc B must have some to send to proc A.
-
-    cellReceivingReferralLists_.setSize(referralProcs.size());
-
-    cellSendingReferralLists_.setSize(referralProcs.size());
-
-    forAll(referralProcs, rPI)
-    {
-        DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
-
-        labelListList translLL(rRL.size());
-
-        forAll(rRL, rRLI)
-        {
-            translLL[rRLI] = rRL[rRLI];
-        }
-
-        cellReceivingReferralLists_[rPI] = receivingReferralList
-        (
-            referralProcs[rPI],
-            translLL
-        );
-    }
-
-    // Send sendingReferralLists to each interacting processor.
-
-    forAll(referralProcs, rPI)
-    {
-
-        DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
-
-        if (referralProcs[rPI] != Pstream::myProcNo())
-        {
-            if (Pstream::parRun())
-            {
-                OPstream toInteractingProc
-                (
-                    Pstream::blocking,
-                    referralProcs[rPI]
-                );
-
-                toInteractingProc << sendingReferralList
-                (
-                    Pstream::myProcNo(),
-                    sRL
-                );
-            }
-        }
-    }
-
-    // Receive sendingReferralLists from each interacting processor.
-
-    forAll(referralProcs, rPI)
-    {
-        if (referralProcs[rPI] != Pstream::myProcNo())
-        {
-            if (Pstream::parRun())
-            {
-                IPstream fromInteractingProc
-                (
-                    Pstream::blocking,
-                    referralProcs[rPI]
-                );
-
-                fromInteractingProc >> cellSendingReferralLists_[rPI];
-            }
-        }
-        else
-        {
-            cellSendingReferralLists_[rPI] = sendingReferralList
-            (
-                Pstream::myProcNo(),
-                cellSendingReferralLists[rPI]
-            );
-        }
-    }
-
-//     Pout << "receiving list: " << nl << cellReceivingReferralLists_ << endl;
-
-//     Pout << "sending list: " << nl << cellSendingReferralLists_ << endl;
-}
-
-
-// ************************************************************************* //

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

@@ -1,188 +0,0 @@
-Info << nl << "Building list of direct interaction neighbours" << endl;
-
-forAll (mesh_.points(), p)
-{
-    forAll(mesh_.faces(), f)
-    {
-        if(testPointFaceDistance(p, f))
-        {
-            const labelList& pCells(mesh_.pointCells()[p]);
-
-            const label cellO(mesh_.faceOwner()[f]);
-
-            const label cellN(mesh_.faceNeighbour()[f]);
-
-            forAll(pCells, pC)
-            {
-                const label cellI(pCells[pC]);
-
-                // cells are not added to their own DIL
-
-                if (cellO > cellI)
-                {
-                    if (findIndex(directInteractionList[cellI], cellO) == -1)
-                    {
-                        directInteractionList[cellI].append(cellO);
-                    }
-                }
-
-                if (cellI > cellO)
-                {
-                    if (findIndex(directInteractionList[cellO], cellI) == -1)
-                    {
-                        directInteractionList[cellO].append(cellI);
-                    }
-                }
-
-                if (mesh_.isInternalFace(f))
-                {
-                    // boundary faces will not have neighbour information
-
-                    if (cellN > cellI)
-                    {
-                        if
-                        (
-                             findIndex(directInteractionList[cellI], cellN)
-                          == -1
-                        )
-                        {
-                            directInteractionList[cellI].append(cellN);
-                        }
-                    }
-
-                    if (cellI > cellN)
-                    {
-                        if
-                        (
-                            findIndex(directInteractionList[cellN], cellI)
-                         == -1
-                        )
-                        {
-                            directInteractionList[cellN].append(cellI);
-                        }
-                    }
-                }
-            }
-        }
-    }
-}
-
-label edgeJIndex;
-
-forAll (mesh_.edges(), edgeIIndex)
-{
-    const edge& eI(mesh_.edges()[edgeIIndex]);
-
-    for
-    (
-        edgeJIndex = edgeIIndex + 1;
-        edgeJIndex != mesh_.edges().size();
-        ++edgeJIndex
-    )
-    {
-        const edge& eJ(mesh_.edges()[edgeJIndex]);
-
-        if (testEdgeEdgeDistance(eI, eJ))
-        {
-            const labelList& eICells(mesh_.edgeCells()[edgeIIndex]);
-
-            const labelList& eJCells(mesh_.edgeCells()[edgeJIndex]);
-
-            forAll(eICells, eIC)
-            {
-                const label cellI(eICells[eIC]);
-
-                forAll(eJCells, eJC)
-                {
-                    const label cellJ(eJCells[eJC]);
-
-                    if (cellJ > cellI)
-                    {
-                        if
-                        (
-                            findIndex(directInteractionList[cellI], cellJ)
-                         == -1
-                        )
-                        {
-                            directInteractionList[cellI].append(cellJ);
-                        }
-                    }
-
-                    if (cellI > cellJ)
-                    {
-                        if
-                        (
-                            findIndex(directInteractionList[cellJ], cellI)
-                         == -1
-                        )
-                        {
-                            directInteractionList[cellJ].append(cellI);
-                        }
-                    }
-                }
-            }
-        }
-    }
-}
-
-// label pointJIndex;
-//
-// forAll (mesh_.points(), pointIIndex)
-// {
-//     const point& ptI
-//     (
-//         mesh_.points()[pointIIndex]
-//     );
-//
-//     for
-//     (
-//         pointJIndex = pointIIndex;
-//         pointJIndex != mesh_.points().size();
-//         ++pointJIndex
-//     )
-//     {
-//         const point& ptJ
-//         (
-//             mesh_.points()[pointJIndex]
-//         );
-//
-//         if (magSqr(ptI - ptJ) <= rCutMaxSqr)
-//         {
-//             const labelList& ptICells
-//             (
-//                 mesh_.pointCells()[pointIIndex]
-//             );
-//
-//             const labelList& ptJCells
-//             (
-//                 mesh_.pointCells()[pointJIndex]
-//             );
-//
-//             forAll(ptICells, pIC)
-//             {
-//                 const label cellI(ptICells[pIC]);
-//
-//                 forAll(ptJCells, pJC)
-//                 {
-//                     const label cellJ(ptJCells[pJC]);
-//
-//                     if (cellJ > cellI)
-//                     {
-//                         if(findIndex(directInteractionList[cellI], cellJ) == -1)
-//                         {
-//                             directInteractionList[cellI].append(cellJ);
-//                         }
-//                     }
-//
-//                     if (cellI > cellJ)
-//                     {
-//                         if(findIndex(directInteractionList[cellJ], cellI) == -1)
-//                         {
-//                             directInteractionList[cellJ].append(cellI);
-//                         }
-//                     }
-//                 }
-//             }
-//         }
-//     }
-// }

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

@@ -1,42 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 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 "moleculeCloud.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-void Foam::moleculeCloud::calculateExternalForce()
-{
-    iterator mol(this->begin());
-
-    for (mol = this->begin(); mol != this->end(); ++mol)
-    {
-        mol().A() += gravity_;
-    }
-}
-
-
-// ************************************************************************* //

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

@@ -1,49 +0,0 @@
-
-vector rIJ;
-
-scalar rIJMag;
-
-scalar rIJMagSq;
-
-vector fIJ;
-
-label idI;
-
-label idJ;
-
-mol = this->begin();
-
-molecule* molI = &mol();
-
-molecule* molJ = &mol();
-
-forAll(directInteractionList_, dIL)
-{
-    forAll(cellOccupancy_[dIL],cellIMols)
-    {
-        molI = cellOccupancy_[dIL][cellIMols];
-
-        forAll(directInteractionList_[dIL], interactingCells)
-        {
-            List< molecule* > cellJ =
-                cellOccupancy_[directInteractionList_[dIL][interactingCells]];
-
-            forAll(cellJ, cellJMols)
-            {
-                molJ = cellJ[cellJMols];
-
-#               include "moleculeCloudCalculatePairForceRealCellsCalculationStep.H"
-            }
-        }
-
-        forAll(cellOccupancy_[dIL],cellIOtherMols)
-        {
-            molJ = cellOccupancy_[dIL][cellIOtherMols];
-
-            if (molJ > molI)
-            {
-#                include "moleculeCloudCalculatePairForceRealCellsCalculationStep.H"
-            }
-        }
-    }
-}

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

@@ -1,33 +0,0 @@
-idI = molI->id();
-
-idJ = molJ->id();
-
-rIJ = molI->position() - molJ->position();
-
-rIJMagSq = magSqr(rIJ);
-
-if(pairPotentials_.rCutSqr(idI, idJ, rIJMagSq))
-{
-    rIJMag = mag(rIJ);
-
-    fIJ = (rIJ/rIJMag)*pairPotentials_.force(idI, idJ, rIJMag);
-
-    // Acceleration increment for molI
-    molI->A() += fIJ/(molI->mass());
-
-    // Acceleration increment for molJ
-    molJ->A() += -fIJ/(molJ->mass());
-
-    scalar potentialEnergy
-    (
-        pairPotentials_.energy(idI, idJ, rIJMag)
-    );
-
-    molI->potentialEnergy() += 0.5*potentialEnergy;
-
-    molJ->potentialEnergy() += 0.5*potentialEnergy;
-
-    molI->rf() += rIJ * fIJ;
-
-    molJ->rf() += rIJ * fIJ;
-}

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

@@ -1,59 +0,0 @@
-vector rKL;
-
-scalar rKLMag;
-
-scalar rKLMagSq;
-
-vector fKL;
-
-label idK;
-
-label idL;
-
-molecule* molK = &mol();
-
-forAll(referredInteractionList_, rIL)
-{
-    const List<label>& realCells =
-        referredInteractionList_[rIL].realCellsForInteraction();
-
-    forAll(referredInteractionList_[rIL], refMols)
-    {
-        referredMolecule* molL = &(referredInteractionList_[rIL][refMols]);
-
-        forAll(realCells, rC)
-        {
-            List<molecule*> cellK = cellOccupancy_[realCells[rC]];
-
-            forAll(cellK, cellKMols)
-            {
-                molK = cellK[cellKMols];
-
-                idK = molK->id();
-
-                idL = molL->id();
-
-                rKL = molK->position() - molL->position();
-
-                rKLMagSq = magSqr(rKL);
-
-                if (pairPotentials_.rCutSqr(idK, idL, rKLMagSq))
-                {
-                    rKLMag = mag(rKL);
-
-                    fKL = (rKL/rKLMag)*pairPotentials_.force(idK, idL, rKLMag);
-
-                    // Acceleration increment for molK
-                    molK->A() += fKL/(molK->mass());
-
-                    // Adding a contribution of 1/2 of the potential energy
-                    // from this interaction
-                    molK->potentialEnergy() +=
-                        0.5*pairPotentials_.energy(idK, idL, rKLMag);
-
-                    molK->rf() += rKL * fKL;
-                }
-            }
-        }
-    }
-}

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

@@ -1,193 +0,0 @@
-// Parallel coding to access boundary information to build up interaction cell info
-
-// See preservePatchTypes for how to read the boundary file.
-
-// Read faceProcAddressing, as per reconstructPar, to get hold of the original,
-// undecomposed face label from a face on a processor mesh.  See email from Mattijs:
-
-// > Is it a case of reading the faceProcAddressing file, in the same way as
-// > something like reconstructPar?
-//         Correct.
-// 
-//         Note that faceProcAddressing is a bit weird since it also includes which side 
-// of an internal face we have. If I remember correctly:
-// 
-//         faceI == 0 illegal
-//         faceI > 0 we have the original owner of  faceI-1 i.e. we have the face in the 
-//         original order.
-//         faceI < 0 we have the original neighbour of -faceI-1 so the face is flipped.
-
-// Use the same functionality as
-// label polyBoundaryMesh::whichPatch(const label faceIndex) const
-// To determine which patch a face was on originally.
-
-if (Pstream::parRun())
-{
-//     if (Pstream::myProcNo() == Pstream::masterNo())
-// //     {
-//         dictionary patchDictionary;
-// 
-//         DynamicList<word> patchNames;
-// 
-//         {
-//             IOobject undecomposedBoundaryHeader
-//             (
-//                 "undecomposedBoundary",
-//                 mesh_.time().constant(),
-//                 polyMesh::meshSubDir,
-//                 mesh_,
-//                 IOobject::MUST_READ,
-//                 IOobject::NO_WRITE,
-//                 false
-//             );
-// 
-//             if (undecomposedBoundaryHeader.headerOk())
-//             {
-//                 polyBoundaryMeshEntries undecomposedPatchEntries
-//                 (
-//                     undecomposedBoundaryHeader
-//                 );
-// 
-//                 forAll(undecomposedPatchEntries, patchi)
-//                 {
-//                     patchNames.append
-//                     (
-//                         undecomposedPatchEntries[patchi].keyword()
-//                     );
-// 
-//                     patchDictionary.add
-//                     (
-//                         undecomposedPatchEntries[patchi]
-//                     );
-//                 }
-//             }
-//             else
-//             {
-//                 FatalErrorIn
-//                 (
-//                     "moleculeCloudBuildCellInteractionLists.C\n"
-//                 )
-//                     << "undecomposedBoundary file not found in "
-//                         "constant/polyMesh"
-//                     << abort(FatalError);
-//             }
-//         }
-// 
-//         labelIOList faceProcAddressing
-//         (
-//             IOobject
-//             (
-//                 "faceProcAddressing",
-//                 mesh_.time().constant(),
-//                 polyMesh::meshSubDir,
-//                 mesh_,
-//                 IOobject::MUST_READ,
-//                 IOobject::NO_WRITE,
-//                 false
-//             )
-//         );
-
-        labelList procPatches(mesh_.globalData().processorPatches());
-
-        forAll(procPatches,pP)
-        {
-            const processorPolyPatch& patch =
-                refCast<const processorPolyPatch>
-            (
-                mesh_.boundaryMesh()[procPatches[pP]]
-            );
-// 
-//             Pout << nl << "name: " << patch.name() << nl
-//                     << "start: " << patch.start() << nl
-//                     << "size: " << patch.size() << nl
-//                     << "separated: " << Switch(patch.separated()) << nl
-//                     << "parallel: " << Switch(patch.parallel()) << nl << endl;
-// 
-//             forAll (patch, pI)
-//             {
-//                 label decomposedMeshFace = patch.start() + pI;
-// 
-//                 label faceProcAdd = faceProcAddressing[decomposedMeshFace];
-// 
-//                 label globalFace = abs(faceProcAdd)-1;
-// 
-//                 Pout << "Patch index: " << pI
-//                         << " " << patch[pI]
-//                         << " Mesh index: " << decomposedMeshFace
-//                         << " faceProcAdd: " << faceProcAdd
-//                         << " globalFace:" << globalFace;
-// 
-//                 label minStart = -1;
-// 
-//                 // Scanning the dictionary each time is a very ugly way of
-//                 // finding out what patch a face originally belonged to, but
-//                 // it proves the concept.  Read the patch info a container
-//                 // class and have a neat way of tell which patch a face is from
-//                 // embedded in that.  Split each processor face down into
-//                 // separate lists for each different originiating patch.
-// 
-//                 forAll(patchNames, patchi)
-//                 {
-//                     if (patchDictionary.found(patchNames[patchi]))
-//                     {
-//                         const dictionary& patchDict = 
-//                                 patchDictionary.subDict(patchNames[patchi]);
-// 
-//                         word faceName(patchNames[patchi]);
-//                         label startFace(readLabel(patchDict.lookup("startFace")));
-//                         label nFaces(readLabel(patchDict.lookup("nFaces")));
-// 
-//                         if
-//                         (
-//                             minStart < 0
-//                          || startFace < minStart
-//                         )
-//                         {
-//                             minStart = startFace;
-//                         }
-// 
-//                         if
-//                         (
-//                             globalFace >= startFace
-//                          && globalFace < startFace + nFaces
-//                         )
-//                         {
-//                             Pout << " original patch: " << faceName  << endl;
-//                         }
-//                     }
-//                 }
-// 
-//                 if (globalFace < minStart)
-//                 {
-//                     Pout << " originally an internal face" << endl;
-//                 }
-//             }
-// 
-            if (patch.separated())
-            {
-                 Pout << patch.separation();
-            }
-
-            if (!patch.parallel())
-            {
-                 Pout << patch.forwardT();
-            }
-        }
-//     }
-//     else
-//     {
-// 
-//     }
-
-        // Get coords of my shared points
-//     vector sharedPoints(vector::one*(Pstream::myProcNo()+1));
-//     label testRedLab(Pstream::myProcNo()+1);
-    
-//     Pout << testRedLab << endl;
-
-        // Append from all processors
-//     combineReduce(sharedPoints, plusEqOp<vector>());
-//     reduce(testRedLab, plusOp<label>());
-
-//     Pout << testRedLab << endl;
-}

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

@@ -24,131 +24,587 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-namespace Foam
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+inline void Foam::moleculeCloud::evaluatePair
+(
+    molecule* molI,
+    molecule* molJ
+)
 {
+    const pairPotentialList& pairPot(pot_.pairPotentials());
+
+    const electrostaticPotential& electrostatic(pot_.electrostatic());
+
+    label idI = molI->id();
+
+    label idJ = molJ->id();
+
+    const molecule::constantProperties& constPropI(constProps(idI));
+
+    const molecule::constantProperties& constPropJ(constProps(idJ));
+
+    List<label> siteIdsI = constPropI.siteIds();
+
+    List<label> siteIdsJ = constPropJ.siteIds();
+
+    List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
+
+    List<bool> electrostaticSitesI = constPropI.electrostaticSites();
+
+    List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
+
+    List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
+
+    forAll(siteIdsI, sI)
+    {
+        label idsI(siteIdsI[sI]);
+
+        forAll(siteIdsJ, sJ)
+        {
+            label idsJ(siteIdsJ[sJ]);
+
+            if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
+            {
+                vector rsIsJ =
+                molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
+
+                scalar rsIsJMagSq = magSqr(rsIsJ);
+
+                if(pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
+                {
+                    scalar rsIsJMag = mag(rsIsJ);
+
+                    vector fsIsJ = (rsIsJ/rsIsJMag)
+                    *pairPot.force(idsI, idsJ, rsIsJMag);
+
+                    molI->siteForces()[sI] += fsIsJ;
+
+                    molJ->siteForces()[sJ] += -fsIsJ;
+
+                    scalar potentialEnergy
+                    (
+                        pairPot.energy(idsI, idsJ, rsIsJMag)
+                    );
+
+                    molI->potentialEnergy() += 0.5*potentialEnergy;
+
+                    molJ->potentialEnergy() += 0.5*potentialEnergy;
+
+                    molI->rf() += rsIsJ * fsIsJ;
+
+                    molJ->rf() += rsIsJ * fsIsJ;
+                }
+            }
+
+            if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
+            {
+                vector rsIsJ =
+                molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
+
+                scalar rsIsJMagSq = magSqr(rsIsJ);
+
+                if(pairPot.rCutMaxSqr(rsIsJMagSq))
+                {
+                    scalar rsIsJMag = mag(rsIsJ);
+
+                    scalar chargeI = constPropI.siteCharges()[sI];
+
+                    scalar chargeJ = constPropJ.siteCharges()[sJ];
+
+                    vector fsIsJ = (rsIsJ/rsIsJMag)
+                    *chargeI*chargeJ*electrostatic.force(rsIsJMag);
+
+                    molI->siteForces()[sI] += fsIsJ;
+
+                    molJ->siteForces()[sJ] += -fsIsJ;
+
+                    scalar potentialEnergy = chargeI*chargeJ
+                    *electrostatic.energy(rsIsJMag);
+
+                    molI->potentialEnergy() += 0.5*potentialEnergy;
+
+                    molJ->potentialEnergy() += 0.5*potentialEnergy;
+
+                    molI->rf() += rsIsJ * fsIsJ;
+
+                    molJ->rf() += rsIsJ * fsIsJ;
+                }
+            }
+        }
+    }
+}
+
+inline void Foam::moleculeCloud::evaluatePair
+(
+    molecule* molReal,
+    referredMolecule* molRef
+)
+{
+    const pairPotentialList& pairPot(pot_.pairPotentials());
+
+    const electrostaticPotential& electrostatic(pot_.electrostatic());
+
+    label idReal = molReal->id();
+
+    label idRef = molRef->id();
+
+    const molecule::constantProperties& constPropReal(constProps(idReal));
+
+    const molecule::constantProperties& constPropRef(constProps(idRef));
+
+    List<label> siteIdsReal = constPropReal.siteIds();
+
+    List<label> siteIdsRef = constPropRef.siteIds();
+
+    List<bool> pairPotentialSitesReal = constPropReal.pairPotentialSites();
+
+    List<bool> electrostaticSitesReal = constPropReal.electrostaticSites();
+
+    List<bool> pairPotentialSitesRef = constPropRef.pairPotentialSites();
+
+    List<bool> electrostaticSitesRef = constPropRef.electrostaticSites();
+
+    forAll(siteIdsReal, sReal)
+    {
+        label idsReal(siteIdsReal[sReal]);
+
+        forAll(siteIdsRef, sRef)
+        {
+            label idsRef(siteIdsRef[sRef]);
+
+            if (pairPotentialSitesReal[sReal] && pairPotentialSitesRef[sRef])
+            {
+                vector rsRealsRef =
+                molReal->sitePositions()[sReal] - molRef->sitePositions()[sRef];
+
+                scalar rsRealsRefMagSq = magSqr(rsRealsRef);
+
+                if(pairPot.rCutSqr(idsReal, idsRef, rsRealsRefMagSq))
+                {
+                    scalar rsRealsRefMag = mag(rsRealsRef);
+
+                    vector fsRealsRef = (rsRealsRef/rsRealsRefMag)
+                    *pairPot.force(idsReal, idsRef, rsRealsRefMag);
+
+                    molReal->siteForces()[sReal] += fsRealsRef;
+
+                    scalar potentialEnergy
+                    (
+                        pairPot.energy(idsReal, idsRef, rsRealsRefMag)
+                    );
+
+                    molReal->potentialEnergy() += 0.5*potentialEnergy;
+
+                    molReal->rf() += rsRealsRef * fsRealsRef;
+
+                }
+            }
+
+            if (electrostaticSitesReal[sReal] && electrostaticSitesRef[sRef])
+            {
+                vector rsRealsRef =
+                molReal->sitePositions()[sReal] - molRef->sitePositions()[sRef];
+
+                scalar rsRealsRefMagSq = magSqr(rsRealsRef);
+
+                if(pairPot.rCutMaxSqr(rsRealsRefMagSq))
+                {
+                    scalar rsRealsRefMag = mag(rsRealsRef);
+
+                    scalar chargeReal = constPropReal.siteCharges()[sReal];
+
+                    scalar chargeRef = constPropRef.siteCharges()[sRef];
+
+                    vector fsRealsRef = (rsRealsRef/rsRealsRefMag)
+                    *chargeReal*chargeRef
+                    *electrostatic.force(rsRealsRefMag);
+
+                    molReal->siteForces()[sReal] += fsRealsRef;
+
+                    scalar potentialEnergy = chargeReal*chargeRef
+                    *electrostatic.energy(rsRealsRefMag);
+
+                    molReal->potentialEnergy() += 0.5*potentialEnergy;
+
+                    molReal->rf() += rsRealsRef * fsRealsRef;
+                }
+            }
+        }
+    }
+}
+
+
+inline bool Foam::moleculeCloud::evaluatePotentialLimit
+(
+    molecule* molI,
+    molecule* molJ
+) const
+{
+    const pairPotentialList& pairPot(pot_.pairPotentials());
+
+    const electrostaticPotential& electrostatic(pot_.electrostatic());
+
+    label idI = molI->id();
+
+    label idJ = molJ->id();
+
+    const molecule::constantProperties& constPropI(constProps(idI));
+
+    const molecule::constantProperties& constPropJ(constProps(idJ));
+
+    List<label> siteIdsI = constPropI.siteIds();
+
+    List<label> siteIdsJ = constPropJ.siteIds();
+
+    List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
+
+    List<bool> electrostaticSitesI = constPropI.electrostaticSites();
+
+    List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
+
+    List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
+
+    forAll(siteIdsI, sI)
+    {
+        label idsI(siteIdsI[sI]);
+
+        forAll(siteIdsJ, sJ)
+        {
+            label idsJ(siteIdsJ[sJ]);
+
+            if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
+            {
+                vector rsIsJ =
+                molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
+
+                scalar rsIsJMagSq = magSqr(rsIsJ);
+
+                if(pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
+                {
+                    scalar rsIsJMag = mag(rsIsJ);
+
+                    // Guard against pairPot.energy being evaluated
+                    // if rIJMag < SMALL. A floating point exception will
+                    // happen otherwise.
+
+                    if (rsIsJMag < SMALL)
+                    {
+                        WarningIn("moleculeCloud::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 molecules"
+                            << " twice. Removing one of the molecules."
+                            << endl;
+
+                        return true;
+                    }
+
+                    // Guard against pairPot.energy being evaluated if rIJMag <
+                    // rMin.  A tabulation lookup error will occur otherwise.
+
+                    if (rsIsJMag < pairPot.rMin(idsI, idsJ))
+                    {
+                        return true;
+                    }
+
+                    if
+                    (
+                        mag(pairPot.energy(idsI, idsJ, rsIsJMag))
+                      > pot_.potentialEnergyLimit()
+                    )
+                    {
+                        return true;
+                    };
+                }
+            }
+
+            if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
+            {
+                vector rsIsJ =
+                molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
+
+                scalar rsIsJMagSq = magSqr(rsIsJ);
+
+                if(pairPot.rCutMaxSqr(rsIsJMagSq))
+                {
+                    scalar rsIsJMag = mag(rsIsJ);
+
+                    // Guard against pairPot.energy being evaluated
+                    // if rIJMag < SMALL. A floating point exception will
+                    // happen otherwise.
+
+                    if (rsIsJMag < SMALL)
+                    {
+                        WarningIn("moleculeCloud::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 molConfig in"
+                            << " parallel or a block filled with molecules "
+                            << " twice. Removing one of the molecules."
+                            << endl;
+
+                        return true;
+                    }
+
+                    scalar chargeI = constPropI.siteCharges()[sI];
+
+                    scalar chargeJ = constPropJ.siteCharges()[sJ];
+
+                    if
+                    (
+                        mag(chargeI*chargeJ*electrostatic.energy(rsIsJMag))
+                      > pot_.potentialEnergyLimit()
+                    )
+                    {
+                        return true;
+                    };
+                }
+            }
+        }
+    }
+
+    return false;
+}
+
+
+inline bool Foam::moleculeCloud::evaluatePotentialLimit
+(
+    molecule* molReal,
+    referredMolecule* molRef
+) const
+{
+    const pairPotentialList& pairPot(pot_.pairPotentials());
+
+    const electrostaticPotential& electrostatic(pot_.electrostatic());
+
+    label idReal = molReal->id();
+
+    label idRef = molRef->id();
+
+    const molecule::constantProperties& constPropReal(constProps(idReal));
+
+    const molecule::constantProperties& constPropRef(constProps(idRef));
+
+    List<label> siteIdsReal = constPropReal.siteIds();
+
+    List<label> siteIdsRef = constPropRef.siteIds();
+
+    List<bool> pairPotentialSitesReal = constPropReal.pairPotentialSites();
+
+    List<bool> electrostaticSitesReal = constPropReal.electrostaticSites();
+
+    List<bool> pairPotentialSitesRef = constPropRef.pairPotentialSites();
+
+    List<bool> electrostaticSitesRef = constPropRef.electrostaticSites();
+
+    forAll(siteIdsReal, sReal)
+    {
+        label idsReal(siteIdsReal[sReal]);
+
+        forAll(siteIdsRef, sRef)
+        {
+            label idsRef(siteIdsRef[sRef]);
+
+            if (pairPotentialSitesReal[sReal] && pairPotentialSitesRef[sRef])
+            {
+                vector rsRealsRef =
+                molReal->sitePositions()[sReal] - molRef->sitePositions()[sRef];
+
+                scalar rsRealsRefMagSq = magSqr(rsRealsRef);
+
+                if(pairPot.rCutSqr(idsReal, idsRef, rsRealsRefMagSq))
+                {
+                    scalar rsRealsRefMag = mag(rsRealsRef);
+
+                    // Guard against pairPot.energy being evaluated
+                    // if rRealRefMag < SMALL. A floating point exception will
+                    // happen otherwise.
+
+                    if (rsRealsRefMag < SMALL)
+                    {
+                        WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
+                        << "Molecule site pair closer than "
+                            << SMALL
+                            << ": mag separation = " << rsRealsRefMag
+                            << ". These may have been placed on top of each"
+                            << " other by a rounding error in molConfig in"
+                            << " parallel or a block filled with molecules "
+                            << " twice. Removing one of the molecules."
+                            << endl;
+
+                        return true;
+                    }
+
+                    // Guard against pairPot.energy being evaluated if
+                    // rRealRefMag < rMin.  A tabulation lookup error will occur
+                    // otherwise.
+
+                    if (rsRealsRefMag < pairPot.rMin(idsReal, idsRef))
+                    {
+                        return true;
+                    }
+
+                    if
+                    (
+                        mag(pairPot.energy(idsReal, idsRef, rsRealsRefMag))
+                      > pot_.potentialEnergyLimit()
+                    )
+                    {
+                        return true;
+                    };
+                }
+            }
+
+            if (electrostaticSitesReal[sReal] && electrostaticSitesRef[sRef])
+            {
+                vector rsRealsRef =
+                molReal->sitePositions()[sReal] - molRef->sitePositions()[sRef];
+
+                scalar rsRealsRefMagSq = magSqr(rsRealsRef);
+
+                if(pairPot.rCutMaxSqr(rsRealsRefMagSq))
+                {
+                    scalar rsRealsRefMag = mag(rsRealsRef);
+
+                    // Guard against pairPot.energy being evaluated
+                    // if rRealRefMag < SMALL. A floating point exception will
+                    // happen otherwise.
+
+                    if (rsRealsRefMag < SMALL)
+                    {
+                        WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
+                        << "Molecule site pair closer than "
+                            << SMALL
+                            << ": mag separation = " << rsRealsRefMag
+                            << ". These may have been placed on top of each"
+                            << " other by a rounding error in molConfig in"
+                            << " parallel or a block filled with molecules "
+                            << " twice. Removing one of the molecules."
+                            << endl;
+
+                        return true;
+                    }
+
+                    scalar chargeReal = constPropReal.siteCharges()[sReal];
+
+                    scalar chargeRef = constPropRef.siteCharges()[sRef];
+
+                    if
+                    (
+                        mag
+                        (
+                            chargeReal*chargeRef
+                            *electrostatic.energy(rsRealsRefMag)
+                        )
+                      > pot_.potentialEnergyLimit()
+                    )
+                    {
+                        return true;
+                    };
+                }
+            }
+        }
+    }
+
+    return false;
+}
+
+
+inline Foam::vector Foam::moleculeCloud::equipartitionLinearVelocity
+(
+    scalar temperature,
+    scalar mass
+)
+{
+    return sqrt(kb*temperature/mass)*vector
+    (
+        rndGen_.GaussNormal(),
+        rndGen_.GaussNormal(),
+        rndGen_.GaussNormal()
+    );
+}
+
+
+inline Foam::vector Foam::moleculeCloud::equipartitionAngularMomentum
+(
+    scalar temperature,
+    const molecule::constantProperties& cP
+)
+{
+    scalar sqrtKbT = sqrt(kb*temperature);
+
+    if (cP.linearMolecule())
+    {
+        return sqrtKbT*vector
+        (
+            0.0,
+            sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal(),
+            sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal()
+        );
+    }
+    else
+    {
+        return sqrtKbT*vector
+        (
+            sqrt(cP.momentOfInertia().xx())*rndGen_.GaussNormal(),
+            sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal(),
+            sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal()
+        );
+    }
+}
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline const polyMesh& moleculeCloud::mesh() const
+inline const Foam::polyMesh& Foam::moleculeCloud::mesh() const
 {
     return mesh_;
 }
 
 
-// MD solution parameters
+inline const Foam::potential& Foam::moleculeCloud::pot() const
-
-inline const moleculeCloud::integrationMethods&
-    moleculeCloud::integrationMethod() const
 {
-    return integrationMethod_;
+    return pot_;
 }
 
 
-inline scalar moleculeCloud::potentialEnergyLimit() const
+inline const Foam::List<Foam::DynamicList<Foam::molecule*> >&
-{
+    Foam::moleculeCloud::cellOccupancy() const
-    return potentialEnergyLimit_;
-}
-
-
-inline label moleculeCloud::nPairPotentials() const
-{
-    return pairPotentials_.size();
-}
-
-
-inline const labelList& moleculeCloud::removalOrder() const
-{
-    return removalOrder_;
-}
-
-
-inline const labelListList& moleculeCloud::directInteractionList() const
-{
-    return directInteractionList_;
-}
-
-
-inline const referredCellList& moleculeCloud::referredInteractionList() const
-{
-    return referredInteractionList_;
-}
-
-
-inline const labelList&
-    moleculeCloud::realCellsWithinRCutMaxOfAnyReferringPatch() const
-{
-    return realCellsWithinRCutMaxOfAnyReferringPatch_;
-}
-
-
-inline const labelList&
-    moleculeCloud::realFacesWithinRCutMaxOfAnyReferringPatch() const
-{
-    return realFacesWithinRCutMaxOfAnyReferringPatch_;
-}
-
-
-inline const labelList&
-    moleculeCloud::realEdgesWithinRCutMaxOfAnyReferringPatch() const
-{
-    return realEdgesWithinRCutMaxOfAnyReferringPatch_;
-}
-
-
-inline const labelList&
-    moleculeCloud::realPointsWithinRCutMaxOfAnyReferringPatch() const
-{
-    return realPointsWithinRCutMaxOfAnyReferringPatch_;
-}
-
-
-inline const List<sendingReferralList>&
-    moleculeCloud::cellSendingReferralLists() const
-{
-    return cellSendingReferralLists_;
-}
-
-
-inline const List<receivingReferralList>&
-    moleculeCloud::cellReceivingReferralLists() const
-{
-    return cellReceivingReferralLists_;
-}
-
-
-inline label moleculeCloud::nInteractingProcs() const
-{
-    return cellReceivingReferralLists_.size();
-}
-
-
-inline const pairPotentialList& moleculeCloud::pairPotentials() const
-{
-    return pairPotentials_;
-}
-
-
-inline const tetherPotentialList& moleculeCloud::tetherPotentials() const
-{
-    return tetherPotentials_;
-}
-
-
-inline const vector& moleculeCloud::gravity() const
-{
-    return gravity_;
-}
-
-
-inline const List< DynamicList<molecule*> >&
-    moleculeCloud::cellOccupancy() const
 {
     return cellOccupancy_;
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+inline const Foam::interactionLists&
+    Foam::moleculeCloud::il() const
+{
+    return il_;
+}
+
+
+inline const Foam::List<Foam::molecule::constantProperties>
+    Foam::moleculeCloud::constProps() const
+{
+    return constPropList_;
+}
+
+
+inline const Foam::molecule::constantProperties&
+    Foam::moleculeCloud::constProps(label id) const
+{
+    return constPropList_[id];
+}
+
+
+inline Foam::Random& Foam::moleculeCloud::rndGen()
+{
+    return rndGen_;
+}
 
-}  // End namespace Foam
 
 // ************************************************************************* //

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

@@ -1,50 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 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 "moleculeCloud.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-void Foam::moleculeCloud::integrateEquationsOfMotion()
-{
-    // Supply trackData to the move method to tell it whether this is the
-    // leapfrog stage 1 or stage 2
-    // or whether it is the
-    // predictor or corrector step.
-
-    molecule::trackData td1(*this, 1);
-
-    Cloud<molecule>::move(td1);
-
-    calculateForce();
-
-    molecule::trackData td2(*this, 2);
-
-    Cloud<molecule>::move(td2);
-}
-
-
-// ************************************************************************* //

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

@@ -1,177 +0,0 @@
-Info<< nl <<  "Reading MD solution parameters:" << endl;
-
-IOdictionary mdSolution
-(
-    IOobject
-    (
-        "mdSolution",
-        mesh_.time().system(),
-        mesh_,
-        IOobject::MUST_READ,
-        IOobject::NO_WRITE
-    )
-);
-
-integrationMethod_ = integrationMethodNames_.read
-(
-    mdSolution.lookup("integrationMethod")
-);
-
-potentialEnergyLimit_ = readScalar
-(
-    mdSolution.lookup("potentialEnergyLimit")
-);
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-IOdictionary potentialDict
-(
-    IOobject
-    (
-        "potentialDict",
-        mesh_.time().system(),
-        mesh_,
-        IOobject::MUST_READ,
-        IOobject::NO_WRITE
-    )
-);
-
-IOdictionary idListDict
-(
-    IOobject
-    (
-        "idList",
-        mesh_.time().constant(),
-        mesh_,
-        IOobject::MUST_READ,
-        IOobject::NO_WRITE
-    )
-);
-
-// ****************************************************************************
-// Pair potentials
-
-if (!potentialDict.found("pair"))
-{
-    FatalErrorIn("moleculeCloudReadPotentials.H") << nl
-            << "pair potential specification subDict not found"
-            << abort(FatalError);
-}
-
-const dictionary& pairDict = potentialDict.subDict("pair");
-
-pairPotentials_.buildPotentials(idListDict, pairDict, mesh_);
-
-if (potentialDict.found("removalOrder"))
-{
-    List<word> remOrd = potentialDict.lookup("removalOrder");
-
-    removalOrder_.setSize(remOrd.size());
-
-    forAll(removalOrder_, rO)
-    {
-        removalOrder_[rO] = findIndex(pairPotentials_.idList(), remOrd[rO]);
-    }
-}
-
-// ****************************************************************************
-// Tether potentials
-
-iterator mol(this->begin());
-
-DynamicList<label> tetherIds;
-
-for
-(
- mol = this->begin();
- mol != this->end();
- ++mol
-)
-{
-    if (mol().tethered())
-    {
-        if (findIndex(tetherIds, mol().id()) == -1)
-        {
-            tetherIds.append
-            (
-                mol().id()
-            );
-        }
-    }
-}
-
-if (Pstream::parRun())
-{
-    List< labelList > allTetherIds(Pstream::nProcs());
-
-    allTetherIds[Pstream::myProcNo()] = tetherIds;
-
-    Pstream::gatherList(allTetherIds);
-
-    if (Pstream::master())
-    {
-        DynamicList<label> globalTetherIds;
-
-        forAll(allTetherIds, procN)
-        {
-            const labelList& procNTetherIds = allTetherIds[procN];
-
-            forAll(procNTetherIds, id)
-            {
-                if (findIndex(globalTetherIds, procNTetherIds[id]) == -1)
-                {
-                    globalTetherIds.append
-                    (
-                        procNTetherIds[id]
-                    );
-                }
-            }
-        }
-
-        globalTetherIds.shrink();
-
-        tetherIds = globalTetherIds;
-    }
-
-    Pstream::scatter(tetherIds);
-}
-
-tetherIds.shrink();
-
-if (tetherIds.size())
-{
-    if (!potentialDict.found("tether"))
-    {
-        FatalErrorIn("moleculeCloudReadPotentials.H") << nl
-                << "tether potential specification subDict not found"
-                << abort(FatalError);
-    }
-
-    const dictionary& tetherDict = potentialDict.subDict("tether");
-
-    tetherPotentials_.buildPotentials(idListDict, tetherDict, tetherIds);
-}
-
-// ****************************************************************************
-// External Forces
-
-gravity_ = vector::zero;
-
-if (potentialDict.found("external"))
-{
-
-    Info << nl << "Reading external forces:" << endl;
-
-    const dictionary& externalDict = potentialDict.subDict("external");
-
-    // ************************************************************************
-    // gravity
-
-    if (externalDict.found("gravity"))
-    {
-        gravity_ = externalDict.lookup("gravity");
-    }
-}
-
-Info << nl << tab << "gravity = " << gravity_ << endl;
-

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

@@ -1,152 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 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 "moleculeCloud.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-const Foam::labelList Foam::moleculeCloud::realCellsInRangeOfSegment
-(
-    const labelList& segmentFaces,
-    const labelList& segmentEdges,
-    const labelList& segmentPoints
-) const
-{
-    DynamicList<label> realCellsFoundInRange;
-
-    forAll(segmentFaces, sF)
-    {
-        const label f = segmentFaces[sF];
-
-        forAll (mesh_.points(), p)
-        {
-            if (testPointFaceDistance(p, f))
-            {
-                const labelList& pCells(mesh_.pointCells()[p]);
-
-                forAll(pCells, pC)
-                {
-                    const label cellI(pCells[pC]);
-
-                    if (findIndex(realCellsFoundInRange, cellI) == -1)
-                    {
-                        realCellsFoundInRange.append(cellI);
-                    }
-                }
-            }
-        }
-    }
-
-    forAll(segmentPoints, sP)
-    {
-        const label p = segmentPoints[sP];
-
-        forAll(mesh_.faces(), f)
-        {
-            if (testPointFaceDistance(p, f))
-            {
-                const label cellO(mesh_.faceOwner()[f]);
-
-                if (findIndex(realCellsFoundInRange, cellO) == -1)
-                {
-                    realCellsFoundInRange.append(cellO);
-                }
-
-                if (mesh_.isInternalFace(f))
-                {
-                    // boundary faces will not have neighbour information
-
-                    const label cellN(mesh_.faceNeighbour()[f]);
-
-                    if (findIndex(realCellsFoundInRange, cellN) == -1)
-                    {
-                        realCellsFoundInRange.append(cellN);
-                    }
-                }
-            }
-        }
-    }
-
-    forAll(segmentEdges, sE)
-    {
-        const edge& eJ(mesh_.edges()[segmentEdges[sE]]);
-
-        forAll (mesh_.edges(), edgeIIndex)
-        {
-            const edge& eI(mesh_.edges()[edgeIIndex]);
-
-            if (testEdgeEdgeDistance(eI, eJ))
-            {
-                const labelList& eICells(mesh_.edgeCells()[edgeIIndex]);
-
-                forAll(eICells, eIC)
-                {
-                    const label cellI(eICells[eIC]);
-
-                    if (findIndex(realCellsFoundInRange, cellI) == -1)
-                    {
-                        realCellsFoundInRange.append(cellI);
-                    }
-                }
-            }
-        }
-    }
-
-//     forAll (points, p)
-//     {
-//         const point& ptI = mesh_.points()[points[p]];
-//
-//         forAll(mesh_.faces(), f)
-//         {
-//             if (testPointFaceDistance(ptI, f))
-//             {
-//                 const label cellO(mesh_.faceOwner()[f]);
-//
-//                 if (findIndex(realCellsFoundInRange, cellO) == -1)
-//                 {
-//                     realCellsFoundInRange.append(cellO);
-//                 }
-//
-//                 if (mesh_.isInternalFace(f))
-//                 {
-//                     // boundary faces will not have neighbour information
-//
-//                     const label cellN(mesh_.faceNeighbour()[f]);
-//
-//                     if (findIndex(realCellsFoundInRange, cellN) == -1)
-//                     {
-//                         realCellsFoundInRange.append(cellN);
-//                     }
-//                 }
-//             }
-//         }
-//     }
-//
-    return realCellsFoundInRange.shrink();
-}
-
-
-// ************************************************************************* //

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

@@ -1,117 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 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 "moleculeCloud.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-const Foam::labelList Foam::moleculeCloud::referredCellsInRangeOfSegment
-(
-    const List<referredCell>& referredInteractionList,
-    const labelList& segmentFaces,
-    const labelList& segmentEdges,
-    const labelList& segmentPoints
-) const
-{
-    DynamicList<label> referredCellsFoundInRange;
-
-    forAll(segmentFaces, sF)
-    {
-        const label f = segmentFaces[sF];
-
-        forAll(referredInteractionList, rIL)
-        {
-            const vectorList& refCellPoints
-                = referredInteractionList[rIL].vertexPositions();
-
-            if (testPointFaceDistance(refCellPoints, f))
-            {
-                if (findIndex(referredCellsFoundInRange, rIL) == -1)
-                {
-                    referredCellsFoundInRange.append(rIL);
-                }
-            }
-        }
-    }
-
-    forAll(segmentPoints, sP)
-    {
-        const label p = segmentPoints[sP];
-
-        forAll(referredInteractionList, rIL)
-        {
-            const referredCell& refCell(referredInteractionList[rIL]);
-
-            if (testPointFaceDistance(p, refCell))
-            {
-                if (findIndex(referredCellsFoundInRange, rIL) == -1)
-                {
-                    referredCellsFoundInRange.append(rIL);
-                }
-            }
-        }
-    }
-
-    forAll(segmentEdges, sE)
-    {
-        const edge& eI(mesh_.edges()[segmentEdges[sE]]);
-
-        forAll(referredInteractionList, rIL)
-        {
-            const vectorList& refCellPoints
-                = referredInteractionList[rIL].vertexPositions();
-
-            const edgeList& refCellEdges
-                = referredInteractionList[rIL].edges();
-
-            forAll(refCellEdges, rCE)
-            {
-                const edge& eJ(refCellEdges[rCE]);
-
-                if
-                (
-                    testEdgeEdgeDistance
-                    (
-                        eI,
-                        refCellPoints[eJ.start()],
-                        refCellPoints[eJ.end()]
-                    )
-                )
-                {
-                    if(findIndex(referredCellsFoundInRange, rIL) == -1)
-                    {
-                        referredCellsFoundInRange.append(rIL);
-                    }
-                }
-            }
-        }
-    }
-
-    return referredCellsFoundInRange.shrink();
-}
-
-
-// ************************************************************************* //

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

@@ -1,55 +0,0 @@
-{ 
-    vector rIJ;
-
-    scalar rIJMag;
-
-    scalar rIJMagSq;
-
-    label idI;
-
-    label idJ;
-
-    mol = this->begin();
-
-    molecule* molI = &mol();
-
-    molecule* molJ = &mol();
-
-    DynamicList<molecule*> molsToDelete;
-
-    forAll(directInteractionList_, dIL)
-    {
-        forAll(cellOccupancy_[dIL],cellIMols)
-        {
-            molI = cellOccupancy_[dIL][cellIMols];
-
-            forAll(directInteractionList_[dIL], interactingCells)
-            {
-                List< molecule* > cellJ =
-                    cellOccupancy_[directInteractionList_[dIL][interactingCells]];
-
-                forAll(cellJ, cellJMols)
-                {
-                    molJ = cellJ[cellJMols];
-
-#                   include "moleculeCloudRemoveHighEnergyOverlapsRealCellsCalculationStep.H"
-                }
-            }
-
-            forAll(cellOccupancy_[dIL],cellIOtherMols)
-            {
-                molJ = cellOccupancy_[dIL][cellIOtherMols];
-
-                if (molJ > molI)
-                {
-#                   include "moleculeCloudRemoveHighEnergyOverlapsRealCellsCalculationStep.H"
-                }
-            }
-        }
-    }
-
-    forAll (molsToDelete, mTD)
-    {
-        deleteParticle(*(molsToDelete[mTD]));
-    }
-}

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

@@ -1,79 +0,0 @@
-idI = molI->id();
-
-idJ = molJ->id();
-
-rIJ = molI->position() - molJ->position();
-
-rIJMagSq = magSqr(rIJ);
-
-if (pairPotentials_.rCutSqr(idI, idJ, rIJMagSq))
-{
-    rIJMag = mag(rIJ);
-
-    bool remove = false;
-
-    // Guard against pairPotentials_.energy being evaluated
-    // if rIJMag < SMALL. A floating point exception will
-    // happen otherwise.
-
-    if (rIJMag < SMALL)
-    {
-        WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
-            << "Real molecule pair "
-            << " idI = " << idI
-            << " at position " << molI->position()
-            << " idJ = " << idJ
-            << " at position " << molJ->position()
-            << " are closer than " << SMALL
-            << ": mag separation = " << rIJMag
-            << ". These may have been placed on top of each"
-            << " other by a rounding error in molConfig in parallel"
-            << " or a block filled with molecules twice."
-            << " Removing one of the molecules."
-            << endl;
-
-        remove = true;
-    }
-    
-    // Guard against pairPotentials_.energy being evaluated
-    // if rIJMag < rMin. A tabulation lookup error will occur otherwise.
-
-    if (rIJMag < pairPotentials_.rMin(idI, idJ))
-    {
-        remove = true;
-    }
-
-    if (!remove)
-    {
-        if
-        (
-            pairPotentials_.energy(idI, idJ, rIJMag) > potentialEnergyLimit_
-        )
-        {
-           
-            remove = true;
-        }
-    }
-
-    if (remove)
-    {
-        if
-        (
-            idI == idJ
-         || findIndex(removalOrder_, idJ) < findIndex(removalOrder_, idI)
-        )
-        {
-            if (findIndex(molsToDelete, molJ) == -1)
-            {
-                molsToDelete.append(molJ);
-            }
-        }
-        else
-        {
-            if (findIndex(molsToDelete, molI) == -1)
-            {
-                molsToDelete.append(molI);
-            }
-        }
-    }
-}

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

@@ -1,152 +0,0 @@
-{
-    vector rKL;
-
-    scalar rKLMag;
-
-    scalar rKLMagSq;
-
-    label idK;
-
-    label idL;
-
-    molecule* molK = &mol();
-
-    DynamicList<molecule*> molsToDelete;
-
-    forAll(referredInteractionList_, rIL)
-    {
-        referredCell& refCell = referredInteractionList_[rIL];
-
-        forAll(refCell, refMols)
-        {
-            referredMolecule* molL = &(refCell[refMols]);
-
-            List <label> realCells = refCell.realCellsForInteraction();
-
-            forAll(realCells, rC)
-            {
-                label cellK = realCells[rC];
-
-                List<molecule*> cellKMols = cellOccupancy_[cellK];
-
-                forAll(cellKMols, cKM)
-                {
-                    molK = cellKMols[cKM];
-
-                    idK = molK->id();
-
-                    idL = molL->id();
-
-                    rKL = molK->position() - molL->position();
-
-                    rKLMagSq = magSqr(rKL);
-
-                    if (pairPotentials_.rCutSqr(idK, idL, rKLMagSq))
-                    {
-                        rKLMag = mag(rKL);
-
-                        bool remove = false;
-
-                        // Guard against pairPotentials_.energy being evaluated
-                        // if rKLMag < SMALL. A floating point exception will
-                        // happen otherwise.
-
-                        if (rKLMag < SMALL)
-                        {
-                            WarningIn
-                            (
-                                "moleculeCloud::removeHighEnergyOverlaps()"
-                            )
-                                << "Real-referred molecule pair "
-                                << " idK = " << idK << " (real)"
-                                << " at position " << molK->position()
-                                << " idL = " << idL << " (referred)"
-                                << " at position " << molL->position()
-                                << " are closer than " << SMALL
-                                << ": mag separation = " << rKLMag
-                                << ". These may have been placed on top of each"
-                                << " other by a rounding error in molConfig in"
-                                << " parallel, a block filled with molecules"
-                                << " twice, or a lattice ending very close"
-                                << " to the edge of the mesh."
-                                << " Removing one of the molecules."
-                                << endl;
-
-                                remove = true;
-                        }
-
-                        // Guard against pairPotentials_.energy being evaluated
-                        // if rIJMag < rMin. A tubulation lookup error will occur otherwise.
-
-                        if (rKLMag < pairPotentials_.rMin(idK, idL))
-                        {
-                            remove = true;
-                        }
-
-                        if (!remove)
-                        {
-                            if
-                            (
-                                pairPotentials_.energy(idK, idL, rKLMag)
-                              > potentialEnergyLimit_
-                            )
-                            {
-                                remove = true;
-                            }
-                        }
-
-                        if (remove)
-                        {
-                            if
-                            (
-                                findIndex(removalOrder_, idK)
-                              < findIndex(removalOrder_, idL)
-                            )
-                            {
-                                if (findIndex(molsToDelete, molK) == -1)
-                                {
-                                    molsToDelete.append(molK);
-                                }
-                            }
-
-                            else if
-                            (
-                                findIndex(removalOrder_, idK)
-                             == findIndex(removalOrder_, idL)
-                            )
-                            {
-                                if
-                                (
-                                    Pstream::myProcNo() == refCell.sourceProc()
-                                 && cellK <= refCell.sourceCell()
-                                )
-                                {
-                                    if (findIndex(molsToDelete, molK) == -1)
-                                    {
-                                        molsToDelete.append(molK);
-                                    }
-                                }
-
-                                else if
-                                (
-                                    Pstream::myProcNo() < refCell.sourceProc()
-                                )
-                                {
-                                    if (findIndex(molsToDelete, molK) == -1)
-                                    {
-                                        molsToDelete.append(molK);
-                                    }
-                                }
-                            }
-                        }
-                    }
-                }
-            }
-        }
-    }
-
-    forAll (molsToDelete, mTD)
-    {
-        deleteParticle(*(molsToDelete[mTD]));
-    }
-}

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

@@ -1,235 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 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 "moleculeCloud.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-bool Foam::moleculeCloud::testPointFaceDistance
-(
-    const label p,
-    const label faceNo
-) const
-{
-    const vector& pointPosition(mesh_.points()[p]);
-
-    return testPointFaceDistance(pointPosition, faceNo);
-}
-
-bool Foam::moleculeCloud::testPointFaceDistance
-(
-    const label p,
-    const referredCell& refCell
-) const
-{
-    const vector& pointPosition(mesh_.points()[p]);
-
-    forAll (refCell.faces(), rCF)
-    {
-        if
-        (
-            testPointFaceDistance
-            (
-                pointPosition,
-                refCell.faces()[rCF],
-                refCell.vertexPositions(),
-                refCell.faceCentres()[rCF],
-                refCell.faceAreas()[rCF]
-            )
-        )
-        {
-            return true;
-        }
-    }
-
-    return false;
-}
-
-bool Foam::moleculeCloud::testPointFaceDistance
-(
-    const vectorList& pointsToTest,
-    const label faceNo
-) const
-{
-    forAll(pointsToTest, pTT)
-    {
-        const vector& p(pointsToTest[pTT]);
-
-        // if any point in the list is in range of the face
-        // then the rest do not need to be tested and
-        // true can be returned
-
-        if (testPointFaceDistance(p, faceNo))
-        {
-            return true;
-        }
-    }
-
-    return false;
-}
-
-bool Foam::moleculeCloud::testPointFaceDistance
-(
-    const vector& p,
-    const label faceNo
-) const
-{
-    const face& faceToTest(mesh_.faces()[faceNo]);
-
-    const vector& faceC(mesh_.faceCentres()[faceNo]);
-
-    const vector& faceA(mesh_.faceAreas()[faceNo]);
-
-    const vectorList& points(mesh_.points());
-
-    return testPointFaceDistance
-    (
-        p,
-        faceToTest,
-        points,
-        faceC,
-        faceA
-    );
-}
-
-bool Foam::moleculeCloud::testPointFaceDistance
-(
-    const vector& p,
-    const labelList& faceToTest,
-    const vectorList& points,
-    const vector& faceC,
-    const vector& faceA
-) const
-{
-    vector faceN(faceA/mag(faceA));
-
-    scalar perpDist((p - faceC) & faceN);
-
-    if (mag(perpDist) > pairPotentials_.rCutMax())
-    {
-        return false;
-    }
-
-    vector pointOnPlane = (p - faceN * perpDist);
-
-    if (magSqr(faceC - pointOnPlane) < pairPotentials_.rCutMaxSqr()*1e-8)
-    {
-        // If pointOnPlane is very close to the face centre
-        // then defining the local axes will be inaccurate
-        // and it is very likely that pointOnPlane will be
-        // inside the face, so return true if the points
-        // are in range to be safe
-
-        return (magSqr(pointOnPlane - p) <= pairPotentials_.rCutMaxSqr());
-    }
-
-    vector xAxis = (faceC - pointOnPlane)/mag(faceC - pointOnPlane);
-
-    vector yAxis =
-        ((faceC - pointOnPlane) ^ faceN)
-       /mag((faceC - pointOnPlane) ^ faceN);
-
-    List<vector2D> local2DVertices(faceToTest.size());
-
-    forAll(faceToTest, fTT)
-    {
-        const vector& V(points[faceToTest[fTT]]);
-
-        if (magSqr(V-p) <= pairPotentials_.rCutMaxSqr())
-        {
-            return true;
-        }
-
-        local2DVertices[fTT] = vector2D
-        (
-            ((V - pointOnPlane) & xAxis),
-            ((V - pointOnPlane) & yAxis)
-        );
-    }
-
-    scalar localFaceCx((faceC - pointOnPlane) & xAxis);
-
-    scalar la_valid = -1;
-
-    forAll(local2DVertices, fV)
-    {
-        const vector2D& va(local2DVertices[fV]);
-
-        const vector2D& vb
-        (
-            local2DVertices[(fV + 1) % local2DVertices.size()]
-        );
-
-        if (mag(vb.y()-va.y()) > SMALL)
-        {
-            scalar la =
-                (
-                    va.x() - va.y()*((vb.x() - va.x())/(vb.y() - va.y()))
-                )
-               /localFaceCx;
-
-            scalar lv = -va.y()/(vb.y() - va.y());
-
-
-            if (la >= 0 && la <= 1 && lv >= 0 && lv <= 1)
-            {
-                la_valid = la;
-
-                break;
-            }
-        }
-    }
-
-    if (la_valid < 0)
-    {
-        // perpendicular point inside face, nearest point is pointOnPlane;
-        return (magSqr(pointOnPlane-p) <= pairPotentials_.rCutMaxSqr());
-    }
-    else
-    {
-        // perpendicular point outside face, nearest point is
-        // on edge that generated la_valid;
-        return
-        (
-            magSqr(pointOnPlane + la_valid*(faceC - pointOnPlane) - p)
-         <= pairPotentials_.rCutMaxSqr()
-        );
-    }
-
-    // if the algorithm hasn't returned anything by now then something has
-    // gone wrong.
-
-    FatalErrorIn("moleculeCloudTestPointAndFaceDistance.C") << nl
-        << "point " << p << " to face " << faceToTest
-        << " comparison did not find a nearest point"
-        << " to be inside or outside face."
-        << abort(FatalError);
-
-    return false;
-}
-
-
-// ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/referredCellList/referredCellList.C

@@ -1,327 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 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 "referredCellList.H"
-#include "moleculeCloud.H"
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-referredCellList::referredCellList(moleculeCloud& molCloud)
-:
-    List<referredCell >(),
-    molCloud_(molCloud)
-{}
-
-
-referredCellList::referredCellList
-(
-    moleculeCloud& molCloud,
-    const List<referredCell>& referredCells,
-    const List<label>& realCells
-)
-:
-    List< referredCell >(referredCells),
-    molCloud_(molCloud)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-referredCellList::~referredCellList()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void referredCellList::setRealCellsInRange()
-{
-    Info<< nl << "Finding real cells in range of referred cells" << endl;
-
-    forAll(*this, rC)
-    {
-        const polyMesh& mesh(molCloud_.mesh());
-
-        referredCell& refCell = (*this)[rC];
-
-        DynamicList<label> realCellsFoundInRange;
-
-        const vectorList& refCellPoints = refCell.vertexPositions();
-
-        forAll(molCloud_.realFacesWithinRCutMaxOfAnyReferringPatch(), rCF)
-        {
-            const label f
-            (
-                molCloud_.realFacesWithinRCutMaxOfAnyReferringPatch()[rCF]
-            );
-
-            if (molCloud_.testPointFaceDistance(refCellPoints,f))
-            {
-                const label cellO(mesh.faceOwner()[f]);
-
-                if (findIndex(realCellsFoundInRange, cellO) == -1)
-                {
-                    realCellsFoundInRange.append(cellO);
-                }
-
-                if (mesh.isInternalFace(f))
-                {
-                    // boundary faces will not have neighbour information
-
-                    const label cellN(mesh.faceNeighbour()[f]);
-
-                    if (findIndex(realCellsFoundInRange, cellN) == -1)
-                    {
-                        realCellsFoundInRange.append(cellN);
-                    }
-                }
-            }
-        }
-
-        forAll(molCloud_.realPointsWithinRCutMaxOfAnyReferringPatch(), rCP)
-        {
-            const label p
-            (
-                molCloud_.realPointsWithinRCutMaxOfAnyReferringPatch()[rCP]
-            );
-
-            if (molCloud_.testPointFaceDistance(p,refCell))
-            {
-                const labelList& pCells(mesh.pointCells()[p]);
-
-                forAll(pCells, pC)
-                {
-                    const label cellI(pCells[pC]);
-
-                    if (findIndex(realCellsFoundInRange, cellI) == -1)
-                    {
-                        realCellsFoundInRange.append(cellI);
-                    }
-                }
-            }
-        }
-
-
-        const edgeList& refCellEdges = refCell.edges();
-
-        forAll(molCloud_.realEdgesWithinRCutMaxOfAnyReferringPatch(), rCE)
-        {
-            const label edgeIIndex
-            (
-                molCloud_.realEdgesWithinRCutMaxOfAnyReferringPatch()[rCE]
-            );
-
-            const edge& eI(mesh.edges()[edgeIIndex]);
-
-            forAll(refCellEdges, rCE)
-            {
-                const edge& eJ(refCellEdges[rCE]);
-
-                if
-                (
-                    molCloud_.testEdgeEdgeDistance
-                    (
-                        eI,
-                        refCellPoints[eJ.start()],
-                        refCellPoints[eJ.end()]
-                    )
-                )
-                {
-                    const labelList& eICells(mesh.edgeCells()[edgeIIndex]);
-
-                    forAll(eICells, eIC)
-                    {
-                        const label cellI(eICells[eIC]);
-
-                        if (findIndex(realCellsFoundInRange, cellI) == -1)
-                        {
-                            realCellsFoundInRange.append(cellI);
-                        }
-                    }
-                }
-            }
-        }
-
-//         scalar rCutMaxSqr = molCloud_.rCutMax()*molCloud_.rCutMax();
-//
-//         forAll (molCloud_.mesh().points(), pointIIndex)
-//         {
-//             const point& ptI
-//             (
-//                 molCloud_.mesh().points()[pointIIndex]
-//             );
-//
-//             forAll(refCellPoints, rCP)
-//             {
-//                 if (magSqr(ptI - refCellPoints[rCP]) <= rCutMaxSqr)
-//                 {
-//                     const labelList& ptICells
-//                     (
-//                         molCloud_.mesh().pointCells()[pointIIndex]
-//                     );
-//
-//                     forAll(ptICells, pIC)
-//                     {
-//                         const label cellI(ptICells[pIC]);
-//
-//                         if(findIndex(realCellsFoundInRange, cellI) == -1)
-//                         {
-//                             realCellsFoundInRange.append(cellI);
-//                         }
-//                     }
-//                 }
-//             }
-//         }
-
-        refCell.realCells() = realCellsFoundInRange.shrink();
-    }
-}
-
-
-void referredCellList::referMolecules()
-{
-    // Create referred molecules for sending using cell occupancy and
-    // cellSendingReferralLists
-
-    const List<DynamicList<molecule*> >& cellOccupancy
-    (
-        molCloud_.cellOccupancy()
-    );
-
-    forAll(molCloud_.cellSendingReferralLists(), cSRL)
-    {
-        const sendingReferralList& sRL
-        (
-            molCloud_.cellSendingReferralLists()[cSRL]
-        );
-
-        List<DynamicList<referredMolecule> > molsToReferOut(sRL.size());
-
-        forAll(sRL, sRLI)
-        {
-            List<molecule*> realMols = cellOccupancy[sRL[sRLI]];
-
-            forAll (realMols, rM)
-            {
-                molecule* mol = realMols[rM];
-
-                molsToReferOut[sRLI].append
-                (
-                    referredMolecule
-                    (
-                        mol->id(),
-                        mol->position()
-                    )
-                );
-            }
-
-            molsToReferOut[sRLI].shrink();
-        }
-
-        // Send lists of referred molecules to other processors
-
-        if (sRL.destinationProc() != Pstream::myProcNo())
-        {
-            if (Pstream::parRun())
-            {
-                OPstream toInteractingProc
-                (
-                    Pstream::blocking,
-                    sRL.destinationProc()
-                );
-
-                toInteractingProc << molsToReferOut;
-            }
-        }
-        else
-        {
-            // Refer molecules directly for referred cells on the same
-            // processor.
-
-            const receivingReferralList& rRL
-            (
-                molCloud_.cellReceivingReferralLists()[cSRL]
-            );
-
-            forAll(rRL, rRLI)
-            {
-                forAll(rRL[rRLI], rC)
-                {
-                    referredCell& refCellToRefMolsTo = (*this)[rRL[rRLI][rC]];
-
-                    refCellToRefMolsTo.referInMols(molsToReferOut[rRLI]);
-                }
-            }
-        }
-    }
-
-    // Receive referred molecule lists to and distribute to referredCells
-    // according tocellReceivingReferralLists, referredCells deal with the
-    // transformations themselves
-
-    forAll(molCloud_.cellReceivingReferralLists(), cRRL)
-    {
-        const receivingReferralList& rRL
-        (
-            molCloud_.cellReceivingReferralLists()[cRRL]
-        );
-
-        List<List<referredMolecule> > molsToReferIn(rRL.size());
-
-        if (rRL.sourceProc() != Pstream::myProcNo())
-        {
-            if (Pstream::parRun())
-            {
-                IPstream fromInteractingProc
-                (
-                    Pstream::blocking,
-                    rRL.sourceProc()
-                );
-
-                fromInteractingProc >> molsToReferIn;
-            }
-
-            forAll(rRL, rRLI)
-            {
-                forAll(rRL[rRLI], rC)
-                {
-                    referredCell& refCellToRefMolsTo = (*this)[rRL[rRLI][rC]];
-
-                    refCellToRefMolsTo.referInMols(molsToReferIn[rRLI]);
-                }
-            }
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

src/lagrangian/molecularDynamics/potential/Make/files

@@ -1,3 +1,7 @@
+potential = potential
+
+$(potential)/potential.C
+
 pairPotential = pairPotential
 
 $(pairPotential)/pairPotentialList/pairPotentialList.C
@@ -9,6 +13,9 @@ $(pairPotential)/basic/newPairPotential.C
 $(pairPotential)/derived/lennardJones/lennardJones.C
 $(pairPotential)/derived/maitlandSmith/maitlandSmith.C
 $(pairPotential)/derived/azizChen/azizChen.C
+$(pairPotential)/derived/exponentialRepulsion/exponentialRepulsion.C
+$(pairPotential)/derived/coulomb/coulomb.C
+$(pairPotential)/derived/noInteraction/noInteraction.C
 
 energyScalingFunction = energyScalingFunction
 
@@ -30,6 +37,10 @@ $(tetherPotential)/basic/newTetherPotential.C
 
 $(tetherPotential)/derived/harmonicSpring/harmonicSpring.C
 $(tetherPotential)/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.C
+$(tetherPotential)/derived/pitchForkRing/pitchForkRing.C
+
+electrostaticPotential = electrostaticPotential
+
+$(electrostaticPotential)/electrostaticPotential.C
 
 LIB = $(FOAM_LIBBIN)/libpotential
-

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

@@ -22,42 +22,32 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-\*----------------------------------------------------------------------------*/
+\*---------------------------------------------------------------------------*/
 
-#include "moleculeCloud.H"
+#include "electrostaticPotential.H"
+#include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-{
 
-void moleculeCloud::buildCellOccupancy()
+Foam::electrostaticPotential::electrostaticPotential()
+:
+    prefactor(1.0/(4.0 * mathematicalConstant::pi * 8.854187817e-12))
+{}
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::scalar Foam::electrostaticPotential::energy(const scalar r) const
 {
-    forAll(cellOccupancy_, cO)
+    return prefactor/r;
-    {
-        cellOccupancy_[cO].clear();
 }
 
-    iterator mol(this->begin());
 
-    for
+Foam::scalar Foam::electrostaticPotential::force(const scalar r) const
-    (
-        mol = this->begin();
-        mol != this->end();
-        ++mol
-    )
 {
-        cellOccupancy_[mol().cell()].append(&mol());
+    return prefactor/(r*r);
 }
 
-    forAll(cellOccupancy_, cO)
-    {
-        cellOccupancy_[cO].shrink();
-    }
-
-    referredInteractionList_.referMolecules();
-}
-
-} // End namespace Foam
 
 // ************************************************************************* //

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

@@ -22,38 +22,66 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-\*----------------------------------------------------------------------------*/
+Class
+    Foam::pairPotentials::electrostaticPotential
 
-#include "moleculeCloud.H"
+Description
+
+
+SourceFiles
+    electrostaticPotential.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef electrostaticPotential_H
+#define electrostaticPotential_H
+
+#include "scalar.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-void Foam::moleculeCloud::calculateForce()
+namespace Foam
 {
-    buildCellOccupancy();
 
-    iterator mol(this->begin());
+/*---------------------------------------------------------------------------*\
+                           Class electrostaticPotential Declaration
+\*---------------------------------------------------------------------------*/
 
-    // Set all accumulated quantities to zero
+class electrostaticPotential
-    for (mol = this->begin(); mol != this->end(); ++mol)
 {
-        mol().A() = vector::zero;
+    // Private data
 
-        mol().potentialEnergy() = 0.0;
+    scalar prefactor;
 
-        mol().rf() = tensor::zero;
-    }
 
-    calculatePairForce();
+public:
 
-    calculateTetherForce();
+    // Constructors
 
-    calculateExternalForce();
+        //- Construct and set prefactor
-}
+        electrostaticPotential();
 
 
+    // Destructor
+
+        ~electrostaticPotential()
+        {}
+
+
+    // Member Functions
+
+        scalar energy(const scalar r) const;
+
+        scalar force(const scalar r) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
 // ************************************************************************* //
-
-
-
-