anchor.diff:    allocate (input_array(num_variables*vector_length, num_points))
     allocate (output_array(vector_length, num_points))
   
  -!  now we have to put in ones and zeroes again to generate derivatives (which then generate matrix elements in ADF)
  -!  density(1, :, :)  unperturbed values
  -!  density(2, :, :)  perturbed values
  -!  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.
  +  ! now we have to put in ones and zeroes again to generate derivatives (which then generate matrix elements in ADF)
  +  ! density(1, :, :)  unperturbed values
  +  ! density(2, :, :)  perturbed values
  +  ! 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, 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
  +    ! 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, num_points
         input_array(:, ipoint) = (/density(:, 1, ipoint), &
  
positive.code:   allocate (input_array(num_variables*vector_length, num_points))
    allocate (output_array(vector_length, num_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, num_variables
      do ipoint = 1, num_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, num_points
        input_array(:, ipoint) = (/density(:, 1, ipoint), &
  
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
    num_variables = 4
    ierr = xc_eval_setup(id, XC_N_NX_NY_NZ, XC_CONTRACTED, order)
  