With[{kmin=99.5,kmax=100.5,dk=4},
Export["/Users/hzhou/Desktop/hz_kf.dat",
Module[{grid,norms,mean},
grid=Ceiling[kmax];
grid={0,Range[-grid,grid,dk]}//Flatten//Union;
grid=Flatten[Outer[List,grid,grid,grid],2];
norms=Norm/@grid//N;
grid=Extract[grid,Position[norms,x_/;kmin<=x<kmax]];
mean=Mean[Norm/@grid]//N;
Flatten@{{Length[grid],mean}//Echo,grid//Transpose}
]
]
]

With[{kdat=Import["/Users/hzhou/Desktop/hz_kf.dat"]},
Module[{n,kave,ki,kx,ky,kz,kvec,kmax},
{{n,kave},ki}={#[[1;;2]],#[[3;;-1]]}&@Flatten[kdat];
{kx,ky,kz}=Partition[#,Length[#]/3]&@ki;
kvec=Transpose[{kx,ky,kz}];
Print["Check n: ",Echo@n==Length[kvec]];
Print["Check kave: ",Echo@kave==Mean[Norm/@N[kvec]]];

Print["All kx: ",Union@kx];
Print["All ky: ",Union@ky];
Print["All kz: ",Union@kz];

kmax=1.5*Max[{kx,ky,kz}//Flatten];
{{ListPointPlot3D[kvec,BoxRatios->{1,1,1},PlotRange->ConstantArray[kmax*{-1,1},3],AxesLabel->{"x","y","z"}],
ListPlot[Cases[kvec,kk_/;kk[[1]]==0]/.{x_,y_,z_}:>{y,z},PlotRange->ConstantArray[kmax*{-1,1},2],AspectRatio->1,FrameLabel->{"y","z"}]},
{ListPlot[Cases[kvec,kk_/;kk[[2]]==0]/.{x_,y_,z_}:>{x,z},AspectRatio->1,PlotRange->ConstantArray[kmax*{-1,1},2],FrameLabel->{"x","z"}],
ListPlot[Cases[kvec,kk_/;kk[[3]]==0]/.{x_,y_,z_}:>{x,y},AspectRatio->1,PlotRange->ConstantArray[kmax*{-1,1},2],FrameLabel->{"x","y"}]}}//Grid
]
]