anchor.diff:          !    ek: single particle eigen value.
           !    mu: chemical potential.
           !    beta: inverse temperature.
  +        ! Returns:
  +        !    fk: Fermi factor for single-particle state ek.
   
           real(p), intent(in) :: ek, mu, beta
  +
           real(p) :: f_k
   
           f_k = 1.0_p / (exp(beta*(ek-mu))+1)
   
       end function fermi_factor
   
  -    function particle_number(sys, beta, mu) result (nav)
  +    pure function particle_number(sys, beta, mu) result (nav)
   
           ! Evaluate the expected number of particles for given chemical potential
           ! mu at temperature T = 1 / beta.
  
positive.code: 
  contains
  
      function fermi_factor(ek, mu, beta) result(f_k)
  
          ! Calculate the Fermi factor, f_k = 1/(e^{beta(e_k-mu)}+1).
  
  
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.,
  