anchor.diff:              end do
               do i = 1, nReac
                   index = nRows - nReac + i
  +                ! Reactants
                   dXij(nIsom+i,j) = dXij(nIsom+i,j) + C(index)
               end do
   
  -            ! The isomers and reactant
  +            ! Convert the dXij information for eigenmode j into phenomenological
  +            ! rate coefficients
               do i = 1, nIsom+nReac+nProd
                   K(i,n) = K(i,n) + V(j) * dXij(i,j)
               end do
  
positive.code:     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
  
negative.code:         end do
      end do
  
      ! Generate set of initial condition vectors, one per isomer and reactant
      ! Each initial condition vector corresponds to a case where a Boltzmann
      ! distribution of only one isomer or reactant is present
      ! The Boltzmann distribution occurs as a result of the fast relaxation of
      ! internal energy modes, which is separable from the chemical modes
      allocate( eqDist(1:nRows, 1:nIsom+nReac) )
      do j = 1, nIsom + nReac
          do r = 1, nRows
  