Closed AleksBL closed 9 months ago
Yes, the matrices that are to be written for TBTrans inputs is that it requires the diagonal elements to be defined. So I think the error makes sense? What did you expect?
How to define them without requiring to make a new Hamiltonian with the fromsp method? even if i do
for i in range(Hd.no):
Hd[i,i] = 0.5
it still wont write
Could you create a complete example. Your example isn't fully complete.
g.xyz with .txt extension instead: g.txt
import numpy as np
import sisl
IP_cutoff=1.6
OP_cutoff=6.
def tbbi(geom,os_0=0.,os_1=None,Vpppi=-2.7,Vpps=0.48,d0_00=1.42,d0_01=3.35,q0_00=2.,q0_01=None,dangling=0.2,field=0.,eps=0.44):
"""Tight-binding Hamiltonian for graphene nanostructure bilayers. The model uses an
orthogonal basis set and the matrix elements between different orbitals are
determined by Slater-Koster-type two-center $\pi$ and $\sigma$ hopping integrals
between $p_z$ orbitals:
Parameters
----------
geom: str, sisl.Geometry
Atomic structure of the system.
os_0: float, optional
On-site energy of the atoms in the bottom layer.
os_1: float, optional
On-site energy of the atoms in the top layer. I will be equal to os_0 by default.
field: float, optional
Electric field perpendicular to the layers. os_1 will be determined according to
this and the interlayer distance.
eps: float, optional
Screening coefficient accounting for the electric field-induced charge redistribution.
Vpppi: float, optional
$pp\pi$ hopping integral.
Vpps: float, optional
$pp\sigma$ hopping integral.
d0_00: float, optional
Interatomic distance between atoms in the same layer.
d0_01: float, optional
Interlayer distance.
q0_00: float, optional
Decay rate for the $pp\pi$ hopping integral.
q0_01: float, optional
Decay rate for the $pp\sigma$ hopping integral. It will be equal
to q0_00 by default, representing isotropic decay rates.
IP_cutoff: float, optional
Cutofff of the hopping between atoms in the same layer.
OP_cutoff: float, optional
Cutoff distance of the hopping between atoms in different layers.
dangling: float, optional
On-site energy shift applied to atoms in edges. These can be atoms in the edges
of the 1D graphene nanostructures or in the edges of pores in nanoporous graphene
nanostructures.
"""
if isinstance(geom,sisl.Geometry):
g = geom.copy()
else:
g = sisl.get_sile(geom).read_geometry()
g = g.remove([i for i in range(g.na) if g.atoms[i].Z == 1])
gtmp = sisl.Geometry(g.xyz,atoms=sisl.Atom(6,R=1.44),sc=g.sc) ### Test
Ham = sisl.Hamiltonian(gtmp) ### Test
Ham.set_nsc((3,3,1))
xyz = Ham.xyz
#d0_01 = g.xyz[:,2].max() - g.xyz[:,2].min()
if os_1 is None:
os_1 = os_0
os_1 += field*eps*d0_01
if q0_01 is None:
q0_01 = q0_00
def cart_to_sph(r):
# Wiki spherical coordinates
Res = np.zeros(r.shape)
R = np.sqrt(np.sum(r**2,axis=1))
theta = np.arctan2( np.sqrt(r[:,0]**2 + r[:,1]**2 ) , r[:,2] )
phi = np.arctan2( r[:,1] , r[:,0] )
Res[:,0] = R
Res[:,1] = theta
Res[:,2] = phi
return Res
def SK_pz(xi,xj):
rij = xj - xi
if len(rij.shape)==1:
rij = np.array([rij])
return _SK_pz(rij)
def _SK_pz(rij, OP_tol = .5):
V = np.zeros(len(rij))
sph = cart_to_sph(rij)
l = np.sin(sph[:,1])
d = sph[:,0]
OP_hop = ( d < OP_cutoff )*( np.abs(rij[:,2])>OP_tol )
idx_inter = OP_hop.nonzero()[0]
l = l[idx_inter] ### Test
d = d[idx_inter] ###
V[idx_inter] = Vpppi*l**2 * np.exp(-(d - d0_00)*q0_00) + Vpps * (1-l**2) * np.exp(-(d-d0_01)*q0_01) ###
return V
Ham.construct(([0.1,1.5],[0.,-2.7])) ### Test
Vals = []
Hops = [(0,0), (0,1),(1,0),(1,1),(-1,0),(0,-1),(-1,1),(1,-1),(-1,-1)]
R1 = Ham.cell[0]; R2 = Ham.cell[1]
for ij in Hops:
T = ij[0]*R1 + ij[1]*R2
V = np.zeros((len(xyz), len(xyz)))
for i in range(len(xyz)):
V [i,:] = SK_pz(xyz[i], xyz+T)
if ij == (0,0):
interz = np.average(xyz[:,2])
idx_0 = np.where(xyz[:,2]<interz)[0]
idx_1 = np.where(xyz[:,2]>interz)[0]
os = np.zeros(Ham.no)
os[idx_0] = os_0
os[idx_1] = os_1
for kk in range(len(xyz)):
xyz33 = g.tile(3,axis=1).tile(3,axis=0).translate(v=-g.cell[0]-g.cell[1]).xyz
if (np.linalg.norm(xyz[kk]-xyz33,axis=1)<1.6).sum()<4:
os[kk] += dangling
idx = np.arange(len(os))
V[idx,idx] = os
Vals += [V]
for i in range(len(Hops)):
ij = Hops[i]
v = Vals[i]
idx = np.where(v!=0)
vals = v[idx[0], idx[1]]
for ii in range(len(idx[0])):
I,J = idx[0][ii], idx[1][ii]
Ham[I, J, ij] = vals[ii]
# Fails with and without these two lines
Ham.eliminate_zeros()
Ham.finalize()
return Ham
H = tbbi(sisl.get_sile('g.xyz').read_geometry())
for i in range(H.no):
H[i,i] = 1e-6 + H[i,i]
H.write('out.TSHS')
@AleksBL it is fixed in the latest commit, some things have changed a bit.
Code to reproduce problem
Gives [1 1 1] SislError: The diagonal elements of your orthogonal Hamiltonian have not been defined, this is a requirement.
Redefining by
makes it write without error.
This bug has appeared at several updates it seems, or it might be some strange stuff going on in the dependencies. I dont know, but i think it was fixed in a 0.13 version.
I didnt understand the origin of it when it was fixed previously, so I wondered if you know where this error comes from?
Best, Aleksander
Version details Run the below code and add to bug-report:
3.10.12 (main, Nov 20 2023, 15:14:05) [GCC 11.4.0] 0.14.2