anchor.diff:      subroutine network_applyApproximateMethod(net, nIsom, nReac, nProd, &
           Elist, nGrains, T, P, method, K)
   
  -        integer, intent(in) :: nIsom, nReac, nProd, nGrains, nT, nP
  -        type(Network), intent(in) :: net
  +        integer, intent(in) :: nIsom, nReac, nProd, nGrains
  +        type(Network), intent(inout) :: net
           real(8), dimension(1:nGrains), intent(in) :: Elist
           real(8), intent(in) :: T, P
  -        character(len=*), intent(in) :: method
  +        integer, intent(in) :: method
           real(8), dimension(1:nIsom+nReac+nProd,1:nIsom+nReac+nProd), intent(out) :: K
   
  -        real(8), dimension(1:nIsom,1:nGrains) densStates
  -        real(8), dimension(1:nIsom+nReac+nProd) Eres
  -        real(8), dimension(1:nIsom,1:nIsom,1:nGrains) Kij
  -        real(8), dimension(1:nReac+nProd,1:nIsom,1:nGrains) Gnj
  -        real(8), dimension(1:nIsom,1:nReac,1:nGrains) Fim
  -
  -        real(8), dimension(1:nIsom) collFreq
  -
  +        real(8), dimension(1:nGrains) :: densStates0
  +        real(8), dimension(1:nIsom,1:nGrains) :: densStates
  +        real(8), dimension(1:nIsom+nReac+nProd) :: Eres
  +        real(8), dimension(1:nIsom) :: E0
  +        real(8), dimension(1:nIsom,1:nGrains,1:nGrains) :: Mcoll
  +        real(8), dimension(1:nIsom,1:nIsom,1:nGrains) :: Kij
  +        real(8), dimension(1:nReac+nProd,1:nIsom,1:nGrains) :: Gnj
  +        real(8), dimension(1:nIsom,1:nReac,1:nGrains) :: Fim
  +
  +        real(8), dimension(1:nIsom) :: collFreq
  +
  +        real(8) dE, dEdown
  +        character(len=128) :: msg
  +
  +        integer invalidRate
  +
  +        integer i, j, r
   
           write(1,*) 'Applying method at', T, 'K,', P/1e5, 'bar...'
   
  -        dE = Elist[1] - Elist[0]
  +        dE = Elist(2) - Elist(1)
   
           ! Density of states per partition function (i.e. normalized density of
           ! states with respect to Boltzmann weighting factor) for each isomer
  
positive.code: 
          ! Active-state energy of each isomer
          do i = 1, nIsom+nReac+nProd
              Eres(i) = isomer_getActiveSpaceEnergy(net%isomers(i), net%reactions)
          end do
  
          ! Isomerization, dissociation, and association microcanonical rate
          ! coefficients, respectively
          do r = 1, size(net%reactions)
              i = net%reactions(r)%reactant
              j = net%reactions(r)%product
              if (i <= nIsom .and j <= nIsom) then ! Isomerization
                  Kij(j,i,:) = net%reactions(r)%kf
                  Kij(i,j,:) = net%reactions(r)%kb
              elseif (i <= nIsom .and j > nIsom) then ! Dissociation
                  Gnj(j-nIsom,i,:) = net%reactions(r)%kf
                  if (j - nIsom <= nReac) Fim(i,j-nIsom,:) = net%reactions(r)%kb
              elseif (i > nIsom .and j <= nIsom) then ! Association
                  if (i - nIsom <= nReac) Fim(j,i-nIsom,:) = net%reactions(r)%kf
                  Gnj(i-nIsom,j,:) = net%reactions(r)%kb
              end if
          end do
  
  
          if (method == 1) then ! Modified strong collision
  
              ! Modified collision frequency of each isomer
              do i = 1, nIsom
                  collFreq(i) = net%isomers(i)%collFreq * &
                      isomer_calculateCollisionEfficiency(net%isomers(i), T, net%reactions, net%dEdown, Elist)
              end do
  
              ! Apply modified strong collision method
  
negative.code:                 Kij, Fim, Gnj, nIsom, nReac, nProd, nGrains, K, msg)
  
          elseif (method == 2) then ! Reservoir state
  
              ! Ground-state energy for each isomer
              do i = 1, nIsom
                  E0(i) = net%isomers(i)%E0
                  collFreq(i) = net%isomers(i)%collFreq
                  densStates0 = net%isomers(i)%densStates
                  call collisionmatrix(T, P, Elist, collFreq, densStates0, E0(i), dEdown, Mcoll(i,:,:), msg)
              end do
  
              call estimateratecoefficients_rs(T, P, Elist, Mcoll, densStates, E0, Eres, &
  