Can Hubbard.plot plot molecule containing oxygen atoms?

[yuanzhangyu]

Thanks very much for reading my issues,

I looked through the hubbard/hubbard/plot/plot.py and found it seems the molecule geometry containing O cannot be plotted by Hubbard. plot. Then I added

elif g.atoms[ia].Z == 8: # O self.axes.add_patch(patches.Circle((np.average(g.xyz[ia, 0]), np.average(g.xyz[ia, 1])), radius=0.7, color='g', lw=2, fill=False)) pi.append(patches.Circle((g.xyz[ia, 0], g.xyz[ia, 1]), radius=0.7))

but still, the molecule geometry containing O cannot be obtained.

I am just wondering is it possible to get the results(geometry plot, spin density) for the molecule containing O?

[sofiasanz]

Dear @yuanzhangyu, thank you very much for your question. On one hand, if you calculated the TB Hamiltonian using our module sp2, the resulting Hamiltonian does not include the O atoms, since we don't have this implemented for the moment. For this reason, the geometry saved in theHubbardHamiltonian object does not include the O atoms. For now, to have the full TB Hamiltonian you should add these atoms by hand, building your own Hamiltonian. I leave here an example of how I might do that with sisl:

import sisl

# Assuming you want a spin-polarized Hamiltonian
Hsp2 = sisl.Hamiltonian(geometry, spin='polarized') 

for ia in Hsp2.geometry:
     idx = Hsp2.geometry.close(ia, R=[0.1, 1.6]) # Look for 1NN atoms around ia

     if Hsp2.geometry.atoms[ia].Z == 6: # Carbon atoms
          Hsp2[ia, ia, :] = eC # set onsite for C sites
          # Set hopping parameters
          # Look for atoms around atom ia and check if they are oxygen or carbon atoms
          for ib in idx[1]:
               if Hsp2.geometry.atoms[ib].Z == 8:
                    Hsp2[ia, ib, :] = -t_CO # Hopping between carbon atom and oxygen
               else: # Assuming your system only has C or O atoms, otherwise one has to add more if statements
                    Hsp2[ia, ib, :] = -t_CC # Hopping between carbon-carbon atoms

     elif Hsp2.geometry.atoms[ia].Z == 8: # Oxygen atoms
          Hsp2[ia, ia, :] = eO # set onsite for O sites
          # Set hopping parameters
          # Look for atoms around atom ia and check if they are oxygen or carbon atoms
          for ib in idx[1]:
               if Hsp2.geometry.atoms[ib].Z == 6:
                    Hsp2[ia, ib, :] = -t_CO # Hopping between oxygen atom and carbon
               else: # Assuming your system only has C or O atoms, otherwise one has to add more if statements
                    Hsp2[ia, ib, :] = -t_OO # Hopping between oxygen-oxygen atoms

Where geometry represents your sisl.Geometry object containing only the sp2 pi-like atoms, ie, it shouldn't contain hydrogen atoms or sp3 atoms, etc. Note tat I'm assuming that there is only one orbital per atom. Maybe there is a better way to do this, but I hope it serves as an example (check that this code gives a hermitian Hamiltonian).

With this sisl Hamiltonian you can now build your HubbardHamiltonian which will contain the Oxygen atoms.

Then you can add the lines you wrote:

elif g.atoms[ia].Z == 8: # O
    self.axes.add_patch(patches.Circle((np.average(g.xyz[ia, 0]), np.average(g.xyz[ia, 1])), radius=0.7, color='g', lw=2, fill=False))
    pi.append(patches.Circle((g.xyz[ia, 0], g.xyz[ia, 1]), radius=0.7))

in the hubbard/hubbard/plot/plot.py, and you should be able to plot the solution for the full system.

[yuanzhangyu]

OK! Thanks a lot for your information!!

Dear @yuanzhangyu https://github.com/yuanzhangyu, thank you very much for your question. On one hand, if you calculated the TB Hamiltonian using our module sp2, the resulting Hamiltonian does not include the O atoms, since we don't have this implemented for the moment. For this reason, the geometry saved in theHubbardHamiltonian object does not include the O atoms. For now, to have the full TB Hamiltonian you should add these atoms by hand, building your own Hamiltonian. I leave here an example of how I might do that with sisl https://sisl.readthedocs.io/en/latest/introduction.html:

import sisl

Assuming you want a spin-polarized HamiltonianHsp2 = sisl.Hamiltonian(geometry, spin='polarized')

