ABINIT Discussion Forums

Hello Abifriends,

I asked a stupid question a few days ago that, after some thinking, I deleted. Sorry about that. Now I have what is a different but maybe still stupid question. Also, I have a report on what is probably a bug. I am using v9.4.0

The bug: If I want to read a WFQ file with optdriver=7 and eph_task=2 or eph_task=-2, Abinit ignores the variable getwfq_filepath. Even with it set in the input file, Abinit complains that getwfq and irdwfq are 0. I can set getwfq and it will read the WFQ file, but getwfq_filepath has no affect.

The other issue:
I am trying to get the gkq electron-phonon matrix elements on the 'arbitrarily dense' grid. We need the k and k+q wave functions. According to section 2 in the respfn help file (https://docs.abinit.org/guide/respfn/), q-points that are the difference between points on the k-grid are already contained in the WFK file. This makes sense. So I pick any k-point in the WFK file as my q-point, which should be a since the file contains Gamma a q=k-gamma=k. I set qpt=... in the input file. Abinit crashes complaining that some k+q point is not contained in the file:

--- !ERROR
src_file: m_bz_mesh.F90
src_line: 2352
mpi_rank: 5
message: |
q = k-kp+G0 not found. kmkp = -0.33333 0.20833 0.12500
...

I think Abinit is adding the q-point to only the irreducible part of k-point grid. I tried calculating the WFK file with kptopt=3 to include all points, but that also didn't work. I can do an nscf calculation with qpt set in the nscf part and use that WFQ file and everything works okay, but having to do ~10000 nscf calculations on a dense k-grid will be impossible. If I am doing something stupid, please berate me. Otherwise, can you please let me know how to fix my issue? I want the gkq elements in the full BZ, which Abinit apparently all ready uses to calculate phonon linewidths on an arbitrary q-path. I just don't seem to be able to access the gkq elements.

Thanks in advance!
Ty

To anyone who is interested: the first thing was 100% a bug. The problem was in m_eph_driver.F90. The variable 'use_wfq' did not check if getwfq_filepath was defined (it only checked getwfq and irdwfq), so it overwrote the input option to the WFK even if getwfq_filepath was in the input file. I fixed it and will try to push the source to github.

I don't think the other thing is a bug. I looked at the abipy workflow since I can't find it in the documentation elsewhere. It looks like doing an nscf calculation at the relevant q-point is the way to do it. It seems unnecessary since the respfn code is able to tell if the q-point is commensurate with the k-point grid. Some comments in the source code look related to this issue, but this is literally my first time looking at the Abinit source so I have no idea whats going on!. That being said, I am going to try to fix it

Ty

Thanks for reporting the problem and fixing the bug.

I have to say that the two options eph_task = -2, 2 are not very well documented.
Let's try to clarify the main purpose of these two options.

eph_task = +2 is mainly used to produce a netcdf file with the matrix elements <k+q,b'|d/d_p V_scf| k, b>
associated to a single atomic perturbation p where p is a shorthand notation for (qpt, idir, ipert).
Although the file produced by this option has extension GKK.nc, this file **cannot** be used by the e-ph part of anaddb.
The main purpose of this option is to generate netcdf files for
the [ElectronPhononCoupling](https://github.com/GkAntonius/ElectronPhononCoupling) python package
that computes the e-ph self-energy and QP corrections.
Some of these features are now implemented in Abinit by eph_task + 4
(see e.g. the [eph4zpr](https://docs.abinit.org/tutorial/eph4zpr/) tutorial).

eph_task = -2, on the other hand, produces a GKQ.nc file with the same e-ph matrix elements in the atomic representation but in a slightly different format.
The main purpose of this task is to generate a netcdf file with e-ph matrix elements with q-points
along a high-symmetry path
so that one can then use abipy to visualize |g(k, q, nu, b, b')| as a function of q for fixed k.
We usually perform N_qpath calculations in which we pass a different WKQ file computed on a shifted k+q mesh.
Note that one has to average the e-ph matrix elements over degenerate phonon modes and degenerate |k+q, b'>, |k, b> bands
(if any) in order to obtain a gauge-independent quantity.
We don't use this option very often because for our recent studies we found that it is much easier to look at the
unit-cell averaged lattice-periodic part of the scattering potential or to Fourier components with G /= 0.
See e.g. Eq 4 and Fig 2 of https://journals.aps.org/prl/pdf/10.110 ... 125.136601 and the SI.
The scattering potential is indeed gauge-independent and important features such as the Frohlich divergence or the jump discontinuity for q --> 0 in semiconductors due to dynamical quadrupoles are clearly visible at the level of the scattering potential
without having to plot |g(q)|^2 directly.
This task corresponds to eph_task +15 and we usually use abipy tools to plot the results.

To summarize, eph_task = +2 is mainly used to interface abinit with a ElectronPhononCoupling for ZPR calculations
while eph_task -2 is used for visualizing the q-dependence of the e-ph matrix elements.
eph_task +15 (analysis of the scattering potentials in q-space) is much easier to use since there's no gauge in the potentials
but this kind of visualization does not give any insight into the contribution to the e-ph matrix elements
coming from the Bloch orbitals. For instance, you may have that particular e-ph matrix elements are zero
due to the symmetry of the scattering potential and of the |k+q>, |k> states but this effect won't show up
if you plot the scattering potential alone.

I looked at the abipy workflow since I can't find it in the documentation elsewhere.
It looks like doing an nscf calculation at the relevant q-point is the way to do it.
It seems unnecessary since the respfn code is able to tell if the q-point is commensurate with the k-point grid.

This step is unecessary if your goal is to have <k',b'|\Delta_{q\nu} V}|k, b> for all k, k' belonging to the ngkpt k-mesh associated to the input WFK. It is needed if you want to compute <k+q,b'|\Delta_{q\nu} V}|k, b> when
q belongs to a high-symmetry path.
The AbiPy workflow was designed for the later case i.e. visualization purposes.

I want the gkq elements in the full BZ,

Why do you need the gkq elements in the full BZ?
For visualization purposes or for computing integrals in the BZ?
The EPH code is in principe able to compute g(k, q) for all k and k'=k+q in the full BZ once you provide an input WFK file
and set eph_ngqp_file == ngkpt but there's no eph_task that is explicitly designed with this kind of operation in mind.
When we compute g(k, q) for all the k- and q-points in the "BZ", we assume the user wants to obtain physical observables
such as QP energies, phonon linewidths, etc but we never precompute and store all the e-ph matrix elements
in memory for performance reasons.
The BZ integration in k/q space of the |g|^2 terms is always done on the fly, possibly by multiple MPI processes.
Only the final results are presented to the user.
On the one hand, I understand that having the full g(k, q, ...) array stored on disk can be useful
if you want to produce cool 2d or 3d plots with |g(k, q)|^2 for fixed k/q or if you need to read the e-ph matrix elements to implement some kind of post-processing step in a high-level language.
On the other hand, writing all this stuff in parallel is not trivial and that's the reason why I would like
to have a better understanding of your problem.

which Abinit apparently all ready uses to calculate phonon linewidths on an arbitrary q-path.
I just don't seem to be able to access the gkq elements.

Actually the phonon linewidths are interpolated using a Fourier-interpolation scheme similar in spirit to the one used for the dynamical matrix.
In a nuthshell, the code computes the phonon linewidths gamma(q) in the IBZ (actually a tensor from which one
can obtain the ph-linewidths but this is not relevant for the discussion).
Then we "rotate" the gamma(q) terms to get the values in the BZ. At this point we can do a q --> R Fourier transform to get gamma(R) in real space and then use the inverse transform to *Fourier interpolate* the signal at arbitrary q-points.
Note that e-ph matrix elements are not needed at this level.
You can interpolate the scattering potentials to obtain gamma(q) on a q-grid that is much denser than the one used at the DFPT level but the q-points still belong to a regular q-mesh so you need an extra level of interpolation
to go from the q-mesh to an *arbitrary* q-path.

I hope this addresses some of the questions raised in your previous posts.
Unfortunately, you are exploring parts of the EPH code that have been mainly used for our internal studies.
Some of these features are not fully documented and the API as well as the netcdf file formats may change in the near future.
We plan to write a tutorial explaining how to perform this kind of analysis but before doing that we need to polish a bit
our python scripts so that they can be easily used by other users.

First of all .. wow. Thanks for such an incredibly thorough response. You answered all of my questions and more and I'm sure this info will be useful for other curious users of Abinit.

Thanks for reporting the problem and fixing the bug.

I am glad I got to help, even for something minor. I will try to help again in the future.

Why do you need the gkq elements in the full BZ?

or if you need to read the e-ph matrix elements to implement some kind of post-processing step in a high-level language.

For now, I want to use the eph matrix elements for post processing. Specifically, I want to compute, analyze, and modify the eph contribution to the force constants (see e.g. equation (1) in https://doi.org/10.1103/PhysRevB.61.12059). I want to do some calculations similar in spirit to that paper, but with different applications in mind. I was going to write my own python code from scratch because I wasn't aware of the ElectronPhononCoupling package. Thanks a lot for making its existence known to me. I can probably extend it for my purposes rather than start from scratch. I am interested in metals (for now) so am not too worried about correctly treating the LR part.

I really only need the eph matrix elements at *arbitrary* q-points for this, but I wanted the dense mesh in the full BZ to use in something like a Holstein model (a possible future project). I don't really know if this is feasible yet, though, so the matrix elements at arbitrary q-points using the WFQ files will suffice.

Actually the phonon linewidths are interpolated using a Fourier-interpolation scheme similar in spirit to the one used for the dynamical matrix.

I actually figured out I was mistaken by looking at the source code and realizing that Abinit was not doing what I thought. I figured that out yesterday. Sorry for not updating my post and saving you the trouble of explaining it.

Note that one has to average the e-ph matrix elements over degenerate phonon modes and degenerate |k+q, b'>, |k, b> bands (if any) in order to obtain a gauge-independent quantity.

Apparently I need to some studying. If convenient, could you suggest a paper(s) relevant to this comment? If not, I will do my own research... which I guess is my job after all ...

Some of these features are not fully documented and the API as well as the netcdf file formats may change in the near future. We plan to write a tutorial explaining how to perform this kind of analysis but before doing that we need to polish a bit our python scripts so that they can be easily used by other users.

I look forward to this!!


Again, thanks a lot for the detailed answers and, in general, for working so hard to provide Abinit to the world.