diff --git a/Make/files b/Make/files
index bfc605f..0bc792e 100644
--- a/Make/files
+++ b/Make/files
@@ -2,8 +2,10 @@ diffusivityModel/Particle/Particle.C
 diffusivityModel/Neutral/Neutral.C
 diffusivityModel/Ion/Ion.C
 
+diffusivityModel/Interaction/Interaction.C
 diffusivityModel/Stockmayer/Stockmayer.C
 diffusivityModel/Coulomb/Coulomb.C
+diffusivityModel/N64/N64.C
 
 diffusivityModel/diffusivityModel/diffusivityModel.C
 
diff --git a/Make/options b/Make/options
index e6da545..2df0871 100644
--- a/Make/options
+++ b/Make/options
@@ -3,8 +3,10 @@ EXE_INC = \
     -IdiffusivityModel/Particle \
     -IdiffusivityModel/Neutral \
     -IdiffusivityModel/Ion \
+    -IdiffusivityModel/Interaction \
     -IdiffusivityModel/Stockmayer \
     -IdiffusivityModel/Coulomb \
+    -IdiffusivityModel/N64 \
     -IdiffusivityModel/diffusivityModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
diff --git a/diffusivityModel/Coulomb/Coulomb.H b/diffusivityModel/Coulomb/Coulomb.H
index 262b2c5..0e9792a 100644
--- a/diffusivityModel/Coulomb/Coulomb.H
+++ b/diffusivityModel/Coulomb/Coulomb.H
@@ -39,6 +39,7 @@ SourceFiles
 #include "scalar.H"
 #include "scalarField.H"
 #include "Ion.H"
+#include "Interaction.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -60,6 +61,7 @@ Ostream& operator<<(Ostream&, const Coulomb&);
 \*---------------------------------------------------------------------------*/
 
 class Coulomb
