Slightly wrong energy written to TBTGF by sisl? #827
Comments
AleksBL
changed the title
|
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! |
|
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 The script im using for this used to work I think. |
|
how are you writing it?? |
|
Share the code please. |
|
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)
´´´ |
|
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 |
|
I have added fmt='.12e' to the write_data function, and the contour.E file now has higher precision: Electrode E1 has no transfer matrix. 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? |
|
ah, you are right, I thought I had fixed this, try and do |
|
Ahh right, it works now, thanks! |
|
Thanks for raising the issue! |
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
Version details
Run the below code and add to bug-report:
The text was updated successfully, but these errors were encountered: