Topological state not appearing

[NikFriedrich]

Trying to calculate the topological edge state of 2B-doped armchair graphene nanoribbon, the TB code seems to fail. The code used to be capable of calculating this state. I don't understand why it fails now.

Here is the script I used to generate the figure:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Feb 01 14:53:13 2022

@author: niklas
"""
import sisl
import numpy as np
import hubbard as h

def agnr2(w, bond=1.42):
    """ AGNR primitive cell of width `w` periodic along x-axis """
    g = sisl.geom.graphene(orthogonal=True).repeat(w//2 + w % 2, axis=1)
    if w % 2:
        g = g.remove(-w//2)
        g = g.remove(0)
    g = g.move(-g.center())
    g.set_nsc([3, 1, 1])
    return g

def agnr2B(w, n, bond=1.42, nB=2):
    """ Create an AGNR with (up to) two B-substitutions """
    g = agnr2(w, bond=bond).tile(n, axis=0)
    g = g.move(-g.center())
    # Find hexagon near origo
    idx = g.close(np.array([0, 0, 0]), R=[1.1*bond])
    # Set first and last atoms to B
    bond = 1.42
    B = sisl.Atom(Z=5, R=bond, orbs=1)
    # N = sisl.Atom(Z=7, R=bond, orbs=1)
    if nB > 0:
        g.atoms[idx[0]] = B
    if nB > 1:
        g.atoms[idx[-1]] = B
    return g

# For the sake of example, let's look at a small 2B-7AGNR
g = agnr2B(7, 15)
g.set_nsc([1, 1, 1])

# Build 1NN TB Hamiltonian
H = sisl.Hamiltonian(g, spin="polarized")
H.construct([[0.1, 1.6], [0, -3]])

# Taken from example sent by Thomas Fredriksen 02.12.2021
mixer = sisl.mixing.PulayMixer(0.5, history=7)
HH = h.HubbardHamiltonian(H, U=3.5)

# Set up spin polarization for topological states at end and 2B-site
HH.set_polarization([1, 113], [96, 208])
HH.converge(h.density.calc_n_insulator, 
            mixer=mixer, 
            steps=1, print_info=True)
p = h.plot.SpinPolarization(HH, colorbar=True, vmax=0.2)
p.savefig("no_topo_state.pdf")

no_topo_state.pdf

[tfrederiksen]

I think the problem lies in the construction of H which does not properly encode the boron atoms. The H.construct call here just sets onsite energies and nearest hoppings as if all sites were carbon.

Try changing the line to H = h.sp2(g).

[NikFriedrich]

Indeed, that did the trick. Thanks a lot.