anchor.diff:  !
   ! For information on the complete list of contributors to the
   ! XCFun library, see: <https://xcfun.readthedocs.io/>
  -
  -
  -subroutine assert(predicate, error_message)
  -
  -  implicit none
  -
  -  logical, intent(in) :: predicate
  -  character(*), intent(in) :: error_message
  -
  -  if (.not. predicate) then
  -    print *, 'ERROR:', error_message
  -    stop 1
  -  endif
  -
  -end subroutine
  -
  -
   program xc_example
   
   ! This example contains calls to f90 interface routines that are needed to "talk
  
positive.code:   ierr = xc_set(id, 'pbe', 1.0d0)
    call assert(ierr == 0, "functional name not recognized")
  
    !-----------------------------------------------------------------------------
    ! in the first example we compute the XC energy ("order 0 derivative")
  
  
negative.code:     density(2, 2, ipoint) = 6.0d0 ! first order
      density(2, 3, ipoint) = 7.0d0 ! first order
      density(2, 4, ipoint) = 8.0d0 ! first order
      density(3, 1, ipoint) = 5.0d0 ! first order
      density(3, 2, ipoint) = 6.0d0 ! first order
      density(3, 3, ipoint) = 7.0d0 ! first order
      density(3, 4, ipoint) = 8.0d0 ! first order
      density(4, 1, ipoint) = 0.0d0 ! first order
      density(4, 2, ipoint) = 0.0d0 ! first order
      density(4, 3, ipoint) = 0.0d0 ! first order
      density(4, 4, ipoint) = 0.0d0 ! first order
    end do
  
    res = derivative(id, &
                     num_points, &
                     num_variables, &
                     vector_length, &
                     density)
    deallocate(density)
  
    ! compare with reference
    call assert((abs(-9.4927931153398468d0 - res) < 1.0d-6), &
                "derivatives do not match reference numbers")
  
  
    !-----------------------------------------------------------------------------
    ! now we try 3nd order, contracted with perturbed densities
  
    order = 3
    ierr = xc_eval_setup(id, XC_N_NX_NY_NZ, XC_CONTRACTED, order)
    call assert(ierr == 0, "xc_eval_setup failed")
  
    vector_length = 2**order
    allocate(density(vector_length, num_variables, num_points))
    density = 0.0d0
    do ipoint = 1, num_points
      ! we use fantasy values here
      density(1, 1, ipoint) = 1.0d0   ! zeroth order
      density(1, 2, ipoint) = 2.0d0   ! zeroth order
      density(1, 3, ipoint) = 3.0d0   ! zeroth order
  