anchor.diff:  
   	! Construct accounting matrix
   	row = 1
  -	do r = minval(nRes)+1, nGrains
  +	do r = 1, nGrains
   		do i = 1, nIsom
   			if (r > nRes(i)) then
   				indices(r,i) = row
  
positive.code: 	Z = 0 * Z
  	
  	! Collisional terms in active-state matrix and RHS vectors
  	do i = 1, nIsom
  		do s = nRes(i)+1, nGrains
  			do r = max(nRes(i)+1, s - halfbandwidth), min(nGrains, s + halfbandwidth)
  
negative.code: 			Z(indices(s,i), i) = sum(Mcoll(i,s,1:nRes(i)) * densStates(i,1:nRes(i)) * exp(-E(1:nRes(i)) / 8.314472 / T))
  		end do
  	end do
  	halfbandwidth = halfbandwidth * nIsom
  
  	! Isomerization terms in active-state matrix and RHS vectors
  	do i = 1, nIsom
  		do j = 1, i-1
  