anchor.diff:      integer, dimension(1:nGrains,1:nIsom) :: indices
   
       ! Intermediate matrices
  -    real(8), dimension(:,:), allocatable :: Mcoll, Mrxn, M, X, Xinv
  +    real(8), dimension(:,:), allocatable :: Mcoll, Mrxn, M, X, Xinv, eqDist
       real(8), dimension(:), allocatable :: S, Sinv, V0, V
   
       ! BLAS and LAPACK temporary variables
  -    real(8), dimension(:), allocatable :: work, iPiv
  +    real(8), dimension(:), allocatable :: work
       integer info
   
       ! Indices
  
positive.code: !        msg = 'Illegal argument passed to DSYEVX.'
  !        return
  !    end if
      
      ! Count the number of chemically distinct eigenvalues
      ! This will be ideally be nIsom+nProd, but might be lower if one or more
      ! chemical eigenmodes is blended with the internal energy eigenmodes
      nCSE = 1
      do i = 1, nIsom+nProd-1
          if (abs(V0(nRows-nIsom-nProd) / V0(nRows-i)) > 3.0) then
              nCSE = nCSE + 1
          end if
      end do
  
negative.code:         msg = 'Zero eigenvalue not found.'
          return
      end if
      
      ! Force largest eigenvalue to be zero (might be small numerical artifact)
      !V0(nRows) = 0.0
  
      allocate( X(1:nRows, 1:nCSE), V(1:nCSE), Xinv(1:nCSE, 1:nRows) )
      do j = 1, nCSE
          ! Find index in full matrix
          count = nRows - (nIsom + nProd) + j
          ! Eigenvalue matrix V
          V(j) = V0(count)
          ! Eigenvector matrix and its inverse/transpose
          do r = 1, nGrains
              do i = 1, nIsom
                  index = indices(r,i)
                  if (index > 0) then
                      ! Unsymmetrize as we go
                      X(index,j) = S(index) * M(index,count)
                      Xinv(j,index) = M(index,count) * Sinv(index)
                  end if
              end do
          end do
          do i = 1, nProd
              index = nRows - nProd + i
              ! Unsymmetrize as we go
              X(index,j) = S(index) * M(index,count)
              Xinv(j,index) = M(index,count) * Sinv(index)
          end do
      end do
  
      ! Generate new M containing only chemically-significant eigenmodes via
      ! multiplication of the above matrices
      do i = 1, nRows
          do j = 1, nRows
              M(i,j) = sum(X(i,:) * (V * Xinv(:,j)))
          end do
      end do
  
      ! Generate set of initial condition vectors (one per isomer and product)
      allocate( eqDist(1:nRows, 1:nIsom+nProd) )
      do j = 1, nIsom + nProd
          do r = 1, nRows
              eqDist(r,j) = 0.0
          end do
      end do
      do r = 1, nGrains
          do i = 1, nIsom
              index = indices(r,i)
              if (index > 0) then
                  eqDist(index,i) = densStates(i,r) * exp(-E(r) / 8.314472 / T)
              end if
          end do
      end do
      do i = 1, nProd
          index = nRows - nProd + i
          eqDist(index,i+nIsom) = 1.0
      end do
  
      ! Construct the rate constant matrix
      do j = 1, nIsom+nProd
          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, nProd
              index = nRows - nProd + i
              K(i+nIsom,j) = sum(M(index,:) * eqDist(:,j))
          end do
      end do
  
      ! Clean up
      deallocate( Mcoll, Mrxn, M, S, Sinv, V0, work, eqDist, X, V, Xinv )
  
  end subroutine
  
  