diff --git a/.gitignore b/.gitignore
index ebb0a8e..e8b9455 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,6 @@
 [._]*.sw[a-p]
 [._]s[a-v][a-z]
 [._]sw[a-p]
+Make
+!Make/file
+!Make/option
diff --git a/diffusivityModel/Neutral/Neutral.H b/diffusivityModel/Neutral/Neutral.H
index f707a2d..4b1d4fb 100644
--- a/diffusivityModel/Neutral/Neutral.H
+++ b/diffusivityModel/Neutral/Neutral.H
@@ -136,6 +136,8 @@ public:
 
         inline scalar Zrot() const;
 
+        inline scalar Zrot(scalar T) const;
+
         // Check
 
         // Edit
diff --git a/diffusivityModel/Neutral/NeutralI.H b/diffusivityModel/Neutral/NeutralI.H
index df32e98..f35f5d4 100644
--- a/diffusivityModel/Neutral/NeutralI.H
+++ b/diffusivityModel/Neutral/NeutralI.H
@@ -84,6 +84,20 @@ inline Foam::scalar Foam::Neutral::Zrot()
     return Zrot_;
 }
 
+inline Foam::scalar Foam::Neutral::Zrot(const scalar T)
+    const
+{
+    scalar a3 = sqrt(pow3(pi));
+    scalar a1 = a3/2.0;
+    scalar a2 = sqr(pi) / 4.0 + 2.0;
+
+    scalar t1 = sqrt(wellDepth_/T);
+    scalar t2 = wellDepth_/T;
+    scalar t3 = t1 * t2;
+
+    return 1.0 + a1 * t1 + a2 * t2 + a3 * t3;
+}
+
 
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
diff --git a/diffusivityModel/Particle/Particle.H b/diffusivityModel/Particle/Particle.H
index 7c06456..a378e72 100644
--- a/diffusivityModel/Particle/Particle.H
+++ b/diffusivityModel/Particle/Particle.H
@@ -68,8 +68,8 @@ public:
     enum Geometry
     {
         ATOM = 0,
-        LINEAR = 1,
-        NONLINEAR = 2
+        LINEAR = 2,
+        NONLINEAR = 3
     };
 
 
@@ -150,6 +150,8 @@ public:
 
         inline scalar z() const;
 
+        inline scalar F() const;
+
         inline Geometry geometry() const;
 
         // Check
diff --git a/diffusivityModel/Particle/ParticleI.H b/diffusivityModel/Particle/ParticleI.H
index 9f76076..702033a 100644
--- a/diffusivityModel/Particle/ParticleI.H
+++ b/diffusivityModel/Particle/ParticleI.H
@@ -57,6 +57,13 @@ inline Foam::scalar Foam::Particle::z()
 }
 
 
