Support for SIESTA-BdG in sisl

[ahkole]

Describe the feature

As was mentioned in a discussion between @zerothi and @rreho #869, it would be good to have support for the BdG Hamiltonian in sisl that is used/produced by SIESTA-BdG (which will be merged into main SIESTA sometime in the future) (see https://journals.aps.org/prb/abstract/10.1103/PhysRevB.110.134505 for implementation paper). Here we will supply some details about how the BdG Hamiltonian is implemented in SIESTA so that we can discuss what is necessary to add this support to sisl.

First of all, we will explain the structure of the Hamiltonian. The structure of the BdG Hamiltonian is (see paper for more details):

   !        | H^ee           Delta |   | H^ee          Delta     |
   ! HBdG = |                      | = |                         |
   !        | Delta^dagger   H^hh  |   | Delta^dagger  -(H^ee)^* |

Within SIESTA-BdG there are two matrices, one for H^ee and one for Delta, both with 8 spin components (same as SOC). The spin components have the following meaning in the SIESTA-BdG code:

   ! 1.  The electron-electron part uses the same structure as in the SOC case
   !     with an additional shift of the chemical potential to zero:
   !
   !                     | H_{jo,iuo}^{u,u}  H_{jo,iuo}^{u,d} |
   !     H^ee_{jo,iuo} = |                                    |
   !                     | H_{jo,iuo}^{d,u}  H_{jo,iuo}^{d,d} |
   !
   !                                            | S_{jo,iuo}  0          |
   !                                 - mu       |                        |
   !                                            | 0           S_{jo,iuo} |
   !
   !         | H(ind,1) + i H(ind,5) - mu S(ind),  H(ind,3) - i H(ind,4)             |
   !       = |                                                                       |
   !         | H(ind,7) + i H(ind,8)            ,  H(ind,2) + i H(ind,6) - mu S(ind) |
   !
   !
   !
   ! 2.  The hole-hole part is -conjg(H^ee) and hence very similar
   !
   !     H^hh_{jo,iuo} = -(H^ee)^*
   !
   !         | - H(ind,1) + i H(ind,5) + mu S(ind),              - H(ind,3) - i H(ind,4) |
   !       = |                                                                           |
   !         | - H(ind,7) + i H(ind,8)            ,  - H(ind,2) + i H(ind,6) + mu S(ind) |
   !
   !
   ! 3.  The structure of the electron-hole part is different from the
   !     usual convention.
   !        - indices 1,2 correspond to the singlet contribution
   !        - indices 3,..,8 to the three triplet contributions  !
   !
   !
   ! The full electron-hole part is the sum of these four contributions
   !
   !     Delta = d_S + d_T(u,u) + d_T(d,d) + d_T(0)
   !
   ! yielding
   !
   !     Delta_{jo,iuo}
   !
   !        |  HBdG(ind,3) + i HBdG(ind,4)  ,   HBdG(ind,1) + i HBdG(ind,2)   |
   !        |                                   + HBdG(ind,7) + i HBdG(ind,8) |
   !      = |                                                                 |
   !        | -HBdG(ind,1) - i HBdG(ind,2)      HBdG(ind,5) + i HBdG(ind,6)   |
   !        |  + HBdG(ind,7) + i HBdG(ind,8),                                 |
   !
   !
   ! The lower off-diagonal block (hole-electron) is related to the
   ! upper-diagonal block by symmetry
   !
   !     (Delta^dagger)_{jo,js,iuo,is} = -conjg(Delta_{jo,js,iuo,is})
   !
   !        | -HBdG(ind,3) + i HBdG(ind,4)  ,  HBdG(ind,1) - i HBdG(ind,2)   |
   !        |                                  + HBdG(ind,7) - i HBdG(ind,8) |
   !      = |                                                                |
   !        | -HBdG(ind,1) + i HBdG(ind,2)    -HBdG(ind,5) + i HBdG(ind,6)   |
   !        |  + HBdG(ind,7) - i HBdG(ind,8),                                |
   !
   !
   !  Singlet:
   !
   !  (d_S)_{jo,iuo}
   !
   !       | 0                           , HBdG(ind,1) + i HBdG(ind,2) |
   !     = |                                                           |
   !       | -HBdG(ind,1) - i HBdG(ind,2), 0                           |
   !
   ! Triplet:
   !
   !     (d_T(u,u))_{jo,iuo}
   !
   !          | +HBdG(ind,3) + i HBdG(ind,4), 0                           |
   !        = |                                                           |
   !          | 0                           , 0                           |
   !
   !     (d_T(d,d))_{jo,iuo}
   !
   !          | 0                           , 0                           |
   !        = |                                                           |
   !          | 0                           , HBdG(ind,5) + i HBdG(ind,6) |
   !
   !
   !     (d_T(0))_{jo,iuo}
   !
   !          | 0                           , HBdG(ind,7) + i HBdG(ind,8) |
   !        = |                                                           |
   !          | +HBdG(ind,7) + i HBdG(ind,8), 0                           |

