[Question/docs] _ncSileTBtrans: pivoting and sorting, device region

[jonaslb]

I had the following happen to me:

In [1]: import sisl as si

In [2]: tbt=si.get_sile("siesta.TBT.SE.nc")

In [3]: l2s=si.utils.list2str

In [4]: l2s(tbt.geom.o2a(tbt.pivot(0), unique=True))
Out[4]: '288-353, 362-373, 382-406, 408-409, 411-415, 418-471'

In [5]: l2s(tbt.geom.sub(tbt.a_dev).o2a(tbt.pivot(0, in_device=True), unique=True))
Out[5]: '3-167'

In [6]: l2s(tbt.a_dev)
Out[6]: '288-490'

In other words, based on Out[4], electrode 0 lies on the orbitals corresponding to, among others, atom 288 which is atom 0 in the device region (as can be seen in Out[6], also this matches my expectation from knowledge of the device geometry). However, when using the in_device parameter, and then using the o2a method of the device geometry, we are told that it lies on atom 3 (from Out[5]). And not only that, but it has become contiguous in the indices even though it shouldn't be.

I did track this down to the pivot function in the _ncSileTBtrans class. I can fix the problem by replacing the call indices(pvt, se_pvt, 0) with indices(self.o_dev, se_pvt, 0), which I believe is what should happen here. However then i found that it would be equivalent to using the sort=True parameter for pivot.

Now the documentation claims that sorting is "mostly useful if the self-energies are returned sorted as well", but I don't see a good reason here to sort the self-energies, unless I've misunderstood the pivoting? In the end I want to do something like H[pvt, pvt.T] += tbt.self_energy(...) where H is a hamiltonian with orbital ordering like you'd get with si.Hamiltonian(device_geometry).Hk(..).

So I wanted to make this a simple PR fixing a bug, but became unsure whether it was intended or not, or if I misunderstood the pivoting. Perhaps that can be cleared up!

[zerothi]

I think it behaves correctly, hopefully.

In [4] you request the orbitals but also pass them to be o2a with unique=True. When you request unique, the atoms are passed through np.unique which also sorts it. So already here it is probably not what you expect.

In [5]. pivot now returns the indices of the first electrode in the pivoted device region. So the orbitals you pass to o2a are not belonging to the geom.sub(tbt.a_dev).

In [6] you get the output of a_dev this is the sorted device region which is not what you have requested in [5].

I hope this clarifies some of the problems, I will amend some documentation by adding notes of what is returned sorted (o_dev and a_dev)

As for your last remark, not exactly what it is supposed to do.

A user may do one of 3 things:

1. Work with the full Hamiltonian of the full thing (i.e. including buffer atoms and downfolding regions). This corresponds to pivot(0, sort=True)
2. Work with the device region Hamiltonian (H.sub(tbt.a_dev)), this corresponds to pivot(0, in_device=True, sort=True)
3. Work with the pivoted device region Hamiltonian (which cannot directly be made in sisl since it may swap orbitals and atoms). However, if you can create this Hamiltonian, one would do pivot(0, in_device=True). Note, I have some scripts that work with the TBtrans BTD matrices (and thus this method is required).

I hope this clarifies.

I have changed the documentation for pivot to this:

in_device : bool, optional
           If ``True`` the pivoting table will be translated to the device region orbitals.
           If `sort` is also true, this would correspond to the orbitals directly translated
           to the geometry ``self.geometry.sub(self.a_dev)``.
        sort : bool, optional
           Whether the returned indices are sorted. Mostly useful if you want to handle
           the device in a non-pivoted manner.

[zerothi]

As a side-note, in_device=False and sort=False is heavily used internally to calculate DOS etc. ;)

[jonaslb]

Thanks for the clarification! The mentioned documentation update sounds good :)

There is also documentation for a sort parameter in the tbtsencSileTBtrans.self_energy method. It says: "if True the returned self-energy will be sorted (equivalent to pivoting the self-energy)". It's not clear to me whether I should use this parameter? I want the sorted, "un-pivoted" situation, but the documentation says that pivoting and sorting are equivalent (perhaps that means 'pivot back to the sorted situation' or some such, but in any case a little unclear).

[zerothi]

I have updated documentation much like in the previous one.

If you see the Examples here: http://zerothi.github.io/sisl/api-generated/sisl.io.tbtrans.tbtsencSileTBtrans.html#sisl.io.tbtrans.tbtsencSileTBtrans

then I think it is more clear what sort does. :) Indeed pivoting and sorting is the same since the pivoting is with respect to the device region, so if you un-pivot it, it is the same as returning the sorted since sorting the pivot table, or sorting the in_device table is the same.

I'll close this for now then :)

Thanks for the question!

[jonaslb]

There is one problem now: The two approaches are not equivalent as documentation says -- as I understand from your previous clarification, using the unsorted approach is relevant when the Hamiltonian is in the pivoted basis, while sorting must be used if the Hamiltonian is in the "regular sorted" basis. See below:

In [1]: import sisl as si, numpy as np                                                                                               

In [2]: tbt = si.get_sile("siesta.TBT.SE.nc")                                                                                        

In [3]: H=np.zeros((tbt.no_d,)*2, dtype=np.complex)                                                                                  

In [4]: H2=H.copy()                                                                                                                  

In [5]: dpvt_unsorted = tbt.pivot(0, in_device=True).reshape(-1, 1)                                                                  

In [6]: dpvt_sorted = tbt.pivot(0, in_device=True, sort=True).reshape(-1, 1)                                                         

In [7]: H[dpvt_unsorted, dpvt_unsorted.T] = tbt.self_energy(0, 0, sort=False)                                                        

In [8]: H2[dpvt_sorted, dpvt_sorted.T] = tbt.self_energy(0, 0, sort=True)                                                            

In [9]: np.allclose(H, H2)                                                                                                           
Out[9]: False

In [10]: devpivot = tbt.pivot(in_device=True).reshape(-1, 1)                                                                         

In [11]: np.allclose(H, H2[devpivot,devpivot.T])                                                                                     
Out[11]: True

[zerothi]

Thanks! You are correct! This is the behaviour intended, I have changed the documentation to reflect what you have shown! ;)

Thanks!!!