mailhexu / TB2J

a python package for computing magnetic interaction parameters
BSD 2-Clause "Simplified" License
67 stars 29 forks source link

sisl and H(k) convention #1

Closed zerothi closed 3 years ago

zerothi commented 4 years ago

Hi He Xu! :) You mentioned something about the convention, now I am intrigued :)

mailhexu commented 4 years ago

Hi Nick,

I calculate the spin density for non-colinear case in this way: $Gk(z)=(zI-Hk)^-1$, where Hk is from sisl. nk=-1/\pi \int_{-\infty}^{Ef} Sk Gk(z) note that I put the Im to later. n is the integration of nk over the BZ. The submatrix ni of each orbital in n(z) is then a 22 matrix. Mi_y = Im (ni[0,1](1j)+ni[1,0](-1j))/2. The M_y calculated from this has the opposite sign to the siesta result. The M_x and M_z has no problem.

In Sisl I find that this function is applied to the Hamiltonian, which seems to be where the convention is defined? Can I find the documentation of the siesta convention somewhere?

Thanks!

def _mat_spin_convert(M, spin=None): """ Conversion of Siesta spin matrices to sisl spin matrices The matrices from Siesta are given in a format adheering to the following concept: A non-colinear calculation has the following entries (in C-index) for the sparse matrix: H[:, [0, 1, 2, 3]] H11 == H[:, 0] H22 == H[:, 1] H12 == H[:, 2] - 1j H[:, 3] # spin-box Hermitian H21 == H[:, 2] + 1j H[:, 3] Although it really does not make sense to change anything, we do change it to adhere to the spin-orbit case (see below). I.e. what Siesta saves is the -Im[H12], which we now store as Im[H12]. A spin-orbit calculation has the following entries (in C-index) for the sparse matrix: H[:, [0, 1, 2, 3, 4, 5, 6, 7]] H11 == H[:, 0] + H[:, 4] H22 == H[:, 1] + H[:, 5] H12 == H[:, 2] - 1j H[:, 3] # spin-box Hermitian H21 == H[:, 6] + 1j H[:, 7] """ if spin is None: if M.spin.is_noncolinear: M._csr._D[:, 3] = -M._csr._D[:, 3] elif M.spin.is_spinorbit: M._csr._D[:, 3] = -M._csr._D[:, 3] elif spin.is_noncolinear: M._D[:, 3] = -M._D[:, 3] elif spin.is_spinorbit: M._D[:, 3] = -M._D[:, 3]

Best wishes, HeXu

zerothi commented 3 years ago

Hi HeXu,

Sorry for the massive delay. I have been struggling my head around this for months, and I think I have finally found a stable and workable solution.

Would you be able to check for me, or share a script for me? That would be great.

I have pushed a branch here: https://github.com/zerothi/sisl/tree/phases-change which seems to fix everything (hopefully)!

mailhexu commented 3 years ago

Thanks a lot! I will test and send the results. I had a workaround for it (by reversing the sign of the y component of the Hamiltonian) and I'll remove this to see if things are still the same. My guess was that there is a sign difference in the convention but from your code it seems more complex than what I have thought. Below are the list of things changed from what I understand:

  1. The y component of spin y component is reversed in _help.py

    • H12 == H[:, 2] - 1j H[:, 3] # spin-box Hermitian H12 == H[:, 2] + 1j H[:, 3] # spin-box Hermitian
  2. in sparse.py: conjugate the spin. I guess this does the same thing.

  3. The H,S, DM, and EDM are transposed. Can you tell me why?

  4. in _phase.pyx and block.py, The phase has been conjugated for Hk=\sum_R H(R) * phase. Is this related to the spin y component or it's just another unrelated change? In TB2J, I define this convention in a variable as each TB/DFT code could use different convention. Can you tell me why this should changed? I should modify TB2J to use different convention for different version of sisl.

Best wishes, HeXu

zerothi commented 3 years ago

Thanks a lot! I will test and send the results. I had a workaround for it (by reversing the sign of the y component of the Hamiltonian) and I'll remove this to see if things are still the same.

Thanks!!

My guess was that there is a sign difference in the convention but from your code it seems more complex than what I have thought. Below are the list of things changed from what I understand:

  1. The y component of spin y component is reversed in _help.py
  • H12 == H[:, 2] - 1j H[:, 3] # spin-box Hermitian H12 == H[:, 2] + 1j H[:, 3] # spin-box Hermitian

