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: !
  !  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.
  
negative.code:   deallocate (input_array)
    deallocate (output_array)
    deallocate (res_reference)
  
    !  Notes:
  
  !  The computed second derivatives are now already contracted with the
  