From c2ef20a5e7f0bb39b8aceeb8d31726c71f1d4aa2 Mon Sep 17 00:00:00 2001 From: Samuel Ponce Date: Mon, 4 Feb 2019 18:02:45 +0000 Subject: [PATCH] Update to EPW USPP to support time-reversal symm. Courtesy of R. Margine. --- EPW/src/adddvscf2.f90 | 12 +- EPW/src/allocate_epwq.f90 | 10 +- EPW/src/close_epw.f90 | 40 +++-- EPW/src/dvanqq2.f90 | 154 +++++++++++-------- EPW/src/dvqpsi_us3.f90 | 252 +++++++++++++++++--------------- EPW/src/dvqpsi_us_only3.f90 | 116 ++++++++------- EPW/src/elphel2_shuffle.f90 | 67 +++++---- EPW/src/elphon_shuffle.f90 | 49 +------ EPW/src/elphon_shuffle_wrap.f90 | 47 +++--- EPW/src/newdq2.f90 | 78 ++++++---- EPW/src/printing.f90 | 4 +- test-suite/extract-epw.x | 2 +- 12 files changed, 435 insertions(+), 396 deletions(-) diff --git a/EPW/src/adddvscf2.f90 b/EPW/src/adddvscf2.f90 index 0e54ab31c..5c9033321 100644 --- a/EPW/src/adddvscf2.f90 +++ b/EPW/src/adddvscf2.f90 @@ -67,9 +67,9 @@ INTEGER :: ijs !! Counter on combined is and js polarization ! - COMPLEX(DP) :: sumA + COMPLEX(kind=DP) :: sum_k !! auxiliary sum variable - COMPLEX(DP) :: sum_nc(npol) + COMPLEX(kind=DP) :: sum_nc(npol) !! auxiliary sum variable non-collinear case ! IF (.not.okvan) RETURN @@ -87,12 +87,12 @@ ! we multiply the integral for the becp term and the beta_n ! DO ibnd = lower_band, upper_band - do ih = 1, nh(nt) + DO ih = 1, nh(nt) ikb = ijkb0 + ih IF (noncolin) THEN sum_nc = czero ELSE - sumA = czero + sum_k = czero ENDIF DO jh = 1, nh(nt) jkb = ijkb0 + jh @@ -106,7 +106,7 @@ ENDDO ENDDO ELSE - sumA = sumA + int3(ih,jh,na,current_spin,ipert) * & + sum_k = sum_k + int3(ih,jh,na,current_spin,ipert) * & becp1(ik)%k(jkb,ibnd) ENDIF ENDDO @@ -114,7 +114,7 @@ CALL zaxpy( npwq, sum_nc(1), vkb(1,ikb), 1, dvpsi(1,ibnd), 1 ) CALL zaxpy( npwq, sum_nc(2), vkb(1,ikb), 1, dvpsi(1+npwx,ibnd), 1 ) ELSE - CALL zaxpy( npwq, sumA, vkb(1,ikb), 1, dvpsi(1,ibnd), 1 ) + CALL zaxpy( npwq, sum_k, vkb(1,ikb), 1, dvpsi(1,ibnd), 1 ) ENDIF ENDDO ENDDO diff --git a/EPW/src/allocate_epwq.f90 b/EPW/src/allocate_epwq.f90 index ada00d25c..db378f447 100644 --- a/EPW/src/allocate_epwq.f90 +++ b/EPW/src/allocate_epwq.f90 @@ -25,10 +25,10 @@ USE gvect, ONLY : ngm USE noncollin_module, ONLY : noncolin, npol, nspin_mag USE spin_orb, ONLY : lspinorb - USE phcom, ONLY : evq, dpsi, vlocq, dmuxc + USE phcom, ONLY : evq, vlocq, dmuxc USE phus, ONLY : int1, int1_nc, int2, int2_so, & - int4, int4_nc, int5, int5_so, becsum_nc, & - alphasum, alphasum_nc, alphap + int4, int4_nc, int5, int5_so, & + alphap USE lr_symm_base, ONLY : rtau USE qpoint, ONLY : eigqts USE lrus, ONLY : becp1 @@ -52,7 +52,6 @@ ! ALLOCATE space for the quantities needed in EPW ! ALLOCATE (evq(npwx*npol, nbnd)) - ALLOCATE (dpsi ( npwx*npol, nbnd)) ALLOCATE (transp_temp(nstemp)) ! ALLOCATE (vlocq(ngm, ntyp)) @@ -82,14 +81,11 @@ IF (noncolin) THEN ALLOCATE (int1_nc(nhm, nhm, 3, nat, nspin)) ALLOCATE (int4_nc(nhm, nhm, 3, 3, nat, nspin)) - ALLOCATE (becsum_nc(nhm*(nhm+1)/2, nat, npol, npol)) - ALLOCATE (alphasum_nc(nhm*(nhm+1)/2, 3, nat, npol, npol)) IF (lspinorb) THEN ALLOCATE (int2_so(nhm, nhm, 3, nat, nat, nspin)) ALLOCATE (int5_so(nhm, nhm, 3, 3, nat, nat, nspin)) ENDIF ENDIF ! noncolin - ALLOCATE (alphasum(nhm * (nhm + 1)/2, 3, nat, nspin_mag)) ENDIF ! ALLOCATE (becp1(nks)) diff --git a/EPW/src/close_epw.f90 b/EPW/src/close_epw.f90 index c3f677825..4fe899405 100644 --- a/EPW/src/close_epw.f90 +++ b/EPW/src/close_epw.f90 @@ -104,9 +104,7 @@ !! Imported the noncolinear case implemented by xlzhang !! !---------------------------------------------------------------------- - USE phcom, ONLY : alphap, alphasum, alphasum_nc, & - becsum_nc, dmuxc, dpsi,& - drc, dpsi, dyn, evq, dvpsi,& + USE phcom, ONLY : alphap, dmuxc, drc, dyn, evq, dvpsi, & int5, vlocq, int2_so, int5_so USE lrus, ONLY : becp1, int3, int3_nc USE phus, ONLY : int1, int1_nc, int2, int4, int4_nc @@ -154,7 +152,6 @@ IF(ASSOCIATED(igkq)) DEALLOCATE(igkq) ! IF(ALLOCATED(dvpsi)) DEALLOCATE (dvpsi) - IF(ALLOCATED(dpsi)) DEALLOCATE ( dpsi) ! IF(ALLOCATED(vlocq)) DEALLOCATE (vlocq) IF(ALLOCATED(dmuxc)) DEALLOCATE (dmuxc) @@ -178,28 +175,25 @@ IF(ALLOCATED(int1_nc)) DEALLOCATE(int1_nc) IF(ALLOCATED(int3_nc)) DEALLOCATE(int3_nc) IF(ALLOCATED(int4_nc)) DEALLOCATE(int4_nc) - IF(ALLOCATED(becsum_nc)) DEALLOCATE(becsum_nc) - IF(ALLOCATED(alphasum_nc)) DEALLOCATE(alphasum_nc) IF(ALLOCATED(int2_so)) DEALLOCATE(int2_so) IF(ALLOCATED(int5_so)) DEALLOCATE(int5_so) - IF(ALLOCATED(alphasum)) DEALLOCATE (alphasum) ! - if(allocated(alphap)) then - do ik=1,nks - do ipol=1,3 - call deallocate_bec_type ( alphap(ipol,ik) ) - enddo - end do - deallocate (alphap) - endif - if(allocated(becp1)) then - do ik=1,size(becp1) - call deallocate_bec_type ( becp1(ik) ) - end do - deallocate(becp1) - end if - call deallocate_bec_type ( becp ) - + IF (allocated(alphap)) THEN + DO ik = 1, nks + DO ipol = 1, 3 + CALL deallocate_bec_type( alphap(ipol,ik) ) + ENDDO + ENDDO + DEALLOCATE(alphap) + ENDIF + IF (allocated(becp1)) THEN + DO ik = 1, size(becp1) + CALL deallocate_bec_type( becp1(ik) ) + ENDDO + DEALLOCATE(becp1) + ENDIF + CALL deallocate_bec_type ( becp ) + ! IF(ALLOCATED(nbnd_occ)) DEALLOCATE(nbnd_occ) IF(ALLOCATED(m_loc)) DEALLOCATE(m_loc) ! diff --git a/EPW/src/dvanqq2.f90 b/EPW/src/dvanqq2.f90 index 04803e52c..2bbf60f63 100644 --- a/EPW/src/dvanqq2.f90 +++ b/EPW/src/dvanqq2.f90 @@ -15,10 +15,16 @@ !! !! New !! This routine calculates two integrals of the Q functions and - !! its derivatives with c V_loc and V_eff which are used + !! its derivatives with V_loc and V_eff which are used !! to compute term dV_bare/dtau * psi in addusdvqpsi. !! The result is stored in int1, int2, int4, int5. The routine is called !! for each q in nqc. + !! int1 -> Eq. B20 of Ref.[1] + !! int2 -> Eq. B21 of Ref.[1] + !! int4 -> Eq. B23 of Ref.[1] + !! int5 -> Eq. B24 of Ref.[1] + !! + !! [1] PRB 64, 235118 (2001). !! !! RM - Nov/Dec 2014 !! Imported the noncolinear case implemented by xlzhang @@ -73,19 +79,27 @@ INTEGER :: is !! counter on spin ! - REAL(DP), ALLOCATABLE :: qmod(:) + REAL(kind=DP), ALLOCATABLE :: qmod(:) !! the modulus of q+G - REAL(DP), ALLOCATABLE :: qmodg(:) + REAL(kind=DP), ALLOCATABLE :: qmodg(:) !! the modulus of G REAL(DP), ALLOCATABLE :: qpg(:,:) !! the q+G vectors - REAL(DP), ALLOCATABLE :: ylmkq(:,:), ylmk0(:,:) - !! the spherical harmonics + REAL(kind=DP), ALLOCATABLE :: ylmkq(:,:) + !! the spherical harmonics at q+G + REAL(kind=DP), ALLOCATABLE :: ylmk0(:,:) + !! the spherical harmonics at G ! - COMPLEX(kind=DP) :: fact, fact1, ZDOTC - COMPLEX(kind=DP), ALLOCATABLE :: aux1(:), aux2(:),& - aux3(:), aux5(:), veff(:,:), sk(:) - ! work space + COMPLEX(kind=DP) :: fact + !! e^{-i q * \tau} * conjg(e^{-i q * \tau}) + COMPLEX(kind=DP) :: fact1 + !! -i * omega + COMPLEX(kind=DP), EXTERNAL :: zdotc + !! the scalar product function + COMPLEX(kind=DP), ALLOCATABLE :: aux1(:), aux2(:), & + aux3(:), aux5(:), sk(:) + COMPLEX(kind=DP), ALLOCATABLE :: veff(:,:) + !! effective potential COMPLEX(kind=DP), ALLOCATABLE, TARGET :: qgm(:) !! the augmentation function at G COMPLEX(kind=DP), POINTER :: qgmq(:) @@ -93,35 +107,35 @@ ! IF (.not.okvan) RETURN ! - CALL start_clock ('dvanqq2') + CALL start_clock('dvanqq2') ! int1(:,:,:,:,:) = czero int2(:,:,:,:,:) = czero int4(:,:,:,:,:) = czero int5(:,:,:,:,:) = czero - ALLOCATE (sk ( ngm)) - ALLOCATE (aux1( ngm)) - ALLOCATE (aux2( ngm)) - ALLOCATE (aux3( ngm)) - ALLOCATE (aux5( ngm)) - ALLOCATE (qmodg( ngm)) - ALLOCATE (ylmk0( ngm, lmaxq * lmaxq)) - ALLOCATE (qgm ( ngm)) - ALLOCATE (ylmkq( ngm, lmaxq * lmaxq)) - ALLOCATE (qmod( ngm)) - ALLOCATE (qgmq( ngm)) + ALLOCATE( sk(ngm) ) + ALLOCATE( aux1(ngm) ) + ALLOCATE( aux2(ngm) ) + ALLOCATE( aux3(ngm) ) + ALLOCATE( aux5(ngm) ) + ALLOCATE( qmodg(ngm) ) + ALLOCATE( qmod(ngm) ) + ALLOCATE( qgmq(ngm) ) + ALLOCATE( qgm(ngm)) + ALLOCATE( ylmk0(ngm, lmaxq * lmaxq) ) + ALLOCATE( ylmkq(ngm, lmaxq * lmaxq) ) ! ! compute spherical harmonics ! - CALL ylmr2 (lmaxq * lmaxq, ngm, g, gg, ylmk0) + CALL ylmr2( lmaxq * lmaxq, ngm, g, gg, ylmk0 ) DO ig = 1, ngm qmodg(ig) = sqrt( gg(ig) ) ENDDO ! - ALLOCATE (qpg(3, ngm)) - CALL setqmod(ngm, xq, g, qmod, qpg) + ALLOCATE( qpg(3, ngm) ) + CALL setqmod( ngm, xq, g, qmod, qpg ) CALL ylmr2(lmaxq * lmaxq, ngm, qpg, qmod, ylmkq) - DEALLOCATE (qpg) + DEALLOCATE(qpg) DO ig = 1, ngm qmod(ig) = sqrt( qmod(ig) ) ENDDO @@ -148,14 +162,15 @@ ! DO ntb = 1, ntyp IF (upf(ntb)%tvanp ) THEN + ! DO ih = 1, nh(ntb) DO jh = ih, nh(ntb) ijh = ijtoh(ih,jh,ntb) ! ! compute the augmentation function ! - CALL qvan2(ngm, ih, jh, ntb, qmodg, qgm, ylmk0) - CALL qvan2(ngm, ih, jh, ntb, qmod, qgmq, ylmkq) + CALL qvan2( ngm, ih, jh, ntb, qmodg, qgm, ylmk0 ) + CALL qvan2( ngm, ih, jh, ntb, qmod, qgmq, ylmkq ) ! ! NB: for this integral the moving atom and the atom of Q ! do not necessarily coincide @@ -167,8 +182,9 @@ * eigts2(mill(2,ig),nb) & * eigts3(mill(3,ig),nb) ENDDO + ! DO na = 1, nat - fact = eigqts(na) * CONJG( eigqts(nb) ) + fact = eigqts(na) * conjg( eigqts(nb) ) ! ! nb is the atom of the augmentation function ! @@ -178,12 +194,14 @@ * eigts2(mill(2,ig),na) & * eigts3(mill(3,ig),na) ENDDO + ! DO ipol = 1, 3 DO ig = 1, ngm aux5(ig) = sk(ig) * ( g(ipol,ig) + xq(ipol) ) ENDDO int2(ih,jh,ipol,na,nb) = fact * fact1 * & - ZDOTC(ngm, aux1, 1, aux5, 1) + zdotc(ngm, aux1, 1, aux5, 1) + ! DO jpol = 1, 3 IF (jpol >= ipol) THEN DO ig = 1, ngm @@ -191,45 +209,51 @@ ( g(jpol,ig) + xq(jpol) ) ENDDO int5(ijh,ipol,jpol,na,nb) = & - CONJG(fact) * tpiba2 * omega * & - ZDOTC(ngm, aux3, 1, aux1, 1) + conjg(fact) * tpiba2 * omega * & + zdotc(ngm, aux3, 1, aux1, 1) ELSE int5(ijh,ipol,jpol,na,nb) = & int5(ijh,jpol,ipol,na,nb) ENDIF ENDDO - ENDDO - ENDDO + ENDDO !ipol + ! + ENDDO !na + ! DO ig = 1, ngm aux1(ig) = qgm(ig) * eigts1(mill(1,ig),nb) & * eigts2(mill(2,ig),nb) & * eigts3(mill(3,ig),nb) ENDDO + ! DO is = 1, nspin_mag DO ipol = 1, 3 DO ig = 1, ngm aux2(ig) = veff(dfftp%nl(ig),is) * g(ipol,ig) ENDDO int1(ih,jh,ipol,nb,is) = - fact1 * & - ZDOTC(ngm, aux1, 1, aux2, 1) + zdotc(ngm, aux1, 1, aux2, 1) DO jpol = 1, 3 IF (jpol >= ipol) THEN DO ig = 1, ngm aux3(ig) = aux2(ig) * g(jpol,ig) ENDDO int4(ijh,ipol,jpol,nb,is) = - tpiba2 * & - omega * ZDOTC(ngm, aux3, 1, aux1, 1) + omega * zdotc(ngm, aux3, 1, aux1, 1) ELSE int4(ijh,ipol,jpol,nb,is) = & int4(ijh,jpol,ipol,nb,is) ENDIF - ENDDO - ENDDO - ENDDO - ENDIF - ENDDO - ENDDO - ENDDO + ENDDO ! jpol + ENDDO ! ipol + ENDDO ! is + ! + ENDIF ! ityp + ENDDO ! nb + ! + ENDDO ! jh + ENDDO ! ih + ! DO ih = 1, nh(ntb) DO jh = ih + 1, nh(ntb) ! @@ -247,14 +271,17 @@ ENDDO ENDIF ENDDO - ENDDO - ENDDO - ENDIF - ENDDO + ! + ENDDO ! jh + ENDDO ! ih + ! + ENDIF ! upf + ENDDO ! ntb CALL mp_sum(int1, intra_pool_comm) CALL mp_sum(int2, intra_pool_comm) CALL mp_sum(int4, intra_pool_comm) CALL mp_sum(int5, intra_pool_comm) + ! IF (noncolin) THEN CALL set_int12_nc(0) int4_nc = czero @@ -276,18 +303,29 @@ ENDDO ENDIF ! - DEALLOCATE (veff) - DEALLOCATE (qgmq) - DEALLOCATE (qmod) - DEALLOCATE (ylmkq) - DEALLOCATE (qgm) - DEALLOCATE (ylmk0) - DEALLOCATE (qmodg) - DEALLOCATE (aux5) - DEALLOCATE (aux3) - DEALLOCATE (aux2) - DEALLOCATE (aux1) - DEALLOCATE (sk) +!DBRM + !write(*,'(a,e20.12)') 'int1 = ', & + !SUM((REAL(REAL(int1(:,:,:,:,:))))**2)+SUM((REAL(AIMAG(int1(:,:,:,:,:))))**2) + !write(*,'(a,e20.12)') 'int2 = ', & + !SUM((REAL(REAL(int2(:,:,:,:,:))))**2)+SUM((REAL(AIMAG(int2(:,:,:,:,:))))**2) + !write(*,'(a,e20.12)') 'int4 = ', & + !SUM((REAL(REAL(int4(:,:,:,:,:))))**2)+SUM((REAL(AIMAG(int4(:,:,:,:,:))))**2) + !write(*,'(a,e20.12)') 'int5 = ', & + !SUM((REAL(REAL(int5(:,:,:,:,:))))**2)+SUM((REAL(AIMAG(int5(:,:,:,:,:))))**2) +!END + ! + DEALLOCATE(sk) + DEALLOCATE(aux1) + DEALLOCATE(aux2) + DEALLOCATE(aux3) + DEALLOCATE(aux5) + DEALLOCATE(qmodg) + DEALLOCATE(qmod) + DEALLOCATE(qgmq) + DEALLOCATE(qgm) + DEALLOCATE(ylmk0) + DEALLOCATE(ylmkq) + DEALLOCATE(veff) ! CALL stop_clock ('dvanqq2') RETURN diff --git a/EPW/src/dvqpsi_us3.f90 b/EPW/src/dvqpsi_us3.f90 index ce1a2d275..7240e2751 100644 --- a/EPW/src/dvqpsi_us3.f90 +++ b/EPW/src/dvqpsi_us3.f90 @@ -10,7 +10,7 @@ ! adapted from PH/dvqpsi_us (QE) ! !---------------------------------------------------------------------- - SUBROUTINE dvqpsi_us3( ik, uact, addnlcc, xxk, xq0 ) + SUBROUTINE dvqpsi_us3( ik, uact, addnlcc, xxkq, xq0 ) !---------------------------------------------------------------------- !! !! This routine calculates dV_bare/dtau * psi for one perturbation @@ -59,10 +59,10 @@ ! REAL(kind=DP), INTENT (in) :: xq0(3) !! Current coarse q-point coordinate - REAL(kind=DP), INTENT (in) :: xxk(3) - !! k-point coordinate + REAL(kind=DP), INTENT (in) :: xxkq(3) + !! k+q point coordinate ! - COMPLEX(kind=DP), INTENT(in) :: uact (3 * nat) + COMPLEX(kind=DP), INTENT(in) :: uact(3 * nat) !! the pattern of displacements ! ! Local variables @@ -86,14 +86,24 @@ INTEGER :: npw !! Number of k+G-vectors inside 'ecut sphere' ! - REAL(DP) :: fac + REAL(kind=DP) :: fac !! spin degeneracy factor ! - COMPLEX(DP) :: gtau, gu, fact, u1, u2, u3, gu0 - COMPLEX(DP), ALLOCATABLE, TARGET :: aux(:) - COMPLEX(DP), ALLOCATABLE :: aux1(:), aux2(:) - COMPLEX(DP), POINTER :: auxs(:) - COMPLEX(DP), ALLOCATABLE :: drhoc(:) + COMPLEX(kind=DP) :: gtau + !! e^{-i G * \tau} + COMPLEX(kind=DP) :: u1, u2, u3 + !! components of displacement pattern u + COMPLEX(kind=DP) :: gu0 + !! scalar product q * u + COMPLEX(kind=DP) :: gu + !! q * u + G * u + COMPLEX(kind=DP) :: fact + !! e^{-i q * \tau} + COMPLEX(kind=DP), ALLOCATABLE, TARGET :: aux(:) + COMPLEX(kind=DP), ALLOCATABLE :: aux1(:), aux2(:) + COMPLEX(kind=DP), POINTER :: auxs(:) + COMPLEX(kind=DP), ALLOCATABLE :: drhoc(:) + !! response core charge density ! CALL start_clock('dvqpsi_us3') ! @@ -115,131 +125,133 @@ dvpsi(:,:) = czero aux1(:) = czero DO na = 1, nat - fact = tpiba * (0.d0, -1.d0) * eigqts(na) - mu = 3 * (na - 1) - IF (abs(uact(mu+1)) + abs(uact(mu+2)) + abs(uact(mu+3)) .gt. eps12) THEN - nt = ityp(na) - u1 = uact(mu + 1) - u2 = uact(mu + 2) - u3 = uact(mu + 3) - gu0 = xq0(1) * u1 + xq0(2) * u2 + xq0(3) * u3 - DO ig = 1, ngms - gtau = eigts1(mill(1,ig),na) * eigts2(mill(2,ig),na) * eigts3(mill(3,ig),na) - gu = gu0 + g(1,ig) * u1 + g(2,ig) * u2 + g(3,ig) * u3 - aux1( dffts%nl(ig) ) = aux1( dffts%nl(ig) ) + vlocq(ig,nt) * gu * fact * gtau - ENDDO - ENDIF + fact = tpiba * (0.d0, -1.d0) * eigqts(na) + mu = 3 * (na - 1) + u1 = uact(mu+1) + u2 = uact(mu+2) + u3 = uact(mu+3) + IF (abs(u1) + abs(u2) + abs(u3) .gt. eps12) THEN + nt = ityp(na) + gu0 = xq0(1) * u1 + xq0(2) * u2 + xq0(3) * u3 + DO ig = 1, ngms + gtau = eigts1(mill(1,ig),na) * & + eigts2(mill(2,ig),na) * & + eigts3(mill(3,ig),na) + gu = gu0 + g(1,ig) * u1 + g(2,ig) * u2 + g(3,ig) * u3 + aux1(dffts%nl(ig)) = aux1(dffts%nl(ig)) + vlocq(ig,nt) * gu * fact * gtau + ENDDO + ENDIF ENDDO ! ! add NLCC when present ! - IF (nlcc_any .AND. addnlcc) THEN - drhoc(:) = czero - DO na = 1, nat - fact = tpiba * (0.d0, -1.d0) * eigqts(na) - mu = 3 * (na - 1) - IF (abs(uact(mu+1)) + abs(uact(mu+2)) + abs(uact(mu+3)) .gt. eps12) THEN - nt = ityp(na) - u1 = uact(mu+1) - u2 = uact(mu+2) - u3 = uact(mu+3) - gu0 = xq0(1) * u1 + xq0(2) * u2 + xq0(3) * u3 - IF (upf(nt)%nlcc) THEN - DO ig = 1,ngm - gtau = eigts1(mill(1,ig),na) * & - eigts2(mill(2,ig),na) * & - eigts3(mill(3,ig),na) - gu = gu0 + g(1,ig) * u1 + g(2,ig) * u2 + g(3,ig) * u3 - drhoc(dfftp%nl(ig)) = drhoc(dfftp%nl(ig)) + drc(ig,nt) * gu * fact * gtau - ENDDO - ENDIF - ENDIF - ENDDO - ! - CALL invfft('Rho', drhoc, dfftp) - ! - IF (.not.lsda) THEN - DO ir = 1, dfftp%nnr - aux(ir) = drhoc(ir) * dmuxc(ir,1,1) - ENDDO - ELSE - is = isk(ik) - DO ir = 1, dfftp%nnr - aux(ir) = drhoc(ir) * 0.5d0 * ( dmuxc(ir,is,1) + dmuxc(ir,is,2) ) - ENDDO + IF (nlcc_any .AND. addnlcc) THEN + drhoc(:) = czero + DO na = 1, nat + fact = tpiba * (0.d0, -1.d0) * eigqts(na) + mu = 3 * (na - 1) + u1 = uact(mu+1) + u2 = uact(mu+2) + u3 = uact(mu+3) + IF (abs(u1) + abs(u2) + abs(u3) .gt. eps12) THEN + nt = ityp(na) + gu0 = xq0(1) * u1 + xq0(2) * u2 + xq0(3) * u3 + IF (upf(nt)%nlcc) THEN + DO ig = 1,ngm + gtau = eigts1(mill(1,ig),na) * & + eigts2(mill(2,ig),na) * & + eigts3(mill(3,ig),na) + gu = gu0 + g(1,ig) * u1 + g(2,ig) * u2 + g(3,ig) * u3 + drhoc(dfftp%nl(ig)) = drhoc(dfftp%nl(ig)) + drc(ig,nt) * gu * fact * gtau + ENDDO + ENDIF ENDIF - ! - fac = 1.d0 / dble(nspin_lsda) - DO is = 1, nspin_lsda - rho%of_r(:,is) = rho%of_r(:,is) + fac * rho_core + ENDDO + ! + CALL invfft('Rho', drhoc, dfftp) + ! + IF (.not.lsda) THEN + DO ir = 1, dfftp%nnr + aux(ir) = drhoc(ir) * dmuxc(ir,1,1) ENDDO - ! - IF ( dft_is_gradient() ) & - CALL dgradcorr( dfftp, rho%of_r, grho, & - dvxc_rr, dvxc_sr, dvxc_ss, dvxc_s, xq0, drhoc, & - 1, nspin_gga, g, aux ) - ! - IF ( dft_is_nonlocc() ) & - CALL dnonloccorr( rho%of_r, drhoc, xq0, aux ) - ! - DO is = 1, nspin_lsda - rho%of_r(:,is) = rho%of_r(:,is) - fac * rho_core - END DO - ! - CALL fwfft('Rho', aux, dfftp) - ! - ! This is needed also when the smooth and the thick grids coincide to - ! cut the potential at the cut-off - ! - auxs(:) = czero - DO ig = 1, ngms - auxs(dffts%nl(ig)) = aux(dfftp%nl(ig)) + ELSE + is = isk(ik) + DO ir = 1, dfftp%nnr + aux(ir) = drhoc(ir) * 0.5d0 * ( dmuxc(ir,is,1) + dmuxc(ir,is,2) ) ENDDO - aux1(:) = aux1(:) + auxs(:) - ENDIF + ENDIF + ! + fac = 1.d0 / dble(nspin_lsda) + DO is = 1, nspin_lsda + rho%of_r(:,is) = rho%of_r(:,is) + fac * rho_core + ENDDO + ! + IF ( dft_is_gradient() ) & + CALL dgradcorr( dfftp, rho%of_r, grho, & + dvxc_rr, dvxc_sr, dvxc_ss, dvxc_s, xq0, drhoc, & + 1, nspin_gga, g, aux ) + ! + IF ( dft_is_nonlocc() ) & + CALL dnonloccorr( rho%of_r, drhoc, xq0, aux ) + ! + DO is = 1, nspin_lsda + rho%of_r(:,is) = rho%of_r(:,is) - fac * rho_core + ENDDO + ! + CALL fwfft('Rho', aux, dfftp) + ! + ! This is needed also when the smooth and the thick grids coincide to + ! cut the potential at the cut-off + ! + auxs(:) = czero + DO ig = 1, ngms + auxs(dffts%nl(ig)) = aux(dfftp%nl(ig)) + ENDDO + aux1(:) = aux1(:) + auxs(:) + ENDIF ! ! Now we compute dV_loc/dtau in real space ! CALL invfft('Rho', aux1, dffts) DO ibnd = lower_band, upper_band - DO ip = 1, npol - aux2(:) = czero - IF ( ip == 1 ) THEN - DO ig = 1, npw - aux2( dffts%nl( igk(ig) ) ) = evc(ig,ibnd) - ENDDO - ELSE - DO ig = 1, npw - aux2( dffts%nl( igk(ig) ) ) = evc(ig+npwx,ibnd) - ENDDO - ENDIF - ! - ! This wavefunction is computed in real space - ! - CALL invfft('Wave', aux2, dffts) - DO ir = 1, dffts%nnr - aux2(ir) = aux2(ir) * aux1(ir) + DO ip = 1, npol + aux2(:) = czero + IF ( ip == 1 ) THEN + DO ig = 1, npw + aux2( dffts%nl( igk(ig) ) ) = evc(ig,ibnd) ENDDO - ! - ! and finally dV_loc/dtau * psi is transformed in reciprocal space - ! - CALL fwfft('Wave', aux2, dffts) - IF ( ip == 1 ) THEN - DO ig = 1, npwq - dvpsi(ig,ibnd) = aux2( dffts%nl( igkq(ig) ) ) - ENDDO - ELSE - DO ig = 1, npwq - dvpsi(ig+npwx,ibnd) = aux2( dffts%nl( igkq(ig) ) ) - ENDDO - ENDIF - ENDDO + ELSE + DO ig = 1, npw + aux2( dffts%nl( igk(ig) ) ) = evc(ig+npwx,ibnd) + ENDDO + ENDIF + ! + ! This wavefunction is computed in real space + ! + CALL invfft('Wave', aux2, dffts) + DO ir = 1, dffts%nnr + aux2(ir) = aux2(ir) * aux1(ir) + ENDDO + ! + ! and finally dV_loc/dtau * psi is transformed in reciprocal space + ! + CALL fwfft('Wave', aux2, dffts) + IF ( ip == 1 ) THEN + DO ig = 1, npwq + dvpsi(ig,ibnd) = aux2( dffts%nl( igkq(ig) ) ) + ENDDO + ELSE + DO ig = 1, npwq + dvpsi(ig+npwx,ibnd) = aux2( dffts%nl( igkq(ig) ) ) + ENDDO + ENDIF + ENDDO ENDDO ! IF (nlcc_any .AND. addnlcc) THEN - DEALLOCATE(drhoc) - DEALLOCATE(aux) - DEALLOCATE(auxs) + DEALLOCATE(drhoc) + DEALLOCATE(aux) + DEALLOCATE(auxs) ENDIF DEALLOCATE(aux1) DEALLOCATE(aux2) @@ -248,7 +260,7 @@ ! First a term similar to the KB case. ! Then a term due to the change of the D coefficients in the perturbat ! - CALL dvqpsi_us_only3( ik, uact, xxk ) + CALL dvqpsi_us_only3( ik, uact, xxkq ) ! CALL stop_clock('dvqpsi_us3') ! diff --git a/EPW/src/dvqpsi_us_only3.f90 b/EPW/src/dvqpsi_us_only3.f90 index b8eb54ce0..674985e28 100644 --- a/EPW/src/dvqpsi_us_only3.f90 +++ b/EPW/src/dvqpsi_us_only3.f90 @@ -10,7 +10,7 @@ ! adapted from PH/dvqpsi_us_only (QE) ! !---------------------------------------------------------------------- - subroutine dvqpsi_us_only3( ik, uact, xxk ) + subroutine dvqpsi_us_only3( ik, uact, xxkq ) !---------------------------------------------------------------------- !! !! This routine calculates dV_bare/dtau * psi for one perturbation @@ -41,11 +41,11 @@ IMPLICIT NONE ! INTEGER, INTENT(in) :: ik - !! Input: the k point - REAL(kind=DP), INTENT(in) :: xxk(3) - !! input: the k point (cartesian coordinates) - COMPLEX(kind=DP), INTENT(in) :: uact (3 * nat) - !! input: the pattern of displacements + !! the k point + REAL(kind=DP), INTENT(in) :: xxkq(3) + !! the k+q point (cartesian coordinates) + COMPLEX(kind=DP), INTENT(in) :: uact(3 * nat) + !! the pattern of displacements ! ! Local variables ! @@ -84,10 +84,10 @@ INTEGER :: ijs !! Counter on combined is and js polarization ! - REAL(DP), ALLOCATABLE :: deff(:,:,:) + REAL(kind=DP), ALLOCATABLE :: deff(:,:,:) ! - COMPLEX(DP), ALLOCATABLE :: ps1(:,:), ps2(:,:,:), aux(:), deff_nc(:,:,:,:) - COMPLEX(DP), ALLOCATABLE :: ps1_nc(:,:,:), ps2_nc(:,:,:,:) + COMPLEX(kind=DP), ALLOCATABLE :: ps1(:,:), ps2(:,:,:), aux(:), deff_nc(:,:,:,:) + COMPLEX(kind=DP), ALLOCATABLE :: ps1_nc(:,:,:), ps2_nc(:,:,:,:) ! LOGICAL :: ok ! @@ -124,16 +124,14 @@ ijkb0 = 0 DO nt = 1, ntyp DO na = 1, nat - IF (ityp (na) .eq. nt) THEN + IF (ityp(na) .eq. nt) THEN mu = 3 * (na - 1) DO ih = 1, nh(nt) ikb = ijkb0 + ih DO jh = 1, nh(nt) jkb = ijkb0 + jh DO ipol = 1, 3 - IF ( abs( uact(mu + 1) ) + & - abs( uact(mu + 2) ) + & - abs( uact(mu + 3) ) > eps12 ) THEN + IF ( abs( uact(mu+1) ) + abs( uact(mu+2) ) + abs( uact(mu+3) ) > eps12 ) THEN IF (noncolin) THEN ijs = 0 DO is = 1, npol @@ -155,52 +153,52 @@ deff(ih,jh,na) * becp1(ik)%k(jkb,ibnd) * & (0.d0,-1.d0) * uact(mu+ipol) * tpiba ENDIF - IF (okvan) THEN - IF (noncolin) THEN - ijs = 0 - DO is = 1, npol - DO js = 1, npol - ijs = ijs + 1 - ps1_nc(ikb,is,ibnd) = ps1_nc(ikb,is,ibnd) + & - int1_nc(ih,jh,ipol,na,ijs) * & - becp1(ik)%nc(jkb,js,ibnd) * uact(mu+ipol) - ENDDO - ENDDO - ELSE - ps1(ikb,ibnd) = ps1(ikb, ibnd) + & - int1(ih,jh,ipol,na,current_spin) * & - becp1(ik)%k(jkb,ibnd) * uact(mu +ipol) - ENDIF - ENDIF +! IF (okvan) THEN +! IF (noncolin) THEN +! ijs = 0 +! DO is = 1, npol +! DO js = 1, npol +! ijs = ijs + 1 +! ps1_nc(ikb,is,ibnd) = ps1_nc(ikb,is,ibnd) + & +! int1_nc(ih,jh,ipol,na,ijs) * & +! becp1(ik)%nc(jkb,js,ibnd) * uact(mu+ipol) +! ENDDO +! ENDDO +! ELSE +! ps1(ikb,ibnd) = ps1(ikb, ibnd) + & +! int1(ih,jh,ipol,na,current_spin) * & +! becp1(ik)%k(jkb,ibnd) * uact(mu+ipol) +! ENDIF +! ENDIF ! okvan ENDIF ! uact>0 - IF (okvan) THEN - DO nb = 1, nat - nu = 3 * (nb - 1) - IF (noncolin) THEN - IF (lspinorb) THEN - ijs = 0 - DO is = 1, npol - DO js = 1, npol - ijs = ijs + 1 - ps1_nc(ikb,is,ibnd) = ps1_nc(ikb,is,ibnd) + & - int2_so(ih,jh,ipol,nb,na,ijs) * & - becp1(ik)%nc(jkb,js,ibnd) * uact(nu+ipol) - ENDDO - ENDDO - ELSE - DO is = 1, npol - ps1_nc(ikb,is,ibnd) = ps1_nc(ikb,is,ibnd) + & - int2(ih,jh,ipol,nb,na) * & - becp1(ik)%nc(jkb,is,ibnd) * uact(nu+ipol) - ENDDO - ENDIF - ELSE - ps1(ikb,ibnd) = ps1(ikb,ibnd) + & - int2(ih,jh,ipol,nb,na) * & - becp1(ik)%k(jkb,ibnd) * uact(nu+ipol) - ENDIF - ENDDO - ENDIF ! okvan +! IF (okvan) THEN +! DO nb = 1, nat +! nu = 3 * (nb - 1) +! IF (noncolin) THEN +! IF (lspinorb) THEN +! ijs = 0 +! DO is = 1, npol +! DO js = 1, npol +! ijs = ijs + 1 +! ps1_nc(ikb,is,ibnd) = ps1_nc(ikb,is,ibnd) + & +! int2_so(ih,jh,ipol,nb,na,ijs) * & +! becp1(ik)%nc(jkb,js,ibnd) * uact(nu+ipol) +! ENDDO +! ENDDO +! ELSE +! DO is = 1, npol +! ps1_nc(ikb,is,ibnd) = ps1_nc(ikb,is,ibnd) + & +! int2(ih,jh,ipol,nb,na) * & +! becp1(ik)%nc(jkb,is,ibnd) * uact(nu+ipol) +! ENDDO +! ENDIF +! ELSE +! ps1(ikb,ibnd) = ps1(ikb,ibnd) + & +! int2(ih,jh,ipol,nb,na) * & +! becp1(ik)%k(jkb,ibnd) * uact(nu+ipol) +! ENDIF +! ENDDO +! ENDIF ! okvan ENDDO ! ipol ENDDO ! jh ENDDO ! ih @@ -241,7 +239,7 @@ DO ig = 1, npwq igg = igkq(ig) !aux(ig) = vkb(ig,ikb) * ( xk(ipol,ikq) + g(ipol,igg) ) - aux(ig) = vkb(ig,ikb) * ( xxk(ipol) + g(ipol,igg) ) + aux(ig) = vkb(ig,ikb) * ( xxkq(ipol) + g(ipol,igg) ) ENDDO DO ibnd = lower_band, upper_band IF (noncolin) THEN diff --git a/EPW/src/elphel2_shuffle.f90 b/EPW/src/elphel2_shuffle.f90 index fc00e1734..4826858f5 100644 --- a/EPW/src/elphel2_shuffle.f90 +++ b/EPW/src/elphel2_shuffle.f90 @@ -154,23 +154,23 @@ !! Absolute index of k+q-point ! ! Local variables for rotating the wavefunctions (in order to use q in the irr wedge) - REAL(DP) :: xktmp(3) + REAL(kind=DP) :: xkqtmp(3) !! Temporary k+q vector for KB projectors - REAL(DP) :: sxk(3) - !! - REAL(DP) :: g0vec_all_r(3,125) + REAL(kind=DP) :: sxk(3) + !! Rotated k-point xk + REAL(kind=DP) :: g0vec_all_r(3,125) !! G_0 vectors needed to fold the k+q grid into the k grid, cartesian coord. - REAL(DP) :: zero_vect(3) + REAL(kind=DP) :: zero_vect(3) !! Temporary zero vector ! - COMPLEX(DP), ALLOCATABLE :: aux1(:,:), aux2(:,:), aux3(:,:) - COMPLEX(DP), ALLOCATABLE :: eptmp(:,:), elphmat(:,:,:) + COMPLEX(kind=DP), ALLOCATABLE :: aux1(:,:), aux2(:,:), aux3(:,:) + COMPLEX(kind=DP), ALLOCATABLE :: eptmp(:,:), elphmat(:,:,:) !! arrays for e-ph matrix elements ! !DBSP - NAG complains ... - COMPLEX(DP), EXTERNAL :: ZDOTC + COMPLEX(kind=DP), EXTERNAL :: zdotc !DBSP -! REAL(kind=DP) :: b,c +! REAL(kind=DP) :: b, c, d !END ! IF ( .not. ALLOCATED(elphmat) ) ALLOCATE( elphmat(nbnd, nbnd, npe) ) @@ -221,8 +221,9 @@ IF (lsda) current_spin = isk(ik) elphmat(:,:,:) = czero !DBSP -! c = 0 -! b = 0 +! b = zero +! c = zero +! d = zero !END ! ! find index, and possibly pool, of k+q @@ -285,7 +286,7 @@ umatq(:,:,ik) = umat_all(:,:,nkq_abs) ! ! the k-vector needed for the KB projectors - xktmp = xkq(:,ik) + xkqtmp = xkq(:,ik) ! ! -------------------------------------------------- ! Fourier translation of the G-sphere igkq @@ -306,7 +307,7 @@ ! (this is needed in the calculation of the KB terms ! for nonlocal pseudos) ! - xktmp = xkq(:,ik) - g0vec_all_r(:,shift(ik+ik0)) + xkqtmp = xkq(:,ik) - g0vec_all_r(:,shift(ik+ik0)) ! ! --------------------------------------------------------------------- ! phase factor arising from fractional traslations @@ -371,8 +372,8 @@ ! now we generate vkb on the igkq() set because dvpsi is needed on that set ! we need S(k)+q_0 in the KB projector: total momentum transfer must be q_0 ! - xktmp = sxk + xq0 - CALL init_us_2( npwq, igkq, xktmp, vkb ) + xkqtmp = sxk + xq0 + CALL init_us_2( npwq, igkq, xkqtmp, vkb ) ! ! -------------------------------------------------- ! Calculation of the matrix element @@ -382,15 +383,15 @@ ! ! recalculate dvbare_q*psi_k ! the call to dvqpsi_us3 differs from the old one to dvqpsi_us - ! only the xktmp passed. + ! only the xkqtmp passed. ! ! we have to use the first q in the star in the dvqpsi_us3 call below (xq0) ! mode = imode0 + ipert IF (timerev) THEN - CALL dvqpsi_us3( ik, CONJG(u(:,mode)), .false., xktmp, xq0 ) + CALL dvqpsi_us3( ik, conjg(u(:,mode)), .false., xkqtmp, xq0 ) ELSE - CALL dvqpsi_us3( ik, u(:,mode), .false., xktmp, xq0 ) + CALL dvqpsi_us3( ik, u(:,mode), .false., xkqtmp, xq0 ) ENDIF !DBSP ! b = b+SUM((REAL(REAL(dvpsi(:,:))))**2)+SUM((REAL(AIMAG(dvpsi(:,:))))**2) @@ -404,7 +405,7 @@ DO ibnd = lower_band, upper_band CALL invfft_wave(npw, igk, evc(:,ibnd), aux1) IF (timerev) THEN - CALL apply_dpot(dffts%nnr, aux1, CONJG(dvscfins(:,:,ipert)), current_spin) + CALL apply_dpot(dffts%nnr, aux1, conjg(dvscfins(:,:,ipert)), current_spin) ELSE CALL apply_dpot(dffts%nnr, aux1, dvscfins(:,:,ipert), current_spin) ENDIF @@ -412,10 +413,13 @@ ENDDO dvpsi = dvpsi + aux3 ! - CALL adddvscf2( ipert, ik ) - ! !DBSP - !c = c+SUM((REAL(REAL(dvpsi(:,:))))**2)+SUM((REAL(AIMAG(dvpsi(:,:))))**2) +! c = c+SUM((REAL(REAL(dvpsi(:,:))))**2)+SUM((REAL(AIMAG(dvpsi(:,:))))**2) +!END + ! + CALL adddvscf2( ipert, ik ) +!DBRM +! d = c+SUM((REAL(REAL(dvpsi(:,:))))**2)+SUM((REAL(AIMAG(dvpsi(:,:))))**2) !END ! ! calculate elphmat(j,i)= for this pertur @@ -424,10 +428,10 @@ DO ibnd =lower_band, upper_band DO jbnd = 1, nbnd elphmat(jbnd,ibnd,ipert) = & - ZDOTC( npwq, evq(1,jbnd), 1, dvpsi(1,ibnd), 1 ) + zdotc( npwq, evq(1,jbnd), 1, dvpsi(1,ibnd), 1 ) IF (noncolin) & elphmat(jbnd,ibnd,ipert) = elphmat(jbnd,ibnd,ipert) + & - ZDOTC( npwq, evq(npwx+1,jbnd), 1, dvpsi(npwx+1,ibnd), 1 ) + zdotc( npwq, evq(npwx+1,jbnd), 1, dvpsi(npwx+1,ibnd), 1 ) ENDDO ENDDO ENDDO @@ -439,7 +443,8 @@ ! IF (ik==2) THEN ! write(*,*)'SUM dvpsi b ', b ! write(*,*)'SUM dvpsi c ', c -! write(*,*)'elphmat(:,:,:)**2',SUM((REAL(REAL(elphmat(:,:,:))))**2)+SUM((REAL(AIMAG(elphmat(:,:,:))))**2) +! write(*,*)'SUM dvpsi d ', d +! write(*,*)'elphmat(:,:,:)**2', SUM((REAL(REAL(elphmat(:,:,:))))**2)+SUM((REAL(AIMAG(elphmat(:,:,:))))**2) ! ENDIF !END ! @@ -498,8 +503,8 @@ IMPLICIT NONE ! INTEGER, INTENT(in) :: npw, igk(npwx) - COMPLEX(DP), INTENT(inout) :: evc(npwx*npol, nbnd) - COMPLEX(DP), INTENT(in) :: eigv1(ngm), eig0v + COMPLEX(kind=DP), INTENT(inout) :: evc(npwx*npol, nbnd) + COMPLEX(kind=DP), INTENT(in) :: eigv1(ngm), eig0v ! INTEGER :: ig !! Counter on G-vectors @@ -529,11 +534,15 @@ ! IMPLICIT NONE ! + REAL(kind=DP), INTENT(in) :: x(3) + !! Input x INTEGER, INTENT(in) :: s(3,3) - REAL(DP), INTENT(in) :: x(3) - REAL(DP), INTENT(out) :: sx(3) + !! Symmetry matrix + REAL(kind=DP), INTENT(out) :: sx(3) + !! Output rotated x ! REAL(DP) :: xcrys(3) + !! x in cartesian coords INTEGER :: i ! xcrys = x diff --git a/EPW/src/elphon_shuffle.f90 b/EPW/src/elphon_shuffle.f90 index 9ccde937b..123161a36 100644 --- a/EPW/src/elphon_shuffle.f90 +++ b/EPW/src/elphon_shuffle.f90 @@ -13,8 +13,6 @@ !! Electron-phonon calculation from data saved in fildvscf !! Shuffle2 mode (shuffle on electrons + load all phonon q's) !! - !! No ultrasoft yet - !! !! RM - Nov/Dec 2014 !! Imported the noncolinear case implemented by xlzhang !! @@ -24,9 +22,7 @@ ! USE kinds, ONLY : DP USE mp, ONLY : mp_barrier, mp_sum - USE mp_global, ONLY : my_pool_id, nproc_pool, npool, kunit, & - inter_pool_comm - USE mp_images, ONLY : nproc_image + USE mp_global, ONLY : my_pool_id, npool, inter_pool_comm USE ions_base, ONLY : nat USE pwcom, ONLY : nbnd, nks, nkstot USE gvect, ONLY : ngm @@ -79,47 +75,14 @@ !! Counter on bands INTEGER :: jbnd !! Counter on bands - INTEGER :: tmp_pool_id - !! temporary pool id - INTEGER :: iks - !! Index of the first k-point block in this pool - INTEGER :: ik0 - !! Index of iks - 1 - INTEGER :: nkl - !! - INTEGER :: nkr - !! ! - COMPLEX(DP), POINTER :: dvscfin(:,:,:) + COMPLEX(kind=DP), POINTER :: dvscfin(:,:,:) !! Change of the scf potential - COMPLEX(DP), POINTER :: dvscfins(:,:,:) + COMPLEX(kind=DP), POINTER :: dvscfins(:,:,:) !! Change of the scf potential (smooth) ! CALL start_clock('elphon_shuffle') ! - ik0 = 0 - tmp_pool_id = 0 - ! - npool = nproc_image / nproc_pool - IF (npool.gt.1) THEN - ! - ! number of kpoint blocks, kpoints per pool and reminder - kunit = 1 - nkl = kunit * ( nkstot / npool ) - nkr = ( nkstot - nkl * npool ) / kunit - ! the reminder goes to the first nkr pools - IF ( my_pool_id < nkr ) nkl = nkl + kunit - ! - iks = nkl * my_pool_id + 1 - IF ( my_pool_id >= nkr ) iks = iks + nkr * kunit - ! - ! the index of the first k point block in this pool - 1 - ! (I will need the index of ik, not ikk) - ! - ik0 = ( iks - 1 ) / kunit - ! - ENDIF - ! ! read Delta Vscf and calculate electron-phonon coefficients ! imode0 = 0 @@ -152,7 +115,7 @@ dvscfins => dvscfin ENDIF ! - CALL newdq2( dvscfin, npe ) + CALL newdq2( dvscfin, npe, xq0, timerev ) CALL elphel2_shuffle( npe, imode0, dvscfins, gmapsym, eigv, isym, xq0, timerev ) ! imode0 = imode0 + npe @@ -170,8 +133,8 @@ ! must be transformed in the cartesian basis ! epmat_{CART} = conjg ( U ) * epmat_{PATTERN} ! - ! note it is not U^\dagger ! Have a look to symdyn_munu.f90 - ! for comparison + ! note it is not U^\dagger but u_pattern! + ! Have a look to symdyn_munu.f90 for comparison ! DO ibnd = 1, nbnd DO jbnd = 1, nbnd diff --git a/EPW/src/elphon_shuffle_wrap.f90 b/EPW/src/elphon_shuffle_wrap.f90 index 97e3a8347..9d5700a56 100644 --- a/EPW/src/elphon_shuffle_wrap.f90 +++ b/EPW/src/elphon_shuffle_wrap.f90 @@ -455,15 +455,12 @@ ! minus_q = (iswitch .gt. -3) ! -! ! Initialize vlocq for the current irreducible q-point -! CALL epw_init(.false.) - ! ! loop over the q points of the star ! DO iq = 1, nq ! SP: First the vlocq needs to be initialized properly with the first ! q in the star - xq = xq0 + xq = xq0 CALL epw_init(.false.) ! ! retrieve the q in the star @@ -655,15 +652,14 @@ ! CALL createkmap( xq ) ! - xq0 = -xq0 - ! CALL loadumat( nbnd, nbndsub, nks, nkstot, xq, cu, cuq, lwin, lwinq, exband, w_centers ) ! ! Calculate overlap U_k+q U_k^\dagger IF (lpolar) CALL compute_umn_c( nbnd, nbndsub, nks, cu, cuq, bmat(:,:,:,nqc) ) ! - CALL elphon_shuffle( iq_irr, nqc_irr, nqc, gmapsym, eigv, isym, xq0, .true. ) + xq0 = -xq0 ! + CALL elphon_shuffle( iq_irr, nqc_irr, nqc, gmapsym, eigv, isym, xq0, .true. ) ! bring epmatq in the mode representation of iq_first, ! and then in the cartesian representation of iq ! @@ -795,20 +791,20 @@ ! IMPLICIT NONE ! - REAL(DP), INTENT(in):: x(3) + REAL(kind=DP), INTENT(in) :: x(3) !! Input x INTEGER, INTENT(in) :: s(3,3) !! Symmetry matrix - REAL(DP), INTENT(out):: sx(3) - !! Output rotated one. + REAL(kind=DP), INTENT(out) :: sx(3) + !! Output rotated x ! ! Local Variable INTEGER :: i ! DO i = 1, 3 - sx (i) = dble( s(i,1) ) * x(1) & - + dble( s(i,2) ) * x(2) & - + dble( s(i,3) ) * x(3) + sx(i) = dble( s(i,1) ) * x(1) & + + dble( s(i,2) ) * x(2) & + + dble( s(i,3) ) * x(3) ENDDO ! RETURN @@ -824,12 +820,12 @@ ! IMPLICIT NONE ! - REAL(DP), INTENT(in) :: x(3) + REAL(kind=DP), INTENT(in) :: x(3) !! input: input vector - REAL(DP), INTENT(in) :: y(3) + REAL(kind=DP), INTENT(in) :: y(3) !! input: second input vector - REAL(DP) :: accep - !! acceptance PARAMETER + REAL(kind=DP) :: accep + !! acceptance parameter PARAMETER (accep = 1.0d-5) ! eqvect_strict = abs( x(1)-y(1) ) .lt. accep .AND. & @@ -853,13 +849,22 @@ ! IMPLICIT NONE ! - INTEGER, INTENT(IN) :: current_iq + INTEGER, INTENT(in) :: current_iq !! Current q-point - INTEGER, INTENT(IN) :: iunpun + INTEGER, INTENT(in) :: iunpun !! Current q-point - INTEGER, INTENT(OUT) :: ierr + INTEGER, INTENT(out) :: ierr !! Error - INTEGER :: imode0, imode, irr, ipert, iq + ! + ! Local variables + INTEGER :: imode0, imode + !! Counter on modes + INTEGER :: irr + !! Counter on irreducible representations + INTEGER :: ipert + !! Counter on perturbations at each irr + INTEGER :: iq + !! Current q-point ! ierr = 0 IF (meta_ionode) THEN diff --git a/EPW/src/newdq2.f90 b/EPW/src/newdq2.f90 index 53474ddd5..0e4f30516 100644 --- a/EPW/src/newdq2.f90 +++ b/EPW/src/newdq2.f90 @@ -11,7 +11,7 @@ ! Adapted from LR_Modules/newdq.f90 (QE) ! !---------------------------------------------------------------------- - SUBROUTINE newdq2( dvscf, npe ) + SUBROUTINE newdq2( dvscf, npe, xq0, timerev ) !---------------------------------------------------------------------- !! !! This routine computes the contribution of the selfconsistent @@ -32,15 +32,19 @@ USE mp_global, ONLY : intra_pool_comm USE mp, ONLY : mp_sum USE lrus, ONLY : int3 - USE qpoint, ONLY : xq, eigqts + USE qpoint, ONLY : eigqts USE constants_epw, ONLY : czero ! IMPLICIT NONE ! INTEGER, INTENT(in) :: npe !! Number of perturbations for this irr representation - COMPLEX(DP), INTENT(in) :: dvscf(dfftp%nnr, nspin_mag, npe) + REAL(kind=DP), INTENT(in) :: xq0(3) + !! The first q-point in the star (cartesian coords.) + COMPLEX(kind=DP), INTENT(in) :: dvscf(dfftp%nnr, nspin_mag, npe) !! Change of the selfconsistent potential + LOGICAL, INTENT(in) :: timerev + !! true if we are using time reversal ! ! Local variables ! @@ -61,16 +65,20 @@ INTEGER :: jh !! Counter on beta functions ! - REAL(DP), ALLOCATABLE :: qmod(:) + REAL(kind=DP), ALLOCATABLE :: qmod(:) !! the modulus of q+G - REAL(DP), ALLOCATABLE :: qg(:,:) + REAL(kind=DP), ALLOCATABLE :: qg(:,:) !! the values of q+G - REAL(DP), ALLOCATABLE :: ylmk0(:,:) - !! the spherical harmonics + REAL(kind=DP), ALLOCATABLE :: ylmk0(:,:) + !! the spherical harmonics at q+G ! - COMPLEX(DP), EXTERNAL :: zdotc + COMPLEX(kind=DP), EXTERNAL :: zdotc !! the scalar product function - COMPLEX(DP), ALLOCATABLE :: aux1(:), aux2(:,:), veff(:), qgm(:) + COMPLEX(kind=DP), ALLOCATABLE :: aux1(:), aux2(:,:) + COMPLEX(kind=DP), ALLOCATABLE :: qgm(:) + !! the augmentation function at G + COMPLEX(kind=DP), ALLOCATABLE :: veff(:) + !! effective potential ! IF (.not.okvan) RETURN ! @@ -88,7 +96,7 @@ ! ! first compute the spherical harmonics ! - CALL setqmod( ngm, xq, g, qmod, qg ) + CALL setqmod( ngm, xq0, g, qmod, qg ) CALL ylmr2( lmaxq * lmaxq, ngm, qg, qmod, ylmk0 ) DO ig = 1, ngm qmod(ig) = sqrt( qmod(ig) ) @@ -102,19 +110,25 @@ ! DO is = 1, nspin_mag DO ir = 1, dfftp%nnr - veff(ir) = dvscf(ir,is,ipert) + IF (timerev) THEN + veff(ir) = conjg(dvscf(ir,is,ipert)) + ELSE + veff(ir) = dvscf(ir,is,ipert) + ENDIF ENDDO CALL fwfft('Rho', veff, dfftp) DO ig = 1, ngm - aux2(ig,is) = veff( dfftp%nl(ig) ) + aux2(ig,is) = veff( dfftp%nl(ig) ) ENDDO ENDDO ! DO nt = 1, ntyp IF (upf(nt)%tvanp ) THEN + ! DO ih = 1, nh(nt) DO jh = ih, nh(nt) CALL qvan2( ngm, ih, jh, nt, qmod, qgm, ylmk0 ) + ! DO na = 1, nat IF (ityp(na) == nt) THEN DO ig = 1, ngm @@ -125,12 +139,14 @@ ENDDO DO is = 1, nspin_mag int3(ih,jh,na,is,ipert) = omega * & - zdotc(ngm,aux1,1,aux2(1,is),1) + zdotc(ngm,aux1,1,aux2(1,is),1) ENDDO ENDIF ENDDO - ENDDO - ENDDO + ! + ENDDO ! jh + ENDDO ! ih + ! DO na = 1, nat IF (ityp(na) == nt) THEN ! @@ -143,23 +159,31 @@ ENDDO ENDDO ENDDO - ENDIF - ENDDO - ENDIF - ENDDO - ENDDO + ! + ENDIF ! ityp + ENDDO ! na + ! + ENDIF ! upf + ENDDO ! nt + ! + ENDDO ! ipert ! CALL mp_sum(int3, intra_pool_comm) ! IF (noncolin) CALL set_int3_nc( npe ) ! - DEALLOCATE (aux1) - DEALLOCATE (aux2) - DEALLOCATE (veff) - DEALLOCATE (ylmk0) - DEALLOCATE (qgm) - DEALLOCATE (qmod) - DEALLOCATE (qg) +!DMRM + !write(*,'(a,e20.12)') 'int3 = ', & + !SUM((REAL(REAL(int3(:,:,:,:,:))))**2)+SUM((REAL(AIMAG(int3(:,:,:,:,:))))**2) +!END + ! + DEALLOCATE(aux1) + DEALLOCATE(aux2) + DEALLOCATE(veff) + DEALLOCATE(ylmk0) + DEALLOCATE(qgm) + DEALLOCATE(qmod) + DEALLOCATE(qg) ! CALL stop_clock('newdq2') ! diff --git a/EPW/src/printing.f90 b/EPW/src/printing.f90 index 5546bf604..624b8f152 100644 --- a/EPW/src/printing.f90 +++ b/EPW/src/printing.f90 @@ -1520,11 +1520,11 @@ CALL print_clock ('epq_init') WRITE( stdout, * ) CALL print_clock ('epq_init') - IF (nlcc_any) call print_clock ('set_drhoc') + IF (nlcc_any) CALL print_clock ('set_drhoc') CALL print_clock ('init_vloc') CALL print_clock ('init_us_1') CALL print_clock ('newd') - CALL print_clock ('dvanqq') + CALL print_clock ('dvanqq2') CALL print_clock ('drho') WRITE( stdout, * ) ! diff --git a/test-suite/extract-epw.x b/test-suite/extract-epw.x index 716da0484..e835ed113 100755 --- a/test-suite/extract-epw.x +++ b/test-suite/extract-epw.x @@ -65,7 +65,7 @@ pi=`grep "Re[Pi]=" $fname | awk '{print $4; print $7; print $10}'` mobvb=`grep "Mobility VB Fermi level" $fname | awk '{print $5}'` mobcb=`grep "Mobility CB Fermi level" $fname | awk '{print $5}'` density=`grep " x-axis" $fname | awk '{print $1; print $2; print $3}'` -mobx=`grep " x-axis" $fname | awk '{print $4}'` +mobx=`grep " x-axis" $fname | awk '{print $5}'` mobav=`grep " avg" $fname | awk '{print $1}'` mobxZ=`grep " x-axis [Z]" $fname | awk '{print $1; print $2; print $3; print $4}'` indabs=`grep " (cm-1)" $fname | awk '{print $1; print $2; print $3; print $4}'`