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

Make [H,x] commutator work for arbitrary direction 

Previously, it assumed crystal coordinate at(:,1) or at(:,2) or at(:,3).
Now, it accepts vpol(1:3) as input.
This commit is contained in:
Jae-Mo Lihm 2021-04-05 11:24:11 +09:00
parent cc1182dcc8
commit 9f3f6384b9
4 changed files with 60 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
PHonon/PH/dvpsi_e.f90
View File

@ -17,7 +17,7 @@ subroutine dvpsi_e (ik, ipol)
! otherwise dvpsi is COMPUTED and WRITTEN on file (vkb and evc must be set)
!
USE kinds, ONLY : DP
USE cell_base, ONLY : tpiba2
USE cell_base, ONLY : tpiba2, at
USE io_global, ONLY : stdout
USE klist, ONLY : xk, ngk, igk_k
USE gvect, ONLY : g
@ -79,7 +79,7 @@ subroutine dvpsi_e (ik, ipol)
! calculate the commutator [H,x_ipol] psi > and store it in dpsi
call commutator_Hx_psi (ikk, nbnd_occ(ikk), becp1(ik), becp2, ipol, dpsi )
call commutator_Hx_psi (ikk, nbnd_occ(ikk), at(:, ipol), becp1(ik), becp2, dpsi )
!
! orthogonalize dpsi to the valence subspace: ps = <evc|dpsi>
! Apply -P^+_c

57
PW/src/commutator_Hx_psi.f90
View File

