makeh[g_]:=Module[{Eval,Evec,U},
{Eval,Evec}=Eigensystem[N[g]];
U=Evec\[Transpose];
Return[U.DiagonalMatrix[Arg[Eval]].U\[ConjugateTranspose]];
];
fjac[g0_,g1_,g2_,g3_]:=Module[{g01,g02,g03,h01,h02,h03},
g01=g1.g0\[ConjugateTranspose];
g02=g2.g0\[ConjugateTranspose];
g03=g3.g0\[ConjugateTranspose];
{h01,h02,h03}=makeh/@{g01,g02,g03}//Chop;
Return[12*I*Tr[h01.(h02.h03-h03.h02)]/(24*Pi^2)/(24)];
];
σ[n_]:=PauliMatrix[n];
m=2.;
lis={};
Do[
q=Sin[kx] σ[1]+Sin[ky] σ[2]+Sin[kz] σ[3]+I (Cos[kx]+Cos[ky]+Cos[kz]+m) σ[0];
q=Table[{u,λ,v}=SingularValueDecomposition[(q/.{kx->nx*2*Pi/L,ky->ny*2*Pi/L,kz->nz*2*Pi/L})];Chop[u.v\[ConjugateTranspose]],{nx,0,L},{ny,0,L},{nz,0,L}];
cube=Table[
{fjac[q[[nx,ny,nz]],q[[nx+1,ny,nz]],q[[nx,ny+1,nz]],q[[nx,ny,nz+1]]](*0124*),
-fjac[q[[nx+1,ny,nz]],q[[nx,ny+1,nz]],q[[nx+1,ny+1,nz]],q[[nx+1,ny,nz+1]]](*1235*),
fjac[q[[nx+1,ny,nz]],q[[nx,ny+1,nz]],q[[nx,ny,nz+1]],q[[nx+1,ny,nz+1]]](*1245*),
-fjac[q[[nx,ny+1,nz]],q[[nx+1,ny+1,nz]],q[[nx+1,ny,nz+1]],q[[nx,ny+1,nz+1]]](*2356*),
fjac[q[[nx,ny+1,nz]],q[[nx,ny,nz+1]],q[[nx+1,ny,nz+1]],q[[nx,ny+1,nz+1]]](*2456*),
-fjac[q[[nx+1,ny+1,nz]],q[[nx+1,ny,nz+1]],q[[nx,ny+1,nz+1]],q[[nx+1,ny+1,nz+1]]](*3567*)}//Chop,{nx,1,L},{ny,1,L},{nz,1,L}];
AppendTo[lis,Total[cube//Flatten]];
Print[{L,lis[[-1]]}];
Print[ListPlot[lis]];
,{L,1,35}];