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