anchor.diff:  
       ! Initialize msg to empty string
       msg = ''
  -
  -    ! Zero collision matrix
  -    do r = 1, Ngrains
  -        do s = 1, Ngrains
  -            Mcoll(r,s) = 0.0
  -        end do
  -    end do
       
       ! Allocate/zero vector of normalization coefficients
       allocate( C(1:Ngrains) )
  
positive.code:     R = 8.314472
  
      ! Evaluate collisional transfer probability - theoretically correct way
      !if (Ej < E0 .or. Ei < E0 .or. rho(j) .eq. 0) then
      !   rate = 0.
      !elseif (Ej >= Ei) then
      !   rate = exp(-(Ej - Ei) / alpha)
      !else
      !   rate = exp(-(Ei - Ej) / alpha) * rho(i) / rho(j) * exp( -(Ei - Ej) / (R * T))
      !end if
  
      ! Evaluate collisional transfer probability - move density of states ratio to other side
      ! so that it is always less than 1
      if (Ej < E0 .or. Ei < E0 .or. rho(j) .eq. 0) then
          rate = 0.
      elseif (Ej >= Ei) then
          rate = exp(-(Ej - Ei) / alpha) * rho(i) / rho(j)
      else
          rate = exp(-(Ei - Ej) / alpha) * exp( -(Ei - Ej) / (R * T))
      end if
  
  end subroutine
  
negative.code:     end do
  
      ! Normalize using detailed balance
      do r = start, Ngrains
          C(r) = (1 - sum(C(start:r-1) * Mcoll(start:r-1,r))) / sum(Mcoll(r:Ngrains,r))
          ! Check for normalization consistency (i.e. all numbers are positive)
          if (C(r) <= 0) then
              msg = 'Error normalizing collisional transfer probabilities matrix.'
              return
          end if
  
      end do
      do r = start, Ngrains
          Mcoll(r,r+1:Ngrains) = Mcoll(r,r+1:Ngrains) * C(r)
          Mcoll(r+1:Ngrains,r) = Mcoll(r+1:Ngrains,r) * C(r)
          Mcoll(r,r) = Mcoll(r,r) * C(r) - 1
      end do
  
  