Clarification on Recent Changes to sisl.physics.Hamiltonian

[rreho]

I am encountering some difficulties understanding the recent changes made to the sisl.physics.Hamiltonian object in the newer versions of sisl. My script, which was built based on the conventions in sisl 0.13.0, no longer behaves as expected (I cannot diagonalize the new Hamiltonian object or the bands give the wrong results).

Specifically, I would like to know if there have been any modifications to how .HSX files are read or how orbitals are specified within the Hamiltonian object. I am working on constructing an HBdG Hamiltonian, which is twice the size of a typical SOC Hamiltonian, and I suspect there might have been changes to either the spin or orbital indexing.

I have attached the relevant portion of my script below for your reference. Any insights or general information about the updates to these objects would be greatly appreciated.

Thank you very much for your help!

def create_atom_h(geom_0,ion_file):
    '''
    Create a new geometry from another geometry replacing the tag to h
    '''
    geom=geom_0.copy()
    atom_ion = sisl.io.siesta.ionxmlSileSiesta(f'{ion_file}')
    atom_ion_basis = atom_ion.read_basis()
    for iatom,atom in enumerate(geom.atoms):
        geom.atoms[iatom] = sisl.Atoms([dict(Z=atom.Z, tag=f'{atom.tag}_h',orbitals=atom_ion_basis.orbitals)])
    geom.reduce()
    return geom
 
def apply_minus_conj(H):
    # Real part of H is stored in [0, 1, 2, 6]
    H[:, :, [0, 1, 2, 6]] = -H[:, :, [0, 1, 2, 6]]
def construct_H_BdG(Hee, Delta_eff, chem_pot, geometry=None):
    """Construct the BdG Hamiltonian from the electron-electron and Delta parts."""
   
    # First shift Hee to Fermi level
    Hee = deepcopy(Hee)
    Hee.shift(-chem_pot)
   
    # Now correct the order of the spin indices for Delta_eff
   
    # In SIESTA the 2x2 spin-box for Delta_eff is constructed like this:
    # !  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,7) + i HBdG(ind,8)             HBdG(ind,5) + i HBdG(ind,6) |
 
    # Meanwhile, sisl constructs the 2x2 spin-box like this:
    # H11 == H[:, 0] + 1j H[:, 4]
    # H22 == H[:, 1] + 1j H[:, 5]
    # H12 == H[:, 2] + 1j H[:, 3] # spin-box Hermitian
    # H21 == H[:, 6] + 1j H[:, 7]
   
    # So we have to map all the spin-components to right one for sisl to construct the Delta spin-box correctly
    # Also note that while reading sisl does H[:, 3] = -H[:, 3]
    Delta_eff = deepcopy(Delta_eff)
    D = Delta_eff._csr._D
    tmp = deepcopy(D)
    tmp[:,0] = D[:, 2]
    tmp[:,1] = D[:, 4]
    tmp[:,2] = D[:, 0] + D[:, 6]
    tmp[:,3] = D[:, 1] + D[:, 7]
    tmp[:,4] = -D[:, 3]  # NB: minus sign because sisl flips sign when reading HSX file
    tmp[:,5] = D[:, 5]
    tmp[:,6] = -D[:, 0] + D[:, 6]
    tmp[:,7] = -D[:, 1] + D[:, 7]
    tmp[:,8] = 0   # Set overlap matrix for Delta part to zero
    Delta_eff._csr._D[:,:] = tmp[:,:]
   
    # Get dense represenations of Hamiltonian blocks
    dense_Hee = Hee._csr.todense()
    dense_Hhh = deepcopy(dense_Hee)
    apply_minus_conj(dense_Hhh)
    dense_Delta_eff = Delta_eff._csr.todense()
    dense_Delta_eff_star = deepcopy(dense_Delta_eff)
    apply_minus_conj(dense_Delta_eff_star)
   
    # Construct HBdG matrix with these blocks
   
    # The structure of the new BdG Hamiltonian will be:
    # | Hee(1, 1, 1) ... Hee(1, no_u, 1) Delta_eff(1, 1, 1) ... Delta_eff(1, no_u, 1) Hee(1, 1, 2) ... Hee(1, no_u, 2) Delta_eff(1, 1, 2) ... Delta_eff(1, no_u, 2) ... ... |
    # |    ...                                                                                                                                                              |
    # | Hee(no_u, 1, 1) ... ...                                                                                                                                             |
    # | Delta_eff_star(1, 1, 1) ... Delta_eff_star(1, no_u, 1) Hhh(1, 1, 1) ... Hhh(1, no_u, 1) Delta_eff_star(1, 1, 2) ... ... ...                                         |
    # |   ...     ...                    ...                                                                                                                                |
    # | Delta_eff_star(no_u, 1, 1) ... ...                                                                                                                                  |
   
    # Here H(io, jo, isc) where io and jo denote orbitals in the unit cell and isc denotes the index of a supercell
   
    no_u, no_s, dim = dense_Hee.shape
    dense_BdG = np.empty((2*no_u, 2*no_s, dim), dtype=dense_Hee.dtype)
    dense_BdG.reshape(2*no_u, -1, 2*no_u, dim)[:no_u, :, :no_u, :] = dense_Hee.reshape(no_u, -1, no_u, dim)[:,:,:,:]
    dense_BdG.reshape(2*no_u, -1, 2*no_u, dim)[:no_u, :, no_u:, :] = dense_Delta_eff.reshape(no_u, -1, no_u, dim)[:,:,:,:]
    dense_BdG.reshape(2*no_u, -1, 2*no_u, dim)[no_u:, :, :no_u, :] = dense_Delta_eff_star.reshape(no_u, -1, no_u, dim)[:,:,:,:]
    dense_BdG.reshape(2*no_u, -1, 2*no_u, dim)[no_u:, :, no_u:, :] = dense_Hhh.reshape(no_u, -1, no_u, dim)[:,:,:,:]
   
    # Convert to sparse matrices needed for building sisl Hamiltonian object
    sp_BdG_matrices = []
    for i in range(dim-1):
        sp_BdG_matrices.append(csr_matrix(dense_BdG[:,:,i]))
    sp_BdG_overlap = csr_matrix(dense_BdG[:,:,-1])
   
    # Construct the new geometry object
    # NB: You can also specify your own geometry to be used (i.e. if you want to plot wavefunctions but the geometry inside the Hamiltonian objects does not contain information about radial part of orbitals
    if geometry is not None:
        soc_geom = geometry
        soc_geom_h = create_atom_h(geom_0=soc_geom,ion_file=f'/Users/Reho0001/workQE/Projects/test-bdg/Pb/Pb.ion.xml')
    else:
        soc_geom = Hee.geometry
    BdG_geom = soc_geom.add(soc_geom_h)  # Effectively doubles the geometry to create both electron and hole atoms/orbitals
   
    # Finally construct new Hamiltonian object
    H_BdG = sisl.Hamiltonian.fromsp(BdG_geom, sp_BdG_matrices, sp_BdG_overlap)
   
    return H_BdG

