anchor.diff:              prod = isomerList(rxn%reac)
   
               ! Calculate forward rate coefficient via inverse Laplace transform
  -            call rateILT(rxn%E0 - reac%E0, reac%densStates, &
  +            call reaction_kineticsILT(rxn%E0 - reac%E0, reac%densStates, &
                   rxn%arrhenius, T, E, rxn%kb)
   
               ! Calculate equilibrium constant
  -            Keq = equilCoeff(reac, prod, T, speciesList)
  +            Keq = reaction_getEquilibriumConstant(reac, prod, T, speciesList)
   
               ! Calculate backward rate coefficient via detailed balance
               ! assuming multimolecular isomer is fully thermalized
  
positive.code:         if (size(reac%species) == 1 .and. size(prod%species) == 1) then
  
              ! Calculate forward rate coefficient via inverse Laplace transform
              call rateILT(rxn%E0 - reac%E0, reac%densStates, &
                  rxn%arrhenius, T, E, rxn%kf)
  
              ! Calculate backward rate coefficient via detailed balance
              do r = 1, nGrains
                  if (prod%densStates(r) /= 0) then
                      rxn%kb(r) = rxn%kf(r) * reac%densStates(r) / prod%densStates(r)
                  end if
              end do
  
          elseif (size(reac%species) == 1 .and. size(prod%species) > 1) then
  
              ! Calculate forward rate coefficient via inverse Laplace transform
              call rateILT(rxn%E0 - reac%E0, reac%densStates, &
                  rxn%arrhenius, T, E, rxn%kf)
  
              ! Calculate equilibrium constant
              Keq = equilCoeff(reac, prod, T, speciesList)
  
              ! Calculate backward rate coefficient via detailed balance
              ! assuming multimolecular isomer is fully thermalized
  
negative.code:         ! Determine free energy of reaction
          dGrxn = 0.0
          do i = 1, size(reac%species)
              dGrxn = dGrxn - species_getFreeEnergy(speciesList(reac%species(i))%thermo, T)
          end do
          do i = 1, size(prod%species)
              dGrxn = dGrxn + species_getFreeEnergy(speciesList(prod%species(i))%thermo, T)
          end do
  
          ! Determine Ka
          Keq = exp(-dGrxn / 8.314472 / T)
          ! Determine Kc
          Keq = Keq * (100000.0 / 8.314472 / T)**(size(prod%species) - size(reac%species))
  
      end function
  
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
      subroutine network_shiftToZeroEnergy(net)
  
          type(Network), intent(inout) :: net
  
          real(8) Emin
          integer i, j
  
          Emin = 1.0e20
  
          do i = 1, size(net%isomers)
              net%isomers(i)%E0 = 0.0
              do j = 1, size(net%isomers(i)%species)
                  net%isomers(i)%E0 = net%isomers(i)%E0 + net%species(net%isomers(i)%species(j))%E0
              end do
              if (net%isomers(i)%E0 < Emin) Emin = net%isomers(i)%E0
          end do
          do i = 1, size(net%isomers)
              net%isomers(i)%E0 = net%isomers(i)%E0 - Emin
          end do
          do i = 1, size(net%reactions)
              net%reactions(i)%E0 = net%reactions(i)%E0 - Emin
          end do
  
      end subroutine
  
      subroutine network_getEnergyGrains(Emin, Emax, dE0, nGrains0, Elist)
  
          real(8), intent(in) :: Emin, Emax, dE0
          integer, intent(in) :: nGrains0
          real(8), dimension(:), allocatable, intent(inout) :: Elist
  
          real(8) dE
          integer nGrains, useGrainSize, r
  
          if (nGrains0 <= 0 .and. dE0 > 0.0) then
              useGrainSize = 0
          elseif (nGrains0 > 0 .and. dE0 <= 0.0) then
              useGrainSize = 1
          else
              ! Choose the tighter constraint
  			useGrainSize = 0
              dE = (Emax - Emin) / (nGrains0 - 1)
              if (dE > dE0) useGrainSize = 1
          end if
  
          if (useGrainSize == 1) then
              nGrains = nGrains0
  			dE = (Emax - Emin) / (nGrains0 - 1)
          else
              nGrains = int((Emax - Emin) / dE0) + 1
  			dE = dE0
          end if
  
          allocate(Elist(1:nGrains))
          do r = 1, nGrains
              Elist(r) = Emin + (r - 1) * dE
          end do
  
      end subroutine
  
  	subroutine network_determineEnergyGrains(net, nIsom, grainSize, nGrains, Tmax, Elist)
  
          type(Network), intent(in) :: net
          integer, intent(in) :: nIsom
          real(8), intent(in) :: grainSize
          integer, intent(in) :: nGrains
          real(8), intent(in) :: Tmax
          real(8), dimension(:), allocatable, intent(inout) :: Elist
  
          integer nE, done, maxIndex
          real(8) Emin, Emax0, Emax, dE, mult, maxValue, tol
          type(Isomer) isom
  
          integer i, r
  
  		! For the purposes of finding the maximum energy we will use 401 grains
  		nE = 401
  		dE = 0.0
  
  		! Determine minimum energy and isomer with minimum ground-state energy
  		isom = net%isomers(1)
          do i = 2, nIsom
              if (isom%E0 > net%isomers(i)%E0) isom = net%isomers(i)
          end do
  		Emin = floor(isom%E0)
  
  		! Determine maximum energy and isomer with maximum ground-state energy
  		isom = net%isomers(1)
          do i = 2, nIsom
              if (isom%E0 < net%isomers(i)%E0) isom = net%isomers(i)
          end do
  		Emax0 = isom%E0
  
  		! (Try to) purposely overestimate Emax using arbitrary multiplier
  		! This is to (hopefully) avoid multiple density of states calculations
  		mult = 50
  		tol = 1e-8
          done = 0
          do while (done == 0)
  
  			Emax = ceiling(Emax0 + mult * 8.314472 * Tmax)
  
  			call network_getEnergyGrains(Emin, Emax, dE, nGrains, Elist)
  			call isomer_getDensityOfStates(net, isom, Elist, nGrains)
  			call isomer_getEqDist(isom, Elist, nGrains, Tmax)
  
  			! Find maximum of distribution
  			maxIndex = 0
  			maxValue = 0.0
  			do r = 1, nGrains
  				if (isom%eqDist(r) > maxValue) then
  					maxValue = isom%eqDist(r)
  					maxIndex = r
                  end if
              end do
  
  			! If tail of distribution is much lower than the maximum, then we've found bounds for Emax
  			if (isom%eqDist(nGrains) / maxValue < tol) then
  				r = nGrains - 1
  				do while (r > 1 .and. done == 0)
  					if (isom%eqDist(r) / maxValue > tol) then
  						done = 1
  					else
  						r = r - 1
                      end if
                  end do
  
                  ! Add difference between isomer ground-state energy and highest
                  ! transition state or isomer energy
                  maxValue = 0.0
                  do r = 1, size(net%isomers)
                      if (net%isomers(r)%E0 > maxValue) maxValue = net%isomers(r)%E0
                  end do
                  do r = 1, size(net%reactions)
                      if (net%reactions(r)%E0 > maxValue) maxValue = net%reactions(r)%E0
                  end do
                  Emax = Elist(r) + maxValue - Emax0
  
                  ! Round Emax up to nearest integer
                  Emax = ceiling(Emax)
  
  			else
  				mult = mult + 50
              end if
  
          end do
  
  		! Return the chosen energy grains
  		call network_getEnergyGrains(Emin, Emax, dE, nE, Elist)
  
      end subroutine
  
      subroutine network_calculateRateCoefficients(net, nIsom, nReac, nProd, &
          Elist, nGrains, Tlist, nT, Plist, nP, method, K)
  
          integer, intent(in) :: nIsom, nReac, nProd, nGrains, nT, nP
          type(Network), intent(inout) :: net
          real(8), dimension(1:nGrains), intent(in) :: Elist
          real(8), dimension(1:nT), intent(in) :: Tlist
          real(8), dimension(1:nP), intent(in) :: Plist
          integer, intent(in) :: method
          real(8), dimension(1:nT,1:nP,1:nIsom+nReac+nProd,1:nIsom+nReac+nProd), intent(out) :: K
  
          real(8) T, P
          integer i, j, u, v, w
  
          do u = 1, nT
  
              T = Tlist(u)
  
              ! Calculate the equilibrium (Boltzmann) distributions for each (unimolecular) well
              do i = 1, nIsom
                  call isomer_getEqDist(net%isomers(i), Elist, nGrains, T)
              end do
  
              ! 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)
              end do
  
              do v = 1, nP
  
                  P = Plist(v)
  
                  ! Calculate collision frequencies
                  do i = 1, nIsom
                      net%isomers(i)%collFreq = isomer_getCollisionFrequency(T, P, net%bathGas, &
                          net%species(net%isomers(i)%species(1))%general)
                  end do
  
  