sisl icon indicating copy to clipboard operation
sisl copied to clipboard

EigenstateElectron.wavefunction fails in some cases

Open pfebrer opened this issue 2 years ago • 19 comments

EigenstateElectron.wavefunction fails at self.change_gauge("R") if:

  • The parent is a geometry
  • Gauge is "r"
  • K is different than gamma

because change_gauge doesn't contemplate that self.parent might be a geometry.

Here's some code to reproduce the error:

import sisl
gr = sisl.geom.graphene()

eigenstate = sisl.EigenstateElectron([[0, 1]], [0], parent=gr, gauge="r", k=[0.5, 0, 0])

grid = sisl.Grid(10, geometry=gr, dtype=complex)
eigenstate[0].wavefunction(grid)

pfebrer avatar Nov 16 '21 12:11 pfebrer

Yes, these corner cases are abundant. My best idea is to force parent to contain the geometry in some on the fly class.

There are simply too many places where it is assumed that parent is some kind of sparse matrix. Even if you fix this in the change_gauge then you'll see the same error in wavefunction. So ;(

zerothi avatar Nov 16 '21 13:11 zerothi

Yes, these corner cases are abundant. My best idea is to force parent to contain the geometry in some on the fly class.

Yes, I would say that's a good idea. Or a method like:

class EigenstateElectron:
    ...
    def get_geometry(self):
        if isinstance(self.parent, Geometry):
            geometry = self.parent
        else:
            geometry = getattr(self.parent, "geometry", None)
       
        return geometry

? Which could even be a property.

Then you wouldn't have to write it everywhere? Just replace *.parent.geometry for *.get_geometry() whenever you want the geometry.

pfebrer avatar Nov 16 '21 13:11 pfebrer

yes, but that will only fix it for State cases... This is where class based stuff is really bad since it really forces stricter rules and limits the use-cases.

In the above, one might as well do:

class State(...)

@property
def geometry(self):
    if isinstance(self.parent, Geometry):
        return self.parent
     return self.parent.geometry

in any case it would be good to explicitly fail (to notify users of bad behaviour), no?

zerothi avatar Nov 16 '21 13:11 zerothi

class GeometryGetter:
    
    def __get__(self, instance, owner):
        if instance is not None and hasattr(instance, 'parent'):
            if isinstance(instance.parent, Geometry):
                return instance.parent
            return instance.parent.geometry

and add it as a class attribute to any class that has a parent?

class ClassWithParent:
    geometry = GeometryGetter()

? You could change "parent" by any other key on initialization I guess too.

pfebrer avatar Nov 16 '21 13:11 pfebrer

in any case it would be good to explicitly fail (to notify users of bad behaviour), no?

Yes, but I would prefer it to be fixed at least in this case, because otherwise you can not plot wavefunctions at k points other than gamma from the WFSX file, and I have some huge systems for which calculating the eigenstates from the hamiltonian in sisl is a really bad idea :(

pfebrer avatar Nov 16 '21 13:11 pfebrer

in any case it would be good to explicitly fail (to notify users of bad behaviour), no?

Yes, but I would prefer it to be fixed at least in this case, because otherwise you can not plot wavefunctions at k points other than gamma from the WFSX file, and I have some huge systems for which calculating the eigenstates from the hamiltonian in sisl is a really bad idea :(

So you say you do not want to raise an error? I don't understand why by what you are saying?

zerothi avatar Nov 16 '21 13:11 zerothi

No, no. Of course if you need parent to be a sparse matrix an error should be raised. But in this case you don't need it, right? You just need the geometry. Or maybe I'm missing something that you need to plot the wavefunction :sweat_smile:

pfebrer avatar Nov 16 '21 13:11 pfebrer

You need the geometry to be able to figure out where the orbitals are, but if you have the geometry elsewhere, you can pass it as a direct argument?

zerothi avatar Nov 16 '21 13:11 zerothi

Hmm not to EigenstateElectron.wavefunction. And you can't call sisl.physics.electron.wavefunction because you need to change the gauge to "R". You also can not pass a geometry to change_gauge so you'd have to copy all the code in EigenstateElectron.wavefunction plus change_gauge passing your geometry instead of self.parent.geometry.

pfebrer avatar Nov 16 '21 13:11 pfebrer

Then just call it your-self :)

        wavefunction(self.state, grid, geometry=geometry, k=k, spinor=spinor, spin=spin, eta=eta)

actually I would like to move away from some of these methods being part of the API because of exactly the reason you have. With many more physical quantities one can calculate, the objects becomes very large and sometimes non-deterministic. Say the state.velocity where degenerate states are calculated.

In fact, I would rather have that velocity has an interface, velocity(state, ...) and then users can adapt their own scripts...

Calling state.DOS(...) or DOS(state, ...) shouldn't matter. And functional programming has other advantages in these aspects, readability for end-users for one. :)

zerothi avatar Nov 16 '21 13:11 zerothi

Then just call it your-self :)

