Allow orbitals with zero occupancy

[Kolaru]

First, thanks a lot for the package! It's great.

I've started using it, and one thing that would be useful for me is to allow the orbitals to be empty.

Typically, when considering x-ray physics, it's not impossible to have empty core orbitals. It is nice to keep track of them, because electrons from other orbitals can decay into the empty ones, and it simplifies the interfacing with external code that actually compute the electronic structure.

Then, probably, it would be nice to have few more utility function e.g. to get the highest occupied orbital, which could no longer be assumed to be config[end].

I can make a PR this week if it sounds like a good idea.

[jagot]

Thanks for your interest in AtomicLevels.jl!

Could you tell me a bit more about what you would like to achieve? I am not sure allowing zero occupancy in a Configuration is a good idea, because it might break things all over the place (also in dependent codes), in more or less subtle ways. Therefore, if you can tell me what you need, maybe we can figure out how to do it, without changing basic assumptions of the code.

For the description of Auger–Meitner decay following core-electron ionization of e.g. neon, I would probably at least include the following CSFs in the CI problem:

julia> csfs([c"1s2 2s2 2p6", c"1s2 2s2 2p5", c"1s 2s2 2p6"])
3-element Vector{NonRelativisticCSF{Orbital{Int64}, Seniority}}:
 1s²(₀¹S|¹S) 2s²(₀¹S|¹S) 2p⁶(₀¹S|¹S)+
 1s²(₀¹S|¹S) 2s²(₀¹S|¹S) 2p⁵(₁²Pᵒ|²Pᵒ)-
 1s(₁²S|²S) 2s²(₀¹S|²S) 2p⁶(₀¹S|²S)+

and if you want to also keep track of the continuum electron:

julia> csfs(Configuration{Orbital}[c"1s2 2s2 2p6", c"1s2 2s2 2p5 ks", c"1s2 2s2 2p5 kd", c"1s 2s2 2p6 kp"])
11-element Vector{NonRelativisticCSF{Orbital, Seniority}}:
 1s²(₀¹S|¹S) 2s²(₀¹S|¹S) 2p⁶(₀¹S|¹S)+
 1s²(₀¹S|¹S) 2s²(₀¹S|¹S) 2p⁵(₁²Pᵒ|²Pᵒ) ks(₁²S|¹Pᵒ)-
 1s²(₀¹S|¹S) 2s²(₀¹S|¹S) 2p⁵(₁²Pᵒ|²Pᵒ) ks(₁²S|³Pᵒ)-
 1s²(₀¹S|¹S) 2s²(₀¹S|¹S) 2p⁵(₁²Pᵒ|²Pᵒ) kd(₁²D|¹Pᵒ)-
 1s²(₀¹S|¹S) 2s²(₀¹S|¹S) 2p⁵(₁²Pᵒ|²Pᵒ) kd(₁²D|¹Dᵒ)-
 1s²(₀¹S|¹S) 2s²(₀¹S|¹S) 2p⁵(₁²Pᵒ|²Pᵒ) kd(₁²D|¹Fᵒ)-
 1s²(₀¹S|¹S) 2s²(₀¹S|¹S) 2p⁵(₁²Pᵒ|²Pᵒ) kd(₁²D|³Pᵒ)-
 1s²(₀¹S|¹S) 2s²(₀¹S|¹S) 2p⁵(₁²Pᵒ|²Pᵒ) kd(₁²D|³Dᵒ)-
 1s²(₀¹S|¹S) 2s²(₀¹S|¹S) 2p⁵(₁²Pᵒ|²Pᵒ) kd(₁²D|³Fᵒ)-
 1s(₁²S|²S) 2s²(₀¹S|²S) 2p⁶(₀¹S|²S) kp(₁²Pᵒ|¹Pᵒ)-
 1s(₁²S|²S) 2s²(₀¹S|²S) 2p⁶(₀¹S|²S) kp(₁²Pᵒ|³Pᵒ)-

[Note that you then have to add Configuration{Orbital} in front of the vector, because k is stored as a "principal" quantum number of type Symbol of the orbital; something I now consider to be a design flaw. Since the principal quantum number is not a group irrep label, except for the very special case of non-relativistic hydrogen, which obeys SO(4), it should not be kept track of. However, it is still useful to have an ordinal quantum number for user input, so you win some you lose some.]

