Problems on energy spectra and MFH-LDOS plotting

[lycheehoo]

Hi everyone, I am working on high-spin nanographene synthesized on metal surfaces using scanning tunneling microscopy (STM), and I am not good at programming, so I wonder can these codes could plot energy spectra and mean-field Hubbard local density of states (MFH-LDOS) maps shown below ([https://dx.doi.org/10.1021/acs.nanolett.0c02939],and then, how to run it? Thank you very much!

[sofiasanz]

Hi @lycheehoo, yes you can do these kind of plots with this module. Regarding the LDOS plots we are just making small adjustments to the code, but it should be solved soon.

For the energy spectra, you can easily obtain them. You can check the examples shown in the documentation webpage. For instance, here I leave you the example of how to solve the self-consistent (SCF) solution of a molecule with the mean-field Hubbard (MFH) model in the following link:

https://dipc-cc.github.io/hubbard/docs/latest/examples/molecules.html

Once you obtained the SCF solution, you can obtain the energy levels simply by diagonalizing the Hamiltonian:

ev, evec = H.eigh(spin=spin, eigvals_only=False)

Where spin=0 represents spin up and spin=1 represents spin down. ev and evec are the eigenvalues are eigenvectors (column vectors), respectively.

Then you can plot them, let's say, as in panel b:

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(4,8))
ax.hlines(ev, xmin=-0.1, xmax=0.1)
ax.set_xlim(-0.2,0.2)
ax.get_xaxis().set_visible(False)
plt.savefig('energy_levels.pdf')

And the density of states (DOS):

import numpy as np
egrid = np.linspace(-10, 10, 400) # Define your energy points at which you want to obtain the DOS
DOS = H.DOS(egrid, spin=spin, eta=eta) # Compute the DOS for the desired spin index along the energy grid

The parameter eta is the level broadening for the Lorentzian distribution. By default is set to 1e-3 eV, but you can change this value.

If you want the total DOS, you can just sum the DOS for spin up and the DOS for spin down.

[lycheehoo]

Hi Sofia, thanks for your detailed reply. I have already got the expected energy spectra following your kind guidance, but I still got 3 questions remained:

1. How to find the suitable parameters set in HH.set_polarization([2,6,10], dn=[24,28,32]) to get a faster convergence? Are there some rules for them?
2. How to plot these DOS data calculated above?
3. If I want to consider systems with B or N atoms, what are the recommed setting of the on-site energy for those atoms (eB and eN)?

[sofiasanz]

Hi Sofia, thanks for your detailed reply. I have already got the expected energy spectra following your kind guidance, but I still got 3 questions remained:

1. How to find the suitable parameters set in `HH.set_polarization([2,6,10], dn=[24,28,32])` to get a faster convergence? Are there some rules for them?

Unfortunately there is no recipe for knowing which is the best spin distribution guess for a system. However, you might get a good insight by looking at the low energy levels (i.e. those occupied/unoccupied close to the Fermi level), as these are the eigenstates that are going to participate more in the emergence of spin polarization. You can for instance look for localized states and where (in which sites of your molecule) they localize. Your spin polarization will most likely look similar to that. If you take for instance the examples molecules in the documentation there you can compare the wavefunctions and the spin polarization, maybe from that you can get an idea of how it works.

To plot the wavefunctions you can just to:

from hubbard import plot
WF = evec[:, int(H.sites/2)] # The second index is for choosing the desired eigenstate. For this index WF stands for the LUMO, int(H.sites/2) -1 would be HOMO, and so forth (for a half-filled molecule)
p = plot.Wavefunction(H, WF*f)
p.savefig('wavefunction.pdf')

The f factor is just for visualization purposes, since probably the weight of the wave function is just too small to see without an amplification factor. I would try with f>1000 (for larger molecules you should use a larger f).

2. How to plot these DOS data calculated above?

You can just simply use matplotlib to plot the DOS against the energy grid:

import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.plot(egrid, DOS)
ax.set_xlabel('Energy [eV]')
ax.set_ylabel('DOS')
plt.savefig('DOS.pdf')

Is this what you are looking for?

3. If I want to consider systems with B or N atoms,  what are the recommed setting of the on-site energy for those atoms (`eB` and `eN`)?

We have a suggestion for these parameters in our module sp2, which you can look in the documentation. I cannot say much about these values since I have not carried many calculations with B or N. If you want to check weather your parameters are a good fit, I would suggest to take a simple system and run both a calculation with DFT and MFH and check if the parameters that you are using look like a good fit.

[lycheehoo]

I got it, thank you very much, Sofia. Hope you everything goes well!