Commit 9ff6acc6 authored by Samuel Poncé's avatar Samuel Poncé

Add possible imposition of gauge through

new variable lphase. Courtesy of Denny Sio.
parent 27cae165
......@@ -47,7 +47,7 @@
scattering, scattering_serta, scattering_0rta, &
int_mob, scissor, carrier, ncarrier, iterative_bte, &
restart, restart_freq, prtgkk, nel, meff, epsiHEG, &
scatread, restart, restart_freq, restart_filq
scatread, restart, restart_freq, restart_filq, lphase
USE elph2, ONLY : elph
USE mp, ONLY : mp_bcast
USE mp_world, ONLY : world_comm
......@@ -130,6 +130,7 @@
CALL mp_bcast (carrier , meta_ionode_id, world_comm)
CALL mp_bcast (restart , meta_ionode_id, world_comm)
CALL mp_bcast (prtgkk , meta_ionode_id, world_comm)
CALL mp_bcast (lphase , meta_ionode_id, world_comm)
!
! integers
!
......
......@@ -57,7 +57,7 @@
title, int_mob, scissor, iterative_bte, scattering, &
ncarrier, carrier, scattering_serta, restart, restart_freq, &
scattering_0rta, longrange, shortrange, scatread, &
restart_filq, prtgkk, nel, meff, epsiHEG
restart_filq, prtgkk, nel, meff, epsiHEG, lphase
USE elph2, ONLY : elph
USE start_k, ONLY : nk1, nk2, nk3
USE constants_epw, ONLY : ryd2mev, ryd2ev, ev2cmm1, kelvin2eV
......@@ -121,7 +121,8 @@
specfun_el, specfun_ph, wmin_specfun, wmax_specfun, nw_specfun, &
delta_approx, scattering, int_mob, scissor, ncarrier, carrier, &
iterative_bte, scattering_serta, scattering_0rta, longrange, shortrange,&
scatread, restart, restart_freq, restart_filq, prtgkk, nel, meff, epsiHEG
scatread, restart, restart_freq, restart_filq, prtgkk, nel, meff, &
epsiHEG, lphase
! tphases, fildvscf0
!
......@@ -278,6 +279,7 @@
! nel : Fractional number of electrons in the unit cell
! meff : Density of state effective mass (in unit of the electron mass)
! epsiHEG : Dielectric constant at zero doping
! lphase : If .true., fix the gauge on the phonon eigenvectors and electronic eigenvectors - DS
!
nk1tmp = 0
nk2tmp = 0
......@@ -461,7 +463,8 @@
prtgkk = .false.
nel = 0.0d0
meff = 1.d0
epsiHEG = 1.d0
epsiHEG = 1.d0
lphase = .false.
!
! reading the namelist inputepw
!
......
......@@ -234,6 +234,8 @@
!! if .true. the system is 2 dimensional (vaccum is in z-direction)
LOGICAL :: prtgkk
!! if .true. print the |g| vertex in [meV].
LOGICAL :: lphase
!! if .true. then fix the gauge when diagonalizing the interpolated dynamical matrix and electronic Hamiltonian.
!
! Superconductivity
LOGICAL :: ephwrite
......
......@@ -38,7 +38,8 @@
!--------------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE constants_epw, ONLY : czero, cone
USE constants_epw, ONLY : czero, cone, eps12
USE epwcom, ONLY : lphase
!
implicit none
!
......@@ -72,6 +73,8 @@
!
INTEGER :: ibnd, jbnd
COMPLEX(kind=DP) :: chf(nbnd, nbnd)
COMPLEX(KIND=DP) :: zdotu
!! Dot product between the two phonon eigenvectors.
!
CALL start_clock('HamW2B')
!----------------------------------------------------------
......@@ -117,7 +120,26 @@
CALL zhpevx ('V', 'A', 'U', nbnd, champ , 0.0, 0.0, &
0, 0,-1.0, neig, w, cz, nbnd, cwork, &
rwork, iwork, ifail, info)
!
! clean noise
DO jbnd = 1, nbnd
DO ibnd = 1, nbnd
IF ( ABS( cz(ibnd,jbnd) ) < eps12 ) cz(ibnd,jbnd) = cmplx(0.d0,0.d0, kind=DP)
ENDDO
ENDDO
!
! DS - Impose phase
IF (lphase) THEN
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) )
EXIT INNER
ENDIF
END DO INNER
ENDDO
ENDIF
!
! rotation matrix and Ham eigenvalues
! [in Ry, mind when comparing with wannier code]
!
......@@ -151,8 +173,8 @@
USE phcom, ONLY : nq1, nq2, nq3
USE ions_base, ONLY : amass, tau, nat, ityp
USE elph2, ONLY : rdw, epsi, zstar
USE epwcom, ONLY : lpolar
USE constants_epw, ONLY : twopi, ci, czero
USE epwcom, ONLY : lpolar, lphase
USE constants_epw, ONLY : twopi, ci, czero, eps12
!
implicit none
!
......@@ -175,16 +197,19 @@
! work variables
! variables for lapack ZHPEVX
integer :: neig, info, ifail( nmodes ), iwork( 5*nmodes )
integer :: imode, jmode, ir, na, nb
real(kind=DP) :: w( nmodes ), rwork( 7*nmodes )
complex(kind=DP) :: champ( nmodes*(nmodes+1)/2 ), &
cwork( 2*nmodes ), cz( nmodes, nmodes)
!
real(kind=DP) :: xq (3)
complex(kind=DP) :: chf(nmodes, nmodes)
! Hamiltonian in Bloch basis, fine mesh
integer :: imode, jmode, ir, na, nb
real(kind=DP) :: rdotk, massfac
complex(kind=DP) :: cfac
COMPLEX(kind=DP) :: chf(nmodes, nmodes)
! Hamiltonian in Bloch basis, fine mesh
COMPLEX(kind=DP) :: champ( nmodes*(nmodes+1)/2 ), &
cwork( 2*nmodes ), cz( nmodes, nmodes)
COMPLEX(kind=DP) :: cfac
!! Complex prefactor for Fourier transform.
COMPLEX(KIND=DP) :: zdotu
!! Dot product between the two phonon eigenvectors.
!
CALL start_clock ( 'DynW2B' )
!----------------------------------------------------------
......@@ -242,15 +267,35 @@
!
DO jmode = 1, nmodes
DO imode = 1, jmode
champ (imode + (jmode - 1) * jmode/2 ) = &
( chf ( imode, jmode) + conjg ( chf ( jmode, imode) ) ) / 2.d0
champ (imode + (jmode - 1) * jmode/2 ) = &
( chf ( imode, jmode) + conjg ( chf ( jmode, imode) ) ) / 2.d0
ENDDO
ENDDO
!
CALL zhpevx ('V', 'A', 'U', nmodes, champ , 0.0, 0.0, &
0, 0,-1.0, neig, w, cz, nmodes, cwork, &
rwork, iwork, ifail, info)
!
!
! clean noise
DO jmode=1,nmodes
DO imode=1,nmodes
IF ( ABS( cz(imode,jmode) ) < eps12 ) cz(imode,jmode) = cmplx(0.d0,0.d0, kind=DP)
ENDDO
ENDDO
!
! DS - Impose phase
IF (lphase) THEN
DO jmode = 1,nmodes
INNER : DO imode = 1,nmodes
IF ( ABS(cz(imode, jmode)) > eps12 ) THEN
cz(:, jmode) = cz(:, jmode) * conjg( cz(imode,jmode) )
cz(:, jmode) = cz(:, jmode)/sqrt( zdotu(nmodes,conjg(cz(:, jmode)),1,cz(:, jmode),1) )
EXIT INNER
ENDIF
END DO INNER
ENDDO
ENDIF
!
! rotation matrix and Ham eigenvalues
! [in Ry, mind when comparing with wannier code]
!
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment