anchor.diff:  !  density(1, :, :)  unperturbed values
   !  density(2, :, :)  perturbed values
   !  density(3, :, :)  ones and zeros
  -!  density(4, :, :)  not used
  -!  in general there will be 2^order elements, 1 ground state, 2^order-2 perturbed and one unused.
  +!  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
  
positive.code:                                  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:   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))
  