! ! Copyright (C) 2001-2008 Quantum ESPRESSO group ! 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 . ! ! !---------------------------------------------------------------------- subroutine dvqpsi_us_only (ik, uact) !---------------------------------------------------------------------- ! ! This routine calculates dV_bare/dtau * psi for one perturbation ! with a given q. The displacements are described by a vector uact. ! The result is stored in dvpsi. The routine is called for each k point ! and for each pattern u. It computes simultaneously all the bands. ! This routine implements Eq. B29 of PRB 64, 235118 (2001). ! Only the contribution of the nonlocal potential is calculated here. ! ! USE kinds, only : DP USE cell_base, ONLY : tpiba USE gvect, ONLY : g USE klist, ONLY : xk USE ions_base, ONLY : nat, ityp, ntyp => nsp USE lsda_mod, ONLY : lsda, current_spin, isk, nspin USE spin_orb, ONLY : lspinorb USE wvfct, ONLY : nbnd, npwx, et USE noncollin_module, ONLY : noncolin, npol USE uspp, ONLY: okvan, nkb, vkb USE uspp_param, ONLY: nh, nhm USE qpoint, ONLY : igkq, npwq, ikks, ikqs USE gwus, ONLY : int1, int1_nc, int2, int2_so, alphap, becp1 USE eqv, ONLY : dvpsi USE control_gw, ONLY : lgamma implicit none ! ! The dummy variables ! integer :: ik ! input: the k point complex(DP) :: uact (3 * nat) ! input: the pattern of displacements ! ! And the local variables ! integer :: na, nb, mu, nu, ikk, ikq, ig, igg, nt, ibnd, ijkb0, & ikb, jkb, ih, jh, ipol, is, js, ijs ! counter on atoms ! counter on modes ! the point k ! the point k+q ! counter on G vectors ! auxiliary counter on G vectors ! counter on atomic types ! counter on bands ! auxiliary variable for counting ! counter on becp functions ! counter on becp functions ! counter on n index ! counter on m index ! counter on polarizations real(DP), parameter :: eps = 1.d-12 complex(DP), allocatable :: ps1 (:,:), ps2 (:,:,:), aux (:), deff_nc(:,:,:,:) real(DP), allocatable :: deff(:,:,:) complex(DP), allocatable :: ps1_nc (:,:,:), ps2_nc (:,:,:,:) ! work space logical :: ok call start_clock ('dvqpsi_us_on') if (noncolin) then allocate (ps1_nc(nkb , npol, nbnd)) allocate (ps2_nc(nkb , npol, nbnd , 3)) allocate (deff_nc(nhm, nhm, nat, nspin)) else allocate (ps1 ( nkb , nbnd)) allocate (ps2 ( nkb , nbnd , 3)) allocate (deff(nhm, nhm, nat)) end if allocate (aux ( npwx)) ikk = ikks(ik) ikq = ikqs(ik) if (lsda) current_spin = isk (ikk) ! ! we first compute the coefficients of the vectors ! if (noncolin) then ps1_nc(:,:,:) = (0.d0, 0.d0) ps2_nc(:,:,:,:) = (0.d0, 0.d0) else ps1(:,:) = (0.d0, 0.d0) ps2(:,:,:) = (0.d0, 0.d0) end if do ibnd = 1, nbnd IF (noncolin) THEN CALL compute_deff_nc(deff_nc,et(ibnd,ikk)) ELSE ! This routine computes the effective value of the D-eS coefficients ! which appear often in many expressions in the US or PAW case. ! This routine is for the collinear case. CALL compute_deff(deff,et(ibnd,ikk)) ENDIF ijkb0 = 0 do nt = 1, ntyp do na = 1, nat 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) ) > eps) 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) + & deff_nc(ih,jh,na,ijs) * & alphap(ipol, ik)%nc(jkb,js,ibnd)* & uact(mu + ipol) ps2_nc(ikb,is,ibnd,ipol)= & ps2_nc(ikb,is,ibnd,ipol)+ & deff_nc(ih,jh,na,ijs) * & becp1(ik)%nc(jkb,js,ibnd) * & (0.d0,-1.d0) * uact(mu+ipol) * tpiba END DO END DO ELSE ps1 (ikb, ibnd) = ps1 (ikb, ibnd) + & deff(ih, jh, na) * & alphap(ipol, ik)%k(jkb, ibnd) * uact (mu + ipol) ps2 (ikb, ibnd, ipol) = ps2 (ikb, ibnd, ipol) +& deff(ih,jh,na)*becp1(ik)%k (jkb, ibnd) * & (0.0_DP,-1.0_DP) * 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) END DO END DO ELSE ps1 (ikb, ibnd) = ps1 (ikb, ibnd) + & (int1 (ih, jh, ipol,na, current_spin) * & becp1(ik)%k (jkb, ibnd) ) * uact (mu +ipol) END IF END IF END IF ! 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) END DO END DO 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) END DO END IF ELSE ps1 (ikb, ibnd) = ps1 (ikb, ibnd) + & (int2 (ih, jh, ipol, nb, na) * & becp1(ik)%k (jkb, ibnd) ) * uact (nu + ipol) END IF enddo endif ! okvan enddo ! ipol enddo ! jh enddo ! ih ijkb0 = ijkb0 + nh (nt) endif enddo ! na enddo ! nt enddo ! nbnd ! ! This term is proportional to beta(k+q+G) ! if (nkb.gt.0) then if (noncolin) then call zgemm ('N', 'N', npwq, nbnd*npol, nkb, & (1.d0, 0.d0), vkb, npwx, ps1_nc, nkb, (1.d0, 0.d0) , dvpsi, npwx) else call zgemm ('N', 'N', npwq, nbnd, nkb, & (1.d0, 0.d0) , vkb, npwx, ps1, nkb, (1.d0, 0.d0) , dvpsi, npwx) end if end if ! ! This term is proportional to (k+q+G)_\alpha*beta(k+q+G) ! do ikb = 1, nkb do ipol = 1, 3 ok = .false. IF (noncolin) THEN do ibnd = 1, nbnd ok = ok.or.(abs (ps2_nc (ikb, 1, ibnd, ipol) ).gt.eps).or. & (abs (ps2_nc (ikb, 2, ibnd, ipol) ).gt.eps) end do ELSE do ibnd = 1, nbnd ok = ok.or. (abs (ps2 (ikb, ibnd, ipol) ) .gt.eps) enddo ENDIF if (ok) then do ig = 1, npwq igg = igkq (ig) aux (ig) = vkb(ig, ikb) * (xk(ipol, ikq) + g(ipol, igg) ) enddo do ibnd = 1, nbnd IF (noncolin) THEN call zaxpy(npwq,ps2_nc(ikb,1,ibnd,ipol),aux,1,dvpsi(1,ibnd),1) call zaxpy(npwq,ps2_nc(ikb,2,ibnd,ipol),aux,1, & dvpsi(1+npwx,ibnd),1) ELSE call zaxpy (npwq, ps2(ikb,ibnd,ipol), aux, 1, dvpsi(1,ibnd), 1) END IF enddo endif enddo enddo deallocate (aux) IF (noncolin) THEN deallocate (ps2_nc) deallocate (ps1_nc) deallocate (deff_nc) ELSE deallocate (ps2) deallocate (ps1) deallocate (deff) END IF call stop_clock ('dvqpsi_us_on') return end subroutine dvqpsi_us_only