cleaned the procedure for locating maxR #575

zerothi · 2023-05-31T07:50:18Z

Instead of searching for a 0 in the radial function we now search for an integration that contians 99.99% of the radial function.
The integration uses the trapezoid integration mechanism with f2*r2 as that should be normalized for a
radial function.

This should make it much more stable against different radial function arguments and produce more accurate ranges for the orbitals.

Now a warning will be issued when the routine cannot determine a correct R (based on inputs).

Ensured that all orbitals with radial functions
can pre-calculate the R in case user did not supply it. One can always forcefully set it by supplying it.

@sofiasanz @tfrederiksen could you please have a look here, this should better adapt the R discovery of your HydrogenicOrbital

Closes Orbital set_radial and friends -- allow custom maximum R finding #574
Tests added
Documentation for functionality
User visible changes (including notable bug fixes) are documented in CHANGELOG

sofiasanz · 2023-05-31T09:26:36Z

Nice!

For the hyrogenic orbital I was using for a carbon pz orbital I get the following numbers using different Zeff:

pz = sisl.HydrogenicOrbital(2, 1, 0, Zeff, q0=1.0, R=-0.999999)

Output (maxR and the radial function evaluated at that point):

****** Zeff = 3.2 ******
pz.R =  3.8761000000000556
pz.radial([pz.R]) =  [0.00057862]
****** Zeff = 4 ******
pz.R =  3.101200000000003
pz.radial([pz.R]) =  [0.00080775]
****** Zeff = 5 ******
pz.R =  2.4820000000000677
pz.radial([pz.R]) =  [0.0011238]

If I don't set the R in pz (i.e. the default value set in sisl), I get:

****** Zeff = 3.2 ******
pz.R =  2.9418000000000886
pz.radial([pz.R]) =  [0.00740379]
****** Zeff = 4 ******
pz.R =  2.3537000000000083
pz.radial([pz.R]) =  [0.0103381]
****** Zeff = 5 ******
pz.R =  1.8838999999999964
pz.radial([pz.R]) =  [0.01439109]

zerothi · 2023-05-31T09:29:01Z

Nice!

For the hyrogenic orbital I was using for a carbon pz orbital I get the following numbers using different Zeff:

pz = sisl.HydrogenicOrbital(2, 1, 0, Zeff, q0=1.0, R=-0.999999)

Output (maxR and the radial function evaluated at that point):

****** Zeff = 3.2 ******
pz.R =  3.8761000000000556
pz.radial([pz.R]) =  [0.00057862]
****** Zeff = 4 ******
pz.R =  3.101200000000003
pz.radial([pz.R]) =  [0.00080775]
****** Zeff = 5 ******
pz.R =  2.4820000000000677
pz.radial([pz.R]) =  [0.0011238]

If I don't set the R in pz (i.e. the default value set in sisl), I get:

****** Zeff = 3.2 ******
pz.R =  2.9418000000000886
pz.radial([pz.R]) =  [0.00740379]
****** Zeff = 4 ******
pz.R =  2.3537000000000083
pz.radial([pz.R]) =  [0.0103381]
****** Zeff = 5 ******
pz.R =  1.8838999999999964
pz.radial([pz.R]) =  [0.01439109]

I would discourage overwriting pz.R like that, that is not intended to work like that. Either

explicitly set R when creating the argument
let sisl handle it and create an R for you

zerothi · 2023-05-31T09:30:16Z

I would have somehow thought the radial function was lower at that percentage integral. I think the normalization is OK. :)

zerothi · 2023-05-31T09:39:49Z

Ah, now I see what you did. :)
Great! Hope the finetuning is OK. :)

Any other suggestions here?

sofiasanz · 2023-05-31T09:40:11Z

I would discourage overwriting pz.R like that, that is not intended to work like that. Either
1. explicitly set R when creating the argument

2. let `sisl` handle it and create an `R` for you

Oh sorry, I'm just printing information, not overwriting (just wanted to see what is the maxR sisl is finding).

I would have somehow thought the radial function was lower at that percentage integral. I think the normalization is OK. :)

Yes, me too, but I guess is just because we are looking for the R that gets these values for the squared radial function, which is much more localized.

sofiasanz · 2023-05-31T09:43:10Z

Ah, now I see what you did. :) Great! Hope the finetuning is OK. :)

Any other suggestions here?

Not for the moment! For me it looks fine! Let's see if @tfrederiksen has any other suggestions.

codecov · 2023-05-31T10:45:48Z

Codecov Report

Merging #575 (07c5743) into main (9e4307b) will decrease coverage by 0.04%.
The diff coverage is 92.10%.

❗ Current head 07c5743 differs from pull request most recent head 37e073b. Consider uploading reports for the commit 37e073b to get more accurate results

