anchor.diff:    deallocate (input_array)
     deallocate (output_array)
   
  -
  -!  now third derivative of the "potential", contracted with perturbed densities.
  +  !  now third derivative of the "potential", contracted with perturbed densities.
     order = 3
     num_variables = 4
     ierr = xc_eval_setup(id, XC_N_NX_NY_NZ, XC_CONTRACTED, order)
  
positive.code:   endif
  
    vector_length = 2**order
  !  allocate density and fill in some fantasy values
    allocate (density(vector_length, num_variables, num_points))
    density = 0.0d0
    do ipoint = 1, num_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
  
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)
  