scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

full lsda rho conversion

This commit is contained in:
fabrizio22 2019-01-21 12:02:37 +01:00
parent 74f7a0c09a
commit 9c87cf7106
48 changed files with 483 additions and 617 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
CPV/src/chargedensity.f90
View File

@ -186,6 +186,14 @@
TRIM(tmp_dir), TRIM(prefix), ndr,postfix
CALL read_rhog ( filename, root_bgrp, intra_bgrp_comm, &
ig_l2g, nspin, rhog )
!
!^^ ... TEMPORARY FIX (newlsda) ...
IF ( nspin==2 ) THEN
rhog(:,1) = ( rhog(:,1) + rhog(:,2) )*0.5d0
rhog(:,2) = rhog(:,1) - rhog(:,2)
ENDIF
!^^.......................
!
CALL rho_g2r ( dfftp, rhog, rhor )
rhopr = rhor
first = .FALSE.

13
CPV/src/cp_restart_new.f90
View File

@ -460,10 +460,21 @@ MODULE cp_restart_new
CALL rho_r2g (dfftp,rho, rhog)
filename = TRIM(dirname) // 'charge-density'
! Only the first band group collects and writes
IF ( my_bgrp_id == root_bgrp_id ) CALL write_rhog &
IF ( my_bgrp_id == root_bgrp_id ) THEN
!
!^^ ... TEMPORARY FIX (newlsda) ...
IF ( lsda ) THEN
rhog(:,1) = rhog(:,1) + rhog(:,2)
rhog(:,2) = rhog(:,1) - rhog(:,2)*2._dp
ENDIF
!^^.......................
!
CALL write_rhog &
( filename, root_bgrp, intra_bgrp_comm, &
tpiba*b1, tpiba*b2, tpiba*b3, gamma_only, &
mill, ig_l2g, rhog, ecutrho )
ENDIF
!
DEALLOCATE ( rhog )
END IF

39
CPV/src/vofrho.f90
View File

@ -117,7 +117,23 @@ SUBROUTINE vofrho_x( nfi, rhor, drhor, rhog, drhog, rhos, rhoc, tfirst, &
CALL start_clock( 'ts_vdw' )
ALLOCATE (stmp(3,nat))
stmp(:,:) = tau0(:,ind_bck(:))
CALL tsvdw_calculate(stmp,rhor)
!
!^^ ... TEMPORARY FIX (newlsda) ...
IF ( nspin==2 ) THEN
rhor(:,1) = rhor(:,1) + rhor(:,2)
rhor(:,2) = rhor(:,1) - rhor(:,2)*2._dp
ENDIF
!^^.......................
!
CALL tsvdw_calculate(stmp,rhor(:,1), nspin)
!
!^^ ... TEMPORARY FIX ...
IF ( nspin==2 ) THEN
rhor(:,1) = ( rhor(:,1) + rhor(:,2) )*0.5_dp
rhor(:,2) = rhor(:,1) - rhor(:,2)
ENDIF
!^^.......................
!
DEALLOCATE (stmp)
CALL stop_clock( 'ts_vdw' )
!
@ -385,14 +401,29 @@ SUBROUTINE vofrho_x( nfi, rhor, drhor, rhog, drhog, rhos, rhoc, tfirst, &
END DO
denlc(:,:) = 0.0_dp
inlc = get_inlc()
!
!^^ ... TEMPORARY FIX (newlsda) ...
IF ( nspin==2 ) THEN
rhosave(:,1) = rhosave(:,1) + rhosave(:,2)
rhosave(:,2) = rhosave(:,1) - rhosave(:,2)*2._dp
ENDIF
!^^.......................
!
if (inlc==1 .or. inlc==2) then
if (nspin>2) call errore('stres_vdW_DF', 'vdW+DF non implemented in spin polarized calculations',1)
CALL stress_vdW_DF(rhosave, rhocsave, nspin, denlc )
CALL stress_vdW_DF(rhosave(:,1), rhocsave, nspin, denlc )
elseif (inlc == 3) then
if (nspin>2) call errore('stress_rVV10', 'rVV10 non implemented with nspin>2',1)
CALL stress_rVV10(rhosave, rhocsave, nspin, denlc )
CALL stress_rVV10(rhosave(:,1), rhocsave, nspin, denlc )
end if
!
!^^ ... TEMPORARY FIX ...
IF ( nspin==2 ) THEN
rhosave(:,1) = ( rhosave(:,1) + rhosave(:,2) )*0.5_dp
rhosave(:,2) = rhosave(:,1) - rhosave(:,2)
ENDIF
!^^.......................
!
dxc(:,:) = dxc(:,:) - omega/e2 * MATMUL(denlc,TRANSPOSE(ainv))
END IF
DEALLOCATE ( rhocsave, rhosave )

6
EPW/src/epw_setup.f90
View File

@ -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, rhoz_or_updw
USE scf, ONLY : v, vrs, vltot, rho, kedtau
USE gvect, ONLY : ngm
USE symm_base, ONLY : nsym, s, irt, t_rev, time_reversal, invs, sr, &
inverse_s
@ -125,12 +125,8 @@
!
! 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()

12
GWW/pw4gww/energies_xc.f90
View File

@ -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, rhoz_or_updw
USE scf, ONLY : rho, vltot, vrs, rho_core,rhog_core, scf_type
USE constants, ONLY :rytoev
USE io_files, ONLY : diropn
USE mp, ONLY : mp_sum, mp_barrier
@ -111,8 +111,6 @@ 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 )
@ -120,9 +118,6 @@ SUBROUTINE energies_xc( lda, n, m, psi, e_xc, e_h,ispin )
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)
if(doublegrid) call fft_interpolate(dfftp, vrs(:,is),dffts,vrs(:,is)) ! interpolate from dense to smooth
@ -362,8 +357,6 @@ 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 )
@ -371,9 +364,6 @@ SUBROUTINE energies_xc( lda, n, m, psi, e_xc, e_h,ispin )
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)

8
HP/src/hp_setup_q.f90
View File

@ -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, rhoz_or_updw
USE scf, ONLY : v, vrs, vltot, rho, kedtau
USE fft_base, ONLY : dfftp
USE gvect, ONLY : ngm
USE gvecs, ONLY : doublegrid
@ -95,14 +95,8 @@ 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()

8
LR_Modules/commutator_Vhubx_psi.f90
View File

@ -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, sgn
REAL(DP) :: nsaux
REAL(DP), ALLOCATABLE :: xyz(:,:), gk(:,:), g2k(:)
COMPLEX(DP), ALLOCATABLE :: dkwfcbessel(:,:), dkwfcylmr(:,:), dkwfcatomk(:,:), &
dpqq26(:,:), dpqq38(:,:), dpqq47(:,:), dkvkbbessel(:,:), &
@ -262,8 +262,7 @@ SUBROUTINE commutator_Vhubx_psi(ik, ipol)
!
trm = (0.d0, 0.d0)
!
sgn= REAL( 2*MOD(current_spin,2)-1 )
nsaux = ( rho%ns(m1, m2, 1, nah) + sgn*rho%ns(m1, m2, nspin, nah) )*0.5d0
nsaux = rho%ns(m1, m2, current_spin, nah)
!
DO ibnd = 1, nbnd_occ(ik)
trm(:,ibnd) = aux_1234(:) * proj1(ibnd,ihubst2) + & ! term_1234
@ -323,8 +322,7 @@ SUBROUTINE commutator_Vhubx_psi(ik, ipol)
!
trm = (0.d0, 0.d0)
!
sgn = REAL( 2*MOD(op_spin,2)-1 )
nsaux = ( rho%ns(m1, m2, 1, nah) + sgn*rho%ns(m1, m2, nspin, nah) )*0.5d0
nsaux = rho%ns(m1, m2, op_spin, nah)
!
DO ibnd = 1, nbnd_occ(ik)
trm(:,ibnd) = aux_1234(:) * proj1(ibnd,ihubst2) + & ! term_1234

9
LR_Modules/dgradcorr.f90
View File

