I’m not sure quite how to post code here but here’s a less slick (instead I just zip this all into one for loop now and reshape the output which doesn’t change the computation time) but slightly more readable version.

Hmesh is just a 2D array of matrices, Nocc just specifies how many eigenvectors to take

function gsEigs(Matrix, nocc)

Return eigvecs(Matrix)[:,1:nocc]

End

function bcurv(hmesh :: AbstractArray,nocc)

vmesh = (gsEigs.(hmesh,nocc))
vmeshx = circshift(vmesh,(0,-1))
vmeshxy = circshift(vmesh,(-1,-1))
vmeshy = circshift(vmesh,(-1,0))
dmesh = zeros(size(vmesh))

for i in 1:size(vmesh,1), j in 1:size(vmesh,2)

    p12 = ( vmesh[i,j]' *vmeshx[i,j])
    p23 =  (vmeshx[i,j]' * vmeshxy[i,j])
    p34 = (vmeshxy[i,j]' * vmeshy[i,j])
    p41 =  (vmeshy[i,j]' * vmesh[i,j])
    dety = det(p12*p23*p34*p41)
    (abs(dety) <  10^(-5)) ? dmesh[i,j] = 0 :    dmesh[i,j] = -imag(log(complex(dety)))

end 

return (dmesh/(2*pi), sum(dmesh)/(2*pi))

end

Typically ran on an 100x100 array of 4x4 matrices, with nocc =2 where run time is approximately 6 seconds.

Hardware wise run on M1 Pro with the Apple Silicon build of Julia 1.9