Calculate force on atoms given the energy of configuration, trajectory of simulation and position of atoms in each frame
tl,dr: Forces are equivalent to the negative of the gradient, use small, arbitrary $\mathrm{d}r$ for each atom

The problem, as you mention, with many neural network methods in chemistry at the moment, is that the $E = f(x,y,z)$ is a black box and thus analytical gradients (and thus forces) are often not available.

First, as you clearly know, $F = -\nabla E(x,y,z)$ so you need to solve $\frac{\mathrm{d}E}{\mathrm{d}r}$ for your neural network. (I indicate this for anyone reading who doesn't realize that the force will inherently be the negative of the gradient of the energy function.

So what you need is a finite difference method. This is, of course, a common problem (e.g., on Math Stack Exchange or Wikipedia):

$$ \frac{\partial E}{\partial x} \approx \frac{f(x+\epsilon)-f(x)}{\epsilon}$$

Picking a good step $\epsilon$ is a bit tricky, but I'd probably pick something like 0.01 Å or 0.001 Å, as a start. In your case, it sounds like you want Cartesian gradients, so you'll take your $f(x,y,z)$ and perform the partial derivatives for all three axes (i.e., the initial point $(x,y,z)$ plus 3 extra energy evaluations per atom.

That is, you'd loop through each atom, move it a bit in the x, y, and z directions and evaluate the system to get the gradient/force on that atom. You'll, therefore, have $3N$ energy evaluations to get all the forces on all the atoms.

You can perform better (and more computationally expensive) finite difference methods if you choose.