Nonlinear core corrections in Gamma-only phonon code fixed

git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@5939 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
giannozz 2009-09-18 09:42:36 +00:00
parent ef4a792bd3
commit 79f2ca2262
4 changed files with 31 additions and 27 deletions

View File

@ -1,5 +1,8 @@
Fixed in version 4.1.1:
* Gamma-only phonon code wasn't working any longer if pseudopotentials
with nonlinear core correction were used
* Check of lspinorb flag consistency between left/right lead and
scattering region in pwcond.x was not working properly; wrong
print-out of E-Ef when Nchannels=0 also fixed.

View File

@ -7,15 +7,15 @@
!
!---------------------------------------------------------------------
subroutine dvb_cc (nlcc,npseu,ngm,nr1,nr2,nr3,nrx1,nrx2,nrx3, &
nl,rho_core,dmuxc,ga,aux,dvb_nlcc)
nl,igtongl,rho_core,dmuxc,ga,aux,dvb_nlcc)
!---------------------------------------------------------------------
! calculate the core-correction contribution to Delta V bare
!
implicit none
integer:: npseu,ngm,nr1,nr2,nr3,nrx1,nrx2,nrx3,np,ng,i
logical :: nlcc(npseu)
integer :: nl(ngm)
real(8) :: rho_core(ngm), dmuxc(nrx1*nrx2*nrx3)
integer :: nl(ngm), igtongl(ngm)
real(8) :: rho_core(*), dmuxc(nrx1*nrx2*nrx3)
complex(8) :: ga(ngm), dvb_nlcc(ngm), aux(nrx1*nrx2*nrx3)
!
do np=1,npseu
@ -26,7 +26,7 @@ subroutine dvb_cc (nlcc,npseu,ngm,nr1,nr2,nr3,nrx1,nrx2,nrx3, &
!
aux(:) = (0.d0, 0.d0)
do ng=1,ngm
aux(nl(ng)) = ga(ng) * rho_core(ng)
aux(nl(ng)) = ga(ng) * rho_core(igtongl(ng))
end do
call cft3(aux,nr1,nr2,nr3,nrx1,nr2,nr3,1)
!

View File

@ -20,7 +20,7 @@ subroutine dvpsi_kb(kpoint,nu)
USE uspp_param, ONLY: upf, nh, nhm
USE uspp, ONLY: dvan, nkb, vkb
USE gvect, ONLY : gstart, nl, nlm, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
nrxx, ngm, g, gg, igtongl
nrxx, ngl, ngm, g, gg, gl, igtongl
USE vlocal, ONLY: vloc
USE wvfct, ONLY: nbnd, npwx, npw, g2kin, igk
USE wavefunctions_module, ONLY: evc, psic
@ -28,16 +28,16 @@ subroutine dvpsi_kb(kpoint,nu)
!
implicit none
integer :: ibnd, ir, ih, jkb, ik, na, nu, ng, mu, nt, kpoint
complex(DP), pointer:: work(:,:), workcc(:), dvloc(:), dvb_cc(:)
complex(DP), pointer:: work(:,:), dvloc(:), dvb_cc(:)
complex(DP) :: exc
real(DP), pointer :: bec1(:,:), bec2(:,:), dv(:)
real(DP), pointer :: bec1(:,:), bec2(:,:), rhocg(:), dv(:)
real(DP) :: gu, gtau
logical :: has_nlcc
!
call start_clock('dvpsi_kb')
!
has_nlcc=.false.
workcc => psic
rhocg => auxr
dv => auxr
dvloc => aux2
dvb_cc => aux3
@ -47,9 +47,9 @@ subroutine dvpsi_kb(kpoint,nu)
mu = 3*(na-1)
if ( u(mu+1,nu)**2+u(mu+2,nu)**2+u(mu+3,nu)**2.gt. 1.0d-12) then
nt=ityp(na)
if (upf(nt)%nlcc) call drhoc (ngm, gg, omega, tpiba2, rgrid(nt)%mesh,&
rgrid(nt)%dx, rgrid(nt)%r, upf(nt)%rho_atc,&
workcc)
if (upf(nt)%nlcc) call drhoc (ngl, gl, omega, tpiba2, rgrid(nt)%mesh,&
rgrid(nt)%r, rgrid(nt)%rab, upf(nt)%rho_atc,&
rhocg )
has_nlcc = has_nlcc .or. upf(nt)%nlcc
do ng = 1,ngm
gtau = tpi * ( g(1,ng)*tau(1,na) + &
@ -60,7 +60,8 @@ subroutine dvpsi_kb(kpoint,nu)
g(3,ng)*u(mu+3,nu) )
exc = gu * CMPLX(-sin(gtau),-cos(gtau),kind=DP)
dvloc (nl(ng))=dvloc (nl(ng)) + vloc(igtongl(ng),nt)*exc
if (upf(nt)%nlcc) dvb_cc(nl(ng)) = dvb_cc(nl(ng)) + workcc(ng) * exc
if (upf(nt)%nlcc) &
dvb_cc(nl(ng)) = dvb_cc(nl(ng)) + rhocg (igtongl(ng)) * exc
end do
end if
end do

View File

@ -18,7 +18,7 @@ subroutine dynmatcc(dyncc)
USE ener, ONLY : etxc, vtxc
USE uspp_param, ONLY : upf
USE gvect, ONLY : nl, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
nrxx, ngm, g, gg
nrxx, ngm, igtongl, ngl, g, gg, gl
USE scf, ONLY : rho, rho_core, rhog_core
USE wavefunctions_module, ONLY: psic
USE wvfct, ONLY: nbnd, npwx, npw, g2kin, igk
@ -32,7 +32,7 @@ subroutine dynmatcc(dyncc)
integer:: i,j,na,nb,nta,ntb,ir,ig,nt, nu_i,nu_j,mu_i,mu_j
complex(DP), pointer:: vxc(:), work1(:), gc(:,:)
complex(DP) :: exc
real(DP), allocatable:: drhocc(:), dyncc1(:,:,:,:)
real(DP), allocatable:: rhocg(:), dyncc1(:,:,:,:)
real(DP) :: exg
logical :: nlcc(ntyp)
!
@ -47,7 +47,7 @@ subroutine dynmatcc(dyncc)
vxc => aux2
allocate ( dyncc1( 3,nat,3,nat))
allocate ( gc ( nrxx, 3))
allocate ( drhocc( nrxx))
allocate ( rhocg( ngl))
!
call v_xc (rho, rho_core, rhog_core, etxc, vtxc, vxc)
!
@ -59,17 +59,17 @@ subroutine dynmatcc(dyncc)
do na=1,nat
nta=ityp(na)
if ( upf(nta)%nlcc ) then
call drhoc (ngm, gg, omega, tpiba2, rgrid(nta)%mesh, rgrid(nta)%dx, &
rgrid(nta)%r, upf(nta)%rho_atc, drhocc)
call drhoc (ngl, gl, omega, tpiba2, rgrid(nta)%mesh, rgrid(nta)%r, &
rgrid(nta)%rab, upf(nta)%rho_atc, rhocg)
do ig=1,ngm
exg = tpi* ( g(1,ig)*tau(1,na) + &
g(2,ig)*tau(2,na) + &
g(3,ig)*tau(3,na) )
exc = CMPLX(cos(exg),-sin(exg),kind=DP)*tpiba2
work1(ig)= drhocc(ig)* exc * CONJG(vxc(nl(ig)))
gc(ig,1) = g(1,ig) * exc * CMPLX(0.0d0,-1.0d0,kind=DP)
gc(ig,2) = g(2,ig) * exc * CMPLX(0.0d0,-1.0d0,kind=DP)
gc(ig,3) = g(3,ig) * exc * CMPLX(0.0d0,-1.0d0,kind=DP)
work1(ig)= rhocg(igtongl(ig))* exc * CONJG(vxc(nl(ig)))
gc(ig,1) = g(1,ig) * exc * (0.0d0,-1.0d0)
gc(ig,2) = g(2,ig) * exc * (0.0d0,-1.0d0)
gc(ig,3) = g(3,ig) * exc * (0.0d0,-1.0d0)
end do
do i=1,3
do j=1,3
@ -81,20 +81,20 @@ subroutine dynmatcc(dyncc)
end do
do i=1,3
call dvb_cc (nlcc,nt,ngm,nr1,nr2,nr3,nrx1,nrx2,nrx3, &
nl,drhocc,dmuxc,gc(1,i),aux3,gc(1,i))
nl,igtongl,rhocg,dmuxc,gc(1,i),aux3,gc(1,i))
end do
do nb=1,nat
ntb=ityp(nb)
if ( upf(ntb)%nlcc ) then
call drhoc (ngm, gg, omega, tpiba2, rgrid(ntb)%mesh, &
rgrid(ntb)%dx, rgrid(ntb)%r, upf(ntb)%rho_atc,&
drhocc)
call drhoc (ngl, gl, omega, tpiba2, rgrid(ntb)%mesh, &
rgrid(ntb)%r, rgrid(ntb)%rab, upf(ntb)%rho_atc,&
rhocg)
do ig=1,ngm
exg = tpi* ( g(1,ig)*tau(1,nb) + &
g(2,ig)*tau(2,nb) + &
g(3,ig)*tau(3,nb) )
exc = -CMPLX(sin(exg),cos(exg),kind=DP)
work1(ig) = exc * drhocc(ig)
work1(ig) = exc * rhocg(igtongl(ig))
end do
do i=1,3
do j=1,3
@ -109,8 +109,8 @@ subroutine dynmatcc(dyncc)
end if
end do
!
deallocate(rhocg)
deallocate(gc)
deallocate(drhocc)
#ifdef __PARA
call mp_sum( dyncc1, intra_pool_comm )
#endif