@@            Coverage Diff             @@
##             main     #575      +/-   ##
==========================================
- Coverage   87.34%   87.31%   -0.04%     
==========================================
  Files         374      374              
  Lines       47334    47420      +86     
==========================================
+ Hits        41345    41403      +58     
- Misses       5989     6017      +28

Impacted Files	Coverage Δ
src/sisl/orbital.py	`90.80% <87.50%> (-1.15%)`	⬇️
src/sisl/tests/test_orbital.py	`100.00% <100.00%> (ø)`

... and 1 file with indirect coverage changes

sofiasanz · 2023-05-31T14:59:46Z

Hi Nick, I think there is something weird now in this branch.

I was plotting the LDOS maps for benzene using this branch (commit f89b533) and I got these following figures which look very weird to me (the top row are the wavefunctions, just to compare with the obtained LDOS):

I am using this orbital

pz = sisl.HydrogenicOrbital(2, 1, 0, Zeff, q0=1.0, R=10)

to expand the lattice in real space.

I compared to what I get with the main branch (using the same grid shape and everything), using version with commit 9e4307b, and I get the following images:

It looks like the LDOS changed a lot and also the resolution, which I don't understand.

zerothi · 2023-05-31T18:03:33Z

Thanks for this @sofiasanz

Could you try the latest commit? It should make R bigger. I believe this branch should be fine to check against if you simply use an explicit R, i.e. if you use R=10 you should get the same as main, that should be a good sanity check, after all.

I believe the problem arose because the integrated quantity was the actual radial integral. Whereas what you see in the wavefunctions are the radial function as is.

Now intsead of integrating and finding R for the integral (radial(r) * r)**2 I now do it for radial(r)*r which is not as quickly approaching 0, hence bigger contributions at further distances. I could also have just used radial(r) however, that might also be too fastly decaying. With the factor r it gets a bit longer. The proper way would be some other scheme, I guess something with log's`, but... That's probably it for another time.

zerothi · 2023-05-31T18:20:40Z

Hmm... Perhaps wait a bit, I need to check some things.

tfrederiksen · 2023-05-31T19:13:04Z

I like the idea that the user can set a cutoff radius through the specification of a negative R (which is then understood as the target volume). My suggestion would be that any R<0 could in principle be acceptable, for instance if the orbital is not normalized.

By the way, one could also have a function that normalizes the orbital on r<=R.

tfrederiksen · 2023-06-01T07:54:31Z

I don't understand the rationale behind radial**2 * r**3. The normalisation condition is int dr radial**2 * r**2 = 1, no?

zerothi · 2023-06-01T08:05:29Z

I don't understand the rationale behind radial**2 * r**3. The normalisation condition is int dr radial**2 * r**2 = 1, no?

Yes, this is where I am currently investigating the problems.
I'll soon have another commit to this branch which will highlight the problem of doing the actual radial integral. Hold tight.

tfrederiksen · 2023-06-01T08:12:10Z

I don't understand the rationale behind radial**2 * r**3. The normalisation condition is int dr radial**2 * r**2 = 1, no?

Yes, this is where I am currently investigating the problems. I'll soon have another commit to this branch which will highlight the problem of doing the actual radial integral. Hold tight.

Isn't it not just to solve for $R$ such that $\int_0^R dr f(r) ^2 r^2 =$ contains?

zerothi · 2023-06-01T08:21:50Z

As @sofiasanz pointed out, the actual function looked a bit weird in the real-space pictures.

This I think can be traced to the short R returned from the radial(r)**2 * r**2 integral. It simply yields a too short R to make sense for a proper real space visualization. This is because the integral converges much faster than the radial function does. Hence some other integral might be more appropriate.

I have tried various things, but I think the simplest and acceptable version would be the abs(radial(r)) function. It is a bit weird, but should suffice.

I have now opted for a different path where users can pass arguments to control the minimization procedure.

This can effectively be done like this:

R = {
"contains": 0.9999, # contain 99.99% of func integrated
"func": lambda radial, r: radial(r)**2 * r**2, # function to integrate
"maxR": 100, # max searched R distance (integral *must* converge within this range).
}
HydrogenicOrbital(..., R=R)

One can still supply a negative number which then translates to {"contains": - R} internally.
The function used is radial_minimize_range which is now also exposed through the API. It will allow users to test and play around with things. Comments on the function name and interface would be appreciated.

Please see the latest commit.

zerothi · 2023-06-01T08:26:17Z

@sofiasanz could you try with the latest commit to see if the new default is more appropriate?

sofiasanz · 2023-06-01T10:07:06Z

I like it :) I think also is much clearer and flexible now.