+    : public Interaction
 {
     // Private data
 
@@ -84,15 +86,6 @@ class Coulomb
         //- Disallow default bitwise assignment
         void operator=(const Coulomb&);
 
-        //- Disallow default bitwise assignment
-        inline scalar power(const scalar a, const scalar b);
-
-        //- Disallow default bitwise assignment
-        inline scalar expon(const scalar a);
-
-        //- Disallow default bitwise assignment
-        inline scalar loge(const scalar a);
-
 
 public:
 
diff --git a/diffusivityModel/Coulomb/CoulombI.H b/diffusivityModel/Coulomb/CoulombI.H
index f3dc1a2..fa842cb 100644
--- a/diffusivityModel/Coulomb/CoulombI.H
+++ b/diffusivityModel/Coulomb/CoulombI.H
@@ -27,22 +27,6 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-inline Foam::scalar Foam::Coulomb::power(const scalar a, const scalar b)
-{
-    return pow(a, b);
-}
-
-
-inline Foam::scalar Foam::Coulomb::expon(const scalar a)
-{
-    return exp(a);
-}
-
-
-inline Foam::scalar Foam::Coulomb::loge(const scalar a)
-{
-    return log(a);
-}
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
diff --git a/diffusivityModel/Interaction/Interaction.C b/diffusivityModel/Interaction/Interaction.C
new file mode 100644
index 0000000..2c047f2
--- /dev/null
+++ b/diffusivityModel/Interaction/Interaction.C
@@ -0,0 +1,77 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     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 "Interaction.H"
+
+#include "error.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::Interaction::Interaction()
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::Interaction::~Interaction()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+
+
+// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
+
+void Foam::Interaction::operator=(const Interaction& rhs)
+{
+    // Check for assignment to self
+    if (this == &rhs)
+    {
+        FatalErrorInFunction
+            << "Attempted assignment to self"
+            << abort(FatalError);
+    }
+}
+
+// * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //
diff --git a/diffusivityModel/Interaction/Interaction.H b/diffusivityModel/Interaction/Interaction.H
new file mode 100644
index 0000000..a1f7761
--- /dev/null
+++ b/diffusivityModel/Interaction/Interaction.H
@@ -0,0 +1,144 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::Interaction
+
+Description
+
+SourceFiles
+    InteractionI.H
+    Interaction.C
+    InteractionIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Interaction_H
+#define Interaction_H
+
+#include "scalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class Istream;
+class Ostream;
+
+// Forward declaration of friend functions and operators
+class Interaction;
+Istream& operator>>(Istream&, Interaction&);
+Ostream& operator<<(Ostream&, const Interaction&);
+
+
+/*---------------------------------------------------------------------------*\
+                         Class Interaction Declaration
+\*---------------------------------------------------------------------------*/
+
+class Interaction
+{
+    // Private data
+
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        Interaction(const Interaction&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const Interaction&);
+
+protected:
+
+        //- Disallow default bitwise assignment
+        inline scalar power(const scalar a, const scalar b);
+
+        //- Disallow default bitwise assignment
+        inline scalar expon(const scalar a);
+
+        //- Disallow default bitwise assignment
+        inline scalar loge(const scalar a);
+
+
+public:
+
+    // Static data members
+
+
+    // Constructors
+
+        //- Construct null
+        Interaction();
+
+        //- Construct from Istream
+        Interaction(Istream&);
+
+        //- Construct as copy
+        // Interaction(const Interaction&);
+
+
+    //- Destructor
+    ~Interaction();
+
+
+    // Member Functions
+
+        // Access
+
+        // Check
+
+        // Edit
+
+        // Write
+
+
+    // Member Operators
+
+
+    // Friend Functions
+
+    // Friend Operators
+
+    // IOstream Operators
+
+        friend Istream& operator>>(Istream&, Interaction&);
+        friend Ostream& operator<<(Ostream&, const Interaction&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "InteractionI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/diffusivityModel/Interaction/InteractionI.H b/diffusivityModel/Interaction/InteractionI.H
new file mode 100644
index 0000000..e7b2406
--- /dev/null
+++ b/diffusivityModel/Interaction/InteractionI.H
@@ -0,0 +1,74 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+inline Foam::scalar Foam::Interaction::power(const scalar a, const scalar b)
+{
+    return pow(a, b);
+}
+
+
+inline Foam::scalar Foam::Interaction::expon(const scalar a)
+{
+    return exp(a);
+}
+
+
+inline Foam::scalar Foam::Interaction::loge(const scalar a)
+{
+    return log(a);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //
diff --git a/diffusivityModel/Ion/Ion.C b/diffusivityModel/Ion/Ion.C
index 487e390..0dd9bbc 100644
--- a/diffusivityModel/Ion/Ion.C
+++ b/diffusivityModel/Ion/Ion.C
@@ -48,13 +48,29 @@ Foam::Ion::Ion(const dictionary& dict)
 :
     Particle(dict),
     alpha_(dict.lookupOrDefault("dipolePolarizability", 0.0)),
-    C6_(dict.lookupOrDefault("dispersionCoef", 0.0)),
+    C6_(0.0),
     Zrot_(dict.lookupOrDefault("rotationalRelaxation", 0.0)),
     n_(dict.lookupOrDefault("nStockmayer", 12)),
     rctTarget_(dict.lookupOrDefault("rctTargets", hashedWordList())),
     Arct_(dict.lookupOrDefault("Arct", scalarList())),
     Brct_(dict.lookupOrDefault("Brct", scalarList()))
-{}
+{
+    if (z() > 0)
+    {
+        C6_ = exp(1.8853 * log(alpha_) + 0.2682); // * sqr(e);
+    }
+    else if (z() < 0)
+    {
+        C6_ = exp(3.2246 * log(alpha_) - 3.2397); // * sqr(e);
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Ion must have nonzero charge number"
+            << abort(FatalError);
+    }
+
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
diff --git a/diffusivityModel/N64/N64.C b/diffusivityModel/N64/N64.C
new file mode 100644
index 0000000..dfeec40
--- /dev/null
+++ b/diffusivityModel/N64/N64.C
@@ -0,0 +1,111 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     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 "N64.H"
+
+#include "error.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+const Foam::scalar Foam::N64::a11[(M+1)*(M+2)/2] = {
+        0., -0.81508875, -0.06484549, 0.31926761, 0.35164219, -0.01484638,
+        -0.02385399, -0.14337633, -0.18493828, 0.05466979, -0.0010004, 0.01936938,
+        -0.0130214, 0.12179535, -0.05237791, -0.00172218, -0.00447904, 0.02211965,
+        0.00964791, -0.07566612, 0.02527412
+    };
+
+const Foam::scalar Foam::N64::a22[(M+1)*(M+2)/2] = {
+        0.00000000e+00, -8.59231154e-01, -7.02815662e-02, 2.79760234e-01,
+        3.59026813e-01, -3.25082644e-01, -8.27621218e-03, -1.12825149e-01,
+        -2.07675933e-01, 1.02433682e+00, -4.59772292e-04, 1.03264160e-02,
+        -1.35152522e-02, 1.69188798e-01, -1.21839345e+00, -1.88957374e-03,
+        -4.45369981e-03, 2.07928092e-02, 7.40762555e-03, -9.71847966e-02,
+        5.07723424e-01
+    };
+
+Foam::scalar Foam::N64::vec[(M+1)*(M+2)/2] = { 0.0 };
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::N64::N64(const Ion& a, const Neutral& b)
+:
+    a_(a),
+    b_(b),
+    ksiij_(a_.alpha()/(sqr(a_.z())*sqrt(b_.alpha())*(1.0+pow(2.0*a_.alpha()/b_.alpha(),2./3.)))),
+    sigmaij_(K1 * (pow(a_.alpha(), 1./3.) + pow(b_.alpha(), 1./3.)) / pow(a_.alpha()*b_.alpha()*(1. + 1./ksiij_), kappa)),
+    eij_(5.2 * b_.alpha() * sqr(a_.z()) * (1.0 + ksiij_) / pow4(sigmaij_)),
+    C6ij_(0.0),
+    gammaij_(0.0),
+    Wij_(1.0/(1.0/a_.W() + 1.0/b_.W()))
+{
+    scalar aij = a_.alpha() / b_.alpha();
+    scalar aji = 1.0 / aij;
+    C6ij_ = 2.0 * a_.C6() * b_.C6() / (aji * a_.C6() + aij * b_.C6());
+    gammaij_ = ((2.0/sqr(a_.z())) * C6ij_ + b_.alphaQ()) / (b_.alpha() * sqr(sigmaij_));
+
+    Info << "N64, gamma = " << gammaij_ << endl;
+
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::N64::~N64()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+
+
+// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
+
+void Foam::N64::operator=(const N64& rhs)
+{
+    // Check for assignment to self
+    if (this == &rhs)
+    {
+        FatalErrorInFunction
+            << "Attempted assignment to self"
+            << abort(FatalError);
+    }
+}
+
+// * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //
diff --git a/diffusivityModel/N64/N64.H b/diffusivityModel/N64/N64.H
new file mode 100644
index 0000000..35da13c
--- /dev/null
+++ b/diffusivityModel/N64/N64.H
@@ -0,0 +1,258 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::N64
+
+Description
+
+SourceFiles
+    N64I.H
+    N64.C
+    N64IO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef N64_H
+#define N64_H
+
+#include "scalar.H"
+#include "scalarField.H"
+#include "Ion.H"
+#include "Neutral.H"
+#include "Interaction.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class Istream;
+class Ostream;
+
+// Forward declaration of friend functions and operators
+class N64;
+Istream& operator>>(Istream&, N64&);
+Ostream& operator<<(Ostream&, const N64&);
+
+
+/*---------------------------------------------------------------------------*\
+                         Class N64 Declaration
+\*---------------------------------------------------------------------------*/
+
+class N64
+    : public Interaction
+{
+    // Private data
+    //
+        //- Highest polynomial order of collision integral fitting formular
+        static constexpr label M = 5;
+
+        //- Collision participant ion a
+        const Ion &a_;
+
+        //- Collision participant neutral b
+        const Neutral &b_;
+
+        //- Combined well depth parameter
+        scalar ksiij_;
+
+        //- Combined collision diameter [A] (Angstrom)
+        scalar sigmaij_;
+
+        //- Combined well depth [eV]
+        scalar eij_;
+
+        //- Combined dispersion coefficient C6/(e^2) [A^5]
+        scalar C6ij_;
+
+        //- Collision integral parameter 2
+        scalar gammaij_;
+
+        //- Combined molecular weight
+        scalar Wij_;
+
+
+        //- Construct null
+        static const scalar a11[(M+1)*(M+2)/2];
+        static const scalar a22[(M+1)*(M+2)/2];
+        static scalar vec[(M+1)*(M+2)/2];
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        N64(const N64&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const N64&);
+
+        // primary template
+        template <label Order>
+            class BivarPolyFeat
+            {
+                public:
+                    static scalar *transform (scalar x1, scalar x2, scalar *vec)
+                    {
+                        scalar *src = BivarPolyFeat<Order-1>::transform(x1,x2,vec);
+
+                        for (label i = 0; i < Order; i++)
+                        {
+                            *(src+Order+i) = *(src + i) * x1;
+                        }
+
+                        *(src+Order+Order) = *(src + (Order - 1)) * x2;
+
+                        return src + Order;
+                    }
+            };
+
+        // convenience function
+        inline void transform (scalar x1, scalar x2, scalar *vec)
+        {
+            BivarPolyFeat<M>::transform(x1,x2,vec);
+        }
+
+        template <label N>
+            struct ScalarProduct
+            {
+                static scalar product (const scalar* first, const scalar* second)
+                {
+                    return ScalarProduct<N-1>::product(first + 1, second + 1)
+                        + *first * *second;
+                }
+            };
+
+        // convenience function
+        inline scalar scalarProduct (const scalar *x1, const scalar *x2)
+        {
+            return ScalarProduct<(M+1)*(M+2)/2>::product(x1,x2);
+        }
+
+
+public:
+
+    // Static data members
+
+        //- Construct null
+        static constexpr scalar K1 = 1.767;
+
+        //- Construct null
+        static constexpr scalar K2 = 0.720;
+
+        //- Construct null
+        static constexpr scalar kappa = 0.0095;
+
+
+    // Constructors
+
+        //- Construct null
+        N64();
+
+        //- Construct from components
+        N64(const Ion& a, const Neutral& b);
+
+        //- Construct from Istream
+        N64(Istream&);
+
+        //- Construct as copy
+        // N64(const N64&);
+
+
+    //- Destructor
+    ~N64();
+
+
+    // Member Functions
+
+        // Access
+
+        // Check
+
+        // Edit
+
+        // Write
+
+        inline scalar Omega11(const scalar T);
+
+
+        inline scalar Omega22(const scalar T);
+
+
+        inline scalar D(const scalar p, const scalar T);
+
+
+        inline scalar mu(const scalar p, const scalar T);
+
+
+        inline scalar k(const scalar p, const scalar T);
+
+
+    // Member Operators
+
+
+    // Friend Functions
+
+    // Friend Operators
+
+    // IOstream Operators
+
+        friend Istream& operator>>(Istream&, N64&);
+        friend Ostream& operator<<(Ostream&, const N64&);
+};
+
+// partial specialization as end criteria
+template <>
+class N64::BivarPolyFeat<0>
+{
+    public:
+        static scalar *transform (scalar x1, scalar x2, scalar *vec)
+        {
+            *vec = 1.0;
+            return vec;
+        }
+};
+
+template <>
+struct N64::ScalarProduct<0>
+{
+    static scalar product (const scalar* first, const scalar* second) {
+        return 0.0;
+    };
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "N64I.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/diffusivityModel/N64/N64I.H b/diffusivityModel/N64/N64I.H
new file mode 100644
index 0000000..f808339
--- /dev/null
+++ b/diffusivityModel/N64/N64I.H
@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+/*
+template <Foam::label I>
+inline Foam::scalar Foam::N64::bivarPolyFeat(const scalar Tstar)
+{
+    const scalar a[21] = {
+        0., -0.81508875, -0.06484549, 0.31926761, 0.35164219, -0.01484638,
+        -0.02385399, -0.14337633, -0.18493828, 0.05466979, -0.0010004, 0.01936938,
+        -0.0130214, 0.12179535, -0.05237791, -0.00172218, -0.00447904, 0.02211965,
+        0.00964791, -0.07566612, 0.02527412
+    };
+
+    const scalar p[21] = { 0. };
+
+    const scalar A = 0.5;
+    const scalar B = -0.129;
+    return (A * loge (Tstar) + B) / sqr(Tstar);
+}
+*/
+
+
+inline Foam::scalar Foam::N64::Omega11(const scalar Tstar)
+{
+    transform (loge(Tstar), gammaij_, vec);
+    return scalarProduct (a11, vec);
+}
+
+
+inline Foam::scalar Foam::N64::Omega22(const scalar Tstar)
+{
+    transform (loge(Tstar), gammaij_, vec);
+    return scalarProduct (a22, vec);
+}
+
+
+inline Foam::scalar Foam::N64::D(const scalar p, const scalar T)
+{
+    return (3.0/8.0) * sqrt(2.0*a_.pi*a_.NA*pow3(a_.k*T)/Wij_)
+        / (a_.pi * sqr(sigmaij_) * p * Omega11(sigmaij_*T/eij_));
+}
+
+
+inline Foam::scalar Foam::N64::mu(const scalar p, const scalar T)
+{
+    return (5.0/16.0) * sqrt(2.0*Wij_*a_.k*T/(a_.pi*a_.NA))
+        / (sqr(sigmaij_) * Omega22(sigmaij_*T/eij_));
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //
diff --git a/diffusivityModel/Neutral/Neutral.C b/diffusivityModel/Neutral/Neutral.C
index d0d927c..a2bb4eb 100644
--- a/diffusivityModel/Neutral/Neutral.C
+++ b/diffusivityModel/Neutral/Neutral.C
@@ -51,7 +51,7 @@ Foam::Neutral::Neutral(const dictionary& dict)
     diameter_(dict.lookupOrDefault("diameter", 0.0)),
     dipoleMoment_(dict.lookupOrDefault("dipoleMoment", 0.0)),
     alpha_(dict.lookupOrDefault("dipolePolarizability", 0.0)),
-    alphaQ_(dict.lookupOrDefault("quadpolePolarizability", 0.0)),
+    alphaQ_(0.0),
     C6_(dict.lookupOrDefault("dispersionCoef", 0.0)),
     Zrot_(dict.lookupOrDefault("rotationalRelaxation", 0.0))
 {
@@ -66,10 +66,27 @@ Foam::Neutral::Neutral(const dictionary& dict)
         diameter_ = tranlib[2];
         dipoleMoment_ = tranlib[3];
         alpha_ = tranlib[4];
-        // alphaQ_ = tranlib[0];
-        // C6_ = tranlib[0];
         Zrot_ = tranlib[5];
+
+        dict.readIfPresent("wellDepth", wellDepth_);
+        dict.readIfPresent("diameter", diameter_);
+        dict.readIfPresent("dipoleMoment", dipoleMoment_);
+        dict.readIfPresent("dipolePolarizability", alpha_);
+        dict.readIfPresent("dispersionCoef", C6_);
+        dict.readIfPresent("rotationalRelaxation", Zrot_);
     }
+    else
+    {
+        wellDepth_ = readScalar(dict.lookup("wellDepth"));
+        diameter_ = readScalar(dict.lookup("diameter"));
+        dipoleMoment_ = readScalar(dict.lookup("dipoleMoment"));
+        alpha_ = readScalar(dict.lookup("dipolePolarizability"));
+        C6_ = readScalar(dict.lookup("dispersionCoef"));
+        Zrot_ = readScalar(dict.lookup("rotationalRelaxation"));
+    }
+
+    C6_ = exp(1.8846 * log(alpha_) - 0.4737); // * sqr(e);
+    alphaQ_ = 2 * C6_;
 
 }
 
diff --git a/diffusivityModel/Stockmayer/Stockmayer.H b/diffusivityModel/Stockmayer/Stockmayer.H
index c01a28b..4a1e970 100644
--- a/diffusivityModel/Stockmayer/Stockmayer.H
+++ b/diffusivityModel/Stockmayer/Stockmayer.H
@@ -39,6 +39,7 @@ SourceFiles
 #include "scalar.H"
 #include "scalarField.H"
 #include "Neutral.H"
+#include "Interaction.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -60,6 +61,7 @@ Ostream& operator<<(Ostream&, const Stockmayer&);
 \*---------------------------------------------------------------------------*/
 
 class Stockmayer
+    : public Interaction
 {
     // Private data
 
@@ -90,12 +92,6 @@ class Stockmayer
         //- Disallow default bitwise assignment
         void operator=(const Stockmayer&);
 
-        //- Disallow default bitwise assignment
-        inline scalar power(const scalar a, const scalar b);
-
-        //- Disallow default bitwise assignment
-        inline scalar expon(const scalar a);
-
 
 public:
 
diff --git a/diffusivityModel/Stockmayer/StockmayerI.H b/diffusivityModel/Stockmayer/StockmayerI.H
index 7e0cb2f..af33882 100644
--- a/diffusivityModel/Stockmayer/StockmayerI.H
+++ b/diffusivityModel/Stockmayer/StockmayerI.H
@@ -27,15 +27,6 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-inline Foam::scalar Foam::Stockmayer::power(const scalar a, const scalar b)
-{
-    return pow(a, b);
-}
-
-inline Foam::scalar Foam::Stockmayer::expon(const scalar a)
-{
-    return exp(a);
-}
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
diff --git a/diffusivityModel/diffusivityModel/diffusivityModel.C b/diffusivityModel/diffusivityModel/diffusivityModel.C
index cb0b7b8..41f3afe 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModel.C
+++ b/diffusivityModel/diffusivityModel/diffusivityModel.C
@@ -34,6 +34,7 @@ License
 #include "Ion.H"
 #include "Stockmayer.H"
 #include "Coulomb.H"
+#include "N64.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -114,7 +115,8 @@ Foam::diffusivityModel::diffusivityModel(const psiReactionThermo& thermo)
     neutrals_(thermo_.composition().species().size()),
     ions_(thermo_.composition().species().size()),
     nns_(thermo_.composition().species().size()*thermo_.composition().species().size()),
-    ccs_(thermo_.composition().species().size()*thermo_.composition().species().size())
+    ccs_(thermo_.composition().species().size()*thermo_.composition().species().size()),
+    cns_(thermo_.composition().species().size()*thermo_.composition().species().size())
 {
 
     const speciesTable &species_(thermo_.composition().species());
@@ -197,8 +199,8 @@ Foam::diffusivityModel::diffusivityModel(const psiReactionThermo& thermo)
     neutrals_.resize(nNeutral);
     nns_.resize(nNeutral*(nNeutral+1)/2);
 
-    // cns_.resize(nNeutral*nIon);
     ccs_.resize((nIon+1)*(nIon+2)/2 - 1);
+    cns_.resize(nNeutral*nIon);
     // ens_.resize(nNeutral);
     // ee_;
 
@@ -260,7 +262,6 @@ Foam::diffusivityModel::diffusivityModel(const psiReactionThermo& thermo)
         }
     }
 
-    /*
     iPair = 0;
     forAll(ions_, i)
     {
@@ -270,6 +271,7 @@ Foam::diffusivityModel::diffusivityModel(const psiReactionThermo& thermo)
             iPair++;
         }
     }
+    /*
     */
 
     iPair = 0;
@@ -498,8 +500,8 @@ void Foam::diffusivityModel::correct()
             {
                 label j = J + 1;
                 // Calculate Dij
-                // Dij(i,j) = nns_[idx1].D(pi,Ti);
-                Dij(i,j) = SMALL;
+                Dij(i,j) = ccs_[idx1].D(pi,Ti,rhoQc2i);
+                // Dij(i,j) = SMALL;
                 Dij(j,i) = Dij(i,j);
                 idx1++;
             }
@@ -508,8 +510,8 @@ void Foam::diffusivityModel::correct()
             for (label j = ions_.size()+1; j < species_.size(); j++)
             {
                 // Calculate Dij
-                // Dij(i,j) = nns_[idx2].D(pi,Ti);
-                Dij(i,j) = SMALL;
+                Dij(i,j) = cns_[idx2].D(pi,Ti);
+                // Dij(i,j) = SMALL;
                 Dij(j,i) = Dij(i,j);
                 idx2++;
             }
diff --git a/diffusivityModel/diffusivityModel/diffusivityModel.H b/diffusivityModel/diffusivityModel/diffusivityModel.H
index 6337762..9c94492 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModel.H
+++ b/diffusivityModel/diffusivityModel/diffusivityModel.H
@@ -50,6 +50,7 @@ class Neutral;
 class Ion;
 class Stockmayer;
 class Coulomb;
+class N64;
 
 // Forward declaration of friend functions and operators
 class diffusivityModel;
@@ -90,6 +91,9 @@ class diffusivityModel
         //- Stockmayer potential list
         PtrList<Coulomb> ccs_;
 
+        //- Stockmayer potential list
+        PtrList<N64> cns_;
+
 
     // Private Member Functions