This was actually a typo, so nothing really changed here.

  1. in sparse.py: conjugate the spin. I guess this does the same thing. change is sparse.py is unrelated, mainly for testing stuff ;)

  2. The H,S, DM, and EDM are transposed. Can you tell me why? This is actually the culprit. Problem was that siesta stores data in CSC format, while I rely on CSR format. This means that I need to transpose+conjugate (Hermitian transpose) data to correct format.

  3. in _phase.pyx and block.py, The phase has been conjugated for Hk=\sum_R H(R) * phase. Is this related to the spin y component or it's just another unrelated change? In TB2J, I define this convention in a variable as each TB/DFT code could use different convention. Can you tell me why this should changed? I should modify TB2J to use different convention for different version of sisl.

This is also because of the above CSR <-> CSC changes.

Basically the phases were probably wrong, but due to sisl never really used spin-orbit+non-collinear I never figured this out, I always had k = -k and thus never saw this difference. I think this is why it pops up now. :)

I think you can safely assume all versions prior to 0.10.0 (current release) has the bug. And once you can assert I am doing things correctly, ;) I can merge the branch and get on with it. :)

Best wishes, HeXu

Thanks for bringing this to my attention HeXu!

zerothi commented 3 years ago

Ok, so I have been thinking a bit more about the convention. And I'll have to work a little more before you test?

So please hold on if you haven't already started.

mailhexu commented 3 years ago

Thanks. I will then wait for your new implementation.

zerothi commented 3 years ago

Thanks. I will then wait for your new implementation.

Thanks, and I have just pushed the fixes. Now I am more sure! ;)

mailhexu commented 3 years ago

Great! I will try it tomorrow.

zerothi commented 3 years ago

So you should check now.

I actually think that your code should stop if you are not using the next sisl version, just to be sure :)

mailhexu commented 3 years ago

Indeed my workaround could be incomplete if the Hamiltonian is transposed. B

mailhexu commented 3 years ago

BTW, do you know if there is a way to stop if people are using old version of TB2J with the new sisl? The results will be wrong in that case.

zerothi commented 3 years ago

Indeed my workaround could be incomplete if the Hamiltonian is transposed. B

Yes, as far as I can tell, the old was using the transposed matrix, plus the transposed spin-box plus a sign change in the phase. So all in all this translates to a sign change in the y-component, as far as I can tell.

But this should hopefully fix it.

BTW, do you know if there is a way to stop if people are using old version of TB2J with the new sisl? The results will be wrong in that case.

No, you can't change people's installations without them doing anything. You could probably make a new release soon that crashes if using old sisl. I will hopefully release pretty soon!

mailhexu commented 3 years ago

Hi, I tested and the spin densities are now good. I use only the Hamiltonian and the overlap matrix so other things are not covered. There is a small issue when reading from the HSX file, there is an error : File "/home/ICN2/xhe/.local/lib/python3.7/site-packages/sisl/io/siesta/binaries.py", line 594, in _xij2system atoms_all = np.stack(atoms) File "<__array_function__ internals>", line 6, in stack File "/home/ICN2/xhe/.local/lib/python3.7/site-packages/numpy/core/shape_base.py", line 425, in stack raise ValueError('all input arrays must have the same shape') seems it is trying to stack the index of orbitals for different species. With the netcdf file it is fine. Thanks a lot! I will make a new version after your release.

Cheers, HeXu

zerothi commented 3 years ago

Thanks! Great, really good!

Yeah, the HSX files are really annoying since they don't really contain anything that lets me decide the geometry. So basically to get correct behaviour one should either pass a geometry hsx.read_hamiltonian(geometry=...) or use the netcdf file. Please note that the HSX file only uses single precision which may be problematic for large BZ calculations. In that case netcdf and double precision is much better!

I have just updated code with more appropriate error messages, could you please have a look and see if it makes sense?
Once this is completed, then I'll merge!

zerothi commented 3 years ago

I am just going to merge in, it is working, and now fine-tuning needs feedback from other 3rd parties. Will be in master asap!

I'll let you close this issue.

mailhexu commented 3 years ago

Thanks!