@ -35,7 +35,7 @@ subroutine dgradcorr (dfft, rho, grho, dvxc_rr, dvxc_sr, dvxc_ss, &
COMPLEX(DP), INTENT(INOUT) :: dvxc (dfft%nnr, nspin)
real(DP), parameter :: epsr = 1.0d-6, epsg = 1.0d-10
real(DP) :: grho2, seg, seg0, amag
real(DP) :: grho2, seg, seg0, amag, sgn(2)
complex(DP) :: s1, fact, term
complex(DP) :: a (2, 2, 2), b (2, 2, 2, 2), c (2, 2, 2), &
ps (2, 2), ps1 (3, 2, 2), ps2 (3, 2, 2, 2)
@ -57,6 +57,8 @@ subroutine dgradcorr (dfft, rho, grho, dvxc_rr, dvxc_sr, dvxc_ss, &
allocate (gdrho( 3, dfft%nnr, nspin0))
allocate (h( 3, dfft%nnr, nspin0))
allocate (dh( dfft%nnr))
sgn(1)=1.d0 ; sgn(2)=-1.d0
h (:, :, :) = (0.d0, 0.d0)
if (noncolin.and.domag) then
@ -98,7 +100,10 @@ subroutine dgradcorr (dfft, rho, grho, dvxc_rr, dvxc_sr, dvxc_ss, &
ELSE
DO is = 1, nspin0
CALL fft_qgradient (dfft, drho(1,is), xq, g, gdrho (1, 1, is) )
rhoout(:,is)=rho(:,is)
!
! rhoout, if LSDA, is in (up,down) format
!
rhoout(:,is)=( rho(:,1) + sgn(is)*rho(:,nspin0) )*0.5_dp
drhoout(:,is)=drho(:,is)
ENDDO
ENDIF

26
LR_Modules/dv_of_drho.f90
View File

@ -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, rhoz_or_updw
USE scf, ONLY : rho, rho_core
USE uspp, ONLY : nlcc_any
USE control_flags, ONLY : gamma_only
USE martyna_tuckerman, ONLY : wg_corr_h, do_comp_mt
@ -96,27 +96,11 @@ 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() ) 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_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_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 ( dft_is_nonlocc() ) call dnonloccorr(rho%of_r, dvscf, xq, dvaux)
!
if (nlcc_any.and.add_nlcc) then
rho%of_r(:, 1) = rho%of_r(:, 1) - rho_core (:)

24
LR_Modules/setup_dgc.f90
View File

@ -20,9 +20,9 @@ subroutine setup_dgc
USE fft_interfaces, ONLY : fwfft
USE gvect, ONLY : ngm, g
USE spin_orb, ONLY : domag
USE scf, ONLY : rho, rho_core, rhog_core
USE scf, ONLY : rho, rho_core, rhog_core, rhoz_or_updw
USE noncollin_module, ONLY : noncolin, ux, nspin_gga, nspin_mag
USE wavefunctions, ONLY : psic
USE wavefunctions, ONLY : psic
USE kinds, ONLY : DP
USE funct, ONLY : dft_is_gradient, gcxc, gcx_spin, &
gcc_spin, dgcxc, dgcxc_spin
@ -36,7 +36,7 @@ subroutine setup_dgc
v1x, v2x, v1c, v2c, vrrx, vsrx, vssx, vrrc, vsrc, vssc, v1xup, &
v1xdw, v2xup, v2xdw, v1cup, v1cdw, vrrxup, vrrxdw, vrsxup, vrsxdw, &
vssxup, vssxdw, vrrcup, vrrcdw, vrscup, vrscdw, vrzcup, vrzcdw, &
amag, seg, seg0
amag, seg, seg0, sgn(2)
COMPLEX(DP), ALLOCATABLE :: rhogout(:,:)
real(DP), allocatable :: rhoout(:,:)
real (DP), parameter :: epsr = 1.0d-6, epsg = 1.0d-10
@ -64,6 +64,7 @@ subroutine setup_dgc
dvxc_ss(:,:,:) = 0.d0
dvxc_s (:,:,:) = 0.d0
grho (:,:,:) = 0.d0
sgn(1)=1.d0 ; sgn(2)=-1.d0
!
! add rho_core
!
@ -86,21 +87,29 @@ subroutine setup_dgc
END DO
DEALLOCATE(rhogout)
ELSE
!
! for convenience, if LSDA, rhoout is kept in (up,down) format
!
do is = 1, nspin_gga
rhoout(:,is) = rho%of_r(:,is)
rhoout(:,is) = ( rho%of_r(:,1) + sgn(is)*rho%of_r(:,nspin_gga) )*0.5d0
enddo
!
! if LSDA rho%of_g is temporarily converted in (up,down) format
!
call rhoz_or_updw(rho, 'only_g', '->updw')
!
if (nlcc_any) then
do is = 1, nspin_gga
rhoout(:,is) = fac * rho_core(:) + rho%of_r(:,is)
rhoout(:,is) = fac * rho_core(:) + rhoout(:,is)
rho%of_g(:,is) = fac * rhog_core(:) + rho%of_g(:,is)
enddo
endif
!
do is = 1, nspin_gga
call fft_gradient_g2r (dfftp, rho%of_g (1, is), g, grho (1, 1, is) )
enddo
END IF
do k = 1, dfftp%nnr
grho2 (1) = grho (1, k, 1) **2 + grho (2, k, 1) **2 + grho (3, k, 1) **2
if (nspin_gga == 1) then
@ -167,6 +176,9 @@ subroutine setup_dgc
rho%of_g(:,is) = rho%of_g(:,is) - fac * rhog_core(:)
enddo
endif
!
CALL rhoz_or_updw(rho, 'only_g', '->rhoz')
!
endif
DEALLOCATE(rhoout)

4
Modules/funct.f90
View File

@ -2897,9 +2897,9 @@ subroutine nlc (rho_valence, rho_core, nspin, enl, vnl, v)
elseif (inlc == 3) then
if(imeta == 0) then
call xc_rVV10 (rho_valence, rho_core, nspin, enl, vnl, v)
call xc_rVV10 (rho_valence(:,1), rho_core, nspin, enl, vnl, v)
else
call xc_rVV10 (rho_valence, rho_core, nspin, enl, vnl, v, 15.7_dp)
call xc_rVV10 (rho_valence(:,1), rho_core, nspin, enl, vnl, v, 15.7_dp)
endif
else
enl = 0.0_DP

51
Modules/tsvdw.f90
View File

@ -591,7 +591,7 @@ PRIVATE :: GetVdWParam
!--------------------------------------------------------------------------------------------------------------
!
!--------------------------------------------------------------------------------------------------------------
SUBROUTINE tsvdw_calculate(tauin, rhor)
SUBROUTINE tsvdw_calculate(tauin, rhor, nspin)
!--------------------------------------------------------------------------------------------------------------
! TS-vdW Management Code: Manages entire calculation of TS-vdW energy, wavefunction forces, and ion forces via
! calls to PRIVATE subroutines below (called in each MD step). The calls to tsvdw_initialize and tsvdw_finalize
@ -602,7 +602,8 @@ PRIVATE :: GetVdWParam
!
! I/O variables
!
REAL(DP), INTENT(IN) :: rhor(:,:)
REAL(DP), INTENT(IN) :: rhor(:)
INTEGER, INTENT(IN) :: nspin
REAL(DP) :: tauin(3,nat)
!
! Parallel initialization...
@ -619,7 +620,7 @@ PRIVATE :: GetVdWParam
!
! Obtain molecular charge density given on the real-space mesh...
!
CALL tsvdw_rhotot( rhor )
CALL tsvdw_rhotot( rhor, nspin )
!
! Determine spherical atomic integration domains and atom overlap (bit array)...
! Compute molecular pro-density (superposition of atomic densities) on the real-space mesh...
@ -969,36 +970,30 @@ PRIVATE :: GetVdWParam
!--------------------------------------------------------------------------------------------------------------
!
!--------------------------------------------------------------------------------------------------------------
SUBROUTINE tsvdw_rhotot( rhor )
SUBROUTINE tsvdw_rhotot( rhor, nspin )
!--------------------------------------------------------------------------------------------------------------
!
IMPLICIT NONE
!
REAL(DP), INTENT(IN) :: rhor(:,:)
REAL(DP), INTENT(IN) :: rhor(:)
INTEGER, INTENT(IN) :: nspin
!
! Local variables
!
INTEGER :: ir,ierr,nspin
REAL(DP), DIMENSION(:), ALLOCATABLE :: rhor_tmp1,rhor_tmp2
INTEGER :: ir,ierr
REAL(DP), DIMENSION(:), ALLOCATABLE :: rhor_tmp
!
CALL start_clock('tsvdw_rhotot')
!
! Initialization of rhotot array (local copy of the real-space charge density)...
!
ALLOCATE(rhotot(nr1*nr2*nr3)); rhotot=0.0_DP
nspin = SIZE(rhor,2)
IF ( nspin < 1 .OR. nspin > 2 ) CALL errore ('tsvdw','invalid nspin',1)
#if defined(__MPI)
!
! Initialization of rhor_tmp temporary buffers...
!
ALLOCATE(rhor_tmp1(nr1*nr2*nr3)); rhor_tmp1=0.0_DP
!
IF (nspin.EQ.2) THEN
!
ALLOCATE(rhor_tmp2(nr1*nr2*nr3)); rhor_tmp2=0.0_DP
!
END IF
ALLOCATE(rhor_tmp(nr1*nr2*nr3)); rhor_tmp=0.0_DP
!
! Collect distributed rhor and broadcast to all processors...
!
@ -1016,36 +1011,20 @@ PRIVATE :: GetVdWParam
!
END DO
!
CALL MPI_ALLGATHERV(rhor(1,1),dfftp%nr3p(me_bgrp+1)*nr1*nr2,&
MPI_DOUBLE_PRECISION,rhor_tmp1(1),recvcount,rdispls,&
CALL MPI_ALLGATHERV(rhor(1),dfftp%nr3p(me_bgrp+1)*nr1*nr2,&
MPI_DOUBLE_PRECISION,rhor_tmp(1),recvcount,rdispls,&
MPI_DOUBLE_PRECISION,intra_bgrp_comm,ierr)
!
IF (nspin.EQ.2) THEN
!
CALL MPI_ALLGATHERV(rhor(1,2),dfftp%nr3p(me_bgrp+1)*nr1*nr2,&
MPI_DOUBLE_PRECISION,rhor_tmp2(1),recvcount,rdispls,&
MPI_DOUBLE_PRECISION,intra_bgrp_comm,ierr)
!
END IF
!
! Transfer rhor temporary arrays to rhotot array...
!
rhotot=rhor_tmp1
!
IF (nspin.EQ.2) THEN
!
rhotot=rhotot+rhor_tmp2
!
END IF
rhotot=rhor_tmp
!
! Clean-up temporary arrays...
!
IF (ALLOCATED(rhor_tmp1)) DEALLOCATE(rhor_tmp1)
IF (ALLOCATED(rhor_tmp2)) DEALLOCATE(rhor_tmp2)
IF (ALLOCATED(rhor_tmp)) DEALLOCATE(rhor_tmp)
!
#else
rhotot(:) = rhor(:,1)
IF (nspin == 2) rhotot(:) = rhotot(:) + rhor(:,2)
rhotot(:) = rhor(:)
#endif
CALL stop_clock('tsvdw_rhotot')
!

27
Modules/xc_rVV10.f90
View File

@ -53,7 +53,7 @@ CONTAINS
!! Local variables
!! ----------------------------------------------------------------------------------
! _
real(dp), intent(IN) :: rho_valence(:,:) !
real(dp), intent(IN) :: rho_valence(:) !
real(dp), intent(IN) :: rho_core(:) ! PWSCF input variables
INTEGER, INTENT(IN) :: nspin !
real(dp), intent(inout) :: etxc, vtxc, v(:,:) !_
@ -129,12 +129,8 @@ CONTAINS
!! ---------------------------------------------------------------------------------------
!! Add together the valence and core charge densities to get the total charge density
!total_rho = rho_valence(:,1) + rho_core(:)
if (nspin == 2) then
total_rho = rho_valence(:,1) + rho_valence(:,2) + rho_core(:)
else
total_rho = rho_valence(:,1) + rho_core(:)
endif
!
total_rho = rho_valence(:) + rho_core(:)
!! -------------------------------------------------------------------------
!! Here we calculate the gradient in reciprocal space using FFT
@ -207,13 +203,8 @@ CONTAINS
grid_cell_volume = omega/(dfftp%nr1*dfftp%nr2*dfftp%nr3)
do i_grid = 1, dfftp%nnr
vtxc = vtxc + grid_cell_volume*rho_valence(i_grid,1)*potential(i_grid)
vtxc = vtxc + grid_cell_volume*rho_valence(i_grid)*potential(i_grid)
end do
if (nspin==2) then
do i_grid = 1, dfftp%nnr
vtxc = vtxc + grid_cell_volume*rho_valence(i_grid,2)*potential(i_grid)
end do
endif
deallocate(potential)
@ -240,7 +231,7 @@ CONTAINS
implicit none
real(dp), intent(IN) :: rho_valence(:,:) !
real(dp), intent(IN) :: rho_valence(:) !
real(dp), intent(IN) :: rho_core(:) ! Input variables
INTEGER, INTENT(IN) :: nspin
real(dp), intent(inout) :: sigma(3,3) !
@ -282,13 +273,7 @@ CONTAINS
!! ---------------------------------------------------------------------------------------
!! Charge
!! ---------------------------------------------------------------------------------------
!total_rho = rho_valence(:,1) + rho_core(:)
if (nspin == 2) then
total_rho = rho_valence(:,1) + rho_valence(:,2) + rho_core(:)
else
total_rho = rho_valence(:,1) + rho_core(:)
endif
total_rho = rho_valence(:) + rho_core(:)
!! -------------------------------------------------------------------------
!! Here we calculate the gradient in reciprocal space using FFT

22
Modules/xc_vdW_DF.f90
View File

@ -477,8 +477,8 @@ CONTAINS
! it is only non-zero for pseudopotentials with non-local core
! corrections.
rho_up = rho_valence(:,1) + 0.5D0*rho_core(:)
rho_down = rho_valence(:,2) + 0.5D0*rho_core(:)
rho_up = ( rho_valence(:,1) + rho_valence(:,2) + rho_core(:) )*0.5D0
rho_down = ( rho_valence(:,1) - rho_valence(:,2) + rho_core(:) )*0.5D0
total_rho = rho_up + rho_down
#if defined (__SPIN_BALANCED)
@ -550,10 +550,10 @@ CONTAINS
grid_cell_volume = omega/(dfftp%nr1*dfftp%nr2*dfftp%nr3)
do i_grid = 1, dfftp%nnr
vtxc = vtxc + e2 * grid_cell_volume * rho_valence(i_grid,1) * potential_up (i_grid) &
+ e2 * grid_cell_volume * rho_valence(i_grid,2) * potential_down(i_grid)
end do
vtxc = vtxc + e2 * grid_cell_volume * (rho_valence(i_grid,1) + rho_valence(i_grid,2)) &
* potential_up (i_grid) &
+ e2 * grid_cell_volume * (rho_valence(i_grid,1) - rho_valence(i_grid,2)) &
* potential_down(i_grid)
deallocate( potential_up, potential_down, q0, grad_rho, grad_rho_up, &
grad_rho_down, dq0_drho_up, dq0_drho_down, thetas, &
@ -1543,7 +1543,7 @@ CONTAINS
implicit none
real(dp), intent(IN) :: rho_valence(:,:) !
real(dp), intent(IN) :: rho_valence(:) !
real(dp), intent(IN) :: rho_core(:) ! Input variables
integer, intent(IN) :: nspin !
real(dp), intent(inout) :: sigma(3,3) !
@ -1597,11 +1597,11 @@ CONTAINS
! --------------------------------------------------------------------
! Charge
total_rho = rho_valence(:,1) + rho_core(:)
total_rho = rho_valence(:) + rho_core(:)
#if defined (__SPIN_BALANCED)
if ( nspin == 2 ) then
total_rho = rho_valence(:,2) + total_rho(:)
end if
!if ( nspin == 2 ) then
! total_rho = rho_valence(:,2) + total_rho(:)
!end if
#endif

48
PHonon/Gamma/cg_setupdgc.f90
View File

@ -13,7 +13,7 @@ SUBROUTINE cg_setupdgc
!
USE kinds, ONLY: dp
USE constants, ONLY: e2
USE scf, ONLY: rho, rho_core, rhog_core, rhoz_or_updw
USE scf, ONLY: rho, rho_core, rhog_core
USE funct, ONLY: gcxc, gcx_spin, gcc_spin, dgcxc, dgcxc_spin, dft_is_gradient
USE fft_base, ONLY: dfftp
USE gvect, ONLY: ngm, g
@ -24,7 +24,7 @@ SUBROUTINE cg_setupdgc
IMPLICIT NONE
INTEGER k, is
real(DP) &
& grho2(2), rh, zeta, grh2, epsr, epsg, &
& grho2(2), rh, zeta, grh2, epsr, epsg, fac(2), sgn(2), &
& sx,sc,v1x,v2x,v1c,v2c,vrrx,vsrx,vssx, &
& vrrc,vsrc,vssc, &
& v1xup,v1xdw,v2xup,v2xdw, &
@ -44,20 +44,30 @@ SUBROUTINE cg_setupdgc
dvxc_s (:,:,:) = 0.d0
grho (:,:,:) = 0.d0
!
! add rho_core
fac(1) = dble(mod(nspin,2)+1) ; fac(2) = 0.d0
sgn(1) = 1.d0 ; sgn(2) = -1.d0
!
! adds rho_core and, if LSDA, directly converts rho in (up,down) format
!
IF (nlcc_any) THEN
rho%of_r(:,1) = rho%of_r(:,1) + rho_core(:)
rho%of_g(:,1) = rho%of_g(:,1) +rhog_core(:)
DO is=1,nspin
rho%of_r(:,is) = ( fac(is)* rho_core(:) + rho%of_r(:,1) + &
sgn(is)*rho%of_r(:,nspin) )*is/2.d0
rho%of_g(:,is) = ( fac(is)*rhog_core(:) + rho%of_g(:,1) + &
sgn(is)*rho%of_g(:,nspin) )*is/2.d0
ENDDO
ENDIF
!
DO is=1,nspin
CALL fft_gradient_g2r (dfftp, rho%of_g(1,is), g, grho(1,1,is))
ENDDO
!
IF (nspin==1) THEN
CALL fft_gradient_g2r (dfftp, rho%of_g(1,1), g, grho(1,1,1))
DO k = 1,dfftp%nnr
grho2(1)=grho(1,k,1)**2+grho(2,k,1)**2+grho(3,k,1)**2
IF (abs(rho%of_r(k,1))>epsr.and.grho2(1)>epsg) THEN
CALL gcxc(rho%of_r(k,1),grho2(1),sx,sc,v1x,v2x,v1c,v2c)
CALL dgcxc(rho%of_r(k,1),grho2(1),vrrx,vsrx,vssx,vrrc,vsrc,vssc)
CALL gcxc(rho%of_r(k,nspin),grho2(1),sx,sc,v1x,v2x,v1c,v2c)
CALL dgcxc(rho%of_r(k,nspin),grho2(1),vrrx,vsrx,vssx,vrrc,vsrc,vssc)
dvxc_rr(k,1,1) = e2 * ( vrrx + vrrc )
dvxc_sr(k,1,1) = e2 * ( vsrx + vsrc )
dvxc_ss(k,1,1) = e2 * ( vssx + vssc )
@ -65,15 +75,9 @@ SUBROUTINE cg_setupdgc
ENDIF
ENDDO
ELSE
CALL rhoz_or_updw(rho, 'r_and_g', 'rhoz_updw')
!! bring (up+down,uo-down) charge to (up,down)
DO is=1,nspin
CALL fft_gradient_g2r (dfftp, rho%of_g(1,is), g, grho(1,1,is))
ENDDO
DO k = 1,dfftp%nnr
rh=rho%of_r(k,1)+rho%of_r(k,2)
grho2(1)=grho(1,k,1)**2+grho(2,k,1)**2+grho(3,k,1)**2
grho2(2)=grho(1,k,2)**2+grho(2,k,2)**2+grho(3,k,2)**2
rh=rho%of_r(k,1)+rho%of_r(k,2)
grh2= (grho(1,k,1)+grho(1,k,2))**2 &
+ (grho(2,k,1)+grho(2,k,2))**2 &
+ (grho(3,k,1)+grho(3,k,2))**2
@ -81,7 +85,7 @@ SUBROUTINE cg_setupdgc
CALL gcx_spin(rho%of_r(k,1),rho%of_r(k,2),grho2(1),grho2(2),sx, &
v1xup,v1xdw,v2xup,v2xdw)
!
CALL dgcxc_spin(rho%of_r(k,1),rho%of_r(k,2),grho(1,k,1),grho(1,k,2), &
CALL dgcxc_spin(rho%of_r(k,1),rho%of_r(k,2),grho(1,k,1),grho(1,k,2), &
vrrxup,vrrxdw,vrsxup,vrsxdw,vssxup,vssxdw, &
vrrcup,vrrcdw,vrscup,vrscdw,vssc,vrzcup,vrzcdw)
!
@ -119,11 +123,17 @@ SUBROUTINE cg_setupdgc
dvxc_ss(k,2,1)=e2*vssc
dvxc_ss(k,2,2)=e2*(vssxdw+vssc)
ENDDO
CALL rhoz_or_updw(rho, 'r_and_g', 'updw_rhoz')
ENDIF
IF (nlcc_any) THEN
rho%of_r(:,1) = rho%of_r(:,1) - rho_core(:)
rho%of_g(:,1) = rho%of_g(:,1) -rhog_core(:)
!
! restore rho to the original values
!
DO is=1,nspin
rho%of_r(:,is) = ( rho%of_r(:,1) + is*sgn(is)*rho%of_r(:,nspin) )*nspin/2.d0 - &
mod(is,2)* rho_core(:)
rho%of_g(:,is) = ( rho%of_g(:,1) + is*sgn(is)*rho%of_g(:,nspin) )*nspin/2.d0 - &
mod(is,2)*rhog_core(:)
ENDDO
ENDIF
CALL stop_clock('setup_dgc')
!

6
PHonon/Gamma/dynmatcc.f90
View File

@ -21,7 +21,7 @@ SUBROUTINE dynmatcc(dyncc)
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, rhoz_or_updw
USE scf, ONLY : rho, rho_core, rhog_core
USE wavefunctions, ONLY: psic
USE cgcom
USE mp_global, ONLY : intra_pool_comm
@ -49,13 +49,9 @@ 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

27
PHonon/PH/addnlcc.f90
View File

@ -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, rhoz_or_updw
use scf, only : rho, rho_core
USE gvect, ONLY : g, ngm
USE fft_base, ONLY : dfftp
USE noncollin_module, ONLY : nspin_lsda, nspin_gga, nspin_mag
@ -96,27 +96,12 @@ 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() ) 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 (nspin_lsda == 2) call rhoz_or_updw(rho, 'only_r', 'updw_rhoz')
!^
endif
if ( dft_is_gradient() ) 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 (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)
enddo

29
PHonon/PH/addnlcc_zstar_eu_us.f90
View File

@ -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, rhoz_or_updw
USE scf, only : rho, rho_core
USE cell_base, ONLY : omega
USE gvect, ONLY : ngm, g
USE fft_base, ONLY : dfftp
@ -74,30 +74,15 @@ SUBROUTINE addnlcc_zstar_eu_us( drhoscf )
! 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() ) 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_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()) 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
IF (dft_is_nonlocc()) CALL dnonloccorr( rho%of_r, drhoscf(1, 1, ipol), &
xq, dvaux )
!
rho%of_r(:,1) = rho%of_r(:,1) - rho_core
!
DO is = 1, nspin_lsda
zstareu0(ipol,mode) = zstareu0(ipol,mode) - &
omega * fac / REAL(nrtot, DP) * &

19
PHonon/PH/dvqhub_barepsi_us.f90
View File

@ -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 : nspin, lsda, current_spin, isk
USE lsda_mod, ONLY : lsda, current_spin, isk
USE wavefunctions, ONLY : evc
USE eqv, ONLY : dvpsi
USE scf, ONLY : rho
@ -67,11 +67,9 @@ 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))
@ -252,14 +250,11 @@ 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) * rhons_m1m2 &
aux2(ig) = dqsphi(ig,ihubst1) * rho%ns(m1,m2,current_spin,nah) &
* proj1(ibnd, ihubst2)
aux4(ig) = swfcatomkpq(ig,ihubst1) * rhons_m1m2 &
aux4(ig) = swfcatomkpq(ig,ihubst1) * rho%ns(m1,m2,current_spin,nah) &
* proj2(ibnd, ihubst2)
aux5(ig) = swfcatomkpq(ig,ihubst1) &
* dnsbare(m1,m2,current_spin,nah,icart,na) &
@ -302,13 +297,11 @@ SUBROUTINE dvqhub_barepsi_us (ik, uact)
DO m2 = 1, 2*Hubbard_l(nt)+1
!
ihubst2 = offsetU(nah) + m2
!
rhons_m1m2 = ( rho%ns(m1,m2,1,nah) + sgn_op * &
rho%ns(m1,m2,nspin,nah) )*0.5d0
!
DO ig = 1, npwq
aux2(ig) = dqsphi(ig, ihubst1) * rhons_m1m2 &
aux2(ig) = dqsphi(ig, ihubst1) * rho%ns(m1,m2,op_spin,nah) &
* proj1(ibnd, ihubst2)
aux4(ig) = swfcatomkpq(ig,ihubst1) * rhons_m1m2 &
aux4(ig) = swfcatomkpq(ig,ihubst1) * rho%ns(m1,m2,op_spin,nah) &
* proj2(ibnd, ihubst2)
aux5(ig) = swfcatomkpq(ig,ihubst1) &
* dnsbare (m1,m2,op_spin,nah,icart,na) &

