diff --git a/doc/manual/parameters.tex b/doc/manual/parameters.tex
index 79b3e3f9d..76c7acd8c 100644
--- a/doc/manual/parameters.tex
+++ b/doc/manual/parameters.tex
@@ -2291,7 +2291,7 @@ \subsection{Parameters in section \tt Post-processing Options}
 {\it Default:} false
 
 
-{\it Description:} [Standard] Write bands for every k-point to an outputfile called 'bands.out' in the units of Ha. This can be used after GS (Ground-state) or NSCF (Non-Self consistent field iteration) modes of solve. This option is by default on for NSCF mode of solve. Outputs a file name 'bands.out'. The first line has 2 entries with first one denoting the number of k-points and second entry denoting the number of eigenvalues(bands) for each k-point. Subsequent lines have 3 columns with first column indicating the k-point index, second column indicating band index and third column indicating corresponding eigenvalue.
+{\it Description:} [Standard] Write bands for every k-point to an outputfile called 'bands.out' in the units of Ha. This can be used after GS (Ground-state) or NSCF (Non-Self consistent field iteration) modes of solve. This option is by default on for NSCF mode of solve. Outputs a file name 'bands.out'. The first line has 3 entries with first one denoting the number of k-points and second entry denoting the number of eigenvalues(bands) for each k-point and third the fermi energy in Ha. Subsequent lines have 4 columns with first column indicating the k-point index, second column indicating band index, third column indicating corresponding eigenvalue and fourth column indicating the corresponding occupation number.
 
 
 {\it Possible values:} A boolean value (true or false)
diff --git a/src/dft/dft.cc b/src/dft/dft.cc
index 6b47c29d9..5e0100a28 100644
--- a/src/dft/dft.cc
+++ b/src/dft/dft.cc
@@ -4472,7 +4472,7 @@ namespace dftfe
       {
         FILE *pFile;
         pFile = fopen("bands.out", "w");
-        fprintf(pFile, "%d %d\n", totkPoints, numberEigenValues);
+        fprintf(pFile, "%d %d %.14g\n", totkPoints, numberEigenValues, FE);
         for (unsigned int kPoint = 0;
              kPoint < totkPoints / (1 + d_dftParamsPtr->spinPolarized);
              ++kPoint)
@@ -4481,16 +4481,33 @@ namespace dftfe
               {
                 if (d_dftParamsPtr->spinPolarized)
                   {
+                    double occupancyUp = dftUtils::getPartialOccupancy(
+                      eigenValuesFlattenedGlobal[2 * kPoint * d_numEigenValues +
+                                                 iWave],
+                      FE,
+                      C_kb,
+                      d_dftParamsPtr->TVal);
+
+                    double occupancyDown = dftUtils::getPartialOccupancy(
+                      eigenValuesFlattenedGlobal[(2 * kPoint + 1) *
+                                                   d_numEigenValues +
+                                                 iWave],
+                      FE,
+                      C_kb,
+                      d_dftParamsPtr->TVal);
+
                     fprintf(
                       pFile,
-                      "%d  %d   %.14g   %.14g\n",
+                      "%d  %d   %.14g   %.14g   %.14g   %.14g\n",
                       kPoint,
                       iWave,
                       eigenValuesFlattenedGlobal[2 * kPoint * d_numEigenValues +
                                                  iWave],
                       eigenValuesFlattenedGlobal[(2 * kPoint + 1) *
                                                    d_numEigenValues +
-                                                 iWave]);
+                                                 iWave],
+                      occupancyUp,
+                      occupancyDown);
                     if (d_dftParamsPtr->reproducible_output &&
                         d_dftParamsPtr->verbosity == 0)
                       {
@@ -4506,20 +4523,33 @@ namespace dftfe
                             (eigenValuesFlattenedGlobal
                                [(2 * kPoint + 1) * d_numEigenValues + iWave])) /
                           1000000000.0;
+                        double occupancyUpTrunc =
+                          std::floor(1000000000 * (occupancyUp)) / 1000000000.0;
+                        double occupancyDownTrunc =
+                          std::floor(1000000000 * (occupancyDown)) /
+                          1000000000.0;
                         pcout << kPoint << "  " << iWave << "  " << std::fixed
                               << std::setprecision(8) << eigenUpTrunc << "  "
-                              << eigenDownTrunc << std::endl;
+                              << eigenDownTrunc << "  " << occupancyUpTrunc
+                              << "  " << occupancyDownTrunc << std::endl;
                       }
                   }
                 else
                   {
+                    double occupancy = dftUtils::getPartialOccupancy(
+                      eigenValuesFlattenedGlobal[kPoint * d_numEigenValues +
+                                                 iWave],
+                      FE,
+                      C_kb,
+                      d_dftParamsPtr->TVal);
                     fprintf(
                       pFile,
-                      "%d  %d %.14g\n",
+                      "%d  %d %.14g %.14g\n",
                       kPoint,
                       iWave,
                       eigenValuesFlattenedGlobal[kPoint * d_numEigenValues +
-                                                 iWave]);
+                                                 iWave],
+                      occupancy);
                     if (d_dftParamsPtr->reproducible_output &&
                         d_dftParamsPtr->verbosity == 0)
                       {
@@ -4528,9 +4558,11 @@ namespace dftfe
                                      (eigenValuesFlattenedGlobal
                                         [kPoint * d_numEigenValues + iWave])) /
                           1000000000.0;
+                        double occupancyTrunc =
+                          std::floor(1000000000 * (occupancy)) / 1000000000.0;
                         pcout << kPoint << "  " << iWave << "  " << std::fixed
-                              << std::setprecision(8) << eigenTrunc
-                              << std::endl;
+                              << std::setprecision(8) << eigenTrunc << " "
+                              << occupancyTrunc << std::endl;
                       }
                   }
               }
diff --git a/utils/dftParameters.cc b/utils/dftParameters.cc
index 174bb530f..e4050f5fc 100644
--- a/utils/dftParameters.cc
+++ b/utils/dftParameters.cc
@@ -158,7 +158,7 @@ namespace dftfe
           "WRITE BANDS",
           "false",
           dealii::Patterns::Bool(),
-          "[Standard] Write bands for every k-point to an outputfile called 'bands.out' in the units of Ha. This can be used after GS (Ground-state) or NSCF (Non-Self consistent field iteration) modes of solve. This option is by default on for NSCF mode of solve. Outputs a file name 'bands.out'. The first line has 2 entries with first one denoting the number of k-points and second entry denoting the number of eigenvalues(bands) for each k-point. Subsequent lines have 3 columns with first column indicating the k-point index, second column indicating band index and third column indicating corresponding eigenvalue.");
+          "[Standard] Write bands for every k-point to an outputfile called 'bands.out' in the units of Ha. This can be used after GS (Ground-state) or NSCF (Non-Self consistent field iteration) modes of solve. This option is by default on for NSCF mode of solve. Outputs a file name 'bands.out'. The first line has 3 entries with first one denoting the number of k-points and second entry denoting the number of eigenvalues(bands) for each k-point and third the fermi energy in Ha. Subsequent lines have 4 columns with first column indicating the k-point index, second column indicating band index, third column indicating corresponding eigenvalue and fourth column indicating the corresponding occupation number.");
       }
       prm.leave_subsection();