System with non-orthogonal basis / Z2Pack on top of Siesta

[nils-wittemeier]

Hey,

I was wondering whether there is a physical reason, why the tight-binding and especially hamilton matrix interfaces only consider orthogonal basis sets. If not do, you believe, that the changes to hm.py (attached) that I have made would be reasonable?

My motivation: I would like to use Z2Pack ontop of SIESTA. SIESTA used localise pseudo-atomic orbitals as the basis set, i.e. the hamiltonian and overlap matrix can be retrieved in a tight-binding-like form and for example read with the sisl package, which also provides interfaces for overlap and Hamiltonian at any k.

I have also been working on extending the functionality of SIESTA to produce wannier90-type .mmn not only for calculations with collinear to non-collinear spins.

Hence, if the above worked it would be a good tool to test.

Thank you for your time.

hm_with_overlap.patch.txt

[greschd]

I don't see a reason why this shouldn't work - seems like a reasonable extension for the code. Please let me know if the test with SIESTA worked as expected.

If we can get a good example and / or test for this functionality, it would be good to merge it into the repository. Please let me know if you are up for creating the corresponding pull request.

[bfocassio]

To test if the sisl Hamiltonians could be used along with z2pack, I tried the Haldane model and Kane-Mele model, written with sisl language, and using Z2pack to compute the Chern number and Z2 invariant, respectively. Everything seems to be working properly. Both examples are attached.   I performed the changes proposed in the issue, and tested it with flat bismuthene computed with SIESTA. The necessary files to run the Z2 calculation are attached.

The problem is that I'm finding it very hard to converge the calculations and get reliable results. The calculation is very sensitive to the 'pos_tol' that needs to be increased at least to 0.05, and it gets worse with increasing the number of lines. Any help is appreciated.

siesta_sisl_z2pack.zip

[greschd]

This problem might be unrelated to the SIESTA implementation - it looks suspiciously like there was a gap closure around ky ~= 0.3. When the direct band gap vanishes anywhere on the surface to be considered, the Z2 invariant becomes ill-defined.

Could you check if this might be the case?

[bfocassio]

The gap does not seem to close anywhere in the BZ.

[greschd]

[Note that I didn't try running the code yet, only looked at the output in the Notebook]

From this figure it seems there are 4 occupied bands, whereas in the notebook len(occupied) is 12 - maybe this is the source of these problems?

[bfocassio]

[I updated the zip because I was using the development version of sisl. Now the code should run with the latest release version, 0.9.8.]

In the image i chose to show only some bands. There is a total of 12 occupied bands.

[greschd]

Ah, I see - so the other occupied bands are below the ones shown here?

I'll have to give this a closer look to see what might be going wrong. Probably this is another place where orthonormality is assumed, and the overlap would need to be taken into account.

[nils-wittemeier]

Indeed, the overlap would have to be taken into account there as well.

I don't see a way of calculating the overlap between two wavefunctions at different k using only the overlap matrix and Hamiltonian matrix. I believe one would also require the matrix elements of e^{ikr}.

Unfortunately, I did not find time recently to finalize the implementation of the Wannier90 modules for spin-orbit coupling in SIESTA either.

[bfocassio]

I think this could be solved by applying the Löwdin transformation, at each k point. With this we get a Hamiltonian in an orthogonal basis.

At each k we have H and S,

We compute the eigenvalues s and eigenvectors U of the overlap matrix, as SU = sU

Then, we compute s^(-1/2) for each s in s, write it as a diagonal matrix s^(-1/2)

Now, transform back to the S basis with S^(-1/2) = U s^(-1/2) U^(dagger)

And finally compute H' = S^(-1/2) H S^(-1/2)

I performed some tests with 2 different systems and several lines. It seems to be stable and working fine. See the attached file for using this with bismuthene.

Let me know if this procedure is correct

siesta_z2pack_bismuthene.zip

[greschd]

Sorry for the late reply - yes, I think this is correct. I'll note how I think about this below, so you can confirm we have the same understanding. It's basically just the derivation of the Loewdin transformation:

Assume B is the matrix containing the non-orthogonal basis vectors as column vectors, with respect to some (arbitrary) orthonormal basis.

S = B^(\dagger) B

We can do singular value decomposition on B:

B = V \Sigma W^(\dagger)

S = W \Sigma^2 W^(\dagger)

S^(-1/2) = W \Sigma^(-1) W^(\dagger)

If H_0 is the Hamiltonian w.r.t. the arbitrary orthonormal basis above,

H = B^(\dagger) H_0 B

H' = S^(-1/2) H S^(-1/2) = = W \Sigma^(-1) W^(\dagger) B^(\dagger) H_0 B W \Sigma^(-1) W^(\dagger) = = W V^(\dagger) H_0 V W^(\dagger)

which is exactly what we want, the Hamiltonian w.r.t a new orthonormal basis.

In principle, this can be used without modifying the Z2Pack code by simply providing H' instead of H and S, correct? Of course, it might pay to include it in the code instead of everyone having to figure it out by themselves. But maybe this means the best way to implement this is not by modifying the hr code itself, but by creating a child class that uses the hr code (similar to what's done for tight-binding models).

[bfocassio]

Yes, that's correct.

Currently I'm calculating H' then passing it to Z2Pack. It would be nice have this implemented within Z2Pack. I can look into this.

[greschd]

Currently I'm calculating H' then passing it to Z2Pack. It would be nice have this implemented within Z2Pack. I can look into this.

Would be much appreciated! Thinking about it some more, I now think adding it directly to z2pack.hm.System might actually be the best idea. Otherwise the class hierarchy becomes a bit complicated, and it would be cumbersome to e.g. add it to z2pack.tb.

So it'd probably amount to adding a basis_overlaps keyword (default None) to that class, and add a property _hamilton_orthogonal that either directly calls the _hamilton or performs the orthogonalization.

Regarding the procedure itself: When the s entries approach zero this becomes numerically unstable - but I guess at that point the basis is also bad for other parts of the calculation, right?

[greschd]

This was merged in #94 (thanks again @bfocassio and @juijan) and released in the 2.2.0 version.

It would be good to add some tests, but I added a separate issue for that (#95).