diff --git a/applications/solvers/cfdemSolverRhoPimpleChem/YEqn.H b/applications/solvers/cfdemSolverRhoPimpleChem/YEqn.H
index cc7f380a..e892eb48 100644
--- a/applications/solvers/cfdemSolverRhoPimpleChem/YEqn.H
+++ b/applications/solvers/cfdemSolverRhoPimpleChem/YEqn.H
@@ -42,7 +42,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
             fvOptions.correct(Yi);
 
-            #include "debugYEqn.H"
+            // #include "debugYEqn.H"
 
             Yi.max(0.0);
             Yt += Yi;
diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.C b/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.C
index 2ed747d9..156de95a 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.C
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.C
@@ -173,12 +173,10 @@ void diffusionCoefficient::execute()
     scalar Pfluid(0);
     scalar molarConcfluid(0);
     scalar Texp(0);
-    List<scalar> dBinarytot_;
-    dBinarytot_.setSize(diffusantGasNames_.size());
-    List<scalar> molecularDiffusion_;
-    molecularDiffusion_.setSize(diffusantGasNames_.size());
-    List<scalar> Xstag_tot_;
-    Xstag_tot_.setSize(diffusantGasNames_.size());
+    List<scalar> MixtureBinaryDiffusion_;
+    MixtureBinaryDiffusion_.setSize(diffusantGasNames_.size());
+    List<scalar> TotalFraction_;
+    TotalFraction_.setSize(diffusantGasNames_.size());
 
     double  **dBinary_   = new double*[diffusantGasNames_.size()];
     double  **molNum_    = new double*[diffusantGasNames_.size()];
@@ -220,7 +218,7 @@ void diffusionCoefficient::execute()
                 }
             }
 
-            Texp            =   pow(Tfluid,1.75);
+            Texp    =   pow(Tfluid,1.75);
 
             for (int i=0; i<diffusantGasNames_.size();i++)
             {
@@ -228,9 +226,8 @@ void diffusionCoefficient::execute()
                 molNum_[i]      =   new double [speciesNames_.size()];
                 volDiff_[i]     =   new double [speciesNames_.size()];
 
-                dBinarytot_[i]  =   0.0;
-                molecularDiffusion_[i]      =   0.0;
-                Xstag_tot_[i]   =   0.0;
+                MixtureBinaryDiffusion_[i]  =   0.0;
+                TotalFraction_[i]   =   0.0;
                 for (int j=0; j < speciesNames_.size();j++)
                 {
                     // get molecular diffusion coefficients if diffusant gas and reactant gas are not equal
@@ -256,33 +253,22 @@ void diffusionCoefficient::execute()
                             Info << "Molecular diffusion for species " << diffusantGasNames_[i] << " in "
                                  << speciesNames_[j] << " is : " << dBinary_[i][j] << nl << endl;
 
-                            // sum of all binary diffusion coefficients except i = j
-                            dBinarytot_[i]  +=  dBinary_[i][j];
-                            Info << "Sum of binary mol. diffusion coefficients for diffusing gas " << diffusantGasNames_[i]
-                                 << " : " << dBinarytot_[i] << nl << endl;
-
                             Info << "Molar fraction of species " << speciesNames_[j] << " : " << Xfluid_[j] << nl << endl;
-                            Info << "What is value of xstag? " << Xstag_tot_[i] << nl << endl;
 
-                            // sum of all molar fractions except the diffusant gas
-                            Xstag_tot_[i]   +=  Xfluid_[j];
-                            Info << "Sum of stagnant gases molar fractions according to diffusing gas " << diffusantGasNames_[i]
-                                 << " : " << Xstag_tot_[i] << nl << endl;
+                            if (Xfluid_[j] != 0.)
+                            {
+                                // sum of all stagnant gases to sum of binary diffusion
+                                TotalFraction_[i]   +=  Xfluid_[j]/dBinary_[i][j];
 
-                            // dCoeff -- diffusion component of diffusant gas
-                            if (Xstag_tot_[i] == 0.0) //(Xstag_tot_[i] == 0. || (1-Xfluid_[i]) == 0.)
-                            {
-                                molecularDiffusion_[i] = 0.0;
-                            } else
-                            {
-                                molecularDiffusion_[i]  =   (1-Xfluid_[i])*dBinarytot_[i]/Xstag_tot_[i];
+                                // dCoeff -- diffusion component of diffusant gas
+                                MixtureBinaryDiffusion_[i]  =   (1.0-Xfluid_[i])/TotalFraction_[i];
                             }
 
                             Info << "Multicomp. mix diffusion for species " << diffusantGasNames_[i]
-                                 << " is: " << molecularDiffusion_[i] << nl << endl;
+                                 << " is: " << MixtureBinaryDiffusion_[i] << nl << endl;
 
                             // pass on dCoeff values to array
-                            diffusionCoefficients_[i][index][0]= molecularDiffusion_[i];
+                            diffusionCoefficients_[i][index][0]= MixtureBinaryDiffusion_[i];
 
                         }else
                         {
@@ -300,7 +286,7 @@ void diffusionCoefficient::execute()
         {
             for(int i =0; i<diffusantGasNames_.size();i++)
             {
-                Info << "effective diffusionCoefficient of species " << diffusantGasNames_[i] << " = " << diffusionCoefficients_[i][index][0] << endl;
+                Info << "diffusionCoefficient of species " << diffusantGasNames_[i] << " = " << diffusionCoefficients_[i][index][0] << endl;
             }
         }
     }
diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C b/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C
index 0c6e9aa3..2386ed47 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C
@@ -262,7 +262,7 @@ void species::execute()
                 for (int i = 0; i<speciesNames_.size();i++)
                 {
                     Xfluid_[i] = X_[i][cellI];
-		    if (Xfluid_[i] <= 0.0) Xfluid_[i] = 0.0;
+                    if (Xfluid_[i] <= 0.0) Xfluid_[i] = 0.0;
                 }
             }
             //fill arrays
@@ -291,7 +291,7 @@ void species::execute()
                 Info << "partTemp_[index][0] = " << partTemp_[index][0] << endl;
                 Info << "Tfluid = " << Tfluid << endl  ;
                 Info << "voidfraction =" << voidfraction << endl;
-                Info << "molarConc_" << molarConc_ << endl;
+                // Info << "molarConc_" << molarConc_ << endl;
             }
         }
     }
