anchor.diff:      end do
   
       ! Calculate the phenomenological rate constants
  -    ! This version is intended to follow the notation of Miller and Klippenstein
  +    ! This version follows the notation of Miller and Klippenstein
       allocate( dXij(1:nIsom+nReac+nProd, 1:nCSE), C(1:nRows) )
  +    ! Iterate over isomers and reactants to determine k(T,P) with each as the
  +    ! reactant
       do n = 1, nIsom+nReac
           do j = 1, nCSE
   
  +            ! dXij contains the change in isomer/reactant/product i as a result 
  +            ! of the jth eigenmode, using isomer/reactant n as the starting 
  +            ! point
               do i = 1, nIsom+nReac+nProd
                   dXij(i,j) = 0.0
               end do
  
positive.code:             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
  
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
  