anchor.diff:      res = output_array(vector_length, 1)
       print *, 'dE/dx_i for i =', ideriv, 'is', res
   
  -    ! test against reference numbers
  +    ! compare with reference
       call assert((abs(res - res_reference(ideriv)) < 1.0d-6), &
                   "derivatives do not match reference numbers")
     end do
  
positive.code:   ierr = xc_set(id, 'pbe', 1.0d0)
    call assert(ierr == 0, "functional name not recognized")
  
    order = 0
    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))
    density = 0.0d0
    do ipoint = 1, num_points
      ! we use fantasy values here
  
negative.code:     density(1, 4, ipoint) = 4.0d0 !nabla_z n
    end do
  
    allocate(input_array(num_variables*vector_length, num_points))
    do ipoint = 1, num_points
      input_array(:, ipoint) = (/density(:, 1, ipoint), &
                                 density(:, 2, ipoint), &
                                 density(:, 3, ipoint), &
                                 density(:, 4, ipoint)/)
    end do
  
    allocate(output_array(vector_length, num_points))
    call xc_eval(id, num_points, input_array, output_array)
    res = output_array(vector_length, 1)
    print *, 'The XC energy density is', res
  
    ! compare with reference
    call assert((abs(-0.86494159400066051d0 - res) < 1.0d-6), &
                "derivatives do not match reference numbers")
  
    deallocate (density)
    deallocate (input_array)
    deallocate (output_array)
  
  
  
  
  
  
  
  
  
    ! now let's compute the first derivatives ('potential')
    ! and contract them with the first order densities
  