Comparing gradients may be tricky sometimes if you want many digits after the decimal to be the same. First thing -- the SCF convergence threshold should be as tight as possible, since the energy is quadratic in delta C and gradients are linear. In DFT, you really want to run the calculations on exactly same grids. This is often impossible, because different packages use different grids, and it's recommended to use ultra fine grids to minimize numerical noise. I heard that some programs (not TeraChem, of course) don't compute the grid derivatives contribution to the gradient, i.e. their gradients are just incorrect, and you don't want to use them in your tests.