@@ -333,8 +333,8 @@ void species::execute()
             Info << "total conversion of species" << speciesNames_[i] << " = " << gSum(changeOfSpeciesMassFields_[i]*1.0*changeOfSpeciesMassFields_[i].mesh().V() * Nevery_ * timestep) << endl;
         }
         massSourceCurr_ = gSum(changeOfGasMassField_*1.0*changeOfGasMassField_.mesh().V() * Nevery_ * timestep);
-	massSourceTot_ += massSourceCurr_;
-	Info << "total conversion of mass:\n\tcurrent source = " << massSourceCurr_ << "\n\ttotal source = " << massSourceTot_ << "\n" << endl;
+        massSourceTot_ += massSourceCurr_;
+        Info << "total conversion of mass:\n\tcurrent source = " << massSourceCurr_ << "\n\ttotal source = " << massSourceTot_ << "\n" << endl;
         Info << "get data done" << endl;
 }
 
diff --git a/tutorials/cfdemSolverRhoPimpleChem/Shrinking_particle_model/DEM/in.liggghts_run b/tutorials/cfdemSolverRhoPimpleChem/Shrinking_particle_model/DEM/in.liggghts_run
index 90b0e00e..e40ad321 100644
--- a/tutorials/cfdemSolverRhoPimpleChem/Shrinking_particle_model/DEM/in.liggghts_run
+++ b/tutorials/cfdemSolverRhoPimpleChem/Shrinking_particle_model/DEM/in.liggghts_run
@@ -51,8 +51,7 @@ fix     cfd2 all couple/cfd/force
 fix     cfd3 all couple/cfd/chemistry n_species 3 species_names O2 CO2 N2
 
 # this should shrink the particle
-#fix    cfd4 all chem/shrink speciesA O2 molMassA 31.99 speciesC CO2 molMassC 44.01 molMassB 12.01 k 2.5e3 rdef 0.2
-fix     cfd4 all chem/shrink speciesA O2 molMassA 31.99 speciesC CO2 molMassC 44.01 molMassB 12.01 k 2.5e-1 rmin 0.005 nevery 1
+fix     cfd4 all chem/shrink speciesA O2 molMassA 0.03199 nuA 1 speciesC CO2 molMassC 0.04401 nuC 1 molMassB 0.01201 nuB 1 k 2.5e-1 rmin 0.005 nevery 1 screen no
 
 # apply nve integration to all particles that are inserted as single particles
 fix     integr all nve/sphere