anchor.diff:          it = 0
   
           ! At zero temperature the chemical potential should be close to the
  -        ! energy of the highest occupied orbital (Fermi Energy). We assume
  -        ! we only deal with completely polarised or unpolarised systems
  +        ! energy of the highest occupied orbital (Fermi Energy).
           ! This is typically a good guess for root finding as the chemical
           ! potential is a monotonically decreasing function of temperature, and
           ! its value at T=0 is the Fermi energy (in the TDL).
           ! The following takes care of the spin polarised and unpolarised case given
  -        ! how we enumerate basis functions.
  +        ! how we enumerate basis functions. We assume we only deal with completely
  +        ! polarised or unpolarised systems
           mu_old = sys%basis%basis_fns(2*sys%nalpha)%sp_eigv
   
           do while (it < max_it)
  
positive.code: module chem_pot
  
  use const, only: p, depsilon
  
  
negative.code: 
      end function particle_number
  
      pure function deriv_particle_number(sys, beta, mu) result(nav_deriv)
  
          ! Derivative of the expected number of particles with respect to the
          ! chemical potential, i.e.,
  