Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Atoms without orbitals? #449

Open
tfrederiksen opened this issue Mar 5, 2022 · 19 comments · May be fixed by #733
Open

Atoms without orbitals? #449

tfrederiksen opened this issue Mar 5, 2022 · 19 comments · May be fixed by #733

Comments

@tfrederiksen
Copy link
Contributor

tfrederiksen commented Mar 5, 2022

I am trying to build a geometry for which I would like to have zero orbitals associated to certain atoms. Can it be done?
As far as I can see, I get at least one orbital associated to an atom:

>>> print(sisl.Atom('H', orbitals=0))
Atom{H, Z: 1, mass(au): 1.00794, maxR: 0.00000,
 Orbital{R: 0.00000, q0: 0.0}
}
>>> print(sisl.Atom('H', orbitals=(0, 1)))
Atom{H, Z: 1, mass(au): 1.00794, maxR: 1.00000,
 Orbital{R: 0.00000, q0: 0.0},
 Orbital{R: 1.00000, q0: 0.0}
}
>>> print(sisl.Atom('H', orbitals=[]))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "~/.local/lib/python3.8/site-packages/sisl/atom.py", line 1013, in __init__
    raise ValueError(f"{self.__class__.__name__}.__init__ got unparseable 'orbitals' argument: {orbitals}")
ValueError: Atom.__init__ got unparseable 'orbitals' argument: []

Motivation: Think of a GNR structure with H-atoms at the edges (it could be a geometry obtained from DFT). Now, if we could simply set zero orbitals for the H and one orbital for C, we would readily arrive at the effective model for the pi-electrions only, without having to remove any physical atoms from the geometry.

@zerothi
Copy link
Owner

zerothi commented Mar 7, 2022

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?

@tfrederiksen
Copy link
Contributor Author

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?

@zerothi
Copy link
Owner

zerothi commented Mar 7, 2022

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 ;)

@tfrederiksen
Copy link
Contributor Author

tfrederiksen commented Mar 7, 2022

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.

@zerothi
Copy link
Owner

zerothi commented Mar 7, 2022

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... :(

@tfrederiksen
Copy link
Contributor Author

tfrederiksen commented Mar 7, 2022

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.

@zerothi
Copy link
Owner

zerothi commented Mar 7, 2022

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?

@tfrederiksen
Copy link
Contributor Author

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.

@zerothi
Copy link
Owner

zerothi commented Mar 7, 2022

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.

@tfrederiksen
Copy link
Contributor Author

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).

@zerothi
Copy link
Owner

zerothi commented Mar 7, 2022

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.

@pfebrer
Copy link
Contributor

pfebrer commented Feb 13, 2024

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?

@zerothi
Copy link
Owner

zerothi commented Feb 13, 2024

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.

@pfebrer
Copy link
Contributor

pfebrer commented Feb 13, 2024

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 :)

@zerothi zerothi linked a pull request Mar 26, 2024 that will close this issue
5 tasks
@zerothi
Copy link
Owner

zerothi commented May 2, 2024

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. :)

@pfebrer
Copy link
Contributor

pfebrer commented May 2, 2024

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.

@pfebrer
Copy link
Contributor

pfebrer commented May 2, 2024

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

@zerothi
Copy link
Owner

zerothi commented May 2, 2024

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.

@zerothi
Copy link
Owner

zerothi commented May 2, 2024

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!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants