scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'polaron_v5.7' into 'EPWv5.7a' 

Polaron v5.7

See merge request epw-code/q-e!37
This commit is contained in:
Hyungjun Lee 2023-03-13 15:17:46 +00:00
parent 4700473d38 48da98913a
commit 2e692a43be
18 changed files with 4806 additions and 1757 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
EPW/CMakeLists.txt
View File

@ -51,9 +51,8 @@ set(sources
src/loadbm.f90
src/bfield.f90
src/io_indabs.f90
src/ephblochkq.f90
src/wfc_elec.f90
src/test_tools.f90)
src/polaron.f90
)
qe_add_library(qe_epw ${sources})
target_link_libraries(

4
EPW/src/Makefile
View File

@ -71,8 +71,8 @@ epw_summary.o \
loadumat.o \
stop_epw.o \
wannierEPW.o \
loadbm.o
EPWOBJS += ephblochkq.o wfc_elec.o test_tools.o
loadbm.o \
polaron.o \
#default : epw

96
EPW/src/bcast_epw_input.f90
View File

@ -78,16 +78,21 @@
USE io_global, ONLY : meta_ionode_id
USE control_flags, ONLY : iverbosity
USE ions_base, ONLY : amass
! ---------------------------------------------------------------------------------
! Added for polaron calculations. Originally by Danny Sio, modified by Chao Lian.
! Shell implementation for future use.
USE epwcom, ONLY : wfcelec, model_vertex , polaron_wf, r01, r02, r03,&
num_cbands, start_band, start_mode, cb_shift, &
polaron_interpol, polaron_bq, polaron_dos, &
electron_dos, phonon_dos, diag_mode, &
restart_polaron_mode, polaron_type, &
emax_plrn, nDOS_plrn, emin_plrn
! -------------------------------------------------------------------------------
! Added for polaron calculations by Chao Lian
USE epwcom, ONLY : plrn, restart_plrn, conv_thr_plrn, end_band_plrn, lrot, &
cal_psir_plrn, start_band_plrn, type_plrn, nstate_plrn, &
interp_Ank_plrn, interp_Bqu_plrn, init_sigma_plrn, &
full_diagon_plrn, mixing_Plrn, init_plrn, niter_plrn, &
nDOS_plrn, edos_max_plrn, edos_min_plrn, edos_sigma_plrn, &
pdos_sigma_plrn, pdos_max_plrn, pdos_min_plrn, &
seed_plrn, ethrdg_plrn, time_rev_A_plrn, nhblock_plrn, &
beta_plrn, Mmn_plrn, recal_Mmn_plrn, r0_plrn, debug_plrn, &
time_rev_U_plrn, g_start_band_plrn, g_end_band_plrn, &
g_start_energy_plrn, g_end_energy_plrn, &
model_vertex_plrn, model_enband_plrn, model_phfreq_plrn, &
kappa_plrn, omega_LO_plrn, m_eff_plrn, step_wf_grid_plrn, &
scell_mat_plrn, scell_mat
!
!
IMPLICIT NONE
!
@ -300,29 +305,54 @@
CALL mp_bcast(scdm_entanglement, meta_ionode_id, world_comm)
!
! ---------------------------------------------------------------------------------
! Added for polaron calculations. Originally by Danny Sio, modified by Chao Lian.
! Shell implementation for future use.
CALL mp_bcast (wfcelec , meta_ionode_id, world_comm)
CALL mp_bcast (model_vertex , meta_ionode_id, world_comm)
CALL mp_bcast (polaron_wf , meta_ionode_id, world_comm)
CALL mp_bcast (polaron_interpol, meta_ionode_id, world_comm)
CALL mp_bcast (polaron_bq , meta_ionode_id, world_comm)
CALL mp_bcast (polaron_dos , meta_ionode_id, world_comm)
CALL mp_bcast (electron_dos , meta_ionode_id, world_comm)
CALL mp_bcast (phonon_dos , meta_ionode_id, world_comm)
CALL mp_bcast (num_cbands , meta_ionode_id, world_comm)
CALL mp_bcast (start_band , meta_ionode_id, world_comm)
CALL mp_bcast (start_mode , meta_ionode_id, world_comm)
CALL mp_bcast (cb_shift , meta_ionode_id, world_comm)
CALL mp_bcast (diag_mode , meta_ionode_id, world_comm)
CALL mp_bcast (restart_polaron_mode, meta_ionode_id, world_comm)
CALL mp_bcast (polaron_type, meta_ionode_id, world_comm)
CALL mp_bcast (r01 , meta_ionode_id, world_comm)
CALL mp_bcast (r02 , meta_ionode_id, world_comm)
CALL mp_bcast (r03 , meta_ionode_id, world_comm)
CALL mp_bcast (nDOS_plrn , meta_ionode_id, world_comm)
CALL mp_bcast (emax_plrn , meta_ionode_id, world_comm)
CALL mp_bcast (emin_plrn , meta_ionode_id, world_comm)
! Added for polaron calculations by Chao Lian.
CALL mp_bcast(lrot , meta_ionode_id, world_comm)
CALL mp_bcast(plrn , meta_ionode_id, world_comm)
CALL mp_bcast(cal_psir_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(interp_Ank_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(interp_Bqu_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(start_band_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(end_band_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(type_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(nDOS_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(edos_max_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(edos_min_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(nstate_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(niter_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(conv_thr_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(restart_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(type_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(init_sigma_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(ethrdg_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(time_rev_A_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(time_rev_U_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(debug_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(full_diagon_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(mixing_Plrn, meta_ionode_id, world_comm)
CAll mp_bcast(init_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(Mmn_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(recal_Mmn_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(r0_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(edos_sigma_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(pdos_sigma_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(pdos_max_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(pdos_min_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(seed_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(nhblock_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(beta_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(g_start_band_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(g_end_band_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(g_start_energy_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(g_end_energy_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(step_wf_grid_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(model_vertex_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(model_enband_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(model_phfreq_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(kappa_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(omega_LO_plrn, meta_ionode_id, world_comm)
CAll mp_bcast(m_eff_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(scell_mat_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(scell_mat, meta_ionode_id, world_comm)
! --------------------------------------------------------------------------------
#endif
!

1010
EPW/src/ephblochkq.f90
View File

File diff suppressed because it is too large Load Diff

206
EPW/src/ephwann_shuffle.f90
View File

@ -118,18 +118,14 @@
USE parallel_include, ONLY : MPI_MODE_RDONLY, MPI_INFO_NULL, MPI_OFFSET_KIND, &
MPI_OFFSET
#endif
! ---------------------------------------------------------------------------------
! Added for polaron calculations. Originally by Danny Sio, modified by Chao Lian.
! Shell implementation for future use.
USE epwcom, ONLY : wfcelec, start_band, polaron_wf, restart_polaron, &
polaron_interpol, polaron_bq, polaron_dos, nPlrn, &
wfcelec_old
USE elph2, ONLY : g2_4
USE ephblochkq, ONLY : interpol_bq, interpol_a_k, compute_a_re
USE polaron, ONLY : wfc_elec, epfall, ufall, Hamil, eigVec, &
interp_plrn_wf, interp_plrn_bq, plot_plrn_wf
USE polaron_old, ONLY : wfc_elec_old
! --------------------------------------------------------------------------------
!!!!!
! Added for polaron calculations.
USE epwcom, ONLY : plrn, time_rev_U_plrn, g_start_band_plrn, g_end_band_plrn
USE epwcom, ONLY : scell_mat_plrn
USE polaron, ONLY : plrn_flow_select, plrn_prepare, plrn_save_g_to_file
USE polaron, ONLY : is_mirror_q, is_mirror_k
USE polaron, ONLY : kpg_map, ikq_all
!!!!!
!
IMPLICIT NONE
!
@ -338,6 +334,9 @@
COMPLEX(KIND = DP), ALLOCATABLE :: eimpmatf(:, :)
!! carrier-ionized impurity matrix in smooth Bloch basis
!!!!!
LOGICAL :: mirror_k, mirror_q, mirror_kpq, ik_global
REAL(KIND = DP) :: xxq_r(3)
!!!!!
!
CALL start_clock('ephwann')
!
@ -975,8 +974,8 @@
!
totq = 0
!
IF (wfcelec) THEN
!
IF (plrn .OR. scell_mat_plrn) THEN
! For polaron calculations, all the q points have to be included
totq = nqf
ALLOCATE(selecq(nqf), STAT = ierr)
IF (ierr /= 0) CALL errore('ephwann_shuffle', 'Error allocating selecq', 1)
@ -984,7 +983,7 @@
selecq(iq) = iq
ENDDO
!
ELSE ! wfcelec
ELSE
! Check if the file has been pre-computed
IF (mpime == ionode_id) THEN
INQUIRE(FILE = 'selecq.fmt', EXIST = exst)
@ -1012,7 +1011,7 @@
WRITE(stdout, '(5x,a,i8,a)')'We only need to compute ', totq, ' q-points'
WRITE(stdout, '(5x,a)')' '
!
ENDIF ! wfcelec
END IF ! plrn
!
! -----------------------------------------------------------------------
! Possible restart during step 1)
@ -1108,38 +1107,8 @@
IF (ierr /= 0) CALL errore('ephwann_shiffle', 'Error allocating epsilon2_abs_lorenz_all', 1)
ENDIF ! indabs
!
! --------------------------------------------------------------------------------------
! Polaron shell implementation for future use
IF (wfcelec) then
IF (polaron_interpol) THEN
ALLOCATE(eigVec(nktotf * nbndfst, nplrn), STAT = ierr)
IF (ierr /= 0) CALL errore('ephwann_shuffle', 'Error allocating eigVec', 1)
eigVec = czero
CALL interp_plrn_wf(nrr_k, ndegen_k, irvec_r, dims)
iq_restart = totq ! Skip the calculation of e-ph element, save the time.
DEALLOCATE(eigVec)
ELSEIF(polaron_bq) THEN
CALL interp_plrn_bq(nrr_q, ndegen_q, irvec_q)
iq_restart = totq ! Skip the calculation of e-ph element, save the time.
ELSEIF(polaron_wf) THEN
CALL plot_plrn_wf()
iq_restart = totq
ELSE
ALLOCATE(eigVec(nktotf * nbndfst, nplrn), STAT = ierr)
IF (ierr /= 0) CALL errore('ephwann_shuffle', 'Error allocating eigVec', 1)
eigVec = czero
ALLOCATE(epfall(nbndfst, nbndfst, nmodes, nkf, nqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('ephwann_shuffle', 'Error allocating epfall', 1)
epfall = czero
ALLOCATE(ufall(nmodes, nmodes, nqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('ephwann_shuffle', 'Error allocating ufall', 1)
ufall = czero
ALLOCATE(Hamil(nkf * nbndfst, nktotf * nbndfst), STAT = ierr)
IF (ierr /= 0) CALL errore('ephwann_shuffle', 'Error allocating Hamil', 1)
Hamil = czero
ENDIF
ENDIF
! -------------------------------------------------------------------------------------
! Polaron calculations: iq_restart may be set to totq in restart/interpolation mode.
IF (plrn) call plrn_prepare(totq, iq_restart)
!
! Restart in SERTA case or self-energy (electron or plasmon) case
IF (restart) THEN
@ -1301,6 +1270,23 @@
ELSE
xxq = xqf(:, iq)
ENDIF
!!!!!!!
! Added by Chao Lian for enforcing the time-rev symmetry of e_q
! for polaron calculations, xxq_r is coordinates of the mirror point of xxq
! for other calculations, xxq_r is xxq.
IF (plrn) THEN
IF (is_mirror_q (iq)) THEN
xxq_r = xqf(:, kpg_map(iq))
mirror_q = .true.
ELSE
xxq_r = xxq
mirror_q = .false.
END IF
ELSE
xxq_r = xxq
mirror_q = .false.
end if
!!!!!!!
! Temporarily commented by H. Lee
! CALL find_gmin(xxq)
!
@ -1309,9 +1295,16 @@
! ------------------------------------------------------
!
IF (.NOT. lifc) THEN
CALL dynwan2bloch(nmodes, nrr_q, irvec_q, ndegen_q, xxq, uf, w2)
!!!!!
! CALL dynwan2bloch(nmodes, nrr_q, irvec_q, ndegen_q, xxq, uf, w2)
CALL dynwan2bloch(nmodes, nrr_q, irvec_q, ndegen_q, xxq_r, uf, w2, mirror_q)
!!!!!
ELSE
CALL dynifc2blochf(nmodes, rws, nrws, xxq, uf, w2)
!!!!!
!TODO: apply degeneracy lift in dynifc2blochf
! CALL dynifc2blochf(nmodes, rws, nrws, xxq, uf, w2)
CALL dynifc2blochf(nmodes, rws, nrws, xxq_r, uf, w2, mirror_q)
!!!!!
ENDIF
!
! ...then take into account the mass factors and square-root the frequencies...
@ -1367,6 +1360,23 @@
!
xkk = xkf(:, ikk)
xkq2 = xkk + xxq
!
IF (plrn .and. time_rev_U_plrn) THEN
mirror_k = is_mirror_k(ik)
ELSE
mirror_k = .false.
END IF
! note that ikq_all return the global index of k+q
! Since ikq2 is global index, we need is_mirror_q instead of is_mirror_k
! is_mirror_q uses global k/q index.
! this is different from is_mirror_k which needs local index k
IF (plrn .and. time_rev_U_plrn) THEN
! kpg_map return global index, thus xkf_all is needed
mirror_kpq = is_mirror_q(ikq_all(ik, iq))
ELSE
mirror_kpq = .false.
END IF
!!!!!
!
CALL DGEMV('t', 3, nrr_k, twopi, irvec_r, 3, xkk, 1, 0.0_DP, rdotk, 1)
CALL DGEMV('t', 3, nrr_k, twopi, irvec_r, 3, xkq2, 1, 0.0_DP, rdotk2, 1)
@ -1393,12 +1403,22 @@
!
! Kohn-Sham first, then get the rotation matricies for following interp.
IF (eig_read) THEN
CALL hamwan2bloch(nbndsub, nrr_k, cufkk, etf_ks(:, ikk), chw_ks, cfac, dims)
CALL hamwan2bloch(nbndsub, nrr_k, cufkq, etf_ks(:, ikq), chw_ks, cfacq, dims)
!!!!!
!CALL hamwan2bloch(nbndsub, nrr_k, cufkk, etf_ks(:, ikk), chw_ks, cfac, dims)
!CALL hamwan2bloch(nbndsub, nrr_k, cufkq, etf_ks(:, ikq), chw_ks, cfacq, dims)
CALL hamwan2bloch(nbndsub, nrr_k, cufkk, etf_ks(:, ikk), chw_ks, cfac, dims, mirror_k)
CALL hamwan2bloch(nbndsub, nrr_k, cufkq, etf_ks(:, ikq), chw_ks, cfacq, dims, mirror_kpq)
!!!!!
ENDIF
!
CALL hamwan2bloch(nbndsub, nrr_k, cufkk, etf(:, ikk), chw, cfac, dims)
CALL hamwan2bloch(nbndsub, nrr_k, cufkq, etf(:, ikq), chw, cfacq, dims)
!!!!!
!
!CALL hamwan2bloch(nbndsub, nrr_k, cufkk, etf(:, ikk), chw, cfac, dims)
!CALL hamwan2bloch(nbndsub, nrr_k, cufkq, etf(:, ikq), chw, cfacq, dims)
CALL hamwan2bloch(nbndsub, nrr_k, cufkk, etf(:, ikk), chw, cfac, dims, mirror_k)
CALL hamwan2bloch(nbndsub, nrr_k, cufkq, etf(:, ikq), chw, cfacq, dims, mirror_kpq)
!!!!!
! Save U^{\dagger}_{mn}(k) matrix to file
!
! Apply a possible scissor shift
etf(icbm:nbndsub, ikk) = etf(icbm:nbndsub, ikk) + scissor
@ -1437,7 +1457,7 @@
! within a Fermi shell of size fsthick
!
IF (((MINVAL(ABS(etf(:, ikk) - ef)) < fsthick) .AND. &
(MINVAL(ABS(etf(:, ikq) - ef)) < fsthick)) .OR. wfcelec) THEN
(MINVAL(ABS(etf(:, ikq) - ef)) < fsthick))) THEN
!
! Compute velocities
!
@ -1581,16 +1601,10 @@
WRITE(stdout, '(7x, a, f12.6, a)' ) 'Adaptative smearing el-impurity = Min: ', DSQRT(2.d0) * MINVAL(valmin) * ryd2mev,' meV'
WRITE(stdout, '(7x, a, f12.6, a)' ) ' Max: ', DSQRT(2.d0) * MAXVAL(valmax) * ryd2mev,' meV'
ENDIF
!
! --------------------------------------------------------------------------------------------------
! Added by Chao Lian for polaron calculations
! Shell implementation for future use.
IF (wfcelec .AND. (.NOT. polaron_bq) .AND. (.NOT. polaron_interpol)) THEN
ufall(1:nmodes, 1:nmodes, iq) = uf(1:nmodes, 1:nmodes)
epfall(1:nbndfst, 1:nbndfst, 1:nmodes, 1:nkf, iq) = epf17(1:nbndfst, 1:nbndfst, 1:nmodes, 1:nkf)
ENDIF
! --------------------------------------------------------------------------------------------------
!
!!!!!
! Added by Chao Lian to save the el-ph matrix element to files
IF (plrn) CALL plrn_save_g_to_file(iq, epf17, wf)
!!!!!
IF (prtgkk ) CALL print_gkk(iq)
IF (phonselfen) CALL selfen_phon_q(iqq, iq, totq)
IF (elecselfen) CALL selfen_elec_q(iqq, iq, totq, first_cycle)
@ -1745,64 +1759,10 @@
ENDIF ! scatread
ENDDO ! end loop over q points
!
! --------------------------------------------------------------------------------
! Added for polaron calculations. Originally by Danny Sio, modified by Chao Lian.
! Shell implementation for future use.
IF (wfcelec .AND. (.NOT. polaron_bq) .AND. (.NOT. polaron_interpol)) THEN
IF (wfcelec_old) then
ALLOCATE(g2_4(ibndmax - ibndmin + 1, ibndmax - ibndmin + 1, nmodes, nkqtotf / 2), STAT = ierr)
IF (ierr /= 0) CALL errore('ephwann_shuffle', 'Error allocating g2_4', 1)
g2_4(:, :, :, :) = czero
CALL wfc_elec_old(nrr_k, nrr_q, nrr_g, irvec_q, irvec_g, &
ndegen_k, ndegen_q, ndegen_g, w2, uf, epmatwef, irvec_r, &
dims, dims2)
ELSE
CALL wfc_elec(nrr_k, ndegen_k, irvec_r, dims)
ENDIF
IF (polaron_wf) THEN ! calculating A(Re) from Ac.txt
CALL compute_a_re (iq, nrr_k, ndegen_k, irvec_r, dims)
RETURN
ENDIF
IF (polaron_interpol) THEN ! interpolate Ak from ar.txt ( A(Re))
CALL interpol_a_k(iq, nrr_k, ndegen_k, irvec_r, dims)
return
ENDIF
DO iq = iq_restart, nqf
xxq = xqf(:, iq)
IF (.NOT. lifc) THEN
CALL dynwan2bloch(nmodes, nrr_q, irvec_q, ndegen_q, xxq, uf, w2)
ELSE
CALL dynifc2blochf(nmodes, rws, nrws, xxq, uf, w2)
ENDIF
!
DO nu = 1, nmodes
!
! wf are the interpolated eigenfrequencies (omega on fine grid)
IF (w2(nu) > 0.d0) THEN
wf(nu, iq) = SQRT(ABS(w2(nu)))
ELSE
wf(nu, iq) = -SQRT(ABS(w2(nu)))
ENDIF
ENDDO
ENDDO
!
IF (polaron_bq) THEN ! interpolate bq from both A(Re) and Ac(k)
DO iq = 1, nqf
CALL interpol_bq(iq, nrr_k, nrr_q, nrr_g, irvec_q, irvec_g, ndegen_k, ndegen_q, ndegen_g, &
w2, uf, epmatwef, irvec_r, dims, dims2)
ENDDO
RETURN
ENDIF
DEALLOCATE(g2_4, STAT = ierr)
IF (ierr /= 0) CALL errore('ephwann_shuffle', 'Error deallocating g2_4', 1)
DEALLOCATE(epfall, STAT = ierr)
IF (ierr /= 0) CALL errore('ephwann_shuffle', 'Error deallocating epfall', 1)
DEALLOCATE(Hamil, STAT = ierr)
IF (ierr /= 0) CALL errore('ephwann_shuffle', 'Error deallocating Hamil', 1)
DEALLOCATE(eigVec, STAT = ierr)
IF (ierr /= 0) CALL errore('ephwann_shuffle', 'Error deallocating eigVec', 1)
ENDIF
! End Polaron Code
!!!!!
! Added for polaron calculations. Originally by Danny Sio, rewritten by Chao Lian.
IF (plrn) CALL plrn_flow_select(nrr_k, ndegen_k, irvec_r, nrr_q, ndegen_q, irvec_q, rws, nrws, dims)
!!!!!
! --------------------------------------------------------------------------------
!
! Check Memory usage

11
EPW/src/epw_init.f90
View File

@ -57,7 +57,6 @@
! --------------------------------------------------------------------------------
! Added for polaron calculations. Originally by Danny Sio, modified by Chao Lian.
! Shell implementation for future use.
USE epwcom, ONLY : polaron_wf
USE grid, ONLY : loadqmesh_serial, loadkmesh_para
! --------------------------------------------------------------------------------
!
@ -251,16 +250,6 @@
CALL dvanqq2()
ENDIF
!
! -------------------------------------------------------------------------------
! Added for polaron calculations. Originally by Danny Sio, modified by Chao Lian.
! Shell implementation for future use.
! IF (polaron_wf) then
! CALL loadqmesh_serial
! CALL loadkmesh_para
! CALL KSstate_extract()
! STOP
! ENDIF
! -------------------------------------------------------------------------------
!
CALL stop_clock('epw_init')
!

176
EPW/src/epw_readin.f90
View File

@ -93,20 +93,23 @@
USE paw_variables, ONLY : okpaw
USE io_epw, ONLY : param_get_range_vector
USE open_close_input_file, ONLY : open_input_file, close_input_file
!
! ---------------------------------------------------------------------------------------
! Added for polaron calculations. Originally by Danny Sio, modified by Chao Lian.
! Shell implementation for future use.
USE epwcom, ONLY : wfcelec, restart_polaron, spherical_cutoff, &
model_vertex, conv_thr_polaron, n_dop, &
polaron_wf, r01, r02, r03, num_cbands, start_band, &
start_mode, cb_shift, polaron_interpol, polaron_bq, &
polaron_dos, electron_dos, phonon_dos, diag_mode, &
restart_polaron_mode, polaron_type, nPlrn, wfcelec_old, &
sigma_plrn, ethr_Plrn, full_diagon_plrn, mixing_Plrn, &
init_plrn_wf, niterPlrn, nDOS_plrn, emax_plrn, emin_plrn, &
sigma_edos_plrn, sigma_pdos_plrn, pmax_plrn, pmin_plrn
!!!!!
! Added for polaron calculations by Chao Lian
USE epwcom, ONLY : plrn, restart_plrn, conv_thr_plrn, end_band_plrn, &
cal_psir_plrn, start_band_plrn, type_plrn, nstate_plrn, &
interp_Ank_plrn, interp_Bqu_plrn, &
init_sigma_plrn, init_k0_plrn, &
full_diagon_plrn, mixing_Plrn, init_plrn, niter_plrn, &
nDOS_plrn, edos_max_plrn, edos_min_plrn, edos_sigma_plrn, &
pdos_sigma_plrn, pdos_max_plrn, pdos_min_plrn, &
seed_plrn, ethrdg_plrn, time_rev_A_plrn, nhblock_plrn, &
beta_plrn, Mmn_plrn, recal_Mmn_plrn, r0_plrn, debug_plrn, &
time_rev_U_plrn, g_start_band_plrn, g_end_band_plrn, &
g_start_energy_plrn, g_end_energy_plrn, lrot, &
model_vertex_plrn, model_enband_plrn, model_phfreq_plrn, &
kappa_plrn, omega_LO_plrn, m_eff_plrn, step_wf_grid_plrn, &
g_power_order_plrn, g_tol_plrn, io_lvl_plrn, &
scell_mat_plrn, scell_mat, init_ntau_plrn, &
adapt_ethrdg_plrn, init_ethrdg_plrn, nethrdg_plrn
!-------------------------------------------------------------------------------------
! SH: Added for tc linearized equation, sparce sampling, and full-bandwidth calculations
USE epwcom, ONLY : gridsamp, griddens, tc_linear, tc_linear_solver, fbw, &
@ -201,16 +204,21 @@
!!!!!
!---------------------------------------------------------------------------------
! Added for polaron calculations. Originally by Danny Sio, modified by Chao Lian.
! Shell implementation for future use.
wfcelec, restart_polaron, spherical_cutoff, model_vertex, start_mode, &
conv_thr_polaron, polaron_wf, r01, r02, r03, num_cbands, start_band, &
cb_shift, polaron_interpol, polaron_bq, polaron_dos, electron_dos , &
phonon_dos, diag_mode, restart_polaron_mode, polaron_type, &
niterPlrn, wfcelec_old, sigma_plrn, ethr_Plrn, full_diagon_plrn, &
mixing_Plrn, init_plrn_wf, nPlrn, nDOS_plrn, emax_plrn, emin_plrn, &
!!!!!
! sigma_edos_plrn, sigma_pdos_plrn, pmax_plrn, pmin_plrn
sigma_edos_plrn, sigma_pdos_plrn, pmax_plrn, pmin_plrn, &
plrn, restart_plrn, conv_thr_plrn, end_band_plrn, lrot, &
cal_psir_plrn, start_band_plrn, type_plrn, nstate_plrn, &
interp_Ank_plrn, interp_Bqu_plrn, init_sigma_plrn, init_k0_plrn, &
full_diagon_plrn, mixing_Plrn, init_plrn, niter_plrn, &
nDOS_plrn, edos_max_plrn, edos_min_plrn, edos_sigma_plrn, &
pdos_sigma_plrn, pdos_max_plrn, pdos_min_plrn, &
seed_plrn, ethrdg_plrn, time_rev_A_plrn, nhblock_plrn, &
beta_plrn, Mmn_plrn, recal_Mmn_plrn, r0_plrn, debug_plrn, &
time_rev_U_plrn, g_start_band_plrn, g_end_band_plrn, &
g_start_energy_plrn, g_end_energy_plrn, &
model_vertex_plrn, model_enband_plrn, model_phfreq_plrn, &
kappa_plrn, omega_LO_plrn, m_eff_plrn, step_wf_grid_plrn, &
g_power_order_plrn, g_tol_plrn, io_lvl_plrn, &
scell_mat_plrn, scell_mat, init_ntau_plrn, &
adapt_ethrdg_plrn, init_ethrdg_plrn, nethrdg_plrn, &
!---------------------------------------------------------------------------------
! SH: Added for tc linearized equation, sparce sampling, and full-bandwidth runs
tc_linear, tc_linear_solver, gridsamp, griddens, fbw, dos_del, muchem
@ -629,6 +637,7 @@
meff = 1.d0
epsiheg = 1.d0
lphase = .FALSE.
lrot = .FALSE.
omegamin = 0.d0 ! eV
omegamax = 10.d0 ! eV
omegastep = 1.d0 ! eV
@ -663,47 +672,63 @@
ii_eps0 = 0.0d0
!!!!!
!
! --------------------------------------------------------------------------------
! Added for polaron calculations. Originally by Danny Sio, modified by Chao Lian.
! Shell implementation for future use.
nPlrn = 1
niterPlrn = 50
n_dop = 0.d0
smear_rpa = 1.d0
wfcelec = .false.
wfcelec_old = .false.
restart_polaron = .false.
spherical_cutoff = 100.d0
model_vertex = .false.
conv_thr_polaron = 1E-5
polaron_wf = .false.
polaron_interpol = .false.
polaron_bq = .false.
polaron_dos = .false.
r01 = 0.d0
r02 = 0.d0
r03 = 0.d0
num_cbands = 1 !2
start_band = 4 !11
start_mode = 1 !1
cb_shift = 0 !0
electron_dos = .false.
phonon_dos = .false.
! Added for polaron calculations by Chao Lian
nstate_plrn = 1
niter_plrn = 50
plrn = .false.
restart_plrn = .false.
model_vertex_plrn = .false.
model_enband_plrn = .false.
model_phfreq_plrn = .false.
kappa_plrn = 0.0
omega_LO_plrn = 0.0
m_eff_plrn = 0.0
conv_thr_plrn = 1E-5
g_power_order_plrn = 1
step_wf_grid_plrn = 1
cal_psir_plrn = .false.
interp_Ank_plrn = .false.
interp_Bqu_plrn = .false.
start_band_plrn = 0
end_band_plrn = 0
g_start_band_plrn = 0
g_end_band_plrn = 0
g_start_energy_plrn = -10.0
g_end_energy_plrn = 10.0
full_diagon_plrn = .false.
mixing_Plrn = 1.0
diag_mode = 1
init_plrn_wf = 2
init_plrn = 1
Mmn_plrn = .false.
recal_Mmn_plrn = .false.
debug_plrn = .false.
r0_plrn = zero
nDOS_plrn = 1000
emin_plrn = zero
pmin_plrn = zero
emax_plrn = 1.d0
pmax_plrn = 1d-2
sigma_edos_plrn = 0.1d0
sigma_pdos_plrn = 1d-3
restart_polaron_mode = 1
polaron_type = -1
sigma_plrn = 4.6
ethr_Plrn = 1E-3
edos_min_plrn = zero ! eV
pdos_min_plrn = zero ! meV
edos_max_plrn = zero ! eV
pdos_max_plrn = zero ! meV
edos_sigma_plrn = 0.01d0 ! eV
pdos_sigma_plrn = 0.1 ! meV
type_plrn = -1
init_sigma_plrn = 4.6
init_k0_plrn = (/1000.d0, 1000.d0, 1000.d0/)
ethrdg_plrn = 1E-6
time_rev_A_plrn = .false.
time_rev_U_plrn = .false.
nhblock_plrn = 1
beta_plrn = 0.0
g_tol_plrn = -0.01
io_lvl_plrn = 0
scell_mat_plrn = .false.
scell_mat(1, 1:3) = (/1, 0, 0/)
scell_mat(2, 1:3) = (/0, 1, 0/)
scell_mat(3, 1:3) = (/0, 0, 1/)
init_ntau_plrn = 1
adapt_ethrdg_plrn = .false.
init_ethrdg_plrn = 1.d-2
nethrdg_plrn = 11
! ---------------------------------------------------------------------------------
!
! Reading the namelist inputepw and check
@ -1180,26 +1205,13 @@
CALL mp_bcast(nk2, meta_ionode_id, world_comm)
CALL mp_bcast(nk3, meta_ionode_id, world_comm)
!
! ---------------------------------------------------------------------------------
! Added for polaron calculations. Originally by Danny Sio, modified by Chao Lian.
! Shell implementation for future use.
CALL mp_bcast(nPlrn, meta_ionode_id, world_comm)
CALL mp_bcast(niterPlrn, meta_ionode_id, world_comm)
CALL mp_bcast(spherical_cutoff, meta_ionode_id, world_comm)
CALL mp_bcast(conv_thr_polaron, meta_ionode_id, world_comm)
CALL mp_bcast(restart_polaron, meta_ionode_id, world_comm)
CALL mp_bcast(polaron_type, meta_ionode_id, world_comm)
CALL mp_bcast(sigma_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(full_diagon_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(mixing_Plrn, meta_ionode_id, world_comm)
CALL mp_bcast(ethr_Plrn, meta_ionode_id, world_comm)
CALL mp_bcast(init_plrn_wf, meta_ionode_id, world_comm)
CALL mp_bcast(sigma_edos_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(sigma_pdos_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(pmax_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(pmin_plrn, meta_ionode_id, world_comm)
! ---------------------------------------------------------------------------------
!
CALL mp_bcast(g_tol_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(io_lvl_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(init_ntau_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(init_k0_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(adapt_ethrdg_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(init_ethrdg_plrn, meta_ionode_id, world_comm)
CALL mp_bcast(nethrdg_plrn, meta_ionode_id, world_comm)
amass = AMU_RY * amass
!
!-----------------------------------------------------------------------

19
EPW/src/epw_setup.f90
View File

@ -43,11 +43,6 @@
USE fft_base, ONLY : dfftp
USE gvecs, ONLY : doublegrid
USE noncollin_module, ONLY : noncolin, domag, m_loc, angle1, angle2, ux, nspin_mag
! ---------------------------------------------------------------------------------
! Added for polaron calculations. Originally by Danny Sio, modified by Chao Lian.
! Shell implementation for future use.
USE epwcom, ONLY: polaron_wf
! ---------------------------------------------------------------------------------
!
IMPLICIT NONE
!
@ -68,20 +63,6 @@
!
CALL start_clock('epw_setup')
!
IF (.NOT. polaron_wf) THEN
DO jk = 1, nkstot
xx_c = xk_cryst(1, jk) * nkc1
yy_c = xk_cryst(2, jk) * nkc2
zz_c = xk_cryst(3, jk) * nkc3
!
! check that the k-mesh was defined in the positive region of 1st BZ
!
IF (xx_c < -eps5 .OR. yy_c < -eps5 .OR. zz_c < -eps5) &
CALL errore('epw_setup', 'coarse k-mesh needs to be strictly positive in 1st BZ', 1)
!
ENDDO
ENDIF ! not polaron_wf
!
! 1) Computes the total local potential (external+scf) on the smooth grid
!
CALL set_vrs(vrs, vltot, v%of_r, kedtau, v%kin_r, dfftp%nnr, nspin, doublegrid)

135
EPW/src/epwcom.f90
View File

@ -118,6 +118,8 @@
!! if .TRUE. print the |g| vertex in [meV].
LOGICAL :: lphase
!! if .TRUE. fix the gauge when diagonalizing the interpolated dynamical matrix and electronic Hamiltonian.
LOGICAL :: lrot
!! if .TRUE. fix the gauge when diagonalizing the interpolated dynamical matrix and electronic Hamiltonian.
LOGICAL :: lindabs
!! if .TRUE. perform phonon-assisted absorption calculations
LOGICAL :: use_ws
@ -422,68 +424,105 @@
! ----------------------------------------------------------------------------------
! Added for polaron calculations. Originally by Danny Sio, modified by Chao Lian.
! Shell implementation for future use.
INTEGER :: num_cbands
!! number of conduction bands accounted in the Hilbert space of polaron Hamilontian
INTEGER :: start_band
!! start band index in matrix element
INTEGER :: nPlrn
!! Number of polaron bands
INTEGER :: nDOS_plrn
INTEGER :: start_band_plrn, end_band_plrn
!! Start and end band index in matrix element
INTEGER :: nstate_plrn
!! Number of polaron states calculated
INTEGER :: ndos_plrn
!! Number of grid in polaron DOS calculation
INTEGER :: start_mode
!! start mode index
INTEGER :: cb_shift
!! CB shifted in polaron calculation
INTEGER :: diag_mode
!! diagonalization mode for polaron solver
INTEGER :: restart_polaron_mode
!! polaron restart mode
INTEGER :: polaron_type
INTEGER :: type_plrn
!! polaron type (electron/hole)
INTEGER :: init_plrn_wf
INTEGER :: init_plrn
!! initial polaron wavefuntion with 1:Gaussian package and 2:Random number
INTEGER :: niterPlrn
INTEGER :: niter_plrn
!! Maximum number of polaron SCF loops.
REAL(KIND = DP) :: spherical_cutoff
!! spherical_cutoff for fast convergence in polaron calculation
REAL(KIND = DP) :: conv_thr_polaron
INTEGER :: seed_plrn
!! The seed number to generate the random initial polaron wavefunction
INTEGER :: nhblock_plrn
!! The seed number to generate the random initial polaron wavefunction
INTEGER :: g_start_band_plrn, g_end_band_plrn
!! Start and end band in saving g matrix
INTEGER :: step_wf_grid_plrn
!! number of grid to skip in output of real space wavefunction
INTEGER :: io_lvl_plrn
!! number of grid to skip in output of real space wavefunction
INTEGER :: scell_mat(3,3)
!! Supercell transformation matrix for polaron
INTEGER :: init_ntau_plrn
!! Number of displacements to be read from file to calculate polaron energy landscape
INTEGER :: nethrdg_plrn
!! Number of steps in the diagonalization if adaptive_ethrdg_plrn=.true.
REAL(KIND = DP) :: conv_thr_plrn
!! convergent threshold for polaron calculation
REAL(KIND = DP) :: r01, r02, r03
!! x,y,z Carsteian coordinate of polaron centre
REAL(KIND = DP) :: emin_plrn, emax_plrn, sigma_edos_plrn
REAL(KIND = DP) :: edos_min_plrn, edos_max_plrn, edos_sigma_plrn
!! Electron Energy range in polaron DOS calculation
REAL(KIND = DP) :: pmin_plrn, pmax_plrn, sigma_pdos_plrn
REAL(KIND = DP) :: pdos_min_plrn, pdos_max_plrn, pdos_sigma_plrn
!! Phonon Energy range in polaron DOS calculation
REAL(KIND = DP) :: n_dop
!! extra added charge per cell (as tot_charge, with opposite sign)
REAL(KIND = DP) :: sigma_plrn
!! decay radius of polaron wavefunction in initialization
REAL(KIND = DP) :: ethr_Plrn
!! decay radius of polaron wavefunction in initialization
REAL(KIND = DP) :: mixing_Plrn
REAL(KIND = DP) :: init_sigma_plrn
!! decay radius of polaron wavefunction in Gaussian initialization
REAL(KIND = DP) :: init_k0_plrn(3)
!! center k of Gaussian initialization
REAL(KIND = DP) :: mixing_plrn
!! mixing weight in plrn iteration
REAL(KIND = DP) :: ethrdg_plrn
!! threshold in diagonalization of the polaron Hamiltonian
REAL(KIND = DP) :: beta_plrn(3)
!! Mixing weight in Self-consistency
LOGICAL :: wfcelec
REAL(KIND = DP) :: r0_plrn(3)
!! Mixing weight in Self-consistency
REAL(KIND = DP) :: g_start_energy_plrn
!! Mixing weight in Self-consistency
REAL(KIND = DP) :: g_end_energy_plrn
!! Mixing weight in Self-consistency
REAL(KIND = DP) :: kappa_plrn
!! Mixing weight in Self-consistency
REAL(KIND = DP) :: g_power_order_plrn
!! kappa
REAL(KIND = DP) :: omega_LO_plrn
!! Higher optical phonon frequency
REAL(KIND = DP) :: m_eff_plrn
!! effective mass in LP model
REAL(KIND = DP) :: g_tol_plrn
!!
REAL(KIND = DP) :: as(3,3)
!! Supercell lattice vectors transformed by Smat, as=S*at
REAL(KIND = DP) :: bs(3,3)
!! Supercell reciprocal lattice vectors transformed by Sbar=transpose(inverse(Smat), bs=Sbar*bg
REAL(KIND = DP) :: init_ethrdg_plrn
!! Initial coarse threshold to be used if adaptive_ethrdg_plrn=.true.
!!
LOGICAL :: plrn
!! if .true. calculates perturbated part of the wavefunction
LOGICAL :: restart_polaron
LOGICAL :: model_vertex_plrn
!! if .true. calculates vertex model
LOGICAL :: model_enband_plrn
!! if .true. calculates vertex model
LOGICAL :: model_phfreq_plrn
!! if .true. calculates vertex model
LOGICAL :: restart_plrn
!! if .true. Using Ack from written outputs
LOGICAL :: model_vertex
!! if .true. Using model vertex and effective mass for polaron calculation
LOGICAL :: wfcelec_old
!! if .true. Using old algorithm by DS
LOGICAL :: full_diagon_plrn
!! if .true. diagonalizing the polaron Hamiltonian with direct diagonalization
LOGICAL :: polaron_wf
!! if .true. diagonalizing the polaron Hamiltonian with direct diagonalization
LOGICAL :: cal_psir_plrn
!! if .true. Generating a 3D-plot for polaron wavefunction
LOGICAL :: polaron_interpol
LOGICAL :: interp_Ank_plrn
!! if .true. interpolating polaron A(k) from A(Re)
LOGICAL :: polaron_bq
LOGICAL :: interp_Bqu_plrn
!! if .true. interpolating polaron phonon component bq from A(Re) and An(k)
LOGICAL :: polaron_dos
!! if .true. calculating polaron dos
LOGICAL :: electron_dos
!! if .true. calculating contributed electron dos
LOGICAL :: phonon_dos
!! if .true. calculating excited phonon dos
LOGICAL :: debug_plrn
!! if .true. interpolating polaron phonon component bq from A(Re) and An(k)
LOGICAL :: time_rev_A_plrn
!! if .true. time_rev_A_plrn
LOGICAL :: time_rev_U_plrn
!! if .true. time_rev_A_plrn
LOGICAL :: Mmn_plrn
!! if .true. time_rev_A_plrn
LOGICAL :: recal_Mmn_plrn
!! if .true. time_rev_A_plrn
LOGICAL :: scell_mat_plrn
!! if .true. activate supercell transformation for polaron
LOGICAL :: adapt_ethrdg_plrn
!! if .true. activate adaptive threshold in polaron davidson diagonalization
! ----------------------------------------------------------------------------------
!
!-----------------------------------------------------------------------

305
EPW/src/grid.f90
View File

@ -29,14 +29,15 @@
USE mp_world, ONLY : mpime
USE kinds, ONLY : DP
USE epwcom, ONLY : filkf, nkf1, nkf2, nkf3, iterative_bte, &
rand_k, rand_nk, mp_mesh_k, system_2d, eig_read, vme
rand_k, rand_nk, mp_mesh_k, system_2d, eig_read, vme, &
scell_mat_plrn, scell_mat, as, bs
USE elph2, ONLY : nkqtotf, nkqf, xkf, wkf, nkf, xkfd, deltaq, &
xkf_irr, wkf_irr, bztoibz, s_bztoibz
USE cell_base, ONLY : at, bg
USE symm_base, ONLY : s, t_rev, time_reversal, nsym
USE io_var, ONLY : iunkf
USE low_lvl, ONLY : init_random_seed
USE constants_epw, ONLY : eps4
USE io_var, ONLY : iunkf, iunRpscell, iunkgridscell
USE low_lvl, ONLY : init_random_seed, matinv3
USE constants_epw, ONLY : eps4, eps8
USE noncollin_module, ONLY : noncolin
# if defined(__MPI)
USE parallel_include, ONLY : MPI_INTEGER2
@ -68,10 +69,32 @@
!! rest from the division of nr of q-points over pools
INTEGER :: ierr
!! Error status
INTEGER :: iRp1, iRp2, iRp3, Rpmax, nRp
!! Number of unit cells within supercell
INTEGER :: Rp_crys_p(3)
!! Unit cell vectors in primitive crystal coordinates
INTEGER, ALLOCATABLE :: Rp(:, :)
!! List of unit cell vectors within supercell in primitive crystal coords
INTEGER :: iGs1, iGs2, iGs3, Gsmax, nGs
!! Number of supercell G-vectors within primitive reciprocal unit cell
INTEGER :: Gs_crys_s(3)
!! Supercell G-vectors in supercell reciprocal coordinates
REAL(KIND = DP) :: ap(3, 3), bp(3, 3)
!! Auxiliary definitions of real and reciprocal primitive cell vector matrix
REAL(KIND = DP), ALLOCATABLE :: xkf_(:, :)
!! coordinates k-points
REAL(KIND = DP), ALLOCATABLE :: wkf_(:)
!! weights k-points
REAL(KIND = DP) :: scell_mat_b(3, 3)
!! Reciprocal lattice transformation matrix
REAL(KIND = DP) :: p2s(3, 3), bs2p(3, 3)
!! Transformation matrix from primitive to supercell crystal coordinates
REAL(KIND = DP) :: Rp_crys_s(3)
!! Unit cell vectors in supercell crystal coordinates
REAL(KIND = DP) :: Gs_crys_p(3)
!! Supercell G-vectors in primitive crystal coordinates
REAL(KIND = DP), ALLOCATABLE :: Gs(:, :)
!! Supercell G-vectors within primitive reciprocal unit cell
!
IF (mpime == ionode_id) THEN
IF (filkf /= '') THEN ! load from file
@ -114,6 +137,144 @@
CALL cryst_to_cart(nkqtotf, xkf_, at, -1)
ENDIF
!
!JLB
ELSEIF (scell_mat_plrn) THEN
!
WRITE(stdout, '(a)') ' '
WRITE(stdout, '(a)') ' Supercell transformation activated (k), as=S*at'
WRITE(stdout, '(a,3i4)') ' S(1, 1:3): ', scell_mat(1, 1:3)
WRITE(stdout, '(a,3i4)') ' S(2, 1:3): ', scell_mat(2, 1:3)
WRITE(stdout, '(a,3i4)') ' S(3, 1:3): ', scell_mat(3, 1:3)
!
ap = TRANSPOSE(at)
as = MATMUL(scell_mat,ap)
!
WRITE(stdout, '(a)') ' Transformed lattice vectors (alat units):'
WRITE(stdout, '(a,3f12.6)') ' as(1, 1:3): ', as(1, 1:3)
WRITE(stdout, '(a,3f12.6)') ' as(2, 1:3): ', as(2, 1:3)
WRITE(stdout, '(a,3f12.6)') ' as(3, 1:3): ', as(3, 1:3)
!
scell_mat_b = matinv3(REAL(scell_mat, DP))
scell_mat_b = TRANSPOSE(scell_mat_b)
!
WRITE(stdout, '(a)') ' Reciprocal lattice transformation matrix, Sbar = (S^{-1})^{t}:'
WRITE(stdout, '(a,3f12.6)') ' Sbar(1, 1:3): ', scell_mat_b(1, 1:3)
WRITE(stdout, '(a,3f12.6)') ' Sbar(2, 1:3): ', scell_mat_b(2, 1:3)
WRITE(stdout, '(a,3f12.6)') ' Sbar(3, 1:3): ', scell_mat_b(3, 1:3)
!
bp = TRANSPOSE(bg)
bs = MATMUL(scell_mat_b, bp)
!
WRITE(stdout, '(a)') ' Transformed reciprocal lattice vectors (2pi/alat units):'
WRITE(stdout, '(a,3f12.6)') ' bs(1, 1:3): ', bs(1, 1:3)
WRITE(stdout, '(a,3f12.6)') ' bs(2, 1:3): ', bs(2, 1:3)
WRITE(stdout, '(a,3f12.6)') ' bs(3, 1:3): ', bs(3, 1:3)
WRITE(stdout, '(a)') ' '
!
! Define transformation matrix from primitive crystal coordinates
! to supercell crystal coordinates Rp_crys_s = ((a_s)^{T})^{-1} (a_p)^{T} Rp_crys_p
p2s = matinv3(TRANSPOSE(as))
p2s = MATMUL(p2s,TRANSPOSE(ap))
!
! Find how many unit cells are contained within the supercell
Rpmax = 5*MAXVAL(scell_mat) ! This should be large enough to find all
ALLOCATE(Rp(3, Rpmax**3), STAT = ierr)
IF (ierr /= 0) CALL errore('loadkmesh_para', 'Error allocating Rp', 1)
Rp = 0
nRp = 0
DO iRp1 = -Rpmax, Rpmax
DO iRp2 = -Rpmax, Rpmax
DO iRp3 = -Rpmax, Rpmax
Rp_crys_p = (/iRp1, iRp2, iRp3/)
Rp_crys_s = MATMUL(p2s, Rp_crys_p)
! Unit cell within supercell if Rp\in(0,1) in supercell crystal coordinates
IF (ALL(Rp_crys_s > -eps8) .AND. ALL(Rp_crys_s < 1.d0-eps8)) THEN
nRp = nRp + 1
Rp(1:3, nRp) = Rp_crys_p
END IF
END DO
END DO
END DO
WRITE(stdout, '(a, 3i6)') ' Number of unit cells within supercell:', nRp
!
! Write Rp-s in supercell to file
IF (mpime == ionode_id) THEN
OPEN(UNIT = iunRpscell, FILE = 'Rp.scell.plrn', ACTION = 'write')
WRITE(iunRpscell, *) nRp
DO iRp1 = 1, nRp
WRITE(iunRpscell, *) Rp(1:3, iRp1)
END DO
CLOSE(iunRpscell)
ENDIF
!
IF (ALLOCATED(Rp)) DEALLOCATE(Rp)
!
! Define transformation matrix from reciprocal supercell crystal coordinates
! to reciprocal primitive crystal coordinates
! Gs_crys_p = ((bp)^{T})^{-1} (bs)^{T} Gs_crys_s
bs2p = matinv3(TRANSPOSE(bp))
bs2p = MATMUL(bs2p, TRANSPOSE(bs))
!
! Find how many k-points lie within primitive reciprocal cell
Gsmax = Rpmax ! This should be large enough to find all
ALLOCATE(Gs(3, Gsmax**3), STAT = ierr)
IF (ierr /= 0) CALL errore('loadqmesh_serial', 'Error allocating Gs', 1)
Gs = 0.d0
nGs = 0
DO iGs1 = -Gsmax, Gsmax
DO iGs2 = -Gsmax, Gsmax
DO iGs3 = -Gsmax, Gsmax
Gs_crys_s = (/iGs1, iGs2, iGs3/)
Gs_crys_p = MATMUL(bs2p, Gs_crys_s)
! Gs within primitive reciprocal unit cell if Gs\in(0,1) in crys_p coords.
IF (ALL(Gs_crys_p > -eps8) .AND. ALL(Gs_crys_p < 1.d0-eps8)) THEN
nGs = nGs + 1
Gs(1:3, nGs) = Gs_crys_p
END IF
END DO
END DO
END DO
WRITE(stdout, '(a, 3i6)') ' Number of k-points needed:', nGs
!DO iGs1 = 1, nGs
! WRITE(stdout, '(3f12.6)') Gs(1:3, iGs1)
!END DO
!
! Write Gs-s within unit cell BZ (k-grid) to file
IF (mpime == ionode_id) THEN
OPEN(UNIT = iunkgridscell, FILE = 'kgrid.scell.plrn', ACTION = 'write')
WRITE(iunkgridscell, *) nGs
DO iGs1 = 1, nGs
WRITE(iunkgridscell, '(3f12.6)') Gs(1:3, iGs1)
END DO
CLOSE(iunkgridscell)
ENDIF
!
! Save list of needed k-points
nkqtotf = nGs
ALLOCATE(xkf_(3, 2 * nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('loadkmesh_para', 'Error allocating xkf_', 1)
ALLOCATE(wkf_(2 * nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('loadkmesh_para', 'Error allocating wkf_', 1)
!
DO ik = 1, nkqtotf
!
ikk = 2 * ik - 1
ikq = ikk + 1
!
xkf_(:, ikk) = Gs(1:3, ik)
wkf_(ikk) = 1.d0 ! weight not important for polaron
!
xkf_(:, ikq) = xkf_(:, ikk)
wkf_(ikq) = 0.d0
!
ENDDO
!
! redefine nkqtotf to include the k+q points
nkqtotf = 2 * nkqtotf
!
IF (ALLOCATED(Gs)) DEALLOCATE(Gs)
!
!JLB
ELSEIF ((nkf1 /= 0) .AND. (nkf2 /= 0) .AND. (nkf3 /= 0)) THEN ! generate grid
IF (mp_mesh_k) THEN
! get size of the mp_mesh in the irr wedge
@ -1716,16 +1877,18 @@
USE mp_global, ONLY : inter_pool_comm
USE mp, ONLY : mp_bcast
USE mp_world, ONLY : mpime
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE epwcom, ONLY : filqf, nqf1, nqf2, nqf3, &
rand_q, rand_nq, mp_mesh_q, system_2d, lscreen, &
plselfen, specfun_pl
plselfen, specfun_pl, &
scell_mat_plrn, scell_mat, as, bs
USE elph2, ONLY : xqf, wqf, nqtotf, nqf
USE cell_base, ONLY : at, bg
USE symm_base, ONLY : s, t_rev, time_reversal, nsym
USE io_var, ONLY : iunqf
USE low_lvl, ONLY : init_random_seed
USE constants_epw, ONLY : eps4
USE low_lvl, ONLY : init_random_seed, matinv3
USE constants_epw, ONLY : eps4, eps8
!
IMPLICIT NONE
!
@ -1741,6 +1904,28 @@
!! Status integer
INTEGER :: ierr
!! Error status
INTEGER :: iRp1, iRp2, iRp3, Rpmax, nRp
!! Number of unit cells within supercell
INTEGER :: Rp_crys_p(3)
!! Unit cell vectors in primitive crystal coordinates
INTEGER, ALLOCATABLE :: Rp(:, :)
!! List of unit cell vectors within supercell in primitive crystal coords
INTEGER :: iGs1, iGs2, iGs3, Gsmax, nGs
!! Number of supercell G-vectors within primitive reciprocal unit cell
INTEGER :: Gs_crys_s(3)
!! Supercell G-vectors in supercell reciprocal coordinates
REAL(KIND = DP) :: ap(3, 3), bp(3, 3)
!! Auxiliary definitions of real and reciprocal primitive cell vector matrix
REAL(KIND = DP) :: scell_mat_b(3, 3)
!! Reciprocal lattice transformation matrix
REAL(KIND = DP) :: p2s(3, 3), bs2p(3, 3)
!! Transformation matrix from primitive to supercell crystal coordinates
REAL(KIND = DP) :: Rp_crys_s(3)
!! Unit cell vectors in supercell crystal coordinates
REAL(KIND = DP) :: Gs_crys_p(3)
!! Supercell G-vectors in primitive crystal coordinates
REAL(KIND = DP), ALLOCATABLE :: Gs(:, :)
!! Supercell G-vectors within primitive reciprocal unit cell
!
IF (mpime == ionode_id) THEN
IF (filqf /= '') THEN ! load from file
@ -1771,6 +1956,112 @@
CALL cryst_to_cart(nqtotf, xqf, at, -1)
ENDIF
!
!JLB
ELSEIF (scell_mat_plrn) THEN
!
WRITE(stdout, '(a)') ' '
WRITE(stdout, '(a)') ' Supercell transformation activated (q), as=S*at'
WRITE(stdout, '(a,3i4)') ' S(1, 1:3): ', scell_mat(1, 1:3)
WRITE(stdout, '(a,3i4)') ' S(2, 1:3): ', scell_mat(2, 1:3)
WRITE(stdout, '(a,3i4)') ' S(3, 1:3): ', scell_mat(3, 1:3)
!
ap = TRANSPOSE(at)
as = MATMUL(scell_mat,ap)
!
WRITE(stdout, '(a)') ' Transformed lattice vectors (alat units):'
WRITE(stdout, '(a,3f12.6)') ' as(1, 1:3): ', as(1, 1:3)
WRITE(stdout, '(a,3f12.6)') ' as(2, 1:3): ', as(2, 1:3)
WRITE(stdout, '(a,3f12.6)') ' as(3, 1:3): ', as(3, 1:3)
!
scell_mat_b = matinv3(REAL(scell_mat, DP))
scell_mat_b = TRANSPOSE(scell_mat_b)
!
WRITE(stdout, '(a)') ' Reciprocal lattice transformation matrix, Sbar = (S^{-1})^{t}:'
WRITE(stdout, '(a,3f12.6)') ' Sbar(1, 1:3): ', scell_mat_b(1, 1:3)
WRITE(stdout, '(a,3f12.6)') ' Sbar(2, 1:3): ', scell_mat_b(2, 1:3)
WRITE(stdout, '(a,3f12.6)') ' Sbar(3, 1:3): ', scell_mat_b(3, 1:3)
!
bp = TRANSPOSE(bg)
bs = MATMUL(scell_mat_b, bp)
!
WRITE(stdout, '(a)') ' Transformed reciprocal lattice vectors (2pi/alat units):'
WRITE(stdout, '(a,3f12.6)') ' bs(1, 1:3): ', bs(1, 1:3)
WRITE(stdout, '(a,3f12.6)') ' bs(2, 1:3): ', bs(2, 1:3)
WRITE(stdout, '(a,3f12.6)') ' bs(3, 1:3): ', bs(3, 1:3)
WRITE(stdout, '(a)') ' '
!
! Define transformation matrix from primitive crystal coordinates
! to supercell crystal coordinates Rp_crys_s = ((a_s)^{T})^{-1} (a_p)^{T} Rp_crys_p
p2s = matinv3(TRANSPOSE(as))
p2s = MATMUL(p2s,TRANSPOSE(ap))
!
! Find how many unit cells are contained within the supercell
Rpmax = 5*MAXVAL(scell_mat) ! This should be large enough to find all
ALLOCATE(Rp(3, Rpmax**3), STAT = ierr)
IF (ierr /= 0) CALL errore('loadqmesh_serial', 'Error allocating Rp', 1)
Rp = 0
nRp = 0
DO iRp1 = -Rpmax, Rpmax
DO iRp2 = -Rpmax, Rpmax
DO iRp3 = -Rpmax, Rpmax
Rp_crys_p = (/iRp1, iRp2, iRp3/)
Rp_crys_s = MATMUL(p2s, Rp_crys_p)
! Unit cell within supercell if Rp\in(0,1) in supercell crystal coordinates
IF (ALL(Rp_crys_s > -eps8) .AND. ALL(Rp_crys_s < 1.d0-eps8)) THEN
nRp = nRp + 1
Rp(1:3, nRp) = Rp_crys_p
END IF
END DO
END DO
END DO
WRITE(stdout, '(a, 3i6)') ' Number of unit cells within supercell:', nRp
!
IF (ALLOCATED(Rp)) DEALLOCATE(Rp)
!
! Define transformation matrix from reciprocal supercell crystal coordinates
! to reciprocal primitive crystal coordinates
! Gs_crys_p = ((bp)^{T})^{-1} (bs)^{T} Gs_crys_s
bs2p = matinv3(TRANSPOSE(bp))
bs2p = MATMUL(bs2p, TRANSPOSE(bs))
!
! Find how many q-points lie within primitive reciprocal cell
Gsmax = Rpmax ! This should be large enough to find all
ALLOCATE(Gs(3, Gsmax**3), STAT = ierr)
IF (ierr /= 0) CALL errore('loadqmesh_serial', 'Error allocating Gs', 1)
Gs = 0.d0
nGs = 0
DO iGs1 = -Gsmax, Gsmax
DO iGs2 = -Gsmax, Gsmax
DO iGs3 = -Gsmax, Gsmax
Gs_crys_s = (/iGs1, iGs2, iGs3/)
Gs_crys_p = MATMUL(bs2p, Gs_crys_s)
! Gs within primitive reciprocal unit cell if Gs\in(0,1) in crys_p coords.
IF (ALL(Gs_crys_p > -eps8) .AND. ALL(Gs_crys_p < 1.d0-eps8)) THEN
nGs = nGs + 1
Gs(1:3, nGs) = Gs_crys_p
END IF
END DO
END DO
END DO
WRITE(stdout, '(a, 3i6)') ' Number of q-points needed:', nGs
!
! Save list of needed q-points
nqtotf = nGs
ALLOCATE(xqf(3, nqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('loadqmesh_serial', 'Error allocating xqf', 1)
ALLOCATE(wqf(nqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('loadqmesh_serial', 'Error allocating wqf', 1)
!
DO iq = 1, nqtotf
!
xqf(:, iq) = Gs(1:3, iq)
wqf(iq) = 1.d0 ! weight not important for polaron
!
ENDDO
!
IF (ALLOCATED(Gs)) DEALLOCATE(Gs)
!
!JLB
ELSEIF ((nqf1 /= 0) .AND. (nqf2 /= 0) .AND. (nqf3 /= 0)) THEN ! generate grid
IF (mp_mesh_q) THEN
IF (lscreen) CALL errore ('loadqmesh', 'If lscreen=.TRUE. do not use mp_mesh_q',1)

7
EPW/src/io_var.f90
View File

@ -44,6 +44,7 @@
iunsparset_merge, iunepmatcb_merge, iunsparseqcb_merge, &
iunsparseicb_merge, iunsparsejcb_merge, iunsparsetcb_merge, &
iunsparsekcb_merge, iunepmat_merge
PUBLIC :: iunRpscell, iunkgridscell, iunpsirscell
!
! Output of physically relevant quantities (60-100)
@ -176,6 +177,12 @@
! Miscellaneous (326-350)
INTEGER :: epwbib = 326 ! EPW bibliographic file.
!
! Output quantities related to polaron (350-400)
!JLB: All the other polaron I/O units should also be defined here for consistency
INTEGER :: iunRpscell = 351 ! Rp unit cell list within polaron supercell
INTEGER :: iunkgridscell = 352 ! Gs k-grid used for transformed supercell
INTEGER :: iunpsirscell = 353 ! Polaron wf in real space for transformed supercell
!
! Merging of files (400-450)
INTEGER :: iunepmat_merge = 400
INTEGER :: iunsparseq_merge = 401

87
EPW/src/low_lvl.f90
View File

@ -864,44 +864,55 @@
!
!-----------------------------------------------------------------------
FUNCTION matinv3(A) RESULT(B)
!-----------------------------------------------------------------------
!!
!! Performs a direct calculation of the inverse of a 3×3 matrix.
!!
USE kinds, ONLY : DP
USE constants_epw, ONLY : eps160
!
REAL(KIND = DP), INTENT(in) :: A(3, 3)
!! Matrix
!
! Local variable
REAL(KIND = DP) :: detinv
!! Inverse of the determinant
REAL(KIND = DP) :: B(3, 3)
!! Inverse matrix
!
! Calculate the inverse determinant of the matrix
detinv = 1 / (A(1, 1) * A(2, 2) * A(3, 3) - A(1, 1) * A(2, 3) * A(3, 2) &
- A(1, 2) * A(2, 1) * A(3, 3) + A(1, 2) * A(2, 3) * A(3, 1) &
+ A(1, 3) * A(2, 1) * A(3, 2) - A(1, 3) * A(2, 2) * A(3, 1))
!
IF (detinv < eps160) THEN
CALL errore('matinv3', 'Inverse does not exist ', 1)
ENDIF
!
! Calculate the inverse of the matrix
B(1, 1) = +detinv * (A(2, 2) * A(3, 3) - A(2, 3) * A(3, 2))
B(2, 1) = -detinv * (A(2, 1) * A(3, 3) - A(2, 3) * A(3, 1))
B(3, 1) = +detinv * (A(2, 1) * A(3, 2) - A(2, 2) * A(3, 1))
B(1, 2) = -detinv * (A(1, 2) * A(3, 3) - A(1, 3) * A(3, 2))
B(2, 2) = +detinv * (A(1, 1) * A(3, 3) - A(1, 3) * A(3, 1))
B(3, 2) = -detinv * (A(1, 1) * A(3, 2) - A(1, 2) * A(3, 1))
B(1, 3) = +detinv * (A(1, 2) * A(2, 3) - A(1, 3) * A(2, 2))
B(2, 3) = -detinv * (A(1, 1) * A(2, 3) - A(1, 3) * A(2, 1))
B(3, 3) = +detinv * (A(1, 1) * A(2, 2) - A(1, 2) * A(2, 1))
!-----------------------------------------------------------------------
END FUNCTION matinv3
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!!
!! Performs a direct calculation of the inverse of a 3×3 matrix.
!!
USE kinds, ONLY : DP
USE constants_epw, ONLY : eps160
!
REAL(KIND = DP), INTENT(in) :: A(3, 3)
!! Matrix
!
! Local variable
REAL(KIND = DP) :: detinv, det
!! Inverse of the determinant
REAL(KIND = DP) :: B(3, 3)
!! Inverse matrix
!
!! Calculate the inverse determinant of the matrix
!detinv = 1 / (A(1, 1) * A(2, 2) * A(3, 3) - A(1, 1) * A(2, 3) * A(3, 2) &
! - A(1, 2) * A(2, 1) * A(3, 3) + A(1, 2) * A(2, 3) * A(3, 1) &
! + A(1, 3) * A(2, 1) * A(3, 2) - A(1, 3) * A(2, 2) * A(3, 1))
!!
!IF (detinv < eps160) THEN
! CALL errore('matinv3', 'Inverse does not exist ', 1)
!ENDIF
!JLB
det = (A(1, 1) * A(2, 2) * A(3, 3) - A(1, 1) * A(2, 3) * A(3, 2) &
- A(1, 2) * A(2, 1) * A(3, 3) + A(1, 2) * A(2, 3) * A(3, 1) &
+ A(1, 3) * A(2, 1) * A(3, 2) - A(1, 3) * A(2, 2) * A(3, 1))
!
IF (ABS(det) < eps160) THEN
CALL errore('matinv3', 'Inverse does not exist ', 1)
END IF
!
detinv = 1 / det
!JLB
!
! Calculate the inverse of the matrix
B(1, 1) = +detinv * (A(2, 2) * A(3, 3) - A(2, 3) * A(3, 2))
B(2, 1) = -detinv * (A(2, 1) * A(3, 3) - A(2, 3) * A(3, 1))
B(3, 1) = +detinv * (A(2, 1) * A(3, 2) - A(2, 2) * A(3, 1))
B(1, 2) = -detinv * (A(1, 2) * A(3, 3) - A(1, 3) * A(3, 2))
B(2, 2) = +detinv * (A(1, 1) * A(3, 3) - A(1, 3) * A(3, 1))
B(3, 2) = -detinv * (A(1, 1) * A(3, 2) - A(1, 2) * A(3, 1))
B(1, 3) = +detinv * (A(1, 2) * A(2, 3) - A(1, 3) * A(2, 2))
B(2, 3) = -detinv * (A(1, 1) * A(2, 3) - A(1, 3) * A(2, 1))
B(3, 3) = +detinv * (A(1, 1) * A(2, 2) - A(1, 2) * A(2, 1))
!-----------------------------------------------------------------------
END FUNCTION matinv3
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
PURE FUNCTION find_minimum(grid, grid_dim) RESULT(minpos)
!-----------------------------------------------------------------------

3792
EPW/src/polaron.f90 Normal file
View File

File diff suppressed because it is too large Load Diff

8
EPW/src/pw2wan2epw.f90
View File

@ -3730,8 +3730,12 @@
OPEN(UNIT = iun_plot, FILE=TRIM(wancube), FORM='formatted', STATUS='unknown')
CALL date_and_tim(cdate, ctime )
! First two lines are comments
WRITE(iun_plot, *) ' Generated by EPW code'
WRITE(iun_plot, *) ' On ', cdate, ' at ', ctime
!!!!!!
!WRITE(iun_plot, *) ' Generated by EPW code'
!WRITE(iun_plot, *) ' On ', cdate, ' at ', ctime
WRITE(iun_plot, *) ' Generated by EPW code On ', cdate, ' at ', ctime
WRITE(iun_plot, '(9I4)') ngx, ngy, ngz, istart(1:3), iend(1:3)
!!!!!!
! Number of atoms, origin of cube (Cartesians) wrt simulation (home) cell
WRITE(iun_plot, '(i4,3f13.5)') icount, orig(1), orig(2), orig(3)
!

201
EPW/src/test_tools.f90
View File

@ -1,201 +0,0 @@
module test_tools
USE kinds, only : DP
INTERFACE para_write
MODULE PROCEDURE &
para_write_i, para_write_i1, para_write_i2, para_write_i3, para_write_i4, para_write_i5, para_write_i6, &
para_write_r, para_write_r1, para_write_r2, para_write_r3, para_write_r4, para_write_r5, para_write_r6,&
para_write_c, para_write_c1, para_write_c2, para_write_c3, para_write_c4, para_write_c5, para_write_c6
END INTERFACE
contains
SUBROUTINE para_write_i(A, filename)
IMPLICIT NONE
INTEGER, INTENT(IN) :: A
CHARACTER(LEN=*), INTENT(IN) :: filename
END SUBROUTINE
SUBROUTINE para_write_i1(A, varName)
USE mp_world, ONLY : mpime
IMPLICIT NONE
INTEGER, INTENT(IN) :: A(:)
CHARACTER(LEN=*), INTENT(IN) :: varName
CHARACTER(LEN=256) :: filename
CHARACTER(LEN=256) :: lineLen
INTEGER :: bounds(1), i, j
bounds = SHAPE(A)
WRITE(filename,'(a, i0)') varName, mpime
OPEN(UNIT=23747806, FILE='./test_out/'//filename, POSITION="append", FORM='formatted')
WRITE(23747806, *) REPEAT('-',10), bounds
WRITE(23747806, '(I5)') (A(j), j=1, bounds(1))
!WRITE(lineLen,'("(",i0,"E12.4)")') bounds(1) ! x2 due to complex
!do i = 1, bounds(2)
!WRITE(23747806, '(E12.4,"+",E12.4,"j")') (A(j,i), j=1, bounds(1))
!WRITE(23747806, '(E12.4,"+",E12.4,"j")', advance="no") A(:,i)
!WRITE(23747806, lineLen) (A(j,i), j=1, bounds(1))
!end do
CLOSE(23747806)
END SUBROUTINE
SUBROUTINE para_write_i2(A)
IMPLICIT NONE
INTEGER, INTENT(IN) :: A(:,:)
END SUBROUTINE
SUBROUTINE para_write_i3(A)
IMPLICIT NONE
INTEGER, INTENT(IN) :: A(:,:,:)
END SUBROUTINE
SUBROUTINE para_write_i4(A)
IMPLICIT NONE
INTEGER, INTENT(IN) :: A(:,:,:,:)
END SUBROUTINE
SUBROUTINE para_write_i5(A)
IMPLICIT NONE
INTEGER, INTENT(IN) :: A(:,:,:,:,:)
END SUBROUTINE
SUBROUTINE para_write_i6(A)
IMPLICIT NONE
INTEGER, INTENT(IN) :: A(:,:,:,:,:,:)
END SUBROUTINE
SUBROUTINE para_write_r(A)
IMPLICIT NONE
REAL(dp), INTENT(IN) :: A
END SUBROUTINE
SUBROUTINE para_write_r1(A, varName)
USE mp_world, ONLY : mpime
IMPLICIT NONE
REAL(dp), INTENT(IN) :: A(:)
CHARACTER(LEN=*), INTENT(IN) :: varName
CHARACTER(LEN=256) :: filename
CHARACTER(LEN=256) :: lineLen
INTEGER :: bounds(1), i, j
bounds = SHAPE(A)
WRITE(filename,'(a, i0)') varName, mpime
OPEN(UNIT=23747806, FILE='./test_out/'//filename, POSITION="append", FORM='formatted')
WRITE(23747806, *) REPEAT('-',10), bounds(1), 1
WRITE(23747806, '(E12.4)') (A(j), j=1, bounds(1))
!WRITE(lineLen,'("(",i0,"E12.4)")') bounds(1) ! x2 due to complex
!do i = 1, bounds(2)
!WRITE(23747806, '(E12.4,"+",E12.4,"j")') (A(j,i), j=1, bounds(1))
!WRITE(23747806, '(E12.4,"+",E12.4,"j")', advance="no") A(:,i)
!WRITE(23747806, lineLen) (A(j,i), j=1, bounds(1))
!end do
CLOSE(23747806)
END SUBROUTINE
SUBROUTINE para_write_r2(A, varName)
USE mp_world, ONLY : mpime
IMPLICIT NONE
REAL(dp), INTENT(IN) :: A(:,:)
CHARACTER(LEN=*), INTENT(IN) :: varName
CHARACTER(LEN=256) :: filename
CHARACTER(LEN=256) :: lineLen
INTEGER :: bounds(2), i, j
bounds = SHAPE(A)
WRITE(filename,'(a, i0)') varName, mpime
OPEN(UNIT=23747806, FILE='./test_out/'//filename, POSITION="append", FORM='formatted')
WRITE(23747806, *) REPEAT('-',10), bounds
WRITE(23747806, '(E12.4)') ((A(j,i), j=1, bounds(1)), i=1, bounds(2))
!WRITE(lineLen,'("(",i0,"E12.4)")') bounds(1) ! x2 due to complex
!do i = 1, bounds(2)
!WRITE(23747806, '(E12.4,"+",E12.4,"j")') (A(j,i), j=1, bounds(1))
!WRITE(23747806, '(E12.4,"+",E12.4,"j")', advance="no") A(:,i)
!WRITE(23747806, lineLen) (A(j,i), j=1, bounds(1))
!end do
CLOSE(23747806)
END SUBROUTINE
SUBROUTINE para_write_r3(A)
IMPLICIT NONE
REAL(dp), INTENT(IN) :: A(:,:,:)
END SUBROUTINE
SUBROUTINE para_write_r4(A)
IMPLICIT NONE
REAL(dp), INTENT(IN) :: A(:,:,:,:)
END SUBROUTINE
SUBROUTINE para_write_r5(A)
IMPLICIT NONE
REAL(dp), INTENT(IN) :: A(:,:,:,:,:)
END SUBROUTINE
SUBROUTINE para_write_r6(A)
IMPLICIT NONE
REAL(dp), INTENT(IN) :: A(:,:,:,:,:,:)
END SUBROUTINE
SUBROUTINE para_write_c(A)
IMPLICIT NONE
COMPLEX(dp), INTENT(IN) :: A
END SUBROUTINE
SUBROUTINE para_write_c1(A, varName)
USE mp_world, ONLY : mpime
IMPLICIT NONE
COMPLEX(dp), INTENT(IN) :: A(:)
CHARACTER(LEN=*), INTENT(IN) :: varName
CHARACTER(LEN=256) :: filename
CHARACTER(LEN=256) :: lineLen
INTEGER :: bounds(1), i, j
bounds = SHAPE(A)
WRITE(filename,'(a, i0)') varName, mpime
OPEN(UNIT=23747806, FILE='./test_out/'//filename, POSITION="append", FORM='formatted')
WRITE(23747806, *) REPEAT('-',10), bounds(1), 1
WRITE(23747806, '("(",E12.4,"+",E12.4,"j",")")') (A(j), j=1, bounds(1))
!WRITE(lineLen,'("(",i0,"E12.4)")') bounds(1) ! x2 due to complex
!do i = 1, bounds(2)
!WRITE(23747806, '(E12.4,"+",E12.4,"j")') (A(j,i), j=1, bounds(1))
!WRITE(23747806, '(E12.4,"+",E12.4,"j")', advance="no") A(:,i)
!WRITE(23747806, lineLen) (A(j,i), j=1, bounds(1))
!end do
CLOSE(23747806)
END SUBROUTINE
SUBROUTINE para_write_c2(A, varName)
USE mp_world, ONLY : mpime
IMPLICIT NONE
COMPLEX(dp), INTENT(IN) :: A(:,:)
CHARACTER(LEN=*), INTENT(IN) :: varName
CHARACTER(LEN=256) :: filename
CHARACTER(LEN=256) :: lineLen
INTEGER :: bounds(2), i, j
bounds = SHAPE(A)
WRITE(filename,'(a, i0)') varName, mpime
OPEN(UNIT=23747806, FILE='./test_out/'//filename, POSITION="append", FORM='formatted')
WRITE(23747806, *) REPEAT('-',10), bounds
WRITE(23747806, '("(",E12.4,"+",E12.4,"j",")")') ((A(j,i), j=1, bounds(1)), i=1, bounds(2))
!WRITE(lineLen,'("(",i0,"E12.4)")') bounds(1) ! x2 due to complex
!do i = 1, bounds(2)
!WRITE(23747806, '(E12.4,"+",E12.4,"j")') (A(j,i), j=1, bounds(1))
!WRITE(23747806, '(E12.4,"+",E12.4,"j")', advance="no") A(:,i)
!WRITE(23747806, lineLen) (A(j,i), j=1, bounds(1))
!end do
CLOSE(23747806)
END SUBROUTINE
SUBROUTINE para_write_c3(A)
IMPLICIT NONE
COMPLEX(dp), INTENT(IN) :: A(:,:,:)
END SUBROUTINE
SUBROUTINE para_write_c4(A)
IMPLICIT NONE
COMPLEX(dp), INTENT(IN) :: A(:,:,:,:)
END SUBROUTINE
SUBROUTINE para_write_c5(A)
IMPLICIT NONE
COMPLEX(dp), INTENT(IN) :: A(:,:,:,:,:)
END SUBROUTINE
SUBROUTINE para_write_c6(A)
IMPLICIT NONE
COMPLEX(dp), INTENT(IN) :: A(:,:,:,:,:,:)
END SUBROUTINE
end module

15
EPW/src/utilities.f90
View File

@ -1434,7 +1434,6 @@
USE pwcom, ONLY : ef
USE mp, ONLY : mp_max, mp_min
USE mp_global, ONLY : inter_pool_comm
USE epwcom, ONLY : wfcelec
USE constants_epw, ONLY : ryd2ev
!
IMPLICIT NONE
@ -1470,20 +1469,6 @@
!
ENDDO
ENDDO
IF (wfcelec) then
DO ik = 1, nkqf
DO ibnd = 1, nbndsub
ebnd = etf(ibnd, ik)
!
IF (ebnd < fsthick + ef .AND. ebnd > ef) THEN
ibndmin = MIN(ibnd, ibndmin)
ibndmax = MAX(ibnd, ibndmax)
ebndmin = MIN(ebnd, ebndmin)
ebndmax = MAX(ebnd, ebndmax)
ENDIF
ENDDO
ENDDO
ENDIF
!
tmp = DBLE(ibndmin)
CALL mp_min(tmp, inter_pool_comm)

350
EPW/src/wan2bloch.f90
View File

@ -17,7 +17,10 @@
CONTAINS
!
!--------------------------------------------------------------------------
SUBROUTINE hamwan2bloch(nbnd, nrr, cuf, eig, chw, cfac, dims)
!!!!!
!SUBROUTINE hamwan2bloch(nbnd, nrr, cuf, eig, chw, cfac, dims)
SUBROUTINE hamwan2bloch(nbnd, nrr, cuf, eig, chw, cfac, dims, is_mirror)
!!!!!
!--------------------------------------------------------------------------
!!
!! From the Hamiltonian in Wannier representation, find the corresponding
@ -31,7 +34,10 @@
!! output : rotation matrix cuf(nbnd, nbnd)
!! interpolated hamiltonian eigenvalues eig(nbnd)
!!
!! 2019: Weng Hong Sio and SP: Lifting of degeneracies.
!! 2021: CL : Replace the random perturbation matrix with prime number matrix
!! in Lifting of degeneracies; control tag : lphase
!! 2021: CL : Rotate the the largest element in eigenvector to real axis. (lrot)
!! 2019: Weng Hong Sio and SP: Lifting of degeneracies. control tag: lrot
!! P_prime = U^dag P U where P is a random perturbation matrix
!! cuf = (eigvector of P_prime) * U
!! P_prime spans the degenenrate subspace.
@ -43,6 +49,9 @@
USE constants_epw, ONLY : czero, cone, zero, one, eps12, eps16
USE epwcom, ONLY : use_ws
USE low_lvl, ONLY : utility_zdotu, degen_sort
!!!!!
USE epwcom, ONLY : debug_plrn, lphase, lrot
!!!!!
!
IMPLICIT NONE
!
@ -62,6 +71,10 @@
!! Hamiltonian in Wannier basis
COMPLEX(KIND = DP), INTENT(out) :: cuf(nbnd, nbnd)
!! Rotation matrix U^\dagger, fine mesh
!!!!!
LOGICAL, INTENT(IN), OPTIONAL :: is_mirror
!! .true. if k-point is a time-reversal invariant point
!!!!!
!
! Local variables
LOGICAL :: duplicates
@ -114,10 +127,22 @@
!! Perturbation matrix made of small complex random number on the full space
COMPLEX(KIND = DP), ALLOCATABLE :: cwork(:)
!! Complex work variable
COMPLEX(KIND = DP), ALLOCATABLE :: P_prime(:, :)
!!!!!
! COMPLEX(KIND = DP), ALLOCATABLE :: P_prime(:, :)
REAL(KIND = DP), ALLOCATABLE :: P_prime(:, :)
!!!!!
!! Perturbation matrix on the subspace
COMPLEX(KIND = DP), ALLOCATABLE :: Uk(:, :)
!! Rotation matrix on the full space
!!!!!
INTEGER :: ibnd_max(1)
!! Index of the maximum element
REAL(KIND = DP) :: norm_vec(nbnd)
!! Real Hamiltonian matrix in Bloch basis for TRI k points
COMPLEX(KIND = DP) :: cfac2(nrr, dims, dims)
!! cfac for the time-reversial point if pointed
COMPLEX(KIND = DP) :: fac_max
!!!!!
!
CALL start_clock('HamW2B')
!----------------------------------------------------------
@ -133,17 +158,31 @@
! H~(k'+q') is chf( nbnd, nbnd, 2*ik )
!
chf(:, :) = czero
!!!!!!!
! Calculate the -k if this is a mirror point
! Do nothing if not specified
cfac2 = cfac
IF (PRESENT(is_mirror)) THEN
IF(is_mirror) cfac2 = CONJG(cfac)
END IF
!!!!!!!
!
IF (use_ws) THEN
DO iw = 1, dims
DO iw2 = 1, dims
DO ir = 1, nrr
chf(iw, iw2) = chf(iw, iw2) + chw(iw, iw2, ir) * cfac(ir, iw, iw2)
!!!!!!!
!chf(iw, iw2) = chf(iw, iw2) + chw(iw, iw2, ir) * cfac(ir, iw, iw2)
chf(iw, iw2) = chf(iw, iw2) + chw(iw, iw2, ir) * cfac2(ir, iw, iw2)
!!!!!!!
ENDDO
ENDDO
ENDDO
ELSE
CALL ZGEMV('n', nbnd**2, nrr, cone, chw, nbnd**2, cfac(:, 1, 1), 1, cone, chf, 1)
!!!!!!!
!CALL ZGEMV('n', nbnd**2, nrr, cone, chw, nbnd**2, cfac(:, 1, 1), 1, cone, chf, 1)
CALL ZGEMV('n', nbnd**2, nrr, cone, chw, nbnd**2, cfac2(:, 1, 1), 1, cone, chf, 1)
!!!!!!!
ENDIF
!
!---------------------------------------------------------------------
@ -199,24 +238,27 @@
degen_group(2, ig) = degen_group(1, ig) + degen_group(2, ig) -1
ENDDO
!
!!!!!!!
! This old random matrix generation is removed
! Generate a pertubation matrix of size (nbnd x nbnd) made of random number
!CALL init_random_seed()
! SP: Using random_number does not work because the perturbation needs to be the
! same when calling hamwan2bloch at k and k+q (see ephwann_shuffle).
! Therefore I fix a "random" number 0.25644832 + 0.01 * ibnd and 0.11584272 + 0.025 * jbnd
P(:, :) = czero
DO ibnd = 1, nbnd
DO jbnd = 1, nbnd
!CALL random_number(rand1)
!CALL random_number(rand2)
rand1 = 0.25644832 + 0.01 * ibnd
rand2 = 0.11584272 + 0.025 * jbnd
P(jbnd, ibnd) = CMPLX(rand1, rand2, KIND = DP)
ENDDO
ENDDO
!P(:, :) = czero
!DO ibnd = 1, nbnd
! DO jbnd = 1, nbnd
! !CALL random_number(rand1)
! !CALL random_number(rand2)
! rand1 = 0.25644832 + 0.01 * ibnd
! rand2 = 0.11584272 + 0.025 * jbnd
! P(jbnd, ibnd) = CMPLX(rand1, rand2, KIND = DP)
! ENDDO
!dENDDO
!
! Hermitize the Perturbation matrix and make it small
P = 0.5d0 * (P + TRANSPOSE(CONJG(P))) * ABS(MINVAL(w)) * 0.1d0
!P = 0.5d0 * (P + TRANSPOSE(CONJG(P))) * ABS(MINVAL(w)) * 0.1d0
!!!!!!!
!
DO ig = 1, ndeg
starting = degen_group(1, ig)
@ -224,7 +266,10 @@
! Size of the degenerate subspace
length = ending - starting + 1
!
ALLOCATE(rwork(1 + 5 * length + 2 * length**2), STAT = ierr)
!!!!!!!
!ALLOCATE(rwork(1 + 5 * length + 2 * length**2), STAT = ierr)
ALLOCATE(rwork(length**2 + 2 * length), STAT = ierr)
!!!!!!!
IF (ierr /= 0) CALL errore('hamwan2bloch', 'Error allocating rwork', 1)
ALLOCATE(iwork(3 + 5 * length), STAT = ierr)
IF (ierr /= 0) CALL errore('hamwan2bloch', 'Error allocating iwork', 1)
@ -238,14 +283,22 @@
IF (ierr /= 0) CALL errore('hamwan2bloch', 'Error allocating wp', 1)
!
Uk(:, :) = cz(:, starting:ending)
P_prime = MATMUL(TRANSPOSE(CONJG(Uk)), MATMUL(P, Uk))
! Create a matrix filled with prime numbers
!!!!!!
! P_prime = MATMUL(TRANSPOSE(CONJG(Uk)), MATMUL(P, Uk))
CALL prime_number_matrix(P_prime, length)
!!!!!!
! Diagonalization of P_prime
CALL ZHEEVD('V', 'L', length, P_prime, length, wp, cwork, &
2 * length + length**2, rwork, 1 + 5 * length + 2 * length**2, &
iwork, 3 + 5 * length, info)
!!!!!!
! CALL ZHEEVD('V', 'L', length, P_prime, length, wp, cwork, &
! 2 * length + length**2, rwork, 1 + 5 * length + 2 * length**2, &
! iwork, 3 + 5 * length, info)
CALL DSYEV('V', 'L', length, P_prime, length, wp, rwork, &
length**2 + 2 * length, info)
!!!!!!
! On exiting P_prime is the eigenvector of the P_prime matrix and wp the eigenvector.
!
cz(:, starting:ending) = MATMUL(Uk, P_prime)
IF (lphase) cz(:, starting:ending) = MATMUL(Uk, P_prime)
!
DEALLOCATE(rwork, STAT = ierr)
IF (ierr /= 0) CALL errore('hamwan2bloch', 'Error deallocating rwork', 1)
@ -261,11 +314,23 @@
IF (ierr /= 0) CALL errore('hamwan2bloch', 'Error deallocating wp', 1)
ENDDO ! ig
!
!!!!!!
! Find the largest element and set it to pure real
IF(lrot) THEN
DO ibnd = 1, nbnd
DO jbnd = 1, nbnd
norm_vec(jbnd) = ABS(cz(jbnd, ibnd))
END DO
ibnd_max(:) = MAXLOC(norm_vec(1:nbnd))
fac_max = cz(ibnd_max(1), ibnd) / norm_vec(ibnd_max(1))
cz(1:nbnd, ibnd) = cz(1:nbnd, ibnd) * CONJG(fac_max)
END DO
END IF
!!!!!!
DO jbnd = 1, nbnd
INNER : DO ibnd = 1, nbnd
IF (ABS(cz(ibnd, jbnd)) > eps12) THEN
cz(:, jbnd) = cz(:, jbnd) * CONJG(cz(ibnd, jbnd))
!cz(:, jbnd) = cz(:, jbnd) / SQRT(zdotu(nbnd, CONJG(cz(:, jbnd)), 1, cz(:, jbnd), 1))
cz(:, jbnd) = cz(:, jbnd) / SQRT(utility_zdotu(CONJG(cz(:, jbnd)), cz(:, jbnd)))
EXIT INNER
ENDIF
@ -277,6 +342,12 @@
!
cuf = CONJG(TRANSPOSE(cz))
eig = w
!!!!!!
! Do the conjugate on eigenvector
IF (PRESENT(is_mirror)) THEN
IF(is_mirror) cuf = TRANSPOSE(cz)
END IF
!!!!!!
!
CALL stop_clock('HamW2B')
!
@ -284,8 +355,49 @@
END SUBROUTINE hamwan2bloch
!--------------------------------------------------------------------------
!
!!!!!!
!--------------------------------------------------------------------------
SUBROUTINE dynwan2bloch(nmodes, nrr_q, irvec_q, ndegen_q, xxq, cuf, eig)
SUBROUTINE prime_number_matrix(A, n)
!--------------------------------------------------------------------------
!!
!! Generating a n x n matrix A filled with prime numbers
!! For example, if n = 4, A =
!! |2 3 5 7|
!! |11 13 17 19|
!! |23 29 31 37|
!! |41 43 47 53|
!! Used to perturb the degenerate eigenstates
!! 2021 Chao Lian
!
USE kinds, ONLY : DP
IMPLICIT NONE
REAL(dp), INTENT(OUT) :: A(:, :)
INTEGER, INTENT(IN) :: n
INTEGER :: prime_numbers(60) = (/2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, &
& 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, &
& 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281/)
INTEGER :: i, j, k
k = 0
DO i = 1, n
DO j = i, n
k = k + 1
A(i, j) = REAL(prime_numbers(k), dp)
END DO
END DO
DO i = 1, n
DO j = 1, i
A(i, j) = (A(j, i))
END DO
END DO
END SUBROUTINE
!
!--------------------------------------------------------------------------
!SUBROUTINE dynwan2bloch(nmodes, nrr_q, irvec_q, ndegen_q, xxq, cuf, eig)
SUBROUTINE dynwan2bloch(nmodes, nrr_q, irvec_q, ndegen_q, xxq, cuf, eig, is_mirror)
!!!!!!
!--------------------------------------------------------------------------
!!
!! From the Hamiltonian in Wannier representation, find the corresponding
@ -296,16 +408,23 @@
!! required to obtain proper phonon dispersion interpolation
!! and corresponds to the reality of the interatomic force constants
!!
!! 2021: CL : Lifting of degeneracies using random perturbation matrix
!! with prime number matrix control tag : lphase
!! Rotate the the largest element in eigenvector to real axis. (lrot)
!
USE kinds, ONLY : DP
USE cell_base, ONLY : at, bg
USE ions_base, ONLY : amass, tau, nat, ityp
USE elph2, ONLY : rdw, epsi, zstar, qrpl
USE epwcom, ONLY : lpolar, lphase, use_ws, nqc1, nqc2, nqc3
USE epwcom, ONLY : lpolar, lphase, lrot, use_ws, nqc1, nqc2, nqc3
USE epwcom, ONLY : debug_plrn
USE constants_epw, ONLY : twopi, ci, czero, zero, one, eps12
USE rigid, ONLY : cdiagh2
USE low_lvl, ONLY : utility_zdotu
USE rigid_epw, ONLY : rgd_blk
!!!!!!
USE low_lvl, ONLY : degen_sort
!!!!!!
!
IMPLICIT NONE
!
@ -323,6 +442,10 @@
!! interpolated dynamical matrix eigenvalues for this kpoint
COMPLEX(KIND = DP), INTENT(out) :: cuf(nmodes, nmodes)
!! Rotation matrix, fine mesh
!!!!!!
LOGICAL, INTENT(IN), OPTIONAL :: is_mirror
!! .true. if q-point is a the mirror point of some original point
!!!!!!
!
! Local variables
INTEGER :: imode
@ -351,6 +474,55 @@
! Dynamical matrix in Bloch basis, fine mesh
COMPLEX(KIND = DP) :: cfac
!! Complex prefactor for Fourier transform.
!!!!!!
INTEGER, ALLOCATABLE :: degen_group(:, :)
!! Index of degenerate subspace
INTEGER :: ierr
!! Error status
INTEGER :: list_dup(nmodes)
!! List of degenerate eigenvalues
INTEGER :: ndeg
!! Number of degeneracies
LOGICAL :: duplicates
!! Returns if the bands contains degeneracices for that k-point.
INTEGER :: ig
!! Counter on real-space index
!JLB
INTEGER :: info
!! "0" successful exit, "<0" i-th argument had an illegal value, ">0" i eigenvectors failed to converge.
REAL(KIND = DP) :: rwork_tri(3*nmodes)
!! Real work array for TRI q case
REAL(KIND = DP) :: rchf(nmodes, nmodes), norm_vec(nmodes)
!! Real Dynamical matrix in Bloch basis for TRI q points
INTEGER, ALLOCATABLE :: iwork(:)
!! IWORK(1) returns the optimal LIWORK.
REAL(KIND = DP) :: rand1
!! Random number
REAL(KIND = DP) :: rand2
!! Random number
REAL(KIND = DP), ALLOCATABLE :: wp(:)
!! Perturbed eigenvalues on the degenerate subspace
REAL(KIND = DP), ALLOCATABLE :: rwork(:)
!! RWORK(1) returns the optimal LRWORK.
COMPLEX(KIND = DP) :: P(nmodes, nmodes)
!! Perturbation matrix made of small complex random number on the full space
COMPLEX(KIND = DP) :: fac_max(1)
!!
COMPLEX(KIND = DP), ALLOCATABLE :: cwork(:)
!! Complex work variable
REAL(KIND = DP), ALLOCATABLE :: P_prime(:, :)
!! Perturbation matrix on the subspace
COMPLEX(KIND = DP), ALLOCATABLE :: Uk(:, :)
!! Rotation matrix on the full space
INTEGER :: starting
!! Starting position
INTEGER :: ending
!! Ending position
INTEGER :: length
!! Size of the degenerate subspace
INTEGER :: imode_max(1)
!! Size of the degenerate subspace
!!!!!!
!
CALL start_clock ('DynW2B')
!----------------------------------------------------------
@ -433,6 +605,78 @@
! 0, 0, -one, neig, w, cz, nmodes, cwork, &
! rwork, iwork, ifail, info)
CALL cdiagh2(nmodes, chf, nmodes, w, cz)
!!!!!!
! Find the degenerate eigenvalues w
CALL degen_sort(w, SIZE(w), duplicates, list_dup)
!
ndeg = MAXVAL(list_dup)
ALLOCATE(degen_group(2, ndeg), STAT = ierr)
IF (ierr /= 0) CALL errore('hamwan2bloch', 'Error allocating degen_group', 1)
degen_group(:, :) = 0
!
! degen_group contains the starting and ending position of each group
! degen_group(1,1) = starting position of group 1
! degen_group(2,1) = ending position of group 1
! degen_group(1,2) = starting position of group 2 ...
DO ig = 1, ndeg
degen_group(2, ig) = 0
DO jmode = 1, nmodes
IF (list_dup(jmode) == ig) THEN
IF (jmode == 1) THEN
degen_group(1, ig) = jmode
ELSE
IF (list_dup(jmode) - list_dup(jmode - 1) /= 0) degen_group(1, ig) = jmode
ENDIF
degen_group(2, ig) = degen_group(2, ig) + 1
ENDIF
ENDDO
degen_group(2, ig) = degen_group(1, ig) + degen_group(2, ig) -1
ENDDO
!
! Generate a pertubation matrix of size (nbnd x nbnd) made of random number
!
DO ig = 1, ndeg
starting = degen_group(1, ig)
ending = degen_group(2, ig)
! Size of the degenerate subspace
length = ending - starting + 1
!
ALLOCATE(rwork(length**2 + 2 * length), STAT = ierr)
IF (ierr /= 0) CALL errore('hamwan2bloch', 'Error allocating rwork', 1)
ALLOCATE(iwork(3 + 5 * length), STAT = ierr)
IF (ierr /= 0) CALL errore('hamwan2bloch', 'Error allocating iwork', 1)
ALLOCATE(cwork(length**2 + 2 * length), STAT = ierr)
IF (ierr /= 0) CALL errore('hamwan2bloch', 'Error allocating cwork', 1)
ALLOCATE(Uk(nmodes, length), STAT = ierr)
IF (ierr /= 0) CALL errore('hamwan2bloch', 'Error allocating Uk', 1)
ALLOCATE(P_prime(length, length), STAT = ierr)
IF (ierr /= 0) CALL errore('hamwan2bloch', 'Error allocating P_prime', 1)
ALLOCATE(wp(length), STAT = ierr)
IF (ierr /= 0) CALL errore('hamwan2bloch', 'Error allocating wp', 1)
!
Uk(:, 1:length) = cz(:, starting:ending)
CALL prime_number_matrix(P_prime, length)
! Diagonalization of P_prime
CALL DSYEV('V', 'L', length, P_prime, length, wp, rwork, &
length**2 + 2 * length, info)
! On exiting P_prime is the eigenvector of the P_prime matrix and wp the eigenvector.
!
IF(lphase) cz(:, starting:ending) = MATMUL(Uk, P_prime)
!
DEALLOCATE(rwork, STAT = ierr)
IF (ierr /= 0) CALL errore('hamwan2bloch', 'Error deallocating rwork', 1)
DEALLOCATE(iwork, STAT = ierr)
IF (ierr /= 0) CALL errore('hamwan2bloch', 'Error deallocating iwork', 1)
DEALLOCATE(cwork, STAT = ierr)
IF (ierr /= 0) CALL errore('hamwan2bloch', 'Error deallocating cwork', 1)
DEALLOCATE(Uk, STAT = ierr)
IF (ierr /= 0) CALL errore('hamwan2bloch', 'Error deallocating Uk', 1)
DEALLOCATE(P_prime, STAT = ierr)
IF (ierr /= 0) CALL errore('hamwan2bloch', 'Error deallocating P_prime', 1)
DEALLOCATE(wp, STAT = ierr)
IF (ierr /= 0) CALL errore('hamwan2bloch', 'Error deallocating wp', 1)
ENDDO ! ig
!!!!!!
!
! clean noise
DO jmode = 1,nmodes
@ -440,6 +684,19 @@
IF (ABS(cz(imode, jmode)) < eps12) cz(imode, jmode) = czero
ENDDO
ENDDO
!!!!!!
! Find the largest element and set it to pure real
IF(lrot) THEN
DO imode = 1, nmodes
DO jmode = 1, nmodes
norm_vec(jmode) = ABS(cz(jmode, imode))
END DO
imode_max(:) = MAXLOC(norm_vec(1:nmodes))
fac_max(1) = cz(imode_max(1), imode)/norm_vec(imode_max(1))
cz(1:nmodes, imode) = cz(1:nmodes, imode) * CONJG(fac_max(1))
END DO
END IF
!!!!!!
!
! DS - Impose phase
IF (lphase) THEN
@ -458,6 +715,11 @@
!
cuf = cz
eig = w
!!!!!!
IF(PRESENT(is_mirror)) THEN
IF(is_mirror) cuf = CONJG(cz)
END IF
!!!!!!
!
CALL stop_clock('DynW2B')
!
@ -466,18 +728,27 @@
!--------------------------------------------------------------------------
!
!--------------------------------------------------------------------------
SUBROUTINE dynifc2blochf(nmodes, rws, nrws, xxq, cuf, eig)
!!!!!
!SUBROUTINE dynifc2blochf(nmodes, rws, nrws, xxq, cuf, eig)
SUBROUTINE dynifc2blochf(nmodes, rws, nrws, xxq, cuf, eig, is_mirror)
!!!!!
!--------------------------------------------------------------------------
!!
!! From the IFCs in the format of q2r, find the corresponding
!! dynamical matrix for a given q point (as in matdyn.x) on the fine grid
!!
!! 2021: CL : Lifting of degeneracies using random perturbation matrix
!! with prime number matrix control tag : lphase
!! Rotate the the largest element in eigenvector to real axis. (lrot)
!
USE kinds, ONLY : DP
USE cell_base, ONLY : at, bg
USE ions_base, ONLY : amass, tau, nat, ityp
USE elph2, ONLY : ifc, epsi, zstar, wscache, qrpl
USE epwcom, ONLY : lpolar, nqc1, nqc2, nqc3, lphase
!!!!!!
USE epwcom, ONLY : debug_plrn, lrot
USE low_lvl, ONLY : degen_sort
!!!!!!
USE io_global, ONLY : stdout
USE rigid_epw, ONLY : rgd_blk
USE low_lvl, ONLY : utility_zdotu
@ -497,6 +768,10 @@
!! interpolated phonon eigenvalues for this qpoint
COMPLEX(KIND = DP), INTENT(out) :: cuf(nmodes, nmodes)
!! Rotation matrix, fine mesh
!!!!!!
LOGICAL, INTENT(IN), OPTIONAL :: is_mirror
!! .true. if q-point is a time-reversal point
!!!!!!
!
! Local variables
LOGICAL, SAVE :: first = .TRUE.
@ -557,6 +832,14 @@
!! Dynamical matrix in Bloch basis, fine mesh
COMPLEX(KIND = DP) :: dyn(3, 3, nat, nat)
!! Dynamical matrix
!!!!!
REAL(KIND = DP) :: norm_vec(nmodes)
!! Vector for one eigenmode
INTEGER :: imode_max(1)
!! index of the max element
COMPLEX(KIND = DP) :: fac_max(1)
!! value of the max element
!!!!!
!
CALL start_clock('DynW2B')
!
@ -686,6 +969,19 @@
CALL zhpevx('V', 'A', 'U', nmodes, champ , zero, zero, &
0, 0, -one, neig, w, cz, nmodes, cwork, rwork, iwork, ifail, info)
!
!!!!!
! Find the largest element and set it to pure real
IF(lrot) THEN
DO imode = 1, nmodes
DO jmode = 1, nmodes
norm_vec(jmode) = ABS(cz(jmode, imode))
END DO
imode_max(:) = MAXLOC(norm_vec(1:nmodes))
fac_max(1) = cz(imode_max(1), imode)/norm_vec(imode_max(1))
cz(1:nmodes, imode) = cz(1:nmodes, imode) * CONJG(fac_max(1))
END DO
END IF
!!!!!
! clean noise
DO jmode = 1,nmodes
DO imode = 1,nmodes

136
EPW/src/wfc_elec.f90
View File

@ -1,136 +0,0 @@
!
! Copyright (C) 2010-2016 Samuel Ponce', Roxana Margine, Carla Verdi, Feliciano Giustino
! Copyright (C) 2007-2009 Jesse Noffsinger, Brad Malone, Feliciano Giustino
!
! This file is distributed under the terms of the GNU General Public
! License. See the file `LICENSE' in the root directory of the
! present distribution, or http://www.gnu.org/copyleft.gpl.txt .
!
!-----------------------------------------------------------------------
MODULE polaron
USE kinds, ONLY : dp
IMPLICIT NONE
! Data block, try to keep it to minimal
COMPLEX(KIND = DP), ALLOCATABLE :: epfall(:, :, :, :, :)
!! el-ph element for all local k and all q
!! epfall need to be filled in ephwann_shuffle
COMPLEX(KIND = DP), ALLOCATABLE :: ufall(:, :, :)
!! el-ph element for all local k and all q
!! epfall need to be filled in ephwann_shuffle
COMPLEX(KIND = DP), ALLOCATABLE :: Hamil(:, :)
!! Hamil need to be passed to h_psi because the parameter space is fixed
!! to meet the requirement of Davidson diagonalization. Ugly but workable.
COMPLEX(KIND = DP), ALLOCATABLE :: eigVec(:, :)
!! polaron eigenvector
REAL(KIND = DP), ALLOCATABLE :: etf_all(:, :)
!! Gather all the eigenvalues
INTEGER, ALLOCATABLE :: ikq_all(:, :), kpg_map(:)
!
PUBLIC :: wfc_elec, interp_plrn_wf, interp_plrn_bq, plot_plrn_wf
CONTAINS
!
!-----------------------------------------------------------------------
SUBROUTINE interp_plrn_wf(nrr_k, ndegen_k, irvec_r, dims)
USE io_global, ONLY : stdout, ionode
IMPLICIT NONE
INTEGER, INTENT (IN) :: nrr_k, dims, ndegen_k(:,:,:) ! ! Added for polaron calculations by Chao Lian.
REAL(DP), INTENT (IN) :: irvec_r(3, nrr_k)
COMPLEX(DP), ALLOCATABLE :: eigvec_wan(:, :)
INTEGER :: nkf1_p, nkf2_p, nkf3_p, nktotf_p, nbndsub_p
END SUBROUTINE
!
!-----------------------------------------------------------------------
SUBROUTINE interp_plrn_bq(nrr_q, ndegen_q, irvec_q)
USE epwcom, ONLY : nkf1, nkf2, nkf3, nbndsub
USE elph2, only : xqf, wf, nqtotf
USE modes, ONLY : nmodes
USE constants_epw, only : eps8, czero, one, two, twopi, ci
USE ions_base, ONLY : nat, amass, ityp, tau
USE wan2bloch, only : dynwan2bloch
IMPLICIT NONE
INTEGER, INTENT (IN) :: nrr_q, ndegen_q(:,:,:) ! ! Added for polaron calculations by Chao Lian.
INTEGER, INTENT (IN) :: irvec_q(3, nrr_q)
INTEGER :: dtau_file
INTEGER :: nkf1_p, nkf2_p, nkf3_p, nktotf_p, nat_p
INTEGER :: iq, inu, ierr, imu, na, iatm, idir
INTEGER :: icount, ix, iy, iz, bmat_file
COMPLEX(DP) :: ctemp, shift(3)
COMPLEX(DP), ALLOCATABLE :: uf(:, :), Bmat(:,:)
COMPLEX(DP), ALLOCATABLE :: dtau(:, :)
REAL(DP), ALLOCATABLE :: w2(:)
REAL(KIND=dp) :: xxq(3)
COMPLEX(KIND=dp) :: expTable(3)
END SUBROUTINE
SUBROUTINE wfc_elec (nrr_k, ndegen_k, irvec_r, dims)
!
! Self consistency calculation of polaron wavefunction.
! Rewritten by Chao Lian based on the implementation by Danny Sio.
!
USE modes, ONLY : nmodes
USE constants_epw, ONLY : ryd2mev, one, ryd2ev, two, zero
USE constants_epw, ONLY : czero, cone, pi, ci, twopi, eps6, eps8, eps5
USE epwcom, ONLY : num_cbands, polaron_type, sigma_plrn, full_diagon_plrn
USE epwcom, ONLY : r01, r02, r03, nPlrn, conv_thr_polaron, cb_shift
USE epwcom, ONLY : mixing_Plrn, init_plrn_wf, niterPlrn
USE epwcom, ONLY : nkf1, nkf2, nkf3, nbndsub
USE io_global, ONLY : stdout, ionode, meta_ionode_id
USE elph2, ONLY : etf, ibndmin, ibndmax, nbndfst
USE elph2, ONLY : nkqf, nkf, nqf, nqtotf, nktotf
USE elph2, ONLY : xkf, xqf, wf, xkq, chw
USE mp_global, ONLY : inter_pool_comm
USE mp_world, ONLY : world_comm
USE cell_base, ONLY : bg
USE mp, ONLY : mp_sum, mp_bcast
USE poolgathering, ONLY : poolgather2
USE test_tools, ONLY : para_write
USE wan2bloch, ONLY : hamwan2bloch
USE ions_base, ONLY : nat
IMPLICIT NONE
! local variables
LOGICAL :: debug
INTEGER :: inu, iq, ik, ikk, jk, ibnd, jbnd, ikq, ik_global, iplrn, ierr
INTEGER :: iter, icount, ix, iy, iz, start_mode, ik_bm, idos, iatm
INTEGER, INTENT (IN) :: nrr_k, dims, ndegen_k(:,:,:) ! ! Added for polaron calculations by Chao Lian.
REAL(DP), INTENT (IN) :: irvec_r(3, nrr_k)
COMPLEX(DP), ALLOCATABLE :: Bmat(:,:), Bmat_save(:,:)
COMPLEX(DP), ALLOCATABLE :: eigvec_wan(:, :), dtau(:, :)
REAL(DP), ALLOCATABLE :: rmat_tmp(:, :)
COMPLEX(KIND=dp) :: cufkk ( nbndsub, nbndsub ), cfac(nrr_k, dims, dims)
!! Rotation matrix, fine mesh, points k
REAL(dp):: estmteRt(nPlrn), eigVal(nPlrn), esterr
REAL(KIND=dp) :: qcart(3), r0(3), xxk(3), xxq(3), prefac, norm
REAL(KIND=dp) :: ef
INTEGER :: band_pos, iqpg, ikpg, ikGamma, iqGamma
INTEGER :: nkf1_p, nkf2_p, nkf3_p, nbndsub_p, nPlrn_p, nktotf_p
REAL(DP) :: eb
REAL(DP) :: xkf_all(3, nktotf), xkf_all_tmp(3, nktotf*2)
REAL(DP) :: EPlrnTot, EPlrnElec, EPlrnPhon
REAL(DP) :: disK, disK_t, shift(3)
COMPLEX(DP) :: ctemp
REAL(DP) :: rtemp
INTEGER :: itemp, jtemp
INTEGER :: dos_file, wan_func_file, bloch_func_file, bmat_file, dtau_file
!LOGICAL :: SCF_run
END SUBROUTINE
SUBROUTINE plot_plrn_wf()
END SUBROUTINE
END MODULE