Commit f3114d67 authored by Fabrizio Ferrari's avatar Fabrizio Ferrari

merge develop with rho_updw_mz

parent c0088709
......@@ -24,7 +24,7 @@
USE io_files, ONLY : tmp_dir
USE klist, ONLY : xk, nks, nkstot
USE lsda_mod, ONLY : nspin, starting_magnetization
USE scf, ONLY : v, vrs, vltot, rho, kedtau
USE scf, ONLY : v, vrs, vltot, rho, kedtau, rhoz_or_updw
USE gvect, ONLY : ngm
USE symm_base, ONLY : nsym, s, irt, t_rev, time_reversal, invs, sr, &
inverse_s
......@@ -125,8 +125,14 @@
!
! 3.1) Setup all gradient correction stuff
!
!^
IF ( nspin == 2 ) CALL rhoz_or_updw( rho, 'r_and_g', 'rhoz_updw' )
!
CALL setup_dgc
!
IF ( nspin == 2 ) CALL rhoz_or_updw( rho, 'r_and_g', 'updw_rhoz' )
!^
!
! 4) Computes the inverse of each matrix of the crystal symmetry group
!
CALL inverse_s()
......
......@@ -33,7 +33,7 @@ SUBROUTINE energies_xc( lda, n, m, psi, e_xc, e_h,ispin )
USE lsda_mod, ONLY : current_spin
USE gvect, ONLY : gstart
USE io_global, ONLY :stdout
USE scf, ONLY : rho, vltot, vrs, rho_core,rhog_core, scf_type
USE scf, ONLY : rho, vltot, vrs, rho_core,rhog_core, scf_type, rhoz_or_updw
USE constants, ONLY :rytoev
USE io_files, ONLY : diropn
USE mp, ONLY : mp_sum, mp_barrier
......@@ -111,12 +111,17 @@ SUBROUTINE energies_xc( lda, n, m, psi, e_xc, e_h,ispin )
! ... the local potential V_Loc psi. First the psi in real space
!set exchange and correlation potential
if(.not.allocated(psic)) write(stdout,*) 'psic not allocated'
!^
IF (nspin == 2) CALL rhoz_or_updw( rho, 'r_and_g', 'rhoz_updw' )
!
if (dft_is_meta()) then
! call v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v%of_r, v%kin_r )
else
CALL v_xc( rho, rho_core, rhog_core, etxc, vtxc, vr )
endif
!
IF (nspin == 2) CALL rhoz_or_updw( rho, 'r_and_g', 'updw_rhoz' )
!^
do is=1,nspin
vrs(:,is)=vr(:,is)
......@@ -357,13 +362,17 @@ SUBROUTINE energies_xc( lda, n, m, psi, e_xc, e_h,ispin )
if(.not.allocated(vr)) write(stdout,*) 'vr not allocated'
allocate(rho_fake_core(dfftp%nnr))
rho_fake_core(:)=0.d0
!^
IF (nspin == 2) CALL rhoz_or_updw( rho, 'r_and_g', 'rhoz_updw' )
!
if (dft_is_meta()) then
! call v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v%of_r, v%kin_r )
else
CALL v_xc( rho, rho_core, rhog_core, etxc, vtxc, vr )
endif
!
IF (nspin == 2) CALL rhoz_or_updw( rho, 'r_and_g', 'updw_rhoz' )
!^
deallocate(rho_fake_core)
......
......@@ -25,13 +25,16 @@ SUBROUTINE hp_ns_trace
IMPLICIT NONE
INTEGER :: is, na, na1, na2, nt, nt1, nt2, m1, m2, ldim, viz
REAL(DP), ALLOCATABLE :: nsaux(:,:) ! auxiliary array for occupations
REAL(DP) :: sfact
!
ALLOCATE(ns(nat))
ALLOCATE(nsaux(nat,nspin))
ns(:) = 0.0d0
nsaux(:,:) = 0.0d0
!
sfact=2.d0
IF (nspin==2) THEN
sfact=1.d0
ALLOCATE(magn(nat))
magn(:) = 0.0d0
ENDIF
......@@ -53,12 +56,8 @@ SUBROUTINE hp_ns_trace
ENDDO
ENDDO
!
IF (nspin==1) THEN
ns(na) = 2.0d0 * nsaux(na,1)
ELSE
ns(na) = nsaux(na,1) + nsaux(na,2)
magn(na) = nsaux(na,1) - nsaux(na,2)
ENDIF
ns(na) = nsaux(na,1) * sfact
IF (nspin /= 1) magn(na) = nsaux(na,2)
!
ENDIF
!
......
......@@ -51,7 +51,7 @@ SUBROUTINE hp_setup_q()
USE cell_base, ONLY : at, bg
USE io_global, ONLY : stdout
USE lsda_mod, ONLY : nspin
USE scf, ONLY : v, vrs, vltot, rho, kedtau
USE scf, ONLY : v, vrs, vltot, rho, kedtau, rhoz_or_updw
USE fft_base, ONLY : dfftp
USE gvect, ONLY : ngm
USE gvecs, ONLY : doublegrid
......@@ -95,8 +95,14 @@ SUBROUTINE hp_setup_q()
!
! 5) Setup gradient correction stuff
!
!^
IF (nspin == 2) CALL rhoz_or_updw( rho, 'r_and_g', 'rhoz_updw' )
!
CALL setup_dgc()
!
IF (nspin == 2) CALL rhoz_or_updw( rho, 'r_and_g', 'updw_rhoz' )
!^
!
! 6) Compute the inverse of each matrix of the crystal symmetry group
!
CALL inverse_s()
......
......@@ -52,7 +52,7 @@ SUBROUTINE commutator_Vhubx_psi(ik, ipol)
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
REAL(DP) :: nsaux
REAL(DP) :: nsaux, sgn
REAL(DP), ALLOCATABLE :: xyz(:,:), gk(:,:), g2k(:)
COMPLEX(DP), ALLOCATABLE :: dkwfcbessel(:,:), dkwfcylmr(:,:), dkwfcatomk(:,:), &
dpqq26(:,:), dpqq38(:,:), dpqq47(:,:), dkvkbbessel(:,:), &
......@@ -262,7 +262,8 @@ SUBROUTINE commutator_Vhubx_psi(ik, ipol)
!
trm = (0.d0, 0.d0)
!
nsaux = rho%ns(m1, m2, current_spin, nah)
sgn= REAL( 2*MOD(current_spin,2)-1 )
nsaux = ( rho%ns(m1, m2, 1, nah) + sgn*rho%ns(m1, m2, nspin, nah) )*0.5d0
!
DO ibnd = 1, nbnd_occ(ik)
trm(:,ibnd) = aux_1234(:) * proj1(ibnd,ihubst2) + & ! term_1234
......@@ -322,7 +323,8 @@ SUBROUTINE commutator_Vhubx_psi(ik, ipol)
!
trm = (0.d0, 0.d0)
!
nsaux = rho%ns(m1, m2, op_spin, nah)
sgn = REAL( 2*MOD(op_spin,2)-1 )
nsaux = ( rho%ns(m1, m2, 1, nah) + sgn*rho%ns(m1, m2, nspin, nah) )*0.5d0
!
DO ibnd = 1, nbnd_occ(ik)
trm(:,ibnd) = aux_1234(:) * proj1(ibnd,ihubst2) + & ! term_1234
......
......@@ -26,7 +26,7 @@ subroutine dv_of_drho (dvscf, add_nlcc, drhoc)
USE cell_base, ONLY : tpiba2, omega
USE noncollin_module, ONLY : nspin_lsda, nspin_mag, nspin_gga
USE funct, ONLY : dft_is_gradient, dft_is_nonlocc
USE scf, ONLY : rho, rho_core
USE scf, ONLY : rho, rho_core, rhoz_or_updw
USE uspp, ONLY : nlcc_any
USE control_flags, ONLY : gamma_only
USE martyna_tuckerman, ONLY : wg_corr_h, do_comp_mt
......@@ -78,8 +78,8 @@ subroutine dv_of_drho (dvscf, add_nlcc, drhoc)
fac = 1.d0 / DBLE (nspin_lsda)
!
if (nlcc_any.and.add_nlcc) then
rho%of_r(:, 1) = rho%of_r(:, 1) + rho_core (:)
do is = 1, nspin_lsda
rho%of_r(:, is) = rho%of_r(:, is) + fac * rho_core (:)
dvscf(:, is) = dvscf(:, is) + fac * drhoc (:)
enddo
endif
......@@ -96,17 +96,31 @@ subroutine dv_of_drho (dvscf, add_nlcc, drhoc)
! NB: If nlcc=.true. we need to add here its contribution.
! grho contains already the core charge
!
if ( dft_is_gradient() ) call dgradcorr &
(dfftp, rho%of_r, grho, dvxc_rr, dvxc_sr, dvxc_ss, dvxc_s, xq, &
dvscf, nspin_mag, nspin_gga, g, dvaux)
if ( dft_is_gradient() ) then
!^
IF (nspin_lsda == 2) CALL rhoz_or_updw( rho, 'only_r', 'rhoz_updw' )
!
call dgradcorr &
(dfftp, rho%of_r, grho, dvxc_rr, dvxc_sr, dvxc_ss, dvxc_s, xq, &
dvscf, nspin_mag, nspin_gga, g, dvaux)
!
IF (nspin_lsda == 2) CALL rhoz_or_updw( rho, 'only_r', 'updw_rhoz' )
!^
endif
!
if (dft_is_nonlocc()) then
!^
IF (nspin_lsda == 2) CALL rhoz_or_updw( rho, 'only_r', 'rhoz_updw' )
!
call dnonloccorr(rho%of_r, dvscf, xq, dvaux)
!
IF (nspin_lsda == 2) CALL rhoz_or_updw( rho, 'only_r', 'updw_rhoz' )
!^
endif
!
if (nlcc_any.and.add_nlcc) then
rho%of_r(:, 1) = rho%of_r(:, 1) - rho_core (:)
do is = 1, nspin_lsda
rho%of_r(:, is) = rho%of_r(:, is) - fac * rho_core (:)
dvscf(:, is) = dvscf(:, is) - fac * drhoc (:)
enddo
endif
......
......@@ -35,8 +35,8 @@ SUBROUTINE setup_dmuxc
!
IF (lsda) THEN
DO ir = 1, dfftp%nnr
rhoup = rho%of_r (ir, 1) + 0.5d0 * rho_core (ir)
rhodw = rho%of_r (ir, 2) + 0.5d0 * rho_core (ir)
rhoup = ( rho%of_r(ir, 1) + rho%of_r(ir, 2) + rho_core(ir) ) * 0.5d0
rhodw = ( rho%of_r(ir, 1) - rho%of_r(ir, 2) + rho_core(ir) ) * 0.5d0
CALL dmxc_spin (rhoup, rhodw, dmuxc(ir,1,1), dmuxc(ir,2,1), &
dmuxc(ir,1,2), dmuxc(ir,2,2) )
ENDDO
......
......@@ -11,7 +11,7 @@
! ... adapted to accept any kind of density by Oliviero Andreussi
!
!--------------------------------------------------------------------
SUBROUTINE compute_dipole( nnr, nspin, rho, r0, dipole, quadrupole )
SUBROUTINE compute_dipole( nnr, rho, r0, dipole, quadrupole )
!--------------------------------------------------------------------
USE kinds, ONLY : DP
USE cell_base, ONLY : at, bg, alat, omega
......@@ -25,8 +25,8 @@ SUBROUTINE compute_dipole( nnr, nspin, rho, r0, dipole, quadrupole )
!
! nnr is passed in input, but nnr should match dfftp%nnr
! for the calculation to be meaningful
INTEGER, INTENT(IN) :: nnr, nspin
REAL(DP), INTENT(IN) :: rho( nnr, nspin )
INTEGER, INTENT(IN) :: nnr
REAL(DP), INTENT(IN) :: rho( nnr )
REAL(DP), INTENT(IN) :: r0(3)
REAL(DP), INTENT(OUT) :: dipole(0:3), quadrupole(3)
!
......@@ -85,9 +85,7 @@ SUBROUTINE compute_dipole( nnr, nspin, rho, r0, dipole, quadrupole )
!
CALL cryst_to_cart( 1, r, at, 1 )
!
rhoir = rho( ir, 1 )
!
IF ( nspin == 2 ) rhoir = rhoir + rho(ir,2)
rhoir = rho( ir )
!
! ... dipole(0) = charge density
!
......
......@@ -561,22 +561,9 @@ MODULE io_base
!
DO ns = 1, nspin
!
! Workaround for LSDA, while waiting for much-needed harmonization:
! we have rhoup and rhodw, we write rhotot=up+dw and rhodif=up-dw
!
IF ( ns == 1 .AND. nspin == 2 ) THEN
DO ig = 1, ngm
rhoaux(ig) = rho(ig,ns) + rho(ig,ns+1)
END DO
ELSE IF ( ns == 2 .AND. nspin == 2 ) THEN
DO ig = 1, ngm
rhoaux(ig) = rho(ig,ns-1) - rho(ig,ns)
END DO
ELSE
DO ig = 1, ngm
DO ig = 1, ngm
rhoaux(ig) = rho(ig,ns)
END DO
END IF
END DO
!
rho_g = 0
CALL mergewf( rhoaux, rho_g, ngm, ig_l2g, me_in_group, &
......@@ -647,7 +634,6 @@ MODULE io_base
!
COMPLEX(dp), ALLOCATABLE :: rho_g(:)
COMPLEX(dp), ALLOCATABLE :: rhoaux(:)
COMPLEX(dp) :: rhoup, rhodw
REAL(dp) :: b1(3), b2(3), b3(3)
INTEGER :: ngm, nspin_, isup, isdw
INTEGER :: iun, ns, ig, ierr
......@@ -789,19 +775,8 @@ MODULE io_base
DO ig = 1, ngm
rho(ig,ns) = rhoaux(ig)
END DO
!
! Workaround for LSDA, while waiting for much-needed harmonization:
! if file contains rhotot=up+dw and rhodif=up-dw (nspin_=2), and
! if we want rhoup and rho down (nspin=2), convert
!
IF ( nspin_ == 2 .AND. nspin == 2 .AND. ns == 2 ) THEN
DO ig = 1, ngm
rhoup = (rho(ig,ns-1) + rhoaux(ig)) / 2.0_dp
rhodw = (rho(ig,ns-1) - rhoaux(ig)) / 2.0_dp
rho(ig,ns-1)= rhoup
rho(ig,ns )= rhodw
END DO
ELSE IF ( nspin_ == 2 .AND. nspin == 4 .AND. ns == 2) THEN
IF ( nspin_ == 2 .AND. nspin == 4 .AND. ns == 2) THEN
DO ig = 1, ngm
rho(ig, 4 ) = rho(ig, 2)
rho(ig, 2 ) = cmplx(0.d0,0.d0, KIND = DP)
......
......@@ -15,7 +15,7 @@ SUBROUTINE A_h(npw,e,h,ah)
USE lsda_mod, ONLY : current_spin, nspin
USE wvfct, ONLY: nbnd, npwx, g2kin
USE wavefunctions, ONLY: evc, psic
USE scf, ONLY : vrs, rho
USE scf, ONLY : vrs, rho, rhoz_or_updw
USE fft_base, ONLY : dffts, dfftp
USE fft_interfaces, ONLY : fwfft, invfft
USE gvect, ONLY : gstart, g, gg
......@@ -120,9 +120,17 @@ SUBROUTINE A_h(npw,e,h,ah)
! add gradient correction contribution (if any)
!
CALL start_clock('dgradcorr')
IF (dft_is_gradient() ) CALL dgradcor1 &
(dfftp, rho%of_r, grho, dvxc_rr, dvxc_sr, dvxc_ss, dvxc_s, &
drho, dpsic, nspin, g, dv)
IF (dft_is_gradient() ) THEN
!^
IF (nspin == 2) CALL rhoz_or_updw(rho, 'only_r', 'rhoz_updw')
!
CALL dgradcor1 &
(dfftp, rho%of_r, grho, dvxc_rr, dvxc_sr, dvxc_ss, dvxc_s, &
drho, dpsic, nspin, g, dv)
!
IF (nspin == 2) CALL rhoz_or_updw(rho, 'only_r', 'updw_rhoz')
!^
ENDIF
CALL stop_clock('dgradcorr')
NULLIFY(dpsic)
NULLIFY (drho)
......
......@@ -13,7 +13,7 @@ SUBROUTINE cg_setup
USE kinds, ONLY: DP
USE cell_base, ONLY: bg
USE ions_base, ONLY: nat, ntyp => nsp, ityp, tau, amass
USE scf, ONLY: rho, rho_core, v, vltot, vrs, kedtau
USE scf, ONLY: rho, rho_core, v, vltot, vrs, kedtau, rhoz_or_updw
USE uspp, ONLY: vkb, nlcc_any
USE uspp_param, ONLY: upf
USE mp_global, ONLY: kunit
......@@ -36,7 +36,7 @@ SUBROUTINE cg_setup
INTEGER :: i, l, nt, ik
LOGICAL :: exst
CHARACTER (len=256) :: filint
REAL(DP) :: rhotot
REAL(DP) :: rhotot, sgn
INTEGER :: ndr, kunittmp, ierr
REAL(DP) :: edum(1,1), wdum(1,1)
!
......@@ -86,17 +86,22 @@ SUBROUTINE cg_setup
!
! derivative of the xc potential
!
sgn = DBLE(2*MOD(current_spin,2)-1)
dmuxc(:) = 0.d0
DO i = 1,dfftp%nnr
rhotot = rho%of_r(i,current_spin)+rho_core(i)
rhotot = ( rho%of_r(i,1) + sgn*rho%of_r(i,nspin) )*0.5d0 + rho_core(i)
IF ( rhotot> 1.d-30 ) dmuxc(i)= dmxc( rhotot)
IF ( rhotot<-1.d-30 ) dmuxc(i)=-dmxc(-rhotot)
ENDDO
!
! initialize data needed for gradient corrections
!^
IF (nspin == 2) CALL rhoz_or_updw(rho, 'r_and_g', 'rhoz_updw')
!
CALL cg_setupdgc
!
IF (nspin == 2) CALL rhoz_or_updw(rho, 'r_and_g', 'updw_rhoz')
!^
iunres=88
!
! open the wavefunction file (already existing)
......
......@@ -15,12 +15,13 @@ SUBROUTINE dynmatcc(dyncc)
USE atom, ONLY : rgrid
USE constants, ONLY : tpi
USE cell_base, ONLY : omega, tpiba2
USE lsda_mod, ONLY : nspin
USE ener, ONLY : etxc, vtxc
USE uspp_param, ONLY : upf
USE fft_base, ONLY : dfftp
USE fft_interfaces, ONLY : fwfft
USE gvect, ONLY : ngm, igtongl, ngl, g, gg, gl
USE scf, ONLY : rho, rho_core, rhog_core
USE scf, ONLY : rho, rho_core, rhog_core, rhoz_or_updw
USE wavefunctions, ONLY: psic
USE cgcom
USE mp_global, ONLY : intra_pool_comm
......@@ -48,9 +49,13 @@ SUBROUTINE dynmatcc(dyncc)
ALLOCATE ( dyncc1( 3,nat,3,nat))
ALLOCATE ( gc ( dfftp%nnr, 3))
ALLOCATE ( rhocg( ngl))
!^
IF (nspin == 2) CALL rhoz_or_updw(rho, 'r_and_g', 'rhoz_updw')
!
CALL v_xc (rho, rho_core, rhog_core, etxc, vtxc, vxc)
!
IF (nspin == 2) CALL rhoz_or_updw(rho, 'r_and_g', 'updw_rhoz')
!^
CALL fwfft ( 'Rho', vxc, dfftp )
!
dyncc1(:,:,:,:) = 0.d0
......
......@@ -418,15 +418,15 @@ SUBROUTINE cg_neweps
USE cell_base, ONLY : omega
USE ions_base, ONLY : nat, tau
USE fft_base, ONLY : dfftp
USE scf, ONLY : rho, rho_core
USE lsda_mod, ONLY : current_spin
USE scf, ONLY : rho, rho_core, rhoz_or_updw
USE lsda_mod, ONLY : nspin, current_spin
USE funct, ONLY : dmxc
USE cgcom
!
IMPLICIT NONE
INTEGER :: i, j
REAL(DP) :: rhotot, chi(3,3)
REAL(DP) :: rhotot, sgn, chi(3,3)
!
! recalculate self-consistent potential etc
!
......@@ -434,17 +434,22 @@ SUBROUTINE cg_neweps
!
! new derivative of the xc potential
!
sgn = DBLE(2*MOD(current_spin,2)-1)
dmuxc(:) = 0.d0
DO i = 1,dfftp%nnr
rhotot = rho%of_r(i,current_spin)+rho_core(i)
rhotot = ( rho%of_r(i,1) + sgn*rho%of_r(i,nspin) )*0.5d0 + rho_core(i)
IF ( rhotot> 1.d-30 ) dmuxc(i)= dmxc( rhotot)
IF ( rhotot<-1.d-30 ) dmuxc(i)=-dmxc(-rhotot)
ENDDO
!
! re-initialize data needed for gradient corrections
!^
IF (nspin == 2) CALL rhoz_or_updw(rho, 'r_and_g', 'rhoz_updw')
!
CALL cg_setupdgc
!
IF (nspin == 2) CALL rhoz_or_updw(rho, 'r_and_g', 'updw_rhoz')
!^
! calculate linear response to macroscopic fields
!
CALL macro
......
......@@ -16,7 +16,7 @@ SUBROUTINE rhod2vkb(dyn0)
USE constants, ONLY: tpi
USE ions_base, ONLY : nat, tau, ityp, ntyp => nsp
USE cell_base, ONLY : tpiba2, tpiba, omega
USE lsda_mod, ONLY : current_spin
USE lsda_mod, ONLY : nspin, current_spin
USE gvect, ONLY : ngm, g, igtongl
USE gvecw, ONLY: gcutw
USE wvfct, ONLY: nbnd, npwx
......@@ -38,7 +38,7 @@ SUBROUTINE rhod2vkb(dyn0)
!
INTEGER :: npw, i, ih, ibnd, na, nt, nu_i,nu_j,mu_i,mu_j, ir, ng, jkb, ik, &
ipol, jpol, ijpol
real(DP) :: weight, fac, gtau
real(DP) :: weight, fac, gtau, sgn
real(DP), ALLOCATABLE :: dynloc(:,:), dynkb(:,:)
COMPLEX(DP), ALLOCATABLE :: dvkb(:,:)
real (DP), ALLOCATABLE :: becp(:,:), becp1(:,:,:), becp2(:,:,:)
......@@ -49,8 +49,9 @@ SUBROUTINE rhod2vkb(dyn0)
!
ALLOCATE ( dynloc( 3*nat, nmodes))
dynloc (:,:) = 0.d0
sgn = DBLE(2*MOD(current_spin,2)-1)
DO ir = 1,dfftp%nnr
psic(ir) = rho%of_r(ir,current_spin)
psic(ir) = ( rho%of_r(ir,1) + sgn*rho%of_r(ir,nspin) )*0.5d0
ENDDO
CALL fwfft ('Rho', psic, dfftp)
DO nu_i = 1,nmodes
......
......@@ -15,7 +15,7 @@ subroutine addnlcc (imode0, drhoscf, npe)
USE ions_base, ONLY : nat
use funct, only : dft_is_gradient, dft_is_nonlocc
USE cell_base, ONLY : omega
use scf, only : rho, rho_core
use scf, only : rho, rho_core, rhoz_or_updw
USE gvect, ONLY : g, ngm
USE fft_base, ONLY : dfftp
USE noncollin_module, ONLY : nspin_lsda, nspin_gga, nspin_mag
......@@ -73,9 +73,7 @@ subroutine addnlcc (imode0, drhoscf, npe)
!
! add core charge to the density
!
DO is=1,nspin_lsda
rho%of_r(:,is) = rho%of_r(:,is) + fac * rho_core(:)
ENDDO
rho%of_r(:,1) = rho%of_r(:,1) + rho_core(:)
!
! Compute the change of xc potential due to the perturbation
!
......@@ -98,11 +96,26 @@ subroutine addnlcc (imode0, drhoscf, npe)
! add gradient correction to xc, NB: if nlcc is true we need to add here
! its contribution. grho contains already the core charge
!
if ( dft_is_gradient() ) &
call dgradcorr (dfftp, rho%of_r, grho, dvxc_rr, dvxc_sr, dvxc_ss, &
if ( dft_is_gradient() ) then
!^
if (nspin_lsda == 2) call rhoz_or_updw(rho, 'only_r', 'rhoz_updw')
!
call dgradcorr (dfftp, rho%of_r, grho, dvxc_rr, dvxc_sr, dvxc_ss, &
dvxc_s, xq, drhoscf (1, 1, ipert), nspin_mag, nspin_gga, g, dvaux)
if (dft_is_nonlocc()) &
call dnonloccorr(rho%of_r, drhoscf (1, 1, ipert), xq, dvaux)
!
if (nspin_lsda == 2) call rhoz_or_updw(rho, 'only_r', 'updw_rhoz')
!^
endif
!
if (dft_is_nonlocc()) then
!^
if (nspin_lsda == 2) call rhoz_or_updw(rho, 'only_r', 'rhoz_updw')
!
call dnonloccorr(rho%of_r, drhoscf (1, 1, ipert), xq, dvaux)
!
if (nspin_lsda == 2) call rhoz_or_updw(rho, 'only_r', 'updw_rhoz')
!^
endif
do is = 1, nspin_lsda
call daxpy (2 * dfftp%nnr, - fac, drhoc, 1, drhoscf (1, is, ipert), 1)
......@@ -120,9 +133,8 @@ subroutine addnlcc (imode0, drhoscf, npe)
enddo
enddo
enddo
DO is=1,nspin_lsda
rho%of_r(:,is) = rho%of_r(:,is) - fac * rho_core(:)
ENDDO
!
rho%of_r(:,1) = rho%of_r(:,1) - rho_core(:)
!
! collect contributions from all r/G points.
!
......
......@@ -12,7 +12,7 @@ SUBROUTINE addnlcc_zstar_eu_us( drhoscf )
USE kinds, ONLY : DP
USE funct, only : dft_is_gradient, dft_is_nonlocc
USE scf, only : rho, rho_core
USE scf, only : rho, rho_core, rhoz_or_updw
USE cell_base, ONLY : omega
USE gvect, ONLY : ngm, g
USE fft_base, ONLY : dfftp
......@@ -59,9 +59,7 @@ SUBROUTINE addnlcc_zstar_eu_us( drhoscf )
dvaux = (0.0_dp,0.0_dp)
CALL addcore (mode, drhoc)
DO is = 1, nspin_lsda
rho%of_r(:,is) = rho%of_r(:,is) + fac * rho_core
END DO
rho%of_r(:,1) = rho%of_r(:,1) + rho_core
DO is = 1, nspin_mag
DO is1 = 1, nspin_mag
......@@ -77,16 +75,28 @@ SUBROUTINE addnlcc_zstar_eu_us( drhoscf )
! its contribution. grho contains already the core charge
!
IF ( dft_is_gradient() ) &
CALL dgradcorr (dfftp, rho%of_r, grho, &
dvxc_rr, dvxc_sr, dvxc_ss, dvxc_s, xq, drhoscf (1,1,ipol),&
nspin_mag, nspin_gga, g, dvaux)
if (dft_is_nonlocc()) &
call dnonloccorr(rho%of_r, drhoscf (1, 1, ipol), xq, dvaux)
DO is = 1, nspin_lsda
rho%of_r(:,is) = rho%of_r(:,is) - fac * rho_core
END DO
IF ( dft_is_gradient() ) THEN
!^
IF (nspin_lsda == 2) CALL rhoz_or_updw(rho, 'only_r', 'rhoz_updw')
!
CALL dgradcorr (dfftp, rho%of_r, grho, dvxc_rr, dvxc_sr, dvxc_ss, &
dvxc_s, xq, drhoscf(1,1,ipol), nspin_mag, nspin_gga, g, dvaux)
!
IF (nspin_lsda == 2) call rhoz_or_updw(rho, 'only_r', 'updw_rhoz')
!^
ENDIF
!
IF (dft_is_nonlocc()) THEN
!^
IF (nspin_lsda == 2) CALL rhoz_or_updw(rho, 'only_r', 'rhoz_updw')
!
CALL dnonloccorr(rho%of_r, drhoscf (1, 1, ipol), xq, dvaux)
!
IF (nspin_lsda == 2) CALL rhoz_or_updw(rho, 'only_r', 'updw_rhoz')
!^
ENDIF
!
rho%of_r(:,1) = rho%of_r(:,1) - rho_core
DO is = 1, nspin_lsda
zstareu0(ipol,mode) = zstareu0(ipol,mode) - &
......
......@@ -47,7 +47,7 @@ SUBROUTINE dvqhub_barepsi_us (ik, uact)
USE control_lr, ONLY : lgamma, ofsbeta
USE units_lr, ONLY : iuatwfc, iuatswfc
USE uspp_param, ONLY : nh
USE lsda_mod, ONLY : lsda, current_spin, isk
USE lsda_mod, ONLY : nspin, lsda, current_spin, isk
USE wavefunctions, ONLY : evc
USE eqv, ONLY : dvpsi
USE scf, ONLY : rho
......@@ -67,9 +67,11 @@ SUBROUTINE dvqhub_barepsi_us (ik, uact)
INTEGER :: i, j, k, icart, counter, na, nt, l, ih, n, mu, ig, &
ihubst, ihubst1, ihubst2, nah, m, m1, m2, ibnd, op_spin, &
ikk, ikq, npw, npwq, ibeta
COMPLEX(DP) :: rhons_m1m2
COMPLEX(DP), ALLOCATABLE :: aux1(:), aux2(:), aux3(:), aux4(:), aux5(:), &
dqsphi(:,:), dmqsphi(:,:), dvqi(:,:), dvqhbar(:,:,:,:), &
vkb_(:,:), dwfcatom_(:)
REAL(DP) :: sgn_cs, sgn_op
COMPLEX(DP), EXTERNAL :: ZDOTC
!
ALLOCATE (proj1(nbnd,nwfcU))
......@@ -250,11 +252,14 @@ SUBROUTINE dvqhub_barepsi_us (ik, uact)
!
ihubst2 = offsetU(nah) + m2
!
rhons_m1m2 = ( rho%ns(m1,m2,1,nah) + sgn_cs * &
rho%ns(m1,m2,nspin,nah) )*0.5d0
!
DO ig = 1, npwq
!
aux2(ig) = dqsphi(ig,ihubst1) * rho%ns(m1,m2,current_spin,nah) &
aux2(ig) = dqsphi(ig,ihubst1) * rhons_m1m2 &
* proj1(ibnd, ihubst2)
aux4(ig) = swfcatomkpq(ig,ihubst1) * rho%ns(m1,m2,current_spin,nah) &
aux4(ig) = swfcatomkpq(ig,ihubst1) * rhons_m1m2 &
* proj2(ibnd, ihubst2)