anchor.diff:                                   density(:, 4, ipoint)/)
       end do
       call xc_eval(id, NR_POINTS, input_array, output_array)
  +
       res = output_array(vector_length, 1)
       print *, 'd^2 E/dx_i dD_pert for i =', ideriv, 'is', res
   
  
positive.code:   res_reference(2) = 2.1489169824673575d-2
    res_reference(3) = 4.1214059228412030d-2
    res_reference(4) = 6.0938948632150471d-2
  
  !  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
  
negative.code:   deallocate (input_array)
    deallocate (output_array)
  
    !  now let's compute the first derivatives ('potential')
    !  and contract them with the first order densities
  !  first set up xcfun for first derivatives
    order = 1
    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
      stop 1
    endif
  
    vector_length = 2**order
  !  allocate density and fill in some phantasy values
    allocate (density(vector_length, nr_variables, NR_POINTS))
    density = 0.0d0
    do ipoint = 1, NR_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
      density(1, 4, ipoint) = 4.0d0 !nabla_z n zeroth order
      density(2, 1, ipoint) = 5.0d0 !        n first order
      density(2, 2, ipoint) = 6.0d0 !nabla_x n first order
      density(2, 3, ipoint) = 7.0d0 !nabla_y n first order
      density(2, 4, ipoint) = 8.0d0 !nabla_z n first order
    end do
  
    allocate (input_array(nr_variables*vector_length, NR_POINTS))
    allocate (output_array(vector_length, NR_POINTS))
  
  !     It is possible to put in other numbers than 1 and 0
  !     this is where perturbed densities go for automatic contraction
  
      do ipoint = 1, NR_POINTS
        input_array(:, ipoint) = (/density(:, 1, ipoint), &
                                   density(:, 2, ipoint), &
                                   density(:, 3, ipoint), &
                                   density(:, 4, ipoint)/)
      end do
      call xc_eval(id, NR_POINTS, input_array, output_array)
  
    deallocate (density)
    deallocate (input_array)
    deallocate (output_array)
    
    
  !  reference values for self test
    allocate (res_reference(nr_variables))
    res_reference(1) = -0.97182634532897016d0
    res_reference(2) = -8.9803044914016916d-3
    res_reference(3) = -1.3470456737102541d-2
    res_reference(4) = -1.7960608982803383d-2  
    
    
  !  now let's compute the first derivatives ('potential')
  !  and do NOT contract them!
  !  first set up xcfun for first derivatives
    order = 1
    nr_variables = 4
  