anchor.diff:    ! in this example we use PBE
     print *, 'Setting up PBE'
     ierr = xc_set(id, 'pbe', 1.0d0)
  -  if (ierr /= 0) then
  -    print *, 'Functional name not recognized, quitting.'
  -    stop 1
  -  endif
  +  call assert(ierr == 0, "functional name not recognized")
   
     ! First we just compute the energy, i.e. the 0-th order integrand.
     ! We have one gridpoint, and four variables, density N and gradient
  
positive.code:   order = 1
    num_variables = 4
    ierr = xc_eval_setup(id, XC_N_NX_NY_NZ, XC_CONTRACTED, order)
    if (ierr /= 0) then
      print *, 'xc_eval_setup failed with error ', ierr
      stop 1
    endif
  
    vector_length = 2**order
    allocate (density(vector_length, num_variables, num_points))
  
negative.code:   order = 1
    num_variables = 4
    ierr = xc_eval_setup(id, XC_N_NX_NY_NZ, XC_CONTRACTED, order)
    call assert(ierr == 0, "xc_eval_setup failed")
  
    vector_length = 2**order
    allocate (density(vector_length, num_variables, num_points))
  