27
PHonon/PH/dvqhub_barepsi_us2.f90
View File

@ -32,7 +32,7 @@ SUBROUTINE dvqhub_barepsi_us2 (ik, dvqhbar, dvqhbar_orth, dvqhbar_orth_lm)
USE control_lr, ONLY : lgamma, ofsbeta
USE units_lr, ONLY : iuatwfc, iuatswfc
USE uspp_param, ONLY : nh
USE lsda_mod, ONLY : nspin, lsda, current_spin, isk
USE lsda_mod, ONLY : lsda, current_spin, isk
USE wavefunctions, ONLY : evc
USE eqv, ONLY : dvpsi
USE scf, ONLY : rho
@ -52,11 +52,8 @@ SUBROUTINE dvqhub_barepsi_us2 (ik, dvqhbar, dvqhbar_orth, dvqhbar_orth_lm)
!
COMPLEX(DP), ALLOCATABLE :: dqsphi(:,:), dmqsphi(:,:), dwfcatom_(:), dvqi(:,:), &
dvqi_orth(:,:), dvqi_orth_lm(:,:), aux1(:), aux2(:)
COMPLEX(DP) :: rhons_m1m2
INTEGER :: i, j, k, icart, na, nt, l, ih, n, mu, ig, npw, npwq, &
ihubst, ihubst1, ihubst2, nah, m, m1, m2, ibnd, op_spin, ikk, ikq, ibeta
REAL(DP) :: sgn_cs, sgn_op
COMPLEX(DP), EXTERNAL :: ZDOTC
!
CALL start_clock( 'dvqhub_barepsi_us2' )
@ -90,8 +87,6 @@ SUBROUTINE dvqhub_barepsi_us2 (ik, dvqhbar, dvqhbar_orth, dvqhbar_orth_lm)
ELSE
op_spin = 1
ENDIF
sgn_cs = DBLE(2*MOD(current_spin,2)-1)
sgn_op = DBLE(2*MOD(op_spin,2)-1)
!
! Compute beta functions at k and k+q
!
@ -250,12 +245,10 @@ SUBROUTINE dvqhub_barepsi_us2 (ik, dvqhbar, dvqhbar_orth, dvqhbar_orth_lm)
DO m2 = 1, 2*Hubbard_l(nt)+1
!
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
!
aux1(ig) = rhons_m1m2 * &
aux1(ig) = rho%ns(m1,m2,current_spin,nah) * &
( dqsphi(ig,ihubst1) * proj1(ibnd,ihubst2) + &
dqsphi(ig,ihubst2) * proj1(ibnd,ihubst1) + &
swfcatomkpq(ig,ihubst1) * proj2(ibnd,ihubst2) + &
@ -271,13 +264,13 @@ SUBROUTINE dvqhub_barepsi_us2 (ik, dvqhbar, dvqhbar_orth, dvqhbar_orth_lm)
!
DO ig = 1, npwq
!
aux1(ig) = rhons_m1m2 * &
aux1(ig) = rho%ns(m1,m2,current_spin,nah) * &
( dqsphi(ig,ihubst2) * proj1(ibnd,ihubst1) + &
swfcatomkpq(ig,ihubst2) * proj2(ibnd,ihubst1) )
!
dvqi_orth(ig,ibnd) = dvqi_orth(ig,ibnd) - aux1(ig)
!
aux1(ig) = rhons_m1m2 * &
aux1(ig) = rho%ns(m1,m2,current_spin,nah) * &
( dqsphi(ig,ihubst1) * proj1(ibnd,ihubst2) + &
swfcatomkpq(ig,ihubst1) * proj2(ibnd,ihubst2) )
!
@ -340,12 +333,10 @@ SUBROUTINE dvqhub_barepsi_us2 (ik, dvqhbar, dvqhbar_orth, dvqhbar_orth_lm)
DO m2 = 1, 2*Hubbard_l(nt)+1
!
ihubst2 = offsetU(nah) + m2
!
rhons_m1m2 = ( rho%ns(m1,m2,1,nah) + sgn_op * &
rho%ns(m1,m2,nspin,nah) )*0.5d0
!
DO ig = 1, npwq
!
aux1(ig) = rhons_m1m2 * &
aux1(ig) = rho%ns(m1,m2,op_spin,nah) * &
( dqsphi(ig,ihubst1) * proj1(ibnd,ihubst2) + &
dqsphi(ig,ihubst2) * proj1(ibnd,ihubst1) + &
swfcatomkpq(ig,ihubst1) * proj2(ibnd,ihubst2) + &
@ -363,13 +354,13 @@ SUBROUTINE dvqhub_barepsi_us2 (ik, dvqhbar, dvqhbar_orth, dvqhbar_orth_lm)
!
DO ig = 1, npwq
!
aux1(ig) = rhons_m1m2 * &
aux1(ig) = rho%ns(m1,m2,op_spin,nah) * &
( dqsphi(ig,ihubst2) * proj1(ibnd,ihubst1) + &
swfcatomkpq(ig,ihubst2) * proj2(ibnd,ihubst1) )
!
dvqi_orth(ig,ibnd) = dvqi_orth(ig,ibnd) + aux1(ig) ! sign change
!
aux1(ig) = rhons_m1m2 * &
aux1(ig) = rho%ns(m1,m2,op_spin,nah) * &
( dqsphi(ig,ihubst1) * proj1(ibnd,ihubst2) + &
swfcatomkpq(ig,ihubst1) * proj2(ibnd,ihubst2) )
!

27
PHonon/PH/dvqpsi_us.f90
View File

@ -29,7 +29,7 @@ subroutine dvqpsi_us (ik, uact, addnlcc)
ngm
USE gvecs, ONLY : ngms, doublegrid
USE lsda_mod, ONLY : lsda, isk
USE scf, ONLY : rho, rho_core, rhoz_or_updw
USE scf, ONLY : rho, rho_core
USE noncollin_module, ONLY : nspin_lsda, nspin_gga, nspin_mag, npol
use uspp_param,ONLY : upf
USE wvfct, ONLY : nbnd, npwx
@ -88,7 +88,6 @@ subroutine dvqpsi_us (ik, uact, addnlcc)
complex(DP) , allocatable, target :: aux (:)
complex(DP) , allocatable :: aux1 (:), aux2 (:)
complex(DP) , pointer :: auxs (:)
REAL(DP) :: fac
COMPLEX(DP), ALLOCATABLE :: drhoc(:)
!
call start_clock ('dvqpsi_us')
@ -174,30 +173,12 @@ subroutine dvqpsi_us (ik, uact, addnlcc)
enddo
endif
fac = 1.d0 / DBLE (nspin_lsda)
rho%of_r(:,1) = rho%of_r(:,1) + rho_core
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, drhoc, 1, nspin_gga, g, aux)
!
IF (nspin_lsda == 2) CALL rhoz_or_updw(rho, 'only_r', 'updw_rhoz')
!^
ENDIF
IF ( dft_is_gradient() ) CALL dgradcorr (dfftp, rho%of_r, grho, dvxc_rr, &
dvxc_sr, dvxc_ss, dvxc_s, xq, drhoc, 1, nspin_gga, g, aux)
IF (dft_is_nonlocc()) THEN
!^
IF (nspin_lsda == 2) CALL rhoz_or_updw(rho, 'only_r', 'rhoz_updw')
!
CALL dnonloccorr(rho%of_r, drhoc, xq, aux)
!
IF (nspin_lsda == 2) CALL rhoz_or_updw(rho, 'only_r', 'updw_rhoz')
!^
ENDIF
IF (dft_is_nonlocc()) CALL dnonloccorr(rho%of_r, drhoc, xq, aux)
rho%of_r(:,1) = rho%of_r(:,1) - rho_core

14
PHonon/PH/dynmat_hub_bare.f90
View File

@ -74,7 +74,7 @@ SUBROUTINE dynmat_hub_bare
na_icart, nap_jcart, na_icar, m1, m2, op_spin, isi, &
imode, jmode, ldim, ik, ikk, ikq, icar, ibnd, &
ios, iund2nsbare, npw, npwq
REAL(DP) :: nsaux, sgn
REAL(DP) :: nsaux
COMPLEX(DP) :: dnsaux1, dnsaux2, d2ns_bare_aux, d2ns_bare_k, work
LOGICAL :: exst
COMPLEX(DP), ALLOCATABLE :: d2ns_bare(:,:,:,:,:,:), dynwrk(:,:)
@ -451,15 +451,12 @@ SUBROUTINE dynmat_hub_bare
work = (0.d0, 0.d0)
!
DO is = 1, nspin
!
sgn = REAL( 2*MOD(is,2) - 1 )
!
DO m1 = 1, 2*Hubbard_l(nt) + 1
!
DO m2 = 1, 2*Hubbard_l(nt) + 1
!
nsaux = ( rho%ns(m1, m2, 1, nah) + sgn * &
rho%ns(m1, m2, nspin, nah) )*0.5d0
nsaux = rho%ns(m1, m2, is, nah)
!
! dnsbare matrix is symmetric i.e. dnsbare(m1,m2) = dnsbare(m2,m1)
! (when keeping only the terms in u and neglecting the ones in u*)
@ -502,14 +499,12 @@ SUBROUTINE dynmat_hub_bare
DO is = 1, nspin
!
isi = nspin/is
sgn = REAL( 2*MOD(isi,2) - 1 )
!
DO m1 = 1, 2*Hubbard_l(nt) + 1
!
DO m2 = 1, 2*Hubbard_l(nt) + 1
!
nsaux = ( rho%ns(m1, m2, 1, nah) + sgn * &
rho%ns(m1, m2, nspin, nah) )*0.5d0
nsaux = rho%ns(m1, m2, isi, nah)
!
! dnsbare matrix is symmetric i.e. dnsbare(m1,m2) = dnsbare(m2,m1)
! (when keeping only the terms in u and neglecting the ones in u*)
@ -521,7 +516,7 @@ SUBROUTINE dynmat_hub_bare
! DO NOT include the delta_m1m2 contribution
! Note the sign change
!
work = work + nsaux * d2ns_bare_aux &
work = work + nsaux * d2ns_bare_aux + &
+ dnsaux1 * CONJG(dnsaux2)
!
ENDDO ! m2
@ -601,4 +596,3 @@ SUBROUTINE dynmat_hub_bare
RETURN
!
END SUBROUTINE dynmat_hub_bare

6
PHonon/PH/dynmatcc.f90
View File

@ -20,7 +20,7 @@ subroutine dynmatcc
USE fft_interfaces, ONLY : fwfft
USE gvect, ONLY : ngm, g
USE lsda_mod, ONLY : nspin
use scf, ONLY : rho, rho_core, rhog_core, rhoz_or_updw
use scf, ONLY : rho, rho_core, rhog_core
USE modes, ONLY : u
USE qpoint, ONLY : xq
USE nlcc_ph, ONLY : drc
@ -49,13 +49,9 @@ subroutine dynmatcc
!
allocate (vxc( dfftp%nnr))
allocate (v ( dfftp%nnr , nspin))
!^
CALL rhoz_or_updw(rho, 'r_and_g', 'rhoz_updw')
!
call v_xc (rho, rho_core, rhog_core, etxcd, vtxcd, v)
!
CALL rhoz_or_updw(rho, 'r_and_g', 'updw_rhoz')
!^
if (nspin == 1 .OR. nspin==4) then
is=1
do ir = 1, dfftp%nnr

7
PHonon/PH/phq_setup.f90
View File

@ -56,7 +56,7 @@ subroutine phq_setup
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, kedtau, rho, rhoz_or_updw
USE scf, ONLY : v, vrs, vltot, kedtau, rho
USE fft_base, ONLY : dfftp
USE gvect, ONLY : ngm
USE gvecs, ONLY : doublegrid
@ -174,14 +174,9 @@ subroutine phq_setup
call setup_dmuxc()
!
! 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()

6
PP/src/add_shift_cc.f90
View File

@ -21,7 +21,7 @@ SUBROUTINE add_shift_cc (shift_cc)
USE gvect, ONLY: ngm, gstart, g, gg, ngl, gl, igtongl
USE ener, ONLY: etxc, vtxc
USE lsda_mod, ONLY: nspin
USE scf, ONLY: rho, rho_core, rhog_core, rhoz_or_updw
USE scf, ONLY: rho, rho_core, rhog_core
USE control_flags, ONLY: gamma_only
USE wavefunctions, ONLY : psic
USE mp_global, ONLY : intra_pool_comm
@ -62,13 +62,9 @@ SUBROUTINE add_shift_cc (shift_cc)
!
ALLOCATE ( vxc(dfftp%nnr,nspin), shift_(nat) )
shift_(:) = 0.d0
!^
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')
!^
IF (nspin==1) THEN
DO ir = 1, dfftp%nnr
psic (ir) = vxc (ir, 1)

16
PP/src/local_dos.f90
View File

@ -34,7 +34,7 @@ SUBROUTINE local_dos (iflag, lsign, kpoint, kband, spin_component, &
USE klist, ONLY : lgauss, degauss, ngauss, nks, wk, xk, &
nkstot, ngk, igk_k
USE lsda_mod, ONLY : lsda, nspin, current_spin, isk
USE scf, ONLY : rho, rhoz_or_updw
USE scf, ONLY : rho
USE symme, ONLY : sym_rho, sym_rho_init, sym_rho_deallocate
USE uspp, ONLY : nkb, vkb, becsum, nhtol, nhtoj, indv
USE uspp_param, ONLY : upf, nh, nhm
@ -372,13 +372,7 @@ SUBROUTINE local_dos (iflag, lsign, kpoint, kband, spin_component, &
!
! Here we add the US contribution to the charge
!
!^
IF ( nspin==2 ) CALL rhoz_or_updw(rho, 'only_r', 'rhoz_updw')
!
CALL addusdens(rho%of_r(:,:))
!
IF ( nspin==2 ) CALL rhoz_or_updw(rho, 'only_r', 'updw_rhoz')
!^
CALL addusdens(rho%of_r(:,:), 'rho-mz')
!
IF (nspin == 1 .or. nspin==4) THEN
is = 1
@ -408,14 +402,8 @@ SUBROUTINE local_dos (iflag, lsign, kpoint, kband, spin_component, &
CALL fwfft ('Rho', psic, dfftp)
rho%of_g(:,1) = psic(dfftp%nl(:))
!
!^
IF ( nspin==2 ) CALL rhoz_or_updw(rho, 'only_g', 'rhoz_updw')
!
CALL sym_rho (1, rho%of_g)
!
IF ( nspin==2 ) CALL rhoz_or_updw(rho, 'only_g', 'updw_rhoz')
!^
!
psic(:) = (0.0_dp, 0.0_dp)
psic(dfftp%nl(:)) = rho%of_g(:,1)
CALL invfft ('Rho', psic, dfftp)

26
PP/src/paw_postproc.f90
View File

@ -148,14 +148,26 @@ SUBROUTINE PAW_make_ae_charge(rho,withcore)
!
! prepare spherical harmonics
CALL ylmr2( i%l**2, 1, posi, distsq, ylm_posi )
DO is = 1,nspin
DO lm = 1, i%l**2
! do interpolation
rho%of_r(ir,is)= rho%of_r(ir,is) + ylm_posi(1,lm) &
* splint(g(i%t)%r(:) , rho_lm(:,lm,is), &
wsp_lm(:,lm,is), sqrt(distsq) )
IF ( nspin/=2 ) THEN
DO is = 1,nspin
DO lm = 1, i%l**2
! do interpolation
rho%of_r(ir,is)= rho%of_r(ir,is) + ylm_posi(1,lm) &
* splint(g(i%t)%r(:) , rho_lm(:,lm,is), &
wsp_lm(:,lm,is), sqrt(distsq) )
ENDDO
ENDDO
ENDDO
ELSE
DO lm = 1, i%l**2
! do interpolation
rho%of_r(ir,1)= rho%of_r(ir,is) + ylm_posi(1,lm) &
* splint(g(i%t)%r(:) , rho_lm(:,lm,is), &
wsp_lm(:,lm,is), sqrt(distsq) )
rho%of_r(ir,2)= rho%of_r(ir,is) + ylm_posi(1,lm)*(2*mod(is,2)-1) &
* splint(g(i%t)%r(:) , rho_lm(:,lm,is), &
wsp_lm(:,lm,is), sqrt(distsq) )
ENDDO
ENDIF
ENDDO rsp_point
!
DEALLOCATE(rho_lm, ylm_posi, wsp_lm)

8
PP/src/ppacf.f90
View File

@ -46,7 +46,7 @@ PROGRAM do_ppacf
USE lsda_mod, ONLY : nspin
USE scf, ONLY : scf_type,create_scf_type,destroy_scf_type
USE scf, ONLY : scf_type_COPY
USE scf, ONLY : rho, rho_core, rhog_core, vltot, rhoz_or_updw
USE scf, ONLY : rho, rho_core, rhog_core, vltot
USE funct, ONLY : xc,xc_spin,gcxc,gcx_spin,gcc_spin,dft_is_nonlocc,nlc
USE funct, ONLY : get_iexch, get_icorr, get_igcx, get_igcc
USE funct, ONLY : set_exx_fraction,set_auxiliary_flags,enforce_input_dft
@ -246,13 +246,9 @@ PROGRAM do_ppacf
etxcccnl=0._DP
vtxcccnl=0._DP
vofrcc=0._DP
!^
IF ( nspin==2 ) CALL rhoz_or_updw( rho, 'only_r', 'rhoz_updw' )
!
CALL nlc( rho%of_r, rho_core, nspin, etxcccnl, vtxcccnl, vofrcc )
!
IF ( nspin==2 ) CALL rhoz_or_updw( rho, 'only_r', 'updw_rhoz' )
!^
CALL mp_sum( etxcccnl , intra_bgrp_comm )
END IF
!
@ -284,7 +280,7 @@ PROGRAM do_ppacf
!
IF (nspin == 1) THEN
tot_rho(:)=rhoout(:,1)
ELSEIF(nspin==2) THEN !^ rhoout is up/dw (for now)
ELSEIF(nspin==2) THEN ! rhoout is (up,down)
tot_rho(:)=rhoout(:,1)+rhoout(:,2)
ELSE
CALL errore ('ppacf','vdW-DF not available for noncollinear spin case',1)

9
PP/src/punch_plot.f90
View File

@ -33,7 +33,7 @@ SUBROUTINE punch_plot (filplot, plot_num, sample_bias, z, dz, &
USE lsda_mod, ONLY : nspin, lsda
USE ener, ONLY : ehart
USE io_global, ONLY : stdout, ionode
USE scf, ONLY : rho, vltot, v, rhoz_or_updw
USE scf, ONLY : rho, vltot, v
USE wvfct, ONLY : nbnd, wg
USE gvecw, ONLY : ecutwfc
USE noncollin_module, ONLY : noncolin
@ -232,13 +232,8 @@ SUBROUTINE punch_plot (filplot, plot_num, sample_bias, z, dz, &
WRITE(stdout, '(7x,a)') "Reconstructing all-electron valence charge."
! code partially duplicate from plot_num=0, should be unified
CALL init_us_1()
!^
IF ( nspin==2 ) CALL rhoz_or_updw(rho, 'r_and_g', 'rhoz_updw')
!^
!
CALL PAW_make_ae_charge(rho,(plot_num==21))
!^
IF ( nspin==2 ) CALL rhoz_or_updw(rho, 'r_and_g', 'updw_rhoz')
!^
!
raux(:) = rho%of_r(:, 1)
IF ( lsda ) THEN

26
PP/src/pw2bgw.f90
View File

@ -1347,7 +1347,7 @@ SUBROUTINE write_rhog ( output_file_name, real_or_complex, symm_type, &
CALL calc_rhog ( rhog_nvmin, rhog_nvmax )
!
IF ( nspin==2 ) CALL rhoz_or_updw(rho, 'r_and_g', 'updw_rhoz')
!^
!
ALLOCATE ( g_g ( nd, ng_g ) )
ALLOCATE ( rhog_g ( ng_g, ns ) )
@ -1508,7 +1508,7 @@ SUBROUTINE write_vxcg ( output_file_name, real_or_complex, symm_type, &
USE lsda_mod, ONLY : nspin
USE mp, ONLY : mp_sum
USE mp_bands, ONLY : intra_bgrp_comm
USE scf, ONLY : rho, rho_core, rhog_core, rhoz_or_updw
USE scf, ONLY : rho, rho_core, rhog_core
USE symm_base, ONLY : s, ftau, nsym
USE wavefunctions, ONLY : psic
USE matrix_inversion
@ -1656,13 +1656,9 @@ SUBROUTINE write_vxcg ( output_file_name, real_or_complex, symm_type, &
rho_core ( : ) = 0.0D0
rhog_core ( : ) = ( 0.0D0, 0.0D0 )
ENDIF
!^
IF ( nspin==2 ) CALL rhoz_or_updw(rho, 'r_and_g', 'rhoz_updw')
!
CALL v_xc ( rho, rho_core, rhog_core, etxc, vtxc, vxcr_g )
!
IF ( nspin==2 ) CALL rhoz_or_updw(rho, 'r_and_g', 'updw_rhoz')
!^
DO is = 1, ns
DO ir = 1, nr
psic ( ir ) = CMPLX ( vxcr_g ( ir, is ), 0.0D0, KIND=dp )
@ -1726,7 +1722,7 @@ SUBROUTINE write_vxc0 ( output_file_name, vxc_zero_rho_core )
USE lsda_mod, ONLY : nspin
USE mp, ONLY : mp_sum
USE mp_bands, ONLY : intra_bgrp_comm
USE scf, ONLY : rho, rho_core, rhog_core, rhoz_or_updw
USE scf, ONLY : rho, rho_core, rhog_core
USE wavefunctions, ONLY : psic
IMPLICIT NONE
@ -1759,13 +1755,9 @@ SUBROUTINE write_vxc0 ( output_file_name, vxc_zero_rho_core )
rho_core ( : ) = 0.0D0
rhog_core ( : ) = ( 0.0D0, 0.0D0 )
ENDIF
!^
IF ( nspin==2 ) CALL rhoz_or_updw(rho, 'r_and_g', 'rhoz_updw')
!
CALL v_xc ( rho, rho_core, rhog_core, etxc, vtxc, vxcr_g )
!
IF ( nspin==2 ) CALL rhoz_or_updw(rho, 'r_and_g', 'updw_rhoz')
!^
DO is = 1, ns
DO ir = 1, nr
psic ( ir ) = CMPLX ( vxcr_g ( ir, is ), 0.0D0, KIND=dp )
@ -1826,7 +1818,7 @@ SUBROUTINE write_vxc_r (output_file_name, diag_nmin, diag_nmax, &
USE mp, ONLY : mp_sum
USE mp_pools, ONLY : inter_pool_comm
USE mp_bands, ONLY : intra_bgrp_comm
USE scf, ONLY : rho, rho_core, rhog_core, rhoz_or_updw
USE scf, ONLY : rho, rho_core, rhog_core
USE wavefunctions, ONLY : evc, psic
USE wvfct, ONLY : nbnd
@ -1892,13 +1884,9 @@ SUBROUTINE write_vxc_r (output_file_name, diag_nmin, diag_nmax, &
rho_core ( : ) = 0.0D0
rhog_core ( : ) = ( 0.0D0, 0.0D0 )
ENDIF
!^
IF ( nspin==2 ) CALL rhoz_or_updw(rho, 'r_and_g', 'rhoz_updw')
!
CALL v_xc (rho, rho_core, rhog_core, etxc, vtxc, vxcr)
!
IF ( nspin==2 ) CALL rhoz_or_updw(rho, 'r_and_g', 'updw_rhoz')
!^
DO ik = iks, ike
npw = ngk ( ik - iks + 1 )
CALL davcio (evc, 2*nwordwfc, iunwfc, ik - iks + 1, -1)
@ -2018,7 +2006,7 @@ SUBROUTINE write_vxc_g (output_file_name, diag_nmin, diag_nmax, &
USE mp, ONLY : mp_sum
USE mp_pools, ONLY : inter_pool_comm
USE mp_bands, ONLY : intra_bgrp_comm
USE scf, ONLY : rho, rho_core, rhog_core, rhoz_or_updw
USE scf, ONLY : rho, rho_core, rhog_core
USE wavefunctions, ONLY : evc, psic
USE wvfct, ONLY : npwx, nbnd
@ -2085,13 +2073,9 @@ SUBROUTINE write_vxc_g (output_file_name, diag_nmin, diag_nmax, &
rho_core ( : ) = 0.0D0
rhog_core ( : ) = ( 0.0D0, 0.0D0 )
ENDIF
!^
IF ( nspin==2 ) CALL rhoz_or_updw(rho, 'r_and_g', 'rhoz_updw')
!
CALL v_xc (rho, rho_core, rhog_core, etxc, vtxc, vxcr)
!
IF ( nspin==2 ) CALL rhoz_or_updw(rho, 'r_and_g', 'updw_rhoz')
!^
DO ik = iks, ike
ikk = ik - iks + 1
npw = ngk ( ik - iks + 1 )

6
PP/src/pw2gw.f90
View File

@ -152,7 +152,7 @@ SUBROUTINE compute_gw( omegamin, omegamax, d_omega, use_gmaps, qplda, vkb, vxcdi
USE mp_world, ONLY : world_comm, mpime, nproc
USE mp_wave, ONLY : mergewf
USE parallel_include
USE scf, ONLY : rho, rho_core, rhog_core, rhoz_or_updw
USE scf, ONLY : rho, rho_core, rhog_core
USE ener, ONLY : etxc, vtxc
USE uspp_param, ONLY : upf, nh
@ -790,13 +790,9 @@ SUBROUTINE compute_gw( omegamin, omegamax, d_omega, use_gmaps, qplda, vkb, vxcdi
OPEN (UNIT=313,FILE="vxcdiag.dat",STATUS="UNKNOWN")
WRITE(313,*) "# K BND <Vxc>"
ALLOCATE ( vxc(dfftp%nnr,nspin) )
!^
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')
!^
DO ik=1,nkpt
npw = ngk(ik)
CALL davcio( evc, 2*nwordwfc, iunwfc, ik, -1 )

26
PW/src/addusdens.f90
View File

@ -5,9 +5,8 @@
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!----------------------------------------------------------------------
SUBROUTINE addusdens(rho)
SUBROUTINE addusdens(rho,lsda_format)
!----------------------------------------------------------------------
!
! ... Add US contribution to the charge density to rho(G)
@ -20,13 +19,19 @@ SUBROUTINE addusdens(rho)
!
IMPLICIT NONE
!
!
COMPLEX(kind=dp), INTENT(inout) :: rho(dfftp%ngm,nspin_mag)
CHARACTER(len=*), OPTIONAL :: lsda_format
INTEGER :: sw_lsda
sw_lsda = 0
IF ( present(lsda_format) .and. nspin_mag==2 ) THEN
IF (lsda_format == 'rho-mz') sw_lsda = 1
ENDIF
!
IF ( tqr ) THEN
CALL addusdens_r(rho)
CALL addusdens_r(rho,sw_lsda)
ELSE
CALL addusdens_g(rho)
CALL addusdens_g(rho,sw_lsda)
ENDIF
!
RETURN
@ -34,7 +39,7 @@ SUBROUTINE addusdens(rho)
END SUBROUTINE addusdens
!
!----------------------------------------------------------------------
SUBROUTINE addusdens_g(rho)
SUBROUTINE addusdens_g(rho, sw_lsda)
!----------------------------------------------------------------------
!
! This routine adds to the charge density rho(G) in reciprocal space
@ -56,6 +61,7 @@ SUBROUTINE addusdens_g(rho)
IMPLICIT NONE
!
COMPLEX(kind=dp), INTENT(inout) :: rho(dfftp%ngm,nspin_mag)
INTEGER, INTENT(in) :: sw_lsda
!
! here the local variables
!
@ -63,7 +69,6 @@ SUBROUTINE addusdens_g(rho)
! starting/ending indices, local number of G-vectors
INTEGER :: ig, na, nt, ih, jh, ijh, is, nab, nb, nij
! counters
REAL(DP), ALLOCATABLE :: tbecsum(:,:,:)
! \sum_kv <\psi_kv|\beta_l><beta_m|\psi_kv> for each species of atoms
REAL(DP), ALLOCATABLE :: qmod (:), ylmk0 (:,:)
@ -159,7 +164,12 @@ SUBROUTINE addusdens_g(rho)
!
! add aux to the charge density in reciprocal space
!
rho(:,:) = rho(:,:) + aux (:,:)
IF (sw_lsda == 0) THEN
rho(:,:) = rho(:,:) + aux(:,:)
ELSE
rho(:,1) = rho(:,1) + aux(:,1) + aux(:,2)
rho(:,2) = rho(:,2) + aux(:,1) - aux(:,2)
ENDIF
!
DEALLOCATE (aux)
!

2
PW/src/atomic_rho.f90
View File

@ -92,7 +92,7 @@ SUBROUTINE atomic_rho_g (rhocg, nspina)
rhoscale = 1.0_dp
ENDIF
!
!^
!
rhocg(:,1) = rhocg(:,1) + &
strf(:,nt) * rhoscale * rhocgnt(igtongl(:)) / omega
!

49
PW/src/electrons.f90
View File

@ -810,17 +810,11 @@ SUBROUTINE electrons_scf ( printout, exxen )
END IF
!
! calculate the xdm energy contribution with converged density
if (lxdm .and. conv_elec) then
!^
IF (nspin == 2) CALL rhoz_or_updw( rho, 'r_and_g', 'rhoz_updw' )
!
IF (lxdm .and. conv_elec) THEN
exdm = energy_xdm()
!
IF (nspin == 2) CALL rhoz_or_updw( rho, 'r_and_g', 'updw_rhoz' )
!^
etot = etot + exdm
hwf_energy = hwf_energy + exdm
end if
END IF
IF (ts_vdw) THEN
! factor 2 converts from Ha to Ry units
etot = etot + 2.0d0*EtsvdW
@ -997,13 +991,13 @@ SUBROUTINE electrons_scf ( printout, exxen )
IMPLICIT NONE
REAL(DP) :: delta_e, delta_e_hub
INTEGER :: ir
!^
!
delta_e = 0._dp
IF ( nspin==2 ) THEN
!
DO ir = 1,dfftp%nnr
delta_e = delta_e - ( rho%of_r(ir,1) + rho%of_r(ir,2) ) * v%of_r(ir,1) & !^up
- ( rho%of_r(ir,1) - rho%of_r(ir,2) ) * v%of_r(ir,2) !^dw
delta_e = delta_e - ( rho%of_r(ir,1) + rho%of_r(ir,2) ) * v%of_r(ir,1) & ! up
- ( rho%of_r(ir,1) - rho%of_r(ir,2) ) * v%of_r(ir,2) ! dw
ENDDO
delta_e = 0.5_dp*delta_e
!
@ -1023,8 +1017,7 @@ SUBROUTINE electrons_scf ( printout, exxen )
delta_e_hub = - SUM (rho%ns_nc(:,:,:,:)*v%ns_nc(:,:,:,:))
delta_e = delta_e + delta_e_hub
else
delta_e_hub = - SUM( (rho%ns(:,:,1,:)+rho%ns(:,:,nspin,:))*v%ns(:,:,1,:) + & !up +
(rho%ns(:,:,1,:)-rho%ns(:,:,nspin,:))*v%ns(:,:,nspin,:) )*0.5d0 !down
delta_e_hub = - SUM (rho%ns(:,:,:,:)*v%ns(:,:,:,:))
if (nspin==1) delta_e_hub = 2.d0 * delta_e_hub
delta_e = delta_e + delta_e_hub
endif
@ -1049,10 +1042,9 @@ SUBROUTINE electrons_scf ( printout, exxen )
!
USE funct, ONLY : dft_is_meta
IMPLICIT NONE
REAL(DP) :: delta_escf, delta_escf_hub
REAL(DP) :: rho_dif(2)
REAL(DP) :: delta_escf, delta_escf_hub, rho_dif(2)
INTEGER :: ir
!^
!
delta_escf=0._dp
IF ( nspin==2 ) THEN
!
@ -1077,24 +1069,19 @@ SUBROUTINE electrons_scf ( printout, exxen )
!
CALL mp_sum( delta_escf, intra_bgrp_comm )
!
if (lda_plus_u) then
if (noncolin) then
IF ( lda_plus_u ) THEN
IF ( noncolin ) THEN
delta_escf_hub = -SUM((rhoin%ns_nc(:,:,:,:)-rho%ns_nc(:,:,:,:))*v%ns_nc(:,:,:,:))
delta_escf = delta_escf + delta_escf_hub
else
!
delta_escf_hub = -SUM( (rhoin%ns(:,:,1,:) + rhoin%ns(:,:,nspin,:) - & !up
( rho%ns(:,:,1,:) + rho%ns(:,:,nspin,:)) )*v%ns(:,:,1,:) + &
(rhoin%ns(:,:,1,:) - rhoin%ns(:,:,nspin,:) - & !down
( rho%ns(:,:,1,:) - rho%ns(:,:,nspin,:)) )*v%ns(:,:,nspin,:) )
!
if (nspin==1) delta_escf_hub = 2.d0 * delta_escf_hub
delta_escf = delta_escf + delta_escf_hub * 0.5d0
endif
end if
ELSE
delta_escf_hub = -SUM((rhoin%ns(:,:,:,:)-rho%ns(:,:,:,:))*v%ns(:,:,:,:))
IF ( nspin==1 ) delta_escf_hub = 2.d0 * delta_escf_hub
delta_escf = delta_escf + delta_escf_hub
ENDIF
ENDIF
IF (okpaw) delta_escf = delta_escf - &
SUM(ddd_paw(:,:,:)*(rhoin%bec(:,:,:)-rho%bec(:,:,:)))
IF ( okpaw ) delta_escf = delta_escf - &
SUM(ddd_paw(:,:,:)*(rhoin%bec(:,:,:)-rho%bec(:,:,:)))
RETURN
!

6
PW/src/force_cc.f90
View File

@ -21,7 +21,7 @@ subroutine force_cc (forcecc)
USE gvect, ONLY : ngm, gstart, g, gg, ngl, gl, igtongl
USE ener, ONLY : etxc, vtxc
USE lsda_mod, ONLY : nspin
USE scf, ONLY : rho, rho_core, rhog_core, rhoz_or_updw
USE scf, ONLY : rho, rho_core, rhog_core
USE control_flags, ONLY : gamma_only
USE noncollin_module, ONLY : noncolin
USE wavefunctions, ONLY : psic
@ -64,12 +64,8 @@ subroutine force_cc (forcecc)
!
allocate ( vxc(dfftp%nnr,nspin) )
!
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')
!
psic=(0.0_DP,0.0_DP)
if (nspin == 1 .or. nspin == 4) then
do ir = 1, dfftp%nnr

13
PW/src/gradcorr.f90
View File

@ -39,7 +39,7 @@ SUBROUTINE gradcorr( rho, rhog, rho_core, rhog_core, etxc, vtxc, v )
COMPLEX(DP), ALLOCATABLE :: rhogsum(:,:)
!
REAL(DP) :: grho2(2), sx, sc, v1x, v2x, v1c, v2c, &
REAL(DP) :: grho2(2), sgn(2), sx, sc, v1x, v2x, v1c, v2c, &
v1xup, v1xdw, v2xup, v2xdw, v1cup, v1cdw , &
etxcgc, vtxcgc, segno, arho, fac, zeta, rh, grh2, amag
REAL(DP) :: v2cup, v2cdw, v2cud, rup, rdw, &
@ -57,6 +57,7 @@ SUBROUTINE gradcorr( rho, rhog, rho_core, rhog_core, etxc, vtxc, v )
if (nspin==4) nspin0=1
if (nspin==4.and.domag) nspin0=2
fac = 1.D0 / DBLE( nspin0 )
sgn(1) = 1.D0 ; sgn(2) = -1.D0
!
ALLOCATE( h( 3, dfftp%nnr, nspin0) )
ALLOCATE( grho( 3, dfftp%nnr, nspin0) )
@ -90,13 +91,17 @@ SUBROUTINE gradcorr( rho, rhog, rho_core, rhog_core, etxc, vtxc, v )
END DO
ELSE
!
rhoout(:,1:nspin0) = rho(:,1:nspin0)
rhogsum(:,1:nspin0) = rhog(:,1:nspin0)
! ... for convenience rhoout and rhogsum are in (up,down) format, if LSDA
!
DO is = 1, nspin0
rhoout(:,is) = ( rho(:,1) + sgn(is) * rho(:,nspin0) ) * 0.5D0
rhogsum(:,is) = ( rhog(:,1) + sgn(is) * rhog(:,nspin0) ) * 0.5D0
ENDDO
!
ENDIF
DO is = 1, nspin0
!
rhoout(:,is) = fac * rho_core(:) + rhoout(:,is)
rhoout(:,is) = fac * rho_core(:) + rhoout(:,is)
rhogsum(:,is) = fac * rhog_core(:) + rhogsum(:,is)
!
CALL fft_gradient_g2r( dfftp, rhogsum(1,is), g, grho(1,1,is) )

18
PW/src/init_ns.f90
View File

@ -27,7 +27,6 @@ subroutine init_ns
implicit none
real(DP) :: totoc
real(DP) :: ns_updw(2)
real(DP), external :: hubbard_occ
integer :: ldim, na, nt, is, m1, majs, mins
@ -55,27 +54,18 @@ subroutine init_ns
if (.not.nm) then
if (totoc.gt.ldim) then
do m1 = 1, ldim
ns_updw(majs) = 1.d0
ns_updw(mins) = (totoc - ldim) / ldim
rho%ns(m1, m1, 1, na) = ns_updw(1) + ns_updw(2)
rho%ns(m1, m1, 2, na) = ns_updw(1) - ns_updw(2)
rho%ns (m1, m1, majs, na) = 1.d0
rho%ns (m1, m1, mins, na) = (totoc - ldim) / ldim
enddo
else
do m1 = 1, ldim
ns_updw(majs) = totoc / ldim
ns_updw(mins) = 0.d0
rho%ns(m1, m1, 1, na) = ns_updw(1) + ns_updw(2)
if (nspin == 2) rho%ns(m1, m1, 2, na) = ns_updw(1) - ns_updw(2)
rho%ns (m1, m1, majs, na) = totoc / ldim
enddo
endif
else
do is = 1,nspin
do m1 = 1, ldim
rho%ns (m1, m1, 1, na) = totoc / 2.d0 / ldim
if (nspin == 2) then
rho%ns (m1, m1, 1, na) = rho%ns (m1, m1, 1, na) * 2.d0
rho%ns (m1, m1, 2, na) = 0.d0
endif
rho%ns (m1, m1, is, na) = totoc / 2.d0 / ldim
enddo
enddo
endif

8
PW/src/ns_adj.f90
View File

@ -14,7 +14,7 @@ subroutine ns_adj
USE kinds, ONLY : DP
USE ions_base, ONLY : nat, ntyp => nsp, ityp
USE ldaU, ONLY : Hubbard_lmax, Hubbard_l, Hubbard_U, starting_ns
USE scf, ONLY : rho, rhoz_or_updw
USE scf, ONLY : rho
USE lsda_mod, ONLY : nspin
USE noncollin_module, ONLY : noncolin, npol
USE io_global, ONLY : stdout
@ -74,8 +74,7 @@ subroutine ns_adj
end do
else
!^
IF (nspin == 2) call rhoz_or_updw(rho, 'hub_ns', 'rhoz_updw')
!
do is = 1, nspin
do m1 = 1, ldim
do m2 = 1, ldim
@ -97,8 +96,7 @@ subroutine ns_adj
enddo
enddo
enddo
IF (nspin == 2) call rhoz_or_updw(rho, 'hub_ns', 'updw_rhoz')
!^
!
endif
endif

23
PW/src/realus.f90
View File

@ -1143,7 +1143,7 @@ MODULE realus
END SUBROUTINE newq_r
!
!------------------------------------------------------------------------
SUBROUTINE addusdens_r( rho )
SUBROUTINE addusdens_r( rho, sw_lsda )
!------------------------------------------------------------------------
!
! ... This routine adds to the charge density the part which is due to
@ -1169,18 +1169,21 @@ MODULE realus
!
IMPLICIT NONE
! The charge density to be augmented (in G-space)
COMPLEX(kind=dp), INTENT(inout) :: rho(dfftp%ngm,nspin_mag)
COMPLEX(kind=dp), INTENT(inout) :: rho(dfftp%ngm,nspin_mag)
INTEGER, INTENT(in) :: sw_lsda
!
INTEGER :: ia, nt, ir, irb, ih, jh, ijh, is, mbia
CHARACTER(len=80) :: msg
REAL(kind=dp), ALLOCATABLE :: rhor(:,:)
REAL(DP) :: charge
REAL(DP) :: charge, sgn(2)
!
!
IF ( .not. okvan ) RETURN
!
CALL start_clock( 'addusdens' )
!
sgn(1)=1._dp ; sgn(2)=-1._dp
!
ALLOCATE ( rhor(dfftp%nnr,nspin_mag) )
rhor(:,:) = 0.0_dp
DO is = 1, nspin_mag
@ -1207,10 +1210,16 @@ MODULE realus
!
ENDDO
!
!
DO is = 1, nspin_mag
psic(:) = rhor(:,is)
CALL fwfft ('Rho', psic, dfftp)
rho(:,is) = rho(:,is) + psic(dfftp%nl(:))
IF (sw_lsda == 0) THEN
rho(:,is) = rho(:,is) + psic(dfftp%nl(:))
ELSE
rho(:,1) = rho(:,1) + psic(dfftp%nl(:))
rho(:,2) = rho(:,2) + sgn(is)*psic(dfftp%nl(:))
ENDIF
END DO
!!! CALL rho_r2g(dfftp, rhor, rho(:,1:nspin_mag) )
!
@ -1220,7 +1229,11 @@ MODULE realus
! ... check the total charge (must not be summed on k-points)
!
IF ( gstart == 2) THEN
charge = SUM(rho(1,1:nspin_lsda) )*omega
IF (sw_lsda == 0) THEN
charge = SUM(rho(1,1:nspin_lsda) )*omega
ELSE
charge = SUM( rho(1,1) )*omega
ENDIF
ELSE
charge = 0.0_dp
ENDIF

19
PW/src/scf_mod.f90
View File

@ -705,9 +705,8 @@ END FUNCTION ns_ddot
!
SUBROUTINE rhoz_or_updw( rho, sp, dir )
!--------------------------------------------------------------------------
! ... Converts rho_up/dw into rho_tot and m_z if dir='updw_rhoz' and
! ... vice versa if dir='rhoz_updw'.
! ... (When conversion will be full, it should become almost useless)
! ... Converts rho(up,dw) into rho(up+dw,up-dw) if dir='->rhoz' and
! ... vice versa if dir='->updw'.
!
USE gvect, ONLY : ngm
!
@ -721,11 +720,11 @@ END FUNCTION ns_ddot
IF ( nspin /= 2 ) RETURN
!
vi = 0._dp
IF (dir == 'rhoz_updw') vi = 0.5_dp
IF (dir == 'updw_rhoz') vi = 1.0_dp
IF (dir == '->updw') vi = 0.5_dp
IF (dir == '->rhoz') vi = 1.0_dp
IF (vi == 0._dp) CALL errore( 'rhoz_or_updw', 'wrong input', 1 )
!
IF ( sp=='only_r' .or. sp=='r_and_g' ) THEN
IF ( sp /= 'only_g' ) THEN
! !$omp parallel do
DO ir = 1, dfftp%nnr
rho%of_r(ir,1) = ( rho%of_r(ir,1) + rho%of_r(ir,nspin) ) * vi
@ -733,18 +732,14 @@ END FUNCTION ns_ddot
ENDDO
! !$omp end parallel do
ENDIF
IF ( sp=='only_g' .or. sp=='r_and_g' ) THEN
! !$omp parallel do
IF ( sp /= 'only_r' ) THEN
! !$omp parallel do
DO ir = 1, ngm
rho%of_g(ir,1) = ( rho%of_g(ir,1) + rho%of_g(ir,nspin) ) * vi
rho%of_g(ir,nspin) = rho%of_g(ir,1) - rho%of_g(ir,nspin) * vi * 2._dp
ENDDO
! !$omp end parallel do
ENDIF
IF ( sp=='hub_ns' ) THEN
rho%ns(:,:,1,:) = ( rho%ns(:,:,1,:) + rho%ns(:,:,2,:) ) * vi
rho%ns(:,:,2,:) = rho%ns(:,:,1,:) - rho%ns(:,:,2,:) * vi * 2._dp
ENDIF
!
RETURN
!

29
PW/src/stres_gradcorr.f90
View File

@ -30,16 +30,16 @@ subroutine stres_gradcorr( rho, rhog, rho_core, rhog_core, kedtau, nspin, &
complex(DP), intent(inout) :: rhog(dfft%ngm, nspin)
! FIXME: should be intent(in)
complex(DP), intent(in) :: rhog_core(dfft%ngm)
real(dp), intent(inout) :: sigmaxc (3, 3)
real(dp), intent(inout) :: sigmaxc(3, 3)
!
integer :: k, l, m, ipol, is, nspin0
integer :: nr1, nr2, nr3, nrxx, ngm
real(DP) , allocatable :: grho (:,:,:)
real(DP), parameter :: epsr = 1.0d-6, epsg = 1.0d-10, e2 = 2.d0
real(DP) :: grh2, grho2 (2), sx, sc, v1x, v2x, v1c, v2c, fac, &
real(DP) :: grh2, grho2(2), sgn(2), fac(2), sx, sc, v1x, v2x, v1c, v2c, &
v1xup, v1xdw, v2xup, v2xdw, v1cup, v1cdw, v2cup, v2cdw, v2cud, &
zeta, rh, rup, rdw, grhoup, grhodw, grhoud, grup, grdw, &
sigma_gradcorr (3, 3)
sigma_gradcorr(3, 3)
logical :: igcc_is_lyp
!dummy variables for meta-gga
real(DP) :: v3x,v3c,v3xup,v3xdw,v3cup,v3cdw
@ -64,18 +64,20 @@ subroutine stres_gradcorr( rho, rhog, rho_core, rhog_core, kedtau, nspin, &
allocate (grho( 3, nrxx, nspin))
nspin0=nspin
if (nspin==4) nspin0=1
fac = 1.d0 / DBLE (nspin0)
fac(1) = dble(mod(nspin0,2)+1) ; fac(2) = 0.d0
sgn(1) = 1.d0 ; sgn(2) = -1.d0
!
! calculate the gradient of rho+rhocore in real space
! in LSDA case rho is temporarily converted in (up,down) format
!
DO is = 1, nspin0
do is = 1, nspin0
!
rho(:,is) = fac * rho_core(:) + rho(:,is)
rhog(:,is) = fac * rhog_core(:) + rhog(:,is)
rho(:,is) = ( fac(is)* rho_core(:) + rho(:,1) + sgn(is)* rho(:,nspin0) )*is/2.d0
rhog(:,is) = ( fac(is)*rhog_core(:) + rhog(:,1) + sgn(is)*rhog(:,nspin0) )*is/2.d0
!
CALL fft_gradient_g2r( dfft, rhog(1,is), g, grho(1,1,is) )
call fft_gradient_g2r( dfft, rhog(1,is), g, grho(1,1,is) )
!
END DO
enddo
!
if (nspin.eq.1) then
!
@ -217,12 +219,13 @@ subroutine stres_gradcorr( rho, rhog, rho_core, rhog_core, kedtau, nspin, &
sigmaxc(:,:) = sigmaxc(:,:) + sigma_gradcorr(:,:) / &
(nr1 * nr2 * nr3)
DO is = 1, nspin0
! rho and rhog are restored to the original values
do is = 1, nspin0
!
rho(:,is) = rho(:,is) - fac * rho_core(:)
rhog(:,is) = rhog(:,is) - fac * rhog_core(:)
rho(:,is) =( rho(:,1) + is*sgn(is)* rho(:,nspin0) )*nspin0/2.d0 - mod(is,2)* rho_core(:)
rhog(:,is)=( rhog(:,1) + is*sgn(is)*rhog(:,nspin0) )*nspin0/2.d0 - mod(is,2)*rhog_core(:)
!
END DO
enddo
!
deallocate(grho)
return

2
PW/src/stres_nonloc_dft.f90
View File

@ -22,7 +22,7 @@ subroutine stres_nonloc_dft( rho, rho_core, nspin, sigma_nonloc_dft )
IMPLICIT NONE
!
integer, intent(in) ::nspin
real(DP), intent(in) :: rho (dfftp%nnr, nspin), rho_core (dfftp%nnr)
real(DP), intent(in) :: rho (dfftp%nnr), rho_core (dfftp%nnr)
real(DP), intent(inout) :: sigma_nonloc_dft (3, 3)
integer :: l, m, inlc

17
PW/src/stress.f90
View File

@ -20,7 +20,7 @@ subroutine stress ( sigma )
USE fft_base, ONLY : dfftp
USE ldaU, ONLY : lda_plus_u, U_projection
USE lsda_mod, ONLY : nspin
USE scf, ONLY : rho, rho_core, rhog_core, rhoz_or_updw
USE scf, ONLY : rho, rho_core, rhog_core
USE control_flags, ONLY : iverbosity, gamma_only, llondon, ldftd3, lxdm, ts_vdw
USE noncollin_module, ONLY : noncolin
USE funct, ONLY : dft_is_meta, dft_is_gradient
@ -76,8 +76,6 @@ subroutine stress ( sigma )
! hartree contribution
!
IF (.not.( do_comp_esm .and. ( esm_bc .ne. 'pbc' ) )) call stres_har(sigmahar)
!^
IF (nspin == 2) CALL rhoz_or_updw(rho, 'r_and_g', 'rhoz_updw')
!
! xc contribution (diagonal)
!
@ -95,10 +93,7 @@ subroutine stress ( sigma )
!
call stres_cc (sigmaxcc)
!
IF (nspin == 2) CALL rhoz_or_updw(rho, 'r_and_g', 'updw_rhoz')
!^
!
IF ( do_comp_esm .and. ( esm_bc .ne. 'pbc' ) ) THEN ! for ESM stress
IF ( do_comp_esm .and. ( esm_bc .ne. 'pbc' ) ) THEN ! for ESM stress
call esm_stres_loclong( sigmaloclong, rho%of_g(:,1) ) ! long range part
sigmaloc(:,:) = sigmaloc(:,:) + sigmaloclong(:,:)
END IF
@ -177,14 +172,8 @@ subroutine stress ( sigma )
!
! DFT-non_local contribution
!
!^
IF (nspin == 2) CALL rhoz_or_updw(rho, 'only_r', 'rhoz_updw')
!
sigma_nonloc_dft (:,:) = 0.d0
call stres_nonloc_dft(rho%of_r, rho_core, nspin, sigma_nonloc_dft)
!
IF (nspin == 2) CALL rhoz_or_updw(rho, 'only_r', 'updw_rhoz')
!^
call stres_nonloc_dft(rho%of_r(:,1), rho_core, nspin, sigma_nonloc_dft)
!
! SUM
!

12
PW/src/sum_band.f90
View File

@ -212,12 +212,12 @@ SUBROUTINE sum_band()
END DO
!
END IF
!^
IF (nspin == 2) THEN
CALL rhoz_or_updw( rho, 'r_and_g', 'updw_rhoz' )
IF (lda_plus_u.AND.(.NOT.noncolin)) CALL rhoz_or_updw( rho, 'hub_ns', 'updw_rhoz' )
ENDIF
!^
!
! ... if LSDA rho%of_r and rho%of_g are converted from (up,dw) to
! ... (up+dw,up-dw) format.
!
IF ( nspin == 2 ) CALL rhoz_or_updw( rho, 'r_and_g', '->rhoz' )
!
CALL stop_clock( 'sum_band' )
!
RETURN

285
PW/src/v_of_rho.f90
View File

@ -22,7 +22,7 @@ SUBROUTINE v_of_rho( rho, rho_core, rhog_core, &
USE ions_base, ONLY : nat, tau
USE ldaU, ONLY : lda_plus_U
USE funct, ONLY : dft_is_meta, get_meta
USE scf, ONLY : scf_type, rhoz_or_updw
USE scf, ONLY : scf_type
USE cell_base, ONLY : alat
USE lsda_mod, ONLY : nspin
USE control_flags, ONLY : ts_vdw
@ -52,9 +52,6 @@ SUBROUTINE v_of_rho( rho, rho_core, rhog_core, &
!
CALL start_clock( 'v_of_rho' )
!
!^
IF (nspin == 2) CALL rhoz_or_updw( rho, 'r_and_g', 'rhoz_updw' )
!
! ... calculate exchange-correlation potential
!
IF (dft_is_meta() .and. (get_meta() /= 4)) then
@ -63,20 +60,6 @@ SUBROUTINE v_of_rho( rho, rho_core, rhog_core, &
CALL v_xc( rho, rho_core, rhog_core, etxc, vtxc, v%of_r )
ENDIF
!
! ... add Tkatchenko-Scheffler potential (factor 2: Ha -> Ry)
!
IF (ts_vdw) THEN
CALL tsvdw_calculate(tau*alat,rho%of_r)
DO is = 1, nspin_lsda
DO ir=1,dfftp%nnr
v%of_r(ir,is)=v%of_r(ir,is)+2.0d0*UtsvdW(ir)
END DO
END DO
END IF
!
IF (nspin == 2) CALL rhoz_or_updw( rho, 'r_and_g', 'updw_rhoz' )
!^
!
! ... add a magnetic field (if any)
!
CALL add_bfield( v%of_r, rho%of_r )
@ -101,6 +84,17 @@ SUBROUTINE v_of_rho( rho, rho_core, rhog_core, &
CALL add_efield(v%of_r(1,is), etotefield, rho%of_r(:,1), .false. )
END DO
!
! ... add Tkatchenko-Scheffler potential (factor 2: Ha -> Ry)
!
IF (ts_vdw) THEN
CALL tsvdw_calculate(tau*alat,rho%of_r(:,1),nspin)
DO is = 1, nspin_lsda
DO ir=1,dfftp%nnr
v%of_r(ir,is)=v%of_r(ir,is)+2.0d0*UtsvdW(ir)
END DO
END DO
END IF
!
CALL stop_clock( 'v_of_rho' )
!
RETURN
@ -141,7 +135,7 @@ SUBROUTINE v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v, kedtaur )
!
! ... local variables
!
REAL(DP) :: zeta, rh
REAL(DP) :: zeta, rh, sgn(2)
INTEGER :: k, ipol, is
REAL(DP) :: ex, ec, v1x, v2x, v3x,v1c, v2c, v3c, &
& v1xup, v1xdw, v2xup, v2xdw, v1cup, v1cdw, &
@ -167,7 +161,8 @@ SUBROUTINE v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v, kedtaur )
vtxc = zero
v(:,:) = zero
rhoneg(:) = zero
!
sgn(1) = 1._dp ; sgn(2) = -1._dp
fac = 1.D0 / DBLE( nspin )
!
ALLOCATE (grho(3,dfftp%nnr,nspin))
ALLOCATE (h(3,dfftp%nnr,nspin))
@ -175,15 +170,12 @@ SUBROUTINE v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v, kedtaur )
ALLOCATE (rhogsum(ngm,nspin))
!
! ... calculate the gradient of rho + rho_core in real space
!
rhoout(:,1:nspin)=rho%of_r(:,1:nspin)
rhogsum(:,1:nspin)=rho%of_g(:,1:nspin)
fac = 1.D0 / DBLE( nspin )
! ... in LSDA case rhoout and rhogsum are defined in (up,down) format
!
DO is = 1, nspin
!
rhoout(:,is) = fac * rho_core(:) + rhoout(:,is)
rhogsum(:,is) = fac * rhog_core(:) + rhogsum(:,is)
rhoout(:,is) =fac*rho_core(:) + ( rho%of_r(:,1) + sgn(is)*rho%of_r(:,nspin) )*0.5D0
rhogsum(:,is)=fac*rhog_core(:) + ( rho%of_g(:,1) + sgn(is)*rho%of_g(:,nspin) )*0.5D0
!
CALL fft_gradient_g2r( dfftp, rhogsum(1,is), g, grho(1,1,is) )
!
@ -229,8 +221,8 @@ SUBROUTINE v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v, kedtaur )
!
! spin-polarised case
!
rhoup=rho%of_r(k, 1)
rhodw=rho%of_r(k, 2)
rhoup = ( rho%of_r(k, 1) + rho%of_r(k, 2) )*0.5d0
rhodw = ( rho%of_r(k, 2) - rho%of_r(k, 2) )*0.5d0
rh = rhoup + rhodw
@ -285,8 +277,8 @@ SUBROUTINE v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v, kedtaur )
end if
if (rho%of_r (k, 1) < zero ) rhoneg(1) = rhoneg(1) - rho%of_r (k, 1)
if (rho%of_r (k, 2) < zero ) rhoneg(2) = rhoneg(2) - rho%of_r (k, 2)
if (rhoup < zero ) rhoneg(1) = rhoneg(1) - rhoup
if (rhodw < zero ) rhoneg(2) = rhoneg(2) - rhodw
end if
end do
@ -357,7 +349,7 @@ SUBROUTINE v_xc( rho, rho_core, rhog_core, etxc, vtxc, v )
!
IMPLICIT NONE
!
TYPE (scf_type), INTENT(IN) :: rho
TYPE (scf_type), INTENT(INOUT) :: rho
REAL(DP), INTENT(IN) :: rho_core(dfftp%nnr)
! the core charge
COMPLEX(DP), INTENT(IN) :: rhog_core(ngm)
@ -429,18 +421,23 @@ SUBROUTINE v_xc( rho, rho_core, rhog_core, etxc, vtxc, v )
!$omp reduction(+:etxc,vtxc), reduction(-:rhoneg)
DO ir = 1, dfftp%nnr
!
rhox = rho%of_r(ir,1) + rho%of_r(ir,2) + rho_core(ir)
rhox = rho%of_r(ir,1) + rho_core(ir)
!
IF ( rho%of_r(ir,1) < 0.D0 ) rhoneg(1) = rhoneg(1) - rho%of_r(ir,1)
!
arhox = ABS( rhox )
!
IF ( arhox > vanishing_charge ) THEN
!
zeta = ( rho%of_r(ir,1) - rho%of_r(ir,2) ) / arhox
zeta = rho%of_r(ir,2) / arhox
!
IF ( ABS( zeta ) > 1.D0 ) zeta = SIGN( 1.D0, zeta )
!
IF ( rho%of_r(ir,1) < 0.D0 ) rhoneg(1) = rhoneg(1) - rho%of_r(ir,1)
IF ( rho%of_r(ir,2) < 0.D0 ) rhoneg(2) = rhoneg(2) - rho%of_r(ir,2)
IF ( ABS( zeta ) > 1.D0 ) THEN
!
rhoneg(2) = rhoneg(2) + 1.D0 / omega
!
zeta = SIGN( 1.D0, zeta )
!
END IF
!
CALL xc_spin( arhox, zeta, ex, ec, vx(1), vx(2), vc(1), vc(2) )
!
@ -448,13 +445,16 @@ SUBROUTINE v_xc( rho, rho_core, rhog_core, etxc, vtxc, v )
!
etxc = etxc + e2*( ex + ec ) * rhox
!
vtxc = vtxc + ( v(ir,1)*rho%of_r(ir,1) + v(ir,2)*rho%of_r(ir,2) )
vtxc = vtxc + ( ( v(ir,1) + v(ir,2) )*rho%of_r(ir,1) + &
( v(ir,1) - v(ir,2) )*rho%of_r(ir,2) )
!
END IF
!
END DO
!$omp end parallel do
!
vtxc = 0.5d0 * vtxc
!
ELSE IF ( nspin == 4 ) THEN
!
! ... noncolinear case
@ -699,85 +699,67 @@ SUBROUTINE v_hubbard(ns, v_hub, eth)
IMPLICIT NONE
!
REAL(DP), INTENT(INOUT) :: ns(2*Hubbard_lmax+1,2*Hubbard_lmax+1,nspin,nat)
REAL(DP), INTENT(IN) :: ns(2*Hubbard_lmax+1,2*Hubbard_lmax+1,nspin,nat)
REAL(DP), INTENT(OUT) :: v_hub(2*Hubbard_lmax+1,2*Hubbard_lmax+1,nspin,nat)
REAL(DP), INTENT(OUT) :: eth
REAL(DP) :: n_tot, n_spin, eth_dc, eth_u, mag2, effU, sgn(2)
REAL(DP) :: n_tot, n_spin, eth_dc, eth_u, mag2, effU, sgn(2)
INTEGER :: is, isop, is1, na, nt, m1, m2, m3, m4
REAL(DP), ALLOCATABLE :: u_matrix(:,:,:,:)
ALLOCATE( u_matrix(2*Hubbard_lmax+1, 2*Hubbard_lmax+1, 2*Hubbard_lmax+1, 2*Hubbard_lmax+1) )
!
eth = 0.d0
eth_dc = 0.d0
eth_u = 0.d0
sgn(1) = 1.d0 ; sgn(2) = -1.d0
!
sgn(1)=1.d0 ; sgn(2)=-1.d0
v_hub(:,:,:,:) = 0.d0
!
IF (lda_plus_u_kind == 0) THEN
!
if (lda_plus_u_kind.eq.0) then
DO na = 1, nat
nt = ityp (na)
IF (Hubbard_U(nt) /= 0.d0 .OR. Hubbard_alpha(nt) /= 0.d0) THEN
!
IF (Hubbard_J0(nt) /= 0.d0) THEN
IF (Hubbard_U(nt).NE.0.d0 .OR. Hubbard_alpha(nt).NE.0.d0) THEN
IF (Hubbard_J0(nt).NE.0.d0) THEN
effU = Hubbard_U(nt) - Hubbard_J0(nt)
ELSE
effU = Hubbard_U(nt)
END IF
!
END IF
DO is = 1, nspin
!
DO m1 = 1, 2 * Hubbard_l(nt) + 1
!
eth = eth + ( Hubbard_alpha(nt) + 0.5D0 * effU ) * &
( ns(m1,m1,1,na) + sgn(is)*ns(m1,m1,nspin,na) )*0.5D0
eth = eth + ( Hubbard_alpha(nt) + 0.5D0*effU )*ns(m1,m1,is,na)
v_hub(m1,m1,is,na) = v_hub(m1,m1,is,na) + &
( Hubbard_alpha(nt) + 0.5D0 * effU )
!
Hubbard_alpha(nt) + 0.5D0*effU
DO m2 = 1, 2 * Hubbard_l(nt) + 1
!
eth = eth - 0.125D0*effU*( ns(m2,m1,1,na) + sgn(is)*ns(m2,m1,nspin,na) ) * &
( ns(m1,m2,1,na) + sgn(is)*ns(m1,m2,nspin,na) )
!
v_hub(m1,m2,is,na) = v_hub(m1,m2,is,na) - effU * &
( ns(m2,m1,1,na) + sgn(is)*ns(m2,m1,nspin,na) )*0.5D0
eth = eth - 0.5D0 * effU * ns(m2,m1,is,na)* ns(m1,m2,is,na)
v_hub(m1,m2,is,na) = v_hub(m1,m2,is,na) - &
effU * ns(m2,m1,is,na)
ENDDO
ENDDO
ENDDO
ENDIF
!
IF (Hubbard_J0(nt) /= 0.d0 .OR. Hubbard_beta(nt) /= 0.d0) THEN
DO is=1, nspin
!
isop = 1
IF ( nspin == 2 .AND. is == 1) isop = 2
!
DO m1 = 1, 2 * Hubbard_l(nt) + 1
!
eth = eth + sgn(is) * Hubbard_beta(nt) * &
( ns(m1,m1,1,na) + sgn(is)*ns(m1,m1,nspin,na) ) * 0.5D0
!
eth = eth + sgn(is)*Hubbard_beta(nt) * ns(m1,m1,is,na)
v_hub(m1,m1,is,na) = v_hub(m1,m1,is,na) + sgn(is)*Hubbard_beta(nt)
!
DO m2 = 1, 2 * Hubbard_l(nt) + 1
!
eth = eth + 0.125D0*Hubbard_J0(nt) * ( ns(m2,m1,1,na) + sgn(is)*ns(m2,m1,nspin,na) )* &
( ns(m1,m2,1,na) + sgn(isop)*ns(m1,m2,nspin,na) )
!
DO m2 = 1, 2*Hubbard_l(nt)+1
eth = eth + 0.5D0*Hubbard_J0(nt)*ns(m2,m1,is,na)*ns(m1,m2,isop,na)
v_hub(m1,m2,is,na) = v_hub(m1,m2,is,na) + Hubbard_J0(nt) * &
( ns(m2,m1,1,na) + sgn(isop)*ns(m2,m1,nspin,na) ) * 0.5D0
ns(m2,m1,isop,na)
END DO
END DO
END DO
END IF
!
END DO
!
IF (nspin == 1) eth = 2.d0 * eth
!
IF (nspin.EQ.1) eth = 2.d0 * eth
!-- output of hubbard energies:
IF ( iverbosity > 0 ) THEN
write(stdout,*) '--- in v_hubbard ---'
@ -785,91 +767,89 @@ SUBROUTINE v_hubbard(ns, v_hub, eth)
write(stdout,*) '-------'
ENDIF
!--
!
ELSE
!
else
DO na = 1, nat
nt = ityp (na)
IF (Hubbard_U(nt) /= 0.d0) THEN
!
! initialize U(m1,m2,m3,m4) matrix
CALL hubbard_matrix (Hubbard_lmax, Hubbard_l(nt), Hubbard_U(nt), &
IF (Hubbard_U(nt).NE.0.d0) THEN
! initialize U(m1,m2,m3,m4) matrix
call hubbard_matrix (Hubbard_lmax, Hubbard_l(nt), Hubbard_U(nt), &
Hubbard_J(1,nt), u_matrix)
!
!--- total N and M^2 for DC (double counting) term
n_tot = 0.d0
!
DO m1 = 1, 2 * Hubbard_l(nt) + 1
n_tot = n_tot + ns(m1,m1,1,na)
ENDDO
IF ( nspin == 1 ) n_tot = 2.d0 * n_tot
!
do is = 1, nspin
do m1 = 1, 2 * Hubbard_l(nt) + 1
n_tot = n_tot + ns(m1,m1,is,na)
enddo
enddo
if (nspin.eq.1) n_tot = 2.d0 * n_tot
mag2 = 0.d0
IF ( nspin == 2 ) THEN
DO m1 = 1, 2 * Hubbard_l(nt) + 1
mag2 = mag2 + ns(m1,m1,2,na)
ENDDO
ENDIF
if (nspin.eq.2) then
do m1 = 1, 2 * Hubbard_l(nt) + 1
mag2 = mag2 + ns(m1,m1,1,na) - ns(m1,m1,2,na)
enddo
endif
mag2 = mag2**2
!--- !
!
!---
!--- hubbard energy: DC term
!
eth_dc = eth_dc + 0.5d0*( Hubbard_U(nt)*n_tot*(n_tot-1.d0) - &
eth_dc = eth_dc + 0.5d0*( Hubbard_U(nt)*n_tot*(n_tot-1.d0) - &
Hubbard_J(1,nt)*n_tot*(0.5d0*n_tot-1.d0) - &
0.5d0*Hubbard_J(1,nt)*mag2 )
!-- !
!--
DO is = 1, nspin
!
!--- n_spin = up/down N
!
n_spin = 0.d0
do m1 = 1, 2 * Hubbard_l(nt) + 1
n_spin = n_spin + ( ns(m1,m1,1,na) + sgn(is)*ns(m1,m1,nspin,na) ) * 0.5D0
n_spin = n_spin + ns(m1,m1,is,na)
enddo
!---
!
DO m1 = 1, 2 * Hubbard_l(nt) + 1
!
! hubbard potential: DC contribution
!
v_hub(m1,m1,is,na) = v_hub(m1,m1,is,na) + Hubbard_J(1,nt)*n_spin + &
0.5d0*(Hubbard_U(nt)-Hubbard_J(1,nt)) - Hubbard_U(nt)*n_tot
!
! +U contributions
!
DO m2 = 1, 2 * Hubbard_l(nt) + 1
DO m3 = 1, 2 * Hubbard_l(nt) + 1
DO m4 = 1, 2 * Hubbard_l(nt) + 1
!
DO is1 = 1, nspin
v_hub(m1,m2,is,na) = v_hub(m1,m2,is,na) &
+ REAL( MOD(nspin,2)+1 ) * u_matrix(m1,m3,m2,m4) * &
( ns(m3,m4,1,na) + sgn(is1)*ns(m3,m4,nspin,na) ) * 0.5D0
ENDDO
!
v_hub(m1,m2,is,na) = v_hub(m1,m2,is,na) - u_matrix(m1,m3,m4,m2) * &
( ns(m3,m4,1,na) + sgn(is)*ns(m3,m4,nspin,na) ) * 0.5D0
!
eth_u = eth_u + 0.125d0*( &
( u_matrix(m1,m2,m3,m4) - u_matrix(m1,m2,m4,m3) ) * &
( ns(m1,m3,1,na) + sgn(is)*ns(m1,m3,nspin,na) ) * &
( ns(m2,m4,1,na) + sgn(is)*ns(m2,m4,nspin,na) ) + u_matrix(m1,m2,m3,m4) * &
( ns(m1,m3,1,na) + sgn(is)*ns(m1,m3,nspin,na) ) * &
( ns(m2,m4,1,na) + sgn(nspin+1-is)*ns(m2,m4,nspin,na) ) )
ENDDO
ENDDO
do m3 = 1, 2 * Hubbard_l(nt) + 1
do m4 = 1, 2 * Hubbard_l(nt) + 1
do is1 = 1, nspin
v_hub(m1,m2,is,na) = v_hub(m1,m2,is,na) + (MOD(nspin,2)+1) * &
u_matrix(m1,m3,m2,m4) * ns(m3,m4,is1,na)
enddo
v_hub(m1,m2,is,na) = v_hub(m1,m2,is,na) - &
u_matrix(m1,m3,m4,m2) * ns(m3,m4,is,na)
eth_u = eth_u + 0.5d0*( ( u_matrix(m1,m2,m3,m4)-u_matrix(m1,m2,m4,m3) ) * &
ns(m1,m3,is,na)*ns(m2,m4,is,na)+u_matrix(m1,m2,m3,m4) * &
ns(m1,m3,is,na)*ns(m2,m4,nspin+1-is,na) )
enddo
enddo
ENDDO
ENDDO
ENDDO
ENDIF
ENDDO
!
IF ( nspin == 1 ) eth_u = 2.d0 * eth_u
endif
enddo
if (nspin.eq.1) eth_u = 2.d0 * eth_u
eth = eth_u - eth_dc
!
!-- output of hubbard energies:
IF ( iverbosity > 0 ) THEN
write(stdout,*) '--- in v_hubbard ---'
@ -877,12 +857,12 @@ SUBROUTINE v_hubbard(ns, v_hub, eth)
write(stdout,*) '-------'
ENDIF
!--
!
ENDIF
!
endif
DEALLOCATE (u_matrix)
RETURN
!
END SUBROUTINE v_hubbard
!-------------------------------------
@ -1098,8 +1078,23 @@ SUBROUTINE v_h_of_rho_r( rhor, ehart, charge, v )
!
! ... compute VH(r) from n(G)
!
CALL v_h( rhog(:,1), ehart, charge, v ) !^ ---wrong input rhog (CPV). Will be adjusted
DEALLOCATE( rhog ) !^ ---before the next merge.
!^^ ... TEMPORARY FIX (newlsda-CPV) ...
IF ( nspin==2 ) THEN
rhog(:,1) = rhog(:,1) + rhog(:,2)
rhog(:,2) = rhog(:,1) - rhog(:,2)*2._dp
ENDIF
!^^.......................
!
CALL v_h( rhog(:,1), ehart, charge, v )
!
!^^ ... TEMPORARY FIX (newlsda) ...
IF ( nspin==2 ) THEN
rhog(:,1) = ( rhog(:,1) + rhog(:,2) )*0.5_dp
rhog(:,2) = rhog(:,1) - rhog(:,2)
ENDIF
!^^.......................
!
DEALLOCATE( rhog )
!
RETURN
!

6
PW/src/write_ns.f90
View File

@ -83,8 +83,7 @@ subroutine write_ns
nsuma = 0.d0
do is = 1, nspin
do m1 = 1, ldim
nsuma(is) = nsuma(is) + ( rho%ns (m1, m1, 1, na) + DBLE(2*MOD(is,2)-1) * &
rho%ns (m1, m1, nspin, na) )*0.5d0
nsuma(is) = nsuma(is) + rho%ns (m1, m1, is, na)
enddo
nsum = nsum + nsuma(is)
enddo
@ -100,8 +99,7 @@ subroutine write_ns
do is = 1, nspin
do m1 = 1, ldim
do m2 = 1, ldim
f (m1, m2) = ( rho%ns(m1,m2,1,na) + DBLE(2*MOD(is,2)-1) * &
rho%ns(m1,m2,nspin,na) )*0.5d0
f (m1, m2) = rho%ns (m1, m2, is, na)
enddo
enddo
call cdiagh(ldim, f, ldmx, lambda, vet)

9
PW/src/xdm_dispersion.f90
View File

@ -184,7 +184,7 @@ CONTAINS
! runs. In addition, forces and stresses are saved for subsequent calls to force_xdm
! and stress_xdm.
USE control_flags, ONLY: lbfgs, lmd
USE scf, ONLY: rho
USE scf, ONLY: rho, rhoz_or_updw
USE io_global, ONLY: stdout, ionode
USE fft_base, ONLY : dfftp
USE cell_base, ONLY : at, alat, omega
@ -230,6 +230,9 @@ CONTAINS
fsave = 0._DP
ssave = 0._DP
atb = at * alat
!
! for convenience rho is converted in (up,down) format, if LSDA
IF (nspin == 2) CALL rhoz_or_updw( rho, 'r_and_g', '->updw' )
! do we need to recalculate the coefficients?
docalc = .NOT.saved .OR. .NOT.(lbfgs .OR. lmd)
@ -559,7 +562,9 @@ CONTAINS
DO nn = 6, 10
CALL mp_sum(ehadd(nn),intra_image_comm)
ENDDO
!
IF (nspin == 2) CALL rhoz_or_updw( rho, 'r_and_g', '->rhoz' )
!
! Convert to Ry
evdw = evdw * 2
for = for * 2

12
PWCOND/src/plus_u_setup.f90
View File

@ -27,7 +27,7 @@ subroutine plus_u_setup(natih, lsr)
natih(2,norbs), lll, kkbeta
integer, allocatable :: ind(:,:), tblms_new(:,:), cross_new(:,:)
real(DP), parameter :: epswfc=1.d-4, eps=1.d-8
REAL(DP) :: r1, beta1, beta2, norm, ledge, redge, sgn
REAL(DP) :: r1, beta1, beta2, norm, ledge, redge
REAL(DP), ALLOCATABLE :: bphi(:,:), rsph(:), taunews_new(:,:), &
gi(:), zpseus_new(:,:,:)
@ -65,7 +65,6 @@ subroutine plus_u_setup(natih, lsr)
bphi(:,:) = 0.d0
rsph(:) = 0.d0
ind(:,:) = 0
sgn = dble(2*mod(iofspin,2)-1)
!--
! Calculate the total number of orbitals (beta + U WF's)
@ -235,8 +234,7 @@ subroutine plus_u_setup(natih, lsr)
if (tblms(3,iwfc).eq.lll) then
do iwfc1 = ind(1,iorb), iorb
if (tblms(3,iwfc1).eq.lll) then
r1 = -( rho%ns(tblms(4,iwfc),tblms(4,iwfc1),1,na) + sgn* &
rho%ns(tblms(4,iwfc),tblms(4,iwfc1),nspin,na) )
r1 = -2.d0*rho%ns(tblms(4,iwfc),tblms(4,iwfc1),iofspin,na)
if (tblms(4,iwfc).eq.tblms(4,iwfc1)) r1 = r1 + 1.d0
zpseus_new(1,iwfc-iorb+iorb1,iwfc1-iorb+iorb1) = &
zpseus_new(1,iwfc-iorb+iorb1,iwfc1-iorb+iorb1) + &
@ -254,8 +252,7 @@ subroutine plus_u_setup(natih, lsr)
do iwfc = ind(1,iorb), iorb
if (tblms(3,iwfc).eq.lll) then
do iwfc1 = 1, ldim
r1 = -( rho%ns(tblms(4,iwfc),iwfc1,1,na) + sgn* &
rho%ns(tblms(4,iwfc),iwfc1,nspin,na) )
r1 = -2.d0*rho%ns(tblms(4,iwfc),iwfc1,iofspin,na)
if (tblms(4,iwfc).eq.iwfc1) r1 = r1 + 1.d0
zpseus_new(1,iwfc-iorb+iorb1,iorb1+iwfc1) = &
0.5d0*Hubbard_U(it)*bphi(tblms(2,iwfc),it)*r1
@ -281,8 +278,7 @@ subroutine plus_u_setup(natih, lsr)
do iwfc = 1, ldim
do iwfc1 = 1, ldim
zpseus_new(1,iorb1+iwfc,iorb1+iwfc1) = &
- Hubbard_U(it) * ( rho%ns(iwfc,iwfc1,1,na) +sgn* &
rho%ns(iwfc,iwfc1,nspin,na) )*0.5d0
- Hubbard_U(it) * rho%ns(iwfc,iwfc1,iofspin,na)
enddo
zpseus_new(1,iorb1+iwfc,iorb1+iwfc) = &
zpseus_new(1,iorb1+iwfc,iorb1+iwfc) + 0.5d0*Hubbard_U(it)