diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.C b/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.C
index ce5cf4cb..d9dc0374 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.C
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.C
@@ -67,28 +67,26 @@ diffusionCoefficient::diffusionCoefficient
     ),
     // create a list from the Species table in the specified diffusionCoefficient dictionary
     speciesNames_(specDict_.lookup("species")),
-    tempFieldName_(propsDict_.lookup("tempFieldName")),
+    tempFieldName_(propsDict_.lookupOrDefault<word>("tempFieldName","T")),
     tempField_(sm.mesh().lookupObject<volScalarField> (tempFieldName_)),
-    pressureFieldName_(propsDict_.lookup("pressureFieldName")),
+    pressureFieldName_(propsDict_.lookupOrDefault<word>("pressureFieldName","p")),
     P_(sm.mesh().lookupObject<volScalarField>(pressureFieldName_)),
-    densityFieldName_(propsDict_.lookup("densityFieldName")),
+    densityFieldName_(propsDict_.lookupOrDefault<word>("densityFieldName","rho")),
     rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
-    totalMoleFieldName_(propsDict_.lookup("totalMoleFieldName")),
+    totalMoleFieldName_(propsDict_.lookupOrDefault<word>("totalMoleFieldName","N")),
     // needed to calculate the mixture diffusion coefficient
     // dcoeff is dependent on molar fraction not mass fraction
     N_(sm.mesh().lookupObject<volScalarField>(totalMoleFieldName_)),
     Y_(speciesNames_.size()),
-    // X_(speciesNames_.size()),
     diffusantGasNames_(propsDict_.lookup("diffusantGasNames")),
-    diffusionCoefficients_(diffusantGasNames_.size(),NULL)
- /*   diffusionCoefficientNames_(propsDict_.lookup("diffusionCoefficientNames")),
-                               //volumeScalarFields created in the ts folders
-    diffusionCoefficients_(diffusionCoefficientNames_.size(),NULL),         //the value of diffusionCoefficient concentration for every diffusionCoefficient 
-    */
+    diffusionCoefficients_(diffusantGasNames_.size(),NULL),
+    X_(diffusantGasNames_.size(), NULL)
 {
     Info << " Reading diffusionCoefficient list: " << diffusantGasNames_ << endl;
     for (int i = 0; i < diffusantGasNames_.size(); i++)
-        Info << " Diffusant names: " << diffusantGasNames_[i] << endl;
+    {
+        Info << "Diffusant names: " << diffusantGasNames_[i] << endl;
+    }
 
     for (int i=0; i<speciesNames_.size(); i++)
     {
@@ -98,28 +96,6 @@ diffusionCoefficient::diffusionCoefficient
         particleCloud_.checkCG(false);
     }
 
-    // if mole fractions field should be generated ??
-    /*forAll(Y, i)
-    {
-        X_.set
-        (
-            i,
-            new volScalarField
-            (
-                IOobject
-                (
-                    "X_" + Y_[i].name(),
-                    mesh_.time().timeName(),
-                    mesh_,
-                    IOobject::NO_READ,
-                    IOobject::AUTO_WRITE
-                ),
-                mesh_,
-                dimensionedScalar("X", dimless, 0)
-            )
-        );
-    } */
-
     allocateMyArrays();
     createCoeffs();
     molWeightTable();
@@ -136,6 +112,7 @@ diffusionCoefficient::~diffusionCoefficient()
     for (int i=0; i<diffusantGasNames_.size(); i++)
     {
         particleCloud_.dataExchangeM().destroy(diffusionCoefficients_[i],nP_);
+        particleCloud_.dataExchangeM().destroy(X_[i], nP_);
     }
 }
 
@@ -148,7 +125,7 @@ diffusionCoefficient::~diffusionCoefficient()
         for (int i=0; i<diffusantGasNames_.size(); i++)
         {
             particleCloud_.dataExchangeM().allocateArray(diffusionCoefficients_[i],initVal,1,"nparticles");
-
+            particleCloud_.dataExchangeM().allocateArray(X_[i],initVal,1,"nparticles");
         }
     }
 }