To avoid confusion, it's important to mention that inside the code, we use HBdG to refer to the matrix object that stores the Delta part of the BdG Hamiltonian. At the end of a calculation, we save H^ee and Delta separately to two HSX files (which both have 8 spin components), namely SystemLabel.HSX and SystemLabel.HBdG.HSX respectively.

We hope that this description is clear enough and a good start for a discussion. If you have any questions let us know.

On another note, @nils-wittemeier do you think we can give @zerothi access to the SIESTA-BdG repo? It might be useful if he can see the code.

[zerothi]

Thanks for this, let me note that I think the entire details should just be stored in one HSX file.
Or, the DELTA matrix should be stored in a different file format, DX or similar (no need to store the S matrix twice).

Is there a reason for not storing everything in one HSX file?

[ahkole]

Is there a reason for not storing everything in one HSX file?

The reason was mostly convenience (or laziness 😅). By doing it this way we could re-use the write_hsx routine from SIESTA directly. And the Delta part could also be read immediately by sisl without having to modify anything. But you are right that maybe it would be better to write everything in a single HSX file. Especially since that avoids writing S twice. I think that should probably also not be too difficult to change, right?

[zerothi]

Is there a reason for not storing everything in one HSX file?

The reason was mostly convenience (or laziness 😅). By doing it this way we could re-use the write_hsx routine from SIESTA directly. And the Delta part could also be read immediately by sisl without having to modify anything. But you are right that maybe it would be better to write everything in a single HSX file. Especially since that avoids writing S twice. I think that should probably also not be too difficult to change, right?

I think the siesta write still works for nspin>8, no? Otherwise it really should! :)

[ahkole]

I think the siesta write still works for nspin>8, no? Otherwise it really should! :)

Hmm, we actually never tried. Maybe we can look into this. That does require us to supply a single object that contains all 16 spin components to the routine (right now the two parts are stored in different objects both with 8 spin components).

[zerothi]

I have created a merge-request so it can be simpler to write a HSX file object https://gitlab.com/siesta-project/siesta/-/merge_requests/400

This should enable you to store two different hamiltonian objects in a single HSX file.

This should be very simple for you! :)

[ahkole]

That's great, thanks! We will update the writing of the HSX file in the SIESTA-BdG code.

[ahkole]

How would you btw handle the writing of Delta and H^ee in other files (such as SystemLabel.nc or DMHS.nc)?

[zerothi]

How would you btw handle the writing of Delta and H^ee in other files (such as SystemLabel.nc or DMHS.nc)?

I would also extend it with more spin dimensions. Unless there are good reasons for not to do it.

[zerothi]

That's great, thanks! We will update the writing of the HSX file in the SIESTA-BdG code.

If you like the solution, I would heavily suggest you to comment in the merge-request, so it can get merged in, then you can rebase ontop of dev!

[ahkole]

If you like the solution, I would heavily suggest you to comment in the merge-request, so it can get merged in, then you can rebase ontop of dev!

I commented on the request to indicate that this would be useful for us.

[zerothi]

Is there a reason why the ordering of the Delta matrix elements are not aligned with that of the Hamiltonian spin-box? It seems a bit odd to have a different ordering of elements, although the matrix spin box looks very different, I think the chances of people mis-understanding this is high.

readability here is likely more important IMHO!

I.e. lets assume one wishes to create a tight-binding model of this thing, would they work in the singlet+triplet basis, or just use the Delta spin-box indices for assigning values?

[ahkole]

Is there a reason why the ordering of the Delta matrix elements are not aligned with that of the Hamiltonian spin-box?

