diff --git a/CPV/src/chargedensity.f90 b/CPV/src/chargedensity.f90 index 58ed48ada..0eebc70c1 100644 --- a/CPV/src/chargedensity.f90 +++ b/CPV/src/chargedensity.f90 @@ -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. diff --git a/CPV/src/cp_restart_new.f90 b/CPV/src/cp_restart_new.f90 index 4a1d5905c..10f3166e0 100644 --- a/CPV/src/cp_restart_new.f90 +++ b/CPV/src/cp_restart_new.f90 @@ -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 diff --git a/CPV/src/vofrho.f90 b/CPV/src/vofrho.f90 index 9fa81cb8f..f41aede05 100644 --- a/CPV/src/vofrho.f90 +++ b/CPV/src/vofrho.f90 @@ -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 ) diff --git a/EPW/src/epw_setup.f90 b/EPW/src/epw_setup.f90 index 4106e4c12..2c2d69717 100644 --- a/EPW/src/epw_setup.f90 +++ b/EPW/src/epw_setup.f90 @@ -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() diff --git a/GWW/pw4gww/energies_xc.f90 b/GWW/pw4gww/energies_xc.f90 index b3daeb8ae..c83cbb24c 100644 --- a/GWW/pw4gww/energies_xc.f90 +++ b/GWW/pw4gww/energies_xc.f90 @@ -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) diff --git a/HP/src/hp_setup_q.f90 b/HP/src/hp_setup_q.f90 index 3c831d783..17345d1bc 100644 --- a/HP/src/hp_setup_q.f90 +++ b/HP/src/hp_setup_q.f90 @@ -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() diff --git a/LR_Modules/commutator_Vhubx_psi.f90 b/LR_Modules/commutator_Vhubx_psi.f90 index ac8484d89..651fd5974 100644 --- a/LR_Modules/commutator_Vhubx_psi.f90 +++ b/LR_Modules/commutator_Vhubx_psi.f90 @@ -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 diff --git a/LR_Modules/dgradcorr.f90 b/LR_Modules/dgradcorr.f90 index 15cc88a10..300bb9a59 100644 --- a/LR_Modules/dgradcorr.f90 +++ b/LR_Modules/dgradcorr.f90 @@ -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 diff --git a/LR_Modules/dv_of_drho.f90 b/LR_Modules/dv_of_drho.f90 index b690492e4..397c422f0 100644 --- a/LR_Modules/dv_of_drho.f90 +++ b/LR_Modules/dv_of_drho.f90 @@ -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 (:) diff --git a/LR_Modules/setup_dgc.f90 b/LR_Modules/setup_dgc.f90 index 8ce328b68..0902b2a45 100644 --- a/LR_Modules/setup_dgc.f90 +++ b/LR_Modules/setup_dgc.f90 @@ -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) diff --git a/Modules/funct.f90 b/Modules/funct.f90 index fc3876d85..a60c85ff0 100644 --- a/Modules/funct.f90 +++ b/Modules/funct.f90 @@ -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 diff --git a/Modules/tsvdw.f90 b/Modules/tsvdw.f90 index 83765585b..26c808704 100644 --- a/Modules/tsvdw.f90 +++ b/Modules/tsvdw.f90 @@ -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') ! diff --git a/Modules/xc_rVV10.f90 b/Modules/xc_rVV10.f90 index 36367f0c6..e35e02ffc 100644 --- a/Modules/xc_rVV10.f90 +++ b/Modules/xc_rVV10.f90 @@ -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 diff --git a/Modules/xc_vdW_DF.f90 b/Modules/xc_vdW_DF.f90 index 752ae9627..33e2c290c 100644 --- a/Modules/xc_vdW_DF.f90 +++ b/Modules/xc_vdW_DF.f90 @@ -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 diff --git a/PHonon/Gamma/cg_setupdgc.f90 b/PHonon/Gamma/cg_setupdgc.f90 index 5c0ca6275..d066ca29a 100644 --- a/PHonon/Gamma/cg_setupdgc.f90 +++ b/PHonon/Gamma/cg_setupdgc.f90 @@ -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') ! diff --git a/PHonon/Gamma/dynmatcc.f90 b/PHonon/Gamma/dynmatcc.f90 index 4eecebc50..e694364d0 100644 --- a/PHonon/Gamma/dynmatcc.f90 +++ b/PHonon/Gamma/dynmatcc.f90 @@ -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 diff --git a/PHonon/PH/addnlcc.f90 b/PHonon/PH/addnlcc.f90 index e55800ddd..6ee7ddea7 100644 --- a/PHonon/PH/addnlcc.f90 +++ b/PHonon/PH/addnlcc.f90 @@ -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 diff --git a/PHonon/PH/addnlcc_zstar_eu_us.f90 b/PHonon/PH/addnlcc_zstar_eu_us.f90 index d43152ea0..27a49d4e5 100644 --- a/PHonon/PH/addnlcc_zstar_eu_us.f90 +++ b/PHonon/PH/addnlcc_zstar_eu_us.f90 @@ -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) * & diff --git a/PHonon/PH/dvqhub_barepsi_us.f90 b/PHonon/PH/dvqhub_barepsi_us.f90 index d9e7417b2..a1e9d7307 100644 --- a/PHonon/PH/dvqhub_barepsi_us.f90 +++ b/PHonon/PH/dvqhub_barepsi_us.f90 @@ -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) & diff --git a/PHonon/PH/dvqhub_barepsi_us2.f90 b/PHonon/PH/dvqhub_barepsi_us2.f90 index 9b0078255..8d289f6f1 100644 --- a/PHonon/PH/dvqhub_barepsi_us2.f90 +++ b/PHonon/PH/dvqhub_barepsi_us2.f90 @@ -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) ) ! diff --git a/PHonon/PH/dvqpsi_us.f90 b/PHonon/PH/dvqpsi_us.f90 index 10df5e672..0c3cd1312 100644 --- a/PHonon/PH/dvqpsi_us.f90 +++ b/PHonon/PH/dvqpsi_us.f90 @@ -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 diff --git a/PHonon/PH/dynmat_hub_bare.f90 b/PHonon/PH/dynmat_hub_bare.f90 index 910434acc..5040cafc6 100644 --- a/PHonon/PH/dynmat_hub_bare.f90 +++ b/PHonon/PH/dynmat_hub_bare.f90 @@ -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 - diff --git a/PHonon/PH/dynmatcc.f90 b/PHonon/PH/dynmatcc.f90 index 38cfb5e30..47ce3e2aa 100644 --- a/PHonon/PH/dynmatcc.f90 +++ b/PHonon/PH/dynmatcc.f90 @@ -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 diff --git a/PHonon/PH/phq_setup.f90 b/PHonon/PH/phq_setup.f90 index 405970169..866c3400d 100644 --- a/PHonon/PH/phq_setup.f90 +++ b/PHonon/PH/phq_setup.f90 @@ -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() diff --git a/PP/src/add_shift_cc.f90 b/PP/src/add_shift_cc.f90 index 80352d48f..203fbb7fc 100644 --- a/PP/src/add_shift_cc.f90 +++ b/PP/src/add_shift_cc.f90 @@ -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) diff --git a/PP/src/local_dos.f90 b/PP/src/local_dos.f90 index d03e98496..3702a70e6 100644 --- a/PP/src/local_dos.f90 +++ b/PP/src/local_dos.f90 @@ -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) diff --git a/PP/src/paw_postproc.f90 b/PP/src/paw_postproc.f90 index 2f047da34..b9f9218f6 100644 --- a/PP/src/paw_postproc.f90 +++ b/PP/src/paw_postproc.f90 @@ -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) diff --git a/PP/src/ppacf.f90 b/PP/src/ppacf.f90 index 43a8ce48f..8cbef7ff9 100644 --- a/PP/src/ppacf.f90 +++ b/PP/src/ppacf.f90 @@ -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) diff --git a/PP/src/punch_plot.f90 b/PP/src/punch_plot.f90 index 48241a27d..7bf7e177d 100644 --- a/PP/src/punch_plot.f90 +++ b/PP/src/punch_plot.f90 @@ -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 diff --git a/PP/src/pw2bgw.f90 b/PP/src/pw2bgw.f90 index 64712ec11..e12724aa7 100644 --- a/PP/src/pw2bgw.f90 +++ b/PP/src/pw2bgw.f90 @@ -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 ) diff --git a/PP/src/pw2gw.f90 b/PP/src/pw2gw.f90 index a3466d557..56435c9a2 100644 --- a/PP/src/pw2gw.f90 +++ b/PP/src/pw2gw.f90 @@ -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 " 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 ) diff --git a/PW/src/addusdens.f90 b/PW/src/addusdens.f90 index 30b16defa..c53a83eba 100644 --- a/PW/src/addusdens.f90 +++ b/PW/src/addusdens.f90 @@ -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> 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) ! diff --git a/PW/src/atomic_rho.f90 b/PW/src/atomic_rho.f90 index a7cc46081..9d5f22b13 100644 --- a/PW/src/atomic_rho.f90 +++ b/PW/src/atomic_rho.f90 @@ -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 ! diff --git a/PW/src/electrons.f90 b/PW/src/electrons.f90 index d609c9735..5bf91ba79 100644 --- a/PW/src/electrons.f90 +++ b/PW/src/electrons.f90 @@ -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 ! diff --git a/PW/src/force_cc.f90 b/PW/src/force_cc.f90 index 05767e596..893fa5d02 100644 --- a/PW/src/force_cc.f90 +++ b/PW/src/force_cc.f90 @@ -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 diff --git a/PW/src/gradcorr.f90 b/PW/src/gradcorr.f90 index 8a3136786..09f1f88ee 100644 --- a/PW/src/gradcorr.f90 +++ b/PW/src/gradcorr.f90 @@ -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) ) diff --git a/PW/src/init_ns.f90 b/PW/src/init_ns.f90 index e915d29c1..5c3993bb2 100644 --- a/PW/src/init_ns.f90 +++ b/PW/src/init_ns.f90 @@ -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 diff --git a/PW/src/ns_adj.f90 b/PW/src/ns_adj.f90 index 6ec5a2285..360e6ec28 100644 --- a/PW/src/ns_adj.f90 +++ b/PW/src/ns_adj.f90 @@ -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 diff --git a/PW/src/realus.f90 b/PW/src/realus.f90 index 5799dbf17..17a285e80 100644 --- a/PW/src/realus.f90 +++ b/PW/src/realus.f90 @@ -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 diff --git a/PW/src/scf_mod.f90 b/PW/src/scf_mod.f90 index 0d689fb83..fec0a37b1 100644 --- a/PW/src/scf_mod.f90 +++ b/PW/src/scf_mod.f90 @@ -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 ! diff --git a/PW/src/stres_gradcorr.f90 b/PW/src/stres_gradcorr.f90 index 90200fc5c..103887eca 100644 --- a/PW/src/stres_gradcorr.f90 +++ b/PW/src/stres_gradcorr.f90 @@ -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 diff --git a/PW/src/stres_nonloc_dft.f90 b/PW/src/stres_nonloc_dft.f90 index b8920db64..3812c1cda 100644 --- a/PW/src/stres_nonloc_dft.f90 +++ b/PW/src/stres_nonloc_dft.f90 @@ -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 diff --git a/PW/src/stress.f90 b/PW/src/stress.f90 index 280adfe9b..201a3ac27 100644 --- a/PW/src/stress.f90 +++ b/PW/src/stress.f90 @@ -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 ! diff --git a/PW/src/sum_band.f90 b/PW/src/sum_band.f90 index a86eeae7b..5d6c3632a 100644 --- a/PW/src/sum_band.f90 +++ b/PW/src/sum_band.f90 @@ -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 diff --git a/PW/src/v_of_rho.f90 b/PW/src/v_of_rho.f90 index 825938bc4..7081a7f31 100644 --- a/PW/src/v_of_rho.f90 +++ b/PW/src/v_of_rho.f90 @@ -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 ! diff --git a/PW/src/write_ns.f90 b/PW/src/write_ns.f90 index 7e8190141..0ef12e344 100644 --- a/PW/src/write_ns.f90 +++ b/PW/src/write_ns.f90 @@ -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) diff --git a/PW/src/xdm_dispersion.f90 b/PW/src/xdm_dispersion.f90 index e9b62ba10..bab6a108f 100644 --- a/PW/src/xdm_dispersion.f90 +++ b/PW/src/xdm_dispersion.f90 @@ -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 diff --git a/PWCOND/src/plus_u_setup.f90 b/PWCOND/src/plus_u_setup.f90 index a66cc416d..999f6ac93 100644 --- a/PWCOND/src/plus_u_setup.f90 +++ b/PWCOND/src/plus_u_setup.f90 @@ -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)