I think it is a really good idea, however it will introduce some corner case problems:

• Geometry this class should be trivially handled
• SparseGeometry (Hamiltonians etc) how should geometries be handled here? For instance Geometry.close will return the non orbital atoms, and this requires extra steps to leave them out.
  Files with 0 orbitals will not work with all file-formats, for instance, Siesta TSHS files requires orbitals on all atoms.

If one takes an atom that has 0 orbitals then H[ia, ...] will be the atom before it. So what should actually happen?

I think it is a really good idea, however it will introduce some corner case problems:

* `Geometry` this class should be trivially handled

* `SparseGeometry` (Hamiltonians etc) how should geometries be handled here? For instance `Geometry.close` will return the _non_ orbital atoms, and this requires extra steps to leave them out.
  Files with 0 orbitals will not work with all file-formats, for instance, Siesta TSHS files requires orbitals on all atoms.

If one takes an atom that has 0 orbitals then H[ia, ...] will be the atom before it. So what should actually happen?

I was worried about the same thing. However, I would say that it would be more logical that a close call on a hamiltonian should return the indices of orbitals (not atoms, which are actually irrelevant). This would solve the issue, no?

Partially, since lots of code is relying on Geometry close calls.

I am also worried about the documentation issue, i.e. which routines uses atoms with orbitals and which are using all atoms. I would intuitively always expect a geometry to behave the same, but that is not what will happen if Hamiltonian removes them (it also requires additional boiler plate code).

What I like about using class objects is that it is very clear how they interact with each other. And in this case it will change the behaviour of Geometry in certain situations. My experience tells me that some may find this intuitive, but new users will not have this intuition... As such I am probably more in favour of not allowing this, because it will break backwards compatibility. Many are using geometry.close to populate Hamiltonians as I have this in the tutorials ;)

Many are using geometry.close to populate Hamiltonians as I have this in the tutorials ;)

Yeah, I understand that there is an issue with backwards compatibility. But to me the current layout is not ideal.

I have been playing with setting up some some multi-orbital TB models and then one ends up with loops like this:

for ia in geom:
    io = geom.a2o(ia, all=True)
    for ja in geom.close(ia, R=(...)):
        jo = geom.a2o(ja, all=True)
        for j in ja:
            H[io, jo] = ...

It starts to look a little ugly and to me it would be nicer (and, and actually, more intuitive I think) to have something like:

for io in H:
    for jo in H.close(io, R=(...)):
        H[io, jo] = ...

It would just boil down to changing the meaning of H.close, no? For the tutorials with one-orbital-per-atom there should be no problem by looping over the geometry itself. But one could also have done it by iterating over the hamiltonian.

I agree, this is why I generally recommend one to create a construct function and pass that, it is a bit easier, but still loong.

The problem of the for io in H is the slow looping since you'll effectively calculate the same thing over and over again (i.e. atoms close to one atom). Also, that is a somewhat different issue that can added without problems.

I agree your suggestion would make sense, but it would require a stateful close function, here what @pfebrer does in #393 should partially solve this. However, for compatibility with other sparse matrix classes an iterator on H will yield row, col values that are populated. So this will be hard to change. This is how it is currently being used:

for row, col in H:

I don't think the above should be changed since that would be very different from other sparse matrix handlers.