[pfebrer]

To clarify, what you do is to read an HSX file from SIESTA with sisl and then convert it into sisl?

Does that mean that the Hamiltonian was simply not correctly read? I mean, sisl should handle converting from SIESTA to sisl conventions, no?

Perhaps the change is that this has been fixed? Just speculating here, I have no facts to back this up😅

[zerothi]

Hmm...

I don't think any of the things you mention have changed.

You are changing the sign of :, 3] which is correct (to be Siesta coherent).

One thing, pre 5, I would not trust HSX files... I would at least check whether the TSHS file format is more complete.
Since 5, the two files are very similar, so that should be ok.

What is not clear from the way you do this, is that your BdG hamiltonian should have twice as many orbitals, but I don't see you incorporating that? I.e. for every orbital you'll have 2 2x2 spin-box matrices, no?

Once #861 is in, I think it would make sense to make a BdG hamiltonian a first-class citizen, that would be better, no?

[zerothi]

I guess what I want, is a more minimal example where these things can be more clarified, i.e. a matrix plot of a 16x16 matrix, or something like that?

[zerothi]

To clarify, what you do is to read an HSX file from SIESTA with sisl and then convert it into sisl?

Does that mean that the Hamiltonian was simply not correctly read? I mean, sisl should handle converting from SIESTA to sisl conventions, no?

Perhaps the change is that this has been fixed? Just speculating here, I have no facts to back this up😅

Yeah, how is sisl reading a Hamiltonian with 16 indices? Does it work?

[rreho]

I will be back with a minimal example. In short, the idea is to create a new Hamiltonian object with 8 spin indices but with double the number of orbitals. H and Delta are two independent .HSX files which are merged to create an object like this
| H Delta|
| -Delta * - H*|

[zerothi]

I will be back with a minimal example. In short, the idea is to create a new Hamiltonian object with 8 spin indices but with double the number of orbitals. H and Delta are two independent .HSX files which are merged to create an object like this | H Delta| | -Delta * - H*|

Indeed that is how I also understood it.

1. Could you clarify how the BdG method stores the Hamiltonian elements. I.e. what does the 16 elements correspond to in the HSX file?
2. Any minimal example would be great, yes! Thanks!!
3. What is your BdG code based on? Is it rebased with dev? If I am not mistaken then the original BdG code originated before 5.0 release, which means it likely used the older HSX file format, which could be erroneous in some situations. So there could be other things in play as well..

[rreho]

