anchor.diff:    deallocate (output_array)
     deallocate (res_reference)
   
  -!  reference values for self test
  +  ! reference values for self test
     allocate (res_reference(num_variables))
     res_reference(1) = -2.0795476275887927d0
     res_reference(2) = 2.1489169824673575d-2
  
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:   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
    num_variables = 4
    ierr = xc_eval_setup(id, XC_N_NX_NY_NZ, XC_CONTRACTED, order)
  