anchor.diff:      density(2, 2, ipoint) = 6.0d0 ! first order
       density(2, 3, ipoint) = 7.0d0 ! first order
       density(2, 4, ipoint) = 8.0d0 ! first order
  -  end do
  -
  -  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.
  -  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
  -
  -    do ipoint = 1, num_points
  -      input_array(:, ipoint) = (/density(:, 1, ipoint), &
  -                                 density(:, 2, ipoint), &
  -                                 density(:, 3, ipoint), &
  -                                 density(:, 4, ipoint)/)
  -    end do
  -    call xc_eval(id, num_points, input_array, output_array)
  -
  -    res = output_array(vector_length, 1)
  -    print *, 'd^2 E/dx_i dD_pert for i =', ideriv, 'is', res
  -
  -    ! compare with reference
  -    call assert((abs(res - res_reference(ideriv)) < 1.0d-6), &
  -                "derivatives do not match reference numbers")
  -  end do
  -
  -  deallocate (density)
  -  deallocate (input_array)
  -  deallocate (output_array)
  -  deallocate (res_reference)
  -
  -!  The computed second derivatives are now already contracted with the
  -!  perturbed density! This works to any order.  In this example we put in ones and
  -!  zeros, which you need to do to compute matrix elements. You can also put in a
  -!  perturbed density instead of ones and zeros, then you generate a response
  -!  contribution d^2E/dD1dD2. Typically this is what you want to do to high order,
  -!  because by the 2N+1 rule you don't have to compute matrix elements to very high
  -!  order. Then you need only one call of xc_eval, not one for each variable. This
  -!  is very efficient, much better than computing all partial derivatives. For
  -!  second derivatives it might be better to compute all partial derivatives and
  -!  reuse them in the response solver.
  -!
  -!  In this example we use the variables that are linear in the density matrix.
  -!  This makes life easier because you can construct them trivially from perturbed
  -!  density matrices. If you use i.e. the square norm of the density gradient
  -!  things get more complicated.
  -!
  -!  For the spin polarized cases (using density and spin density or up- and down
  -!  densities) we would use 8 variables instead of 4.
  -!
  -!  This will be exemplified in the following.
  -
  -
  -
  -  ! Now contract the 2nd derivative completely
  -  !  allocate density and fill in some fantasy values
  -  !  Now the shape of zdensity is different from the last example.
  -  !  Typically, we want to contract the first kernel with e.g.
  -  !  first and/or second-order perturbed densities and again, we have
  -  !  2^order, i.e. 4 elements:
  -  !  1.) ground state density (zeroth order)
  -  !  2.) perturbed density a) (first order)
  -  !  3.) perturbed density b) (first order)
  -  !  4.) perturbed density ab) (second order)
  -  allocate (density(vector_length, num_variables, num_points))
  -  density = 0.0d0
  -  do ipoint = 1, num_points
  -    density(1, 1, ipoint) = 1.0d0    ! zeroth order
  -    density(1, 2, ipoint) = 2.0d0    ! zeroth order
  -    density(1, 3, ipoint) = 3.0d0    ! zeroth order
  -    density(1, 4, ipoint) = 4.0d0    ! zeroth order
  -    density(2, 1, ipoint) = 5.0d0    ! first order
  -    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
  -
  -  allocate (input_array(num_variables*vector_length, num_points))
  -  allocate (output_array(vector_length, num_points))
  -
  -  do ipoint = 1, num_points
  -    input_array(:, ipoint) = (/density(:, 1, ipoint), &
  -                               density(:, 2, ipoint), &
  -                               density(:, 3, ipoint), &
  -                               density(:, 4, ipoint)/)
  -  end do
  -  call xc_eval(id, num_points, input_array, output_array)
  -
  -  !  The output array will then contain:
  -  ! 1.) Product of functional with ground state density
  -  ! 2.) Product of first-order functional derivative with first-order density a)
  -  ! 3.) Product of first-order functional derivative with first-order density b)
  -  ! 4.) Product of first-order functional derivative with second-order density ab) +
  -  !     Product of second-order functional derivative with first-order densities a) and b)
  -
  -  deallocate (density)
  -  deallocate (input_array)
  -  deallocate (output_array)
  -
  -  !  now third derivative of the "potential", contracted with perturbed densities.
  +    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 and fill in some fantasy values
  -  allocate (density(vector_length, num_variables, num_points))
  -  density = 0.0d0
  -  do ipoint = 1, num_points
  +  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
  
positive.code: !
  ! 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
  
negative.code:   integer, parameter :: num_variables = 4
  
    character(1000) :: text
    integer :: id, order, ierr, ipoint
    integer :: vector_length
    real(8) :: res
  
    real(8), allocatable :: density(:, :, :)
  
  
    ! print some info and copyright about the library
  