anchor.diff:    ! and do NOT contract them!
     ! first set up xcfun for first derivatives
     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")
   
  
positive.code:   implicit none
  
    integer, parameter :: num_points = 1
  
    character(1000) :: text
    integer :: id, order, ierr, ideriv, ipoint
    integer :: vector_length
    integer :: num_variables
    real(8) :: res
  
    real(8), allocatable :: density(:, :, :)
  
negative.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
    ierr = xc_eval_setup(id, XC_N_NX_NY_NZ, XC_CONTRACTED, order)
    call assert(ierr == 0, "xc_eval_setup failed")
  
  