anchor.diff:          eqDist(index,i+nIsom) = 1.0
       end do
   
  -    ! Construct the rate constant matrix
  -    do j = 1, nIsom+nReac
  -        do r = 1, nGrains
  -            do i = 1, nIsom
  -                index = indices(r,i)
  -                if (index > 0) then
  -                    K(i,j) = K(i,j) + sum(M(index,:) * eqDist(:,j))
  -                end if
  -            end do
  -        end do
  -        do i = 1, nReac
  -            index = nRows - nReac + i
  -            K(i+nIsom,j) = sum(M(index,:) * eqDist(:,j))
  -        end do
  +    ! Calculate the phenomenological rate constants
  +    ! This version is intended to follow the notation of Miller and Klippenstein
  +    allocate( dXij(1:nIsom+nReac+nProd, 1:nCSE), C(1:nRows) )
  +    do n = 1, nIsom+nReac
  +        do j = 1, nCSE
  +
  +            do i = 1, nIsom+nReac+nProd
  +                dXij(i,j) = 0.0
  +            end do
  +
  +            do r = 1, nRows
  +                C(r) = X(r,j) * sum(Xinv(j,:) * eqDist(:,n))
  +            end do
  +
  +            do r = 1, nGrains
  +                do i = 1, nIsom
  +                    index = indices(r,i)
  +                    if (index > 0) then
  +                        ! Isomers
  +                        dXij(i,j) = dXij(i,j) + C(index)
  +                        ! Products
  +                        do h = nReac+1, nReac+nProd
  +                            dXij(nIsom+h,j) = dXij(nIsom+h,j) + C(index) * Gnj(h,i,r) / V(j)
  +                        end do
  +                    end if
  +                end do
  +            end do
  +            do i = 1, nReac
  +                index = nRows - nReac + i
  +                dXij(nIsom+i,j) = dXij(nIsom+i,j) + C(index)
  +            end do
  +
  +            ! The isomers and reactant
  +            do i = 1, nIsom+nReac+nProd
  +                K(i,n) = K(i,n) + V(j) * dXij(i,j)
  +            end do
  +
  +        end do
  +
       end do
   
       ! Clean up
  -    deallocate( Mcoll, Mrxn, M, S, Sinv, V0, work, eqDist, X, V, Xinv )
  +    deallocate( Mcoll, Mrxn, M, S, Sinv, V0, work, eqDist, X, V, Xinv, dXij, C )
   
   end subroutine
   
  
positive.code: 
      ! The size of the master equation matrix
      integer nRows
      
      ! The number of chemically-significant eigenvalues
      integer nCSE
  
  
negative.code: !   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  !   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  !   DEALINGS IN THE SOFTWARE.
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
  subroutine estimateRateCoefficients_CSE(T, P, E, Mcoll0, E0, densStates, &
  