1. generate many k-points in the first Brillouin zone using a very dense k-mesh;
2. calculate the eigenvalues for those k-points generated in the 1st step;
3. arrange the calculated data into a proper format, e.g, bxsf for Fermi surface; and then visualize it by XCRYSDEN .
or in the following format
kx_1 ky_1 eigenvalue_1 (highest occupied \pi band)
.....
kx_n ky_n eigenvalue_n (highest occupied \pi band)
kx_1 ky_1 eigenvalue_1 (lowest unoccupied \pi band)
....
kx_n ky_n eigenvalue_n (lowest unoccupied \pi band)
The data in the above format can by plotted by the GNUPLOT, python matplotlib, or MATLAB.