+inline Foam::scalar Foam::Particle::F()
+    const
+{
+    return (geometry_);
+}
+
+
 inline Foam::Particle::Geometry Foam::Particle::geometry()
     const
 {
diff --git a/diffusivityModel/diffusivityModel/diffusivityModel.C b/diffusivityModel/diffusivityModel/diffusivityModel.C
index 15c48ee..2cff8c3 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModel.C
+++ b/diffusivityModel/diffusivityModel/diffusivityModel.C
@@ -43,6 +43,62 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+Foam::scalar Foam::diffusivityModel::mixAvgD(const label k, const UList<scalar>& Di, const UList<scalar>& X, const UList<scalar>& Y)
+{
+    scalarField XD((X)/Di);
+    scalarField YD((Y)/Di);
+
+    scalar XY = X[k] / (1. - Y[k]);
+
+    XD[k] = 0.0;
+    YD[k] = 0.0;
+    
+    return 1. / (sum(XD) + XY*sum(YD));
+}
+
+
+Foam::scalar Foam::diffusivityModel::mixAvgMu(const UList<scalar>& muI, const UList<scalar>& X, const UList<scalar>& W)
+{
+    const label K = muI.size();
+
+    const scalar a = 1. / sqrt(8.);
+
+    scalarSquareMatrix Phi(K);
+
+    scalarField den(K);
+
+    for (label k = 0; k < K; k++)
+    {
+        for (label j = 0; j < K; j++)
+        {
+            const scalar Wj = W[j];
+            const scalar Wk = W[k];
+
+            const scalar muj = muI[j];
+            const scalar muk = muI[k];
+
+            Phi(k,j) = a * sqrt(1. + Wk/Wj) * sqr(1. + sqrt(muk/muj)*pow(Wj/Wk, 1./4.));
+        }
+    }
+
+
+    for (label k = 0; k < K; k++)
+    {
+        den[k] = sum(X * UList<scalar>(Phi[k], Phi.n()));
+    }
+
+    return sum(X * muI / den);
+
+}
+
+
+Foam::scalar Foam::diffusivityModel::mixAvgK(const UList<scalar>& k, const UList<scalar>& X)
+{
+    scalar sum1 = sum(k*X);
+    scalar sum2 = 1.0 / sum(X/k);
+    return (sum1 + sum2) / 2.0;
+}
+
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
@@ -190,6 +246,40 @@ Foam::diffusivityModel::diffusivityModel(const psiReactionThermo& thermo)
         }
     }
 
