Energy differentiation with cut-off function

I am a new learner of molecular dynamics (MD) simulations methods and has a simple question regarding handling of cutoff functions. In MD, pairwise energy between two atoms is assumed to be a function, $$f(r)$$ (eg Lennard Jones potential) which is multiplied by a cut off function, $$t(r)$$ to make interactions beyond cutoff distance to be zero. So technically pairwise energy of a system is written as

$$E = \sum_{i,j} f(r_{ij})*t(r_{ij})$$

How do we handle force term then? Is it the differentiation of the whole term as written above or is it the differentiation of $$\sum_{i.j}f(r_{ij})$$ and then multiplied separately by the same cutoff function.

It is necessary to differentiate the whole thing. Monitoring the energy conservation is an essential aspect of MD simulation algorithms, and this will not even be possible in principle if the forces are not derived from the actual potential energy. So you need both terms: $$f'(r)t(r) + f(r)t'(r)$$ where $$'$$ represents the derivative (and then you multiply by the appropriate vector to give the gradient and hence the force).
Even then, there will be some numerical issues arising from any lack of smoothness of the modified energy function at the cutoff distance. Usually $$t(r)$$ is chosen to vanish at the cutoff, or alternatively an additive term is also included (a vertical shift) so that the pairwise energy has no step discontinuity at the cutoff distance: such a discontinuity would cause a delta-function in the forces, which would be tricky to handle. Even if this is done, though, the existence of a discontinuity in the higher derivatives of the potential will have a more subtle effect on the accuracy of the algorithm and the energy conservation; in practice, the cutoff distance is chosen to be large enough to make these effects small. Also, quite often $$t(r)$$ is chosen to make the first derivative of the product $$t(r)f(r)$$ vanish at the cutoff distance.
There is a bit of formal mathematics in the literature regarding these small discontinuities in higher derivatives. Symplectic MD algorithms are known to exactly conserve a modified Hamiltonian which is within a small amount, $$\delta t^n$$, of the true Hamiltonian, where $$\delta t$$ is the time step and the exponent $$n$$ depends on the algorithm. But this nice result relies on being able to expand the potential energy in a series in $$r$$, i.e. it breaks down if the potential is not smooth.