electronic-structure / SIRIUS

Domain specific library for electronic structure calculations
BSD 3-Clause "New" or "Revised" License
126 stars 40 forks source link

Hubbard orbitals #578

Closed toxa81 closed 1 year ago

toxa81 commented 3 years ago

This is a placeholder to discuss the indexing of Hubbard orbitals. In QE Hubbard orbitals can be generated either from pseudo-atomic radial wave-functions (stored in the UPF) or from the Wannier functions.

Atomic orbitals

UPF files can be of two kinds: "normal" and "relativistic". In the first case pseudo-atomic orbitals are labeled by orbital quantum number ℓ={s,p,d,f}. In the second case pseudo-atomic wave functions should be 2-component spinors with "major" and "minor" components (in reality they are not, only one components is stored. which? "major"?) and are labeled by the total orbital quantum number j. For example, for d-shell one will find j=3/2 and j=5/2 radial wave-functions. Later, the QE code will create an average radial wave-function out of two (supposedly "major") 3/2 and 5/2 components. So basically, in SO case the code still deals with scalar {s,p,d,f}-like orbitals.

For non-magnetic calculation the Hubbard orbital basis is just |nℓm> localized pseudo-atomic functions. For magnetic collinear it is basically the same, just transformed into pure spinors (|nℓm>, 0) and (0, |nℓm>) but we can think of them as of scalar functions without a spin index. The non-collinear case is not clear. Are localized functions the same as in collinear case?