@ -6,16 +6,19 @@
! or http://www.gnu.org/copyleft/gpl.txt .
!
!----------------------------------------------------------------------
subroutine commutator_Hx_psi (ik, nbnd_occ, becp1, becp2, ipol, dpsi)
subroutine commutator_Hx_psi (ik, nbnd_occ, vpol, becp1, becp2, dpsi)
!----------------------------------------------------------------------
!
! On output: dpsi contains [H,x_ipol] | psi_ik > in crystal axis
! (projected on at(*,ipol) )
! On output: dpsi contains [H, r \dot vpol] | psi_ik >
!
! vpol is the polarization vector in Cartesian coordinates.
! For crystal coordinate, use vpol = at(:, ipol).
! For Cartesian coordinate, use vpol = (1.0, 0.0, 0.0) or other permutations.
!
! vkb and evc must be properly set for the appropriate k-point
! in addition becp1 must be set equal to becp1 = <vkb|evc>
! as it is done in PH/phq_init.f90 for the k-point ik
! NB: here the last index of becp1 is missing, hence it refers
! NB: here the last index of becp1 is missing, hence it refers
! to a single k-point
!
! CALL calbec (npw, vkb, evc, becp1(:,:) )
@ -41,43 +44,44 @@ subroutine commutator_Hx_psi (ik, nbnd_occ, becp1, becp2, ipol, dpsi)
TYPE(bec_type), INTENT(IN) :: becp1 ! dimensions ( nkb, nbnd )
TYPE(bec_type), INTENT(INOUT) :: becp2 ! dimensions ( nkb, nbnd )
!
INTEGER, INTENT(IN) :: ik, nbnd_occ, ipol
INTEGER, INTENT(IN) :: ik, nbnd_occ
REAL(DP), INTENT(IN) :: vpol(3)
!! polarization vector in Cartesian coordinates
!
! Local variables
!
integer :: npw, ig, na, ibnd, jbnd, ikb, jkb, nt, lter, ih, jh, ijkb0, &
nrec, is, js, ijs
! counters
real(DP), allocatable :: gk (:,:), g2k(:)
real(DP) :: gk_ig(3)
real(DP), allocatable :: g2k(:), gk_vpol(:)
! the derivative of |k+G|
complex(DP), allocatable :: ps2(:,:,:), dvkb (:,:), dvkb1 (:,:), &
work (:,:), psc(:,:,:,:), deff_nc(:,:,:,:)
REAL(DP), allocatable :: deff(:,:,:)
!
CALL start_clock ('commutator_Hx_psi')
dpsi=(0.d0, 0.d0)
!
npw = ngk(ik)
allocate (gk(3, npw), g2k(npw) )
allocate (gk_vpol(npw), g2k(npw) )
do ig = 1, npw
gk (1:3, ig) = (xk (1:3, ik) + g (1:3, igk_k(ig,ik) ) ) * tpiba
g2k (ig) = SUM(gk (1:3, ig) **2 )
gk_ig(1:3) = (xk (1:3, ik) + g (1:3, igk_k(ig,ik) ) ) * tpiba
g2k (ig) = SUM(gk_ig**2)
!
! Take the component along the vpol vector
gk_vpol(ig) = SUM(vpol * gk_ig(:))
enddo
!
! this is the kinetic contribution to [H,x]: -2i (k+G)_ipol * psi
!
do ibnd = 1, nbnd_occ
do ig = 1, npw
dpsi(ig,ibnd) = SUM(at(1:3,ipol)*gk(1:3,ig))*(0.d0,-2.d0)*evc (ig,ibnd)
dpsi(ig,ibnd) = gk_vpol(ig)*(0.d0,-2.d0)*evc (ig,ibnd)
enddo
IF (noncolin) THEN
do ig = 1, npw
dpsi (ig+npwx, ibnd) = (at(1, ipol) * gk(1, ig) + &
at(2, ipol) * gk(2, ig) + &
at(3, ipol) * gk(3, ig) ) &
*(0.d0,-2.d0)*evc (ig+npwx, ibnd)
dpsi (ig+npwx, ibnd) = gk_vpol(ig)*(0.d0,-2.d0)*evc (ig+npwx, ibnd)
end do
END IF
enddo
@ -101,18 +105,14 @@ subroutine commutator_Hx_psi (ik, nbnd_occ, becp1, becp2, ipol, dpsi)
allocate (dvkb (npwx, nkb), dvkb1(npwx, nkb))
dvkb (:,:) = (0.d0, 0.d0)
dvkb1(:,:) = (0.d0, 0.d0)
call gen_us_dj (ik, dvkb)
call gen_us_dy (ik, at (1, ipol), dvkb1)
call gen_us_dy (ik, vpol, dvkb1)
do ig = 1, npw
if (g2k (ig) < 1.0d-10) then
gk (1, ig) = 0.d0
gk (2, ig) = 0.d0
gk (3, ig) = 0.d0
gk_vpol(ig) = 0.d0
else
gk (1, ig) = gk (1, ig) / sqrt (g2k (ig) )
gk (2, ig) = gk (2, ig) / sqrt (g2k (ig) )
gk (3, ig) = gk (3, ig) / sqrt (g2k (ig) )
gk_vpol(ig) = gk_vpol(ig) / sqrt (g2k (ig) )
endif
enddo
@ -124,16 +124,13 @@ subroutine commutator_Hx_psi (ik, nbnd_occ, becp1, becp2, ipol, dpsi)
do ikb = 1, nh (nt)
jkb = jkb + 1
do ig = 1, npw
work (ig,jkb) = dvkb1 (ig, jkb) + dvkb (ig, jkb) * &
(at (1, ipol) * gk (1, ig) + &
at (2, ipol) * gk (2, ig) + &
at (3, ipol) * gk (3, ig) )
work (ig,jkb) = dvkb1 (ig, jkb) + dvkb (ig, jkb) * gk_vpol(ig)
enddo
enddo
endif
enddo
enddo
deallocate (g2k, gk)
deallocate (g2k, gk_vpol)
! In the case of gamma point systems becp2 is real
! so we have to include a factor of i before calling
@ -230,7 +227,7 @@ subroutine commutator_Hx_psi (ik, nbnd_occ, becp1, becp2, ipol, dpsi)
! Compute the commutator between the non-local Hubbard potential
! and the position operator
!
IF (lda_plus_u) CALL commutator_Vhubx_psi(ik, ipol, nbnd_occ, dpsi)
IF (lda_plus_u) CALL commutator_Vhubx_psi(ik, nbnd_occ, vpol, dpsi)
!
111 continue
!

56
PW/src/commutator_Vhubx_psi.f90
View File