I found the commit that gives the correct plot. I share a folder with a notebook, files and plot of correct and wrong superconducting bands.
Pb.zip

I exclude it is related with siesta version before or after 5.0 as I tested for both cases.

I share here a minimal example that works with this commit of sisl

commit be29de5f72fde34f6791aec8227ec45aabcbf125
Author: Nick Papior <[email protected]>
Date:   Wed May 22 12:10:33 2024 +0200

    fixed hsxSileSiesta reading from fdf siles

and it does not with this specifications

[sys]
3.11.9 | packaged by conda-forge | (main, Apr 19 2024, 18:45:13) [Clang 16.0.6 ]
[sisl]
version                       0.1.dev4906+g0181f9f
path                          /Users/Reho0001/opt/anaconda3/envs/aiida-2024/lib/python3.11/site-packages/sisl
CC                            /Library/Developer/CommandLineTools/usr/bin/clang
CFLAGS                        -O3 -DNDEBUG
C version                     12.0.0.12000032
FC                            /usr/local/bin/gfortran
FFLAGS                        -O3
FC version                    11.3.0
cython build version          Cython version 3.0.11
numpy build version           2.1.2
[sisl.definitions]
NPY_NO_DEPRECATED_API         NPY_1_20_API_VERSION
CYTHON_NO_PYINIT_EXPORT       1
[sisl.cmake_args]
CMAKE_BUILD_TYPE              Release
WITH_FORTRAN                  ON
F2PY_REPORT_ON_ARRAY_COPY     10
WITH_F2PY_REPORT_COPY         OFF
WITH_F2PY_REPORT_EXIT         OFF
WITH_COVERAGE                 OFF
WITH_LINE_DIRECTIVES          OFF
WITH_ANNOTATE                 OFF
WITH_GDB                      OFF
NO_COMPILATION                
[sisl.env]
SISL_LOG_FILE                 /Users/Reho0001/workQE/Projects/test-bdg/Pb
SISL_LOG_LEVEL                INFO
SISL_NUM_PROCS                1
SISL_PAR_CHUNKSIZE            0.1
SISL_TMP                      /Users/Reho0001/workQE/Projects/test-bdg/Pb/.sisl_tmp
SISL_CONFIGDIR                /Users/Reho0001/.config/sisl
SISL_FILES_TESTS              /Users/Reho0001/workQE/Projects/test-bdg/Pb/_THIS_DIRECTORY_DOES_NOT_EXIST_
SISL_SHOW_PROGRESS            False
SISL_IO_DEFAULT               
SISL_UNIT_SIESTA              codata2018
[runtime]
numpy                         1.26.4
scipy                         1.14.0
xarray                        2024.7.0
netCDF4                       1.7.1.post2
pandas                        2.2.2
matplotlib                    3.9.2
dill                          0.3.8
pathos                        0.3.2
skimage                       0.24.0
plotly                        5.23.0
ase                           3.23.0
pymatgen                      not found
                              module 'pymatgen' has no attribute '__version__'
[install]
pip                           numpy==1.26.4 scipy==1.1
4.0 xarray==2024.7.0 netCDF4==1.7.1.post2 pandas==2.2.2 matplotlib==3.9.2 dill==0.3.8 pathos==0.3.2 scikit-image==0.24.0 plotly==5.23.0 ase==3.23.0
conda                         numpy==1.26.4 scipy==1.14.0 xarray==2024.7.0 netCDF4==1.7.1.post2 pandas==2.2.2 matplotlib==3.9.2 dill==0.3.8 pathos==0.3.2 scikit-image==0.24.0 plotly==5.23.0 ase==3.23.0

The BdG code is a code that we developed in SIESTA and we are merging it right now with dev

[zerothi]

I am having difficulties understanding everything. It could be as simple as reading the geometry, because I don't think anything else has changed since that commit.

So could you try and get a better understanding of where this happens, i.e. by comparing different values of what you are reading? That would be very helpful to get better at pin-pointing this.

You could even do: https://git-scm.com/docs/git-bisect
to get the exact commit that introduces this (no manual search, just re-installing and running).
I really don't know what I am looking for here ;)

[pfebrer]

I found the commit that gives the correct plot.

You mean that with the next commit it doesn't work anymore?

[zerothi]

I have amended the test suite of sisl, see 443c51f which adds checks for eigenvalue spectrums of sisl and siesta.

It basically reads the .EIG and .KP files and diagonalizes the Hamiltonian at the same k-points and compares the eigenvalues. For spin-orbit coupling it succeeds. So it reads the HSX file in a consistent manner. As far as I can see.

But I would be interested in knowing how your problem occurs, because it seems that this would be good to have in sisl as well.

[zerothi]

And could you share the BdG branch? Then we can get that in as 1st class citizen!