def get_uniformity_kpoint(n):
    kPoints = []
    for i in range(n):
        for j in range(n):
            kPoints.append(i/n * rec[0] + j/n * rec[1])
    np.save("kPoints.npy", np.array(kPoints))
    return kPoints

Take points evenly along the inverted parallelogram. Carry out the previous calculation.

The whole Brillouin zone is restored by translation.

import numpy as np
import matplotlib.cm as cm
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
from matplotlib.ticker import LinearLocator, FormatStrFormatter

basis_vector = [[1.37287871,1.37287871,-2.74575742],[-2.74575742,1.37287871,1.37287871],[13.36629497,13.36629497,13.36629497]]
V = np.dot(basis_vector[0], np.cross(basis_vector[1], basis_vector[2]) )
rec = [np.cross(basis_vector[1], basis_vector[2]) * 2 * np.pi/V,
       np.cross(basis_vector[2], basis_vector[0]) * 2 * np.pi/V,
       np.cross(basis_vector[0], basis_vector[1]) * 2 * np.pi/V]

print(rec)
nk = 200
Ek = np.load("Ek_mesh.npy")[:, 0]
Ek = Ek.reshape(nk, nk)
print(Ek.shape)

kPoints = np.load("kPoints.npy")
print(kPoints.shape)
kx = kPoints[:,0]
ky = kPoints[:,1]
print(kx.shape)
kx = kx.reshape(nk, nk)
ky = ky.reshape(nk, nk)



plt.contour(kx, ky, Ek, [0.0])
plt.contour(kx + rec[0][0] -rec[1][0] , ky + rec[0][1] -rec[1][1] , Ek, [0.0])
plt.contour(kx - rec[0][0], ky - rec[0][1] , Ek, [0.0])
plt.contour(kx - rec[1][0], ky - rec[1][1] , Ek, [0.0])
plt.contour(kx  +rec[0][0], ky + rec[0][1] , Ek, [0.0])
plt.show()