for ia in Hsp2.geometry: idx = Hsp2.geometry.close(ia, R=[0.1, 1.6]) # Look for 1NN atoms around ia

 if Hsp2.geometry.atoms[ia].Z == 6: # Carbon atoms
      Hsp2[ia, ia, :] = eC # set onsite for C sites
      # Set hopping parameters
      # Look for atoms around atom ia and check if they are oxygen or carbon atoms
      for ib in idx[1]:
           if Hsp2.geometry.atoms[ib].Z == 8:
                Hsp2[ia, ib, :] = -t_CO # Hopping between carbon atom and oxygen
           else: # Assuming your system only has C or O atoms, otherwise one has to add more if statements
                Hsp2[ia, ib, :] = -t_CC # Hopping between carbon-carbon atoms

 elif Hsp2.geometry.atoms[ia].Z == 8: # Oxygen atoms
      Hsp2[ia, ia, :] = eO # set onsite for O sites
      # Set hopping parameters
      # Look for atoms around atom ia and check if they are oxygen or carbon atoms
      for ib in idx[1]:
           if Hsp2.geometry.atoms[ib].Z == 6:
                Hsp2[ia, ib, :] = -t_CO # Hopping between oxygen atom and carbon
           else: # Assuming your system only has C or O atoms, otherwise one has to add more if statements
                Hsp2[ia, ib, :] = -t_OO # Hopping between oxygen-oxygen atoms

Where geometry represents your sisl.Geometry object containing only the sp2 pi-like atoms, ie, it shouldn't contain hydrogen atoms or sp3 atoms, etc. Note tat I'm assuming that there is only one orbital per atom. Maybe there is a better way to do this, but I hope it serves as an example (check that this code gives a hermitian Hamiltonian).

With this sisl Hamiltonian you can now build your HubbardHamiltonian which will contain the Oxygen atoms.

Then you can add the lines you wrote:

elif g.atoms[ia].Z == 8: # O self.axes.add_patch(patches.Circle((np.average(g.xyz[ia, 0]), np.average(g.xyz[ia, 1])), radius=0.7, color='g', lw=2, fill=False)) pi.append(patches.Circle((g.xyz[ia, 0], g.xyz[ia, 1]), radius=0.7))

in the hubbard/hubbard/plot/plot.py, and you should be able to plot the solution for the full system.

— Reply to this email directly, view it on GitHub https://github.com/dipc-cc/hubbard/issues/131#issuecomment-1380193156, or unsubscribe https://github.com/notifications/unsubscribe-auth/A45LDZGKSOO7ZZYYH4CGFKTWR7TRVANCNFSM6AAAAAATWN3EVI . You are receiving this because you were mentioned.Message ID: @.***>

[yuanzhangyu]

Thanks very much for your information!

I modified the sp2.py like this, ...

Determine pz sites

aux = []
sp3 = []
for ia, atom in enumerate(ext_geom.atoms.iter()):
    # Append non C-type atoms in aux list
    if atom.Z not in [6,8]:
        aux.append(ia)
    idx = ext_geom.close(ia, R=[0.1, 1.6])
    if len(idx[1]) == 4: # Search for atoms with 4 neighbors
        if atom.Z == 6:
            sp3.append(ia)

... Hsp2 = sisl.Hamiltonian(pi_geom, orthogonal=orthogonal,spin=spin)

for ia in pi_geom:
    idx = pi_geom.close(ia, R=[0.1, 1.6]) # Look for 1NN atoms around ia
    # NB: I found that ':' is necessary in the following lines, but I don't understand why...
    if pi_geom.atoms[ia].Z == 6: # Carbon atoms
        Hsp2[ia, ia, :] = eC # set onsite for C sites
        # Set hopping parameters
        # Look for atoms around atom ia and check if they are oxygen or carbon atoms
        for ib in idx[1]:
            if pi_geom.atoms[ib].Z == 8:
                Hsp2[ia, ib, :] = -t1 # Hopping between carbon atom and oxygen
            else: # Assuming your system only has C or O atoms, otherwise one has to add more if statements
                Hsp2[ia, ib, :] = -t2 # Hopping between carbon-carbon atoms
    elif pi_geom.atoms[ia].Z == 8: # Oxygen atoms
        Hsp2[ia, ia, :] = eO # set onsite for O sites
        # Set hopping parameters
        # Look for atoms around atom ia and check if they are oxygen or carbon atoms
        for ib in idx[1]:
            if pi_geom.atoms[ib].Z == 6:
                Hsp2[ia, ib, :] = -t1 # Hopping between oxygen atom and carbon
            else: # Assuming your system only has C or O atoms, otherwise one has to add more if statements
                Hsp2[ia, ib, :] = -t3 # Hopping between oxygen-oxygen atoms
    if not Hsp2.orthogonal:
        Hsp2[ia, ia]=1.0
        Hsp2[ia, idx[1]] = s1
return Hsp2

and also modified the plot.py, but it doesn't seem to work. Did I do some wrong things?

