anchor.diff:  
   !  now third derivative of the "potential", contracted with perturbed densities.
     order = 3
  -  nr_variables = 4
  +  num_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
  
positive.code: !  density(3, :, :)  ones and zeros
  !  density(4, :, :)  not used in this case
  !  For the second order case, there will be 2^order elements, 1 ground state, 2^order-2 perturbed and one unused.
    do ideriv = 1, nr_variables
      do ipoint = 1, NR_POINTS
        density(3, :, ipoint) = 0.0d0
        density(3, ideriv, ipoint) = 1.0d0
      end do
  !     It is possible to put in other numbers than 1 and 0 in the density2 array
  !     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)
  
      res = output_array(vector_length, 1)
      print *, 'd^2 E/dx_i dD_pert for i =', ideriv, 'is', res
  
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
  