From a20caaa2affd3fa8bd0aad15e01c6605dc29c9f9 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 22 Jul 2011 10:34:34 +0100
Subject: [PATCH 1/2] ENH:molecule: reverted to old molecule.

Removed MoleculeCloud, controllers, monoatomic, polyatomic molecules
---
 .../Templates/MoleculeCloud/MoleculeCloud.C   | 1265 -----------------
 .../Templates/MoleculeCloud/MoleculeCloud.H   |  252 ----
 .../Templates/MoleculeCloud/MoleculeCloudI.H  |  416 ------
 .../baseClasses/moleculeCloud/moleculeCloud.C |   48 -
 .../baseClasses/moleculeCloud/moleculeCloud.H |   83 --
 .../derived/monoatomicCloud/monoatomicCloud.H |   52 -
 .../derived/polyatomicCloud/polyatomicCloud.H |   52 -
 .../basic/controllers/controllers.C           |  573 --------
 .../basic/controllers/controllers.H           |  163 ---
 .../basic/controllers/controllersI.H          |   64 -
 .../basic/fluxController/fluxController.C     |  524 -------
 .../basic/fluxController/fluxController.H     |  272 ----
 .../basic/stateController/stateController.C   |  490 -------
 .../basic/stateController/stateController.H   |  270 ----
 .../molecules/constPropSite/constPropSite.C   |  114 --
 .../molecules/constPropSite/constPropSite.H   |  185 ---
 .../molecules/constPropSite/constPropSiteI.H  |  139 --
 .../molecules/monoatomic/monoatomic.C         |  199 ---
 .../molecules/monoatomic/monoatomic.H         |  418 ------
 .../molecules/monoatomic/monoatomicI.H        |  328 -----
 .../molecules/monoatomic/monoatomicIO.C       |  305 ----
 .../molecules/polyatomic/polyatomic.C         |  320 -----
 .../molecules/polyatomic/polyatomic.H         |  483 -------
 .../molecules/polyatomic/polyatomicI.H        |  653 ---------
 .../molecules/polyatomic/polyatomicIO.C       |  414 ------
 25 files changed, 8082 deletions(-)
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.C
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.H
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloudI.H
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/clouds/baseClasses/moleculeCloud/moleculeCloud.C
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/clouds/baseClasses/moleculeCloud/moleculeCloud.H
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/clouds/derived/monoatomicCloud/monoatomicCloud.H
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/clouds/derived/polyatomicCloud/polyatomicCloud.H
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.C
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.H
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllersI.H
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.C
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.H
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.C
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.H
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.C
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.H
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSiteI.H
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomic.C
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomic.H
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicI.H
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicIO.C
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.C
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.H
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomicI.H
 delete mode 100644 src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomicIO.C

