VASP Forum

Hello! Recently, I have been performing electron-phonon renormalization calculations using VASP 6.5.1. I input the k-point coordinates in the KPOINTS_ELPH file in the following format:

k-points along high symmetry lines
40
Line-mode
rec
0.333333333333 0.333333333333 0.0 ! K
0.0 0.0 0.0 ! \Gamma

0.0 0.0 0.0 ! \Gamma
0.5 0.0 0.0 ! M

Unfortunately, I encountered the following error:

forrtl: severe (71): integer divide by zero
Image              PC                Routine            Line        Source             
vasp_std           0000000001FFB681  tbk_trace_stack_i     Unknown  Unknown
vasp_std           0000000001FF97BB  tbk_string_stack_     Unknown  Unknown
vasp_std           0000000001F61034  Unknown               Unknown  Unknown
vasp_std           0000000001F60E46  tbk_stack_trace       Unknown  Unknown
vasp_std           0000000001EE1AE7  for__issue_diagno     Unknown  Unknown
vasp_std           0000000001EEBA90  for__signal_handl     Unknown  Unknown
libpthread-2.17.s  00002B9AD11565E0  Unknown               Unknown  Unknown
vasp_std           000000000092AF02  Unknown               Unknown  Unknown
vasp_std           000000000092F935  Unknown               Unknown  Unknown
vasp_std           0000000001D102B1  Unknown               Unknown  Unknown

Could you please advise on how to resolve this issue? Thank you!

I cannot reproduce your error. I would need additional files like phelel_params.hdf5. So please upload all files from your directory where the crash occured.

https://github.com/Zhaoyihang-PRC/vasp6.5.1-example.git

As the file is relatively large, I have shared it via the link provided here. Thank you!

Thank you for reporting the error and providing the example.
I can reproduce your behaviour, a division by zero in a modulo function leads to the crash.

I'm not working on this part of the code so I will give it to my colleague. He is now on christmas vacation, so it will take some time until he can look at the problem.

Hello,

I'm from the electron-phonon team and I've had a look at your example. It is true that the code crashes in a most unpleasant way with this KPOINTS_ELPH file. Unfortunately, we do not support line-mode for this type of calculation at the moment. We should definitely update the code so that the crash is at least avoided and produces a meaningful error message. Thank you for bringing this to your attention.

May I ask you what you want to calculate with these settings? The KPOINTS_ELPH file is used to integrate the electron self-energy and with k-points in line mode, this will not give any meaningful results. For standard electron-phonon calculations, it is recommended to use an automatic regular k-point mesh in the KPOINTS_ELPH file to sample the first Brillouin zone.

A line mode is usually only reasonable for plotting/writing the electron-phonon matrix elements using ELPH_DRIVER = mels, but not for integrating the self-energy (ELPH_DRIVER=el).

Thank you for your reply! In fact, we intend to calculate the overall variation characteristics of the band structure along the high-symmetry k-point path by computing the real part of the electronic self-energy, and further plot the temperature-dependent band structures accordingly. With the conventional KPOINTS settings, we have no alternative but to manually identify the k-points that lie on the high-symmetry path from the IBZKPT file, manually extract the variation of their real parts, and then calculate and superimpose such variations onto the original band structure. This practice will not only reduce the sampling density of the band structure, but also compromise the computational efficiency significantly.

Thank you for your explanation. Unfortunately, it is currently not possible to select the k-points for self-energy calculations on a different mesh than the q-points, meaning that \(\mathbf{k}\) in the self-energy needs to exist on a regular k/q-point mesh (from KPOINTS_ELPH file), exactly as you describe.

However, things are not as dire as they might seem: In order to obtain a reliable self-energy, it needs to be converged with respect to the q-point density. This usually means going to very high sampling densities in the KPOINTS_ELPH file (for example 64x64x64). At these high densities, you will have enough k-points lying on the high-symmetry path that can be extracted from the calculation. This will not incur any additional computational cost as these points are already part of the grid anyway. Therefore, the procedure you describe is currently the intended way of doing what you want.

As an additional note, it is also necessary to increase the number of intermediate bands to a very large number in order to converge the real part of the self-energy. The default settings from ELPH_MODE = renorm already provide a convenient way to do this via ELPH_NBANDS = -2 and KPOINTS_OPT_MODE = 2. This computes many additional bands for electron-phonon interactions by diagonalizing the Hamiltonian at each k-point in KPOINTS_ELPH after the electronic convergence has been reached. Your current value of ELPH_NBANDS = 24 will most certainly not be enough to reach convergence. I encourage you to make a convergence study via ELPH_NBANDS_SUM to see if the number of bands is sufficient (see if the self-energy still changes when you reach ELPH_NBANDS).

Let me know if I was able to answer your question or if I can still help you with something.

Thank you for your reply! I have another question. Since I indeed want to obtain the real part of the electronic self-energy, theoretically, for accuracy, I need to set ELPH_NBANDS = -2 and KPOINTS_OPT_MODE = 2. However, the current calculation accuracy combined with a dense k-mesh seems to greatly affect computational efficiency and even consume a great deal of memory. Because I hope to get the self-energy for the full-space k-mesh, is there any solution to accelerate our calculations? Looking forward to your reply!

Meanwhile, I would also like to consult with you: if I need to extract the imaginary part of the electronic self-energy in the full space as the relaxation time for subsequent transport calculations, what essential parameters do we need to ensure accuracy? Thanks！

I understand that computing the self-energy at many points in the Brillouin zone can become quite expensive. However, I still suspect that there is something in your settings that causes the computation to be so expensive. Could you maybe share the INCAR and OUTCAR files of such a heavy computation with me so that I can have a look?

The biggest time saver in this case is to limit the number of bands that need to be computed for the self-energy (the index \(n\) in \(\Sigma_{n \mathbf{k}}\)). You can limit this with ELPH_SELFEN_BAND_START and ELPH_SELFEN_BAND_STOP. I have found these settings in your INCAR on github but they were commented out.