If you instead work with Slater determinants, you can get those as

julia> spin_configurations([c"1s2 2s2 2p6", c"1s2 2s2 2p5 ks", c"1s 2s2 2p6"])
15-element Vector{SpinConfiguration{SpinOrbital{<:Orbital, Tuple{Int64, HalfIntegers.Half{Int64}}}}}:
 1s₀α 1s₀β 2s₀α 2s₀β 2p₋₁α 2p₋₁β 2p₀α 2p₀β 2p₁α 2p₁β
 1s₀α 1s₀β 2s₀α 2s₀β 2p₋₁α 2p₋₁β 2p₀α 2p₀β 2p₁α ks₀α
 1s₀α 1s₀β 2s₀α 2s₀β 2p₋₁α 2p₋₁β 2p₀α 2p₀β 2p₁α ks₀β
 1s₀α 1s₀β 2s₀α 2s₀β 2p₋₁α 2p₋₁β 2p₀α 2p₀β 2p₁β ks₀α
 1s₀α 1s₀β 2s₀α 2s₀β 2p₋₁α 2p₋₁β 2p₀α 2p₀β 2p₁β ks₀β
 1s₀α 1s₀β 2s₀α 2s₀β 2p₋₁α 2p₋₁β 2p₀α 2p₁α 2p₁β ks₀α
 1s₀α 1s₀β 2s₀α 2s₀β 2p₋₁α 2p₋₁β 2p₀α 2p₁α 2p₁β ks₀β
 1s₀α 1s₀β 2s₀α 2s₀β 2p₋₁α 2p₋₁β 2p₀β 2p₁α 2p₁β ks₀α
 1s₀α 1s₀β 2s₀α 2s₀β 2p₋₁α 2p₋₁β 2p₀β 2p₁α 2p₁β ks₀β
 1s₀α 1s₀β 2s₀α 2s₀β 2p₋₁α 2p₀α 2p₀β 2p₁α 2p₁β ks₀α
 1s₀α 1s₀β 2s₀α 2s₀β 2p₋₁α 2p₀α 2p₀β 2p₁α 2p₁β ks₀β
 1s₀α 1s₀β 2s₀α 2s₀β 2p₋₁β 2p₀α 2p₀β 2p₁α 2p₁β ks₀α
 1s₀α 1s₀β 2s₀α 2s₀β 2p₋₁β 2p₀α 2p₀β 2p₁α 2p₁β ks₀β
 1s₀α 2s₀α 2s₀β 2p₋₁α 2p₋₁β 2p₀α 2p₀β 2p₁α 2p₁β
 1s₀β 2s₀α 2s₀β 2p₋₁α 2p₋₁β 2p₀α 2p₀β 2p₁α 2p₁β

When you have the CSFs/Slater determinants, it is a little bit up to you what you do with them. Typically, I would collect all unique orbitals into a vector, thereby establishing a canonical ordering, and then create a matrix of integers (for CSFs) or booleans (for Slater determinants), where the element O[i,j] corresponds to the occupancy of the ith orbital in jth configuration. Note that for CSFs, you would also need to keep track of the coupling tree. For Slater determinants, I use something similar to https://arxiv.org/abs/1311.6244v1

[Kolaru]

Thanks a lot for your answer!

Currently, I'm doing something quite simple: I am using AtomicLevels.jl to keep track of the electronic configuration of several atoms. I calculate rates of atomic processes with an external software (XATOM), and I include them in a Monte-Carlo scheme to simulate cascades of processes following one or multiples core ionizations (I am not myself that deep into the development of electronic structure calculations).

