Slightly wrong energy written to TBTGF by sisl?

[AleksBL]

When writing a custom tbtgf file using the write_header of "sisl.io.tbtgfSileTBtrans" the energies in read by TBtrans will differ by x1.00000785 from the input value given inside python.

The calculation of the greens function is going on in the python code below:

There doesnt seem to be any descrepancy when iterating over the f:

The fail happens in TBtrans, in the output:

sisl v. 0.15
tbtrans from siesta-rel-5.0, but also happens with tbtrans from siesta v.4.1.5

The script below shows that it seems to be consistent across python, but something it seems to be some intermediate step

import sisl
import numpy as np
import tqdm

g      = sisl.geom.graphene()
H_elec = sisl.Hamiltonian(g)
H_elec.construct([[0.1,1.5], [0.0, -2.7]])
gamma = sisl.MonkhorstPack(H_elec, [1] * 3)
E    = np.linspace(-1,1)
eta  = 0.001
mu   = None

with sisl.io.tbtgfSileTBtrans('Test.TBTGF') as f:
    val = E + 1j * eta
    # print(val[0])
    # print(val.dtype)
    if mu is not None:f.write_header(gamma, E + 1j * eta, mu = mu)
    else:             f.write_header(gamma, E + 1j * eta)
    
    for ispin, new_k, k, e in tqdm.tqdm(f):
        if new_k:
            f.write_hamiltonian(H_elec.Hk(format='array', dtype=np.complex128),
                                H_elec.Sk(format='array', dtype=np.complex128))
        
        # else:          SeHSE = SRSSE.self_energy(e + 1j*eta, bulk=True, coupling=True)
        rand_mat  = np.random.random((H_elec.no, H_elec.no))
        rand_mat += rand_mat.T
        
        f.write_self_energy(rand_mat)

gf = sisl.get_sile('Test.TBTGF')

Version details
Run the below code and add to bug-report:

import sisl._debug_info as sd
sd.print_debug_info()

[zerothi]

And how is the input file for tbtrans? Because that is where the real problem enters.

I.e. the energies used should be precise enough, if you write with limited precision, then thats the error!

[AleksBL]

These are the first couple of energies in the file containing the custom energies (my_energies from the tbtrans manual i think)

-3.9999 0.01 1.0 eV
-3.933233333333333 0.01 1.0 eV
-3.8665666666666665 0.01 1.0 eV
-3.7998999999999996 0.01 1.0 eV
-3.733233333333333 0.01 1.0 eV
-3.6665666666666663 0.01 1.0 eV
-3.5999 0.01 1.0 eV
-3.533233333333333 0.01 1.0 eV
-3.4665666666666666 0.01 1.0 eV

The script im using for this used to work I think.

[zerothi]

how are you writing it??

[zerothi]

Share the code please.

[AleksBL]

The script:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Sep 12 16:23:57 2024

