anchor.diff:    implicit none
   
     integer, parameter :: num_points = 1
  +  integer, parameter :: num_variables = 4  ! we will use XC_N_NX_NY_NZ
   
     character(1000) :: text
     integer :: id, order, ierr, ideriv, ipoint
     integer :: vector_length
  -  integer :: num_variables
     real(8) :: res
   
     real(8), allocatable :: density(:, :, :)
  
positive.code:   ! now second derivative of the "potential", contracted with ONE perturbed density.
    ! hopefully the strange input/output_array format will start to make sense
    order = 2
    num_variables = 4
    ierr = xc_eval_setup(id, XC_N_NX_NY_NZ, XC_CONTRACTED, order)
    call assert(ierr == 0, "xc_eval_setup failed")
  
  
negative.code:   ! and contract them with the first order densities
    ! first set up xcfun for first derivatives
    order = 1
    ierr = xc_eval_setup(id, XC_N_NX_NY_NZ, XC_CONTRACTED, order)
    call assert(ierr == 0, "xc_eval_setup failed")
  
  