@@ -162,6 +139,7 @@ void diffusionCoefficient::reAllocMyArrays() const
         for (int i=0; i<diffusantGasNames_.size(); i++)
         {
             particleCloud_.dataExchangeM().allocateArray(diffusionCoefficients_[i],initVal,1);
+            particleCloud_.dataExchangeM().allocateArray(X_[i],initVal,1);
         }
     }
 }
@@ -241,7 +219,7 @@ void diffusionCoefficient::execute()
                 // and for every reactant gas
                 for (int j=0; j < speciesNames_.size();j++)
                 {
-                    // if diffusant gas ans reactant gas are not equal get diffusion coefficients
+                    // if diffusant gas and reactant gas are not equal get diffusion coefficients
                     if (diffusantGasNames_[i] != speciesNames_[j])
                     {
                         Info << "molar weights diffuser gases: " << molWeight(speciesNames_[j]) << nl << endl;
@@ -277,6 +255,9 @@ void diffusionCoefficient::execute()
                             Info << "molar fraction diffusing gases - 2:" << Xfluid_[i] << nl << endl;
                             Info << "molar fraction stagnant gases - 2:" << Xfluid_[j] << nl << endl;
 
+                            // fill X array with diffusant gas molar fraction to be transferred to DEM
+                            X_[i][index][0] = Xfluid_[i];
+
                             // for now considering only one species is diffusing, the rest of the molar fractions are from
                             // stagnangt gases we can get the total stagnant gas molar fractions by
                             Xstag_tot[i]   =   1-Xfluid_[i];
@@ -346,6 +327,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 << "molar fraction of species" << diffusantGasNames_[i] << " = " << X_[i][index][0] << endl;
             }
         }
     }
@@ -354,6 +336,8 @@ void diffusionCoefficient::execute()
     {
         word pushName = diffusantGasNames_[i] + "_diffCoeff";
         particleCloud_.dataExchangeM().giveData(pushName,"scalar-atom",diffusionCoefficients_[i]);
+        word moleFractionName = "X_" + diffusantGasNames_[i];
+        particleCloud_.dataExchangeM().giveData(moleFractionName, "scalar-atom", X_[i]);
     };
 
     Info << "give data done" << endl;
diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.H b/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.H
index 588dbbae..33a30e0f 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.H
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/diffusionCoefficients/diffusionCoefficients.H
@@ -83,12 +83,12 @@ private:
 
     UPtrList<volScalarField> Y_;
 
-    // PtrList<volScalarField> X_;
-
     wordList diffusantGasNames_;
     
     mutable List<double**> diffusionCoefficients_;
 
+    mutable List<double**> X_;
+
     HashTable<scalar, word> coeffs;
 
     HashTable<scalar, word> molWeight;
diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.C b/src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.C
index 1341cf33..f1623605 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.C
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/massTransferCoeff/massTransferCoeff.C
@@ -56,15 +56,15 @@ massTransferCoeff::massTransferCoeff
     verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
     interpolation_(propsDict_.lookupOrDefault<bool>("interpolation",false)),
     mesh_(sm.mesh()),
-    velFieldName_(propsDict_.lookup("velFieldName")),
+    velFieldName_(propsDict_.lookupOrDefault<word>("velFieldName","U")),
     U_(mesh_.lookupObject<volVectorField>(velFieldName_)),
-    voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")),
+    voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
     voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
     densityFieldName_(propsDict_.lookupOrDefault<word>("densityFieldName","rho")),
     rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
-    partNuName_(propsDict_.lookup("partNu")),
+    partNuName_(propsDict_.lookupOrDefault<word>("partViscos","partNu")),
     partNu_(NULL),
-    partReynolds_(propsDict_.lookup("partReynolds")),
+    partReynolds_(propsDict_.lookupOrDefault<word>("partReynolds","partRe")),
     partRe_(NULL)
 {}
 
diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C b/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C
index e1cea815..551b5994 100644
--- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C
+++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C
@@ -86,20 +86,20 @@ species::species
         mesh_,
         dimensionedScalar("zero",dimMass/(dimVol*dimTime),0.0)
     ),
-    tempFieldName_(propsDict_.lookup("tempFieldName")),
+    tempFieldName_(propsDict_.lookupOrDefault<word>("tempFieldName","T")),
     tempField_(sm.mesh().lookupObject<volScalarField> (tempFieldName_)),
-    partTempName_(propsDict_.lookup("partTempName")),
+    partTempName_(propsDict_.lookupOrDefault<word>("partTempName","partTemp")),
     partTemp_(NULL),
-    densityFieldName_(propsDict_.lookup("densityFieldName")),
+    densityFieldName_(propsDict_.lookupOrDefault<word>("densityFieldName","rho")),
     rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
-    partRhoName_(propsDict_.lookup("partRhoName")),
+    partRhoName_(propsDict_.lookupOrDefault<word>("partRhoName","partRho")),
     partRho_(NULL),
-    voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")),
+    voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
     voidfraction_(sm.mesh().lookupObject<volScalarField>(voidfractionFieldName_)),
     // total mole field
-    totalMoleFieldName_(propsDict_.lookup("totalMoleFieldName")),
+    totalMoleFieldName_(propsDict_.lookupOrDefault<word>("totalMoleFieldName","N")),
     N_(sm.mesh().lookupObject<volScalarField>(totalMoleFieldName_)),
-    partMoleName_(propsDict_.lookup("partMoleName")),
+    partMoleName_(propsDict_.lookupOrDefault<word>("partMoleName","partN")),
     partN_(NULL),
     loopCounter_(-1),
     Nevery_(propsDict_.lookupOrDefault<label>("Nevery",1)),
@@ -216,7 +216,7 @@ void species::execute()
     // realloc the arrays
     reAllocMyArrays();
 
-  // get Y_i, T, rho at particle positions, fill arrays with them and push to LIGGGHTS
+    // get Y_i, T, rho at particle positions, fill arrays with them and push to LIGGGHTS
 
     label  cellI=0;
     scalar Tfluid(0);
@@ -226,7 +226,6 @@ void species::execute()
     Yfluid_.setSize(speciesNames_.size());
     scalar Nfluid(0);
 
-
     // defining interpolators for T, rho, voidfraction, N
     interpolationCellPoint <scalar> TInterpolator_(tempField_);
     interpolationCellPoint <scalar> rhoInterpolator_(rho_);
@@ -260,8 +259,8 @@ void species::execute()
             }
             //fill arrays
             partTemp_[index][0] =   Tfluid;
-	    // partRho was filled with rhofluid*voidfraction before
-	    // probably wrong: need actual gas density, not averaged one
+            // partRho was filled with rhofluid*voidfraction before
+            // probably wrong: need actual gas density, not averaged one
             partRho_[index][0]  =   rhofluid;
             partN_[index][0]    =   Nfluid;
 
@@ -294,7 +293,7 @@ void species::execute()
 
         // give DEM data
         particleCloud_.dataExchangeM().giveData(partTempName_, "scalar-atom", partTemp_);
-        particleCloud_.dataExchangeM().giveData(partRhoName_, "scalar-atom", partRho_);
+        particleCloud_.dataExchangeM().giveData(partRhoName_,  "scalar-atom", partRho_);
         particleCloud_.dataExchangeM().giveData(partMoleName_, "scalar-atom", partN_);
 
         for (int i=0; i<speciesNames_.size();i++)
@@ -305,7 +304,7 @@ void species::execute()
         Info << "give data done" << endl;
 
         // pull changeOfSpeciesMass_, transform onto fields changeOfSpeciesMassFields_, add them up on changeOfGasMassField_
-	scalar timestep = mesh_.time().deltaTValue();
+        scalar timestep = mesh_.time().deltaTValue();
         changeOfGasMassField_.primitiveFieldRef() = 0.0;
         changeOfGasMassField_.boundaryFieldRef() = 0.0;
         for (int i=0; i<speciesNames_.size();i++)
