Return of missing diagonal elements error

[AleksBL]

Code to reproduce problem

print(Hd.nsc)
Hd.write('out.TSHS')

Gives
[1 1 1]
SislError: The diagonal elements of your orthogonal Hamiltonian have not been defined, this is a requirement.

Redefining by

Hd2 = sisl.Hamiltonian(sisl.Geometry(Hd.xyz, atoms=Hd.atoms))
import scipy.sparse as sp
Hd2 = Hd2.fromsp(sisl.Geometry(Hd.xyz, atoms=Hd.atoms), Hd.Hk(), S = sp.eye(Hd.na))
Hd2.write('out.TSHS')

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

[zerothi]

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?

[AleksBL]

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

[zerothi]

Could you create a complete example. Your example isn't fully complete.

[AleksBL]

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')

[zerothi]

@AleksBL it is fixed in the latest commit, some things have changed a bit.