+    mu_.set
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "mu",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimensionedScalar("zero", dimArea/dimTime, 0.0)
+        )
+    );
+
+    k_.set
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "lambda",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimensionedScalar("zero", dimArea/dimTime, 0.0)
+        )
+    );
+
     forAll(species_, i)
     {
         D_.set
@@ -240,12 +330,10 @@ void Foam::diffusivityModel::correct()
 {
     const volScalarField &T = thermo_.T();
     const volScalarField &p = thermo_.p();
+    const volScalarField rho(thermo_.rho());
 
     const speciesTable &species_(thermo_.composition().species());
 
-    scalarSymmetricSquareMatrix Dij(species_.size());
-
-    const volScalarField Cv(thermo_.Cv());
     const volScalarField Wbar(thermo_.composition().W());
 
     const PtrList<volScalarField> &Y(thermo_.composition().Y());
@@ -256,17 +344,29 @@ void Foam::diffusivityModel::correct()
         Wpure[i] = thermo_.composition().W(i);
     }
 
+    scalarSymmetricSquareMatrix Dij(species_.size());
+
+    scalarField localX(species_.size());
+    scalarField localY(species_.size());
+    scalarField localCv(species_.size());
+
+    scalarField localD(species_.size());
+
+    scalarField Dii(species_.size());
+    scalarField muI(species_.size());
+    scalarField kI(species_.size());
+
     forAll (p, celli)
     {
+        const scalar rhoi = rho[celli];
         const scalar pi = p[celli];
         const scalar Ti = T[celli];
-        const scalar Cvi = Cv[celli];
         const scalar WbarI = Wbar[celli];
 
-        scalarField localX(species_.size());
-        scalarField localY(species_.size());
-
-        scalarField localD(species_.size());
+        forAll (species_, i)
+        {
+            localCv[i] = thermo_.composition().Cv(i,pi,Ti);
+        }
 
         forAll (species_, i)
         {
@@ -278,25 +378,75 @@ void Foam::diffusivityModel::correct()
 
         forAll (species_, i)
         {
-            for (label j = i; j < species_.size(); j++)
+            muI[i] = nns_[idx].mu(pi,Ti);
+
+            Dii[i] = nns_[idx].D(pi,Ti);
+            Dij(i,i) = Dii[i];
+            idx++;
+
+            for (label j = i+1; j < species_.size(); j++)
             {
                 // Calculate Dij
-                // Dij[idx] = interactions[idx].D(pi,Ti);
                 Dij(i,j) = nns_[idx].D(pi,Ti);
                 Dij(j,i) = Dij(i,j);
                 idx++;
             }
-
-            // DMat = 0;
-            //
         }
 
         forAll (species_, i)
         {
-            // localD[i] = mixAvgD(Dij.row(i));
+            const scalar R = Neutral::RR / Wpure[i];
+            const scalar CvTrans = (3./2.)*R;
+            const scalar CvRot = (neutrals_[i].F()/2.)*R;
+            const scalar CvVib = localCv[i] - CvTrans - CvRot;
+
+            const scalar rSc = rhoi * Dii[i] / muI[i];
+
+            const scalar A = 5./2. - rSc;
+            const scalar Zrot = neutrals_[i].Zrot(Ti);
+            const scalar B = Zrot + (2./Neutral::pi) * ((5./3.)*CvRot + rSc);
+            const scalar AB = (2./Neutral::pi)*(A/B);
+
+            const scalar fTrans = (5./2.) * (1.0 - AB * CvRot / CvTrans);
+            const scalar fRot = rSc*(1.0+AB);
+            const scalar fVib = rSc;
+
+            kI[i] = muI[i]*(fTrans*CvTrans + fRot*CvRot + fVib*CvVib);
+
         }
 
-        Info << localY << endl;
+        Switch pure = false;
+        label pureSpecieI = -1;
+        forAll (species_, i)
+        {
+            if (1. - localX[i] < 1e-12)
+            {
+                pureSpecieI = i;
+                pure = true;
+                break;
+            }
+        }
+
+        if (pure)
+        {
+            UList<scalar> Dpure(Dij[pureSpecieI], Dij.n());
+
+            forAll (species_, i)
+            {
+                D_[i][celli] = Dpure[i];
+            }
+        }
+        else
+        {
+            forAll (species_, i)
+            {
+                UList<scalar> Di(Dij[i], Dij.n());
+                D_[i][celli] = mixAvgD(i, Di, localX, localY);
+            }
+        }
+
+        mu_()[celli] = mixAvgMu(muI, localX, Wpure);
+        k_()[celli] = mixAvgK(kI, localX);
     }
 
     return;
diff --git a/diffusivityModel/diffusivityModel/diffusivityModel.H b/diffusivityModel/diffusivityModel/diffusivityModel.H
index f7fc2fa..67e8bfa 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModel.H
+++ b/diffusivityModel/diffusivityModel/diffusivityModel.H
@@ -67,6 +67,12 @@ class diffusivityModel
         //- Mixture averaged diffusivity fields
         PtrList<volScalarField> D_;
 
+        //- Mixture averaged viscosity field
+        autoPtr<volScalarField> mu_;
+
+        //- Mixture averaged thermal conductivity field
+        autoPtr<volScalarField> k_;
+
         //- Electron object
         // autoPtr<Electron> electron_;
 
@@ -88,6 +94,15 @@ class diffusivityModel
         //- Disallow default bitwise assignment
         void operator=(const diffusivityModel&);
 
+        //- Disallow default bitwise assignment
+        scalar mixAvgD(const label k, const UList<scalar>& Di, const UList<scalar>& X, const UList<scalar>& Y);
+
+        //- Disallow default bitwise assignment
+        scalar mixAvgMu(const UList<scalar>& muI, const UList<scalar>& X, const UList<scalar>& W);
+
+        //- Disallow default bitwise assignment
+        scalar mixAvgK(const UList<scalar>& k, const UList<scalar>& X);
+
 
 public:
 
diff --git a/eReactingFoam.C b/eReactingFoam.C
index 70f506f..584372d 100644
--- a/eReactingFoam.C
+++ b/eReactingFoam.C
@@ -91,6 +91,8 @@ int main(int argc, char *argv[])
 
         while (pimple.loop())
         {
+            diff.correct();
+
             #include "UEqn.H"
             #include "YEqn.H"
             #include "EEqn.H"