anchor.diff:  
           ! Active-state energy of each isomer
           do i = 1, nIsom+nReac+nProd
  -            Eres(i) = isomer_getActiveSpaceEnergy(net%isomers(i), net%reactions)
  +            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)%reactant
  -            j = net%reactions(r)%product
  -            if (i <= nIsom .and j <= nIsom) then ! Isomerization
  +            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
  +            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
  +            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_calculateCollisionEfficiency(net%isomers(i), T, net%reactions, net%dEdown, Elist)
  +                    isomer_getCollisionEfficiency(Elist, nGrains, T, i, &
  +                    net%isomers(i)%E0, net%isomers(i)%densStates, net%reactions, dEdown)
               end do
   
               ! Apply modified strong collision method
  
positive.code:         ! Determine free energy of reaction
          dGrxn = 0.0
          do i = 1, size(reac%species)
              dGrxn = dGrxn - freeEnergy(speciesList(reac%species(i))%thermo, T)
          end do
          do i = 1, size(prod%species)
              dGrxn = dGrxn + freeEnergy(speciesList(prod%species(i))%thermo, T)
          end do
  
          ! Determine Ka
          reaction_equilCoeff = exp(-dGrxn / 8.314472 / T)
          ! Determine Kc
          reaction_equilCoeff = reaction_equilCoeff * (100000.0 / 8.314472 / T)**(size(prod%species) - size(reac%species))
  
      end function
  
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
      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(in) :: 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
          character(len=*), 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, size(net%Tlist)
  
              T = Tlist(u)
  
              ! Calculate the equilibrium (Boltzmann) distributions for each (unimolecular) well
              do i = 1, nIsom
                  eqDist(i,:) = isomer_eqDist(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), Elist, nGrains, T, net%species)
              end do
  
              do v = 1, size(net%Plist)
  
                  P = Plist(v)
  
                  ! Calculate collision frequencies
                  do i = 1, nIsom
                      net%isomers(i)%collFreq = collisionFrequency(T, P, net%bathGas, &
                          net%species(net%isomers(i)%species(1))%general)
                  end do
  
  
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, &
  