@ -7,14 +7,18 @@
!
!
!-----------------------------------------------------------------------
SUBROUTINE commutator_Vhubx_psi(ik, ipol, nbnd_calc, dpsi)
SUBROUTINE commutator_Vhubx_psi(ik, nbnd_calc, vpol, dpsi)
!-----------------------------------------------------------------------
!
! This routine computes the commutator between the non-local
! Hubbard potential and the position operator, applied to psi
! of the current k point, i.e. [V_hub,r]|psi_nk> .
! of the current k point, i.e. [V_hub, r \dot vpol]|psi_nk>.
! The result is added to dpsi.
!
! vpol is the polarization vector in Cartesian coordinates.
! For crystal coordinate, use vpol = at(:, ipol).
! For Cartesian coordinate, use vpol = (1.0, 0.0, 0.0) or other permutations.
!
! Some insights about the formulas here can be found e.g.
! in I. Timrov's PhD thesis, Sec. 6.1.3,
! https://pastel.archives-ouvertes.fr/pastel-00823758
@ -33,7 +37,7 @@ SUBROUTINE commutator_Vhubx_psi(ik, ipol, nbnd_calc, dpsi)
USE uspp_param, ONLY : nh, upf
USE lsda_mod, ONLY : lsda, current_spin, isk, nspin
USE klist, ONLY : xk, ngk, igk_k
USE cell_base, ONLY : tpiba, at
USE cell_base, ONLY : tpiba
USE gvect, ONLY : g
USE scf, ONLY : rho
USE mp, ONLY : mp_sum
@ -46,18 +50,20 @@ SUBROUTINE commutator_Vhubx_psi(ik, ipol, nbnd_calc, dpsi)
!
INTEGER, INTENT(IN) :: ik
!! k point index
INTEGER, INTENT(IN) :: ipol
!! polarization in crystal coordinates
INTEGER, INTENT(IN) :: nbnd_calc
!! Number of bands to calculate [V_hub, x_ipol]|psi_ik>
REAL(DP), INTENT(IN) :: vpol(3)
!! polarization vector in Cartesian coordinates
COMPLEX(DP), INTENT(OUT) :: dpsi(npwx*npol, nbnd)
!! Output wavefunction where [V_hub, x_ipol]|psi_ik> is added
!
REAL(DP), PARAMETER :: eps = 1.0d-8
INTEGER :: na, n ,l, nt, nah, ikb , m, m1, m2, ibnd, ib, ig, jkb, i, &
ihubst, ihubst1, ihubst2, icart, op_spin, npw, offpm, offpmU
REAL(DP) :: nsaux
REAL(DP), ALLOCATABLE :: xyz(:,:), gk(:,:), g2k(:)
REAL(DP) :: nsaux, g2k, gk_ig(3)
REAL(DP), ALLOCATABLE :: xyz(:,:)
REAL(DP), ALLOCATABLE :: gk_vpol(:)
!! ipol component of the G/|G| vector, in crystal coordinates
COMPLEX(DP), ALLOCATABLE :: dkwfcbessel(:,:), dkwfcylmr(:,:), dkwfcatomk(:,:), &
dpqq26(:,:), dpqq38(:,:), dpqq47(:,:), dkvkbbessel(:,:), &
dkvkbylmr(:,:), dkvkb(:,:), aux_1234(:), termi(:,:), trm(:,:), &
@ -87,8 +93,7 @@ SUBROUTINE commutator_Vhubx_psi(ik, ipol, nbnd_calc, dpsi)
ALLOCATE (termi(npwx,nbnd))
ALLOCATE (trm(npwx,nbnd))
ALLOCATE (xyz(3,3))
ALLOCATE (gk(3,npw))
ALLOCATE (g2k(npw))
ALLOCATE (gk_vpol(npw))
!
dpqq26 = (0.d0, 0.d0)
dpqq38 = (0.d0, 0.d0)
@ -124,7 +129,7 @@ SUBROUTINE commutator_Vhubx_psi(ik, ipol, nbnd_calc, dpsi)
! Derivative of Bessel functions and spherical harmonics wrt to crystal axis ipol
!
CALL gen_at_dj (ik, dkwfcbessel)
CALL gen_at_dy (ik, at(:,ipol), dkwfcylmr)
CALL gen_at_dy (ik, vpol, dkwfcylmr)
!
DO ig = 1, npw
!
@ -133,16 +138,22 @@ SUBROUTINE commutator_Vhubx_psi(ik, ipol, nbnd_calc, dpsi)
! the cartesian component and then to crystal component ipol
! gk_icart= d|k+G|/d(k+G)_icart
!
gk(:,ig) = (xk(:,ik) + g(:,igk_k(ig,ik))) * tpiba
gk_ig = (xk(:,ik) + g(:,igk_k(ig,ik))) * tpiba
g2k = SUM( gk_ig**2 )
!
g2k(ig) = SUM( gk(1:3,ig)**2 )
! Take the component along the vpol vector
gk_vpol(ig) = SUM(vpol(:) * gk_ig(:))
!
IF (g2k(ig) < 1.0d-10) THEN
gk(:,ig) = 0.d0
IF (g2k < 1.0d-10) THEN
gk_vpol(ig) = 0.d0
ELSE
gk(:,ig) = gk(:,ig) / SQRT(g2k(ig))
gk_vpol(ig) = gk_vpol(ig) / SQRT(g2k)
ENDIF
!
ENDDO
!
DO ig = 1, npw
!
! Derivative wrt crystal axis ipol
! d|k+G|/d(k+G)_ipol = \sum_{icart} d|k+G|/d(k+G)_icart * at (icart,ipol)
! The derivative is done for all the atomic wfc and for each ig
@ -154,10 +165,7 @@ SUBROUTINE commutator_Vhubx_psi(ik, ipol, nbnd_calc, dpsi)
offpm = oatwfc(na)
DO m1 = 1, 2*Hubbard_l(nt)+1
dkwfcatomk(ig,offpmU+m1) = dkwfcylmr(ig,offpm+m1) &
+ dkwfcbessel(ig,offpm+m1) * &
( at (1, ipol) * gk (1, ig) + &
at (2, ipol) * gk (2, ig) + &
at (3, ipol) * gk (3, ig) )
+ dkwfcbessel(ig,offpm+m1) * gk_vpol(ig)
ENDDO
ENDIF
ENDDO
@ -172,7 +180,7 @@ SUBROUTINE commutator_Vhubx_psi(ik, ipol, nbnd_calc, dpsi)
! Here the derivative of the Bessel functions and the spherical harmonics
!
CALL gen_us_dj (ik, dkvkbbessel)
CALL gen_us_dy (ik, at(:,ipol), dkvkbylmr)
CALL gen_us_dy (ik, vpol, dkvkbylmr)
!
jkb = 0
DO nt = 1, ntyp
@ -181,10 +189,7 @@ SUBROUTINE commutator_Vhubx_psi(ik, ipol, nbnd_calc, dpsi)
DO ikb = 1, nh(nt)
jkb = jkb + 1
DO ig = 1, npw
dkvkb(ig,jkb) = dkvkbylmr(ig,jkb) + dkvkbbessel(ig,jkb) * &
(at (1, ipol) * gk (1, ig) + &
at (2, ipol) * gk (2, ig) + &
at (3, ipol) * gk (3, ig) )
dkvkb(ig,jkb) = dkvkbylmr(ig,jkb) + dkvkbbessel(ig,jkb) * gk_vpol(ig)
ENDDO
ENDDO
ENDIF
@ -388,8 +393,7 @@ SUBROUTINE commutator_Vhubx_psi(ik, ipol, nbnd_calc, dpsi)
DEALLOCATE (termi)
DEALLOCATE (trm)
DEALLOCATE (xyz)
DEALLOCATE (gk)
DEALLOCATE (g2k)
DEALLOCATE (gk_vpol)
!
CALL stop_clock ('commutator_Vhubx_psi')
!

2
TDDFPT/src/lr_dvpsi_e.f90
View File

@ -84,7 +84,7 @@ SUBROUTINE lr_dvpsi_e(ik,ipol,dvpsi)
!
CALL allocate_bec_type ( nkb, nbnd, becp2 )
!
CALL commutator_Hx_psi (ik, nbnd_occ(ik), becp1, becp2, ipol, d0psi )
CALL commutator_Hx_psi (ik, nbnd_occ(ik), at(:, ipol), becp1, becp2, d0psi )
!
IF (okvan) CALL calbec ( npw, vkb, evc, becp, nbnd)
!