[Pw_forum] modified explicit exchange potential

[degironc]

Dear Helen,

Helen wrote:
> Hi,
>
> I sent a previous message which wasn't answered, but which I've since 
> been working on, so I think I have found some answers myself. I'd 
> really appreciate Dr. de Gironcoli's input as I'm sure he is a lot 
> more knowledgable than I am.
>
> 1. I asked a question about band-structure calculations when using 
> exx.f90, as in the examples it states that only k-points used in the 
> original scf calculation can be used to calculate band structure. I 
> think a way to get round this would be to in the scf input file, first 
> list the k-points with normal weights for scf calculation, followed by 
> a list of k-points for band structure calculation with very small 
> weights. This way the values at these points are calculated but they 
> won't effect any BZ integrations. (For direct band gaps just use 
> Monkhorst-pack k-point generator which includes k=(0, 0, 0)). 

In the present status of the EXX part calculating a "simple" 
band-structure is not so simple.
It is certainly possible (as you say) to perform a calculation with a 
regular grid that contains gamma and
most high symmetry points (specifying nbnd>nelec/2 in order to have 
access to some empty bands).
This will also contains some point along (001) and can give an idea of 
the dispersion.
I do not see any easy way (actually I dont see any way) by which one can 
calculate the band structure in a point that does not belong to a 
regular grid without implementing something new.

The problem is that Vx is defined as a sum over the BZ

<r|Vx|psi_k> = C * sum_q f(k+q) sum_G <r|exp(-i(q+G)r|phi_k+q> 
<phi_k+q|exp(+i(q+G)r'|psi_k>/(q+G)2

where phi_k+q are the self-consistently calculated occupied bands in k+q 
and q is a point in a regular (unshifted) grid.
This implies that in order to calculate the bands in k the occupied 
bands in all k+q are needed..which you tipically don't have

A possibility would be to change the definition of Vx to

<r|Vx|psi_k> = C * sum_k' f(k') sum_G <r|exp(-i(k'-k+G)r|phi_k'> 
<phi_k'|exp(+i(k'-k+G)r'|psi_k>/(k'-k+G)2
where now k' is a regular (shifted or unshifted) grid.
The point here is that the grid of k'-points is always the same as k 
changes, as opposed to the grid of k+q that is always different.. this 
would allow to do a band calculation in a simpler way because one could 
keep the occupied phi_k'  from the scf calculation.

This requires some coding in exx and probably some reorganization in the 
management of files in order not to overwrite the files containing the 
phi_k' in a nscf calculation...
I'm not sure if the two calculations (sum over q or sum over k') are 
strictly equivalent in the scf case when one tries to reduce the BZ sum 
in Vx to a subset of the points... I think they are not equivalent but I 
don't know which one is "better".

Another option could be: 1) calculate the (occupied and some empty) 
bands in a regular grid of points, 2) use the wannier90 code to extract 
the maximally localized wannier-functions, 3) use them to build a 
tight-binding representation of the hamiltonian that allows you to 4) 
fourier interpolate the bandstructure throughout the BZ...

It looks complicated but it could be the simplest option...
and if automated and shown to be accurate it could actually be used also 
to speed up the EXX part in general...maybe
> 2. The error involved in discarding the q=0 term
> goes as 1/(NQS*Omega). This fact is utilized in removing the error. I 
> think this will still be valid if the e- interaction is not U(r)=1/r but
> U(r)=efr(gam r)/r ? [Fourier Transform=exp(-q2/4*gam2)*4Pi/q2, so
> divergence as q->0 is same as U(r)=1/r, so I think they should be the 
> same)
the q=0 therm is problematic for two reasons: i) it contains an 
integrable divergence and this is taken care of by the Gygi-Baldereschi 
trick and ii) the remaining non-divergent term in the q->0 cannot be 
calculated
in q=0 because it is a 0/0 limit. This is why there is the extrapolation 
trick. I think this still remain a problem even if U(r) is erf(gam r)/r 
because the problem is the leading q2 dependence of 
|<psi_k|exp(-iqr)|phi_k+q>|^2 which is not easy to extract...

If instead one would  like to replace the 1/r with a short range 
interaction (yukawa or other) this limit should be well defined and one 
would not need extrapolation. For this and other reasons hybrid functionals
with a screened short-range non-local exchange might have some 
computational advantages.

> 3. I think that in exx.f90 all I need to do to account for a 
> descreened potential is to change fac(ig):
>                 fac(ig)=e2*fpi*dexp(-qq/(4d0*gam))/(tpiba2*qq + yukawa 
> ) * grid_factor
> I did this in the routines exxenergy2 and in vexx. I haven't touched 
> exx_divergence, do I need to? Are there any other places which need to 
> be changed?
I'm not sure... exx_divergence is basically the Ewald energy of a set of 
unit point charges located at the Bravais lattice vectors+uniform 
compensating background. If you change the interaction I expect this 
number to change... try to see how your exchange energy converge for an 
isolated system (an atom or a molecule) as a function on the box-size, 
w/o the extrapolation trick. If the divergence is correctly calculated
it should converge as 1/omega, if not I would expect something like 
1/boxsize.

best regards,

Stefano de Gironcoli, SISSA & DEMOCRITOS, Trieste (Italy)