anchor.diff:          real(dp) :: qD, qi, ktf, beta, qs, wf, sf, L2, eta_n_0, G_k, Tp, eta_n
           real(dp) :: isospin, n_0, n_2, n_l, n_in
           real(dp), dimension(:), pointer :: y
  -        integer, intent(out) :: ierr
           integer ,dimension(:), pointer :: iwork => null()
           real(dp), dimension(:), pointer :: work => null()
  -        real(dp), dimension(1) :: rtol, atol
  +        real(dp), dimension(num_integration_variables) :: rtol, atol
           integer, dimension(:), pointer :: ipar => null()
           real(dp), dimension(:), pointer :: rpar => null()
           real(dp) :: h, kn_end, kn_start, kn, kfi
  -        integer :: liwork, lwork, itol, lipar, lrpar
  -        integer :: n, idid, lout, iout, npts, num_integration_variables, num_ipar, num_rpar
  -
  -        allocate(y(1))
  -        num_ipar = 0
  -        num_rpar = 6    
  +        integer :: liwork, lwork, itol, lipar, lrpar,idid, lout, iout, npts
  +
  +        allocate(y(num_integration_variables))
  +
           kFn = (threepisquare*nn)**onethird
           kfi = (threepisquare*nn/ionic%Yn*(1.0-ionic%Yn)/ionic%A)**onethird
  -        num_integration_variables = 1
  +
           call dop853_work_sizes( &
           &   num_integration_variables,num_integration_variables,liwork,lwork)
           allocate(iwork(liwork), work(lwork),ipar(num_ipar), rpar(num_rpar))
  
positive.code:         V = hbar**2*(threepisquare*n_in)**(2.0/3.0)/2.0/Mneutron &
          &   *(1.0 - (nn/n_in)**(2.0/3.0))
      end function neutron_potential
  
  
      subroutine structure_factor_phonon(n,x,h,y,dy,lrpar,rpar,lipar,ipar,ierr)
          ! integrand for eq. (16), Deibel et al. (2017)
  
negative.code:         rpar(5) = kfn
          rpar(6) = R_a
  
          call dop853(num_integration_variables, &
          &   structure_factor_phonon,kn_start,y,kn_end,h, &
          &   int_default_max_step_size,int_default_max_steps, &
          &   rtol,atol,itol, null_solout, iout, work, lwork, iwork, liwork,  &
          &   num_rpar, rpar, num_ipar, ipar, lout, idid)
  