Given the input dimension in the function you’re taking a gradient of is really small (3), this is probably better suited for ForwardDiff. Try with
funct(x, y, b) = ForwardDiff.gradient(x -> pair_energy(x, y, b), x)
You’ll also need to remove the ::Float64 annotations on vector1D and change the pair_energy ones to just SVector{3} (or remove them completely). With that I get:
julia> @btime test_analytical_force = analytical_total_force($r, $box_size)
924.744 μs (1 allocation: 2.50 KiB)