@author: investigator
"""
import os
os.environ['OMP_NUM_THREADS']='1'
os.environ['SISL_NUM_PROCS']='1'

from siesta_python.siesta_python import SiP
from siesta_python.funcs import zigzag_g
import sisl
import matplotlib.pyplot as plt
import numpy as np
from Zandpack.TimedependentTransport import TD_Transport as TDT

RUN = True
dk  = 1000.0

g = zigzag_g()
tx, ty, td = 2,3,10

gem   = g.tile(tx,0)
gep   = g.tile(tx,0).move(g.cell[0] * (td- tx))
slab  = g.tile(td,0).tile(ty,1)
emmax = gem.xyz[:,0].max()+.1
epmin = gep.xyz[:,0].min()-.1

#y1,y2,y3,y4,y5,y6 = 2.0, 7.0, 6.0, 15.0, 15,32
b11,b12 = 2,  10
b21,b22 = 12, 18
b31,b32 = 25, 31
x1,x2   = 22, 26

def rcond(r):
    r0,r1,r2 = r
    if r0<emmax or r0>epmin:
        return False
    if b11<r1 and b12>r1:
        return False
    if b21<r1 and b22>r1:
        return False
    if b31<r1 and b32>r1:
        return False
    if x1<r0 and x2 > r0 and r1>b11 and r1<b32:
        return False
    return True

dev = slab.remove([i for i in range(slab.na) if rcond(slab.xyz[i])])


C = dev.center()
Dij = np.linalg.norm(dev.xyz - C,axis=1)
idx= np.where(Dij == Dij.min())[0][0]
H  = sisl.geom.sc(1.0, 'H').move(dev.xyz[idx]).move([0,0,3.0])
dev = dev.add(H)


e1  = SiP(gem.cell, gem.xyz, gem.atoms.Z,
         directory_name ='E1', sl = 'E1', sm = 'E1',
         semi_inf       = '-a1',
         elec_RSSE      = True,
         elec_SurRSSE   = True,
         )

e2  = SiP(gep.cell, gep.xyz, gep.atoms.Z,
         directory_name='E2', sl = 'E2', sm = 'E2',
         semi_inf      = '+a1',
         elec_RSSE     = True,
         elec_SurRSSE  = True,
         )


Dev = SiP(dev.cell, dev.xyz, dev.atoms.Z,
          directory_name = 'Device',
          elecs = [e1,e2], Chem_Pot = [0.0, 0.0],
          kp_tbtrans = [1,1,1],
          save_SE = True,
          )
# In[]
vals = [[0.1, 1.5, 3.2], [0.0, -2.7, 0.0]]
Hem  = sisl.Hamiltonian(gem); Hem.construct(vals)
Hep  = sisl.Hamiltonian(gep); Hep.construct(vals)
Hd   = sisl.Hamiltonian(dev); Hd. construct(vals)
dHj  = np.linalg.norm(Hd.xyz[-1] - Hd.xyz, axis=1)
idx  = np.where((dHj>0.1)*(dHj<3.2))[0]
print(idx)
for i in idx:
    Hd[Hd.no-1,i] = 2.0
    Hd[i,Hd.no-1] = 2.0

e1.manual_H(Hem)
e2.manual_H(Hep)
# In[]
R  = TDT ([e1.to_sisl(), e2.to_sisl()],  Dev.to_sisl(), kT_i = [0.025, 0.025], mu_i=[0.25, 0.25])
line = np.linspace(-4,4,121)+1j*1e-2 + 1e-4
line = np.vstack([line]*2)
R.Make_Contour(line, 17, pole_mode = 'JieHu2011')
Dev.custom_tbtrans_contour = R.Contour
if RUN == False:
    assert 1 == 0

e1.Real_space_SI(1, (1,ty,1), 0.0, R.Contour,(1,3,1),
                  parallel_E = False,num_procs = 4,
                  dk = dk)

e2.Real_space_SI(1, (1,ty,1), 0.0, R.Contour,(1,3,1),
                  parallel_E = False,num_procs = 4,
                  dk = dk)

Dev.find_elec_inds()
Dev.fdf(eta = 1e-3)

# assert 1 == 0
Dev.manual_H(Hd.sub(Dev._rearange_indices))

Dev.run_tbtrans_in_dir(DOS_GF = True)
R.Device = Dev
R.read_data()
R.pickle('1Rib_U')
´´

Real_space_SI (More or less the sisl real space SE tutorial):
```python
def Real_space_SI(self, ax_integrate, supercell, eta, Contour,
                      nsc = (1,3,1), 
                      ending = 'TSHS',
                      Hsurf_in = None,
                      only_couplings = False,
                      dk = 1000.0,
                      mu = None,
                      parallel_E = False,
                      num_procs  = 4,
                      sisl_patch_1 = False,
                      keep_unpicklable_objects = False,
                      ):
        """ Similar to Real_space_SE, refer to tutorials for how to use this.
            
        """
        
        import tqdm
        from sisl.physics import RecursiveSI, RealSpaceSI
        if self.elec_RSSE == False:
            print('set elec_RSSE to True!')
            error()
        if self.elec_SurRSSE == False:
            print('set elec_SurRSSE to True!')
            error()
        
        if hasattr(self,'SurRSSE_dict'):
            ax_integrate = self.SurRSSE_dict['ax_integrate']
            supercell    = self.SurRSSE_dict['supercell']
            eta          = self.SurRSSE_dict['eta']
            try:
                Contour  = self.SurRSSE_dict['Contour']
            except:
                pass
            nsc          = self.SurRSSE_dict['nsc']
            Hsurf_in     = self.SurRSSE_dict['Hsurf_in']
        
        E = Contour
        def translate(s):
            return s.replace('a1','A').replace('a2','B').replace('a3','C')
        H_minimal = sisl.get_sile(self.dir + '/' + self.sl + '.'+ending).read_hamiltonian()
        SE        = RecursiveSI(H_minimal, translate(self.semi_inf))
        if Hsurf_in is None:
            Hsurf = H_minimal.copy()
        else:
            Hsurf = Hsurf_in.copy()
        Hsurf.set_nsc(nsc)
        if isinstance(ax_integrate, int):
            if not nsc[ax_integrate]>1:
                print('Wrong nsc or ax_integrate passed!')
                error()
        if isinstance(ax_integrate, list):
            if not (nsc[ax_integrate]>1).any():
                print('Wrong nsc or ax_integrate passed!')
                error()
        
        SRSSE = RealSpaceSI(SE, Hsurf, 
                            ax_integrate, 
                            unfold = supercell,
                            dk = dk)
        if sisl_patch_1:
            print('---Warning---')
            print('Youre now using a custom real_space_coupling method in the RealSpaceSI class')
            
            from siesta_python.sisl_patches import Alt_real_space_coupling
            SRSSE.real_space_coupling = Alt_real_space_coupling.__get__(SRSSE, RealSpaceSI)
            SRSSE.initialize()
        if keep_unpicklable_objects:
            self._RSSI_SE = SE
            self._RSSI_SRSSE = SRSSE
        
        H_elec, elec_indices = SRSSE.real_space_coupling(ret_indices=True)
        # Before we overwrite the electrode Hamiltonian for TBTrans we move the 
        # minimal electrode Hamiltonian
        _ori_name = self.dir + '/' + self.sl + '.'+ending
        _new_name = self.dir + '/' + self.sl + '_minimal_old.'+ending
        os.system('cp '+_ori_name + ' ' + _new_name)
        
        H_elec.write(self.dir + '/'+ self.sl + '.'+ending)
        H = SRSSE.real_space_parent()
        
        #indices = np.arange(len(H))
        #indices = np.delete(indices, elec_indices)
        #indices = np.concatenate([elec_indices, indices])
        np.save(self.dir + '/RS_Coupling_pos', H.xyz[elec_indices])
        np.save(self.dir + '/RS_Coupling_specie', H.atoms.Z[elec_indices])
        if only_couplings:
            return
        
        gamma = sisl.MonkhorstPack(H_elec, [1] * 3)
        sisl.io.tableSile(self.dir + '/contour.E', 'w').write_data(E, np.zeros(E.size) + 0.)
        if parallel_E:
            import joblib as Jl
            global SE_func # for multiprocessing backend
            def SE_func(e):
                return e, SRSSE.self_energy(e, bulk=True, coupling=True)
            results = Jl.Parallel(n_jobs  = num_procs,
                                  backend = joblib_backend,
                                  verbose = 10)(Jl.delayed(SE_func)(e + 1j*eta) for e in E)
            del SE_func
            ResDic = {}
            for r in results:
                ResDic.update({r[0]:r[1]})
        
        with sisl.io.tbtgfSileTBtrans(self.dir +'/' + self.sl +'.TBTGF') as f:
            if mu is not None:f.write_header(gamma, E + 1j * eta, mu = mu)
            else:             f.write_header(gamma, E + 1j * eta)
            
            for ispin, new_k, k, e in tqdm.tqdm(f):
                
                if new_k:
                    f.write_hamiltonian(H_elec.Hk(format='array', dtype=np.complex128),
                                        H_elec.Sk(format='array', dtype=np.complex128))
                if parallel_E: SeHSE = ResDic[e+1j*eta]
                else:          SeHSE = SRSSE.self_energy(e + 1j*eta, bulk=True, coupling=True)
                f.write_self_energy(SeHSE)
´´´

[zerothi]

That isn't surprising:

...
        sisl.io.tableSile(self.dir + '/contour.E', 'w').write_data(E, np.zeros(E.size) + 0.)
...

you are writing with the default precision, which is 1e-5, you have to explicitly set this to something higher! Just do 1e-12 to be sure.

[AleksBL]

I have added fmt='.12e' to the write_data function, and the contour.E file now has higher precision:
-3.999900000000e+00+1.000000000000e-02j 0.000000000000e+00+0.000000000000e+00j
-3.933233333333e+00+1.000000000000e-02j 0.000000000000e+00+0.000000000000e+00j
-3.866566666667e+00+1.000000000000e-02j 0.000000000000e+00+0.000000000000e+00j
-3.799900000000e+00+1.000000000000e-02j 0.000000000000e+00+0.000000000000e+00j
-3.733233333333e+00+1.000000000000e-02j 0.000000000000e+00+0.000000000000e+00j
-3.666566666667e+00+1.000000000000e-02j 0.000000000000e+00+0.000000000000e+00j
-3.599900000000e+00+1.000000000000e-02j 0.000000000000e+00+0.000000000000e+00j
-3.533233333333e+00+1.000000000000e-02j 0.000000000000e+00+0.000000000000e+00j
-3.466566666667e+00+1.000000000000e-02j 0.000000000000e+00+0.000000000000e+00j
-3.399900000000e+00+1.000000000000e-02j 0.000000000000e+00+0.000000000000e+00j
-3.333233333333e+00+1.000000000000e-02j 0.000000000000e+00+0.000000000000e+00j
.....

Electrode E1 has no transfer matrix.
Assuming ../E1/E1.TBTGF contains H, S and E*S - H - \Sigma for calculating the correct self-energies and scattering matrices.
E1 principal cell is perfect!
Electrode Green's function file: '../E1/E1.TBTGF' already exist.
Warning: contours differ by > 9.9999999999999995E-007
ERROR : contours differ by > 9.9999999999999991E-006
TS E.real [eV] TS E.imag [eV] File E.real [eV] File E.imag [eV]
-0.39999000000000002E+01 0.10000000000000002E-01 -0.39998685797720466E+01 0.99999214474663038E-02
-0.39332333333333334E+01 0.10000000000000002E-01 -0.39332024367889376E+01 0.99999214474663038E-02
-0.38665666666666665E+01 0.10000000000000002E-01 -0.38665362938058294E+01 0.99999214474663038E-02
-0.37999000000000001E+01 0.10000000000000002E-01 -0.37998701508227199E+01 0.99999214474663038E-02
-0.37332333333333332E+01 0.10000000000000002E-01 -0.37332040078396123E+01 0.99999214474663038E-02
-0.36665666666666663E+01 0.10000000000000002E-01 -0.36665378648565028E+01 0.99999214474663038E-02
-0.35998999999999999E+01 0.10000000000000002E-01 -0.35998717218733947E+01 0.99999214474663038E-02
-0.35332333333333330E+01 0.10000000000000002E-01 -0.35332055788902856E+01 0.99999214474663038E-02
-0.34665666666666666E+01 0.10000000000000002E-01 -0.34665394359071775E+01 0.99999214474663038E-02
-0.33998999999999997E+01 0.10000000000000002E-01 -0.33998732929240680E+01 0.99999214474663038E-02
-0.33332333333333333E+01 0.10000000000000002E-01 -0.33332071499409599E+01 0.99999214474663038E-02

It seems the constants for converting between Rydberg and eV was updated to this CODATA2018(?) recently, does it have anything to do with this? maybe there is a legacy constant hanging out somewhere in the code that writes the TBTGF files?

[zerothi]

ah, you are right, I thought I had fixed this, try and do tbtGfSileSiesta(..., version="new")
I'll fix this so that it is handled by an env.

[AleksBL]

Ahh right, it works now, thanks!

[zerothi]

Thanks for raising the issue!