Your code contains a small bug, in the leapfrog HMC function the line pnext = -p should be pnext = -pnext. Also this line is only of theoretical importance, it

Your code contains a small bug, in the leapfrog HMC function the line pnext = -p should be pnext = -pnext. Also this line is only of theoretical importance, it has no effect on the final result (unless you use a weird asymmetric kinetic energy).