Change energy reference

[tfrederiksen]

After discussions with @ROrtiz1993 I realize that it may be more natural to define the diagonal elements in the MFH Hamiltonian relative to the average charge per spin channel, i.e., q0 / 2.

As a consequence of this, the midgap energy for a graphene nanoflake remains at zero for a first-nearest neighbour model (1NN), independent of the onsite repulsion parameter U.

[tfrederiksen]

It might also be a good idea to update the documentation to reflect our choice of energy reference: https://dipc-cc.github.io/hubbard/docs/latest/index.html

[sofiasanz]

Hi, thank you! Have you tested it? Maybe the tests files should be updated too. I just run the test-hubbard-quick.py, test-hubbard-quantitative.py and they seem to not be working properly.

[zerothi]

doesn't the energy factor depend on whether you have 2 spins or not? I.e. are the occupations doubled so that the energy is the same, or not? I.e. that would be a good test ;)

[sofiasanz]

I tested the change of this branch for the case of the Clar goblet molecule.

With the 1NN model I get the following eigenvalues:

HOMO-1 -1.0562566111567124
HOMO -0.3635487720214219
LUMO 0.3635487720214266
LUMO+1 1.0562566111567155

With the 3NN model I get:

HOMO-1 -0.48132864342771287
HOMO 0.06883450639151843
LUMO 0.8078526711122002
LUMO+1 1.3925001558617789

Without this change (i.e. in the original implementation), for the 1NN I get:

HOMO-1 -2.806256611174714
HOMO -2.113548772030526
LUMO -1.3864512279707735
LUMO+1 -0.6937433888242578

Al values are given in eV. The Coulomb repulsion parameter is U=3.5 eV. I also ran a similar calculation for the spin-degenerate case and I got similar results.

However, I think that the fact that for the 1NN model the eigeinvalues are referenced to the midgap using the energy shift U*q0/2 instead shift U*q0 is a characteristic for bipartite lattices since I also tested our molecular junction from Nat. Commun. volume 10, Article number: 200 (2019) and for this molecule, using the 1NN I get:

HOMO-1 -0.4709273766507112
HOMO -0.2819595063689108
LUMO -0.09374137556195279
LUMO+1 0.4462153433372981

With this I think this change is OK, however we should update the test scripts since I quickly ran them and they don't work.

[tfrederiksen]

Thanks @sofiasanz for having a look at this.

What about the total energy? Did you check that we get the same result with spin-degenerate and polarised settings for a system that is unpolarized (as one should)?

Could you update tests and documentation before merge?

[sofiasanz]

Did you check that we get the same result with spin-degenerate and polarised settings for a system that is unpolarized (as one should)?

Yes, that works properly.

The thing is that for our test molecule, which is composed by carbon, boron and nitrogen atoms, the resulting densities obtained before and after this change are different since the charge associated to each atom in the geometry (self.q0) is different, and therefore the added quantity is not just a constant term added to the diagonal.

[zerothi]

Did you check that we get the same result with spin-degenerate and polarised settings for a system that is unpolarized (as one should)?

Yes, that works properly.

The thing is that for our test molecule, which is composed by carbon, boron and nitrogen atoms, the resulting densities obtained before and after this change are different since the charge associated to each atom in the geometry (self.q0) is different, and therefore the added quantity is not just a constant term added to the diagonal.

Is this because you essentially want to explicitly denote how many electrons in each of the spin channels? Or? Perhaps it would make sense if atoms/orbitals could contain spin charges that are different?

[tfrederiksen]

The thing is that for our test molecule, which is composed by carbon, boron and nitrogen atoms, the resulting densities obtained before and after this change are different since the charge associated to each atom in the geometry (self.q0) is different, and therefore the added quantity is not just a constant term added to the diagonal.

This sounds like a bug (one way or the other). We do not want the solution to depend on the chosen energy reference.

[sofiasanz]

Is this because you essentially want to explicitly denote how many electrons in each of the spin channels? Or? Perhaps it would make sense if atoms/orbitals could contain spin charges that are different?

Do you mean different for each spin channel? Hmmm I'm not sure if that makes much sense?

The problem is that we subtract self.q0 to the diagonal of the Hamiltonian, and the discussion of this branch is that this term doesn't make much sense because we are subtracting two times (one per spin Hamiltonian), so there was a factor 1/2 missing.

On the other hand, this factor (self.q0) was not only just an energy reference since it is a vector (is the charge associated to each orbital at each atomic site), so for heteroatomic structures self.q0[i] is different.

This sounds like a bug (one way or the other). We do not want the solution to depend on the chosen energy reference.

Yes, we should fix this! Maybe we should change self.q0 (charge associated to each orbital) by the total charge, no?

[sofiasanz]

I have updated this branch using the average charge per site as the energy reference. I also updated the density file containing the density of the test molecule (tests/mol-ref/density.nc)

[tfrederiksen]

Looks good, Sofia.

Note to @NikFriedrich , this fix will affect your calculations with B-atoms!

[sofiasanz]

Looks good, Sofia.

If you agree then I'll merge this branch :)

[sofiasanz]

Actually I have found one of our test files that is not working properly. In particular tests/test-hubbard-multi-orbitals.py. In this test we compare the results for a periodic ZGNR with only pz orbitals, with the case of a ZGNR where the carbon atoms have pz orbitals and a fictitious "s" orbital that have associated charge equal to zero and U_s = 0 (like a void orbital).

Result of the main branch:

Results of this branch:

[tfrederiksen]

Why is this new result incorrect? Aren't electrons simply transferring to the $s$-band? I guess you will recover the pristine ZGNR bands if you put the onsite $s$ well above that of $p_z$.

[sofiasanz]

The idea of this test was to have non interacting s-orbitals, like to have just more zeros in the Hamiltonian. I think that this should give the same result since the hoppings are only set between pz orbitals so the s orbitals are non interacting.

[tfrederiksen]

The idea of this test was to have non interacting s-orbitals, like to have just more zeros in the Hamiltonian. I think that this should give the same result since the hoppings are only set between pz orbitals so the s orbitals are non interacting.

Did you try what I suggest? Since states are filled from below I do not see why the electrons would not end up in the $s$-orbitals.

[sofiasanz]

Yes you were correct, now it works! What I don't understand is why it was working before...

This is what I get now:

The flat line position depends on the on-set value I set for the s-orbitals.