anchor.diff:  
           elseif (rxn%reac <= nIsom .and. rxn%prod > nIsom) then
   
  -            ! Calculate forward rate coefficient via inverse Laplace transform
  -            call reaction_kineticsILT(rxn%E0 - reac%E0, reac%densStates, &
  -                rxn%arrhenius, T, E, rxn%kf)
  -
               ! Calculate equilibrium constant
               Keq = reaction_getEquilibriumConstant(reac, prod, T, speciesList)
   
  
positive.code:             write (1,*) '    Expected kf at', T, 'K =', kf0 
              write (1,*) '      Actual kf at', T, 'K =', kf
          end if
          
          !write (*,*) rxn%reac, rxn%prod, rxn%kf(nGrains), rxn%kb(nGrains)
  
      end subroutine
  
negative.code: 
          elseif (rxn%reac > nIsom .and. rxn%prod <= nIsom) then
  
              ! Calculate equilibrium constant
              Keq = reaction_getEquilibriumConstant(reac, prod, T, speciesList)
  
              ! Calculate backward rate coefficient via detailed balance
              do r = 1, nGrains
                  if (prod%densStates(r) /= 0) then
                      rxn%kb(r) = rxn%kf(r) / Keq * &
                          (reac%densStates(r) / reac%Q) / (prod%densStates(r) / prod%Q)
                  end if
              end do
  
              ! Include equilibrium distribution in forward direction
              do r = 1, nGrains
                  rxn%kf(r) = rxn%kf(r) * reac%densStates(r) * exp(-E(r) / 8.314472 / T) / reac%Q * dE
              end do
  
              ! Check that forward and reverse rates integrate to give the proper k(E) values
              kf0 = rxn%arrhenius%A * T ** rxn%arrhenius%n * exp(-rxn%arrhenius%Ea / 8.314472 / T)
              kf = sum(rxn%kf)
              kb = sum(rxn%kb * prod%densStates * exp(-E / 8.314472 / T) / prod%Q * dE)           
              Keq0 = reaction_getEquilibriumConstant(reac, prod, T, speciesList)
              Keq = kf / kb  
  
          end if
          
  