anchor.diff:    order = 2
     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
  +  call assert(ierr == 0, "xc_eval_setup failed")
   
     vector_length = 2**order
     allocate (density(vector_length, num_variables, num_points))
  
positive.code:   order = 0
    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:     print *, 'd^2 E/dx_i dD_pert for i =', ideriv, 'is', res
  
      ! test against reference numbers
      call assert((abs(res - res_reference(ideriv)) < 1.0d-6), &
                  "derivatives do not match reference numbers")
    end do
  
    deallocate (density)
  