scf_pot = np.abs(potential.energy(xyz))
scf_pot = scf_pot.value # get numerical value from Astropy Quantity

# log-spaced contour levels
levels = np.logspace(np.log10(scf_pot.min()), np.log10(scf_pot.max()), 16)

plt.figure(figsize=(6,6))

plt.contour(x, z, scf_pot.reshape(x.shape), cmap='inferno_r',
            levels=levels, locator=ticker.LogLocator())

plt.title("Equipotential")
plt.xlabel("$x$", fontsize=22)
plt.ylabel("$z$", fontsize=22)
plt.tight_layout()