Anomalous Hall tensor conductivity

[apezol]

Hi, I was playing with sisl and noticed that it provides a way to calculate the AHC tensor. However when calculating this quantity for a non-orthogonal basis it gives a warning message. I think the non-orthogonality of the basis is taken into account in the definition of the velocity matrix operator but would It be possible to make a transformation in the Hamiltonian itself instead?, something like: H' = S^(-1/2) H S^(-1/2) where S is the overlap matrix. Since the calculation of the AHC involves only velocity matrix elements, these may come from the now 'orthogonalized' Hamiltonian.

[zerothi]

Hi this is the same problem as outlined in #356.

So I will close this issue, please track the other one.

[zerothi]

Erhm, actually the velocity matrix can easily be calculated (and sisl do this) even for non-orthogonal basis sets.

As far as I can see then sisl can calculate the conductivity without problems? Did you try?

I.e. Was it a warning or error?

[zerothi]

The warning is there because I haven't tested this against known values and hence put it in. If you can assert it works, then I can remove the warning?

[apezol]

Hi, I haven't had time to construct or simulate a particular model that displays a large value to compare the results. But I have tested a simulation I had for MoS2(a TR invariant system) for which the AHC is vanishing. So, what I did is to calculate the AHC tensor with sisl using:

cond = si.physics.electron.conductivity(bz)

which gives: [ 2225.66563072, 1969.53093094, -0. ] [-1923.70036413, -822.65225741, -0. ] [ -0. , -0. , -0. ]

Then I calculated the same quantity by using a simple routine that diagonalizes a shifted Hamiltonian to emulate the derivative. When omitting the Overlap matrix, the ahc is ~ 2630, while if I use the Lowdin transformed Hamiltonian the AHC is ~ -0.02, so it's a zero. In all cases the conductivity was calculated at E=0, being a gapped material with TRS, I think the only value calculated that makes sense is the last one, i.e., that one corresponding to the Lowdin transformation.

I think that physics.electron.conductivity() has some other arguments, in particular one which takes into account degeneracies, but since for this particular system the zero energy corresponds to the gap, we don't need to care about it right?.

[zerothi]

Could you pass the input and python script for this?

[apezol]

sisl_test.zip I'm sharing the archives from siesta and a notebook having the lines. Maybe it's a little time consuming apart from the fact the these calculations usually need an enough fine grid, hopefully this helps (the numbers within the loop are particular for this system since it has 41 states). Overcoming this issue might lead to extending the AHC function to the SHC for instance (or orbital Hall for that matter).

[zerothi]

Could you try out the latest commits? I think there were some problems with the decoupling of states. Since the velocity decoupling is important to get correct velocities for matrix multiplications.
So I think this was a problem of degenerate states, and now your code also seems to do exactly this (i.e. degenerate states are not handled). Whether this is the true cause is not immediately apparent to me, but also I don't think that the units in your approach and in mine are the same. I use the explicit derivatives which has / epsilon explicitly made in correct units.

cond = array([[-0.00095338,  0.00085425, -0.        ],
       [-0.00087331,  0.00033775, -0.        ],
       [-0.        , -0.        , -0.        ]])

ahc(0.) = -2.0240379398032258

ahc_prime(0.) = -2.135904152515741

Again, I don't think that the Lowdin is necessary since the velocity matrix can be calculated exactly in the non-orthogonal basis.

Could you test and report back whether it looks good now?

[apezol]

Hi Nick, I've obtained the same value you posted and it makes more sense now. You're right, there might be some factor I'm missing since I've only applied that script for dimensionless toy models. Now, looking at the function that calculates the velocity matrix is clearer how the Overlap matrix enters into the calculation, a certain reminiscence on how the Green's function is obtained. On the other hand, If I need to calculate the conductivity for an arbitrary energy, what would be the argument?.

[zerothi]

Ok, great. :)
I must admit I am no expert in this calculation. And in fact I think that the sisl unit is wrong, so I would really like some insight and some tests that indicate to me an exact number + checks on units etc. ;)

As for which energy, I think you should just pass the occupation function to conductivity if so desired. The default distribution function is the Fermi-Dirac distribution, to change it to something different, do:

# kT = 0.025 eV
# Ef = 0
dist = si.get_distribution("fermi-dirac", 0.025, 0.)
# or kT = 0.05 eV; Ef = -0.1 eV
dist = si.get_distribution("fermi-dirac", 0.05, -0.1)
cond = conductivity(bz, distribution=dist)

or similar.

[apezol]

haldane.zip Hi Nick, thanks for the suggestion. Attached is a simple test for the Haldane model. Within the gap the value of the Hall conductance is -9.072277938746234e-05 and apart from the units it seems to reproduce(qualitatively) the same results I get using other codes. Since this model is well known to display the value of conductance quantum within the gap(regardless its value), maybe it can help to define the units. I think is customary to express the conductivity in terms of S/cm(\Omega^-1 cm^-1) , then e^2/h is roughly ~ 3874 \Omega^-1 cm^-1.

