anchor.diff:    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
  +  nr_variables = 4
  +  ierr = xc_eval_setup(id, XC_N_NX_NY_NZ, XC_CONTRACTED, order)
  +  if (ierr /= 0) then
  +    print *, 'xc_eval_setup failed with error ', ierr
  +    stop 1
  +  endif
  +
  +  vector_length = 2**order
  +!  allocate density and fill in some phantasy values
  +  allocate (density(vector_length, nr_variables, NR_POINTS))
  +  density = 0.0d0
  +  do ipoint = 1, NR_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
  +    density(1, 4, ipoint) = 4.0d0 !nabla_z n zeroth order
  +    density(2, 1, ipoint) = 5.0d0 !        n first order
  +    density(2, 2, ipoint) = 6.0d0 !nabla_x n first order
  +    density(2, 3, ipoint) = 7.0d0 !nabla_y n first order
  +    density(2, 4, ipoint) = 8.0d0 !nabla_z n first order
  +  end do
  +
  +  allocate (input_array(nr_variables*vector_length, NR_POINTS))
  +  allocate (output_array(vector_length, NR_POINTS))
  +
  +!     It is possible to put in other numbers than 1 and 0
  +!     this is where perturbed densities go for automatic contraction
  +
  +    do ipoint = 1, NR_POINTS
  +      input_array(:, ipoint) = (/density(:, 1, ipoint), &
  +                                 density(:, 2, ipoint), &
  +                                 density(:, 3, ipoint), &
  +                                 density(:, 4, ipoint)/)
  +    end do
  +    call xc_eval(id, NR_POINTS, input_array, output_array)
  +
  +  deallocate (density)
  +  deallocate (input_array)
  +  deallocate (output_array)
  +  
  +  
   !  reference values for self test
     allocate (res_reference(nr_variables))
     res_reference(1) = -0.97182634532897016d0
     res_reference(2) = -8.9803044914016916d-3
     res_reference(3) = -1.3470456737102541d-2
  -  res_reference(4) = -1.7960608982803383d-2
  -
  +  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
     nr_variables = 4
  
positive.code:   res_reference(2) = 2.1489169824673575d-2
    res_reference(3) = 4.1214059228412030d-2
    res_reference(4) = 6.0938948632150471d-2
  
  !  now second derivative of the "potential", contracted with one perturbed density.
  !  hopefully the strange input/output_array format will start to make sense
    order = 2
    nr_variables = 4
  
negative.code: !  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 phantasy 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, nr_variables, NR_POINTS))
    density = 0.0d0
    do ipoint = 1, NR_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(nr_variables*vector_length, NR_POINTS))
    allocate (output_array(vector_length, NR_POINTS))
  
    do ipoint = 1, NR_POINTS
      input_array(:, ipoint) = (/density(:, 1, ipoint), &
                                 density(:, 2, ipoint), &
                                 density(:, 3, ipoint), &
                                 density(:, 4, ipoint)/)
    end do
    call xc_eval(id, NR_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.
    order = 3
    nr_variables = 4
    ierr = xc_eval_setup(id, XC_N_NX_NY_NZ, XC_CONTRACTED, order)
    if (ierr /= 0) then
      print *, 'xc_eval_setup failed with error ', ierr
      stop 1
    endif
  
    vector_length = 2**order
  !  allocate density and fill in some phantasy values
    allocate (density(vector_length, nr_variables, NR_POINTS))
    density = 0.0d0
    do ipoint = 1, NR_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 (1)
      density(2, 2, ipoint) = 6.0d0   ! first order (1)
      density(2, 3, ipoint) = 7.0d0   ! first order (1)
      density(2, 4, ipoint) = 8.0d0   ! first order (1)
      density(3, 1, ipoint) = 9.0d0   ! first order (2)
      density(3, 2, ipoint) = 10.0d0  ! first order (2)
      density(3, 3, ipoint) = 11.0d0  ! first order (2)
      density(3, 4, ipoint) = 12.0d0  ! first order (2)
      density(4, 1, ipoint) = 5.0d0   ! second order (depending on (1) and (2))
      density(4, 2, ipoint) = 6.0d0   ! second order (depending on (1) and (2))
      density(4, 3, ipoint) = 7.0d0   ! second order (depending on (1) and (2))
      density(4, 4, ipoint) = 8.0d0   ! second order (depending on (1) and (2))
      density(5, 1, ipoint) = 9.0d0   ! first order (3)
      density(5, 2, ipoint) = 10.0d0  ! first order (3)
      density(5, 3, ipoint) = 11.0d0  ! first order (3)
      density(5, 4, ipoint) = 12.0d0  ! first order (3)
      density(6, 1, ipoint) = 5.0d0   ! second order (depending on (1) and (3))
      density(6, 2, ipoint) = 6.0d0   ! second order (depending on (1) and (3))
      density(6, 3, ipoint) = 7.0d0   ! second order (depending on (1) and (3))
      density(6, 4, ipoint) = 8.0d0   ! second order (depending on (1) and (3))
      density(7, 1, ipoint) = 9.0d0   ! second order (depending on (2) and (3))
      density(7, 2, ipoint) = 10.0d0  ! second order (depending on (2) and (3))
      density(7, 3, ipoint) = 11.0d0  ! second order (depending on (2) and (3))
      density(7, 4, ipoint) = 12.0d0  ! second order (depending on (2) and (3))
      density(8, 1, ipoint) = 0.0d0   ! third order (depending on (1), (2) and (3))
      density(8, 2, ipoint) = 0.0d0   ! third order (depending on (1), (2) and (3))
      density(8, 3, ipoint) = 0.0d0   ! third order (depending on (1), (2) and (3))
      density(8, 4, ipoint) = 0.0d0   ! third order (depending on (1), (2) and (3))
    end do
  
    allocate (input_array(nr_variables*vector_length, NR_POINTS))
    allocate (output_array(vector_length, NR_POINTS))
  
      do ipoint = 1, NR_POINTS
        input_array(:, ipoint) = (/density(:, 1, ipoint), &
                                   density(:, 2, ipoint), &
                                   density(:, 3, ipoint), &
                                   density(:, 4, ipoint)/)
      end do
      
      call xc_eval(id, NR_POINTS, input_array, output_array)
  
      
    deallocate (density)
    deallocate (input_array)
    deallocate (output_array) 
    
    write(*,*)'Kernel test has ended properly!'
    
  
  end program
  