diff --git a/diffusivityModel/Interaction/InteractionI.H b/diffusivityModel/Interaction/InteractionI.H
index 97cc6fb..3726675 100644
--- a/diffusivityModel/Interaction/InteractionI.H
+++ b/diffusivityModel/Interaction/InteractionI.H
@@ -53,7 +53,7 @@ inline Foam::scalar Foam::Interaction::D(const scalar mij, const scalar p, const
 
 inline Foam::scalar Foam::Interaction::mu(const scalar mk, const scalar T, const scalar OmegaD)
 {
-    return (5.0/16.0) * sqrt(Foam::Particle::pi*mk*Foam::Particle::k*T) / OmegaD;
+    return (5.0/16.0) * sqrt(Foam::Particle::pi*mk*Foam::Particle::k*T/Foam::Particle::NA) / OmegaD;
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
diff --git a/diffusivityModel/Particle/Particle.C b/diffusivityModel/Particle/Particle.C
index 5362536..be1ece6 100644
--- a/diffusivityModel/Particle/Particle.C
+++ b/diffusivityModel/Particle/Particle.C
@@ -49,8 +49,8 @@ const Foam::scalar Foam::Particle::k = constant::physicoChemical::k.value();
 //- 1 Angstom in meter
 const Foam::scalar Foam::Particle::Angstrom = 1e-10;
 
-//- 1 Debye in [ Coulomb x meter ]
-const Foam::scalar Foam::Particle::Debye = 1.0 / 2.99792458e29;
+//- 1 Debye in [ sqrt( J x meter ^ 3 ) ]
+const Foam::scalar Foam::Particle::Debye = 1e-18 * Foam::pow(1./100., 5./2) * Foam::pow(1./1000., 1./2.);
 
 //- Circumference / Diameter
 const Foam::scalar Foam::Particle::pi = constant::mathematical::pi;
diff --git a/diffusivityModel/Stockmayer/Stockmayer.C b/diffusivityModel/Stockmayer/Stockmayer.C
index 93b5838..d2ba155 100644
--- a/diffusivityModel/Stockmayer/Stockmayer.C
+++ b/diffusivityModel/Stockmayer/Stockmayer.C
@@ -68,15 +68,17 @@ Foam::Stockmayer::Stockmayer(const Neutral& a, const Neutral& b)
     Wij_(1.0/(1.0/a_.W() + 1.0/b_.W())),
     deltaij_(a_.dipoleMoment()*b_.dipoleMoment()/2.0/(eij_*a_.k)/pow3(sigmaij_))
 {
-    if (b_.dipoleMoment() > 0.0 && deltaij_ > 0.0)
+    if (b_.dipoleMoment() > 0.0 && deltaij_ == 0.0)
     {
-        scalar alphaStar = a_.alpha() / a_.diameter() / a_.diameter() / a_.diameter();
 
-        scalar muStar = b_.dipoleMoment() / Foam::sqrt(b_.wellDepth()*b_.diameter()*b_.diameter()*b_.diameter());
+        const scalar alphaStar = a_.alpha() / pow3(a_.diameter()/a.Angstrom);
 
-        scalar ksi = 1.0 + alphaStar * muStar * Foam::sqrt(b_.wellDepth()/a_.wellDepth());
+        const scalar muStar = b_.dipoleMoment() / Foam::sqrt(b_.wellDepth()*b_.k*pow3(b_.diameter()));
+
+        const scalar ksi = 1.0 + alphaStar * muStar * Foam::sqrt(b_.wellDepth()/a_.wellDepth()) / 4.0;
+
+        eij_ *= sqr(ksi);
 
-        eij_ *= ksi * ksi;
         sigmaij_ *= Foam::pow(ksi, -1.0/6.0);
     }
 }
diff --git a/diffusivityModel/Stockmayer/StockmayerI.H b/diffusivityModel/Stockmayer/StockmayerI.H
index 8f3f5c3..d40b93f 100644
--- a/diffusivityModel/Stockmayer/StockmayerI.H
+++ b/diffusivityModel/Stockmayer/StockmayerI.H
@@ -40,14 +40,14 @@ License
 inline Foam::scalar Foam::Stockmayer::Omega11(const scalar Tstar)
 {
     constexpr scalar a[6] = { 1.0548, 0.15504, 0.55909, 2.1705, 0.093193, 1.5 };
-    return a[0] * power(Tstar, -a[1]) + power(Tstar + a[2], -a[3]);
+    return (a[0] * power(Tstar, -a[1]) + power(Tstar + a[2], -a[3]))*f11(Tstar);
 }
 
 
 inline Foam::scalar Foam::Stockmayer::Omega22(const scalar Tstar)
 {
     constexpr scalar b[6] = { 1.0413, 0.11930, 0.43628, 1.6041, 0.095661, 2.0 };
-    return b[0] * power(Tstar, -b[1]) + power(Tstar + b[2], -b[3]);
+    return (b[0] * power(Tstar, -b[1]) + power(Tstar + b[2], -b[3]))*f22(Tstar);
 }
 
 
diff --git a/diffusivityModel/diffusivityModel/diffusivityModel.C b/diffusivityModel/diffusivityModel/diffusivityModel.C
index 658afde..43f45c3 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModel.C
+++ b/diffusivityModel/diffusivityModel/diffusivityModel.C
@@ -139,7 +139,7 @@ Foam::scalar Foam::diffusivityModel::mixAvgMu(const UList<scalar>& muI, const UL
             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.));
+            Phi(k,j) = a * sqr(1. + sqrt(muk/muj)*pow(Wj/Wk, 1./4.)) / sqrt(1. + Wk/Wj);
         }
     }