anchor.diff:  
     vector_length = 2**order
   !  allocate density and fill in some fantasy values
  -  allocate (density(vector_length, nr_variables, NR_POINTS))
  -  density = 0.0d0
  -  do ipoint = 1, NR_POINTS
  +  allocate (density(vector_length, num_variables, num_points))
  +  density = 0.0d0
  +  do ipoint = 1, num_points
       density(1, 1, ipoint) = 1.0d0 !        n zeroth order
       density(1, 2, ipoint) = 2.0d0 !nabla_x n zeroth order
       density(1, 3, ipoint) = 3.0d0 !nabla_y n zeroth order
  
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
    nr_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
  
negative.code:     density(4, 4, ipoint) = 0.0d0    ! first order
    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)
  
    !  The output array will then contain:
    ! 1.) Product of functional with ground state density
  