anchor.diff:    deallocate (density)
     deallocate (input_array)
     deallocate (output_array)
  -  
  -  
  +
  +
   !  now third derivative of the "potential", contracted with perturbed densities.
     order = 3
     nr_variables = 4
  
positive.code:   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  
    
    
  !  now let's compute the first derivatives ('potential')
  !  and do NOT contract them!
  !  first set up xcfun for first derivatives
  
negative.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
  