Non-physical solution of the polaron calculation of 2D blue-P

[jack]

Dear,
I am attempting to calculate the hole polaron in two-dimensional blue phosphorus (blue P). However, during the EPW calculation, an error message appears: "Non-physical solution, check initial guess and convergence."

The monolayer blue P structure contains two equivalent P atoms. The valence band comprises 5 bands without spin polarization.
First, I performed Wannierization using the projection: proj(1) = 'f=0.3333333330,0.6666666670,0.4705074384:p', with nbndsub = 3 and bands_skipped = 'exclude_bands = 1-2, 6-20'. The fitted bands agree well with the DFT results. It should be noted that wannier_plot_supercell was set to 6 6 2, while the coarse q-point and k-point grids were 6 6 1. If wannier_plot_supercell is set to 6 6 1 or 6*n 6*n 1, the Wannierization fails with the error: "Wrong range in plotting WF cube".

Next, I proceeded to calculate the hole polaron. Despite using the parameters listed below, EPW reported a non-physical solution. After reducing niter_plrn to 12, the calculation ran without error and yielded a formation energy with an error of approximately 0.05 eV. However, the resulting polaron wavefunction appears anomalous: it is localized around 14 Å along the Z-axis, whereas the atoms are positioned near 9 Å. This indicates that the wavefunction is centered outside the atomic sites.
How can I obtain the correct polaron wavefunction? Thank you.
The version of EPW is 6.0.
*****************epw1.in
prefix = 'blue'
outdir = './'
elph = .true.
epwwrite = .true.
lpolar = .true.
nbndsub = 3
dvscf_dir = './save'
etf_mem = 0
bands_skipped = 'exclude_bands = 1-2, 6-20'
wannierize = .true.
num_iter = 500
iprint = 2
proj(1) = 'f=0.3333333330,0.6666666670,0.4705074384:p'
wdata(1) = 'use_ws_distance = T'
dis_win_min = -17
dis_win_max = -5
!dis_froz_min = -3
dis_froz_max = -5
!system_2d = 'dipole_sp'
wannier_plot = .true.
wannier_plot_supercell = 6 6 1
!wannier_plot_radius = 2
nk1 = 6
nk2 = 6
nk3 = 1
nq1 = 6
nq2 = 6
nq3 = 1
band_plot = .true.
filkf = './kpoints'
filqf = './kpoints'
/
******************************epw2.in
&inputepw
prefix = 'blue'
outdir = './'
elph = .true.
epwread = .true.
lpolar = .true.
nbndsub = 3
dvscf_dir = './save/'
!system_2d = 'dipole_sp'
etf_mem = 0
wannierize = .false.
plrn = .true.
restart_plrn = .false.
type_plrn = 1
init_plrn = 1
!init_ntau_plrn = 1
init_sigma_plrn = 1
niter_plrn = 100
conv_thr_plrn = 5E-4
ethrdg_plrn = 5E-4
!full_diagon_plrn = .true.
!adapt_ethrdg_plrn = .true.
nk1 = 6
nk2 = 6
nk3 = 1
nq1 = 6
nq2 = 6
nq3 = 1
nkf1 = 12
nkf2 = 12
nkf3 = 1
nqf1 = 12
nqf2 = 12
nqf3 = 1
/

[jlb]

Dear jack,

Thank you for your post.

The error "Non-physical solution, check initial guess and convergence" is issued as a safeguard to prevent unphysical results when there are atomic displacements that exceed half of the lattice parameter. This usually indicates that there’s a problem in your calculation.

Can you please share the scf, ph, and nscf input files you used prior to the EPW calculation, so that we can reproduce the error?

Best,
Jon Lafuente-Bartolome

[jack]

Dear,
Thank you for your response. I have uploaded all input files.
The wannierization is implied with the "wannier_plot_supercell = 6 6 2", and the proj is set to P's p orbit. The fitting bands agree well with the DFT.
The polaron I calcualted is hole polaron.

[jlb]

Dear jack,

Thanks for sharing the files. Looking at the results obtained with your inputs, it seems that the spread of the Wannier functions is way too large.

In this case, I found that it is more convenient to build a Wannier function set comprising 8 bands (5 valence+3 conduction), so that P:s and P:p can be chosen as projections.

Then, if you only want to consider the topmost three valence bands in the polaron calculation as in the input files you shared, you can limit the bands that enter the solution of the polaron equations and for which the matrix elements are to be interpolated with the following input variables in epw2.in and epw3.in:

start_band_plrn = 3
end_band_plrn = 5
g_start_band_plrn = 3
g_end_band_plrn = 5

In any case, it is always good to make sure that the Wannier interpolation is correct not only by looking at the bands, but also looking at the electron-phonon matrix elements obtained by EPW (proper decay in Wannier representation, and comparison to direct DFPT calculations). You can find more information at the "Hands-On tutorial 2: Interpolation of electron-phonon matrix elements" of the 2024 summer school (https://docs.epw-code.org/doc/School2024.html).

If you are sure that all your calculations are well converged, and still you get the "non-physical solution" error, you can tune the maximum displacement allowed by modifying the input parameter "dtau_max_plrn", which is given in units of alat, and by default is set to 0.5d0. But please, use this with caution, as very large displacements might signal errors in the calculation, or that the solution is beyond the harmonic approximation.

Best regards,
Jon Lafuente-Bartolome

[jack]

Dear,
Thank you for your helpful response. I've noticed the larger spread of the Wannier functions and was attempting to generate them by including both valence and conduction bands. However, I was unsure how to properly exclude the conduction bands. Your reply has provided me with a clear guideline. Thank you for your guidance.
Best.