diff --git a/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.C b/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.C
deleted file mode 100644
index c6d54036f7..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.C
+++ /dev/null
@@ -1,1265 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*----------------------------------------------------------------------------*/
-
-#include "MoleculeCloud.H"
-#include "fvMesh.H"
-#include "mathematicalConstants.H"
-
-using namespace Foam::constant::mathematical;
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template<class MoleculeType>
-void Foam::MoleculeCloud<MoleculeType>::buildConstProps()
-{
-    const List<word>& idList(pot_.idList());
-
-    constPropList_.setSize(idList.size());
-
-    const List<word>& siteIdList(pot_.siteIdList());
-
-    IOdictionary molPropertiesDict
-    (
-        IOobject
-        (
-            this->name() + "Properties",
-            mesh_.time().constant(),
-            mesh_,
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE,
-            false
-        )
-    );
-
-    Info<< nl << "Reading " << molPropertiesDict.name() << endl;
-
-    forAll(idList, i)
-    {
-        const word& id = idList[i];
-
-        const dictionary& molDict(molPropertiesDict.subDict(id));
-
-        List<word> siteIdNames = molDict.lookup("siteIds");
-
-        List<label> siteIds(siteIdNames.size());
-
-        forAll(siteIdNames, sI)
-        {
-            const word& siteId = siteIdNames[sI];
-
-            siteIds[sI] = findIndex(siteIdList, siteId);
-
-            if (siteIds[sI] == -1)
-            {
-                FatalErrorIn
-                (
-                    "void Foam::MoleculeCloud<MoleculeType>::buildConstProps()"
-                )
-                    << siteId << " site not found."
-                    << nl << abort(FatalError);
-            }
-        }
-
-        constPropList_[i] = typename MoleculeType::constantProperties
-        (
-            molDict,
-            siteIds
-        );
-    }
-}
-
-
-template<class MoleculeType>
-void Foam::MoleculeCloud<MoleculeType>::setSiteSizesAndPositions()
-{
-    forAllIter(typename MoleculeCloud<MoleculeType>, *this, mol)
-    {
-        const typename MoleculeType::constantProperties& cP
-        (
-            constProps(mol().id())
-        );
-
-        mol().setSiteSizes(cP.nSites());
-
-        mol().setSitePositions(cP);
-    }
-}
-
-
-template<class MoleculeType>
-void Foam::MoleculeCloud<MoleculeType>::buildCellOccupancy()
-{
-    forAll(cellOccupancy_, cO)
-    {
-        cellOccupancy_[cO].clear();
-    }
-
-    forAllIter(typename MoleculeCloud<MoleculeType>, *this, mol)
-    {
-        cellOccupancy_[mol().cell()].append(&mol());
-    }
-
-    forAll(cellOccupancy_, cO)
-    {
-        cellOccupancy_[cO].shrink();
-    }
-}
-
-
-template<class MoleculeType>
-void Foam::MoleculeCloud<MoleculeType>::calculatePairForce()
-{
-    PstreamBuffers pBufs(Pstream::nonBlocking);
-
-    // Start sending referred data
-    il_.sendReferredData(cellOccupancy(), pBufs);
-
-    MoleculeType* molI = NULL;
-    MoleculeType* molJ = NULL;
-
-    {
-        // Real-Real interactions
-
-        const labelListList& dil = il_.dil();
-
-        forAll(dil, d)
-        {
-            forAll(cellOccupancy_[d],cellIMols)
-            {
-                molI = cellOccupancy_[d][cellIMols];
-
-                forAll(dil[d], interactingCells)
-                {
-                    List<MoleculeType*> cellJ =
-                        cellOccupancy_[dil[d][interactingCells]];
-
-                    forAll(cellJ, cellJMols)
-                    {
-                        molJ = cellJ[cellJMols];
-
-                        evaluatePair(*molI, *molJ);
-                    }
-                }
-
-                forAll(cellOccupancy_[d], cellIOtherMols)
-                {
-                    molJ = cellOccupancy_[d][cellIOtherMols];
-
-                    if (molJ > molI)
-                    {
-                        evaluatePair(*molI, *molJ);
-                    }
-                }
-            }
-        }
-    }
-
-    // Receive referred data
-    il_.receiveReferredData(pBufs);
-
-    {
-        // Real-Referred interactions
-
-        const labelListList& ril = il_.ril();
-
-        List<IDLList<MoleculeType> >& referredMols = il_.referredParticles();
-
-        forAll(ril, r)
-        {
-            const List<label>& realCells = ril[r];
-
-            IDLList<MoleculeType>& refMols = referredMols[r];
-
-            forAllIter
-            (
-                typename IDLList<MoleculeType>,
-                refMols,
-                refMol
-            )
-            {
-                forAll(realCells, rC)
-                {
-                    List<MoleculeType*> cellI = cellOccupancy_[realCells[rC]];
-
-                    forAll(cellI, cellIMols)
-                    {
-                        molI = cellI[cellIMols];
-
-                        evaluatePair(*molI, refMol());
-                    }
-                }
-            }
-        }
-    }
-}
-
-
-template<class MoleculeType>
-void Foam::MoleculeCloud<MoleculeType>::calculateTetherForce()
-{
-    const tetherPotentialList& tetherPot(pot_.tetherPotentials());
-
-    forAllIter(typename MoleculeCloud<MoleculeType>, *this, mol)
-    {
-        if (mol().tethered())
-        {
-            vector rIT = mol().position() - mol().specialPosition();
-
-            label idI = mol().id();
-
-            scalar massI = constProps(idI).mass();
-
-            vector fIT = tetherPot.force(idI, rIT);
-
-            mol().a() += fIT/massI;
-
-            mol().potentialEnergy() += tetherPot.energy(idI, rIT);
-
-            // What to do here?
-            // mol().rf() += rIT*fIT;
-        }
-    }
-}
-
-
-template<class MoleculeType>
-void Foam::MoleculeCloud<MoleculeType>::calculateExternalForce()
-{
-    forAllIter(typename MoleculeCloud<MoleculeType>, *this, mol)
-    {
-        mol().a() += pot_.gravity();
-    }
-}
-
-
-template<class MoleculeType>
-void Foam::MoleculeCloud<MoleculeType>::removeHighEnergyOverlaps()
-{
-    Info<< nl << "Removing high energy overlaps, limit = "
-        << pot_.potentialEnergyLimit()
-        << nl << "Removal order:";
-
-    forAll(pot_.removalOrder(), rO)
-    {
-        Info<< ' ' << pot_.idList()[pot_.removalOrder()[rO]];
-    }
-
-    Info<< nl ;
-
-    label initialSize = this->size();
-
-    buildCellOccupancy();
-
-    // Real-Real interaction
-
-    MoleculeType* molI = NULL;
-    MoleculeType* molJ = NULL;
-
-    {
-        DynamicList<MoleculeType*> molsToDelete;
-
-        const labelListList& dil(il_.dil());
-
-        forAll(dil, d)
-        {
-            forAll(cellOccupancy_[d], cellIMols)
-            {
-                molI = cellOccupancy_[d][cellIMols];
-
-                forAll(dil[d], interactingCells)
-                {
-                    List<MoleculeType*> cellJ =
-                        cellOccupancy_[dil[d][interactingCells]];
-
-                    forAll(cellJ, cellJMols)
-                    {
-                        molJ = cellJ[cellJMols];
-
-                        if (evaluatePotentialLimit(*molI, *molJ))
-                        {
-                            label idI = molI->id();
-
-                            label idJ = molJ->id();
-
-                            if
-                            (
-                                idI == idJ
-                             || findIndex(pot_.removalOrder(), idJ)
-                              < findIndex(pot_.removalOrder(), idI)
-                            )
-                            {
-                                if (findIndex(molsToDelete, molJ) == -1)
-                                {
-                                    molsToDelete.append(molJ);
-                                }
-                            }
-                            else if (findIndex(molsToDelete, molI) == -1)
-                            {
-                                molsToDelete.append(molI);
-                            }
-                        }
-                    }
-                }
-
-                forAll(cellOccupancy_[d], cellIOtherMols)
-                {
-                    molJ = cellOccupancy_[d][cellIOtherMols];
-
-                    if (molJ > molI)
-                    {
-                        if (evaluatePotentialLimit(*molI, *molJ))
-                        {
-                            label idI = molI->id();
-
-                            label idJ = molJ->id();
-
-                            if
-                            (
-                                idI == idJ
-                             || findIndex(pot_.removalOrder(), idJ)
-                              < findIndex(pot_.removalOrder(), idI)
-                            )
-                            {
-                                if (findIndex(molsToDelete, molJ) == -1)
-                                {
-                                    molsToDelete.append(molJ);
-                                }
-                            }
-                            else if (findIndex(molsToDelete, molI) == -1)
-                            {
-                                molsToDelete.append(molI);
-                            }
-                        }
-                    }
-                }
-            }
-        }
-
-        forAll(molsToDelete, mTD)
-        {
-            deleteParticle(*(molsToDelete[mTD]));
-        }
-    }
-
-    buildCellOccupancy();
-
-    PstreamBuffers pBufs(Pstream::nonBlocking);
-
-    // Start sending referred data
-    il_.sendReferredData(cellOccupancy(), pBufs);
-
-        // Receive referred data
-    il_.receiveReferredData(pBufs);
-
-    // Real-Referred interaction
-
-    {
-        DynamicList<MoleculeType*> molsToDelete;
-
-        const labelListList& ril(il_.ril());
-
-        List<IDLList<MoleculeType> >& referredMols = il_.referredParticles();
-
-        forAll(ril, r)
-        {
-            IDLList<MoleculeType>& refMols = referredMols[r];
-
-            forAllIter
-            (
-                typename IDLList<MoleculeType>,
-                refMols,
-                refMol
-            )
-            {
-                molJ = &refMol();
-
-                const List<label>& realCells = ril[r];
-
-                forAll(realCells, rC)
-                {
-                    label cellI = realCells[rC];
-
-                    List<MoleculeType*> cellIMols = cellOccupancy_[cellI];
-
-                    forAll(cellIMols, cIM)
-                    {
-                        molI = cellIMols[cIM];
-
-                        if (evaluatePotentialLimit(*molI, *molJ))
-                        {
-                            label idI = molI->id();
-
-                            label idJ = molJ->id();
-
-                            if
-                            (
-                                findIndex(pot_.removalOrder(), idI)
-                              < findIndex(pot_.removalOrder(), idJ)
-                            )
-                            {
-                                if (findIndex(molsToDelete, molI) == -1)
-                                {
-                                    molsToDelete.append(molI);
-                                }
-                            }
-                            else if
-                            (
-                                findIndex(pot_.removalOrder(), idI)
-                             == findIndex(pot_.removalOrder(), idJ)
-                            )
-                            {
-                                // Remove one of the molecules
-                                // arbitrarily, assuring that a
-                                // consistent decision is made for
-                                // both real-referred pairs.
-
-                                if (molI->origId() > molJ->origId())
-                                {
-                                    if (findIndex(molsToDelete, molI) == -1)
-                                    {
-                                        molsToDelete.append(molI);
-                                    }
-                                }
-                            }
-                        }
-                    }
-                }
-            }
-        }
-
-        forAll(molsToDelete, mTD)
-        {
-            deleteParticle(*(molsToDelete[mTD]));
-        }
-    }
-
-    buildCellOccupancy();
-
-    // Start sending referred data
-    il_.sendReferredData(cellOccupancy(), pBufs);
-
-    // Receive referred data
-    il_.receiveReferredData(pBufs);
-
-    label molsRemoved = initialSize - this->size();
-
-    if (Pstream::parRun())
-    {
-        reduce(molsRemoved, sumOp<label>());
-    }
-
-    Info<< tab << molsRemoved << " molecules removed" << endl;
-}
-
-
-template<class MoleculeType>
-void Foam::MoleculeCloud<MoleculeType>::initialiseMolecules
-(
-    const dictionary& mdInitialiseDict
-)
-{
-    Info<< nl
-        << "Initialising molecules in each zone specified in "
-        << mdInitialiseDict.name()
-        << endl;
-
-    const cellZoneMesh& cellZones = mesh_.cellZones();
-
-    if (!cellZones.size())
-    {
-        FatalErrorIn
-        (
-            "void Foam::MoleculeCloud<MoleculeType>::initialiseMolecules"
-        )
-            << "No cellZones found in the mesh."
-            << abort(FatalError);
-    }
-
-    forAll(cellZones, z)
-    {
-        const cellZone& zone(cellZones[z]);
-
-        if (zone.size())
-        {
-            if (!mdInitialiseDict.found(zone.name()))
-            {
-                Info<< "No specification subDictionary for zone "
-                    << zone.name() << " found, skipping." << endl;
-            }
-            else
-            {
-                const dictionary& zoneDict =
-                    mdInitialiseDict.subDict(zone.name());
-
-                const scalar temperature
-                (
-                    readScalar(zoneDict.lookup("temperature"))
-                );
-
-                const vector bulkVelocity(zoneDict.lookup("bulkVelocity"));
-
-                List<word> latticeIds
-                (
-                    zoneDict.lookup("latticeIds")
-                );
-
-                List<vector> latticePositions
-                (
-                    zoneDict.lookup("latticePositions")
-                );
-
-                if (latticeIds.size() != latticePositions.size())
-                {
-                    FatalErrorIn
-                    (
-                        "Foam::MoleculeCloud<MoleculeType>::initialiseMolecules"
-                    )
-                        << "latticeIds and latticePositions must be the same "
-                        << " size." << nl
-                        << abort(FatalError);
-                }
-
-                diagTensor latticeCellShape
-                (
-                    zoneDict.lookup("latticeCellShape")
-                );
-
-                scalar latticeCellScale = 0.0;
-
-                if (zoneDict.found("numberDensity"))
-                {
-                    scalar numberDensity = readScalar
-                    (
-                        zoneDict.lookup("numberDensity")
-                    );
-
-                    if (numberDensity < VSMALL)
-                    {
-                        WarningIn
-                        (
-                            "MoleculeCloud<MoleculeType>::initialiseMolecules"
-                        )
-                            << "numberDensity too small, not filling zone "
-                            << zone.name() << endl;
-
-                        continue;
-                    }
-
-                    latticeCellScale = pow
-                    (
-                        latticeIds.size()/(det(latticeCellShape)*numberDensity),
-                        (1.0/3.0)
-                    );
-                }
-                else if (zoneDict.found("massDensity"))
-                {
-                    scalar unitCellMass = 0.0;
-
-                    forAll(latticeIds, i)
-                    {
-                        label id = findIndex(pot_.idList(), latticeIds[i]);
-
-                        const typename MoleculeType::constantProperties& cP
-                        (
-                            constProps(id)
-                        );
-
-                        unitCellMass += cP.mass();
-                    }
-
-                    scalar massDensity = readScalar
-                    (
-                        zoneDict.lookup("massDensity")
-                    );
-
-                    if (massDensity < VSMALL)
-                    {
-                        WarningIn
-                        (
-                            "MoleculeCloud<MoleculeType>::initialiseMolecules"
-                        )
-                            << "massDensity too small, not filling zone "
-                            << zone.name() << endl;
-
-                        continue;
-                    }
-
-                    latticeCellScale = pow
-                    (
-                        unitCellMass/(det(latticeCellShape)*massDensity),
-                        (1.0/3.0)
-                    );
-                }
-                else
-                {
-                    FatalErrorIn
-                    (
-                        "Foam::MoleculeCloud<MoleculeType>::initialiseMolecules"
-                    )
-                        << "massDensity or numberDensity not specified " << nl
-                        << abort(FatalError);
-                }
-
-                latticeCellShape *= latticeCellScale;
-
-                vector anchor(zoneDict.lookup("anchor"));
-
-                bool tethered = false;
-
-                if (zoneDict.found("tetherSiteIds"))
-                {
-                    tethered = bool
-                    (
-                        List<word>(zoneDict.lookup("tetherSiteIds")).size()
-                    );
-                }
-
-                const vector orientationAngles
-                (
-                    zoneDict.lookup("orientationAngles")
-                );
-
-                scalar phi(orientationAngles.x()*pi/180.0);
-
-                scalar theta(orientationAngles.y()*pi/180.0);
-
-                scalar psi(orientationAngles.z()*pi/180.0);
-
-                const tensor R
-                (
-                    cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi),
-                    cos(psi)*sin(phi) + cos(theta)*cos(phi)*sin(psi),
-                    sin(psi)*sin(theta),
-                    - sin(psi)*cos(phi) - cos(theta)*sin(phi)*cos(psi),
-                    - sin(psi)*sin(phi) + cos(theta)*cos(phi)*cos(psi),
-                    cos(psi)*sin(theta),
-                    sin(theta)*sin(phi),
-                    - sin(theta)*cos(phi),
-                    cos(theta)
-                );
-
-                // Find the optimal anchor position.  Finding the approximate
-                // mid-point of the zone of cells and snapping to the nearest
-                // lattice location.
-
-                vector zoneMin = VGREAT*vector::one;
-
-                vector zoneMax = -VGREAT*vector::one;
-
-                forAll(zone, cell)
-                {
-                    const point cellCentre = mesh_.cellCentres()[zone[cell]];
-
-                    if (cellCentre.x() > zoneMax.x())
-                    {
-                        zoneMax.x() = cellCentre.x();
-                    }
-                    if (cellCentre.x() < zoneMin.x())
-                    {
-                        zoneMin.x() = cellCentre.x();
-                    }
-                    if (cellCentre.y() > zoneMax.y())
-                    {
-                        zoneMax.y() = cellCentre.y();
-                    }
-                    if (cellCentre.y() < zoneMin.y())
-                    {
-                        zoneMin.y() = cellCentre.y();
-                    }
-                    if (cellCentre.z() > zoneMax.z())
-                    {
-                        zoneMax.z() = cellCentre.z();
-                    }
-                    if (cellCentre.z() < zoneMin.z())
-                    {
-                        zoneMin.z() = cellCentre.z();
-                    }
-                }
-
-                point zoneMid = 0.5*(zoneMin + zoneMax);
-
-                point latticeMid = inv(latticeCellShape) & (R.T()
-                & (zoneMid - anchor));
-
-                point latticeAnchor
-                (
-                    label(latticeMid.x() + 0.5*sign(latticeMid.x())),
-                    label(latticeMid.y() + 0.5*sign(latticeMid.y())),
-                    label(latticeMid.z() + 0.5*sign(latticeMid.z()))
-                );
-
-                anchor += (R & (latticeCellShape & latticeAnchor));
-
-                // Continue trying to place molecule as long as at
-                // least one molecule is placed in each iteration.
-                // The "|| totalZoneMols == 0" condition means that the
-                // algorithm will continue if the origin is outside the
-                // zone.
-
-                label n = 0;
-
-                label totalZoneMols = 0;
-
-                label molsPlacedThisIteration = 0;
-
-                while
-                (
-                    molsPlacedThisIteration != 0
-                 || totalZoneMols == 0
-                )
-                {
-                    label sizeBeforeIteration = this->size();
-
-                    bool partOfLayerInBounds = false;
-
-                    if (n == 0)
-                    {
-                        // Special treatment is required for the first position,
-                        // i.e. iteration zero.
-
-                        labelVector unitCellLatticePosition(0,0,0);
-
-                        forAll(latticePositions, p)
-                        {
-                            label id = findIndex(pot_.idList(), latticeIds[p]);
-
-                            const vector& latticePosition =
-                                vector
-                                (
-                                    unitCellLatticePosition.x(),
-                                    unitCellLatticePosition.y(),
-                                    unitCellLatticePosition.z()
-                                )
-                              + latticePositions[p];
-
-                            point globalPosition =
-                                anchor
-                              + (R & (latticeCellShape & latticePosition));
-
-                            partOfLayerInBounds = mesh_.bounds().contains
-                            (
-                                globalPosition
-                            );
-
-                            label cell = -1;
-                            label tetFace = -1;
-                            label tetPt = -1;
-
-                            mesh_.findCellFacePt
-                            (
-                                globalPosition,
-                                cell,
-                                tetFace,
-                                tetPt
-                            );
-
-                            if (findIndex(zone, cell) != -1)
-                            {
-                                createMolecule
-                                (
-                                    globalPosition,
-                                    cell,
-                                    tetFace,
-                                    tetPt,
-                                    id,
-                                    tethered,
-                                    temperature,
-                                    bulkVelocity
-                                );
-                            }
-                        }
-                    }
-                    else
-                    {
-                        // Place top and bottom caps.
-
-                        labelVector unitCellLatticePosition(0,0,0);
-
-                        for
-                        (
-                            unitCellLatticePosition.z() = -n;
-                            unitCellLatticePosition.z() <= n;
-                            unitCellLatticePosition.z() += 2*n
-                        )
-                        {
-                            for
-                            (
-                                unitCellLatticePosition.y() = -n;
-                                unitCellLatticePosition.y() <= n;
-                                unitCellLatticePosition.y()++
-                            )
-                            {
-                                for
-                                (
-                                    unitCellLatticePosition.x() = -n;
-                                    unitCellLatticePosition.x() <= n;
-                                    unitCellLatticePosition.x()++
-                                )
-                                {
-                                    forAll(latticePositions, p)
-                                    {
-                                        label id = findIndex
-                                        (
-                                            pot_.idList(),
-                                            latticeIds[p]
-                                        );
-
-                                        const vector& latticePosition =
-                                            vector
-                                            (
-                                                unitCellLatticePosition.x(),
-                                                unitCellLatticePosition.y(),
-                                                unitCellLatticePosition.z()
-                                            )
-                                          + latticePositions[p];
-
-                                        point globalPosition =
-                                            anchor
-                                          + (
-                                                R
-                                              & (
-                                                    latticeCellShape
-                                                  & latticePosition
-                                                )
-                                            );
-
-                                        partOfLayerInBounds =
-                                            mesh_.bounds().contains
-                                            (
-                                                globalPosition
-                                            );
-
-                                        label cell = -1;
-                                        label tetFace = -1;
-                                        label tetPt = -1;
-
-                                        mesh_.findCellFacePt
-                                        (
-                                            globalPosition,
-                                            cell,
-                                            tetFace,
-                                            tetPt
-                                        );
-
-                                        if (findIndex(zone, cell) != -1)
-                                        {
-                                            createMolecule
-                                            (
-                                                globalPosition,
-                                                cell,
-                                                tetFace,
-                                                tetPt,
-                                                id,
-                                                tethered,
-                                                temperature,
-                                                bulkVelocity
-                                            );
-                                        }
-                                    }
-                                }
-                            }
-                        }
-
-                        for
-                        (
-                            unitCellLatticePosition.z() = -(n-1);
-                            unitCellLatticePosition.z() <= (n-1);
-                            unitCellLatticePosition.z()++
-                        )
-                        {
-                            for (label iR = 0; iR <= 2*n -1; iR++)
-                            {
-                                unitCellLatticePosition.x() = n;
-
-                                unitCellLatticePosition.y() = -n + (iR + 1);
-
-                                for (label iK = 0; iK < 4; iK++)
-                                {
-                                    forAll(latticePositions, p)
-                                    {
-                                        label id = findIndex
-                                        (
-                                            pot_.idList(),
-                                            latticeIds[p]
-                                        );
-
-                                        const vector& latticePosition =
-                                            vector
-                                            (
-                                                unitCellLatticePosition.x(),
-                                                unitCellLatticePosition.y(),
-                                                unitCellLatticePosition.z()
-                                            )
-                                          + latticePositions[p];
-
-                                        point globalPosition =
-                                            anchor
-                                          + (
-                                                R
-                                              & (
-                                                   latticeCellShape
-                                                 & latticePosition
-                                                )
-                                            );
-
-                                        partOfLayerInBounds =
-                                            mesh_.bounds().contains
-                                            (
-                                                globalPosition
-                                            );
-
-                                        label cell = -1;
-                                        label tetFace = -1;
-                                        label tetPt = -1;
-
-                                        mesh_.findCellFacePt
-                                        (
-                                            globalPosition,
-                                            cell,
-                                            tetFace,
-                                            tetPt
-                                        );
-
-                                        if (findIndex(zone, cell) != -1)
-                                        {
-                                            createMolecule
-                                            (
-                                                globalPosition,
-                                                cell,
-                                                tetFace,
-                                                tetPt,
-                                                id,
-                                                tethered,
-                                                temperature,
-                                                bulkVelocity
-                                            );
-                                        }
-                                    }
-
-                                    unitCellLatticePosition =
-                                        labelVector
-                                        (
-                                          - unitCellLatticePosition.y(),
-                                            unitCellLatticePosition.x(),
-                                            unitCellLatticePosition.z()
-                                        );
-                                }
-                            }
-                        }
-                    }
-
-                    if
-                    (
-                        totalZoneMols == 0
-                     && !partOfLayerInBounds
-                    )
-                    {
-                        WarningIn
-                        (
-                            "Foam::MoleculeCloud<MoleculeType>::"
-                            "initialiseMolecules()"
-                        )
-                            << "A whole layer of unit cells was placed "
-                            << "outside the bounds of the mesh, but no "
-                            << "molecules have been placed in zone '"
-                            << zone.name()
-                            << "'.  This is likely to be because the zone "
-                            << "has few cells ("
-                            << zone.size()
-                            << " in this case) and no lattice position "
-                            << "fell inside them.  "
-                            << "Aborting filling this zone."
-                            << endl;
-
-                        break;
-                    }
-
-                    molsPlacedThisIteration =
-                        this->size() - sizeBeforeIteration;
-
-                    totalZoneMols += molsPlacedThisIteration;
-
-                    n++;
-                }
-            }
-        }
-    }
-}
-
-
-template<class MoleculeType>
-void Foam::MoleculeCloud<MoleculeType>::createMolecule
-(
-    const point& position,
-    label cell,
-    label tetFace,
-    label tetPt,
-    label id,
-    bool tethered,
-    scalar temperature,
-    const vector& bulkVelocity
-)
-{
-    if (cell == -1)
-    {
-        mesh_.findCellFacePt(position, cell, tetFace, tetPt);
-    }
-
-    if (cell == -1)
-    {
-        FatalErrorIn("Foam::MoleculeCloud<MoleculeType>::createMolecule")
-            << "Position specified does not correspond to a mesh cell." << nl
-            << abort(FatalError);
-    }
-
-    point specialPosition(vector::zero);
-
-    label special = 0;
-
-    if (tethered)
-    {
-        specialPosition = position;
-
-        special = MoleculeType::SPECIAL_TETHERED;
-    }
-
-    const typename MoleculeType::constantProperties& cP(constProps(id));
-
-    typename MoleculeType::trackingData td
-    (
-        *this,
-        MoleculeType::trackingData::tpAccess
-    );
-
-    addParticle
-    (
-        new MoleculeType
-        (
-            mesh_,
-            position,
-            cell,
-            tetFace,
-            tetPt,
-            temperature,
-            bulkVelocity,
-            specialPosition,
-            cP,
-            td,
-            special,
-            id
-        )
-    );
-}
-
-
-template<class MoleculeType>
-Foam::label Foam::MoleculeCloud<MoleculeType>::nSites() const
-{
-    label n = 0;
-
-    forAllConstIter(typename MoleculeCloud<MoleculeType>, *this, mol)
-    {
-        n += constProps(mol().id()).nSites();
-    }
-
-    return n;
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class MoleculeType>
-Foam::MoleculeCloud<MoleculeType>::MoleculeCloud
-(
-    const word& cloudName,
-    const polyMesh& mesh,
-    const potential& pot,
-    bool readFields
-)
-:
-    Cloud<MoleculeType>(mesh, cloudName, false),
-    moleculeCloud(),
-    mesh_(mesh),
-    pot_(pot),
-    cellOccupancy_(mesh_.nCells()),
-    il_(mesh_, pot_.pairPotentials().rCutMax(), false),
-    constPropList_(),
-    rndGen_(label(971501) + 1526*Pstream::myProcNo())
-{
-    if (readFields)
-    {
-        MoleculeType::readFields(*this);
-    }
-
-    buildConstProps();
-
-    setSiteSizesAndPositions();
-
-    removeHighEnergyOverlaps();
-
-    calculateForce();
-}
-
-
-template<class MoleculeType>
-Foam::MoleculeCloud<MoleculeType>::MoleculeCloud
-(
-    const word& cloudName,
-    const polyMesh& mesh,
-    const potential& pot,
-    const dictionary& mdInitialiseDict,
-    bool readFields
-)
-:
-    Cloud<MoleculeType>(mesh, cloudName, false),
-    moleculeCloud(),
-    mesh_(mesh),
-    pot_(pot),
-    il_(mesh_),
-    constPropList_(),
-    rndGen_(label(971501) + 1526*Pstream::myProcNo())
-{
-    if (readFields)
-    {
-        MoleculeType::readFields(*this);
-    }
-
-    this->clear();
-
-    buildConstProps();
-
-    initialiseMolecules(mdInitialiseDict);
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class MoleculeType>
-void Foam::MoleculeCloud<MoleculeType>::evolve()
-{
-    typename MoleculeType::trackingData td0
-    (
-        *this,
-        MoleculeType::trackingData::tpFirstVelocityHalfStep
-    );
-    Cloud<MoleculeType>::move(td0, mesh_.time().deltaTValue());
-
-    typename MoleculeType::trackingData td1
-    (
-        *this,
-        MoleculeType::trackingData::tpLinearTrack
-    );
-    Cloud<MoleculeType>::move(td1, mesh_.time().deltaTValue());
-
-    typename MoleculeType::trackingData td2
-    (
-        *this,
-        MoleculeType::trackingData::tpRotationalTrack
-    );
-    Cloud<MoleculeType>::move(td2, mesh_.time().deltaTValue());
-
-    calculateForce();
-
-    typename MoleculeType::trackingData td3
-    (
-        *this,
-        MoleculeType::trackingData::tpSecondVelocityHalfStep
-    );
-    Cloud<MoleculeType>::move(td3, mesh_.time().deltaTValue());
-
-    info();
-}
-
-
-template<class MoleculeType>
-void Foam::MoleculeCloud<MoleculeType>::calculateForce()
-{
-    buildCellOccupancy();
-
-    // Set accumulated quantities to zero
-    forAllIter(typename MoleculeCloud<MoleculeType>, *this, mol)
-    {
-        mol().siteForces() = vector::zero;
-
-        mol().potentialEnergy() = 0.0;
-
-        mol().rf() = tensor::zero;
-    }
-
-    calculatePairForce();
-
-    calculateTetherForce();
-
-    calculateExternalForce();
-}
-
-
-template<class MoleculeType>
-void Foam::MoleculeCloud<MoleculeType>::info()
-{
-    // Calculates and prints the mean momentum and energy in the system
-    // and the number of molecules.
-
-    typename MoleculeType::trackingData td
-    (
-        *this,
-        MoleculeType::trackingData::tpAccess
-    );
-
-    MoleculeType::info(td);
-}
-
-
-template<class MoleculeType>
-void Foam::MoleculeCloud<MoleculeType>::writeXYZ(const fileName& fName) const
-{
-    OFstream os(fName);
-
-    os  << nSites() << nl
-        << "MoleculeCloud<MoleculeType> site positions in angstroms" << nl;
-
-    forAllConstIter(typename MoleculeCloud<MoleculeType>, *this, mol)
-    {
-        const typename MoleculeType::constantProperties& cP
-        (
-            constProps(mol().id())
-        );
-
-        forAll(mol().sitePositions(), i)
-        {
-            const point& sP = mol().sitePositions()[i];
-
-            os  << pot_.siteIdList()[cP.sites()[i].siteId()]
-                << ' ' << sP.x()*1e10
-                << ' ' << sP.y()*1e10
-                << ' ' << sP.z()*1e10
-                << nl;
-        }
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.H b/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.H
deleted file mode 100644
index 38c912a17d..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.H
+++ /dev/null
@@ -1,252 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::MoleculeCloud
-
-Description
-
-SourceFiles
-    MoleculeCloudI.H
-    MoleculeCloud.C
-
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef MoleculeCloud_H
-#define MoleculeCloud_H
-
-#include "Cloud.H"
-#include "moleculeCloud.H"
-#include "IOdictionary.H"
-#include "potential.H"
-#include "InteractionLists.H"
-#include "labelVector.H"
-#include "Random.H"
-#include "fileName.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                       Class MoleculeCloud Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class MoleculeType>
-class MoleculeCloud
-:
-    public Cloud<MoleculeType>,
-    public moleculeCloud
-{
-
-private:
-
-    // Private data
-
-        //-
-        const polyMesh& mesh_;
-
-        //-
-        const potential& pot_;
-
-        //-
-        List<DynamicList<MoleculeType*> > cellOccupancy_;
-
-        //-
-        InteractionLists<MoleculeType> il_;
-
-        //-
-        List<typename MoleculeType::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
-        (
-            MoleculeType& molI,
-            MoleculeType& molJ
-        );
-
-        //-
-        inline bool evaluatePotentialLimit
-        (
-            MoleculeType& molI,
-            MoleculeType& molJ
-        ) const;
-
-        //-
-        void calculateTetherForce();
-
-        //-
-        void calculateExternalForce();
-
-        //-
-        void removeHighEnergyOverlaps();
-
-        //-
-        void initialiseMolecules(const dictionary& mdInitialiseDict);
-
-        //-
-        void createMolecule
-        (
-            const point& position,
-            label cell,
-            label tetFace,
-            label tetPt,
-            label id,
-            bool tethered,
-            scalar temperature,
-            const vector& bulkVelocity
-        );
-
-        //-
-        label nSites() const;
-
-        //- Disallow default bitwise copy construct
-        MoleculeCloud(const MoleculeCloud&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const MoleculeCloud&);
-
-
-public:
-
-    // Constructors
-
-        //- Construct given mesh and potential references
-        MoleculeCloud
-        (
-            const word& cloudName,
-            const polyMesh& mesh,
-            const potential& pot,
-            bool readFields = true
-        );
-
-        //- Construct given mesh, potential and mdInitialiseDict
-        MoleculeCloud
-        (
-            const word& cloudName,
-            const polyMesh& mesh,
-            const potential& pot,
-            const dictionary& mdInitialiseDict,
-            bool readFields = true
-        );
-
-
-    // Member Functions
-
-        //- Evolve the molecules (move, calculate forces, control state etc)
-        void evolve();
-
-        //-
-        void calculateForce();
-
-        //- Print cloud information
-        void info();
-
-
-        // Access
-
-            //-
-            inline const polyMesh& mesh() const;
-
-            //-
-            inline const potential& pot() const;
-
-            //-
-            inline const List<DynamicList<MoleculeType*> >&
-            cellOccupancy() const;
-
-            //-
-            inline const InteractionLists<MoleculeType>& il() const;
-
-            //-
-            inline const List<typename MoleculeType::constantProperties>
-            constProps() const;
-
-            //-
-            inline const typename MoleculeType::constantProperties&
-            constProps(label id) const;
-
-            //-
-            inline Random& rndGen();
-
-            //-
-            inline vector equipartitionLinearVelocity
-            (
-                scalar temperature,
-                scalar mass
-            );
-
-            //-
-            inline vector equipartitionAngularMomentum
-            (
-                scalar temperature,
-                const typename MoleculeType::constantProperties& cP
-            );
-
-
-    // Member Operators
-
-        //- Write molecule sites in XYZ format
-        void writeXYZ(const fileName& fName) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "MoleculeCloudI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "MoleculeCloud.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloudI.H b/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloudI.H
deleted file mode 100644
index 42c506024f..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloudI.H
+++ /dev/null
@@ -1,416 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "constants.H"
-
-using namespace Foam::constant;
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template<class MoleculeType>
-inline void Foam::MoleculeCloud<MoleculeType>::evaluatePair
-(
-    MoleculeType& molI,
-    MoleculeType& molJ
-)
-{
-    const pairPotentialList& pairPot = pot_.pairPotentials();
-
-    const pairPotential& electrostatic = pairPot.electrostatic();
-
-    label idI = molI.id();
-
-    label idJ = molJ.id();
-
-    const typename MoleculeType::constantProperties& constPropI
-    (
-        constProps(idI)
-    );
-
-    const typename MoleculeType::constantProperties& constPropJ
-    (
-        constProps(idJ)
-    );
-
-    forAll(constPropI.pairPotSites(), pI)
-    {
-        label sI = constPropI.pairPotSites()[pI];
-
-        label idsI = constPropI.sites()[sI].siteId();
-
-        forAll(constPropJ.pairPotSites(), pJ)
-        {
-            label sJ = constPropJ.pairPotSites()[pJ];
-
-            label idsJ = constPropJ.sites()[sJ].siteId();
-
-            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;
-
-                vector rIJ = molI.position() - molJ.position();
-
-                tensor virialContribution =
-                (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
-
-                molI.rf() += virialContribution;
-
-                molJ.rf() += virialContribution;
-            }
-        }
-    }
-
-
-    forAll(constPropI.electrostaticSites(), pI)
-    {
-        label sI = constPropI.electrostaticSites()[pI];
-
-        forAll(constPropJ.electrostaticSites(), pJ)
-        {
-            label sJ = constPropJ.electrostaticSites()[pJ];
-
-            vector rsIsJ =
-                molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
-
-            scalar rsIsJMagSq = magSqr(rsIsJ);
-
-            if (rsIsJMagSq <= electrostatic.rCutSqr())
-            {
-                scalar rsIsJMag = mag(rsIsJ);
-
-                scalar chargeI = constPropI.sites()[sI].siteCharge();
-
-                scalar chargeJ = constPropJ.sites()[sJ].siteCharge();
-
-                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;
-
-                vector rIJ = molI.position() - molJ.position();
-
-                tensor virialContribution =
-                    (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
-
-                molI.rf() += virialContribution;
-
-                molJ.rf() += virialContribution;
-            }
-        }
-    }
-}
-
-
-template<class MoleculeType>
-inline bool Foam::MoleculeCloud<MoleculeType>::evaluatePotentialLimit
-(
-    MoleculeType& molI,
-    MoleculeType& molJ
-) const
-{
-    const pairPotentialList& pairPot = pot_.pairPotentials();
-
-    const pairPotential& electrostatic = pairPot.electrostatic();
-
-    label idI = molI.id();
-
-    label idJ = molJ.id();
-
-    const typename MoleculeType::constantProperties& constPropI
-    (
-        constProps(idI)
-    );
-
-    const typename MoleculeType::constantProperties& constPropJ
-    (
-        constProps(idJ)
-    );
-
-    forAll(constPropI.pairPotSites(), pI)
-    {
-        label sI = constPropI.pairPotSites()[pI];
-
-        label idsI = constPropI.sites()[sI].siteId();
-
-        forAll(constPropJ.pairPotSites(), pJ)
-        {
-            label sJ = constPropJ.pairPotSites()[pJ];
-
-            label idsJ = constPropJ.sites()[sJ].siteId();
-
-            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<MoleculeType>::"
-                        "removeHighEnergyOverlaps()"
-                    )
-                        << "Molecule site pair closer than "
-                        << SMALL
-                        << ": mag separation = " << rsIsJMag
-                        << ". These may have been placed on top of each"
-                        << " other by a rounding error in mdInitialise in"
-                        << " parallel or a block filled with moleculess"
-                        << " 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;
-                };
-            }
-        }
-    }
-
-    forAll(constPropI.electrostaticSites(), pI)
-    {
-        label sI = constPropI.electrostaticSites()[pI];
-
-        forAll(constPropJ.electrostaticSites(), pJ)
-        {
-            label sJ = constPropJ.electrostaticSites()[pJ];
-
-            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<MoleculeType>::"
-                        "removeHighEnergyOverlaps()"
-                    )
-                        << "Molecule site pair closer than "
-                        << SMALL
-                        << ": mag separation = " << rsIsJMag
-                        << ". These may have been placed on top of each"
-                        << " other by a rounding error in mdInitialise in"
-                        << " parallel or a block filled with molecules"
-                        << " twice. Removing one of the molecules."
-                        << endl;
-
-                    return true;
-                }
-
-                if (rsIsJMag < electrostatic.rMin())
-                {
-                    return true;
-                }
-
-                scalar chargeI = constPropI.sites()[sI].siteCharge();
-
-                scalar chargeJ = constPropJ.sites()[sJ].siteCharge();
-
-                if
-                (
-                    mag(chargeI*chargeJ*electrostatic.energy(rsIsJMag))
-                  > pot_.potentialEnergyLimit()
-                )
-                {
-                    return true;
-                };
-            }
-        }
-    }
-
-    return false;
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class MoleculeType>
-inline const Foam::polyMesh& Foam::MoleculeCloud<MoleculeType>::mesh() const
-{
-    return mesh_;
-}
-
-
-template<class MoleculeType>
-inline const Foam::potential& Foam::MoleculeCloud<MoleculeType>::pot() const
-{
-    return pot_;
-}
-
-
-template<class MoleculeType>
-inline const Foam::List<Foam::DynamicList<MoleculeType*> >&
-Foam::MoleculeCloud<MoleculeType>::cellOccupancy() const
-{
-    return cellOccupancy_;
-}
-
-
-template<class MoleculeType>
-inline const Foam::InteractionLists<MoleculeType>&
-Foam::MoleculeCloud<MoleculeType>::il() const
-{
-    return il_;
-}
-
-
-template<class MoleculeType>
-inline const Foam::List<typename MoleculeType::constantProperties>
-Foam::MoleculeCloud<MoleculeType>::constProps() const
-{
-    return constPropList_;
-}
-
-
-template<class MoleculeType>
-inline const typename MoleculeType::constantProperties&
-Foam::MoleculeCloud<MoleculeType>::constProps(label id) const
-{
-    return constPropList_[id];
-}
-
-
-template<class MoleculeType>
-inline Foam::Random& Foam::MoleculeCloud<MoleculeType>::rndGen()
-{
-    return rndGen_;
-}
-
-
-template<class MoleculeType>
-inline Foam::vector
-Foam::MoleculeCloud<MoleculeType>::equipartitionLinearVelocity
-(
-    scalar temperature,
-    scalar mass
-)
-{
-    return sqrt(physicoChemical::k.value()*temperature/mass)*vector
-    (
-        rndGen_.GaussNormal(),
-        rndGen_.GaussNormal(),
-        rndGen_.GaussNormal()
-    );
-}
-
-
-template<class MoleculeType>
-inline Foam::vector
-Foam::MoleculeCloud<MoleculeType>::equipartitionAngularMomentum
-(
-    scalar temperature,
-    const typename MoleculeType::constantProperties& cP
-)
-{
-    scalar sqrtKbT = sqrt(physicoChemical::k.value()*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()
-        );
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/clouds/baseClasses/moleculeCloud/moleculeCloud.C b/src/lagrangian/molecularDynamics/molecule/clouds/baseClasses/moleculeCloud/moleculeCloud.C
deleted file mode 100644
index 4fda04df7d..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/clouds/baseClasses/moleculeCloud/moleculeCloud.C
+++ /dev/null
@@ -1,48 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "moleculeCloud.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(moleculeCloud, 0);
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::moleculeCloud::moleculeCloud()
-{}
-
-
-// * * * * * * * * * * * * * * * *  Destructors  * * * * * * * * * * * * * * //
-
-Foam::moleculeCloud::~moleculeCloud()
-{}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/clouds/baseClasses/moleculeCloud/moleculeCloud.H b/src/lagrangian/molecularDynamics/molecule/clouds/baseClasses/moleculeCloud/moleculeCloud.H
deleted file mode 100644
index c6ba2096d3..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/clouds/baseClasses/moleculeCloud/moleculeCloud.H
+++ /dev/null
@@ -1,83 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::moleculeCloud
-
-Description
-    Virtual abstract base class for templated moleculeCloud
-
-SourceFiles
-    moleculeCloud.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef moleculeCloud_H
-#define moleculeCloud_H
-
-#include "volFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                       Class moleculeCloud Declaration
-\*---------------------------------------------------------------------------*/
-
-class moleculeCloud
-{
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        moleculeCloud(const moleculeCloud&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const moleculeCloud&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("moleculeCloud");
-
-    // Constructors
-
-        //- Null constructor
-        moleculeCloud();
-
-    //- Destructor
-    virtual ~moleculeCloud();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/clouds/derived/monoatomicCloud/monoatomicCloud.H b/src/lagrangian/molecularDynamics/molecule/clouds/derived/monoatomicCloud/monoatomicCloud.H
deleted file mode 100644
index 727bfe226e..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/clouds/derived/monoatomicCloud/monoatomicCloud.H
+++ /dev/null
@@ -1,52 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::monoatomicCloud
-
-Description
-    Cloud class to simulate monoatomic molecules
-
-SourceFiles
-    monoatomicCloud.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef monoatomicCloud_H
-#define monoatomicCloud_H
-
-#include "MoleculeCloud.H"
-#include "monoatomic.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    typedef MoleculeCloud<monoatomic> monoatomicCloud;
-}
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/clouds/derived/polyatomicCloud/polyatomicCloud.H b/src/lagrangian/molecularDynamics/molecule/clouds/derived/polyatomicCloud/polyatomicCloud.H
deleted file mode 100644
index bb902fe3a4..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/clouds/derived/polyatomicCloud/polyatomicCloud.H
+++ /dev/null
@@ -1,52 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::polyatomicCloud
-
-Description
-    Cloud class to simulate polyatomic molecules
-
-SourceFiles
-    polyatomicCloud.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef polyatomicCloud_H
-#define polyatomicCloud_H
-
-#include "MoleculeCloud.H"
-#include "polyatomic.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    typedef MoleculeCloud<polyatomic> polyatomicCloud;
-}
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.C b/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.C
deleted file mode 100644
index f0bad336ee..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.C
+++ /dev/null
@@ -1,573 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "controllers.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::controllers::controllers
-(
-    const polyMesh& mesh
-)
-:
-    time_(mesh.time()),
-    controllersDict_
-    (
-        IOobject
-        (
-            "controllersDict",
-            time_.system(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        )
-    ),
-    stateControllersList_(),
-    sCNames_(),
-    sCIds_(),
-    sCFixedPathNames_(),
-    stateControllers_(),
-    fluxControllersList_(),
-    fCNames_(),
-    fCIds_(),
-    fCFixedPathNames_(),
-    fluxControllers_()
-{}
-
-
-Foam::controllers::controllers
-(
-    const polyMesh& mesh,
-    polyatomicCloud& cloud
-)
-:
-    time_(mesh.time()),
-    controllersDict_
-    (
-        IOobject
-        (
-            "controllersDict",
-            time_.system(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    ),
-    stateControllersList_(controllersDict_.lookup("stateControllers")),
-    sCNames_(stateControllersList_.size()),
-    sCIds_(stateControllersList_.size()),
-    sCFixedPathNames_(stateControllersList_.size()),
-    stateControllers_(stateControllersList_.size()),
-    fluxControllersList_(controllersDict_.lookup("fluxControllers")),
-    fCNames_(fluxControllersList_.size()),
-    fCIds_(fluxControllersList_.size()),
-    fCFixedPathNames_(fluxControllersList_.size()),
-    fluxControllers_(fluxControllersList_.size())
-{
-
-    Info << nl << "Creating controllers" << nl << endl;
-
-    // state controllers
-
-    if (!stateControllers_.empty())
-    {
-        forAll(stateControllers_, sC)
-        {
-            const entry& controllersI = stateControllersList_[sC];
-
-            const dictionary& controllersIDict = controllersI.dict();
-
-            stateControllers_[sC] = autoPtr<stateController>
-            (
-                stateController::New(time_, cloud, controllersIDict)
-            );
-
-            sCNames_[sC] = stateControllers_[sC]->type();
-
-            sCIds_[sC] = sC;
-        }
-    }
-
-    //- flux controllers
-
-    if (!fluxControllers_.empty())
-    {
-        forAll(fluxControllers_, fC)
-        {
-            const entry& controllersI = fluxControllersList_[fC];
-
-            const dictionary& controllersIDict = controllersI.dict();
-
-            fluxControllers_[fC] = autoPtr<fluxController>
-            (
-                fluxController::New(time_, cloud, controllersIDict)
-            );
-
-            fCNames_[fC] = fluxControllers_[fC]->type();
-            fCIds_[fC] = fC;
-        }
-    }
-
-    // creating directories for state controllers
-    if (!nStateControllers_.empty())
-    {
-        // case/controllers
-        fileName controllersPath(time_.path()/"controllers");
-
-        if (!isDir(controllersPath))
-        {
-            mkDir(controllersPath);
-        }
-
-        // case/controllers/<cloudName>
-        fileName controllersPath(controllersPath/cloud.name());
-
-        if (!isDir(controllersPath))
-        {
-            mkDir(controllersPath);
-        }
-
-        // case/controllers/<cloudName>/stateControllers
-        fileName stateControllersPath(controllersPath/"stateControllers");
-
-        if (!isDir(stateControllersPath))
-        {
-            mkDir(stateControllersPath);
-        }
-
-        forAll(stateControllers_, sC)
-        {
-            if (stateControllers_[sC]->writeInCase())
-            {
-                // case/controllers/<cloudName>/
-                //     stateControllers/<stateControllerModel>
-                fileName stateControllerPath(stateControllersPath/sCNames_[sC]);
-
-                if (!isDir(stateControllerPath))
-                {
-                    mkDir(stateControllerPath);
-                }
-
-                const word& regionName = stateControllers_[sC]->regionName();
-
-                // case/controllers/<cloudName>/
-                //     stateControllers/<stateControllerModel>/<cellZoneName>
-                fileName zonePath(stateControllerPath/regionName);
-
-                if (!isDir(zonePath))
-                {
-                    mkDir(zonePath);
-                }
-
-                sCFixedPathNames_[sC] = zonePath;
-            }
-        }
-    }
-
-    // creating directories for flux controllers
-    if (nFluxControllers_ > 0)
-    {
-        // case/controllers
-        fileName controllersPath(time_.path()/"controllers");
-
-        if ( !isDir(controllersPath) )
-        {
-            mkDir(controllersPath);
-        }
-
-        // case/controllers/<cloudName>
-        fileName controllersPath(time_.path()/cloud.name());
-
-        if ( !isDir(controllersPath) )
-        {
-            mkDir(controllersPath);
-        }
-
-        // case/controllers/<cloudName>/fluxControllers
-        fileName fluxControllersPath(controllersPath/"fluxControllers");
-
-        if (!isDir(fluxControllersPath))
-        {
-            mkDir(fluxControllersPath);
-        }
-
-        forAll(fluxControllers_, fC)
-        {
-            if (fluxControllers_[fC]->writeInCase())
-            {
-                // case/controllers/<cloudName>/
-                //    fluxControllers/<fluxControllerModel>
-                fileName fluxControllerPath(fluxControllersPath/fCNames_[fC]);
-
-                if (!isDir(fluxControllerPath))
-                {
-                    mkDir(fluxControllerPath);
-                }
-
-                const word& regionName = fluxControllers_[fC]->regionName();
-
-                // case/controllers/<cloudName>/
-                //    fluxControllers/<fluxControllerModel>/<faceZoneName>
-                fileName zonePath(fluxControllerPath/regionName);
-
-                if (!isDir(zonePath))
-                {
-                    mkDir(zonePath);
-                }
-
-                fCFixedPathNames_[fC] = zonePath;
-            }
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-controllers::~controllers()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void controllers::initialConfig()
-{
-    forAll(stateControllers_, sC)
-    {
-        stateControllers_[sC]->initialConfiguration();
-    }
-
-    forAll(fluxControllers_, fC)
-    {
-        fluxControllers_[fC]->initialConfiguration();
-    }
-}
-
-
-void controllers::updateTimeInfo()
-{
-    forAll(stateControllers_, sC)
-    {
-        stateControllers_[sC]->updateTime();
-    }
-
-    forAll(fluxControllers_, fC)
-    {
-        fluxControllers_[fC]->updateTime();
-    }
-}
-
-
-void controllers::controlState()
-{
-    forAll(stateControllers_, sC)
-    {
-        stateControllers_[sC]->controlMols();
-    }
-}
-
-
-void controllers::controlVelocitiesI()
-{
-    forAll(stateControllers_, sC)
-    {
-        stateControllers_[sC]->controlMolsBeg();
-    }
-}
-
-void controllers::controlVelocitiesII()
-{
-    forAll(stateControllers_, sC)
-    {
-        stateControllers_[sC]->controlMolsEnd();
-    }
-}
-
-
-void  controllers::controlPriorToForces()
-{
-    forAll(stateControllers_, sC)
-    {
-        stateControllers_[sC]->controlBeforeForces();
-    }
-}
-
-
-void controllers::calculateStateProps()
-{
-    forAll(stateControllers_, sC)
-    {
-        stateControllers_[sC]->calculateProperties();
-    }
-
-    forAll(fluxControllers_, fC)
-    {
-        fluxControllers_[fC]->calculateProperties();
-    }
-}
-
-
-void controllers::outputStateResults()
-{
-    const Time& runTime = time_;
-
-    if (runTime.outputTime())
-    {
-        // creating a set of directories in the current time directory
-        {
-            List<fileName> timePathNames(sCFixedPathNames_.size());
-
-            if (nStateControllers_ > 0)
-            {
-                if (Pstream::master())
-                {
-                    // case/<timeDir>/uniform
-                    fileName uniformTimePath
-                    (
-                        runTime.path()/runTime.timeName()/"uniform"
-                    );
-
-                    if (!isDir(uniformTimePath))
-                    {
-                        mkDir(uniformTimePath);
-                    }
-
-                    if (!stateControllers_.empty())
-                    {
-                        // case/<timeDir>/uniform/controllers
-                        fileName controllersTimePath
-                        (
-                            uniformTimePath/"controllers"
-                        );
-
-                        if (!isDir(controllersTimePath))
-                        {
-                            mkDir(controllersTimePath);
-                        }
-
-                        // case/<timeDir>/uniform/controllers/<cloudName>
-                        fileName cloudTimePath
-                        (
-                            controllersTimePath/cloud.name()
-                        );
-
-                        if (!isDir(cloudTimePath))
-                        {
-                            mkDir(cloudTimePath);
-                        }
-
-                        // case/<timeDir>/uniform/controllers/<cloudName>/
-                        fileName stateControllersTimePath
-                        (
-                            cloudTimePath/"stateControllers"
-                        );
-
-                        if (!isDir(stateControllersTimePath))
-                        {
-                            mkDir(stateControllersTimePath);
-                        }
-
-                        forAll(stateControllers_, sC)
-                        {
-                            if (stateControllers_[sC]->writeInTimeDir())
-                            {
-                                // case/<timeDir>/uniform/controllers/
-                                //     <cloudName>/<stateControllerModel>
-                                fileName sCTimePath
-                                (
-                                    stateControllersTimePath/sCNames_[sC]
-                                );
-
-                                if (!isDir(sCTimePath))
-                                {
-                                    mkDir(sCTimePath);
-                                }
-
-                                // Creating directory for different zones but
-                                // of the same model
-                                const word& regionName =
-                                    stateControllers_[sC]->regionName();
-
-                                // case/<timeDir>/uniform/controllers/
-                                //    <cloudName>/<stateControllerModel>/
-                                //    <cellZoneName>
-                                fileName zoneTimePath(sCTimePath/regionName);
-
-                                if (!isDir(zoneTimePath))
-                                {
-                                    mkDir(zoneTimePath);
-                                }
-
-                                timePathNames[sC] = zoneTimePath;
-                            }
-                        }
-                    }
-                }
-            }
-
-            // write out data
-            forAll(stateControllers_, sC)
-            {
-                stateControllers_[sC]->output
-                (
-                    sCFixedPathNames_[sC],
-                    timePathNames[sC]
-                );
-            }
-        }
-
-        {
-            List<fileName> timePathNames(fCFixedPathNames_.size());
-
-            if (nFluxControllers_ > 0)
-            {
-                if (Pstream::master())
-                {
-                    // case/<timeDir>/uniform
-                    fileName uniformTimePath
-                    (
-                        runTime.path()/runTime.timeName()/"uniform"
-                    );
-
-                    if (!isDir(uniformTimePath))
-                    {
-                        mkDir(uniformTimePath);
-                    }
-
-                    if (!fluxControllers_.empty())
-                    {
-                        // case/<timeDir>/uniform/controllers
-                        fileName controllersTimePath
-                        (
-                            uniformTimePath/"controllers"
-                        );
-
-                        if (!isDir(controllersTimePath))
-                        {
-                            mkDir(controllersTimePath);
-                        }
-
-                        // case/<timeDir>/uniform/controllers/<cloudName>
-                        fileName cloudTimePath
-                        (
-                            controllersTimePath/cloud.name()
-                        );
-
-                        if (!isDir(cloudTimePath))
-                        {
-                            mkDir(cloudTimePath);
-                        }
-
-                        // case/<timeDir>/uniform/fluxControllers
-                        fileName controllersTimePath
-                        (
-                            cloudTimePath/"fluxControllers"
-                        );
-
-                        if (!isDir(controllersTimePath))
-                        {
-                            mkDir(controllersTimePath);
-                        }
-
-                        forAll(fluxControllers_, fC)
-                        {
-                            if (stateControllers_[fC]->writeInTimeDir())
-                            {
-                                // case/<timeDir>/uniform/controllers/
-                                //     <cloudName>/<fluxControllerModel>
-                                fileName fCTimePath
-                                (
-                                    controllersTimePath/fCNames_[fC]
-                                );
-
-                                if (!isDir(fCTimePath))
-                                {
-                                    mkDir(fCTimePath);
-                                }
-
-                                const word& regionName =
-                                    fluxControllers_[fC]->regionName();
-
-                                // case/<timeDir>/uniform/controllers/
-                                //     <cloudName>/<fluxControllerModel>/
-                                //     <faceZoneName>
-                                fileName zoneTimePath(fCTimePath/regionName);
-
-                                if (!isDir(zoneTimePath))
-                                {
-                                    mkDir(zoneTimePath);
-                                }
-
-                                timePathNames[fC] = zoneTimePath;
-                            }
-                        }
-                    }
-                }
-            }
-
-            // write out data
-            forAll(fluxControllers_, fC)
-            {
-                fluxControllers_[fC]->output
-                (
-                    fCFixedPathNames_[fC],
-                    timePathNames[fC]
-                );
-            }
-        }
-
-        // Re-read dictionaries for modified properties (run-time selection)
-        {
-            stateControllersList_.clear();
-
-            stateControllersList_ = controllersDict_.lookup("stateControllers");
-
-            forAll(stateControllers_, sC)
-            {
-                const entry& controllersI = stateControllersList_[sC];
-                const dictionary& controllersIDict = controllersI.dict();
-
-                stateControllers_[sC]->updateProperties(controllersIDict);
-            }
-        }
-
-        {
-            fluxControllersList_.clear();
-
-            fluxControllersList_ = controllersDict_.lookup("fluxControllers");
-
-            forAll(fluxControllers_, fC)
-            {
-                const entry& controllersI = fluxControllersList_[fC];
-                const dictionary& controllersIDict = controllersI.dict();
-
-                fluxControllers_[fC]->updateProperties(controllersIDict);
-            }
-        }
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.H b/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.H
deleted file mode 100644
index 22e9ca7cca..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.H
+++ /dev/null
@@ -1,163 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-
-    controllers
-
-Description
-
-    Stores all the information for the controllers models defined within
-    the controllersDict, and selects & builds the models automatically.
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef controllers_H
-#define controllers_H
-
-#include "List.H"
-#include "IOdictionary.H"
-#include "autoPtr.H"
-#include "polyMesh.H"
-#include "timeData.H"
-#include "stateController.H"
-#include "fluxController.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                        Class controllers Declaration
-\*---------------------------------------------------------------------------*/
-
-class controllers
-{
-    // Private data
-
-        Time& time_;
-
-        //- The entire dictionary (containing multiple subDictionaries)
-        IOdictionary controllersDict_;
-
-        //- state controllers
-        PtrList<entry> stateControllersList_;
-        List<word> sCNames_;
-        List<label> sCIds_;
-        List<fileName> sCFixedPathNames_;
-        List< autoPtr<stateController> > stateControllers_;
-
-        //- flux controllers
-        PtrList<entry> fluxControllersList_;
-        List<word> fCNames_;
-        List<label> fCIds_;
-        List<fileName> fCFixedPathNames_;
-        List< autoPtr<fluxController> > fluxControllers_;
-
-
-public:
-
-    // Constructors
-
-        //- Null Constructor
-        controllers
-        (
-            const polyMesh& mesh
-        );
-
-        //- Constructor for with cloud
-        controllers
-        (
-            const polyMesh& mesh,
-            polyatomicCloud& cloud
-        );
-
-    //- Destructor
-    ~controllers();
-
-    // Member Functions
-
-        //- Initial configuration call this function after the polyatomicCloud
-        //  is completely initialised
-        void initialConfig();
-
-        //- this function is to be called at the beginning of the MD time-step.
-        //  since we have placed a non-referenced time-data class in the
-        //  state-controller class.
-        void updateTimeInfo();
-
-        //- control molecular state -- call this after the intermolecular force
-        //  calulation
-        void controlState();
-
-        //-
-        void controlVelocitiesI();
-
-        //-
-        void controlVelocitiesII();
-
-        //-
-        void controlPriorToForces();
-
-        //- calculate properties -- call this at the end of the MD time-step.
-        void calculateStateProps();
-
-        //- output -- call this function at the end of the MD time-step
-        void outputStateResults();
-
-        // Access
-
-            //-
-            inline List< autoPtr<stateController> >& stateControllers();
-
-            //-
-            inline const List< autoPtr<stateController> >&
-            stateControllers() const;
-
-            //-
-            inline List< autoPtr<fluxController> >& fluxControllers();
-
-            //-
-            inline const List< autoPtr<fluxController> >&
-            fluxControllers() const;
-
-            //-
-            inline const List<word>& stateControllersNames() const;
-
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "controllersI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllersI.H b/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllersI.H
deleted file mode 100644
index 66c10f0784..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllersI.H
+++ /dev/null
@@ -1,64 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-\*---------------------------------------------------------------------------*/
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-Foam::List<Foam::autoPtr<Foam::stateController> >&
-Foam::controllers::stateControllers()
-{
-    return stateControllers_;
-}
-
-
-const Foam::List<Foam::autoPtr<Foam::stateController> >&
-Foam::controllers::stateControllers() const
-{
-    return stateControllers_;
-}
-
-
-Foam::List<Foam::autoPtr<Foam::fluxController> >&
-Foam::controllers::fluxControllers()
-{
-    return fluxControllers_;
-}
-
-
-const Foam::List< autoPtr<fluxController> >&
-Foam::controllers::fluxControllers() const
-{
-    return fluxControllers_;
-}
-
-
-const Foam::List<Foam::word>& Foam::controllers::stateControllersNames() const
-{
-    return sCNames_;
-}
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.C b/src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.C
deleted file mode 100644
index 711b661758..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.C
+++ /dev/null
@@ -1,524 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "waterFluxController.H"
-#include "IFstream.H"
-#include "graph.H"
-#include "polyatomicCloud.H"
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-defineTypeNameAndDebug(waterFluxController, 0);
-
-defineRunTimeSelectionTable(waterFluxController, dictionary);
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-// Construct from components
-waterFluxController::waterFluxController
-(
-    Time& t,
-    polyatomicCloud& cloud,
-    const dictionary& dict
-)
-:
-    mesh_(refCast<const fvMesh>(cloud.mesh())),
-    cloud_(cloud),
-    rndGen_(clock::getTime()),
-    controllerDict_(dict.subDict("controllerProperties")),
-    timeDict_(controllerDict_.subDict("timeProperties")),
-    time_(t, timeDict_),
-    regionName_(controllerDict_.lookup("zoneName")),
-    regionId_(-1),
-    zoneSurfaceArea_(0.0),
-    internalFaces_(),
-    processorFaces_(),
-    control_(true),
-    readStateFromFile_(true),
-    singleValueController_(false),
-    density_(0.0),
-    velocity_(vector::zero),
-    temperature_(0.0),
-    pressure_(0.0),
-    strainRate_(tensor::zero),
-    tempGradient_(vector::zero),
-    fieldController_(false),
-    densities_(),
-    velocities_(),
-    temperatures_(),
-    pressures_(),
-    writeInTimeDir_(true),
-    writeInCase_(true)
-{
-    const faceZoneMesh& faceZones = mesh_.faceZones();
-    regionId_ = faceZones.findZoneID(regionName_);
-
-    if (regionId_ == -1)
-    {
-        FatalErrorIn("waterFluxController::waterFluxController()")
-        << "Cannot find region (faceZone): " << regionName_ << nl << "in: "
-        << t.time().system()/"controllersDict"
-        << exit(FatalError);
-    }
-
-    control_ = Switch(controllerDict_.lookup("controlSwitch"));
-    readStateFromFile_ = Switch(controllerDict_.lookup("readStateFromFile"));
-
-    setFacesInfo();
-}
-
-
-// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
-
-autoPtr<waterFluxController> waterFluxController::New
-(
-    Time& t,
-    polyatomicCloud& cloud,
-    const dictionary& dict
-)
-{
-    word waterFluxControllerName
-    (
-        dict.lookup("fluxControllerModel")
-    );
-
-    Info<< "Selecting fluxController "
-         << waterFluxControllerName << endl;
-
-    dictionaryConstructorTable::iterator cstrIter =
-        dictionaryConstructorTablePtr_->find(waterFluxControllerName);
-
-    if (cstrIter == dictionaryConstructorTablePtr_->end())
-    {
-        FatalError
-            << "waterFluxController::New(const dictionary&) : " << endl
-            << "    unknown waterFluxController type "
-            << waterFluxControllerName
-            << ", constructor not in hash table" << endl << endl
-            << "    Valid injector types are :" << endl;
-        Info<< dictionaryConstructorTablePtr_->toc() << abort(FatalError);
-    }
-
-    return autoPtr<waterFluxController>
-        (
-                cstrIter()(t, cloud, dict)
-        );
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-waterFluxController::~waterFluxController()
-{}
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-// void waterFluxController::updateTime()
-// {
-//     time_++;
-//
-//     const scalar& t = time_.time().timeOutputValue();
-//
-//     if ((t - initialTime_) < timePeriod_)
-//     {
-//         time_.controlTimeInterval().endTime() = false;
-// //         control_ = false;
-//     }
-//     else
-//     {
-// //         control_ = true;
-//     }
-// }
-
-
-void waterFluxController::setFacesInfo()
-{
-    const labelList& faces = controlZone();
-
-    if (Pstream::parRun())
-    {
-        DynamicList<label> processorFaces(0);
-
-        forAll(mesh_.boundaryMesh(), patchI)
-        {
-            const polyPatch& patch = mesh_.boundaryMesh()[patchI];
-
-            if (isA<processorPolyPatch>(patch))
-            {
-                for (label p = 0; p < patch.size(); p++)
-                {
-                    label patchFaceI = p + patch.start();
-                    label faceId = findIndex (faces, patchFaceI);
-
-                    if (faceId != -1)
-                    {
-                        processorFaces.append(patchFaceI);
-                    }
-                }
-            }
-        }
-
-        processorFaces.shrink();
-
-        processorFaces_.setSize(processorFaces.size(), -1);
-
-        forAll(processorFaces, f)
-        {
-            processorFaces_[f] = processorFaces[f];
-        }
-
-        label nInternalFaces = faces.size() - processorFaces.size();
-        internalFaces_.setSize(nInternalFaces, -1);
-
-        label counter = 0;
-
-        forAll(faces, f)
-        {
-            const label& faceI = faces[f];
-
-            if (findIndex(processorFaces, faceI) == -1)
-            {
-                internalFaces_[counter] = faceI;
-                counter++;
-            }
-        }
-
-//         Pout << "processorFaces: " << processorFaces_ << endl;
-//         Pout << "internalFaces: " << internalFaces_ << endl;
-
-        forAll(internalFaces_, f)
-        {
-            const label& faceI = internalFaces_[f];
-            zoneSurfaceArea_ += mag(mesh_.faceAreas()[faceI]);
-        }
-
-        // faces on a zone located on a processor cut belong to both processors
-        // (hence the 0.5)
-
-        forAll(processorFaces_, f)
-        {
-            const label& faceI = processorFaces_[f];
-            zoneSurfaceArea_ += 0.5*mag(mesh_.faceAreas()[faceI]);
-        }
-
-
-        if (Pstream::parRun())
-        {
-            for (int p = 0; p < Pstream::nProcs(); p++)
-            {
-                if (p != Pstream::myProcNo())
-                {
-                    const int proc = p;
-                    {
-                        OPstream toNeighbour(Pstream::blocking, proc);
-                        toNeighbour << zoneSurfaceArea_;
-                    }
-                }
-            }
-
-            //- receiving
-            for (int p = 0; p < Pstream::nProcs(); p++)
-            {
-                if (p != Pstream::myProcNo())
-                {
-                    scalar zoneSurfaceAreaProc;
-
-                    const int proc = p;
-                    {
-                        IPstream fromNeighbour(Pstream::blocking, proc);
-                        fromNeighbour >> zoneSurfaceAreaProc;
-                    }
-
-                    zoneSurfaceArea_ += zoneSurfaceAreaProc;
-                }
-            }
-        }
-    }
-    else
-    {
-        forAll(faces, f)
-        {
-            const label& faceI = faces[f];
-
-            zoneSurfaceArea_ += mag(mesh_.faceAreas()[faceI]);
-        }
-    }
-}
-
-
-
-void waterFluxController::updateTime()
-{
-    time_++;
-
-//     const scalar& t = time_.time().timeOutputValue();
-//
-//     if ((t - initialTime_) < timePeriod_)
-//     {
-//         time_.controlTimeInterval().endTime() = false;
-// //         control_ = false;
-//     }
-//     else
-//     {
-// //         control_ = true;
-//     }
-}
-
-void waterFluxController::updateFluxControllerProperties
-(
-    const dictionary& newDict
-)
-{
-    controllerDict_ = newDict.subDict("controllerProperties");
-
-    //- you can reset the controlling zone from here. This essentially
-    //  means that the coupling zone can infact move arbitrarily. To make
-    //  this happen we probably need to devise a technique for automatically
-    //  changing the cellZone else where, and then calling this function to
-    //  reset the controlling zone in which the controller operates in.
-
-    if (controllerDict_.found("controlSwitch"))
-    {
-        control_ = Switch(controllerDict_.lookup("controlSwitch"));
-    }
-
-    if (controllerDict_.found("readStateFromFile"))
-    {
-        readStateFromFile_ = Switch
-        (
-            controllerDict_.lookup("readStateFromFile")
-        );
-    }
-}
-
-const labelList& waterFluxController::controlZone() const
-{
-    return mesh_.faceZones()[regionId_];
-}
-
-label waterFluxController::isFaceOnControlZone(const label& faceI)
-{
-    const label f = findIndex(controlZone(), faceI);
-
-    return f;
-}
-
-const word& waterFluxController::regionName() const
-{
-    return regionName_;
-}
-
-const scalar& waterFluxController::density() const
-{
-    return density_;
-}
-
-scalar& waterFluxController::density()
-{
-    return density_;
-}
-
-const vector& waterFluxController::velocity() const
-{
-    return velocity_;
-}
-
-vector& waterFluxController::velocity()
-{
-    return velocity_;
-}
-
-const scalar& waterFluxController::temperature() const
-{
-    return temperature_;
-}
-
-scalar& waterFluxController::temperature()
-{
-    return temperature_;
-}
-
-const scalar& waterFluxController::pressure() const
-{
-    return pressure_;
-}
-
-scalar& waterFluxController::pressure()
-{
-    return pressure_;
-}
-
-const tensor& waterFluxController::strainRate() const
-{
-    return strainRate_;
-}
-
-tensor& waterFluxController::strainRate()
-{
-    return strainRate_;
-}
-
-const vector& waterFluxController::tempGradient() const
-{
-    return tempGradient_;
-}
-
-vector& waterFluxController::tempGradient()
-{
-    return tempGradient_;
-}
-
-
-const scalarField& waterFluxController::densityField() const
-{
-    return densities_;
-}
-
-scalarField& waterFluxController::densityField()
-{
-    return densities_;
-}
-
-const vectorField& waterFluxController::velocityField() const
-{
-    return velocities_;
-}
-vectorField& waterFluxController::velocityField()
-{
-    return velocities_;
-}
-
-const scalarField& waterFluxController::temperatureField() const
-{
-    return temperatures_;
-}
-
-scalarField& waterFluxController::temperatureField()
-{
-    return temperatures_;
-}
-
-const scalarField& waterFluxController::pressureField() const
-{
-    return pressures_;
-}
-
-scalarField& waterFluxController::pressureField()
-{
-    return pressures_;
-}
-
-
-const bool& waterFluxController::singleValueController() const
-{
-    return singleValueController_;
-}
-
-bool& waterFluxController::singleValueController()
-{
-    return singleValueController_;
-}
-
-const bool& waterFluxController::fieldController() const
-{
-    return fieldController_;
-}
-
-bool& waterFluxController::fieldController()
-{
-    return fieldController_;
-}
-
-
-const bool& waterFluxController::writeInTimeDir() const
-{
-    return writeInTimeDir_;
-}
-
-const bool& waterFluxController::writeInCase() const
-{
-    return writeInCase_;
-}
-
-
-// const scalar waterFluxController::avReqDensity() const
-// {
-//     scalar totalDensity = 0.0;
-//
-//     forAll(densities_, c)
-//     {
-//         totalDensity += densities_[c];
-//     }
-//
-//     if (cells_.size() > 0)
-//     {
-//         totalDensity /= scalar(cells_.size());
-//     }
-//
-//     return totalDensity;
-// }
-//
-// const vector waterFluxController::avReqVelocity() const
-// {
-//     vector totalVel = vector::zero;
-//
-//     forAll(velocities_, c)
-//     {
-//         totalVel += velocities_[c];
-//     }
-//
-//     if (cells_.size() > 0)
-//     {
-//         totalVel /= scalar(cells_.size());
-//     }
-//
-//     return totalVel;
-// }
-//
-// const scalar waterFluxController::avReqTemperature() const
-// {
-//     scalar totalTemp = 0.0;
-//
-//     forAll(densities_, c)
-//     {
-//         totalTemp += temperatures_[c];
-//     }
-//
-//     if (cells_.size() > 0)
-//     {
-//         totalTemp /= scalar(cells_.size());
-//     }
-//
-//     return totalTemp;
-// }
-
-
-
-} // End namespace Foam
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.H b/src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.H
deleted file mode 100644
index 9b8a375277..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.H
+++ /dev/null
@@ -1,272 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    waterFluxController
-
-Description
-
-SourceFiles
-    waterFluxControllerI.H
-    waterFluxController.C
-    waterFluxControllerIO.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef waterFluxController_H
-#define waterFluxController_H
-
-#include "IOdictionary.H"
-#include "Time.H"
-#include "autoPtr.H"
-#include "runTimeSelectionTables.H"
-#include "vector.H"
-#include "volFields.H"
-#include "Random.H"
-#include "polyatomic.H"
-#include "timeData.H"
-#include "writeTimeData.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class waterFluxController Declaration
-\*---------------------------------------------------------------------------*/
-
-class waterFluxController
-{
-
-protected:
-
-    // Protected data
-
-//         Time& time_;
-
-        const fvMesh& mesh_;
-
-        polyatomicCloud& cloud_;
-
-        Random rndGen_;
-
-        //- subDictionary containing the properties
-        dictionary controllerDict_;
-
-
-        dictionary timeDict_;
-
-        timeData time_;
-
-        //- name of face zone
-        word regionName_;
-        label regionId_;
-//         labelList faces_;
-
-        scalar zoneSurfaceArea_;
-
-        labelList internalFaces_;
-        labelList processorFaces_;
-
-        bool control_;
-        bool readStateFromFile_;
-
-        //- set all the properties below from model if required
-
-            bool singleValueController_;
-
-            // target values
-            scalar density_;
-            vector velocity_;
-            scalar temperature_;
-            scalar pressure_;
-
-            tensor strainRate_;
-            vector tempGradient_;
-
-            bool fieldController_;
-
-            //- targeted fields
-            scalarField densities_;
-            vectorField velocities_;
-            scalarField temperatures_;
-            scalarField pressures_;
-
-            bool writeInTimeDir_;
-            bool writeInCase_;
-
-
-
-
-    // Private Member Functions
-
-        void setFacesInfo();
-
-
-public:
-
-    //- Runtime type information
-    TypeName("waterFluxController");
-
-    // Declare runtime constructor selection table
-        declareRunTimeSelectionTable
-        (
-            autoPtr,
-            waterFluxController,
-            dictionary,
-            (
-                Time& t,
-                polyatomicCloud& cloud,
-                const dictionary& dict
-            ),
-            (t, cloud, dict)
-        );
-
-    // Constructors
-
-        //- Construct from components
-        waterFluxController
-        (
-                        Time& t,
-            polyatomicCloud& cloud,
-            const dictionary& dict
-        );
-
-
-    // Selectors
-
-        static autoPtr<waterFluxController> New
-        (
-                        Time& t,
-            polyatomicCloud& cloud,
-            const dictionary& dict
-        );
-
-
-    // Destructor
-
-        virtual ~waterFluxController();
-
-
-    // Member Functions
-
-        void updateTime();
-
-        //- create an initial configuration
-        virtual void initialConfiguration() = 0;
-
-
-        //- calculate any required properties
-        virtual void calculateProperties() = 0;
-
-
-        //- control the polyatomic from the tracking function
-        virtual void controlMol
-        (
-            polyatomic& mol,
-            polyatomic::trackData& td
-        ) = 0;
-
-        //- output data
-        virtual void output
-        (
-            const fileName& fixedPathName,
-            const fileName& timePath
-        ) = 0;
-
-
-        //- E. update properties from a modified dictionary
-        virtual void updateProperties(const dictionary&) = 0;
-
-
-        void updateFluxControllerProperties(const dictionary&);
-
-        // Access
-
-            //- return the control zone cells
-            const labelList& controlZone() const;
-
-
-            label isFaceOnControlZone(const label& faceI);
-
-            //- return the control zone name
-            const word& regionName() const;
-
-           //- return the targeted values
-            const scalar& density() const;
-            scalar& density();
-
-            const vector& velocity() const;
-            vector& velocity();
-
-            const scalar& temperature() const;
-            scalar& temperature();
-
-            const scalar& pressure() const;
-            scalar& pressure();
-
-            const tensor& strainRate() const;
-            tensor& strainRate();
-
-            const vector& tempGradient() const;
-            vector& tempGradient();
-
-            //- return the targeted fields
-            const scalarField& densityField() const;
-            scalarField& densityField();
-
-            const vectorField& velocityField() const;
-            vectorField& velocityField();
-
-            const scalarField& temperatureField() const;
-            scalarField& temperatureField();
-
-            const scalarField& pressureField() const;
-            scalarField& pressureField();
-
-
-            const bool& singleValueController() const;
-            bool& singleValueController();
-
-            const bool& fieldController() const;
-            bool& fieldController();
-
-            const bool& writeInTimeDir() const;
-            const bool& writeInCase() const;
-
-//             const scalar avReqDensity() const;
-//             const vector avReqVelocity() const;
-//             const scalar avReqTemperature() const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.C b/src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.C
deleted file mode 100644
index b98aa126a3..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.C
+++ /dev/null
@@ -1,490 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "stateController.H"
-#include "IFstream.H"
-#include "polyatomicCloud.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-defineTypeNameAndDebug(stateController, 0);
-
-defineRunTimeSelectionTable(stateController, dictionary);
-
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::stateController::stateController
-(
-    polyatomicCloud& cloud,
-    const dictionary& dict
-)
-:
-    mesh_(refCast<const fvMesh>(cloud.mesh())),
-    cloud_(cloud),
-    rndGen_(clock::getTime()),
-    controllerDict_(dict.subDict("controllerProperties")),
-    timeDict_(controllerDict_.subDict("timeProperties")),
-    time_(mesh_.time(), timeDict_),
-    timePeriod_(readScalar(timeDict_.lookup("initialTimePeriod"))), //temp
-    initialTime_(time_.time().startTime().value()),
-    regionName_(controllerDict_.lookup("zoneName")),
-    regionId_(-1),
-    control_(true),
-    readStateFromFile_(true),
-    singleValueController_(false),
-    density_(0.0),
-    velocity_(vector::zero),
-    temperature_(0.0),
-    pressure_(0.0),
-    strainRate_(tensor::zero),
-    tempGradient_(vector::zero),
-    fieldController_(false),
-    densities_(),
-    velocities_(),
-    temperatures_(),
-    pressures_(),
-    writeInTimeDir_(true),
-    writeInCase_(true)
-{
-    const cellZoneMesh& cellZones = mesh_.cellZones();
-
-    regionId_ = cellZones.findZoneID(regionName_);
-
-    if (regionId_ == -1)
-    {
-        FatalErrorIn("stateController::stateController()")
-            << "Cannot find region: " << regionName_ << nl << "in: "
-            << time_.time().system()/"controllersDict"
-            << exit(FatalError);
-    }
-
-    control_ = Switch(controllerDict_.lookup("controlSwitch"));
-
-    readStateFromFile_ = Switch(controllerDict_.lookup("readStateFromFile"));
-
-    const scalar& avTimeInterval = time_.averageTimeInterval().deltaT();
-
-    if ((timePeriod_ < avTimeInterval) && (timePeriod_ > 0.0))
-    {
-        timePeriod_ = avTimeInterval;
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
-
-Foam::autoPtr<Foam::stateController> Foam::stateController::New
-(
-    polyatomicCloud& cloud,
-    const dictionary& dict
-)
-{
-    word stateControllerName
-    (
-        dict.lookup("stateControllerModel")
-    );
-
-    Info<< "Selecting stateController "
-         << stateControllerName << endl;
-
-    dictionaryConstructorTable::iterator cstrIter =
-        dictionaryConstructorTablePtr_->find(stateControllerName);
-
-    if (cstrIter == dictionaryConstructorTablePtr_->end())
-    {
-        FatalError
-            << "stateController::New(const dictionary&) : " << endl
-            << "    unknown stateController type "
-            << stateControllerName
-            << ", constructor not in hash table" << endl << endl
-            << "    Valid types are :" << endl;
-        Info<< dictionaryConstructorTablePtr_->toc() << abort(FatalError);
-    }
-
-    return autoPtr<stateController>
-    (
-        cstrIter()(cloud, dict)
-    );
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::stateController::~stateController()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::stateController::updateTime()
-{
-    time_++;
-
-    const scalar& t = time_.time().timeOutputValue();
-
-    if ((t - initialTime_) < timePeriod_)
-    {
-        time_.controlTimeInterval().endTime() = false;
-    }
-}
-
-
-void Foam::stateController::updateStateControllerProperties
-(
-    const dictionary& newDict
-)
-{
-    controllerDict_ = newDict.subDict("controllerProperties");
-
-    if (controllerDict_.found("controlSwitch"))
-    {
-        control_ = Switch(controllerDict_.lookup("controlSwitch"));
-    }
-
-    if (controllerDict_.found("readStateFromFile"))
-    {
-        readStateFromFile_ = Switch
-        (
-            controllerDict_.lookup("readStateFromFile")
-        );
-    }
-
-    timeDict_ = controllerDict_.subDict("timeProperties");
-
-    if (timeDict_.found("resetAtOutput"))
-    {
-        time_.resetFieldsAtOutput() = Switch(timeDict_.lookup("resetAtOutput"));
-    }
-}
-
-
-const Foam::labelList& Foam::stateController::controlZone() const
-{
-    return mesh_.cellZones()[regionId_];
-}
-
-const Foam::word& Foam::stateController::regionName() const
-{
-    return regionName_;
-}
-
-
-Foam::scalar Foam::stateController::density() const
-{
-    return density_;
-}
-
-
-Foam::scalar& Foam::stateController::density()
-{
-    return density_;
-}
-
-
-const Foam::vector& Foam::stateController::velocity() const
-{
-    return velocity_;
-}
-
-
-Foam::vector& Foam::stateController::velocity()
-{
-    return velocity_;
-}
-
-
-Foam::scalar Foam::stateController::temperature() const
-{
-    return temperature_;
-}
-
-
-Foam::scalar& Foam::stateController::temperature()
-{
-    return temperature_;
-}
-
-const Foam::scalar& Foam::stateController::pressure() const
-{
-    return pressure_;
-}
-
-
-Foam::scalar& Foam::stateController::pressure()
-{
-    return pressure_;
-}
-
-
-const Foam::tensor& Foam::stateController::strainRate() const
-{
-    return strainRate_;
-}
-
-
-Foam::tensor& Foam::stateController::strainRate()
-{
-    return strainRate_;
-}
-
-
-const Foam::vector& Foam::stateController::tempGradient() const
-{
-    return tempGradient_;
-}
-
-
-Foam::vector& Foam::stateController::tempGradient()
-{
-    return tempGradient_;
-}
-
-
-const Foam::scalarField& Foam::stateController::densityField() const
-{
-    return densities_;
-}
-
-Foam::scalarField& Foam::stateController::densityField()
-{
-    return densities_;
-}
-
-
-const Foam::vectorField& Foam::stateController::velocityField() const
-{
-    return velocities_;
-}
-
-
-Foam::vectorField& Foam::stateController::velocityField()
-{
-    return velocities_;
-}
-
-
-const Foam::scalarField& Foam::stateController::temperatureField() const
-{
-    return temperatures_;
-}
-
-
-Foam::scalarField& Foam::stateController::temperatureField()
-{
-    return temperatures_;
-}
-
-
-const Foam::scalarField& Foam::stateController::pressureField() const
-{
-    return pressures_;
-}
-
-
-Foam::scalarField& Foam::stateController::pressureField()
-{
-    return pressures_;
-}
-
-
-bool Foam::stateController::singleValueController() const
-{
-    return singleValueController_;
-}
-
-
-bool& Foam::stateController::singleValueController()
-{
-    return singleValueController_;
-}
-
-
-bool Foam::stateController::fieldController() const
-{
-    return fieldController_;
-}
-
-
-bool& Foam::stateController::fieldController()
-{
-    return fieldController_;
-}
-
-
-bool Foam::stateController::writeInTimeDir() const
-{
-    return writeInTimeDir_;
-}
-
-
-bool Foam::stateController::writeInCase() const
-{
-    return writeInCase_;
-}
-
-
-Foam::scalar Foam::stateController::avReqDensity()
-{
-    scalar totalDensity = 0.0;
-
-    if (singleValueController_)
-    {
-        totalDensity = density_;
-    }
-    else if (fieldController_)
-    {
-        label controlCells = controlZone().size();
-
-        forAll(densities_, c)
-        {
-            totalDensity += densities_[c];
-        }
-
-        if (Pstream::parRun())
-        {
-            reduce(totalDensity, sumOp<scalar>());
-
-            reduce(controlCells, sumOp<label>());
-        }
-
-        if (controlCells > 0)
-        {
-            totalDensity /= scalar(controlCells);
-        }
-    }
-
-    return totalDensity;
-}
-
-
-Foam::vector Foam::stateController::avReqVelocity()
-{
-    vector totalVel = vector::zero;
-
-    if (singleValueController_)
-    {
-        totalVel = velocity_;
-    }
-    else if (fieldController_)
-    {
-        label controlCells = controlZone().size();
-
-        forAll(velocities_, c)
-        {
-            totalVel += velocities_[c];
-        }
-
-        if (Pstream::parRun())
-        {
-            reduce(totalVel, sumOp<vector>());
-
-            reduce(controlCells, sumOp<label>());
-        }
-
-        if (controlCells > 0)
-        {
-            totalVel /= scalar(controlCells);
-        }
-    }
-
-    return totalVel;
-}
-
-
-Foam::scalar Foam::stateController::avReqTemperature()
-{
-    scalar totalTemp = 0.0;
-
-    if (singleValueController_)
-    {
-        totalTemp = temperature_;
-    }
-    else if (fieldController_)
-    {
-        label controlCells = controlZone().size();
-
-        forAll(temperatures_, c)
-        {
-            totalTemp += temperatures_[c];
-        }
-
-        if (Pstream::parRun())
-        {
-            reduce(totalTemp, sumOp<scalar>());
-
-            reduce(controlCells, sumOp<label>());
-        }
-
-        if (controlCells > 0)
-        {
-            totalTemp /= scalar(controlCells);
-        }
-    }
-
-    return totalTemp;
-}
-
-
-Foam::scalar Foam::stateController::avReqPressure()
-{
-    scalar totalPressure = 0.0;
-
-    if (singleValueController_)
-    {
-        totalPressure = pressure_;
-    }
-    else if (fieldController_)
-    {
-        label controlCells = controlZone().size();
-
-        forAll(pressures_, c)
-        {
-            totalPressure += pressures_[c];
-        }
-
-        if (Pstream::parRun())
-        {
-            reduce(totalPressure, sumOp<scalar>());
-
-            reduce(controlCells, sumOp<label>());
-        }
-
-        if (controlCells > 0)
-        {
-            totalPressure /= scalar(controlCells);
-        }
-    }
-
-    return totalPressure;
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.H b/src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.H
deleted file mode 100644
index 3e142c3962..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.H
+++ /dev/null
@@ -1,270 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    stateController
-
-Description
-
-    Basic/abstract class of a state controller
-
-SourceFiles
-    stateControllerI.H
-    stateController.C
-    stateControllerIO.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef stateController_H
-#define stateController_H
-
-#include "IOdictionary.H"
-#include "autoPtr.H"
-#include "runTimeSelectionTables.H"
-#include "vector.H"
-#include "volFields.H"
-#include "Random.H"
-#include "polyatomic.H"
-#include "timeData.H"
-#include "writeTimeData.H"
-#include "selectIds.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                      Class stateController Declaration
-\*---------------------------------------------------------------------------*/
-
-class stateController
-{
-
-protected:
-
-    // Protected data
-
-        //-
-        const fvMesh& mesh_;
-
-        //-
-        polyatomicCloud& cloud_;
-
-        //-
-        Random rndGen_;
-
-        //- subDictionary containing the properties
-        dictionary controllerDict_;
-
-        //-
-        dictionary timeDict_;
-
-        //-
-        timeData time_;
-
-        //-
-        scalar timePeriod_;
-
-        //-
-        scalar initialTime_;
-
-        //- name of control zone
-        word regionName_;
-
-        //-
-        label regionId_;
-
-        //-
-        bool control_;
-
-        //-
-        bool readStateFromFile_;
-
-        //- set all the properties below from model if required
-
-            //-
-            bool singleValueController_;
-
-            //- target values
-            scalar density_;
-            vector velocity_;
-            scalar temperature_;
-            scalar pressure_;
-            tensor strainRate_;
-            vector tempGradient_;
-
-            //- set this in model
-            bool fieldController_;
-
-            //- targeted fields
-            scalarField densities_;
-            vectorField velocities_;
-            scalarField temperatures_;
-            scalarField pressures_;
-
-            // set these in model
-            bool writeInTimeDir_;
-            bool writeInCase_;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("stateController");
-
-    //- Declare runtime constructor selection table
-        declareRunTimeSelectionTable
-        (
-            autoPtr,
-            stateController,
-            dictionary,
-            (
-                polyatomicCloud& cloud,
-                const dictionary& dict
-            ),
-            (t, cloud, dict)
-        );
-
-    // Constructors
-
-        //- Construct from components
-        stateController
-        (
-            polyatomicCloud& cloud,
-            const dictionary& dict
-        );
-
-
-    // Selectors
-
-        static autoPtr<stateController> New
-        (
-            polyatomicCloud& cloud,
-            const dictionary& dict
-        );
-
-
-    // Destructor
-
-        virtual ~stateController();
-
-
-    // Member Functions
-
-
-        void updateTime();
-
-        //- create an initial configuration
-        virtual void initialConfiguration() = 0;
-
-        //- calculate any required properties
-        virtual void calculateProperties() = 0;
-
-        //- control molecules at different stages of the integration time-step
-        virtual void controlMolsBeg() = 0;
-
-        virtual void controlBeforeForces() = 0;
-
-        virtual void controlMols() = 0;
-
-        virtual void controlMolsEnd() = 0;
-
-
-        //- output data
-        virtual void output
-        (
-            const fileName& fixedPathName,
-            const fileName& timePath
-        ) = 0;
-
-        //- update properties from a modified dictionary
-        virtual void updateProperties(const dictionary&) = 0;
-
-        void updateStateControllerProperties(const dictionary&);
-
-        // Access
-
-            //- return the control zone cells
-            const labelList& controlZone() const;
-
-            //- return the control zone name
-            const word& regionName() const;
-
-            //- return the targeted fields
-            scalar density() const;
-            scalar& density();
-
-            const vector& velocity() const;
-            vector& velocity();
-
-            scalar temperature() const;
-            scalar& temperature();
-
-            scalar pressure() const;
-            scalar& pressure();
-
-            const tensor& strainRate() const;
-            tensor& strainRate();
-
-            const vector& tempGradient() const;
-            vector& tempGradient();
-
-            //- return the targeted fields
-            const scalarField& densityField() const;
-            scalarField& densityField();
-
-            const vectorField& velocityField() const;
-            vectorField& velocityField();
-
-            const scalarField& temperatureField() const;
-            scalarField& temperatureField();
-
-            const scalarField& pressureField() const;
-            scalarField& pressureField();
-
-            bool singleValueController() const;
-            bool& singleValueController();
-
-            bool fieldController() const;
-            bool& fieldController();
-
-            bool writeInTimeDir() const;
-            bool writeInCase() const;
-
-            scalar avReqDensity();
-            vector avReqVelocity();
-            scalar avReqTemperature();
-            scalar avReqPressure();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.C b/src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.C
deleted file mode 100644
index 08ea263dd6..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.C
+++ /dev/null
@@ -1,114 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "constPropSite.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::constPropSite::constPropSite()
-:
-    siteReferencePosition_(vector::zero),
-    siteMass_(0.0),
-    siteCharge_(0.0),
-    siteId_(0),
-    name_(),
-    pairPotentialSite_(false),
-    electrostaticSite_(false)
-{}
-
-
-Foam::constPropSite::constPropSite
-(
-    const vector& siteReferencePosition,
-    const scalar& siteMass,
-    const scalar& siteCharge,
-    const label& siteId,
-    const word& name,
-    const bool& pairPotentialSite,
-    const bool& electrostaticSite
-)
-:
-    siteReferencePosition_(siteReferencePosition),
-    siteMass_(siteMass),
-    siteCharge_(siteCharge),
-    siteId_(siteId),
-    name_(name),
-    pairPotentialSite_(pairPotentialSite),
-    electrostaticSite_(electrostaticSite)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::constPropSite::~constPropSite()
-{}
-
-
-// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
-
-Foam::Istream& Foam::operator>>(Istream& is, constPropSite& cPS)
-{
-    is  >> cPS.siteReferencePosition_
-        >> cPS.siteMass_
-        >> cPS.siteCharge_
-        >> cPS.siteId_
-        >> cPS.name_
-        >> cPS.pairPotentialSite_
-        >> cPS.electrostaticSite_;
-
-    // Check state of Istream
-    is.check
-    (
-        "Foam::Istream& Foam::operator>>"
-        "(Foam::Istream&, Foam::constPropSite&)"
-    );
-
-    return is;
-}
-
-
-Foam::Ostream& Foam::operator<<(Ostream& os, const constPropSite& cPS)
-{
-
-    os  << token::SPACE << cPS.siteReferencePosition()
-        << token::SPACE << cPS.siteMass()
-        << token::SPACE << cPS.siteCharge()
-        << token::SPACE << cPS.siteId()
-        << token::SPACE << cPS.name()
-        << token::SPACE << cPS.pairPotentialSite()
-        << token::SPACE << cPS.electrostaticSite();
-
-    // Check state of Ostream
-    os.check
-    (
-        "Foam::Ostream& Foam::operator<<(Foam::Ostream&, "
-        "const Foam::constPropSite&)"
-    );
-
-    return os;
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.H b/src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.H
deleted file mode 100644
index 910c0ccaef..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.H
+++ /dev/null
@@ -1,185 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-
-
-Description
-
-
-SourceFiles
-    constPropSiteI.H
-    constPropSite.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef constPropSite_H
-#define constPropSite_H
-
-#include "vector.H"
-#include "IOstreams.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-class Istream;
-class Ostream;
-
-// Forward declaration of friend functions and operators
-class constPropSite;
-Istream& operator>>(Istream&, constPropSite&);
-Ostream& operator<<(Ostream&, const constPropSite&);
-
-
-/*---------------------------------------------------------------------------*\
-                         Class constPropSite Declaration
-\*---------------------------------------------------------------------------*/
-
-class constPropSite
-{
-    // Private data
-
-        //-
-        vector siteReferencePosition_;
-
-        //-
-        scalar siteMass_;
-
-        //-
-        scalar siteCharge_;
-
-        //-
-        label siteId_;
-
-        //-
-        word name_;
-
-        //-
-        bool pairPotentialSite_;
-
-        //-
-        bool electrostaticSite_;
-
-
-public:
-
-    // Constructors
-
-        //- Construct null
-        constPropSite();
-
-        //- Construct from components
-        constPropSite
-        (
-            const vector& siteReferencePosition,
-            const scalar& siteMass,
-            const scalar& siteCharge,
-            const label& siteId,
-            const word& name,
-            const bool& pairPotentialSite,
-            const bool& electrostaticSite
-        );
-
-
-    // Selectors
-
-    // Destructor
-
-        ~constPropSite();
-
-    // Member Functions
-
-        // Access
-
-            //-
-            inline const vector& siteReferencePosition() const;
-
-            //-
-            inline vector& siteReferencePosition();
-
-            //-
-            inline const scalar& siteMass() const;
-
-            //-
-            inline scalar& siteMass();
-
-            //-
-            inline const scalar& siteCharge() const;
-
-            //-
-            inline scalar& siteCharge();
-
-            //-
-            inline const label& siteId() const;
-
-            //-
-            inline label& siteId();
-
-            //-
-            inline const word& name() const;
-
-            //-
-            inline word& name();
-
-            //-
-            inline const bool& pairPotentialSite() const;
-
-            //-
-            inline bool& pairPotentialSite();
-
-            //-
-            inline const bool& electrostaticSite() const;
-
-            //-
-            inline bool& electrostaticSite();
-
-
-    // Member Operators
-
-        inline bool operator==(const constPropSite&) const;
-        inline bool operator!=(const constPropSite&) const;
-
-
-    // IOstream Operators
-
-        friend Istream& operator>>(Istream&, constPropSite&);
-        friend Ostream& operator<<(Ostream&, const constPropSite&);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "constPropSiteI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSiteI.H b/src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSiteI.H
deleted file mode 100644
index 9c648c94c8..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSiteI.H
+++ /dev/null
@@ -1,139 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-inline const Foam::vector& Foam::constPropSite::siteReferencePosition() const
-{
-    return siteReferencePosition_;
-}
-
-
-inline Foam::vector& Foam::constPropSite::siteReferencePosition()
-{
-    return siteReferencePosition_;
-}
-
-
-inline const Foam::scalar& Foam::constPropSite::siteMass() const
-{
-    return siteMass_;
-}
-
-
-inline Foam::scalar& Foam::constPropSite::siteMass()
-{
-    return siteMass_;
-}
-
-
-inline const Foam::scalar& Foam::constPropSite::siteCharge() const
-{
-    return siteCharge_;
-}
-
-
-inline Foam::scalar& Foam::constPropSite::siteCharge()
-{
-    return siteCharge_;
-}
-
-
-inline const Foam::label& Foam::constPropSite::siteId() const
-{
-    return siteId_;
-}
-
-
-inline Foam::label& Foam::constPropSite::siteId()
-{
-    return siteId_;
-}
-
-
-inline const Foam::word& Foam::constPropSite::name() const
-{
-    return name_;
-}
-
-
-inline Foam::word& Foam::constPropSite::name()
-{
-    return name_;
-}
-
-
-inline const bool& Foam::constPropSite::pairPotentialSite() const
-{
-    return pairPotentialSite_;
-}
-
-
-inline bool& Foam::constPropSite::pairPotentialSite()
-{
-    return pairPotentialSite_;
-}
-
-
-inline const bool& Foam::constPropSite::electrostaticSite() const
-{
-    return electrostaticSite_;
-}
-
-
-inline bool& Foam::constPropSite::electrostaticSite()
-{
-    return electrostaticSite_;
-}
-
-
-// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
-
-bool Foam::constPropSite::operator==
-(
-    const constPropSite& rhs
-) const
-{
-    return
-       siteReferencePosition_ == rhs.siteReferencePosition_
-    && siteMass_ == rhs.siteMass_
-    && siteCharge_ == rhs.siteCharge_
-    && siteId_ == rhs.siteId_
-    && name_ == rhs.name_
-    && pairPotentialSite_ == rhs.pairPotentialSite_
-    && electrostaticSite_ == rhs.electrostaticSite_;
-}
-
-
-bool Foam::constPropSite::operator!=
-(
-    const constPropSite& rhs
-) const
-{
-    return !(*this == rhs);
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomic.C b/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomic.C
deleted file mode 100644
index f5b0316487..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomic.C
+++ /dev/null
@@ -1,199 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*----------------------------------------------------------------------------*/
-
-#include "monoatomic.H"
-#include "Random.H"
-#include "Time.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(monoatomic, 0);
-    defineTemplateTypeNameAndDebug(Cloud<monoatomic>, 0);
-}
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::monoatomic::move
-(
-    monoatomic::trackingData& td,
-    const scalar trackTime
-)
-{
-    td.switchProcessor = false;
-    td.keepParticle = true;
-
-    if (special_ != SPECIAL_FROZEN)
-    {
-        return td.keepParticle;
-    }
-
-    const constantProperties& constProps(td.cloud().constProps(id_));
-
-    if (td.part() == trackingData::tpFirstVelocityHalfStep)
-    {
-        // First leapfrog velocity adjust part, required before tracking+force
-        // part
-
-        v_ += 0.5*trackTime*a_;
-    }
-    else if (td.part() == trackingData::tpLinearTrack)
-    {
-        // Leapfrog tracking part
-
-        scalar tEnd = (1.0 - stepFraction())*trackTime;
-        scalar dtMax = tEnd;
-
-        while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
-        {
-            // set the lagrangian time-step
-            scalar dt = min(dtMax, tEnd);
-
-            dt *= trackToFace(position() + dt*v_, td);
-
-            tEnd -= dt;
-            stepFraction() = 1.0 - tEnd/trackTime;
-        }
-
-        setSitePositions(constProps);
-    }
-    else if (td.part() == trackingData::tpSecondVelocityHalfStep)
-    {
-        // Second leapfrog velocity adjust part, required after tracking+force
-        // part
-
-        a_ = siteForces_[0]/constProps.mass();
-
-        v_ += 0.5*trackTime*a_;
-    }
-    else if (td.part() != trackingData::tpRotationalTrack)
-    {
-        FatalErrorIn("monoatomic::move(trackingData&, const scalar)") << nl
-            << td.part() << " is an invalid part of the integration method."
-            << abort(FatalError);
-    }
-
-    return td.keepParticle;
-}
-
-
-void Foam::monoatomic::transformProperties(const tensor& T)
-{
-    particle::transformProperties(T);
-
-    v_ = transform(T, v_);
-
-    a_ = transform(T, a_);
-
-    rf_ = transform(T, rf_);
-
-    sitePositions_[0] = position_ + (T & (sitePositions_[0] - position_));
-
-    siteForces_[0] = T & siteForces_[0];
-}
-
-
-void Foam::monoatomic::transformProperties(const vector& separation)
-{
-    particle::transformProperties(separation);
-
-    if (special_ == SPECIAL_TETHERED)
-    {
-        specialPosition_ += separation;
-    }
-
-    sitePositions_[0] += separation;
-}
-
-
-void Foam::monoatomic::setSitePositions(const constantProperties& constProps)
-{
-    sitePositions_[0] = position_;
-}
-
-
-void Foam::monoatomic::setSiteSizes(label size)
-{
-    // Nothing required, size controlled internally
-}
-
-
-bool Foam::monoatomic::hitPatch
-(
-    const polyPatch&,
-    trackingData&,
-    const label,
-    const scalar,
-    const tetIndices&
-)
-{
-    return false;
-}
-
-
-void Foam::monoatomic::hitProcessorPatch
-(
-    const processorPolyPatch&,
-    trackingData& td
-)
-{
-    td.switchProcessor = true;
-}
-
-
-void Foam::monoatomic::hitWallPatch
-(
-    const wallPolyPatch& wpp,
-    trackingData& td,
-    const tetIndices& tetIs
-)
-{
-    // Use of the normal from tetIs is not required as
-    // hasWallImpactDistance is false.
-    vector nw = normal();
-    nw /= mag(nw);
-
-    scalar vn = v_ & nw;
-
-    // Specular reflection
-    if (vn > 0)
-    {
-        v_ -= 2*vn*nw;
-    }
-}
-
-
-void Foam::monoatomic::hitPatch
-(
-    const polyPatch&,
-    trackingData& td
-)
-{
-    td.keepParticle = false;
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomic.H b/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomic.H
deleted file mode 100644
index d1e999435b..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomic.H
+++ /dev/null
@@ -1,418 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::monoatomic
-
-Description
-    Foam::monoatomic
-
-SourceFiles
-    monoatomicI.H
-    monoatomic.C
-    monoatomicIO.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef monoatomic_H
-#define monoatomic_H
-
-#include "particle.H"
-#include "IOstream.H"
-#include "autoPtr.H"
-#include "diagTensor.H"
-#include "constPropSite.H"
-#include "MoleculeCloud.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class monoatomic Declaration
-\*---------------------------------------------------------------------------*/
-
-class monoatomic
-:
-    public particle
-{
-
-public:
-
-    // 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 monoatomic constant properties
-    class constantProperties
-    {
-        // Private data
-
-            //- Sites of mass, charge or interaction
-            FixedList<constPropSite, 1> sites_;
-
-            //- Which sites require electrostatic interactions
-            FixedList<label, 1> electrostaticSites_;
-
-            //- Which sites require pair interactions
-            FixedList<label, 1> pairPotSites_;
-
-            //-
-            scalar mass_;
-
-
-    public:
-
-        //-
-        inline constantProperties();
-
-        //- Construct from dictionary
-        inline constantProperties
-        (
-            const dictionary& dict,
-            const List<label>& siteIds
-        );
-
-        // Member functions
-
-        //-
-        inline const FixedList<constPropSite, 1>& sites() const;
-
-        //-
-        inline const FixedList<label, 1>& pairPotSites() const;
-
-        //-
-        inline const FixedList<label, 1>& electrostaticSites() const;
-
-        //-
-        inline label degreesOfFreedom() const;
-
-        //-
-        inline scalar mass() const;
-
-        //-
-        inline label nSites() const;
-    };
-
-
-    //- Class used to pass tracking data to the trackToFace function
-    class trackingData
-    :
-        public particle::TrackingData<MoleculeCloud<monoatomic> >
-    {
-    public:
-
-        enum trackPart
-        {
-            tpFirstVelocityHalfStep,
-            tpLinearTrack,
-            tpRotationalTrack,
-            tpSecondVelocityHalfStep,
-            tpAccess
-        };
-
-    private:
-
-        // label specifying which part of the integration algorithm is taking
-        label part_;
-
-
-    public:
-
-        // Constructors
-
-            trackingData(MoleculeCloud<monoatomic>& cloud, trackPart part)
-            :
-                particle::TrackingData<MoleculeCloud<monoatomic> >(cloud),
-                part_(part)
-            {}
-
-        // Member functions
-
-            inline label part() const
-            {
-                return part_;
-            }
-    };
-
-
-private:
-
-    // Private data
-
-          //- Linear velocity of monoatomic
-        vector v_;
-
-        //- Total linear acceleration of monoatomic
-        vector a_;
-
-            //-
-        vector specialPosition_;
-
-        //-
-        scalar potentialEnergy_;
-
-        // - r_ij f_ij, stress dyad
-        tensor rf_;
-
-//         // - r_ij outer product f_ij: virial contribution
-//         tensor rDotf_;
-
-        //-
-        label special_;
-
-        //-
-        label id_;
-
-        //-
-        FixedList<vector, 1> siteForces_;
-
-        //-
-        FixedList<vector, 1> sitePositions_;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("monoatomic");
-
-    friend class Cloud<monoatomic>;
-
-    // Constructors
-
-        //- Construct with macroscopic description
-        inline monoatomic
-        (
-            const polyMesh& mesh,
-            const vector& position,
-            const label cellI,
-            const label tetFaceI,
-            const label tetPtI,
-            const scalar temperature,
-            const vector& bulkVelocity,
-            const vector& specialPosition,
-            const constantProperties& constProps,
-            trackingData& td,
-            const label special,
-            const label id
-        );
-
-        //- Construct from all components
-        inline monoatomic
-        (
-            const polyMesh& mesh,
-            const vector& position,
-            const label cellI,
-            const label tetFaceI,
-            const label tetPtI,
-            const vector& v,
-            const vector& a,
-            const vector& specialPosition,
-            const constantProperties& constProps,
-            const label special,
-            const label id
-        );
-
-        //- Construct from Istream
-        monoatomic
-        (
-            const polyMesh& mesh,
-            Istream& is,
-            bool readFields = true
-        );
-
-        //- Construct and return a clone
-        autoPtr<particle> clone() const
-        {
-            return autoPtr<particle>(new monoatomic(*this));
-        }
-
-        //- Factory class to read-construct particles used for
-        //  parallel transfer
-        class iNew
-        {
-            const polyMesh& mesh_;
-
-        public:
-
-            iNew(const polyMesh& mesh)
-            :
-                mesh_(mesh)
-            {}
-
-            autoPtr<monoatomic> operator()(Istream& is) const
-            {
-                return autoPtr<monoatomic>(new monoatomic(mesh_, is, true));
-            }
-        };
-
-
-    // Member Functions
-
-        // Tracking
-
-            //-
-            bool move(trackingData&, const scalar trackTime);
-
-            //-
-            virtual void transformProperties(const tensor& T);
-
-            //-
-            virtual void transformProperties(const vector& separation);
-
-            //-
-            void setSitePositions(const constantProperties& constProps);
-
-            //-
-            void setSiteSizes(label size);
-
-
-        // Access
-
-            //-
-            inline const vector& v() const;
-
-            //-
-            inline vector& v();
-
-            //-
-            inline const vector& a() const;
-
-            //-
-            inline vector& a();
-
-            //-
-            inline const FixedList<vector, 1>& siteForces() const;
-
-            //-
-            inline FixedList<vector, 1>& siteForces();
-
-            //-
-            inline const FixedList<vector, 1>& sitePositions() const;
-
-            //-
-            inline FixedList<vector, 1>& sitePositions();
-
-            //-
-            inline const vector& specialPosition() const;
-
-            //-
-            inline vector& specialPosition();
-
-            //-
-            inline scalar potentialEnergy() const;
-
-            //-
-            inline scalar& potentialEnergy();
-
-            //-
-            inline const tensor& rf() const;
-
-            //-
-            inline tensor& rf();
-
-            //-
-            inline label special() const;
-
-            //-
-            inline bool tethered() const;
-
-            //-
-            inline label id() const;
-
-
-    // Member Operators
-
-        //- Overridable function to handle the particle hitting a patch
-        //  Executed before other patch-hitting functions
-        bool hitPatch
-        (
-            const polyPatch&,
-            trackingData& td,
-            const label patchI,
-            const scalar trackFraction,
-            const tetIndices& tetIs
-        );
-
-        //- Overridable function to handle the particle hitting a processorPatch
-        void hitProcessorPatch
-        (
-            const processorPolyPatch&,
-            trackingData& td
-        );
-
-        //- Overridable function to handle the particle hitting a wallPatch
-        void hitWallPatch
-        (
-            const wallPolyPatch&,
-            trackingData& td,
-            const tetIndices&
-        );
-
-        //- Overridable function to handle the particle hitting a polyPatch
-        void hitPatch
-        (
-            const polyPatch&,
-            trackingData& td
-        );
-
-
-    // I-O
-
-        //- Read
-        static void readFields(Cloud<monoatomic>& mC);
-
-        //- Write
-        static void writeFields(const Cloud<monoatomic>& mC);
-
-        //- Show info
-        static void info(trackingData& td);
-
-
-    // IOstream Operators
-
-        friend Ostream& operator<<(Ostream&, const monoatomic&);
-};
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "monoatomicI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicI.H b/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicI.H
deleted file mode 100644
index 4485d53b8d..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicI.H
+++ /dev/null
@@ -1,328 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-inline Foam::monoatomic::constantProperties::constantProperties()
-:
-    sites_(),
-    electrostaticSites_(-1),
-    pairPotSites_(-1),
-    mass_(0)
-{}
-
-
-inline Foam::monoatomic::constantProperties::constantProperties
-(
-    const dictionary& dict,
-    const List<label>& siteIds
-)
-:
-    sites_(),
-    electrostaticSites_(-1),
-    pairPotSites_(-1),
-    mass_(0)
-{
-    if (siteIds.size() != 1)
-    {
-        FatalErrorIn
-        (
-            "Foam::monoatomic::constantProperties::constantProperties"
-            "("
-                "const dictionary& dict, "
-                "const List<label>& siteIds"
-            ")"
-        )
-            << "monoatomic, single site only, given: " << dict
-            << nl << abort(FatalError);
-    }
-
-    FixedList<word, 1> siteIdNames = dict.lookup("siteIds");
-
-    FixedList<vector, 1> siteReferencePositions
-    (
-        dict.lookup("siteReferencePositions")
-    );
-
-    FixedList<scalar, 1> siteMasses(dict.lookup("siteMasses"));
-
-    FixedList<scalar, 1> siteCharges(dict.lookup("siteCharges"));
-
-    FixedList<word, 1> pairPotentialIds(dict.lookup("pairPotentialSiteIds"));
-
-    constPropSite site = sites_[0];
-
-    site = constPropSite
-    (
-        siteReferencePositions[0],
-        siteMasses[0],
-        siteCharges[0],
-        siteIds[0],
-        siteIdNames[0],
-        (findIndex(pairPotentialIds, siteIdNames[0]) != -1),    // pair
-        (mag(siteCharges[0]) > VSMALL)                          // charge
-    );
-
-    mass_ = site.siteMass();
-
-    if (site.pairPotentialSite())
-    {
-        pairPotSites_[0] = 0;
-    }
-
-    if (site.electrostaticSite())
-    {
-        electrostaticSites_[0] = 0;
-    }
-
-    if (!site.pairPotentialSite() && !site.electrostaticSite())
-    {
-        WarningIn
-        (
-            "Foam::monoatomic::constantProperties::constantProperties"
-            "("
-            "const dictionary& dict"
-            ")"
-        )
-        << siteIdNames[0] << " is a non-interacting site." << endl;
-    }
-
-    // Single site monoatomic - no rotational motion.
-    site.siteReferencePosition() = vector::zero;
-}
-
-
-inline Foam::monoatomic::monoatomic
-(
-    const polyMesh& mesh,
-    const vector& position,
-    const label cellI,
-    const label tetFaceI,
-    const label tetPtI,
-    const scalar temperature,
-    const vector& bulkVelocity,
-    const vector& specialPosition,
-    const constantProperties& cP,
-    monoatomic::trackingData& td,
-    const label special,
-    const label id
-)
-:
-    particle(mesh, position, cellI, tetFaceI, tetPtI),
-    v_(bulkVelocity),
-    a_(vector::zero),
-    specialPosition_(specialPosition),
-    potentialEnergy_(0.0),
-    rf_(tensor::zero),
-    special_(special),
-    id_(id),
-    siteForces_(vector::zero),
-    sitePositions_()
-{
-    setSitePositions(cP);
-
-    v_ += td.cloud().equipartitionLinearVelocity(temperature, cP.mass());
-}
-
-
-inline Foam::monoatomic::monoatomic
-(
-    const polyMesh& mesh,
-    const vector& position,
-    const label cellI,
-    const label tetFaceI,
-    const label tetPtI,
-    const vector& v,
-    const vector& a,
-    const vector& specialPosition,
-    const constantProperties& cP,
-    const label special,
-    const label id
-
-)
-:
-    particle(mesh, position, cellI, tetFaceI, tetPtI),
-    v_(v),
-    a_(a),
-    specialPosition_(specialPosition),
-    potentialEnergy_(0.0),
-    rf_(tensor::zero),
-    special_(special),
-    id_(id),
-    siteForces_(vector::zero),
-    sitePositions_()
-{
-    setSitePositions(cP);
-}
-
-// * * * * * * * constantProperties Member Functions * * * * * * * * * * * * //
-
-inline const Foam::FixedList<Foam::constPropSite, 1>&
-Foam::monoatomic::constantProperties::sites() const
-{
-    return sites_;
-}
-
-
-inline const Foam::FixedList<Foam::label, 1>&
-Foam::monoatomic::constantProperties::pairPotSites() const
-{
-    return pairPotSites_;
-}
-
-
-inline const Foam::FixedList<Foam::label, 1>&
-Foam::monoatomic::constantProperties::electrostaticSites() const
-{
-    return electrostaticSites_;
-}
-
-
-inline Foam::label
-Foam::monoatomic::constantProperties::degreesOfFreedom() const
-{
-    return 3;
-}
-
-
-inline Foam::scalar Foam::monoatomic::constantProperties::mass() const
-{
-    return mass_;
-}
-
-
-inline Foam::label Foam::monoatomic::constantProperties::nSites() const
-{
-    return 1;
-}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-inline const Foam::vector& Foam::monoatomic::v() const
-{
-    return v_;
-}
-
-
-inline Foam::vector& Foam::monoatomic::v()
-{
-    return v_;
-}
-
-
-inline const Foam::vector& Foam::monoatomic::a() const
-{
-    return a_;
-}
-
-
-inline Foam::vector& Foam::monoatomic::a()
-{
-    return a_;
-}
-
-
-inline const Foam::FixedList<Foam::vector, 1>&
-Foam::monoatomic::siteForces() const
-{
-    return siteForces_;
-}
-
-
-inline Foam::FixedList<Foam::vector, 1>& Foam::monoatomic::siteForces()
-{
-    return siteForces_;
-}
-
-
-inline const Foam::FixedList<Foam::vector, 1>&
-Foam::monoatomic::sitePositions() const
-{
-    return sitePositions_;
-}
-
-
-inline Foam::FixedList<Foam::vector, 1>& Foam::monoatomic::sitePositions()
-{
-    return sitePositions_;
-}
-
-
-inline const Foam::vector& Foam::monoatomic::specialPosition() const
-{
-    return specialPosition_;
-}
-
-
-inline Foam::vector& Foam::monoatomic::specialPosition()
-{
-    return specialPosition_;
-}
-
-
-inline Foam::scalar Foam::monoatomic::potentialEnergy() const
-{
-    return potentialEnergy_;
-}
-
-
-inline Foam::scalar& Foam::monoatomic::potentialEnergy()
-{
-    return potentialEnergy_;
-}
-
-
-inline const Foam::tensor& Foam::monoatomic::rf() const
-{
-    return rf_;
-}
-
-
-inline Foam::tensor& Foam::monoatomic::rf()
-{
-    return rf_;
-}
-
-
-inline Foam::label Foam::monoatomic::special() const
-{
-    return special_;
-}
-
-
-inline bool Foam::monoatomic::tethered() const
-{
-    return special_ == SPECIAL_TETHERED;
-}
-
-
-inline Foam::label Foam::monoatomic::id() const
-{
-    return id_;
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicIO.C b/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicIO.C
deleted file mode 100644
index a9a27bb992..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicIO.C
+++ /dev/null
@@ -1,305 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "monoatomic.H"
-#include "IOstreams.H"
-#include "Cloud.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::monoatomic::monoatomic
-(
-    const polyMesh& mesh,
-    Istream& is,
-    bool readFields
-)
-:
-    particle(mesh, is, readFields),
-    v_(vector::zero),
-    a_(vector::zero),
-    specialPosition_(vector::zero),
-    potentialEnergy_(0.0),
-    rf_(tensor::zero),
-    special_(0),
-    id_(0),
-    siteForces_(),
-    sitePositions_()
-{
-    if (readFields)
-    {
-        if (is.format() == IOstream::ASCII)
-        {
-            is  >> v_;
-            is  >> a_;
-            is  >> siteForces_;
-            potentialEnergy_ = readScalar(is);
-            is  >> rf_;
-            special_ = readLabel(is);
-            id_ = readLabel(is);
-            is  >> sitePositions_;
-            is  >> specialPosition_;
-        }
-        else
-        {
-            is.read
-            (
-                reinterpret_cast<char*>(&v_),
-                sizeof(v_)
-              + sizeof(a_)
-              + sizeof(specialPosition_)
-              + sizeof(potentialEnergy_)
-              + sizeof(rf_)
-              + sizeof(special_)
-              + sizeof(id_)
-            );
-
-            is  >> siteForces_ >> sitePositions_;
-        }
-    }
-
-    // Check state of Istream
-    is.check
-    (
-        "Foam::monoatomic::monoatomic"
-        "("
-            "const polyMesh& mesh,"
-            "Istream& is,"
-            "bool readFields"
-        ")"
-    );
-}
-
-
-void Foam::monoatomic::readFields(Cloud<monoatomic>& mC)
-{
-    if (!mC.size())
-    {
-        return;
-    }
-
-    particle::readFields(mC);
-
-    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> 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);
-
-    label i = 0;
-
-    forAllIter(typename Cloud<monoatomic>, mC, iter)
-    {
-        monoatomic& mol = iter();
-
-        mol.v_ = v[i];
-        mol.a_ = a[i];
-        mol.specialPosition_ = specialPosition[i];
-        mol.special_ = special[i];
-        mol.id_ = id[i];
-        i++;
-    }
-}
-
-
-void Foam::monoatomic::writeFields(const Cloud<monoatomic>& mC)
-{
-    particle::writeFields(mC);
-
-    label np = mC.size();
-
-    IOField<vector> v(mC.fieldIOobject("v", IOobject::NO_READ), np);
-    IOField<vector> a(mC.fieldIOobject("a", 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);
-
-    label i = 0;
-    forAllConstIter(typename Cloud<monoatomic>, mC, iter)
-    {
-        const monoatomic& mol = iter();
-
-        v[i] = mol.v_;
-        a[i] = mol.a_;
-        specialPosition[i] = mol.specialPosition_;
-        special[i] = mol.special_;
-        id[i] = mol.id_;
-        i++;
-    }
-
-    v.write();
-    a.write();
-    specialPosition.write();
-    special.write();
-    id.write();
-}
-
-
-void Foam::monoatomic::info(monoatomic::trackingData& td)
-{
-    vector totalLinearMomentum(vector::zero);
-    scalar maxVelocityMag = 0.0;
-    scalar totalMass = 0.0;
-    scalar totalLinearKE = 0.0;
-    scalar totalPE = 0.0;
-    scalar totalrDotf = 0.0;
-    label nMols = td.cloud().size();
-
-    forAllConstIter(typename Cloud<monoatomic>, td.cloud(), mol)
-    {
-        const label molId = mol().id();
-        scalar molMass(td.cloud().constProps(molId).mass());
-        totalMass += molMass;
-    }
-
-    forAllConstIter(typename Cloud<monoatomic>, td.cloud(), mol)
-    {
-        const label molId = mol().id();
-        const monoatomic::constantProperties cP
-        (
-            td.cloud().constProps(molId)
-        );
-        scalar molMass(cP.mass());
-        const vector& molV(mol().v());
-        totalLinearMomentum += molV * molMass;
-        if (mag(molV) > maxVelocityMag)
-        {
-            maxVelocityMag = mag(molV);
-        }
-        totalLinearKE += 0.5*molMass*magSqr(molV);
-        totalPE += mol().potentialEnergy();
-        totalrDotf += tr(mol().rf());
-    }
-
-    scalar meshVolume = sum(td.cloud().mesh().cellVolumes());
-
-    if (Pstream::parRun())
-    {
-        reduce(totalLinearMomentum, sumOp<vector>());
-        reduce(maxVelocityMag, maxOp<scalar>());
-        reduce(totalMass, sumOp<scalar>());
-        reduce(totalLinearKE, sumOp<scalar>());
-        reduce(totalPE, sumOp<scalar>());
-        reduce(totalrDotf, sumOp<scalar>());
-        reduce(nMols, sumOp<label>());
-        reduce(meshVolume, sumOp<scalar>());
-    }
-
-    if (nMols)
-    {
-        Info<< nl << "Number of molecules in " << td.cloud().name() << " = "
-            << nMols << nl
-            << "    Overall number density = "
-            << nMols/meshVolume << nl
-            << "    Overall mass density = "
-            << totalMass/meshVolume << nl
-            << "    Average linear momentum per molecule = "
-            << totalLinearMomentum/nMols << ' '
-            << mag(totalLinearMomentum)/nMols << nl
-            << "    maximum |velocity| = "
-            << maxVelocityMag << nl
-            << "    Average linear KE per molecule = "
-            << totalLinearKE/nMols << nl
-            << "    Average angular KE per molecule = "
-            << totalPE/nMols << nl
-            << "    Average TE per molecule = "
-            <<
-            (
-                  totalLinearKE
-                + totalPE
-            )
-            /nMols
-            << nl << endl;
-    }
-    else
-    {
-        Info<< nl << "No molecules in " << td.cloud().name() << endl;
-    }
-}
-
-
-// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
-
-Foam::Ostream& Foam::operator<<(Ostream& os, const monoatomic& mol)
-{
-    if (os.format() == IOstream::ASCII)
-    {
-        os  << token::SPACE << static_cast<const particle&>(mol)
-            << token::SPACE << mol.face()
-            << token::SPACE << mol.stepFraction()
-            << token::SPACE << mol.v_
-            << token::SPACE << mol.a_
-            << token::SPACE << mol.specialPosition_
-            << token::SPACE << mol.potentialEnergy_
-            << token::SPACE << mol.rf_
-            << token::SPACE << mol.special_
-            << token::SPACE << mol.id_
-            << token::SPACE << mol.siteForces_
-            << token::SPACE << mol.sitePositions_;
-    }
-    else
-    {
-        os  << static_cast<const particle&>(mol);
-        os.write
-        (
-            reinterpret_cast<const char*>(&mol.v_),
-            sizeof(mol.v_)
-          + sizeof(mol.a_)
-          + sizeof(mol.specialPosition_)
-          + sizeof(mol.potentialEnergy_)
-          + sizeof(mol.rf_)
-          + sizeof(mol.special_)
-          + sizeof(mol.id_)
-        );
-        os  << mol.siteForces_ << mol.sitePositions_;
-    }
-
-    // Check state of Ostream
-    os.check
-    (
-        "Foam::Ostream& Foam::operator<<"
-        "(Foam::Ostream&, const Foam::monoatomic&)"
-    );
-
-    return os;
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.C b/src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.C
deleted file mode 100644
index 18fbfbb5c9..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.C
+++ /dev/null
@@ -1,320 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*----------------------------------------------------------------------------*/
-
-#include "polyatomic.H"
-#include "Random.H"
-#include "Time.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(polyatomic, 0);
-    defineTemplateTypeNameAndDebug(Cloud<polyatomic>, 0);
-}
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-Foam::tensor Foam::polyatomic::rotationTensorX(scalar phi) const
-{
-    return tensor
-    (
-        1, 0, 0,
-        0, Foam::cos(phi), -Foam::sin(phi),
-        0, Foam::sin(phi), Foam::cos(phi)
-    );
-}
-
-
-Foam::tensor Foam::polyatomic::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::polyatomic::rotationTensorZ(scalar phi) const
-{
-    return tensor
-    (
-        Foam::cos(phi), -Foam::sin(phi), 0,
-        Foam::sin(phi), Foam::cos(phi), 0,
-        0, 0, 1
-    );
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::polyatomic::move
-(
-    polyatomic::trackingData& td,
-    const scalar trackTime
-)
-{
-    td.switchProcessor = false;
-    td.keepParticle = true;
-
-    if (special_ != SPECIAL_FROZEN)
-    {
-        return td.keepParticle;
-    }
-
-    const constantProperties& constProps(td.cloud().constProps(id_));
-
-    if (td.part() == trackingData::tpFirstVelocityHalfStep)
-    {
-        // First leapfrog velocity adjust part, required before tracking+force
-        // part
-
-        v_ += 0.5*trackTime*a_;
-
-        pi_ += 0.5*trackTime*tau_;
-    }
-    else if (td.part() == trackingData::tpLinearTrack)
-    {
-        // Leapfrog tracking part
-
-        scalar tEnd = (1.0 - stepFraction())*trackTime;
-        scalar dtMax = tEnd;
-
-        while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
-        {
-            // set the lagrangian time-step
-            scalar dt = min(dtMax, tEnd);
-
-            dt *= trackToFace(position() + dt*v_, td);
-
-            tEnd -= dt;
-            stepFraction() = 1.0 - tEnd/trackTime;
-        }
-    }
-    else if (td.part() == trackingData::tpRotationalTrack)
-    {
-        // 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*trackTime*pi_.x()/momentOfInertia.xx());
-                pi_ = pi_ & R;
-                Q_ = Q_ & R;
-            }
-
-            R = rotationTensorY(0.5*trackTime*pi_.y()/momentOfInertia.yy());
-            pi_ = pi_ & R;
-            Q_ = Q_ & R;
-
-            R = rotationTensorZ(trackTime*pi_.z()/momentOfInertia.zz());
-            pi_ = pi_ & R;
-            Q_ = Q_ & R;
-
-            R = rotationTensorY(0.5*trackTime*pi_.y()/momentOfInertia.yy());
-            pi_ = pi_ & R;
-            Q_ = Q_ & R;
-
-            if (!constProps.linearMolecule())
-            {
-                R = rotationTensorX(0.5*trackTime*pi_.x()/momentOfInertia.xx());
-                pi_ = pi_ & R;
-                Q_ = Q_ & R;
-            }
-        }
-
-        setSitePositions(constProps);
-    }
-    else if (td.part() == trackingData::tpSecondVelocityHalfStep)
-    {
-        // Second leapfrog velocity adjust part, required after tracking+force
-        // part
-
-        scalar m = constProps.mass();
-
-        a_ = vector::zero;
-
-        tau_ = vector::zero;
-
-        forAll(siteForces_, sI)
-        {
-            const vector& f = siteForces_[sI];
-
-            a_ += f/m;
-
-            tau_ +=
-                constProps.sites()[sI].siteReferencePosition()
-              ^ (Q_.T() & f);
-        }
-
-        v_ += 0.5*trackTime*a_;
-
-        pi_ += 0.5*trackTime*tau_;
-
-        if (constProps.pointMolecule())
-        {
-            tau_ = vector::zero;
-
-            pi_ = vector::zero;
-        }
-
-        if (constProps.linearMolecule())
-        {
-            tau_.x() = 0.0;
-
-            pi_.x() = 0.0;
-        }
-    }
-    else
-    {
-        FatalErrorIn("polyatomic::move(trackingData&, const scalar)") << nl
-            << td.part() << " is an invalid part of the integration method."
-            << abort(FatalError);
-    }
-
-    return td.keepParticle;
-}
-
-
-void Foam::polyatomic::transformProperties(const tensor& T)
-{
-    particle::transformProperties(T);
-
-    Q_ = T & Q_;
-
-    v_ = transform(T, v_);
-
-    a_ = transform(T, a_);
-
-    pi_ = Q_.T() & transform(T, Q_ & pi_);
-
-    tau_ = Q_.T() & transform(T, Q_ & tau_);
-
-    rf_ = transform(T, rf_);
-
-    sitePositions_ = position_ + (T & (sitePositions_ - position_));
-
-    siteForces_ = T & siteForces_;
-}
-
-
-void Foam::polyatomic::transformProperties(const vector& separation)
-{
-    particle::transformProperties(separation);
-
-    if (special_ == SPECIAL_TETHERED)
-    {
-        specialPosition_ += separation;
-    }
-
-    sitePositions_ = sitePositions_ + separation;
-}
-
-
-void Foam::polyatomic::setSitePositions(const constantProperties& constProps)
-{
-    forAll(constProps.sites(), sI)
-    {
-        sitePositions_[sI] =
-            position_
-          + (Q_ & constProps.sites()[sI].siteReferencePosition());
-    }
-}
-
-
-void Foam::polyatomic::setSiteSizes(label size)
-{
-    sitePositions_.setSize(size);
-
-    siteForces_.setSize(size);
-}
-
-
-bool Foam::polyatomic::hitPatch
-(
-    const polyPatch&,
-    trackingData&,
-    const label,
-    const scalar,
-    const tetIndices&
-)
-{
-    return false;
-}
-
-
-void Foam::polyatomic::hitProcessorPatch
-(
-    const processorPolyPatch&,
-    trackingData& td
-)
-{
-    td.switchProcessor = true;
-}
-
-
-void Foam::polyatomic::hitWallPatch
-(
-    const wallPolyPatch& wpp,
-    trackingData& td,
-    const tetIndices& tetIs
-)
-{
-    // Use of the normal from tetIs is not required as
-    // hasWallImpactDistance is false.
-    vector nw = normal();
-    nw /= mag(nw);
-
-    scalar vn = v_ & nw;
-
-    // Specular reflection
-    if (vn > 0)
-    {
-        v_ -= 2*vn*nw;
-    }
-}
-
-
-void Foam::polyatomic::hitPatch
-(
-    const polyPatch&,
-    trackingData& td
-)
-{
-    td.keepParticle = false;
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.H b/src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.H
deleted file mode 100644
index 0a7797516b..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.H
+++ /dev/null
@@ -1,483 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::polyatomic
-
-Description
-    Foam::polyatomic
-
-SourceFiles
-    polyatomicI.H
-    polyatomic.C
-    polyatomicIO.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef polyatomic_H
-#define polyatomic_H
-
-#include "particle.H"
-#include "IOstream.H"
-#include "autoPtr.H"
-#include "diagTensor.H"
-#include "constPropSite.H"
-#include "MoleculeCloud.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class polyatomic Declaration
-\*---------------------------------------------------------------------------*/
-
-class polyatomic
-:
-    public particle
-{
-
-public:
-
-    // 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 polyatomic constant properties
-    class constantProperties
-    {
-        // Private data
-
-            //- Sites of mass, charge or interaction
-            List<constPropSite> sites_;
-
-            //- Which sites require electrostatic interactions
-            List<label> electrostaticSites_;
-
-            //- Which sites require pair interactions
-            List<label> pairPotSites_;
-
-            //- Moment of intertia (in principal axis configiration)
-            diagTensor momentOfInertia_;
-
-            //-
-            scalar mass_;
-
-
-        // Private Member Functions
-
-            //-
-            bool linearMoleculeTest() const;
-
-
-    public:
-
-        //-
-        inline constantProperties();
-
-        //- Construct from dictionary
-        inline constantProperties
-        (
-            const dictionary& dict,
-            const labelList& siteIds
-        );
-
-        // Member functions
-
-        //-
-        inline const List<constPropSite>& sites() const;
-
-        //-
-        inline const List<label>& pairPotSites() const;
-
-        //-
-        inline const List<label>& electrostaticSites() const;
-
-        //-
-        inline const diagTensor& momentOfInertia() const;
-
-        //-
-        inline bool linearMolecule() const;
-
-        //-
-        inline bool pointMolecule() const;
-
-        //-
-        inline label degreesOfFreedom() const;
-
-        //-
-        inline scalar mass() const;
-
-        //-
-        inline label nSites() const;
-    };
-
-
-    //- Class used to pass tracking data to the trackToFace function
-    class trackingData
-    :
-        public particle::TrackingData<MoleculeCloud<polyatomic> >
-    {
-    public:
-
-        enum trackPart
-        {
-            tpFirstVelocityHalfStep,
-            tpLinearTrack,
-            tpRotationalTrack,
-            tpSecondVelocityHalfStep,
-            tpAccess
-        };
-
-
-    private:
-
-        // label specifying which part of the integration algorithm is taking
-        label part_;
-
-
-    public:
-
-        // Constructors
-
-            trackingData(MoleculeCloud<polyatomic>& cloud, trackPart part)
-            :
-                particle::TrackingData<MoleculeCloud<polyatomic> >(cloud),
-                part_(part)
-            {}
-
-        // Member functions
-
-            inline label part() const
-            {
-                return part_;
-            }
-    };
-
-
-private:
-
-    // Private data
-
-        //- Orientation, stored as the rotation tensor to transform
-        //  from the polyatomic to the global reference frame, i.e.:
-        //  globalVector = Q_ & bodyLocalVector
-        //  bodyLocalVector = Q_.T() & globalVector
-        tensor Q_;
-
-        //- Linear velocity of polyatomic
-        vector v_;
-
-        //- Total linear acceleration of polyatomic
-        vector a_;
-
-        //- Angular momentum of polyatomic, in body local reference frame
-        vector pi_;
-
-        //- Total torque on polyatomic, in body local reference frame
-        vector tau_;
-
-        //-
-        vector specialPosition_;
-
-        //-
-        scalar potentialEnergy_;
-
-        // - r_ij f_ij, stress dyad
-        tensor rf_;
-
-//         // - r_ij outer product f_ij: virial contribution
-//         tensor rDotf_;
-
-        //-
-        label special_;
-
-        //-
-        label id_;
-
-        //-
-        List<vector> siteForces_;
-
-        //-
-        List<vector> sitePositions_;
-
-
-    // Private Member Functions
-
-
-        //-
-        tensor rotationTensorX(scalar deltaT) const;
-
-        //-
-        tensor rotationTensorY(scalar deltaT) const;
-
-        //-
-        tensor rotationTensorZ(scalar deltaT) const;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("polyatomic");
-
-    friend class Cloud<polyatomic>;
-
-    // Constructors
-
-        //- Construct with macroscopic description
-        inline polyatomic
-        (
-            const polyMesh& mesh,
-            const vector& position,
-            const label cellI,
-            const label tetFaceI,
-            const label tetPtI,
-            const scalar temperature,
-            const vector& bulkVelocity,
-            const vector& specialPosition,
-            const constantProperties& constProps,
-            polyatomic::trackingData& td,
-            const label special,
-            const label id
-        );
-
-        //- Construct from all components
-        inline polyatomic
-        (
-            const polyMesh& mesh,
-            const vector& position,
-            const label cellI,
-            const label tetFaceI,
-            const label tetPtI,
-            const tensor& Q,
-            const vector& v,
-            const vector& a,
-            const vector& pi,
-            const vector& tau,
-            const vector& specialPosition,
-            const constantProperties& constProps,
-            const label special,
-            const label id
-        );
-
-        //- Construct from Istream
-        polyatomic
-        (
-            const polyMesh& mesh,
-            Istream& is,
-            bool readFields = true
-        );
-
-        //- Construct and return a clone
-        autoPtr<particle> clone() const
-        {
-            return autoPtr<particle>(new polyatomic(*this));
-        }
-
-        //- Factory class to read-construct particles used for
-        //  parallel transfer
-        class iNew
-        {
-            const polyMesh& mesh_;
-
-        public:
-
-            iNew(const polyMesh& mesh)
-            :
-                mesh_(mesh)
-            {}
-
-            autoPtr<polyatomic> operator()(Istream& is) const
-            {
-                return autoPtr<polyatomic>(new polyatomic(mesh_, is, true));
-            }
-        };
-
-
-    // Member Functions
-
-        // Tracking
-
-            //-
-            bool move(trackingData&, const scalar trackTime);
-
-            //-
-            virtual void transformProperties(const tensor& T);
-
-            //-
-            virtual void transformProperties(const vector& separation);
-
-            //-
-            void setSitePositions(const constantProperties& constProps);
-
-            //-
-            void setSiteSizes(label size);
-
-
-        // Access
-
-            //-
-            inline const tensor& Q() const;
-
-            //-
-            inline tensor& Q();
-
-            //-
-            inline const vector& v() const;
-
-            //-
-            inline vector& v();
-
-            //-
-            inline const vector& a() const;
-
-            //-
-            inline vector& a();
-
-            //-
-            inline const vector& pi() const;
-
-            //-
-            inline vector& pi();
-
-            //-
-            inline const vector& tau() const;
-
-            //-
-            inline vector& tau();
-
-            //-
-            inline const List<vector>& siteForces() const;
-
-            //-
-            inline List<vector>& siteForces();
-
-            //-
-            inline const List<vector>& sitePositions() const;
-
-            //-
-            inline List<vector>& sitePositions();
-
-            //-
-            inline const vector& specialPosition() const;
-
-            //-
-            inline vector& specialPosition();
-
-            //-
-            inline scalar potentialEnergy() const;
-
-            //-
-            inline scalar& potentialEnergy();
-
-            //-
-            inline const tensor& rf() const;
-
-            //-
-            inline tensor& rf();
-
-            //-
-            inline label special() const;
-
-            //-
-            inline bool tethered() const;
-
-            //-
-            inline label id() const;
-
-
-    // Member Operators
-
-        //- Overridable function to handle the particle hitting a patch
-        //  Executed before other patch-hitting functions
-        bool hitPatch
-        (
-            const polyPatch&,
-            trackingData& td,
-            const label patchI,
-            const scalar trackFraction,
-            const tetIndices& tetIs
-        );
-
-        //- Overridable function to handle the particle hitting a processorPatch
-        void hitProcessorPatch
-        (
-            const processorPolyPatch&,
-            trackingData& td
-        );
-
-        //- Overridable function to handle the particle hitting a wallPatch
-        void hitWallPatch
-        (
-            const wallPolyPatch&,
-            trackingData& td,
-            const tetIndices&
-        );
-
-        //- Overridable function to handle the particle hitting a polyPatch
-        void hitPatch
-        (
-            const polyPatch&,
-            trackingData& td
-        );
-
-
-    // I-O
-
-        //- Read
-        static void readFields(Cloud<polyatomic>& mC);
-
-        //- Write
-        static void writeFields(const Cloud<polyatomic>& mC);
-
-        //- Show info
-        static void info(trackingData& td);
-
-
-    // IOstream Operators
-
-        friend Ostream& operator<<(Ostream&, const polyatomic&);
-};
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "polyatomicI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomicI.H b/src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomicI.H
deleted file mode 100644
index 57149e5fa9..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomicI.H
+++ /dev/null
@@ -1,653 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-inline Foam::polyatomic::constantProperties::constantProperties()
-:
-    sites_(),
-    electrostaticSites_(),
-    pairPotSites_(),
-    momentOfInertia_(diagTensor(0, 0, 0)),
-    mass_(0)
-{}
-
-
-inline Foam::polyatomic::constantProperties::constantProperties
-(
-    const dictionary& dict,
-    const labelList& siteIds
-)
-:
-    sites_(),
-    electrostaticSites_(),
-    pairPotSites_(),
-    momentOfInertia_(),
-    mass_(0)
-{
-    List<word> siteIdNames = dict.lookup("siteIds");
-
-    List<vector> siteReferencePositions(dict.lookup("siteReferencePositions"));
-
-    List<scalar> siteMasses(dict.lookup("siteMasses"));
-
-    List<scalar> siteCharges(dict.lookup("siteCharges"));
-
-    List<word> pairPotentialIds(dict.lookup("pairPotentialSiteIds"));
-
-    sites_.setSize(siteReferencePositions.size());
-
-    if
-    (
-        (siteIdNames.size() != sites_.size())
-     || (siteCharges.size() != sites_.size())
-    )
-    {
-            FatalErrorIn
-            (
-                "Foam::polyatomic::constantProperties::constantProperties"
-                "("
-                    "const dictionary& dict"
-                ")"
-            )
-            << "Sizes of site id, charge and "
-            << "referencePositions are not the same: " << nl
-            << siteIdNames  << nl
-            << siteReferencePositions << nl
-            << siteCharges << nl
-            << abort(FatalError);
-    }
-
-    electrostaticSites_.setSize(sites_.size(), -1);
-    pairPotSites_.setSize(sites_.size(), -1);
-
-    label pairI = 0;
-    label electrostaticI = 0;
-
-    forAll(sites_, sI)
-    {
-        sites_[sI] = constPropSite
-        (
-            siteReferencePositions[sI],
-            siteMasses[sI],
-            siteCharges[sI],
-            siteIds[sI],
-            siteIdNames[sI],
-            (findIndex(pairPotentialIds, siteIdNames[sI]) != -1),   // pair
-            (mag(siteCharges[sI]) > VSMALL)                         // charge
-        );
-
-        if (sites_[sI].pairPotentialSite())
-        {
-            pairPotSites_[pairI++] = sI;
-        }
-
-        if (sites_[sI].electrostaticSite())
-        {
-            electrostaticSites_[electrostaticI++] = sI;
-        }
-
-        if (!sites_[sI].pairPotentialSite() && !sites_[sI].electrostaticSite())
-        {
-            WarningIn
-            (
-                "Foam::polyatomic::constantProperties::constantProperties"
-                "("
-                    "const dictionary& dict"
-                ")"
-            )
-                << siteIdNames[sI] << " is a non-interacting site." << endl;
-        }
-    }
-
-    pairPotSites_.setSize(pairI);
-    electrostaticSites_.setSize(electrostaticI);
-
-    // Calculate the centre of mass of the body and subtract it from each
-    // position
-
-    vector centreOfMass(vector::zero);
-
-    forAll(sites_, sI)
-    {
-        mass_ += sites_[sI].siteMass();
-
-        centreOfMass +=
-            sites_[sI].siteReferencePosition()*sites_[sI].siteMass();
-    }
-
-    centreOfMass /= mass_;
-
-    forAll(sites_, sI)
-    {
-        sites_[sI].siteReferencePosition() -= centreOfMass;
-    }
-
-    if (sites_.size() == 1)
-    {
-        // Single site polyatomic - no rotational motion.
-
-        sites_[0].siteReferencePosition() = vector::zero;
-
-        momentOfInertia_ = diagTensor(-1, -1, -1);
-    }
-    else if (linearMoleculeTest())
-    {
-        // Linear polyatomic.
-
-        Info<< nl << "Linear polyatomic." << endl;
-
-        vector dir =
-            sites_[1].siteReferencePosition()
-          - sites_[0].siteReferencePosition();
-
-        dir /= mag(dir);
-
-        tensor Q = rotationTensor(dir, vector(1,0,0));
-
-        forAll(sites_, sI)
-        {
-            sites_[sI].siteReferencePosition() =
-                (Q & sites_[sI].siteReferencePosition());
-        }
-
-        // 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(sites_, sI)
-        {
-            centreOfMass +=
-                sites_[sI].siteReferencePosition()*sites_[sI].siteMass();
-        }
-
-        centreOfMass /= mass_;
-
-        forAll(sites_, sI)
-        {
-            sites_[sI].siteReferencePosition() -= centreOfMass;
-        }
-
-        diagTensor momOfInertia = diagTensor::zero;
-
-        forAll(sites_, sI)
-        {
-            const vector& p(sites_[sI].siteReferencePosition());
-
-            momOfInertia +=
-                sites_[sI].siteMass()*diagTensor(0, p.x()*p.x(), p.x()*p.x());
-        }
-
-        momentOfInertia_ = diagTensor
-         (
-             -1,
-             momOfInertia.yy(),
-             momOfInertia.zz()
-         );
-    }
-    else
-    {
-        // Fully 6DOF polyatomic
-
-        // Calculate the inertia tensor in the current orientation
-
-        tensor momOfInertia(tensor::zero);
-
-        forAll(sites_, sI)
-        {
-            const vector& p(sites_[sI].siteReferencePosition());
-
-            momOfInertia += sites_[sI].siteMass()*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("polyatomic::constantProperties::constantProperties")
-                << "An eigenvalue of the inertia tensor is zero, the "
-                << "polyatomic " << 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
-
-        forAll(sites_, sI)
-        {
-            sites_[sI].siteReferencePosition() =
-                (Q & sites_[sI].siteReferencePosition());
-        }
-
-        // 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(sites_, sI)
-        {
-            centreOfMass +=
-                sites_[sI].siteReferencePosition()*sites_[sI].siteMass();
-        }
-
-        centreOfMass /= mass_;
-
-        forAll(sites_, sI)
-        {
-            sites_[sI].siteReferencePosition() -= centreOfMass;
-        }
-
-        // Calculate the moment of inertia in the principle component
-        // reference frame
-
-        momOfInertia = tensor::zero;
-
-        forAll(sites_, sI)
-        {
-            const vector& p(sites_[sI].siteReferencePosition());
-
-            momOfInertia += sites_[sI].siteMass()*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::polyatomic::polyatomic
-(
-    const polyMesh& mesh,
-    const vector& position,
-    const label cellI,
-    const label tetFaceI,
-    const label tetPtI,
-    const scalar temperature,
-    const vector& bulkVelocity,
-    const vector& specialPosition,
-    const constantProperties& cP,
-    polyatomic::trackingData& td,
-    const label special,
-    const label id
-)
-:
-    particle(mesh, position, cellI, tetFaceI, tetPtI),
-    Q_(I),
-    v_(bulkVelocity),
-    a_(vector::zero),
-    pi_(vector::zero),
-    tau_(vector::zero),
-    specialPosition_(specialPosition),
-    potentialEnergy_(0.0),
-    rf_(tensor::zero),
-    special_(special),
-    id_(id),
-    siteForces_(cP.nSites(), vector::zero),
-    sitePositions_(cP.nSites())
-{
-    setSitePositions(cP);
-
-    v_ += td.cloud().equipartitionLinearVelocity(temperature, cP.mass());
-
-    if (!cP.pointMolecule())
-    {
-        pi_ = td.cloud().equipartitionAngularMomentum(temperature, cP);
-
-        scalar phi(td.cloud().rndGen().scalar01()*twoPi);
-
-        scalar theta(td.cloud().rndGen().scalar01()*twoPi);
-
-        scalar psi(td.cloud().rndGen().scalar01()*twoPi);
-
-        Q_ = tensor
-        (
-            cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi),
-            cos(psi)*sin(phi) + cos(theta)*cos(phi)*sin(psi),
-            sin(psi)*sin(theta),
-            - sin(psi)*cos(phi) - cos(theta)*sin(phi)*cos(psi),
-            - sin(psi)*sin(phi) + cos(theta)*cos(phi)*cos(psi),
-            cos(psi)*sin(theta),
-            sin(theta)*sin(phi),
-            - sin(theta)*cos(phi),
-            cos(theta)
-        );
-    }
-}
-
-
-Foam::polyatomic::polyatomic
-(
-    const polyMesh& mesh,
-    const vector& position,
-    const label cellI,
-    const label tetFaceI,
-    const label tetPtI,
-    const tensor& Q,
-    const vector& v,
-    const vector& a,
-    const vector& pi,
-    const vector& tau,
-    const vector& specialPosition,
-    const constantProperties& cP,
-    const label special,
-    const label id
-)
-:
-    particle(mesh, position, cellI, tetFaceI, tetPtI),
-    Q_(Q),
-    v_(v),
-    a_(a),
-    pi_(pi),
-    tau_(tau),
-    specialPosition_(specialPosition),
-    potentialEnergy_(0.0),
-    rf_(tensor::zero),
-    special_(special),
-    id_(id),
-    siteForces_(cP.nSites(), vector::zero),
-    sitePositions_(cP.nSites())
-{
-    setSitePositions(cP);
-}
-
-
-// * * * * * * * constantProperties Private Member Functions * * * * * * * * //
-
-inline bool Foam::polyatomic::constantProperties::linearMoleculeTest() const
-{
-    if (sites_.size() == 2)
-    {
-        return true;
-    }
-
-    vector refDir =
-        sites_[1].siteReferencePosition()
-      - sites_[0].siteReferencePosition();
-
-    refDir /= mag(refDir);
-
-    for
-    (
-        label i = 2;
-        i < sites_.size();
-        i++
-    )
-    {
-        vector dir =
-            sites_[i].siteReferencePosition()
-          - sites_[i - 1].siteReferencePosition();
-
-        dir /= mag(dir);
-
-        if (mag(refDir & dir) < 1 - SMALL)
-        {
-            return false;
-        }
-    }
-
-    return true;
-}
-
-
-// * * * * * * * constantProperties Member Functions * * * * * * * * * * * * //
-
-inline const Foam::List<Foam::constPropSite>&
-Foam::polyatomic::constantProperties::sites() const
-{
-    return sites_;
-}
-
-
-inline const Foam::List<Foam::label>&
-Foam::polyatomic::constantProperties::pairPotSites() const
-{
-    return pairPotSites_;
-}
-
-
-inline const Foam::List<Foam::label>&
-Foam::polyatomic::constantProperties::electrostaticSites() const
-{
-    return electrostaticSites_;
-}
-
-
-inline const Foam::diagTensor&
-Foam::polyatomic::constantProperties::momentOfInertia() const
-{
-    return momentOfInertia_;
-}
-
-
-inline bool Foam::polyatomic::constantProperties::linearMolecule() const
-{
-    return ((momentOfInertia_.xx() < 0) && (momentOfInertia_.yy() > 0));
-}
-
-
-inline  bool Foam::polyatomic::constantProperties::pointMolecule() const
-{
-    return (momentOfInertia_.zz() < 0);
-}
-
-
-inline Foam::label
-Foam::polyatomic::constantProperties::degreesOfFreedom() const
-{
-    if (linearMolecule())
-    {
-        return 5;
-    }
-    else if (pointMolecule())
-    {
-        return 3;
-    }
-    else
-    {
-        return 6;
-    }
-}
-
-
-inline Foam::scalar Foam::polyatomic::constantProperties::mass() const
-{
-    return mass_;
-}
-
-
-inline Foam::label Foam::polyatomic::constantProperties::nSites() const
-{
-    return sites_.size();
-}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-inline const Foam::tensor& Foam::polyatomic::Q() const
-{
-    return Q_;
-}
-
-
-inline Foam::tensor& Foam::polyatomic::Q()
-{
-    return Q_;
-}
-
-
-inline const Foam::vector& Foam::polyatomic::v() const
-{
-    return v_;
-}
-
-
-inline Foam::vector& Foam::polyatomic::v()
-{
-    return v_;
-}
-
-
-inline const Foam::vector& Foam::polyatomic::a() const
-{
-    return a_;
-}
-
-
-inline Foam::vector& Foam::polyatomic::a()
-{
-    return a_;
-}
-
-
-inline const Foam::vector& Foam::polyatomic::pi() const
-{
-    return pi_;
-}
-
-
-inline Foam::vector& Foam::polyatomic::pi()
-{
-    return pi_;
-}
-
-
-inline const Foam::vector& Foam::polyatomic::tau() const
-{
-    return tau_;
-}
-
-
-inline Foam::vector& Foam::polyatomic::tau()
-{
-    return tau_;
-}
-
-
-inline const Foam::List<Foam::vector>& Foam::polyatomic::siteForces() const
-{
-    return siteForces_;
-}
-
-
-inline Foam::List<Foam::vector>& Foam::polyatomic::siteForces()
-{
-    return siteForces_;
-}
-
-
-inline const Foam::List<Foam::vector>& Foam::polyatomic::sitePositions() const
-{
-    return sitePositions_;
-}
-
-
-inline Foam::List<Foam::vector>& Foam::polyatomic::sitePositions()
-{
-    return sitePositions_;
-}
-
-
-inline const Foam::vector& Foam::polyatomic::specialPosition() const
-{
-    return specialPosition_;
-}
-
-
-inline Foam::vector& Foam::polyatomic::specialPosition()
-{
-    return specialPosition_;
-}
-
-
-inline Foam::scalar Foam::polyatomic::potentialEnergy() const
-{
-    return potentialEnergy_;
-}
-
-
-inline Foam::scalar& Foam::polyatomic::potentialEnergy()
-{
-    return potentialEnergy_;
-}
-
-
-inline const Foam::tensor& Foam::polyatomic::rf() const
-{
-    return rf_;
-}
-
-
-inline Foam::tensor& Foam::polyatomic::rf()
-{
-    return rf_;
-}
-
-
-inline Foam::label Foam::polyatomic::special() const
-{
-    return special_;
-}
-
-
-inline bool Foam::polyatomic::tethered() const
-{
-    return special_ == SPECIAL_TETHERED;
-}
-
-
-inline Foam::label Foam::polyatomic::id() const
-{
-    return id_;
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomicIO.C b/src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomicIO.C
deleted file mode 100644
index 86dcbfe49f..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomicIO.C
+++ /dev/null
@@ -1,414 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "polyatomic.H"
-#include "IOstreams.H"
-#include "Cloud.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::polyatomic::polyatomic
-(
-    const polyMesh& mesh,
-    Istream& is,
-    bool readFields
-)
-:
-    particle(mesh, is, readFields),
-    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)
-        {
-            is  >> Q_;
-            is  >> v_;
-            is  >> a_;
-            is  >> pi_;
-            is  >> tau_;
-            is  >> siteForces_;
-            potentialEnergy_ = readScalar(is);
-            is  >> rf_;
-            special_ = readLabel(is);
-            id_ = readLabel(is);
-            is  >> sitePositions_;
-            is  >> specialPosition_;
-        }
-        else
-        {
-            is.read
-            (
-                reinterpret_cast<char*>(&Q_),
-                sizeof(Q_)
-              + sizeof(v_)
-              + sizeof(a_)
-              + sizeof(pi_)
-              + sizeof(tau_)
-              + sizeof(specialPosition_)
-              + sizeof(potentialEnergy_)
-              + sizeof(rf_)
-              + sizeof(special_)
-              + sizeof(id_)
-            );
-
-            is  >> siteForces_ >> sitePositions_;
-        }
-    }
-
-    // Check state of Istream
-    is.check
-    (
-        "Foam::polyatomic::polyatomic"
-        "("
-            "const polyMesh& mesh,"
-            "Istream& is,"
-            "bool readFields"
-        ")"
-    );
-}
-
-
-void Foam::polyatomic::readFields(Cloud<polyatomic>& mC)
-{
-    if (!mC.size())
-    {
-        return;
-    }
-
-    particle::readFields(mC);
-
-    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);
-
-    label i = 0;
-
-    forAllIter(typename Cloud<polyatomic>, mC, iter)
-    {
-        polyatomic& 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];
-        i++;
-    }
-}
-
-
-void Foam::polyatomic::writeFields(const Cloud<polyatomic>& mC)
-{
-    particle::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);
-
-    // Post processing fields
-
-    IOField<vector> piGlobal
-    (
-        mC.fieldIOobject("piGlobal", IOobject::NO_READ),
-        np
-    );
-
-    IOField<vector> tauGlobal
-    (
-        mC.fieldIOobject("tauGlobal", IOobject::NO_READ),
-        np
-    );
-
-    IOField<vector> orientation1
-    (
-        mC.fieldIOobject("orientation1", IOobject::NO_READ),
-        np
-    );
-
-    IOField<vector> orientation2
-    (
-        mC.fieldIOobject("orientation2", IOobject::NO_READ),
-        np
-    );
-
-    IOField<vector> orientation3
-    (
-        mC.fieldIOobject("orientation3", IOobject::NO_READ),
-        np
-    );
-
-    label i = 0;
-    forAllConstIter(typename Cloud<polyatomic>, mC, iter)
-    {
-        const polyatomic& 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_;
-
-        piGlobal[i] = mol.Q_ & mol.pi_;
-        tauGlobal[i] = mol.Q_ & mol.tau_;
-
-        orientation1[i] = mol.Q_ & vector(1,0,0);
-        orientation2[i] = mol.Q_ & vector(0,1,0);
-        orientation3[i] = mol.Q_ & vector(0,0,1);
-
-        i++;
-    }
-
-    Q.write();
-    v.write();
-    a.write();
-    pi.write();
-    tau.write();
-    specialPosition.write();
-    special.write();
-    id.write();
-
-    piGlobal.write();
-    tauGlobal.write();
-
-    orientation1.write();
-    orientation2.write();
-    orientation3.write();
-}
-
-
-void Foam::polyatomic::info(polyatomic::trackingData& td)
-{
-    vector totalLinearMomentum(vector::zero);
-    vector totalAngularMomentum(vector::zero);
-    scalar maxVelocityMag = 0.0;
-    scalar totalMass = 0.0;
-    scalar totalLinearKE = 0.0;
-    scalar totalAngularKE = 0.0;
-    scalar totalPE = 0.0;
-    scalar totalrDotf = 0.0;
-    //vector CentreOfMass(vector::zero);
-    label nMols = td.cloud().size();
-    label dofs = 0;
-
-    forAllConstIter(typename Cloud<polyatomic>, td.cloud(), mol)
-    {
-        const label molId = mol().id();
-        scalar molMass(td.cloud().constProps(molId).mass());
-        totalMass += molMass;
-        //CentreOfMass += mol().position()*molMass;
-    }
-
-    // if (nMols)
-    // {
-    //     CentreOfMass /= totalMass;
-    // }
-
-    forAllConstIter(typename Cloud<polyatomic>, td.cloud(), mol)
-    {
-        const label molId = mol().id();
-        const polyatomic::constantProperties cP
-        (
-            td.cloud().constProps(molId)
-        );
-        scalar molMass(cP.mass());
-        const diagTensor& molMoI(cP.momentOfInertia());
-        const vector& molV(mol().v());
-        const vector& molOmega(inv(molMoI) & mol().pi());
-        vector molPiGlobal = mol().Q() & mol().pi();
-        totalLinearMomentum += molV * molMass;
-        totalAngularMomentum += molPiGlobal;
-            //+((mol().position() - CentreOfMass) ^ (molV * molMass));
-        if (mag(molV) > maxVelocityMag)
-        {
-            maxVelocityMag = mag(molV);
-        }
-        totalLinearKE += 0.5*molMass*magSqr(molV);
-        totalAngularKE += 0.5*(molOmega & molMoI & molOmega);
-        totalPE += mol().potentialEnergy();
-        totalrDotf += tr(mol().rf());
-        dofs += cP.degreesOfFreedom();
-    }
-
-    scalar meshVolume = sum(td.cloud().mesh().cellVolumes());
-
-    if (Pstream::parRun())
-    {
-        reduce(totalLinearMomentum, sumOp<vector>());
-        reduce(totalAngularMomentum, sumOp<vector>());
-        reduce(maxVelocityMag, maxOp<scalar>());
-        reduce(totalMass, sumOp<scalar>());
-        reduce(totalLinearKE, sumOp<scalar>());
-        reduce(totalAngularKE, sumOp<scalar>());
-        reduce(totalPE, sumOp<scalar>());
-        reduce(totalrDotf, sumOp<scalar>());
-        reduce(nMols, sumOp<label>());
-        reduce(dofs, sumOp<label>());
-        reduce(meshVolume, sumOp<scalar>());
-    }
-
-    if (nMols)
-    {
-        Info<< nl << "Number of molecules in " << td.cloud().name() << " = "
-            << nMols << nl
-            << "    Overall number density = "
-            << nMols/meshVolume << nl
-            << "    Overall mass density = "
-            << totalMass/meshVolume << nl
-            << "    Average linear momentum per molecule = "
-            << totalLinearMomentum/nMols << ' '
-            << mag(totalLinearMomentum)/nMols << nl
-            << "    Average angular momentum per molecule = "
-            << totalAngularMomentum << ' '
-            << mag(totalAngularMomentum)/nMols << nl
-            << "    maximum |velocity| = "
-            << maxVelocityMag << nl
-            << "    Average linear KE per molecule = "
-            << totalLinearKE/nMols << nl
-            << "    Average angular KE per molecule = "
-            << totalAngularKE/nMols << nl
-            << "    Average PE per molecule = "
-            << totalPE/nMols << nl
-            << "    Average TE per molecule = "
-            <<
-            (
-                  totalLinearKE
-                + totalAngularKE
-                + totalPE
-            )
-            /nMols
-            << nl << endl;
-    }
-    else
-    {
-        Info<< nl << "No molecules in " << td.cloud().name() << endl;
-    }
-}
-
-
-// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
-
-Foam::Ostream& Foam::operator<<(Ostream& os, const polyatomic& mol)
-{
-    if (os.format() == IOstream::ASCII)
-    {
-        os  << token::SPACE << static_cast<const particle&>(mol)
-            << token::SPACE << mol.face()
-            << token::SPACE << mol.stepFraction()
-            << token::SPACE << mol.Q_
-            << 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.special_
-            << token::SPACE << mol.id_
-            << token::SPACE << mol.siteForces_
-            << token::SPACE << mol.sitePositions_;
-    }
-    else
-    {
-        os  << static_cast<const particle&>(mol);
-        os.write
-        (
-            reinterpret_cast<const char*>(&mol.Q_),
-            sizeof(mol.Q_)
-          + sizeof(mol.v_)
-          + sizeof(mol.a_)
-          + sizeof(mol.pi_)
-          + sizeof(mol.tau_)
-          + sizeof(mol.specialPosition_)
-          + sizeof(mol.potentialEnergy_)
-          + sizeof(mol.rf_)
-          + sizeof(mol.special_)
-          + sizeof(mol.id_)
-        );
-        os  << mol.siteForces_ << mol.sitePositions_;
-    }
-
-    // Check state of Ostream
-    os.check
-    (
-        "Foam::Ostream& Foam::operator<<"
-        "(Foam::Ostream&, const Foam::polyatomic&)"
-    );
-
-    return os;
-}
-
-
-// ************************************************************************* //

From 6be1a5d417ed922ce36ce2625d7cd758666b8eff Mon Sep 17 00:00:00 2001
From: sergio <s.ferraris@opencfd.co.uk>
Date: Fri, 22 Jul 2011 12:44:23 +0100
Subject: [PATCH 2/2] Adding directionality to porous media

---
 .../porousBafflePressureFvPatchFields.C                     | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/turbulenceModels/derivedFvPatchFields/porousBafflePressure/porousBafflePressureFvPatchFields.C b/src/turbulenceModels/derivedFvPatchFields/porousBafflePressure/porousBafflePressureFvPatchFields.C
index f82551589d..1652ab878e 100644
--- a/src/turbulenceModels/derivedFvPatchFields/porousBafflePressure/porousBafflePressureFvPatchFields.C
+++ b/src/turbulenceModels/derivedFvPatchFields/porousBafflePressure/porousBafflePressureFvPatchFields.C
@@ -66,6 +66,8 @@ void Foam::porousBafflePressureFvPatchField<Foam::scalar>::updateCoeffs()
 
     scalarField Un(phip/patch().magSf());
 
+    scalarField magUn(mag(Un));
+
     if (phi.dimensions() == dimensionSet(0, 3, -1, 0, 0))
     {
         const incompressible::turbulenceModel& model =
@@ -76,7 +78,7 @@ void Foam::porousBafflePressureFvPatchField<Foam::scalar>::updateCoeffs()
 
         const scalarField nuEffw = model.nuEff()().boundaryField()[patchI];
 
-        jump_ = -(I_*nuEffw*mag(Un) + D_*0.5*magSqr(Un)*length_);
+        jump_ = -sign(Un)*(I_*nuEffw + D_*0.5*magUn*length_)*magUn;
     }
     else
     {
@@ -93,7 +95,7 @@ void Foam::porousBafflePressureFvPatchField<Foam::scalar>::updateCoeffs()
 
         Un /= rhow;
 
-        jump_ = -(I_*muEffw*mag(Un) + D_*0.5*rhow*magSqr(Un)*length_);
+        jump_ = -sign(Un)*(I_*muEffw + D_*0.5*rhow*magUn*length_)*magUn;
     }
 
     if (debug)