@@ -322,7 +321,7 @@ void species::execute()
             );
 
             // take care for implementation in LIGGGHTS: species produced from particles defined positive
-	    // changeOf...Fields need to be mass per volume per timestep
+            // changeOf...Fields need to be mass per volume per timestep
             changeOfSpeciesMassFields_[i].primitiveFieldRef() /= (changeOfSpeciesMassFields_[i].mesh().V() * Nevery_ * timestep);
             changeOfSpeciesMassFields_[i].correctBoundaryConditions();
             changeOfGasMassField_ += changeOfSpeciesMassFields_[i];
diff --git a/tutorials/cfdemSolverRhoPimpleChem/test_case/CFD/constant/couplingProperties b/tutorials/cfdemSolverRhoPimpleChem/test_case/CFD/constant/couplingProperties
index 5f7a73a0..e0cb2c69 100644
--- a/tutorials/cfdemSolverRhoPimpleChem/test_case/CFD/constant/couplingProperties
+++ b/tutorials/cfdemSolverRhoPimpleChem/test_case/CFD/constant/couplingProperties
@@ -105,26 +105,14 @@ reactionHeatProps
 speciesProps
 {
 	ChemistryFile 		"$casePath/CFD/constant/foam.inp";
-	tempFieldName 		"T";
-	partTempName  		"partTemp";
-	densityFieldName 	"rho";
-	partRhoName		"partRho";
-        voidfractionFieldName   "voidfraction";
-        totalMoleFieldName      "N";
-        partMoleName            "partN";
-        pressureFieldName       "p";
-        partPName               "partP";
-        // verbose true;
+        // verbose	true;
+	// Nevery	1;
 }
 
 diffusionCoefficientsProps
 {
 	verbose			true;        
 	ChemistryFile 		"$casePath/CFD/constant/foam.inp";
-	tempFieldName 		"T";
-	totalMoleFieldName      "N";
-        pressureFieldName       "p";
-        densityFieldName 	"rho";
         diffusantGasNames       (    CO
                                      // H2
                                 ); 
@@ -132,10 +120,6 @@ diffusionCoefficientsProps
 
 massTransferCoeffProps
 {
-    velFieldName    		"U";
-    voidfractionFieldName   	"voidfraction";
-    partNu          		"partNu";
-    partReynolds   		"partRe";
     verbose 			true;
 }
 
diff --git a/tutorials/cfdemSolverRhoPimpleChem/test_case/CFD/system/controlDict b/tutorials/cfdemSolverRhoPimpleChem/test_case/CFD/system/controlDict
index 82d6e262..946debef 100644
--- a/tutorials/cfdemSolverRhoPimpleChem/test_case/CFD/system/controlDict
+++ b/tutorials/cfdemSolverRhoPimpleChem/test_case/CFD/system/controlDict
@@ -23,13 +23,13 @@ startTime       0;
 
 stopAt          endTime;
 
-endTime         5;
+endTime         4000;
 
 deltaT          0.1;
 
 writeControl    adjustableRunTime;  //runTime;//// timeStep; // 
 
-writeInterval   0.5;
+writeInterval   100;
 
 purgeWrite      0;
 
diff --git a/tutorials/cfdemSolverRhoPimpleChem/test_case/DEM/in.liggghts_run b/tutorials/cfdemSolverRhoPimpleChem/test_case/DEM/in.liggghts_run
index 2ab52e36..e1af20fa 100644
--- a/tutorials/cfdemSolverRhoPimpleChem/test_case/DEM/in.liggghts_run
+++ b/tutorials/cfdemSolverRhoPimpleChem/test_case/DEM/in.liggghts_run
@@ -64,7 +64,7 @@ fix     Ea all property/global Ea_cfd5 vector 69488 73674 113859
 # particle porosities adn tortuosity
 fix     porosity        all property/global porosity_       scalar 0.15
 fix     tortuosity      all property/global tortuosity_     scalar 3
-fix     pore_diameter   all property/global pore_diameter_  scalar 0.00001 #3.91e-9
+fix     pore_diameter   all property/global pore_diameter_  scalar 3.91e-5 #0.00001 #3.91e-9
 
 # Material properties for unreacted chemical shrink core (activate only when chem/shrink/core is active)
 fix     density all property/global density_all vector 7870 5740 5170 5240 
@@ -78,38 +78,38 @@ fix     integr all nve/sphere
 
 ###############################################
 ## DEBUG ##
-#compute masschange all property/atom mass
-#compute massreduce all reduce sum c_masschange
-#fix     rmass all ave/time 100 1 100 c_massreduce file rmass.dat 
+compute masschange all property/atom mass
+compute massreduce all reduce sum c_masschange
+fix     rmass all ave/time 100 1 100 c_massreduce file rmass.dat 
 
-#compute layerRad1 all reduce sum f_layerRelRad[1]
-#fix     redRad1 all ave/time 100 1 100 c_layerRad1 file relRad1.dat
+compute layerRad1 all reduce sum f_layerRelRad[1]
+fix     redRad1 all ave/time 100 1 100 c_layerRad1 file relRad1.dat
 
-#compute layerRad2 all reduce sum f_layerRelRad[2]
-#fix     redRad2 all ave/time 100 1 100 c_layerRad2 file relRad2.dat 
+compute layerRad2 all reduce sum f_layerRelRad[2]
+fix     redRad2 all ave/time 100 1 100 c_layerRad2 file relRad2.dat 
 
-#compute layerRad3 all reduce sum f_layerRelRad[3]
-#fix     redRad3 all ave/time 100 1 100 c_layerRad3 file relRad3.dat 
+compute layerRad3 all reduce sum f_layerRelRad[3]
+fix     redRad3 all ave/time 100 1 100 c_layerRad3 file relRad3.dat 
 
-#compute layerRad4 all reduce sum f_layerRelRad[4]
-#fix     redRad4 all ave/time 100 1 100 c_layerRad4 file relRad4.dat 
+compute layerRad4 all reduce sum f_layerRelRad[4]
+fix     redRad4 all ave/time 100 1 100 c_layerRad4 file relRad4.dat 
 
 
-#variable    rad_fe  atom   "f_layerRelRad[1]*radius"
-#compute r_fe    all reduce  sum v_rad_fe
-#fix         r_fe1   all     ave/time    100  1   100  c_r_fe  file    rfe.dat 
+variable    rad_fe  atom   "f_layerRelRad[1]*radius"
+compute r_fe    all reduce  sum v_rad_fe
+fix         r_fe1   all     ave/time    100  1   100  c_r_fe  file    rfe.dat 
 
-#variable    rad_w   atom   "f_layerRelRad[2]*radius"
-#compute     r_w    all     reduce      sum v_rad_w
-#fix         r_w1   all     ave/time    100  1   100  c_r_w  file    rw.dat 
+variable    rad_w   atom   "f_layerRelRad[2]*radius"
+compute     r_w    all     reduce      sum v_rad_w
+fix         r_w1   all     ave/time    100  1   100  c_r_w  file    rw.dat 
 
-#variable    rad_m  atom   "f_layerRelRad[3]*radius"
-#compute     r_m    all     reduce      sum v_rad_m
-#fix         r_m1   all     ave/time    100  1   100  c_r_m  file    rm.dat 
+variable    rad_m  atom   "f_layerRelRad[3]*radius"
+compute     r_m    all     reduce      sum v_rad_m
+fix         r_m1   all     ave/time    100  1   100  c_r_m  file    rm.dat 
 
-#variable    rad_h  atom   "f_layerRelRad[4]*radius"
-#compute     r_h    all     reduce      sum v_rad_h
-#fix         r_h1   all     ave/time    100  1   100  c_r_h  file    rh.dat 
+variable    rad_h  atom   "f_layerRelRad[4]*radius"
+compute     r_h    all     reduce      sum v_rad_h
+fix         r_h1   all     ave/time    100  1   100  c_r_h  file    rh.dat 
 
 ###############################################
 # screen output