But I can't, because I have to change the gauge. And I can't simply call state.change_gauge() because this looks for self.parent.geometry. So I would have to calculate the phase and do the conversion in user code.

pfebrer avatar Nov 16 '21 15:11 pfebrer

Hmm well I'm not sure which approach is better (object oriented vs functional), but to me it is not that clear that functional is better, specially for interactive programming.

If you are doing some analysis in a jupyter notebook I see more advantages than disadvantages in having methods. It is much easier to find what you can do with a given object and explore the possibilities. With the other approach, you have to import everything that you want to use.

It may be good for the readability of scripts or library code, but not good for a dynamic environment like IPython. Plus, all that improvement in readability is lost when people do things like:

from a import *
from b import *
from c import *
from d import *

some_function_from_god_knows_where()

which I've seen is done all across SIESTA, for example.

Also, classes can be used as modules if you wish, you don't need to call methods like obj.method(). These two calls are equivalent:

class A:
    def some_method(self, ...):
        ...

a = A()

a.some_method()
A.some_method(a)

You can not do from module.A import some_method, but that may even be a good thing if you are discussing readability.

And if you want to do it purely functional, as in "the result of the function depends only on explicit arguments", then I don't think that's very suitable for interactive programming.

To me, what you have done with the functions in sisl.physics.electron is perfect. That is, the functions that actually do the job are much more declarative, which is very flexible and clear. But still you build a wrapper (in the form of a method) that hides some details in favour of a faster usage and less expertise requirements from the user side. The method allows you to play, while the function underneath allows you to tweak as much as you want.

pfebrer avatar Nov 16 '21 15:11 pfebrer

I actually agree on all points... :) My main worry is that classes become so complex with too many routines. This is the feedback I get on say the geometry class which have ridiculously many methods...

zerothi avatar Nov 16 '21 18:11 zerothi

I'm curious, who gave you that feedback? Was it from users? What is the problem of having many methods? Is it maintainability, usage or some other thing?

pfebrer avatar Nov 16 '21 22:11 pfebrer

I'm curious, who gave you that feedback? Was it from users? What is the problem of having many methods? Is it maintainability, usage or some other thing?

Yes, users and experienced users. The problem is probably more with using the program. One can do soo many things, and it is incomprehensible with everything. Some of this could simply be solved by grouping documentation https://github.com/zerothi/sisl/issues/203 and improving it in general.

zerothi avatar Nov 17 '21 07:11 zerothi

Hmm there are a lot of methods to convert indices which maybe could be in a single method called indices? It could work much like a unit converter:

geometry.indices([1,2,3], "atom", "sc")

Then there are "operations", which modify the coordinates of the system. They could be under an op attribute:

geometry.op.tile(2, 1)

although that is bad for chaining operations in the same line :(

And then there are "measurements" (e.g. distances, angles, center, finding neighbors...).

This splitting would be great for a GUI as well (I'm always thinking about graphical interfaces hehe).

pfebrer avatar Nov 17 '21 10:11 pfebrer

Yes, one way would be to group them like you suggest. But as you say this has the problem of verbosity.

After your #393 suggestion I have thought about something like this:

geometry = ...

neighbours = GeometryNeighbours(geometry, **flags)

neighbours.close(...)

The same idea could be applied to measurements etc. However, I can't understand whether this is just a convenience for documentation, or will actually help end-users...

zerothi avatar Nov 17 '21 10:11 zerothi

although that is bad for chaining operations in the same line :(

You could always add a __getattr__ to Geometry to get the methods from the operations, measurements and indices managers. That would keep backwards compatibility and still allow easy usage while making the documentation easier to read.

pfebrer avatar Nov 17 '21 10:11 pfebrer

making the documentation easier to read.

And easier for the GUI to present methods to the user of course :)

pfebrer avatar Nov 17 '21 10:11 pfebrer