Wannier orbitals Wannier orbitals are created from the unitary transformation of Bloch KS wave-functions. They look similar to the atomic orbitals (let's think of them as 2ℓ+1 orthogonal functions). In magnetic collinear case the wave-functions have spin-up and spin-down channels and correspondingly, Wannier functions become spin-dependent. In the non-collinear case the KS wave-functions are two-component spinors. The Wannier functions, in principle, should also be two-component spinors. The number of such functions should be (2ℓ+1)

To summarize: in the case of atomic orbitals the number of independent Hubbard orbitals is 2ℓ+1 because, basically, this the only information about pseudo-atomic wave-functions, obtained from the UPF files. In case of Wannier orbitals the number of independent Hubbard orbitals is N_spin * (2ℓ+1) and in case of non-collinear magnetism each of the 2(2ℓ+1) Wannier functions has two components.

When we index Hubbard orbitals, we can stay in "spin-up", "spin-down" notation only if Hubbard orbitals are the pure spinors of (|phi>, 0) (0, |phi>) kind. Is it the case in current LDA+U implementation of QE?

toxa81 commented 3 years ago

For SIRIUS we should change the way Hubbard orbitals are loaded. We either set the atomic orbital for a particular ℓ-channel (allow multiple ℓ-channels) or we set plane-wave coefficients of the Wannier orbital. In the former case we give either a scalar orbital and provide ℓ or spinor orbital and provide j. The m-index will be expanded later by SIRIUS. In the later case we should set ℓ and then provide 2ℓ+1 sets of plane-wave coefficients for each k-point.

The averaging of the Hubbard radial functions should happen on the QE side.

Possible API:

! assume that atomic orbitals psi is already prepared by QE
! loop over atom types
do iat=1,nat
  ! loop over Hubbard orbitals of a given atom type
  do iorb=1,norb(iat)
    ! get orbital quantum number
    l = l_of_orbital(iorb, iat)
    ! get principal quantum number
    n = n_of_orbital(iorb, iat)
    ! get occupancy
    occ = o_of_orbital(iorb, iat)
    ! add Hubbard orbital and specify radial atomic part
    call sirius_add_hubbard_orbital(sim_ctx, label(iat), n, l, occ, U(iorb, iat), J(irob, iat), f=psi(:, iorb, iat))
  enddo
enddo

! for Wannier orbitals
call sirius_add_hubbard_orbital(sim_ctx, 'Ni', n=3, l=2, occupancy, U, J)
! and then later for k-points
do ik=1,nk
  i = 0
  do ia=1,nat
    do orb=1,norb
      l = l_of_psi(orb, ia)
      do m = -l, l
        i = i + 1
        call sirius_add_hubbard_orbital_k(ks_handler, ik, ia, l, 
enddo 
mtaillefumier commented 3 years ago

This is a placeholder to discuss the indexing of Hubbard orbitals. In QE Hubbard orbitals can be generated either from pseudo-atomic radial wave-functions (stored in the UPF) or from the Wannier functions.

Atomic orbitals

UPF files can be of two kinds: "normal" and "relativistic". In the first case pseudo-atomic orbitals are labeled by orbital quantum number ℓ={s,p,d,f}. In the second case pseudo-atomic wave functions should be 2-component spinors with "major" and "minor" components (in reality they are not, only one components is stored. which? "major"?) and are labeled by the total orbital quantum number j. For example, for d-shell one will find j=3/2 and j=5/2 radial wave-functions. Later, the QE code will create an average radial wave-function out of two (supposedly "major") 3/2 and 5/2 components. So basically, in SO case the code still deals with scalar {s,p,d,f}-like orbitals.

I reopened a book of quantum mechanics (unfortunately not Landau Lipschitz because I do not where the book is) to check how the orbital should look like in presence of SO. So according to the theory the radial part of the wave functions depend only n and j not j_z. So your description of the pseudo-wave functions is reasonable. The pseudo potential should have two radial wave functions for each (n,l) one major l + 1/2 and one minor for l-1/2. QE and sirius just create an average of these two orbitals when we use hubbard in the relativistic case.

For non-magnetic calculation the Hubbard orbital basis is just |nℓm> localized pseudo-atomic functions. For magnetic collinear it is basically the same, just transformed into pure spinors (|nℓm>, 0) and (0, |nℓm>) but we can think of them as of scalar functions without a spin index. The non-collinear case is not clear. Are localized functions the same as in collinear case?

Indeed they are the same (colinear, non-colinear) for pseudo that are not relativistic because we do not have the information about the level splitting induced by e-e interactions and/or spin-orbit coupling.

The non-colinear case is tricky when relativistic pseudo are used. Hubbard should a priori be different for the groups l+1/2 and l-1/2 but QE and sirius do not make the distinction.

Wannier orbitals Wannier orbitals are created from the unitary transformation of Bloch KS wave-functions. They look similar to the atomic orbitals (let's think of them as 2ℓ+1 orthogonal functions). In magnetic collinear case the wave-functions have spin-up and spin-down channels and correspondingly, Wannier functions become spin-dependent. In the non-collinear case the KS wave-functions are two-component spinors. The Wannier functions, in principle, should also be two-component spinors. The number of such functions should be (2ℓ+1)

To summarize: in the case of atomic orbitals the number of independent Hubbard orbitals is 2ℓ+1 because, basically, this the only information about pseudo-atomic wave-functions, obtained from the UPF files. In case of Wannier orbitals the number of independent Hubbard orbitals is N_spin * (2ℓ+1) and in case of non-collinear magnetism each of the 2(2ℓ+1) Wannier functions has two components.

When we index Hubbard orbitals, we can stay in "spin-up", "spin-down" notation only if Hubbard orbitals are the pure spinors of (|phi>, 0) (0, |phi>) kind. Is it the case in current LDA+U implementation of QE?

The only problem comes from using relativistic pseudo with hubbard.

timrov commented 3 years ago

Hi guys, A few points:

When we index Hubbard orbitals, we can stay in "spin-up", "spin-down" notation only if Hubbard orbitals are the pure spinors of (|phi>, 0) (0, |phi>) kind. Is it the case in current LDA+U implementation of QE?

As I wrote above, so far we have been dealing with the 'atomic' (and 'ortho-atomic') orbitals, which do not depend on the spin index. I haven't thought about the generalization to the case when they depend on the spin index for the relativistic case. But this is something that we want to do. One PhD student wants to work on the generalization to the non-collinear case.

A general question about SIRIUS. As far as I understand, the logic of SIRIUS does not follow exactly the same logic of QE, correct? So basically you reimplement everything in a different way, correct? So SIRIUS is just a different code, not a GPU-enabled mirror of QE, right?

mtaillefumier commented 3 years ago

Hi Iurii,

On 10/27/20 11:43 AM, Iurii Timrov wrote:

Hi guys, A few points:

  • Wannier functions are not currently supported. Recently I have started working on this. I made progress on the spin-unpolarized case. Now I am working on the collinear spin-polarized case, and indeed I started thinking that now WFs depend on the spin index and so this is different from the 'atomic' case which does not depend on spin indices.
  • So far I have been working on DFT+U for the non-spin-polarized and collinear spin-polarized cases. So I do not know the answers to your questions about the relativistic PPs, spin-orbit coupling, etc. etc.

    When we index Hubbard orbitals, we can stay in "spin-up", "spin-down" notation only if Hubbard orbitals are the pure spinors of (|phi>, 0) (0, |phi>) kind. Is it the case in current LDA+U implementation of QE?

In that case, we can only store the up or down component of the spinor but it only works (from a theory standpoint) when the pps are non relativistic. QE uses the full spinor for non-colinear magnetism in general. I do not remember what implementation they chose for the colinear one. As for SIRIUS, we have

As I wrote above, so far we have been dealing with the 'atomic' (and 'ortho-atomic') orbitals, which do not depend on the spin index. I haven't thought about the generalization to the case when they depend on the spin index for the relativistic case. But this is something that we want to do. One PhD student wants to work on the generalization to the non-collinear case.

It is a complex problem in the sense that atomic wave functions are indexed with the total angular momentum. I do not know about the beta projectors but I suspect something similar. if so, properly treating relativistic pps require working in the generalized spherical harmonics basis etc... But it is an interesting scientific question. I do not think that DFT should have a problem treating this but it might not be straightforward to implement a full relativistic DFT.

A general question about SIRIUS. As far as I understand, the logic of SIRIUS does not follow exactly the same logic of QE, correct? So basically you reimplement everything in a different way, correct? So SIRIUS is just a different code, not a GPU-enabled mirror of QE, right?

SIRIUS and QE are two different implementations of the same theory. SIRIUS can also treat other variants of potentials that QE does not support but the underlying theory is the same. The implementation however is (completely) different that's why you do not get the exact same numbers when you run QE and QE+SIRIUS. Moreover SIRIUS was implemented with GPUs in mind and this reflected in the structure of the code. Anton can correct me if I say something wrong.

haampie commented 3 years ago

As far as I understand, the logic of SIRIUS does not follow exactly the same logic of QE, correct?

To add to @mtaillefumier's comment; a few differences I know of: the mixer code in SIRIUS is different form QE, and in the Davidson solver SIRIUS implements eigenvector locking, which is not implemented in QE. Probably there's more differences.

mtaillefumier commented 3 years ago

On 10/27/20 1:30 PM, Harmen Stoppels wrote:

As far as I understand, the logic of SIRIUS does not follow
exactly the same logic of QE, correct?

To add to @mtaillefumier https://github.com/mtaillefumier's comment; a few differences I know of: the mixer code in SIRIUS is different form QE, and in the Davidson solver SIRIUS implements eigenvector locking, which is not implemented in QE. Probably there's more differences.

Indeed, it is something I omitted in my reply. But it is true. There are technical differences (sometimes major) between the two codes but they try to solve essentially the same problem. :)

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/electronic-structure/SIRIUS/issues/578#issuecomment-717210450, or unsubscribe https://github.com/notifications/unsubscribe-auth/AHAE5JN4JEP4RNOWFZBNOF3SM24NPANCNFSM4S7NI3QQ.

timrov commented 3 years ago

Thanks for your replies!

I do not think that DFT should have a problem treating this but it might not be straightforward to implement a full relativistic DFT.

As far as I know, there is no need in relativistic DFT. QE solves Pauli equations (for two-component spinors) and using fully-relativistic pseudopotentials. See e.g. sections 2.2.2(c) and 2.2.2(d) in my PhD thesis: https://pastel.archives-ouvertes.fr/pastel-00823758/

mtaillefumier commented 3 years ago

On 10/27/20 5:12 PM, Iurii Timrov wrote:

Thanks for your replies!

I do not think that DFT should have a problem treating this but it
might not be straightforward to implement a full relativistic DFT.

As far as I know, there is no need in relativistic DFT. QE solves Pauli equations (for two-component spinors) and using fully-relativistic pseudopotentials. See e.g. sections 2.2.2(c) and 2.2.2(d) in my PhD thesis: https://pastel.archives-ouvertes.fr/pastel-00823758/ https://pastel.archives-ouvertes.fr/pastel-00823758/

Hopefully not. Solving Pauli equations should be enough in many cases. The reading was interesting actually. I am not sure to understand what you call fully-relativistic pseudopotentials. If these pps generated from the Pauli equation for atoms or from the Dirac equation.

timrov commented 3 years ago

Fully-relativistic PPs are generated by solving a Dirac equation for an atom. But we do not need to solve a Dirac equation in a solid (instead we solve a Pauli-type equation + fully-relativistic PPs, which is much cheaper and accurate enough). I recommend to read the introduction of this paper: A. Dal Corso, Phys. Rev. B 82, 075116 (2010).

mtaillefumier commented 3 years ago

On 10/27/20 5:26 PM, Iurii Timrov wrote:

Fully-relativistic PPs are generated by solving a Dirac equation for an atom. But we do not need to solve a Dirac equation in a solid (instead we solve a Pauli-type equation + fully-relativistic PPs, which is much cheaper and accurate enough). I recommend to read the introduction of this paper: A. Dal Corso, Phys. Rev. B 82, 075116 (2010).

thanks, I will.

toxa81 commented 3 years ago

Iurii @timrov , can I assume that the case of spin-orbit + Hubbard is not going to be used in the nearest future? If it is important, how the localized atomic orbitals look like in this case? Are they pure spinors, constructed from the pseudo-atomic radial functions? Or something else?

timrov commented 3 years ago

can I assume that the case of spin-orbit + Hubbard is not going to be used in the nearest future?

Hard to say. There is a PhD student in our group who is planning to work on that. So, in the best case scenario, it can be implemented in ~1-2 years (but this is not guaranteed, because I do not know how it will go).

If it is important, how the localized atomic orbitals look like in this case? Are they pure spinors, constructed from the pseudo-atomic radial functions? Or something else?

I do not know, I haven't though about this yet. I would suggest to postpone this case for later.

toxa81 commented 3 years ago

@mtaillefumier I'm trying to formalize the JSON input schema for the U,V correction parts. We need to provide 1) on-site iteration terms 2) off-site interactions. My first attempt looks like this

{
    "U": [{
            "atom_type": "Fe",
            "l": 2,
            "U": 4.0
        },
        {
            "atom_type": "Fe",
            "l": 1,
            "U": 1.0
        }
    ],
    "V": [{
            "atom_pair": [1, 4],
            "l": [2, 1],
            "V": 0.5
        },
        {
            "atom_pair": [2, 4],
            "l": [1, 1],
            "V": 0.1
        }
    ]
}

Does it look reasonable/obvious?

mtaillefumier commented 3 years ago

It is reasonable I think. I doubt anyone will use this directly except us. One thing I would add is the displacement vectors because we may want to have atom pairs that extend further than the unit cell. if not given then we can maybe assume [0,0,0], if not we just have to indicate for instance

"displacement" : [0, 0, 1]

timrov commented 3 years ago

In the new version of QE we want to do something similar: to have [i,j,k] where i,j,k=-1,0,+1 to indicate neighbouring cells.

mtaillefumier commented 3 years ago

I think it is a reasonable way to describe interactions with neighboring cells. It also avoid computation of the atom coordinates inside the input file.

toxa81 commented 3 years ago

of course, I completely forgot about the the displacements!

{
    "U": [{
            "atom_type": "Fe",
            "l": 2,
            "U": 4.0
        },
        {
            "atom_type": "Fe",
            "l": 1,
            "U": 1.0
        }
    ],
    "V": [{
        "atoms": [1, 4],
        "T": [0, 0, 1],
        "V": 0.5,
        "l": [2, 1]
    }]

}
toxa81 commented 3 years ago

U is a list of on-site interactions, V is a list of off-site interactions. For V, is "l": [2, 1] self-evident?

timrov commented 3 years ago

"l": [2, 1]

Well, to me it is evident

mtaillefumier commented 3 years ago

to me as well.

On 2/1/21 2:05 PM, Iurii Timrov wrote:

"l": [2, 1]

Well, to me it is evident

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/electronic-structure/SIRIUS/issues/578#issuecomment-770842655, or unsubscribe https://github.com/notifications/unsubscribe-auth/AHAE5JL7FBRRTZRUNCPXSN3S42RKPANCNFSM4S7NI3QQ.

mtaillefumier commented 1 year ago

Should we close this issue ?

toxa81 commented 1 year ago
toxa81 commented 1 year ago

agree