[zerothi]

Ok, so I have looked a bit more at this. With the following script and the latest commit I get.

#!/usr/bin/env python3
import tqdm
import numpy as np
import sisl as si

uc = si.SuperCell([[1, 0, 0], [0.5, 3 ** .5 / 2, 0], [0, 0, 1]], nsc=[3, 3, 1])
gr = si.Geometry(np.dot([[1/3, 1/3, 0], [2/3, 2/3, 0]], uc.cell), si.Atom(6, R=0.58), uc)

m = 0.2
t1 = 1.
t2 = m/3.0/np.sqrt(3)*2.0
phi = np.pi / 2

H = si.Hamiltonian(gr, dtype=np.complex128)
H.construct([[0.1, .6], [0, t1]])

H[0, 0] = -m
H[1, 1] = m

H[0, 0, (1, 0)] = t2 * np.exp(1j * phi)
H[0, 0, (-1, 0)] = np.conjugate(t2 * np.exp(1j * phi))

H[1, 1, (1,-1)] = t2 * np.exp(1j * phi)
H[1, 1, (-1,1)] = np.conjugate(t2 * np.exp(1j * phi))

H[1, 1, (0, 1)] = t2 * np.exp(1j * phi)
H[1, 1, (0, -1)] = np.conjugate(t2 * np.exp(1j * phi))

H[1, 1, (1, 0)] = np.conjugate(t2 * np.exp(1j * phi))
H[1, 1, (-1, 0)] = t2 * np.exp(1j * phi)

H[0, 0, (1, -1)] = np.conjugate(t2 * np.exp(1j * phi))
H[0, 0, (-1, 1)] = t2 * np.exp(1j * phi)

H[0, 0, (0, 1)] = np.conjugate(t2 * np.exp(1j * phi))
H[0, 0, (0, -1)] = t2 * np.exp(1j * phi)

E = np.linspace(-3, 3, 3)
E = [0.]

def conductance(e):
    dist = si.get_distribution("step_function", e)
    cond = si.physics.electron.conductivity(bz, distribution=dist)
    return cond[0][1]

nk = [81, 81, 1]
#nk = [3, 3, 1]

bz = si.MonkhorstPack(H, nk, trs=False)
values = np.asarray(list(map(conductance, tqdm.tqdm(E))))
print(values[0])

bz.set_parent(H.tile(2, 0))
values = np.asarray(list(map(conductance, tqdm.tqdm(E))))
print(values[0])

The values I get are (I think the units are now S/cm^-1):

-4251.441264431725

for both because I am now dividing by the BZ volume. Would you agree this is ok close to the expected value? Is there some other way to test this?

[zerothi]

I will leave it as it is and add a warning until someone tells me it is correct or suggests a change. :)

[apezol]

Hi Nick, sorry for the missed reply. Indeed, the units seem correct for this one. On the other, I've been testing a routine for extending the AHC to other quantities related, like the spin and orbital Hall conductivities. As a benchmark I'm attaching a picture comparing the results using sisl (left) to those in a publication (shown in the picture), for the Spin Hall conductivity in Pt. That is, changing the velocity operator for the spin/orbital current operators will lead to the desired transport quantity. I guess this isn't implemented in sisl yet, but it would be useful.

[zerothi]

Great, thanks, could you show the code, or the equations that change? Thanks!

[apezol]

Yes!, here it is. I made comments for some lines in the python file, it isn't optimized (of course) but for this case, when time reversal symmetry is preserved, the time consumption might be reduced by choosing the Brillouin zone with this symmetry. The main difference with respect to the AHC, is related to the following line: '' ... np.dot(np.dot(np.conjugate(v[:,i]), 0.5((dHdy-w[j]dSdy)@Sz_total+Sz_total@(dHdy-w[j]*dSdy))), v[:,j]) ... ''

where one of the current operators (the velocity matrix) has been changed to the spin current operator, which is the anti commutator of the spin and velocity ( {v,S} ). In general the S operator can be the spin/orbital matrices depending on the type of response. I'm making some tests for the orbital part and I can share them later.

In short, for the implementation in sisl, I guess the function for the velocity can take an arbitrary operator (O), which turns out to be the identity in the simple case of the AHC.

sisl_test.zip

[zerothi]

So what you are saying is that all that is needed is to change the operator in the velocity matrix calculation?

Could you provide a paper that describes the operators and terminology behind the anti-commutator of v and S?

[apezol]

Yes, here is for the spin part:

https://arxiv.org/pdf/cond-mat/0502351.pdf

and for the orbital is given in this one:

https://arxiv.org/pdf/1804.02118.pdf

Although, maybe I need to check this one: https://arxiv.org/pdf/1810.07637.pdf I've only tested the SHC for Pt and seems to work fine, but the formula applies for the periodic part of the Bloch states only. I'll check Eq. 23 in the last paper since I wasn't able to obtain the right result for the orbital part.