So while one can fix the Hamiltonian.close it still has many problems related to files that contain orbital information. I.e. reading/writing to files TSHS, nc files, if one has 0-orbital atoms, they are not parseable by siesta/tbtrans but people will request this... :(

OK, I see. Another approach then could then be to allow Geometry to hold some kind of external atoms, such as atoms without orbitals? They could then be included in visualizations and geom.close/geom.write calls etc. with some optional argument external_atoms=True.

Hmm.. Ok, what about allowing 0 orbitals, and when creating SparseGeometries with geometries having 0 orbitals an error will be raised? I can't immediately see all implications, but at least this will be annoying enough that people will see it! And it will force users to change the geometry before creating Hamiltonians etc?

Hmm.. Ok, what about allowing 0 orbitals, and when creating SparseGeometries with geometries having 0 orbitals an error will be raised? I can't immediately see all implications, but at least this will be annoying enough that people will see it! And it will force users to change the geometry before creating Hamiltonians etc?

Sounds interesting. Indeed I think it would be intuitive to allow for 0 orbitals. But I do not see how that would help to tackle the situation where one would like to work with a physical geometry and some effective Hamiltonian for a subset of the atoms.

I think there wouldn't be too much of a problem, I am just very worried that the 0-orbital atom will cause havoc in the SparseGeometry classes.

The workflow for 0-orbital geometries would be:

1. create geometry
2. remove 0-orbital atoms
3. pass to SparseGeometry
4. Plotting functionality should then plot from the SparseGeometry + the geometry with only 0-orbitals

Not pretty, but I think this will not break anything...
I still cannot see immediately how easy they are to manage behind the scenes. If you want to try it out in a simple PR, it would be great to gain some experience.

What if such a geometry passed to SparseGeometry is internally filtered to move zero-orbital atoms into a separate placeholder? Then it could be simple flag if they should be included or not in various calls (visualizations, file io, close calls etc).

What if such a geometry passed to SparseGeometry is internally filtered to move zero-orbital atoms into a separate placeholder? Then it could be simple flag if they should be included or not in various calls (visualizations, file io, close calls etc).

I have done this with the Geometry._names and it turned out to be horrible to maintain. It requires changes all around which I don't think is justifying its use. Basically all calls which alters a SparseGeometry needs some way of defining how it should handle these atoms.

Also, a class should do one thing, and it should be clear to the end user. I think this is crossing the "do one thing" without providing much aid for the end user. I really try to make things simple in sisl, and I think this is a borderline case where things are being needlessly more complicated without gain (especially for the SparseGeometry since they intrinsically deals with orbitals). And there are simple ways around it by having a geometry with the atoms having 0-orbitals only. That can easily be plotted on top of another geometry.

It is not that I don't understand how it can aid in these plotting functionalities, but... I would still think that having a separate geometry that gets plotted as a subsequent step is not too much of a deal compared to the internal restructure and unsureness of how a Geometry interacts with SparseGeometry and when a 0-orbital atom is used vs. when it is not used. I am really worried that it would be really useful for TB implementations, but that it will be confusing to actually use.

I would also like this feature for QM/MM setups :)

Perhaps a first step could be to allow atoms to have 0 orbitals but raise an error when trying to create a SparseOrbital with a geometry containing one?

For SparseAtom there is no problem, right?

Hmm...

How should a default geometry be created? Currently a base of 1 orbital per atom is the default.

Lets get #665 in first (restructuring). I have just returned from vacation. So I should be more available now.

How should a default geometry be created? Currently a base of 1 orbital per atom is the default.

Probably keeping the current default is the only sane option :)

I would also like this feature for QM/MM setups :)

Perhaps a first step could be to allow atoms to have 0 orbitals but raise an error when trying to create a SparseOrbital with a geometry containing one?

For SparseAtom there is no problem, right?

I think you should start here @pfebrer , just die when creating matrices with 0-orbital atoms. :)

Is there something to discuss really? Or should I just try and see what happens?

One option would be to put all 0-orbital atoms at the end or at the beggining, but I don't think it is a good option.

Another option would be to keep a mask with len(orbitals) > 0, or an array with the indices of the orbital-free atoms.

Also could we try to merge #496 and #765 first? Otherwise my head will explode with so many different branches to think about 😅

Is there something to discuss really? Or should I just try and see what happens?
I am sure something will pop up ;)
Lets discuss those details here, instead of in the PR.

One option would be to put all 0-orbital atoms at the end or at the beggining, but I don't think it is a good option.
Agreed, we shouldn't try to force these things.

Another option would be to keep a mask with len(orbitals) > 0, or an array with the indices of the orbital-free atoms.

Lets raise an error first on sparse matrices, then see what needs fixing.

Also could we try to merge #496 and #765 first? Otherwise my head will explode with so many different branches to think about 😅

Definitely #765 , I have to wrap my head around #496 and figure out how it all fits together!