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:     density(1, 4, ipoint) = 4.0d0 !nabla_z n
    end do
  
    allocate (input_array(num_variables*vector_length, num_points))
    allocate (output_array(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
    call xc_eval(id, num_points, input_array, output_array)
    res = output_array(vector_length, 1)
    print *, 'The XC energy density is', res
  
    deallocate (density)
    deallocate (input_array)
    deallocate (output_array)
  
    ! now let's compute the first derivatives ('potential')
    ! and contract them with the first order densities
  
negative.code:     res = output_array(vector_length, 1)
      print *, 'd^2 E/dx_i dD_pert for i =', ideriv, 'is', res
  
      ! compare with reference
      call assert((abs(res - res_reference(ideriv)) < 1.0d-6), &
                  "derivatives do not match reference numbers")
    end do
  