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
/