I'll try this commit in the afternoon and plot again the benzene LDOS to compare to what we get with the main branch. What is the default value for R now?

zerothi · 2023-06-01T11:47:05Z

I like it :) I think also is much clearer and flexible now.

I'll try this commit in the afternoon and plot again the benzene LDOS to compare to what we get with the main branch. What is the default value for R now?

for Z in [3.2, 4, 5]:
    orb = HydrogenicOrbital(2, 1, 0, Z)#, R=-0.999999)
    R = np.arange(0, orb.R, 0.001)
    r = orb.radial(R)
    print("99.99", orb.R, r[-1] / r.max()*100)

    orb = HydrogenicOrbital(2, 1, 0, Z, R=-0.999999)
    R = np.arange(0, orb.R, 0.001)
    r = orb.radial(R)
    print("99.9999", orb.R, r[-1] / r.max()*100)

yields:

99.99   3.888400000000018 0.025069398892111586
99.9999 5.519599999999978 0.0002567943876730368
99.99   3.110700000000002 0.025104120475539306
99.9999 4.415699999999987 0.0002569771567260434
99.99   2.4886000000000186 0.025104120475539272
99.9999 3.5326000000000057 0.00025697715672604347

Hope this should be fine enough.

The last value shows the percentage of the maximum radial value at the maximum R distance. I.e. 0.025 means that orb.radial([orb.R]) == orb.radial(<peak>) * 0.00025 (0.025%)

sofiasanz · 2023-06-01T15:29:15Z

I produced this image (using pz = sisl.HydrogenicOrbital(2, 1, 0, Zeff, q0=1.0, R=10)) with this branch (f62c12d) which looks like what I got with the main branch yesterday, so it seems to be working just fine:

zerothi · 2023-06-01T17:30:39Z

But what if you don't specify R? That was the idea.

zerothi · 2023-06-01T18:21:19Z

Or perhaps in another way, at what minimum R would you find the plots acceptable? :)

zerothi · 2023-06-06T10:20:36Z

@sofiasanz thanks for your check, I have now added a new commit which should fix the problem, now R is copied explicitly for all orbitals.

sofiasanz · 2023-06-06T11:31:33Z

Yes, I wrote a message before saying that when printing the Grid and Geometry objects the maxR didn't match, but I deleted it since I realized that the problem was because I was copying the geometry and then the information about the maxR was somehow lost. I'll check again with the new commits :)

sofiasanz · 2023-06-06T11:56:02Z

Here I have the images I get with this last commit for Zeff=3.2 and different R values:

I start to get a little weird images in the vicinity of R~2.5 Ang for this effective nuclear charge.
I think now is working fine :)

And BTW now, even though I use a copy of the geometry the information of the maxR is preserved 👍

zerothi · 2023-06-07T07:24:04Z

Thanks for checking.

Just for completeness:

import numpy as np
from sisl import HydrogenicOrbital

r = np.arange(0, 10, 0.001)
for Z in [3.2, 4, 5]:
    orb = HydrogenicOrbital(2, 1, 0, Z)
    R = np.arange(0, orb.R, 0.001)
    r = orb.radial(R)
    print("99.99  ", Z, orb.R, r[-1] / r.max()*100)

    orb = HydrogenicOrbital(2, 1, 0, Z, R=-0.999999)
    R = np.arange(0, orb.R, 0.001)
    r = orb.radial(R)
    print("99.9999", Z, orb.R, r[-1] / r.max()*100)

which produces

99.99   3.2 3.888400000000018 0.025069398892111586
99.9999 3.2 5.519599999999978 0.0002567943876730368
99.99   4 3.110700000000002 0.025104120475539306
99.9999 4 4.415699999999987 0.0002569771567260434
99.99   5 2.4886000000000186 0.025104120475539272
99.9999 5 3.5326000000000057 0.00025697715672604347

meaning that the default value produces a R=3.888 which should be fine then. I'll merge this then. :)

Instead of searching for a 0 in the radial function we now search for an integration that contians 99.99% of the radial function. The integration uses the trapezoid integration mechanism with f(r)**2*r**3 as that should be normalized for a radial function. This should make it much more stable against different radial function arguments and produce more accurate ranges for the orbitals. Now a warning will be issued when the routine cannot determine a correct R (based on inputs). Ensured that all orbitals with radial functions can pre-calculate the R in case user did not supply it. One can always forcefully set it by supplying it. Signed-off-by: Nick Papior <[email protected]>

This way we can converge faster to a better approach than the current absolute value of the functional integral. Tried functions: 1. f(r) -> problematic when f turns negative 2. f(r) ** 2 -> yields somewhat short R 3. f(r) * r -> problematic when f turns negative 4. (f(r) * r)**2 -> yields too short R 5. f(r)**2 * r**3 -> much better 6. abs(f(r)) -> yields a pretty long tail, but should be fine With the current approach users can manually control the behaviour by supplying a dictionary which passed the arguments to the radial_minimize_range method. Signed-off-by: Nick Papior <[email protected]>

