[WIP] Hybrid functionals

[BKaperick]

Started work at the DFTK Summer School with @RashidRafeek and @ChRickert in response to issue #530.

• [x] Initialize fock exchange term and operator
• [x] Compute energy for 1 kpoint
• [x] Compute Hamiltonian
• [ ] Extend LibXc and DftFunctionals interface to allow hybrid families
• [ ] Add standard models like PBE0 and add test case
• [ ] Generalize to fractional occupations
• [ ] Make energy term work with arbitrary number of kpoints
• [ ] Generalize SCF to mix orbitals
• [ ] Implement ACE method for SCF for better performance

[antoine-levitt]

Looks super, great work!

Extend energy code to work for arbitrary number of kpoints

Should be OK to do: when i = (n,k) and j = (n',k'), rho_ij = psi_i^* psi_j is exp(i(k'-k)x) u_nk*(x) u_mk(x), so V_ij can be computed in Fourier space just as you do now but with a |k'-k+G|^2 instead of G^2. That should possibly be the only change needed, except for normalizations maybe.

For testing, we should be able to compare hartree-fock against eg abinit.

additional TODOs (not necessarily in the scope of this PR, I'm just adding for future reference) is to generalize SCF to mix orbitals, and to implement the ACE method for performance.

[antoine-levitt]

also as a TODO: generalize to fractional occupations

[antoine-levitt]

Note that the BZ integral is singular and needs to be regularized, either by short-ranging the Coulomb (this is what abinit does, apparently) or by a correction technique, see eg https://journals.aps.org/prb/pdf/10.1103/PhysRevB.34.4405 (this is what QE does by default). There's also https://journals.aps.org/prb/pdf/10.1103/PhysRevB.73.205119 which I haven't read yet, and probably more.

[antoine-levitt]

See also https://www.vasp.at/wiki/index.php/Coulomb_singularity

[zhubonan]

Hi @antoine-levitt @mfherbst,

I am interested in implementing the hybrid functional with k-point sampling. Actually have tried a bit, but there are quite a few things to be done such as getting the coulomb singularity correction computed and generating the q point set and corresponding wave functions. I am wondering if this is already being actively worked on and what is the best way to contribute.

[antoine-levitt]

Hi,

As far as I know nobody's working on it actively, feel free to pick it up (the simplest is probably to start another PR; unless @BKaperick objects I think it's fine to use some of the code in this PR if you want, with proper attribution in a commit message). Getting this to work would be awesome!

Regarding the Coulomb the simplest solution to get something working is probably to just do a damped Coulomb.

For the set of q points, I'm not sure a specific set needs to be computed, can't you just use the existing kpoint grid? Or do you mean in the presence of symmetries? Note also that phonons are being actively worked on (cc @epolack), they might share some of the technology?

Feel free to get in touch (here, by email or on zulip)

[zhubonan]

Thanks for the reply! I will try to pick it up and extend the functionality to multiple k points.

For the set of q points, I'm not sure a specific set needs to be computed, can't you just use the existing kpoint grid? Or do you mean in the presence of symmetries?

Yes, the existing grid can be used, but if the actual kpoints are reduced by symmetry only those in the irreducible BZ are calculated. However, the double summation requires the full grid to be used for the q points. So for each q we need to get its wave function from the symmetry related k point. This requires some transformation of the wave function (see https://docs.abinit.org/theory/wavefunctions/#symmetry-properties).

Note also that phonons are being actively worked on (cc @epolack), they might share some of the technology?

Good call, I will take a look there.

[antoine-levitt]

I have not looked closely if symmetries can be used for exchange computations, presumably they can and then the symmetric structure can be used? We use things like apply_symop (in symmetries.jl) to convert when needed. If they can't (or as a first step before tackling the symmetries), it's reasonable to just do a fully non-symmetry-adapted calculation (symmetries=false in Model()). See also unfold_bz(scfres) in symmetry.jl.

[BKaperick]

Hi,

As far as I know nobody's working on it actively, feel free to pick it up (the simplest is probably to start another PR; unless @BKaperick objects I think it's fine to use some of the code in this PR if you want, with proper attribution in a commit message). Getting this to work would be awesome!

Regarding the Coulomb the simplest solution to get something working is probably to just do a damped Coulomb.

For the set of q points, I'm not sure a specific set needs to be computed, can't you just use the existing kpoint grid? Or do you mean in the presence of symmetries? Note also that phonons are being actively worked on (cc @epolack), they might share some of the technology?

Feel free to get in touch (here, by email or on zulip)

Yes, sorry for letting this drag on @mfherbst and @antoine-levitt , I just could not find the time to get it working. Thank you for taking over the topic, @zhubonan. And of course you can use whatever code is useful from this PR.

[mfherbst]

Superseeded by #942.