As @ptrblck  mentioned, you might have to put the network into eval mode to fix BN.
There is a small thing in your observation code which actually modifies the gradients. .numpy() doesn’t make a copy but stores a reference, and then you += more to it. Also, if you accumulate the gradients across batches (by not calling zero_grad), you don’t have to accumulate them manually or you count the earlier ones multiple times. If you comment out the accumulation code and just compare the p.grad of parameters afterwards, it seems to be the same.