Inquiry about 2D Polar Correction and Negative Phonon

[balabala12s]

Dear Nakib,

I hope this message finds you well.

I am using elphbolt to calculate the phonon drag of 2D materials.When I set polar = .true. and twod = .true., I encounter the error: "Large negative phonon energy found! Stopping!" The code works fine if I turn the polar flag off.

I noticed the source code uses 3d_q functions (such as spectral_head_polarizability_3d_q). Does elphbolt support 2D-specific polar correction yet? If not, are there any other methods to resolve this numerical instability?

Thanks for your help!

Best regards,
Lee

[nakib]

Dear Lee,

Thank you for your message.

Have you checked if the phonons are stable when you generate them through EPW? We can proceed from there to identify the issue.

Those other subroutines you mentioned are not related to polar corrections.

Best,
Nakib

[balabala12s]

Dear Nakib,

Thank you for your guidance.

Regarding your question, I have double-checked the phonon dispersions from both ph.x and EPW, and they are fully stable with no imaginary frequencies throughout the Brillouin zone.

However, the crash in elphbolt only occurs when polar is set to .true. (with born and epsilon provided). If I disable the polar corrections, the code runs smoothly.

I suspect the issue might be related to the Q-point density: my current ph.x q-grid might be too coarse. Could this lead to numerical instabilities or unphysical (negative) values during the long-range dipole-dipole interpolation in elphbolt?

I would be very grateful for any suggestions or advice you may have on how to proceed.

Best regards,
Lee

[nakib]

Dear Lee,

EPW and elphbolt's phonon dispersions should be identical apart from very small numerical differences. This is what we always check. And it is a strong check that everything is consistent. Would you be able to tell me more about what acoustic sum rule corrections you are using in EPW? In elphbolt, we don't have anything other than the simplest acoustic sum rule correction in the style of ShengBTE. This could be one source of difference. The other source of difference could be how polar corrections are treated. If in EPW you are using a polar correction specifically designed for 2D systems (Coulomb cutoff, etc.) then that is also not something currently available in elphbolt (but it's on the list of issues that we want to implement in the future). Could you please check this first?

Best,
Nakib

[balabala12s]

Hi Nakib,

Thanks for the suggestions. Regarding your points:

ASR Correction: I am indeed using the simplest acoustic sum rule in EPW (asr_typ = 'simple'), so this should be consistent with the approach in Elphbolt.

Polar Correction: This is likely where the difference originates. In my EPW calculations, even after including the polar correction (lpolar = .true.), I noticed that the LO-TO splitting at the Gamma point for the optical branches is not explicitly shown in the phonon dispersion (This is different from the QE). For example, as seen in this forum case (https://forum.epw-code.org/viewtopic.php?p=3696&hilit=LO+TO+splitting#p3696), the results show that the LO-TO splitting of the optical branches does not appear in the phonon dispersion even after considering the polar correction.

Could you let me know if there are any methods to correct this error? Could this error cause Elphbolt to show huge negative phonon energies?

Best,
Lee

[nakib]

Hi Lee,

Regarding the simple ASR correction, I agree.

Which version of EPW are you using? If you have questions regarding the EPW interpolation of phonon energies, it would be best if you discuss it in their forum.

I am aware of one recent work by Quan and Liao (https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.134.226301) that used the EPW+elphbolt workflow on a polar 2D material. So it seems like it is possible to treat these types of systems with the current version of elphbolt. But I do not have any specific tips at the moment apart from trying making the DFTP q-mesh finer and/or playing around with QE's recent ASRs that strictly enforce the rotational invariance.

We will work on the special 2D polar corrections and will announce when these are ready for public use.

Best,
Nakib

[balabala12s]

Dear Nakib,

I apologize for my very late reply, and thank you very much for your helpful response.

Following your suggestion, I contacted the authors of the recent work by Quan and Liao (Phys. Rev. Lett. 134, 226301). They confirmed that they did not modify the EPW part; instead, they modified elphbolt to include the Fröhlich e‑ph interaction for 2D MoS₂ using Eq. 62 of Phys. Rev. B 105, 115414, with the material parameters provided in that paper.

I have now implemented the same modification in elphbolt. Specifically, I modified wannier.f90 by adding the 2D Fröhlich correction (the long_range_prefac subroutine) and modified crystal.f90 to read the LO phonon frequency and the effective dielectric thickness from the input file. I also removed the vacuum contribution from the static and high‑frequency dielectric constants using the capacitor‑stack model (Eqs. 64‑66 of the same PRB paper). My input file and the modified source files are attached.

Unfortunately, after performing the coupled electron‑phonon transport calculations, my results for MoS₂ are still very far from those reported in the PRL paper.

I would greatly appreciate it if you could have a look at my modifications and input file. Do you see any obvious mistake in my implementation? Could the problem be related to other parts of the code? Or is there any additional parameter that I might have missed?

Thank you so much for your time and help. I look forward to any advice you might have.

Best regards,
code and result.zip