anchor.diff:  
               ! Calculate microcanonical rate coefficients at the current conditions
               do w = 1, size(net%reactions)
  +                allocate(net%reactions(w)%kf(1:nGrains))
  +                allocate(net%reactions(w)%kb(1:nGrains))
                   call reaction_calculateMicrocanonicalRates(net%reactions(w), T, Elist, nGrains, net%isomers, net%species, &
                       nIsom, nReac, nProd)
               end do
  
positive.code:                 ! Determine phenomenological rate coefficients
                  call network_applyApproximateMethod(net, nIsom, nReac, nProd, Elist, nGrains, T, P, method, K(u,v,:,:))
  
              end do
  
          end do
  
negative.code:         integer, intent(in) :: method
          real(8), dimension(1:nT,1:nP,1:nIsom+nReac+nProd,1:nIsom+nReac+nProd), intent(out) :: K
  
          real(8), dimension(:), allocatable :: Elist
          real(8), dimension(1:nIsom+nReac+nProd,1:nGrains0) :: densStates0
          real(8) T, P
          integer i, j, u, v, w, nGrains
  
          ! Save the original densities of states for each isomer
          ! (We'll be overwriting the ones on each isomer at each temperature)
          do i = 1, nIsom
              densStates0(i,:) = net%isomers(i)%densStates
              deallocate(net%isomers(i)%densStates)
          end do
  
          do u = 1, nT
  
              T = Tlist(u)
  
              ! Choose energy grains
              call network_determineEnergyGrains(net, nIsom, grainSize, numGrains, T, Elist)
              nGrains = size(Elist)
              write (1,*) '    Using', nGrains, 'grains of size', Elist(2) - Elist(1), 'J/mol in range', minval(Elist), &
                  ' to', maxval(Elist), 'J/mol'
  
              do i = 1, nIsom
                  allocate( net%isomers(i)%densStates(1:nGrains) )
                  call network_mapDensityOfStates(net%isomers(i)%E0, densStates0(i,:), Elist0, nGrains0, &
                      net%isomers(i)%densStates, Elist, nGrains)
              end do
  
              ! Calculate the equilibrium (Boltzmann) distributions for each (unimolecular) well
              do i = 1, nIsom
  