anchor.diff:    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
  
positive.code:                                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)
  
negative.code: !
  !  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.
  