diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
index bfad310038..9eddbe4efd 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
@@ -268,14 +268,9 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const
 
     phaseSystem::heatTransferTable& eqns = eqnsPtr();
 
-    forAllConstIter
-    (
-        phaseSystem::phaseModelList,
-        this->phaseModels_,
-        phaseModelIter
-    )
+    forAll(this->phaseModels_, phasei)
     {
-        const phaseModel& phase(phaseModelIter());
+        const phaseModel& phase = this->phaseModels_[phasei];
 
         eqns.insert
         (
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C
index a0881ee33a..f4a1284db8 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C
@@ -120,14 +120,9 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const
 
     phaseSystem::heatTransferTable& eqns = eqnsPtr();
 
-    forAllConstIter
-    (
-        phaseSystem::phaseModelList,
-        this->phaseModels_,
-        phaseModelIter
-    )
+    forAll(this->phaseModels_, phasei)
     {
-        const phaseModel& phase(phaseModelIter());
+        const phaseModel& phase = this->phaseModels_[phasei];
 
         eqns.insert
         (
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C
index 801ceb8caa..fe0343b2ac 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C
@@ -69,14 +69,9 @@ massTransfer() const
 
     phaseSystem::massTransferTable& eqns = eqnsPtr();
 
-    forAllConstIter
-    (
-        phaseSystem::phaseModelList,
-        this->phaseModels_,
-        phaseModelIter
-    )
+    forAll(this->phaseModels_, phasei)
     {
-        const phaseModel& phase(phaseModelIter());
+        const phaseModel& phase = this->phaseModels_[phasei];
 
         const PtrList<volScalarField>& Yi = phase.Y();
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
index 0f5470d26e..dc9335f7ea 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
@@ -413,14 +413,9 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
 
     phaseSystem::momentumTransferTable& eqns = eqnsPtr();
 
-    forAllConstIter
-    (
-        phaseSystem::phaseModelList,
-        this->phaseModels_,
-        phaseModelIter
-    )
+    forAll(this->phaseModels_, phasei)
     {
-        const phaseModel& phase(phaseModelIter());
+        const phaseModel& phase = this->phaseModels_[phasei];
 
         eqns.insert
         (
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
index 1eb96b3b26..f0453940c7 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
@@ -64,14 +64,9 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::massTransfer() const
 
     phaseSystem::massTransferTable& eqns = eqnsPtr();
 
-    forAllConstIter
-    (
-        phaseSystem::phaseModelList,
-        this->phaseModels_,
-        phaseModelIter
-    )
+    forAll(this->phaseModels_, phasei)
     {
-        const phaseModel& phase(phaseModelIter());
+        const phaseModel& phase = this->phaseModels_[phasei];
 
         const PtrList<volScalarField>& Yi = phase.Y();
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
index ae23455d88..aa9380758d 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
@@ -46,24 +46,19 @@ Foam::tmp<Foam::surfaceScalarField> Foam::phaseSystem::calcPhi
     const phaseModelList& phaseModels
 ) const
 {
-    phaseModelList::const_iterator phaseModelIter = phaseModels.begin();
-
     tmp<surfaceScalarField> tmpPhi
     (
         new surfaceScalarField
         (
             "phi",
-            fvc::interpolate(phaseModelIter())*phaseModelIter().phi()
+            fvc::interpolate(phaseModels[0])*phaseModels[0].phi()
         )
     );
 
-    ++phaseModelIter;
-
-    for (; phaseModelIter != phaseModels.end(); ++ phaseModelIter)
+    for (label phasei=1; phasei<phaseModels.size(); phasei++)
     {
         tmpPhi() +=
-            fvc::interpolate(phaseModelIter())
-           *phaseModelIter().phi();
+            fvc::interpolate(phaseModels[phasei])*phaseModels[phasei].phi();
     }
 
     return tmpPhi;
@@ -196,16 +191,14 @@ Foam::phaseSystem::~phaseSystem()
 
 Foam::tmp<Foam::volScalarField> Foam::phaseSystem::rho() const
 {
-    phaseModelList::const_iterator phaseModelIter = phaseModels_.begin();
-
     tmp<volScalarField> tmpRho
     (
-        phaseModelIter()*phaseModelIter().rho()
+        phaseModels_[0]*phaseModels_[0].rho()
     );
 
-    for (; phaseModelIter != phaseModels_.end(); ++ phaseModelIter)
+    for (label phasei=1; phasei<phaseModels_.size(); phasei++)
     {
-        tmpRho() += phaseModelIter()*phaseModelIter().rho();
+        tmpRho() += phaseModels_[phasei]*phaseModels_[phasei].rho();
     }
 
     return tmpRho;
@@ -214,16 +207,14 @@ Foam::tmp<Foam::volScalarField> Foam::phaseSystem::rho() const
 
 Foam::tmp<Foam::volVectorField> Foam::phaseSystem::U() const
 {
-    phaseModelList::const_iterator phaseModelIter = phaseModels_.begin();
-
     tmp<volVectorField> tmpU
     (
-        phaseModelIter()*phaseModelIter().U()
+        phaseModels_[0]*phaseModels_[0].U()
     );
 
-    for (; phaseModelIter != phaseModels_.end(); ++ phaseModelIter)
+    for (label phasei=1; phasei<phaseModels_.size(); phasei++)
     {
-        tmpU() += phaseModelIter()*phaseModelIter().U();
+        tmpU() += phaseModels_[phasei]*phaseModels_[phasei].U();
     }
 
     return tmpU;
@@ -266,9 +257,9 @@ void Foam::phaseSystem::solve()
 
 void Foam::phaseSystem::correct()
 {
-    forAllIter(phaseModelList, phaseModels_, phaseModelIter)
+    forAll(phaseModels_, phasei)
     {
-        phaseModelIter().correct();
+        phaseModels_[phasei].correct();
     }
 }
 
@@ -277,11 +268,11 @@ void Foam::phaseSystem::correctKinematics()
 {
     bool updateDpdt = false;
 
-    forAllIter(phaseModelList, phaseModels_, phaseModelIter)
+    forAll(phaseModels_, phasei)
     {
-        phaseModelIter().correctKinematics();
+        phaseModels_[phasei].correctKinematics();
 
-        updateDpdt = updateDpdt || phaseModelIter().thermo().dpdt();
+        updateDpdt = updateDpdt || phaseModels_[phasei].thermo().dpdt();
     }
 
     // Update the pressure time-derivative if required
@@ -294,27 +285,27 @@ void Foam::phaseSystem::correctKinematics()
 
 void Foam::phaseSystem::correctThermo()
 {
-    forAllIter(phaseModelList, phaseModels_, phaseModelIter)
+    forAll(phaseModels_, phasei)
     {
-        phaseModelIter().correctThermo();
+        phaseModels_[phasei].correctThermo();
     }
 }
 
 
 void Foam::phaseSystem::correctTurbulence()
 {
-    forAllIter(phaseModelList, phaseModels_, phaseModelIter)
+    forAll(phaseModels_, phasei)
     {
-        phaseModelIter().correctTurbulence();
+        phaseModels_[phasei].correctTurbulence();
     }
 }
 
 
 void Foam::phaseSystem::correctEnergyTransport()
 {
-    forAllIter(phaseModelList, phaseModels_, phaseModelIter)
+    forAll(phaseModels_, phasei)
     {
-        phaseModelIter().correctEnergyTransport();
+        phaseModels_[phasei].correctEnergyTransport();
     }
 }
 
@@ -325,9 +316,9 @@ bool Foam::phaseSystem::read()
     {
         bool readOK = true;
 
-        forAllIter(phaseModelList, phaseModels_, phaseModelIter)
+        forAll(phaseModels_, phasei)
         {
-            readOK &= phaseModelIter().read();
+            readOK &= phaseModels_[phasei].read();
         }
 
         // models ...
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemTemplates.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemTemplates.H
index 37a3627e1c..eee3fcd3b3 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemTemplates.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemTemplates.H
@@ -151,12 +151,12 @@ void Foam::phaseSystem::generatePairsAndSubModels
         HashTable<autoPtr<modelType>, phasePairKey, phasePairKey::hash>
         modelTypeTable;
 
-    forAllConstIter(phaseModelList, phaseModels_, phaseModelIter)
+    forAll(phaseModels_, phasei)
     {
         modelTypeTable tempModels;
         generatePairsAndSubModels
         (
-            IOobject::groupName(modelName, phaseModelIter().name()),
+            IOobject::groupName(modelName, phaseModels_[phasei].name()),
             tempModels
         );
 
@@ -175,7 +175,7 @@ void Foam::phaseSystem::generatePairsAndSubModels
 
             models[tempModelIter.key()].insert
             (
-                phaseModelIter().name(),
+                phaseModels_[phasei].name(),
                 *tempModelIter
             );
         }