anchor.diff:          do i = 1, nIsom
               K(nIsom+n,nIsom+n) = K(nIsom+n,nIsom+n) - sum(Fim(i,n,:))
           end do
  +        ! Reaction from isomer or reactant j to reactant n
           do j = 1, nIsom+nReac
               do i = 1, nIsom
                   K(n+nIsom,j) = K(n+nIsom,j) + sum(Gnj(n,i,nRes(i)+1:nGrains) * pa(nRes(i)+1:nGrains,j,i))
  
positive.code:     do i = 1, nIsom
          bres(i) = sum(eqDist(i,1:nRes(i)))
      end do
      
      ! Determine pseudo-steady state populations of active state
      if (nIsom == 1) then
          ! Not worth it to call banded solver when only one well
  
negative.code:         end do
      end do
  
      ! Renormalize
      do i = 1, nIsom
          pa(:,i,:) = pa(:,i,:) / bres(i)
          eqDist(i,:) = eqDist(i,:) / bres(i)
      end do
  
      ! Initialize phenomenological rate coefficient matrix
      do i = 1, nIsom+nReac+nProd
          do n = 1, nIsom+nReac+nProd
  