anchor.diff:  !  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
  +  do ideriv = 1, num_variables
  +    do ipoint = 1, num_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
  +    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, NR_POINTS, input_array, output_array)
  +    call xc_eval(id, num_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: 
    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
      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
    end do
  
    allocate (input_array(nr_variables*vector_length, NR_POINTS))
    allocate (output_array(vector_length, NR_POINTS))
  
  !  here is a tricky part, we need a loop and we need to put
  !  in a '1' where we want the derivative
  !  we loop over all variables we want a derivative with respect to
    do ideriv = 1, nr_variables
      do ipoint = 1, NR_POINTS
        density(2, :, ipoint) = 0.0d0
        density(2, ideriv, ipoint) = 1.0d0
      end do
  !     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)
  
      res = output_array(vector_length, 1)
      print *, 'dE/dx_i for i =', ideriv, 'is', res
  
negative.code:   !  2.) perturbed density a) (first order)
    !  3.) perturbed density b) (first order)
    !  4.) perturbed density ab) (second order)
    allocate (density(vector_length, num_variables, num_points))
    density = 0.0d0
    do ipoint = 1, num_points
      density(1, 1, ipoint) = 1.0d0    ! zeroth order
      density(1, 2, ipoint) = 2.0d0    ! zeroth order
      density(1, 3, ipoint) = 3.0d0    ! zeroth order
  