Maybe @nils-wittemeier can also comment on this but as far as I understood the reason is that we store the singlet and three triplet components as four complex numbers sequentially in HBdG, i.e. Singlet = HBdG(1) + i HBdG(2), Triplet_{u,u} = HBdG(3) + i HBdG(4), etc., see also the last part of the comment above that explains the matrix form of the singlet and triplet contributions. This is inherited from Gabor's initial implementation I believe and this has probably been done this way for easy mapping between input and HBdG. It's always possible of course to internally combine the singlet and triplet components to store the Delta spin-box in HBdG in the same order as for the electron-electron Hamiltonian spin-box. I don't know which one would have better readability?

I.e. lets assume one wishes to create a tight-binding model of this thing, would they work in the singlet+triplet basis, or just use the Delta spin-box indices for assigning values?

For users it would be easier to work in the singlet-triplet basis I think. This is also how they specify the input Delta now in SIESTA-BdG. I.e., they would specify that they want singlet pairing between atom 1 and 2 or something. Otherwise you would have to manually figure out which components of the Delta spin-box to set to get the correct singlet or triplet pairing that you want.

[zerothi]

Is there a reason why the ordering of the Delta matrix elements are not aligned with that of the Hamiltonian spin-box?

Maybe @nils-wittemeier can also comment on this but as far as I understood the reason is that we store the singlet and three triplet components as four complex numbers sequentially in HBdG, i.e. Singlet = HBdG(1) + i HBdG(2), Triplet_{u,u} = HBdG(3) + i HBdG(4), etc., see also the last part of the comment above that explains the matrix form of the singlet and triplet contributions. This is inherited from Gabor's initial implementation I believe and this has probably been done this way for easy mapping between input and HBdG. It's always possible of course to internally combine the singlet and triplet components to store the Delta spin-box in HBdG in the same order as for the electron-electron Hamiltonian spin-box. I don't know which one would have better readability?

NOTE: this should be discussed in a future MR against siesta!

IMHO, it would make sense to have the internal representation to be just as the siesta hamiltonian, right now the Hamiltonian indices are difficult enough to grasp, and with this additional complexity it becames very hard.
Unless, the singlet/triplet components are used separately I would argue it would make sense to have them "understandable" ;)

I.e. lets assume one wishes to create a tight-binding model of this thing, would they work in the singlet+triplet basis, or just use the Delta spin-box indices for assigning values?

For users it would be easier to work in the singlet-triplet basis I think. This is also how they specify the input Delta now in SIESTA-BdG. I.e., they would specify that they want singlet pairing between atom 1 and 2 or something. Otherwise you would have to manually figure out which components of the Delta spin-box to set to get the correct singlet or triplet pairing that you want.

Yes, but in sisl? Would you still do it like this?

[nils-wittemeier]

Hmmm ... I have to say that the seperation in singlet and triplet component seems to be much easier to understand for users compares to the abstract uu ud etc.

How things should be handled internally is ofcourse another matter and I understand your point as well, Nick. For maintenance of the code it might be better to stick to the normal convention. As you pointed out this should be discussed on the siesta MR. I do quickly want to noted that this would require a decomposition into singlet triplet for grid operations.

[zerothi]

Ok, seems like the single/triplet separation makes sense internally as well.

[rreho]

I would like to comment that singlet and triplet components are easier to grasp/commonly used in the superconductors community and they are associated to specific physical phenomena/interpretations.

[ahkole]

I do quickly want to noted that this would require a decomposition into singlet triplet for grid operations.

@nils-wittemeier this is no different from what has to be done for the diagonal (electron-electron) quantities right? There you also have to combine/decompose the DM spin-box to create the X and Y spin components for rho in rhoofd.F90 for the SOC case.

So I don't think that would prevent us from using the same convention for the spin-box of Delta internally.

[nils-wittemeier]

That's a fair point

[zerothi]

How does the density matrix look like in BdG? Does the DM have the same elements, singlet + triplet?

[zerothi]

And could you comment on how PDOS etc looks like with wavefunctions from BdG simulations?

[zerothi]

Could you please have a look in #877? It should implement everything for nambu code style in sisl (except pdos and some minor details).

[ahkole]

