anchor.diff:      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
  
positive.code:     ! Determine phenomenological rate coefficients from pa(1:nGrains, 1:nIsom+nReac, 1:nIsom)
      ! Rows relating to isomers
      do i = 1, nIsom
          do r = 1, nRes(i)
              K(i,i) = K(i,i) + sum(Mcoll(i,r,1:nRes(i)) * eqDist(i,1:nRes(i)))
          end do
          do j = 1, nIsom
              do r = 1, nRes(i)
                  K(i,j) = K(i,j) + sum(Mcoll(i,r,nRes(i)+1:nGrains) * pa(nRes(i)+1:nGrains,j,i))
              end do
          end do
          do n = nIsom+1, nIsom+nReac
              do r = 1, nRes(i)
                  K(i,n) = K(i,n) + sum(Mcoll(i,r,nRes(i)+1:nGrains) * pa(nRes(i)+1:nGrains,n,i))
  
negative.code:     real(8), dimension(1:nGrains, 1:nIsom+nReac, 1:nIsom)  ::  pa
  
      ! Indices
      integer i, j, m, n, r
  
      real(8), dimension(1:nIsom,1:nGrains) :: eqDist
      real(8), dimension(1:nIsom) :: bres
  