[sofiasanz]

What error are you getting? Which values for the needed parameters are you using?

[yuanzhangyu]

Thanks for your reply,

I didn’t get an error, actually, I just used hubbard.plot to draw the structure in Ref. Nano Lett. 2022, 22, 1, 164–171 https://pubs.acs.org/doi/full/10.1021/acs.nanolett.1c03578. But I didn’t get the O atom.

The parameters in the SI of the reference are used in sp2.py.

What error are you getting? Which values for the needed parameters are you using?

— Reply to this email directly, view it on GitHub https://github.com/dipc-cc/hubbard/issues/131#issuecomment-1381769871, or unsubscribe https://github.com/notifications/unsubscribe-auth/A45LDZBXKHPFDEWFJCEIU43WSFBA7ANCNFSM6AAAAAATWN3EVI . You are receiving this because you were mentioned.Message ID: @.***>

[sofiasanz]

Can you pass your geometry? In format .xyz for example?

[yuanzhangyu]

OK! Of course! Attached please find it, thanks so much!

Geometry_Test.xyz.zip

[sofiasanz]

Dear @yuanzhangyu, sorry for the delay in my answer. Here I leave something that works (it sets correctly the tight-binding Hamiltonian and it plots the spin polarization for the self-consistent solution for that geometry taking into account the O atoms):

import numpy as np
import sisl
from hubbard import HubbardHamiltonian, density, plot

geometry = sisl.get_sile('Geometry_Test.xyz').read_geometry()

# Set tight-binding Hamiltonian for the sp2 system
Hsp2 = sisl.Hamiltonian(geometry, spin='polarized') 
eC, eO, t_CC, t_CO = 0, -1.5, 2.7, 3.8 # eV (parameters taken from ref. Ref. Nano Lett. 2022, 22, 1, 164–171 <https://pubs.acs.org/doi/full/10.1021/acs.nanolett.1c03578>)

for ia in Hsp2.geometry:
     idx = Hsp2.geometry.close(ia, R=[0.1, 1.6]) # Look for 1NN atoms around ia

     if Hsp2.geometry.atoms[ia].Z == 6: # Carbon atoms
          Hsp2[ia, ia, :] = eC # set onsite for C sites
          # Set hopping parameters
          # Look for atoms around atom ia and check if they are oxygen or carbon atoms
          for ib in idx[1]:
               if Hsp2.geometry.atoms[ib].Z == 8:
                    Hsp2[ia, ib, :] = -t_CO # Hopping between carbon atom and oxygen
               else: # Assuming your system only has C or O atoms, otherwise one has to add more if statements
                    Hsp2[ia, ib, :] = -t_CC # Hopping between carbon-carbon atoms

     elif Hsp2.geometry.atoms[ia].Z == 8: # Oxygen atoms
          Hsp2[ia, ia, :] = eO # set onsite for O sites
          # Set hopping parameters
          # Look for atoms around atom ia and check if they are oxygen or carbon atoms
          for ib in idx[1]:
               if Hsp2.geometry.atoms[ib].Z == 6:
                    Hsp2[ia, ib, :] = -t_CO # Hopping between oxygen atom and carbon

# Plot the Hamiltonian
p = plot.BondHoppings(Hsp2, off_diagonal_only=False)
p.legend()
p.savefig('bonds.pdf')

# Build HubbardHamiltonian object and converge
H = HubbardHamiltonian(Hsp2, U=3)
H.set_polarization([44], [27])
H.converge(density.calc_n, tol=1e-10)

p = plot.SpinPolarization(H, vmin=-0.2, vmax=0.2, colorbar=True)
p.savefig('polarization.pdf')

And by adding the line that you mentioned hubbard/hubbard/plot/plot.py :

elif g.atoms[ia].Z == 8: # O
    self.axes.add_patch(patches.Circle((np.average(g.xyz[ia, 0]), np.average(g.xyz[ia, 1])), radius=0.7, color='g', lw=2, fill=False))
    pi.append(patches.Circle((g.xyz[ia, 0], g.xyz[ia, 1]), radius=0.7))

Beware that this is an example script and you should check the quantitative and qualitative result that is giving, I haven't checked that this is physically correct, just making sure that the program works for your purposes.

In this case, I would also recommend that you check that the tight-binding Hamiltonian is hermitian, for instance by checking the following difference:

Hsp2.Hk(format='array') - np.conjugate(Hsp2.Hk(format='array')).T

By the way, you should re-install your local hubbard module if you make changes in the package (for instance if you change the lines in hubbard/hubbard/plot/plot.py you should do again python -m pip install . to update the the package with your changes).

[yuanzhangyu]

OK, I see ! Thanks very much for your information ! This issue can be closed now !