How does the density matrix look like in BdG? Does the DM have the same elements, singlet + triplet?

Correct me if I'm wrong @nils-wittemeier: The structure is similar as for the Hamiltonian. So we have two separate matrix objects in the SIESTA-BdG code, DM and ADM (anomalous density matrix), where the ADM is analogous to Delta (HBdG in the code), so we also store the singlet and triplet components in the same order in the ADM matrix, i.e. Singlet = ADM(1) + i ADM(2), Triplet_{u,u} = ADM(3) + i ADM(4), etc. We also store the ADM in a separate file at the moment (with extension .ADM). I guess for consistency we should then also instead store everything in a single DM file with 16 spin components? I guess that also requires a MR for SIESTA that allows writing of two DM objects in a single DM file?

And could you comment on how PDOS etc looks like with wavefunctions from BdG simulations?

The wavefunctions that come out of diagonalizing the BdG Hamiltonian are 4-component Nambu spinors of the form,

           | e_{u} |
psi =      | e_{d} |
           | h_{u} |
           | h_{d} |

where e are the electron coefficients and h are the hole coefficients. In the orbital basis the output of cdiag would have the shape

psi = | e^{s1,orb1,u} e^{s2,orb1,u} .... e^{sn,orb1,u} |
      | e^{s1,orb1,d} e^{s2,orb1,d} .... e^{sn,orb1,d} |
      | h^{s1,orb1,u} h^{s2,orb1,u} .... h^{sn,orb1,u} |
      | h^{s1,orb1,d} h^{s2,orb1,d} .... h^{sn,orb1,d} |
      | e^{s1,orb2,u} e^{s2,orb2,u} .... e^{sn,orb2,u} |
      | e^{s1,orb2,d} e^{s2,orb2,d} .... e^{sn,orb2,d} |
      | h^{s1,orb2,u} h^{s2,orb2,u} .... h^{sn,orb2,u} |
      | h^{s1,orb2,d} h^{s2,orb2,d} .... h^{sn,orb2,d} |
      |     .                .                 .       |
      |     .                .                 .       |
      |     .                .                 .       |
      | h^{s1,orbn,d} h^{s2,orbn,d} .... h^{sn,orbn,d} |

Here s1 (etc.) labels the states (bands) and orb1 (etc.) labels the orbitals.

For the PDOS the number of spin components is doubled and we basically store first the electron PDOS and then the hole PDOS with for both the same spin order as for SOC, so PDOS = e^{uu} e^{dd} e^{x} e^{y} h^{uu} h^{dd} h^{x} h^{y}. We compute the different spin components in the PDOS in the same way as in the SOC case, just using the electron coefficients for the electron PDOS and the hole coefficients for the hole PDOS. In the normal PDOS the electron and hole coefficients don't mix.

I hope the above is clear enough.

[ahkole]

Could you please have a look in #877? It should implement everything for nambu code style in sisl (except pdos and some minor details).

That is fast! On first glance things look good. I'll try to have a closer look later (also comparing with how we do things in the SIESTA-BdG code).

[zerothi]

PDOS is now implemented, and one should also be able to do all other kinds of stuff, I think the only thing missing is the trs operation of the Matrix....

[rreho]

I have a general question on this thread. Would it also be a good idea to make sure that we can exploit sisl.viz functionalities to plot the quantities that we compute within the SIESTA-BdG code?

Would those changes done for the PDOS ensure compatibility with sisl.viz ?
As an additional note, we also compute and plotgrid quantities such as the anomalous density in real space, the initial and final pairing potential $\Delta$.
However, we did implement those changes in sisl but used our own python functions.

[zerothi]

I have a general question on this thread. Would it also be a good idea to make sure that we can exploit sisl.viz functionalities to plot the quantities that we compute within the SIESTA-BdG code?

Would those changes done for the PDOS ensure compatibility with sisl.viz ? As an additional note, we also compute and plotgrid quantities such as the anomalous density in real space, the initial and final pairing potential $\Delta$. However, we did implement those changes in sisl but used our own python functions.

Of course these things would be ideal to have there as well.

I don't know what is involved in the grid-operations you have. But I can't see that there should be a huge limitation here, so of course that should exist in sisl. I have taken a first stab, and then I would advice you to test that backend, once in, you can extend sisl to your liking and give a PR.