diff --git a/applications/utilities/preProcessing/mdInitialise/mdInitialise.C b/applications/utilities/preProcessing/mdInitialise/mdInitialise.C
index 5d06f35fb6..10a6602181 100644
--- a/applications/utilities/preProcessing/mdInitialise/mdInitialise.C
+++ b/applications/utilities/preProcessing/mdInitialise/mdInitialise.C
@@ -2,7 +2,7 @@
  =========                   |
  \\      /   F ield          | OpenFOAM: The Open Source CFD Toolbox
   \\    /    O peration      |
-   \\  /     A nd            | Copyright (C) 2008-2010 OpenCFD Ltd.
+   \\  /     A nd            | Copyright (C) 2008-2011 OpenCFD Ltd.
     \\/      M anipulation   |
 -------------------------------------------------------------------------------
 License
diff --git a/src/lagrangian/molecularDynamics/molecule/Make/files b/src/lagrangian/molecularDynamics/molecule/Make/files
index 7f8abcb0f3..94c3ffd8e8 100644
--- a/src/lagrangian/molecularDynamics/molecule/Make/files
+++ b/src/lagrangian/molecularDynamics/molecule/Make/files
@@ -5,4 +5,8 @@ polyatomic/polyatomic/polyatomicIO.C
 
 polyatomic/polyatomicCloud/polyatomicCloud.C
 
+/* controllers/basic/controllers/controllers.C
+controllers/basic/stateController/stateController.C
+controllers/basic/fluxController/fluxController.C */
+
 LIB = $(FOAM_LIBBIN)/libmolecule
diff --git a/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.C b/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.C
index 5647869233..08ea263dd6 100644
--- a/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.C
+++ b/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.C
@@ -21,12 +21,7 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Class
-    constPropSite
-
-Description
-
-\*----------------------------------------------------------------------------*/
+\*---------------------------------------------------------------------------*/
 
 #include "constPropSite.H"
 
diff --git a/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.H b/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.H
index c53eba2ee8..910c0ccaef 100644
--- a/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.H
+++ b/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.H
@@ -20,11 +20,13 @@ License
 
     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
diff --git a/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSiteI.H b/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSiteI.H
index b4f1ad9177..9c648c94c8 100644
--- a/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSiteI.H
+++ b/src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSiteI.H
@@ -21,11 +21,6 @@ License
     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
-
 \*---------------------------------------------------------------------------*/
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.C b/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.C
new file mode 100644
index 0000000000..f0bad336ee
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.C
@@ -0,0 +1,573 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
new file mode 100644
index 0000000000..22e9ca7cca
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.H
@@ -0,0 +1,163 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
new file mode 100644
index 0000000000..66c10f0784
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllersI.H
@@ -0,0 +1,64 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
new file mode 100644
index 0000000000..711b661758
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.C
@@ -0,0 +1,524 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
new file mode 100644
index 0000000000..9b8a375277
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.H
@@ -0,0 +1,272 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
new file mode 100644
index 0000000000..b98aa126a3
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.C
@@ -0,0 +1,490 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
new file mode 100644
index 0000000000..3e142c3962
--- /dev/null
+++ b/src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.H
@@ -0,0 +1,270 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/polyatomic/polyatomicCloud/polyatomicCloud.C b/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloud.C
index 0ad96b90c8..2bfb744eb0 100644
--- a/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloud.C
+++ b/src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloud.C
@@ -1340,7 +1340,7 @@ void Foam::polyatomicCloud::info() const
 
     if (nMols)
     {
-        Info<< "Number of molecules in " << this->name() << "  = "
+        Info<< "Number of molecules in " << this->name() << " = "
             << nMols << nl
             << "Overall number density = "
             << nMols/meshVolume << nl