The issue is that if, for example, I photoionize the 1s orbital twice (I'm considering very intense x-ray pulses), it just disappears:

julia> c"1s2 2s2 2p6" - o"1s" - o"1s"
2s² 2p⁶

I believe (I only weakly believe it as I am a user of electronic structure code, but not an expert ^^') that you cannot decay to an orbital that doesn't explicitly exist.

But really, for me, it is a convenience thing: having explicit empty orbitals makes it easier to track what configurations are available in principle (I just move electrons around in the orbitals given), as well as making it easier to interface with other codes.

As a side note, I have to admit that I am not familiar with the the k-notation that you are using for continuum electrons (eg c"1s2 2s2 2p5 ks"), and I don't know if it could be useful to me (it would not for the purpose of continuum electrons, as I'm modelling them as point particle in the Monte-Carlo scheme -- I'm doing quite crude approximations, but it has be working okay so far :D).

[jagot]

I think I understand; you want to generate input to XATOM (say hello to Sang-Kil Son from me!), or maybe parse its output. I have never used this code myself, so I am not familiar with exactly what it requires, but I can roughly guess.

If you already know this, I apologize, but it can perhaps still be useful:

It is important to understand that an orbital does not exist as a physical entity. It is merely a convenient mathematical object, that gives a good description of the state of one electron in a many-body wavefunction. Formally, an orbital is a vector in the vector space of possible states of one electron. The many-body wavefunction is formed as an antisymmetrized outer product of such orbitals, this is commonly referred to as a Slater determinant. Even better, the many-body wavefunction is expanded a linear combination of such Slater determinants, and this is referred to as multiconfigurational. A configuration is a short-hand notation for which orbitals are occupied in a given wavefunction.

In Auger–Meitner decay (and similar processes such as autoionization), the system is prepared in a highly excited state, which through correlation (typically the Coulomb interaction between electrons) will after some time decay to a lower-lying state. The highly excited state is typically characterized by a inner-shell vacancy, or equivalently a hole in a lower-lying orbital, whereas the lower-lying state which the initial state decays to, this vacancy is filled. In a pictographic language, we say that the electron "jumps" from one orbital to another, to fill the vacancy. The probability for this to happen is proportional to the absolute square of the Hamiltonian matrix element between the states: $|\langle f | \hat{H} | i \rangle|^2$. In that sense, the decay can absolutely happen, even if the "orbital does not exist" in $|i\rangle$; it is not occupied in $|i\rangle$, but it is in $|f\rangle$.

$|k \ell\rangle$ is the typical notation for a free electron with linear momentum magnitude $k$ and orbital angular momentum $\ell$; it is similar to $|n \ell\rangle$ for bound electrons, where the principal quantum number $n$ is related to the number of nodes of the radial wavefunction ($n-\ell-1$ nodes).

A word of warning though: e.g. the 1s orbital is not necessarily the same in the ground state of neon 1s^2 2s^2 2p^6 and in singly ionized neon 1s 2s^2 2p^6; the remaining electrons will have adjusted when one electron is missing.

If you let me know some more details of what you need to feed into XATOM, or extract from it, in terms of configurations from AtomicLevels.jl, I will be happy to help. However, changing Configuration to hold empty orbitals is not an option, I am afraid, since that can probably (most likely will) break other things.

[Kolaru]

Please, assume that I know nothing :D

More seriously, thanks a lot for the explanation. The key part that I missed is indeed that you need the full initial and final state: the probability of the process is not computed from the initial state alone.

Now regarding what I need for XATOM: I only need a string describing the configuration itself, that I generate from a Configuration. I am really using AtomicLevels as only book keeping device for now, which may not use its full potential, but which is extremely convenient.

I am convince that zero-occupancy is indeed neither necessary, nor desirable.

Thanks a lot for your time!

PS: You are likely to hear from me again soon, as I am extending the package on my side with small utilities that I'd like to share e.g. generating a configuration directly from a PeriodicTables elements or indexing a Configuration with an Orbital.
PS2: I'll say hi to Sang-Kil for you ;)

[jagot]

Be sure to check out https://github.com/JuliaAtoms/AtomicStructure.jl/blob/master/src/potentials/table_of_elements.jl; maybe this does what you want? In that case, we can probably move this functionality from AtomicStructure.jl to AtomicLevels.jl.

[Kolaru]

I am a simple person, I was really just thinking to propose the following as a package extension for PeriodicTables.jl:

function Configuration(element::Element)
    return parse(Configuration{Orbital}, element.el_config)
end

That would allow to do the following:

julia> using PeriodicTable

julia> using AtomicLevels

julia> Configuration(elements[:I])
1s² 2s² 2p⁶ 3s² 3p⁶ 4s² 3d¹⁰ 4p⁶ 5s² 4d¹⁰ 5p⁵