diff --git a/PH/adddvepsi_us.f90 b/PH/adddvepsi_us.f90 index c0714b099..3ebd65962 100644 --- a/PH/adddvepsi_us.f90 +++ b/PH/adddvepsi_us.f90 @@ -8,7 +8,7 @@ subroutine adddvepsi_us(becp2,ipol,kpoint) ! This subdoutine adds to dvpsi the terms which depend on the augmentation ! charge. It assumes that the variable dpqq, has been set and it is in - ! the crystal base. + ! the crystal basis. ! It calculates the last two terms of Eq.10 in JCP 21, 9934 (2004). ! P^+_c is applied in solve_e. ! @@ -23,16 +23,17 @@ subroutine adddvepsi_us(becp2,ipol,kpoint) USE noncollin_module, ONLY : noncolin, npol USE uspp_param, only: nh USE phus, ONLY : becp1, dpqq, dpqq_so + USE becmod, ONLY : bec_type USE control_ph, ONLY: nbnd_occ USE eqv, ONLY : dvpsi implicit none integer, intent(in) :: ipol, kpoint - complex(DP), intent(in) :: becp2(nkb,npol,nbnd) + TYPE(bec_type), intent(in) :: becp2 complex(DP), allocatable :: ps(:), ps_nc(:,:) - integer:: ijkb0, nt, na, ih, jh, ikb, jkb, ibnd, ip, is + integer:: ijkb0, nt, na, ih, jh, ikb, jkb, ibnd, ip, is, js, ijs IF (noncolin) THEN allocate (ps_nc(nbnd,npol)) @@ -55,28 +56,30 @@ subroutine adddvepsi_us(becp2,ipol,kpoint) jkb = ijkb0 + jh do ibnd=1, nbnd_occ(kpoint) IF (noncolin) THEN - DO ip=1,npol - IF (lspinorb) THEN - ps_nc(ibnd,ip)=ps_nc(ibnd,ip) + & - (0.d0,1.d0)*(becp2(jkb,1,ibnd)* & - qq_so(ih,jh,1+(ip-1)*2,nt) + & - becp2(jkb,2,ibnd) * & - qq_so(ih,jh,2+(ip-1)*2,nt) ) & - + becp1(kpoint)%nc(jkb,1,ibnd)* & - dpqq_so(ih,jh,1+(ip-1)*2,ipol,nt) & - + becp1(kpoint)%nc(jkb,2,ibnd)* & - dpqq_so(ih,jh,2+(ip-1)*2,ipol,nt) - ELSE - ps_nc(ibnd,ip)=ps_nc(ibnd,ip)+ & - becp2(jkb,ip,ibnd)*(0.d0,1.d0)* & - qq(ih,jh,nt)+becp1(kpoint)%nc(jkb,ip,ibnd) & - * dpqq(ih,jh,ipol,nt) - END IF - END DO + IF (lspinorb) THEN + ijs=0 + do is=1,npol + do js=1,npol + ijs=ijs+1 + ps_nc(ibnd,is)=ps_nc(ibnd,is) + & + qq_so(ih,jh,ijs,nt)* & + (0.d0,1.d0)*becp2%nc(jkb,js,ibnd) & + + becp1(kpoint)%nc(jkb,js,ibnd)* & + dpqq_so(ih,jh,ijs,ipol,nt) + enddo + enddo + ELSE + DO is=1,npol + ps_nc(ibnd,is)=ps_nc(ibnd,is)+ & + qq(ih,jh,nt)*becp2%nc(jkb,is,ibnd)*(0.d0,1.d0) & + + dpqq(ih,jh,ipol,nt)* & + becp1(kpoint)%nc(jkb,is,ibnd) + END DO + END IF ELSE - ps(ibnd) = ps(ibnd) & - + becp2(jkb,1,ibnd)*(0.d0,1.d0)*qq(ih,jh,nt)+ & - becp1(kpoint)%k(jkb,ibnd)*dpqq(ih,jh,ipol,nt) + ps(ibnd) = ps(ibnd)+qq(ih,jh,nt)*becp2%k(jkb,ibnd) & + *(0.d0,1.d0) + & + dpqq(ih,jh,ipol,nt)* becp1(kpoint)%k(jkb,ibnd) END IF enddo enddo diff --git a/PH/adddvscf.f90 b/PH/adddvscf.f90 index 97036b561..67c6ef118 100644 --- a/PH/adddvscf.f90 +++ b/PH/adddvscf.f90 @@ -38,7 +38,7 @@ subroutine adddvscf (ipert, ik) ! ! And the local variables ! - integer :: na, nt, ibnd, ih, jh, ijkb0, ikk, ikb, jkb, ip + integer :: na, nt, ibnd, ih, jh, ijkb0, ikk, ikb, jkb, is, js, ijs ! counter on atoms ! counter on atomic types ! counter on bands @@ -74,19 +74,18 @@ subroutine adddvscf (ipert, ik) do jh = 1, nh (nt) jkb = ijkb0 + jh IF (noncolin) THEN - sum_nc(1)=sum_nc(1)+ & - int3_nc(ih,jh,ipert,na,1)* & - becp1(ik)%nc (jkb, 1, ibnd)+ & - int3_nc(ih,jh,ipert,na,2)* & - becp1(ik)%nc (jkb, 2, ibnd) - sum_nc(2)=sum_nc(2)+ & - int3_nc(ih,jh,ipert,na,3)* & - becp1(ik)%nc (jkb, 1, ibnd)+ & - int3_nc(ih,jh,ipert,na,4)* & - becp1(ik)%nc (jkb, 2, ibnd) + ijs=0 + do is=1,npol + do js=1,npol + ijs=ijs+1 + sum_nc(is)=sum_nc(is)+ & + int3_nc(ih,jh,ipert,na,ijs)* & + becp1(ik)%nc(jkb, js, ibnd) + enddo + enddo ELSE sum = sum + int3 (ih, jh, ipert, na, current_spin)*& - becp1(ik)%k (jkb, ibnd) + becp1(ik)%k(jkb, ibnd) END IF enddo IF (noncolin) THEN diff --git a/PH/ch_psi_all.f90 b/PH/ch_psi_all.f90 index b2afb4f4f..a010d8f47 100644 --- a/PH/ch_psi_all.f90 +++ b/PH/ch_psi_all.f90 @@ -111,11 +111,7 @@ subroutine ch_psi_all (n, h, ah, e, ik, m) ! ! And apply S again ! - IF (noncolin) THEN - call calbec (n, vkb, hpsi, becp, m) - ELSE - call calbec (n, vkb, hpsi, becp, m) - END IF + call calbec (n, vkb, hpsi, becp, m) call s_psi (npwx, n, m, hpsi, spsi) do ibnd = 1, m do ig = 1, n diff --git a/PH/dvpsi_e.f90 b/PH/dvpsi_e.f90 index 4a11ff202..1282ba1df 100644 --- a/PH/dvpsi_e.f90 +++ b/PH/dvpsi_e.f90 @@ -131,19 +131,11 @@ subroutine dvpsi_e (ik, ipol) call davcio (dvpsi, lrcom, iucom, nrec, 1) ! allocate (spsi ( npwx*npol, nbnd)) - IF (noncolin) THEN - CALL calbec (npw, vkb, dvpsi, becp ) - ELSE - CALL calbec (npw, vkb, dvpsi, becp ) - END IF + CALL calbec (npw, vkb, dvpsi, becp ) CALL s_psi(npwx,npw,nbnd,dvpsi,spsi) call dcopy(2*npwx*npol*nbnd,spsi,1,dvpsi,1) deallocate (spsi) - IF (noncolin) THEN - call adddvepsi_us(becp2%nc,ipol,ik) - ELSE - call adddvepsi_us(becp2%k,ipol,ik) - END IF + call adddvepsi_us(becp2,ipol,ik) endif IF (nkb > 0) call deallocate_bec_type (becp2)