diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90
index 64dc908..6b64c1b 100644
--- a/src/HF/HFB.f90
+++ b/src/HF/HFB.f90
@@ -1,6 +1,6 @@
 subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc,     & 
                nBas,nOrb,nO,S,T,V,Hc,ERI,dipole_int,X,EHFB,eHF,c,P,Panom,F,Delta, &
-               temperature,sigma)
+               temperature,sigma,chem_pot_hf)
 
 ! Perform Hartree-Fock Bogoliubov calculation
 
@@ -34,6 +34,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc,
 
 ! Local variables
 
+  logical                       :: chem_pot_hf
   integer                       :: nBas2
   integer                       :: iorb
   integer                       :: nSCF
@@ -135,6 +136,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc,
    Occ(:)     = 0d0
    Occ(1:nO)  = 1d0
    call fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF)
+   if(chem_pot_hf) chem_pot = 0.5d0*(eHF(nO)+eHF(nO+1))
    P(:,:)      = 0d0
    Panom(:,:)  = 0d0
    do iorb=1,nOrb
diff --git a/src/QuAcK/BQuAcK.f90 b/src/QuAcK/BQuAcK.f90
index 572118b..5bf5637 100644
--- a/src/QuAcK/BQuAcK.f90
+++ b/src/QuAcK/BQuAcK.f90
@@ -1,6 +1,6 @@
 subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc,                                    &
                   S,T,V,Hc,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift,                                  &
-                  guess_type,mix,temperature,sigma)
+                  guess_type,mix,temperature,sigma,chem_pot_hf)
 
 ! Restricted branch of QuAcK
 
@@ -13,6 +13,7 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,
 
   logical,intent(in)            :: doHFB
 
+  logical,intent(in)            :: chem_pot_hf
   integer,intent(in)            :: nNuc,nBas,nOrb
   integer,intent(in)            :: nC
   integer,intent(in)            :: nO
@@ -94,7 +95,7 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,
     call wall_time(start_HF)
     call HFB(dotest,maxSCF_HF,thresh_HF,max_diis_HF,level_shift,nNuc,ZNuc,rNuc,ENuc, &
              nBas,nOrb,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHFB,eHF,cHF,PHF,PanomHF,  &
-             FHF,Delta,temperature,sigma)
+             FHF,Delta,temperature,sigma,chem_pot_hf)
     call wall_time(end_HF)
 
     t_HF = end_HF - start_HF
diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90
index af50d21..5081fe3 100644
--- a/src/QuAcK/QuAcK.f90
+++ b/src/QuAcK/QuAcK.f90
@@ -73,6 +73,7 @@ program QuAcK
 
   logical                       :: dotest,doRtest,doUtest,doGtest
 
+  logical                       :: chem_pot_hf
   double precision              :: temperature,sigma
 
   character(len=256)            :: working_dir
@@ -137,7 +138,7 @@ program QuAcK
                     maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T,                 & 
                     doACFDT,exchange_kernel,doXBS,                                              &
                     dophBSE,dophBSE2,doppBSE,dBSE,dTDA,                                         &
-                    temperature,sigma)
+                    temperature,sigma,chem_pot_hf)
 
 !------------------!
 ! Hardware         !
@@ -294,7 +295,7 @@ program QuAcK
   if(doBQuAcK) & 
     call BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc,                           &
                 S,T,V,Hc,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,          &
-                temperature,sigma)
+                temperature,sigma,chem_pot_hf)
 
 !-----------!
 ! Stop Test !
diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90
index 4fb3e4d..778e305 100644
--- a/src/QuAcK/read_options.f90
+++ b/src/QuAcK/read_options.f90
@@ -8,7 +8,7 @@ subroutine read_options(working_dir,
                         maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T,                 &
                         doACFDT,exchange_kernel,doXBS,                                              &
                         dophBSE,dophBSE2,doppBSE,dBSE,dTDA,                                         &
-                        temperature,sigma)
+                        temperature,sigma,chem_pot_hf)
 
 ! Read desired methods 
 
@@ -73,6 +73,7 @@ subroutine read_options(working_dir,
   logical,intent(out)           :: dBSE
   logical,intent(out)           :: dTDA
 
+  logical,intent(out)           :: chem_pot_hf
   double precision,intent(out)  :: temperature
   double precision,intent(out)  :: sigma
 
@@ -224,9 +225,13 @@ subroutine read_options(working_dir,
       ! Options for Hartree-Fock Bogoliubov
     
       temperature  = 0d0
+      sigma        = 1d0
+      chem_pot_hf  = .false.
     
       read(1,*) 
-      read(1,*) temperature,sigma
+      read(1,*) temperature,sigma,ans1
+
+      if(ans1 == 'T') chem_pot_hf  = .true.
     
     endif