Hi,
Unfortunately the autograd does not compute Jacobians, but vector Jacobian products.
Note that here I did not give the output of NN to the .grad() but the loss itself. There are two main reasons here:

I would argue that the gradient you want to enforce is not on the Jacobian of the function but the gradient of the loss wrt the input.
The loss is scalar, so a single backward pass actually gives the full Jacobian (gradients) and you can do the penalty directly without having to do a for-loop of .grad()