dr argument controls the precision of the integrals. Signed-off-by: Nick Papior <[email protected]>

This ensures that copying does not destroy the R value. Added tests for all orbital classes to ensure it is copied over correctly. Signed-off-by: Nick Papior <[email protected]>

zerothi force-pushed the 574-orbital-range-rework branch from 4d2833b to 7f1aab0 Compare May 31, 2023 10:34

zerothi force-pushed the 574-orbital-range-rework branch 2 times, most recently from 7dc13ae to f89b533 Compare May 31, 2023 12:06

zerothi force-pushed the 574-orbital-range-rework branch from 8f5fe0c to 8ad7d80 Compare May 31, 2023 20:27

zerothi force-pushed the 574-orbital-range-rework branch 2 times, most recently from 07c5743 to 781abbd Compare June 6, 2023 10:41

zerothi added 4 commits June 7, 2023 09:24

amended typing and added dr argument

5f42b60

dr argument controls the precision of the integrals. Signed-off-by: Nick Papior <[email protected]>

Ensured orbital copy passes over the R value

37e073b

This ensures that copying does not destroy the R value. Added tests for all orbital classes to ensure it is copied over correctly. Signed-off-by: Nick Papior <[email protected]>

zerothi force-pushed the 574-orbital-range-rework branch from 781abbd to 37e073b Compare June 7, 2023 07:24

zerothi merged commit 92d6a87 into main Jun 7, 2023

zerothi deleted the 574-orbital-range-rework branch June 7, 2023 08:58

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

cleaned the procedure for locating maxR #575

cleaned the procedure for locating maxR #575

zerothi commented May 31, 2023

sofiasanz commented May 31, 2023

zerothi commented May 31, 2023

zerothi commented May 31, 2023

zerothi commented May 31, 2023

sofiasanz commented May 31, 2023

sofiasanz commented May 31, 2023

codecov bot commented May 31, 2023 •

edited

Loading

sofiasanz commented May 31, 2023

zerothi commented May 31, 2023 •

edited

Loading

zerothi commented May 31, 2023

tfrederiksen commented May 31, 2023

tfrederiksen commented Jun 1, 2023

zerothi commented Jun 1, 2023

tfrederiksen commented Jun 1, 2023 •

edited

Loading

zerothi commented Jun 1, 2023

zerothi commented Jun 1, 2023

sofiasanz commented Jun 1, 2023

zerothi commented Jun 1, 2023 •

edited

Loading

sofiasanz commented Jun 1, 2023

zerothi commented Jun 1, 2023

zerothi commented Jun 1, 2023

zerothi commented Jun 6, 2023

sofiasanz commented Jun 6, 2023 •

edited

Loading

sofiasanz commented Jun 6, 2023 •

edited

Loading

zerothi commented Jun 7, 2023

cleaned the procedure for locating maxR #575

cleaned the procedure for locating maxR #575

Conversation

zerothi commented May 31, 2023

sofiasanz commented May 31, 2023

zerothi commented May 31, 2023

zerothi commented May 31, 2023

zerothi commented May 31, 2023

sofiasanz commented May 31, 2023

sofiasanz commented May 31, 2023

codecov bot commented May 31, 2023 • edited Loading

Codecov Report

sofiasanz commented May 31, 2023

zerothi commented May 31, 2023 • edited Loading

zerothi commented May 31, 2023

tfrederiksen commented May 31, 2023

tfrederiksen commented Jun 1, 2023

zerothi commented Jun 1, 2023

tfrederiksen commented Jun 1, 2023 • edited Loading

zerothi commented Jun 1, 2023

zerothi commented Jun 1, 2023

sofiasanz commented Jun 1, 2023

zerothi commented Jun 1, 2023 • edited Loading

sofiasanz commented Jun 1, 2023

zerothi commented Jun 1, 2023

zerothi commented Jun 1, 2023

zerothi commented Jun 6, 2023

sofiasanz commented Jun 6, 2023 • edited Loading

sofiasanz commented Jun 6, 2023 • edited Loading

zerothi commented Jun 7, 2023

codecov bot commented May 31, 2023 •

edited

Loading

zerothi commented May 31, 2023 •

edited

Loading

tfrederiksen commented Jun 1, 2023 •

edited

Loading

zerothi commented Jun 1, 2023 •

edited

Loading

sofiasanz commented Jun 6, 2023 •

edited

Loading

sofiasanz commented Jun 6, 2023 •

edited

Loading