anchor.diff:  
       end subroutine
   
  +    function isomer_getActiveSpaceEnergy(isom, reactions) result(Eres)
  +
  +        integer, intent(in) :: isom
  +        type(Reaction), dimension(:), intent(in) :: reactions
  +        real(8) :: Eres
  +
  +        integer i, r
  +
  +        Eres = 1.0e20
  +
  +        do r = 1, size(reactions)
  +            if (reactions(r)%reac == isom .or. reactions(r)%prod == isom) then
  +                if (reactions(r)%E0 < Eres) Eres = reactions(r)%E0
  +            end if
  +        end do
  +
  +    end function
  +
  +    subroutine isomer_getDensityOfStates(net, isom, Elist, nGrains)
  +
  +        integer, intent(in) :: nGrains
  +        type(Network), intent(in) :: net
  +        type(Isomer), intent(inout) :: isom
  +        real(8), dimension(1:nGrains), intent(in) :: Elist
  +
  +        real(8), dimension(1:nGrains) :: densStates0, densStates
  +
  +        integer index, i, r
  +
  +        ! Initialize density of states arrays
  +        do r = 1, nGrains
  +            isom%densStates(r) = 0.0
  +            densStates(r) = 0.0
  +            densStates0(r) = 0.0
  +        end do
  +
  +        ! Calculate the density of states for each species, convolving when
  +		! multiple species are present
  +        do i = 1, size(isom%species)
  +            call species_getDensityOfStates(net%species(isom%species(i)), Elist, nGrains, densStates0)
  +            call convolve(densStates, densStates0, Elist, nGrains)
  +        end do
  +
  +		! Shift to appropriate energy grain using E0
  +		index = 0
  +        do r = 1, nGrains
  +            if (isom%E0 < Elist(r) .and. index == 0) index = r
  +        end do
  +        isom%densStates(index:nGrains) = densStates(1:nGrains-index)
  +
  +    end subroutine
  +
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   
  -    subroutine reaction_calculateMicrocanonicalRates(rxn, T, E, isomerList, speciesList)
  -
  +    subroutine reaction_calculateMicrocanonicalRates(rxn, T, E, nE, isomerList, speciesList)
  +
  +        integer, intent(in) :: nE
           type(Reaction), intent(inout) :: rxn
           real(8), intent(in) :: T
  -        real(8), dimension(:), intent(in) :: E
  +        real(8), dimension(1:nE), intent(in) :: E
           type(Isomer), dimension(:), intent(in) :: isomerList
           type(Species), dimension(:), intent(in) :: speciesList
   
  
positive.code: 
      end function
  
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
      function isomer_collisionFrequency(T, P, bathGas, spec)
  
          ! Provide parameter type-checking
          real(8), intent(in)	:: T
          real(8), intent(in)	:: P
          type(GeneralData), intent(in) :: bathGas
          type(GeneralData), intent(in) :: spec
          real(8) :: collisionFrequency
  
          real(8)		::	kB, collisionIntegral, sigma, eps, mu, gasConc
  
  
negative.code: 
          ! Active-state energy of each isomer
          do i = 1, nIsom+nReac+nProd
              Eres(i) = isomer_getActiveSpaceEnergy(i, net%reactions)
          end do
  
          ! Isomerization, dissociation, and association microcanonical rate
          ! coefficients, respectively
          do r = 1, size(net%reactions)
              i = net%reactions(r)%reac
              j = net%reactions(r)%prod
              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
  
  
          ! Average energy transferred in a deactivating collision
          dEdown = net%bathGas%dEdown
  
          if (method == 1) then ! Modified strong collision
  
              ! Modified collision frequency of each isomer
              do i = 1, nIsom
                  collFreq(i) = net%isomers(i)%collFreq * &
                      isomer_getCollisionEfficiency(Elist, nGrains, T, i, &
                      net%isomers(i)%E0, net%isomers(i)%densStates, net%reactions, dEdown)
              end do
  
              ! Apply modified strong collision method
  