fix: use correct form of pz orbital

[sofiasanz]

I think we were using the incorrect pz orbital for carbon atoms. With this fix we use sisl.HydrogenicOrbital instead.

[sofiasanz]

When we expand the wavefunctions in the real space grid we should take into account that depending on the definition of the basis orbitals they may overlap with each other. Since we typically assume orthogonal basis, maybe we should chose a reasonable Zeff by default that gives little overlap between the orbitals so we can still assume orthogonality.

Here I leave a comparison of LDOS plots for the three first states of benzene using three different Zeff (which in the end will result in more or less localized pz orbitals).

The Hamiltonian corresponds to the 3NN model where the hopping between first, second and third nearest neighbors are $t_1=2.7, t_2 =0.2, t_3=0.18$ eV.

Here I also test what is the norm of the expanded wavefunction by computing the absolute square of the real space grid for each eigenstate. We can see that the more confined the pz orbitals are, the closer to 1.0 the norm of the real space wavefunction is.

Z_eff = 3.2
state 0, Grid ** 2 =  1.716084494207556
state 1, Grid ** 2 =  0.9848432517472101
state 2, Grid ** 2 =  1.2098811567433554
state 3, Grid ** 2 =  0.7235034541768871
state 4, Grid ** 2 =  0.8579891599673729
state 5, Grid ** 2 =  0.5141571166414545

Z_eff = 4.0
state 0, Grid ** 2 =  1.3426258715268453
state 1, Grid ** 2 =  1.0187585113982105
state 2, Grid ** 2 =  1.122284396475522
state 3, Grid ** 2 =  0.85941107930204
state 4, Grid ** 2 =  0.9387608837177124
state 5, Grid ** 2 =  0.7173649340722447

Z_eff = 5.0
state 0, Grid ** 2 =  1.1355326400701464
state 1, Grid ** 2 =  1.0137233632460871
state 2, Grid ** 2 =  1.0537180125542505
state 3, Grid ** 2 =  0.9418538694872627
state 4, Grid ** 2 =  0.9773910904344912
state 5, Grid ** 2 =  0.8732008627789625

[zerothi]

Would it make sense to calculate an effective Z based on a radius? In this way you sort of go around the problem of fine-tuning Z, but instead you explicitly set R?

Have you checked what maxR is after using these Zeff?

[sofiasanz]

Would it make sense to calculate an effective Z based on a radius? In this way you sort of go around the problem of fine-tuning Z, but instead you explicitly set R?

Yes that could be an option, to take into account the distance between the orbitals in the lattice and try to make them "orthogonal".

Have you checked what maxR is after using these Zeff?

But isn't maxR an external parameter?

[zerothi]

Would it make sense to calculate an effective Z based on a radius? In this way you sort of go around the problem of fine-tuning Z, but instead you explicitly set R?

Yes that could be an option, to take into account the distance between the orbitals in the lattice and try to make them "orthogonal".

I don't know exactly the implications, it was just an idea to have a "better" control of the visual appearance.

Have you checked what maxR is after using these Zeff?

But isn't maxR an external parameter?

Once you have created an AtomicOrbital or an HydrogenicOrbital the maxR will be defined based on the function. It all depends on the arguments at the instantiation point. For instance:

h1 = HydrogenicOrbital(..., R=1)
h10 = HydrogenicOrbital(..., R=10) # the default
h = HydrogenicOrbital(..., R=-1) # let sisl optimize (see below)

Actually, @tfrederiksen hard-coded the default max-R which might be a bit unfortunate. Should we change to something else? Or just go back to sisl's default?

def radial(r):
    z = 2 * Z / (n * Bohr)
    zr = z * r
    pref = (z ** 3 * factorial(n - l - 1) / (2 * n * factorial(n + l)) ) ** 0.5
    L = eval_genlaguerre(n - l - 1, 2 * l + 1, zr)
    return pref * np.exp(-zr / 2) * zr ** l * L

The default of sisl (when R is not specified) is to find the "largest R where the function goes to 0 in a range of 0<50 Angstrom". The 50 Ang is just emperical, and should probably be an argument. So it depends.

[zerothi]

Actually, I think the current default of the functional is a bit unfortunate, so I think we need to change it for consistency. I.e. the radial function is based only on 0<10 Ang (if R isn't specified). This may not be good enough for very long range functions, the above should partly fix this by allowing sisl to optimize R up to 50 Ang.

[sofiasanz]

OK I see, I just tried setting R=-1 and printing the max(r) (which if I understood it correctly gives the maximum radius that sisl optimized) for the different Zeff, and I always get max(r)=10.0, but it might give a larger value if the default was larger? Or?

[zerothi]

OK I see, I just tried setting R=-1 and printing the max(r) (which if I understood it correctly gives the maximum radius that sisl optimized) for the different Zeff, and I always get max(r)=10.0, but it might give a larger value if the default was larger? Or?

Ah, I see. This is because the R field is maximally 10 Ang, it definitely needs changing. :) I'll change it to adapt based on Z only. But I can see that the actual check is also hard. We might need to change some logic in sisl for this to actually be properly contained.

[zerothi]

I will try and fix this in sisl, the way it is done is perhaps not the best one.

I will probably implement something like this:

# forcefully set the length to 10 Ang
HydrogenicOrbital(..., R=10)
# determine R based on the integrand contained 99.99% of the total integrand
HydrogenicOrbital(..., R=-0.9999)

how would you think about this?

[sofiasanz]

I will try and fix this in sisl, the way it is done is perhaps not the best one.

I will probably implement something like this:

# forcefully set the length to 10 Ang
HydrogenicOrbital(..., R=10)
# determine R based on the integrand contained 99.99% of the total integrand
HydrogenicOrbital(..., R=-0.9999)

how would you think about this?

What happens if one uses another negative radius? Would it determine R based on the integrand contained, let's say, 70% of the total integrand (using R=-0.7)?

[zerothi]

What happens if one uses another negative radius? Would it determine R based on the integrand contained, let's say, 70% of the total integrand (using R=-0.7)?

Yes, then it is 70%. The main idea would be to allow users to fine-tune this. Instead of hard-coding this in sisl. For instance a very slowly decaying function might be annoying to have a long tail, yeah, you'll have an abrupt discontinuity, but...

[sofiasanz]

What happens if one uses another negative radius? Would it determine R based on the integrand contained, let's say, 70% of the total integrand (using R=-0.7)?

Yes, then it is 70%. The main idea would be to allow users to fine-tune this. Instead of hard-coding this in sisl. For instance a very slowly decaying function might be annoying to have a long tail, yeah, you'll have an abrupt discontinuity, but...

Yes, totally agreed. I think then it looks nice! In this way it allows the user to set the desired minimum value for this function, which I think is nice. I have no better suggestions for this flag for the moment ;)

[sofiasanz]

In my last commit (bc85fae) I set an argument that allows the user to control the atoms in the geometry when building the sp2 Hamiltonian, and therefore the orbitals inside these objects.

[sofiasanz]

Yes, I think this is better. :)

Excellent, thanks :)

I'll clean and merge this PR then :)