anchor.diff:                  ! Determine phenomenological rate coefficients
                   call network_applyApproximateMethod(net, nIsom, nReac, nProd, Elist, nGrains, T, P, method, K(u,v,:,:))
   
  +            end do
  +
  +            do i = 1, nIsom
  +                deallocate(net%isomers(i)%densStates)
  +                deallocate(net%isomers(i)%eqDist)
  +            end do
  +            do w = 1, size(net%reactions)
  +                deallocate(net%reactions(w)%kf)
  +                deallocate(net%reactions(w)%kb)
               end do
   
           end do
  
positive.code: 
              ! Calculate microcanonical rate coefficients at the current conditions
              do w = 1, size(net%reactions)
                  call reaction_calculateMicrocanonicalRates(net%reactions(w), T, Elist, nGrains, net%isomers, net%species, &
                      nIsom, nReac, nProd)
              end do
  
negative.code:         ! Return the chosen energy grains
          call network_getEnergyGrains(Emin, Emax, grainSize, nGrains, Elist)
  
      end subroutine
  
      subroutine network_mapDensityOfStates(E0, densStates0, Elist0, nGrains0, densStates, Elist, nGrains)
  
          real(8), intent(in) :: E0        
          integer, intent(in) :: nGrains0
          real(8), dimension(1:nGrains0), intent(in) :: densStates0
          real(8), dimension(1:nGrains0), intent(in) :: Elist0
          integer, intent(in) :: nGrains
          real(8), dimension(1:nGrains), intent(out) :: densStates
          real(8), dimension(1:nGrains), intent(in) :: Elist
  
          integer :: r, s
  
          do r = 1, nGrains
              densStates(r) = 0.
          end do
  
          do r = 1, nGrains
              if (Elist(r) >= E0) then
                  do s = 2, nGrains0
                      if (E0 + Elist0(s) >= Elist(r)) then
                          if (densStates0(s-1) > 0 .and. densStates0(s) > 0) then
                              densStates(r) = densStates0(s) * (densStates0(s-1) / densStates0(s)) ** &
                                  ((Elist(r) - E0 - Elist0(s)) / (Elist0(s-1) - Elist0(s)))
                          else
                              densStates(r) = densStates0(s) + (densStates0(s-1) - densStates0(s)) * &
                                  (Elist(r) - E0 - Elist0(s)) / (Elist0(s-1) - Elist0(s))
                          end if
                          exit
                      end if
                  end do
              end if
          end do
  
      end subroutine
  
  