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

also the scf potential is defined using the scf_type. 

then call to v_of_r has been simplified


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@4442 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
degironc 2007-11-18 20:25:11 +00:00
parent bf74ac09d1
commit ae8be06438
34 changed files with 159 additions and 174 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
D3/d3_setup.f90
View File

@ -47,7 +47,7 @@ SUBROUTINE d3_setup()
USE io_global, ONLY : stdout
USE kinds, ONLY : DP
USE pwcom
USE scf, only : rho, rho_core, vr, vltot, vrs
USE scf, only : rho, rho_core, v, vltot, vrs, kedtau
USE uspp_param, ONLY : upf
USE control_flags, ONLY : iverbosity, modenum
USE constants, ONLY : degspin
@ -84,7 +84,7 @@ SUBROUTINE d3_setup()
!
! 1) Computes the total local potential (external+scf) on the smoot grid
!
CALL set_vrs (vrs, vltot, vr, nrxx, nspin, doublegrid)
CALL set_vrs (vrs, vltot, v%of_r, kedtau, v%kin_r, nrxx, nspin, doublegrid)
!
! 2) Computes the derivative of the xc potential
!

2
GIPAW/gipaw_module.f90
View File

@ -724,7 +724,7 @@ CONTAINS
! computes the total local potential (external+scf) on the smooth grid
call setlocal
call set_vrs (vrs, vltot, vr, nrxx, nspin, doublegrid)
call set_vrs (vrs, vltot, vr, kedtau, kedtaur, nrxx, nspin, doublegrid)
! compute the D for the pseudopotentials
call newd

4
Gamma/cg_setup.f90
View File

@ -14,7 +14,7 @@ SUBROUTINE cg_setup
USE kinds, ONLY: DP
USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau, amass
USE pwcom
USE scf, only : rho, rho_core, vr, vltot, vrs
USE scf, only : rho, rho_core, v, vltot, vrs, kedtau
USE uspp_param, ONLY: upf
USE mp_global, ONLY : kunit
USE wavefunctions_module, ONLY: evc
@ -39,7 +39,7 @@ SUBROUTINE cg_setup
!
! sum self-consistent part (vr) and local part (vltot) of potential
!
CALL set_vrs(vrs,vltot,vr,nrxx,nspin,doublegrid)
CALL set_vrs(vrs,vltot,v%of_r,kedtau, v%kin_r, nrxx,nspin,doublegrid)
!
! allocate memory for various arrays
!

6
PH/dvanqq.f90
View File

@ -20,7 +20,7 @@ subroutine dvanqq
!
USE ions_base, ONLY : nat, ityp, ntyp => nsp
use pwcom
use scf, only : vr, vltot
use scf, only : v, vltot
use noncollin_module, ONLY : noncolin, npol
USE kinds, only : DP
use phcom
@ -100,11 +100,11 @@ subroutine dvanqq
do is = 1, nspin
if (nspin.ne.4.or.is==1) then
do ir = 1, nrxx
veff (ir, is) = CMPLX (vltot (ir) + vr (ir, is), 0.d0)
veff (ir, is) = CMPLX (vltot (ir) + v%of_r (ir, is), 0.d0)
enddo
else
do ir = 1, nrxx
veff (ir, is) = CMPLX (vr (ir, is), 0.d0)
veff (ir, is) = CMPLX (v%of_r (ir, is), 0.d0)
enddo
endif
call cft3 (veff (1, is), nr1, nr2, nr3, nrx1, nrx2, nrx3, - 1)

4
PH/phq_setup.f90
View File

@ -51,7 +51,7 @@ subroutine phq_setup
USE klist, ONLY : xk, lgauss, degauss, ngauss, nks, nelec
USE ktetra, ONLY : ltetra, tetra
USE lsda_mod, ONLY : nspin, lsda, starting_magnetization
USE scf, ONLY : vr, vrs, vltot, rho, rho_core
USE scf, ONLY : v, vrs, vltot, rho, rho_core, kedtau
USE gvect, ONLY : nrxx, ngm
USE gsmooth, ONLY : doublegrid
USE symme, ONLY : nsym, s, ftau, irt, t_rev
@ -105,7 +105,7 @@ subroutine phq_setup
!
! 1) Computes the total local potential (external+scf) on the smooth grid
!
call set_vrs (vrs, vltot, vr, nrxx, nspin, doublegrid)
call set_vrs (vrs, vltot, v%of_r, kedtau, v%kin_r, nrxx, nspin, doublegrid)
!
! 2) Set non linear core correction stuff
!

10
PP/punch_plot.f90
View File

@ -53,7 +53,7 @@ SUBROUTINE punch_plot (filplot, plot_num, sample_bias, z, dz, &
USE lsda_mod, ONLY : nspin, current_spin
USE ener, ONLY : ehart
USE io_global, ONLY : stdout, ionode
USE scf, ONLY : rho, vltot, vr
USE scf, ONLY : rho, vltot, v
USE wvfct, ONLY : npw, nbnd, wg, igk, gamma_only
USE noncollin_module, ONLY : noncolin
@ -115,17 +115,17 @@ SUBROUTINE punch_plot (filplot, plot_num, sample_bias, z, dz, &
! The total self-consistent potential V_H+V_xc on output
!
IF (noncolin) THEN
call DCOPY (nrxx, vr, 1, raux, 1)
call DCOPY (nrxx, v%of_r, 1, raux, 1)
ELSE
IF (spin_component == 0) THEN
CALL DCOPY (nrxx, vr, 1, raux, 1)
CALL DCOPY (nrxx, v%of_r, 1, raux, 1)
DO is = 2, nspin
CALL DAXPY (nrxx, 1.0d0, vr (1, is), 1, raux, 1)
CALL DAXPY (nrxx, 1.0d0, v%of_r (1, is), 1, raux, 1)
ENDDO
CALL DSCAL (nrxx, 1.d0 / nspin, raux, 1)
ELSE
IF (nspin == 2) current_spin = spin_component
CALL DCOPY (nrxx, vr (1, current_spin), 1, raux, 1)
CALL DCOPY (nrxx, v%of_r (1, current_spin), 1, raux, 1)
ENDIF
ENDIF
CALL DAXPY (nrxx, 1.0d0, vltot, 1, raux, 1)

14
PP/pw2casino.f90
View File

@ -69,9 +69,9 @@ SUBROUTINE compute_casino
nrxx, g, gg, ecutwfc, gcutm, nl, igtongl
USE klist , ONLY: nks, nelec, xk
USE lsda_mod, ONLY: lsda, nspin
USE scf, ONLY: rho, rho_core, rhog_core
USE scf, ONLY: rho, rho_core, rhog_core, vnew
USE ldaU, ONLY : lda_plus_u, eth, Hubbard_lmax
USE vlocal, ONLY: vloc, vnew, strf
USE vlocal, ONLY: vloc, strf
USE wvfct, ONLY: npw, npwx, nbnd, gamma_only, igk, g2kin, wg, et
USE uspp, ONLY: nkb, vkb, dvan
USE uspp_param, ONLY: nh
@ -79,6 +79,7 @@ SUBROUTINE compute_casino
USE io_global, ONLY: stdout
USE io_files, ONLY: nd_nmbr, nwordwfc, iunwfc
USE wavefunctions_module, ONLY : evc
USE funct, ONLY : dft_is_meta
IMPLICIT NONE
INTEGER :: ig, ibnd, ik, io, na, j, ispin, nbndup, nbnddown, &
nk, ngtot, ig7, ikk, nt, ijkb0, ikb, ih, jh, jkb, at_num
@ -86,7 +87,7 @@ SUBROUTINE compute_casino
LOGICAL :: exst, found
REAL(DP) :: ek, eloc, enl, charge, etotefield
COMPLEX(DP), ALLOCATABLE :: aux(:), hpsi(:,:)
REAL(DP), allocatable :: v_h_new(:,:,:,:)
REAL(DP), allocatable :: v_h_new(:,:,:,:), kedtaur_new(:,:)
INTEGER :: ios
INTEGER, EXTERNAL :: atomic_number
REAL (DP), EXTERNAL :: ewald
@ -106,6 +107,11 @@ SUBROUTINE compute_casino
ALLOCATE (idx (4*npwx) )
ALLOCATE (igtog (4*npwx) )
if (lda_plus_u) allocate(v_h_new(2*Hubbard_lmax+1,2*Hubbard_lmax+1,nspin,nat))
if (dft_is_meta()) then
allocate (kedtaur_new(nrxx,nspin))
else
allocate (kedtaur_new(1,nspin))
endif
hpsi (:,:) = (0.d0, 0.d0)
idx(:) = 0
igtog(:) = 0
@ -222,7 +228,7 @@ SUBROUTINE compute_casino
! compute hartree and xc contribution
!
CALL v_of_rho( rho, rho_core, rhog_core, &
ehart, etxc, vtxc, eth, etotefield, charge, vnew, v_h_new )
ehart, etxc, vtxc, eth, etotefield, charge, vnew )
!
etot=(ek + (etxc-etxcc)+ehart+eloc+enl+ewld)
!

6
PP/work_function.f90
View File

@ -15,7 +15,7 @@ SUBROUTINE work_function (wf)
USE io_global, ONLY : stdout, ionode, ionode_id
USE ener, ONLY : ef
USE lsda_mod, ONLY : nspin, current_spin
USE scf, ONLY : rho, vltot, vr, rho_core, rhog_core
USE scf, ONLY : rho, vltot, v, rho_core, rhog_core
USE gvect
USE cell_base, ONLY : omega, alat
USE mp, ONLY : mp_bcast
@ -60,12 +60,12 @@ SUBROUTINE work_function (wf)
#endif
!
#ifdef __PARA
aux(:) = vltot(:) + vr(:,current_spin)
aux(:) = vltot(:) + v%of_r(:,current_spin)
CALL gather (aux, vaux1)
aux(:) = aux(:) - vxc(:,current_spin)
CALL gather (aux, vaux2)
#else
vaux1(1:nrxx) = vltot(1:nrxx) + vr(1:nrxx,current_spin)
vaux1(1:nrxx) = vltot(1:nrxx) + v%of_r(1:nrxx,current_spin)
vaux2(1:nrxx) = vaux1(1:nrxx) -vxc(1:nrxx,current_spin)
#endif
!

6
PW/addusforce.f90
View File

@ -21,7 +21,7 @@ subroutine addusforce (forcenl)
USE gvect, ONLY : nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, ngm, &
nl, nlm, gg, g, eigts1, eigts2, eigts3, ig1, ig2, ig3
USE lsda_mod, ONLY : nspin
USE scf, ONLY : vr, vltot
USE scf, ONLY : v, vltot
USE uspp, ONLY : becsum, okvan
USE uspp_param, ONLY : upf, lmaxq, nh, nhm
USE wvfct, ONLY : gamma_only
@ -52,9 +52,9 @@ subroutine addusforce (forcenl)
allocate (vg(nrxx))
do is = 1, nspin
if (nspin.eq.4.and.is.ne.1) then
vg (:) = vr(:,is)
vg (:) = v%of_r(:,is)
else
vg (:) = vltot (:) + vr (:, is)
vg (:) = vltot (:) + v%of_r (:, is)
endif
call cft3 (vg, nr1, nr2, nr3, nrx1, nrx2, nrx3, - 1)
aux (:, is) = vg (nl (:) ) * tpiba * (0.d0, -1.d0)

6
PW/addusstress.f90
View File

@ -19,7 +19,7 @@ subroutine addusstres (sigmanlc)
USE gvect, ONLY : nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, ngm, &
nl, nlm, gg, g, eigts1, eigts2, eigts3, ig1, ig2, ig3
USE lsda_mod, ONLY : nspin
USE scf, ONLY : vr, vltot
USE scf, ONLY : v, vltot
USE uspp, ONLY : becsum, okvan
USE uspp_param, ONLY : upf, lmaxq, nh, nhm
USE wvfct, ONLY : gamma_only
@ -73,11 +73,11 @@ subroutine addusstres (sigmanlc)
do is = 1, nspin
if ( nspin == 4 .and. is /= 1 ) then
!
vg(:) = vr(:,is)
vg(:) = v%of_r(:,is)
!
ELSE
!
vg(:) = vltot(:) + vr(:,is)
vg(:) = vltot(:) + v%of_r(:,is)
!
END IF
call cft3 (vg, nr1, nr2, nr3, nrx1, nrx2, nrx3, - 1)

11
PW/allocate_fft.f90
View File

@ -20,9 +20,8 @@ subroutine allocate_fft
USE gsmooth, ONLY : nr1s,nr2s,nr3s,nrxxs,ngms, nls, nlsm, doublegrid
USE ions_base, ONLY : nat
USE lsda_mod, ONLY : nspin
USE scf, ONLY : rho, vr, vltot, vrs, rho_core, rhog_core, &
kedtau, kedtaur, create_scf_type
USE vlocal, ONLY : vnew
USE scf, ONLY : rho, v, vnew, vltot, vrs, rho_core, rhog_core, &
kedtau, create_scf_type
USE wvfct, ONLY : gamma_only
USE noncollin_module, ONLY : pointlist, factlist, r_loc, &
report, i_cons, noncolin, npol
@ -64,16 +63,14 @@ subroutine allocate_fft
allocate (ig3( ngm))
call create_scf_type(rho)
allocate (vr( nrxx,nspin))
call create_scf_type(v)
call create_scf_type(vnew)
allocate (vltot( nrxx))
allocate (vnew ( nrxx, nspin))
allocate (rho_core( nrxx))
if (dft_is_meta() ) then
allocate ( kedtau(nrxxs,nspin) )
allocate ( kedtaur(nrxx,nspin) )
else
allocate ( kedtau(1,nspin) )
allocate ( kedtaur(1,nspin) )
end if
ALLOCATE( rhog_core( ngm ) )
allocate (psic( nrxx))

20
PW/allocate_nlpot.f90
View File

@ -31,7 +31,7 @@ subroutine allocate_nlpot
USE gvect, ONLY : ngm, gcutm, ecutwfc, g
USE klist, ONLY : xk, wk, ngk, nks, xqq
USE lsda_mod, ONLY : nspin
USE ldaU, ONLY : Hubbard_lmax, v_hub
USE ldaU, ONLY : Hubbard_lmax
USE scf, ONLY :
USE noncollin_module, ONLY : noncolin
USE wvfct, ONLY : npwx, npw, igk, g2kin
@ -46,7 +46,7 @@ subroutine allocate_nlpot
!
! a few local variables
!
integer :: nt, na, nb, ldim, nwfcm
integer :: nt, na, nb, nwfcm
! counters on atom type, atoms, beta functions
!
! calculate number of PWs for all kpoints
@ -110,20 +110,8 @@ subroutine allocate_nlpot
if (nkb > 0) allocate (vkb( npwx, nkb))
allocate (becsum( nhm * (nhm + 1)/2, nat, nspin))
!
! ... Allocate space for Hubbard potential
! ... These arrays are allocated ALWAYS even if lda_plus_u = .FALSE.
! ... This is needed since they are passed as arguments of mix_rho
! ... no matter lda_plus_u is .TRUE. or .FALSE. ( 23/10/2003 C.S. )
!
! if (lda_plus_u) then
!
ldim = 2 * Hubbard_lmax + 1
ALLOCATE( v_hub( ldim, ldim, nspin, nat ) )
!
! endif
!
! Calculate dimensions for array tab (including a possible factor
! coming from cell contraction during variable cell relaxation/MD)
! Calculate dimensions for array tab (including a possible factor
! coming from cell contraction during variable cell relaxation/MD)
!
nqx = INT( (sqrt (ecutwfc) / dq + 4) * cell_factor )

14
PW/clean_pw.f90
View File

@ -20,15 +20,15 @@ SUBROUTINE clean_pw( lflag )
USE klist, ONLY : ngk
USE reciprocal_vectors, ONLY : ig_l2g
USE symme, ONLY : irt
USE vlocal, ONLY : strf, vloc, vnew
USE vlocal, ONLY : strf, vloc
USE wvfct, ONLY : igk, g2kin, et, wg, btype, gamma_only
USE force_mod, ONLY : force
USE scf, ONLY : rho, vr, vltot, rho_core, rhog_core, &
vrs, kedtau, kedtaur, destroy_scf_type
USE scf, ONLY : rho, v, vltot, rho_core, rhog_core, &
vrs, kedtau, destroy_scf_type, vnew
USE wavefunctions_module, ONLY : evc, psic, psic_nc
USE us, ONLY : qrad, tab, tab_at, tab_d2y, spline_ps
USE uspp, ONLY : deallocate_uspp
USE ldaU, ONLY : swfcatom, v_hub
USE ldaU, ONLY : swfcatom
USE extfield, ONLY : forcefield
USE sticks, ONLY : dfftp, dffts
USE stick_base, ONLY : sticks_deallocate
@ -80,11 +80,10 @@ SUBROUTINE clean_pw( lflag )
IF ( ALLOCATED( ig2 ) ) DEALLOCATE( ig2 )
IF ( ALLOCATED( ig3 ) ) DEALLOCATE( ig3 )
call destroy_scf_type(rho)
call destroy_scf_type(v)
call destroy_scf_type(vnew)
IF ( ALLOCATED( kedtau ) ) DEALLOCATE( kedtau )
IF ( ALLOCATED( kedtaur ) ) DEALLOCATE( kedtaur )
IF ( ALLOCATED( vr ) ) DEALLOCATE( vr )
IF ( ALLOCATED( vltot ) ) DEALLOCATE( vltot )
IF ( ALLOCATED( vnew ) ) DEALLOCATE( vnew )
IF ( ALLOCATED( rho_core ) ) DEALLOCATE( rho_core )
IF ( ALLOCATED( rhog_core ) ) DEALLOCATE( rhog_core )
IF ( ALLOCATED( psic ) ) DEALLOCATE( psic )
@ -114,7 +113,6 @@ SUBROUTINE clean_pw( lflag )
IF ( ALLOCATED( igk ) ) DEALLOCATE( igk )
IF ( ALLOCATED( g2kin ) ) DEALLOCATE( g2kin )
IF ( ALLOCATED( qrad ) ) DEALLOCATE( qrad )
IF ( ALLOCATED( v_hub ) ) DEALLOCATE( v_hub )
IF ( ALLOCATED( tab ) ) DEALLOCATE( tab )
IF ( ALLOCATED( tab_at ) ) DEALLOCATE( tab_at )
IF ( lspinorb ) THEN

32
PW/electrons.f90
View File

@ -30,14 +30,14 @@ SUBROUTINE electrons()
USE gsmooth, ONLY : doublegrid, ngms
USE klist, ONLY : xk, wk, nelec, ngk, nks, nkstot, lgauss
USE lsda_mod, ONLY : lsda, nspin, magtot, absmag, isk
USE vlocal, ONLY : strf, vnew
USE vlocal, ONLY : strf
USE wvfct, ONLY : nbnd, et, gamma_only, npwx
USE ener, ONLY : etot, hwf_energy, eband, deband, ehart, &
vtxc, etxc, etxcc, ewld, demet
USE scf, ONLY : scf_type, &
create_scf_type, destroy_scf_type, &
rho, rho_core, rhog_core, &
vr, vltot, vrs, kedtau, kedtaur
v, vltot, vrs, kedtau, vnew
USE control_flags, ONLY : mixing_beta, tr2, ethr, niter, nmix, &
iprint, istep, lscf, lmd, conv_elec, &
restart, io_level, assume_isolated
@ -45,7 +45,7 @@ SUBROUTINE electrons()
iunefield
USE buffers, ONLY : save_buffer
USE ldaU, ONLY : eth, Hubbard_U, Hubbard_lmax, &
niter_with_fixed_ns, lda_plus_u, v_hub
niter_with_fixed_ns, lda_plus_u
USE extfield, ONLY : tefield, etotefield
USE wavefunctions_module, ONLY : evc, psic
USE noncollin_module, ONLY : noncolin, magtot_nc, i_cons, bfield, &
@ -324,7 +324,7 @@ SUBROUTINE electrons()
! ... charge density (i.e. the new estimate)
!
CALL v_of_rho( rhoin, rho_core, rhog_core, &
ehart, etxc, vtxc, eth, etotefield, charge, vr, v_hub)
ehart, etxc, vtxc, eth, etotefield, charge, v)
!
! ... estimate correction needed to have variational energy:
! ... T + E_ion (eband + deband) are calculated in sum_band
@ -365,12 +365,12 @@ SUBROUTINE electrons()
! ... 1) the output HXC-potential is saved in vr
! ... 2) vnew contains V(out)-V(in) ( used to correct the forces ).
!
vnew(:,:) = vr(:,:)
vnew%of_r(:,:) = v%of_r(:,:)
!
CALL v_of_rho( rho,rho_core,rhog_core, &
ehart, etxc, vtxc, eth, etotefield, charge, vr, v_hub)
ehart, etxc, vtxc, eth, etotefield, charge, v)
!
vnew(:,:) = vr(:,:) - vnew(:,:)
vnew%of_r(:,:) = v%of_r(:,:) - vnew%of_r(:,:)
!
! CHECKME: is it becsum or becstep??
IF (okpaw) &
@ -406,7 +406,7 @@ SUBROUTINE electrons()
!
! ... define the total local potential (external + scf)
!
CALL set_vrs( vrs, vltot, vr, nrxx, nspin, doublegrid )
CALL set_vrs( vrs, vltot, v%of_r, kedtau, v%kin_r, nrxx, nspin, doublegrid )
!
! ... in the US case we have to recompute the self-consistent
! ... term in the nonlocal potential
@ -475,9 +475,9 @@ SUBROUTINE electrons()
!
fock0 = exxenergy2()
CALL v_of_rho( rho, rho_core,rhog_core, &
ehart, etxc, vtxc, eth, etotefield, charge, vr, v_hub)
ehart, etxc, vtxc, eth, etotefield, charge, v)
!
CALL set_vrs( vrs, vltot, vr, nrxx, nspin, doublegrid )
CALL set_vrs( vrs, vltot, v%of_r, kedtau, v%kin_r, nrxx, nspin, doublegrid )
!
WRITE( stdout, * ) " NOW GO BACK TO REFINE HYBRID CALCULATION"
WRITE( stdout, * ) fock0
@ -783,13 +783,13 @@ SUBROUTINE electrons()
!
DO is = 1, nspin
!
delta_e = delta_e - SUM( rho%of_r(:,is)*vr(:,is) )
delta_e = delta_e - SUM( rho%of_r(:,is)*v%of_r(:,is) )
!
END DO
!
IF ( dft_is_meta() ) THEN
DO is = 1, nspin
delta_e = delta_e - SUM( rho%kin_r(:,is)*kedtaur(:,is) )
delta_e = delta_e - SUM( rho%kin_r(:,is)*v%kin_r(:,is) )
END DO
END IF
!
@ -798,7 +798,7 @@ SUBROUTINE electrons()
CALL mp_sum( delta_e, intra_pool_comm )
!
if (lda_plus_u) then
delta_e_hub = - SUM (rho%ns(:,:,:,:)*v_hub(:,:,:,:))
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
end if
@ -828,14 +828,14 @@ SUBROUTINE electrons()
DO is = 1, nspin
!
delta_escf = delta_escf - &
SUM( ( rhoin%of_r(:,is)-rho%of_r(:,is) )*vr(:,is) )
SUM( ( rhoin%of_r(:,is)-rho%of_r(:,is) )*v%of_r(:,is) )
!
END DO
!
IF ( dft_is_meta() ) THEN
DO is = 1, nspin
delta_escf = delta_escf - &
SUM( (rhoin%kin_r(:,is)-rho%kin_r(:,is) )*kedtaur(:,is))
SUM( (rhoin%kin_r(:,is)-rho%kin_r(:,is) )*v%kin_r(:,is))
END DO
END IF
!
@ -845,7 +845,7 @@ SUBROUTINE electrons()
!
!
if (lda_plus_u) then
delta_escf_hub = - SUM ((rhoin%ns(:,:,:,:)-rho%ns(:,:,:,:))*v_hub(:,:,:,:))
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
end if

6
PW/force_corr.f90
View File

@ -27,7 +27,7 @@ subroutine force_corr (forcescc)
USE gvect, ONLY : ngm, gstart, nr1, nr2, nr3, nrx1, nrx2, &
nrx3, nrxx, nl, g, ngl, gl, igtongl
USE lsda_mod, ONLY : nspin
USE vlocal, ONLY : vnew
USE scf, ONLY : vnew
USE wvfct, ONLY : gamma_only
USE wavefunctions_module, ONLY : psic
!
@ -45,11 +45,11 @@ subroutine force_corr (forcescc)
! vnew is V_out - V_in, psic is the temp space
!
if (nspin == 1 .or. nspin == 4) then
psic(:) = vnew (:, 1)
psic(:) = vnew%of_r (:, 1)
else
isup = 1
isdw = 2
psic(:) = (vnew (:, isup) + vnew (:, isdw)) * 0.5d0
psic(:) = (vnew%of_r (:, isup) + vnew%of_r (:, isdw)) * 0.5d0
end if
!
ndm = MAXVAL ( msh(1:ntyp) )

5
PW/force_hub.f90
View File

@ -21,11 +21,12 @@ SUBROUTINE force_hub(forceh)
USE cell_base, ONLY : at, bg
USE ldaU, ONLY : hubbard_lmax, hubbard_l, hubbard_u, &
hubbard_alpha, U_projection, &
swfcatom, v_hub
swfcatom
USE symme, ONLY : s, nsym, irt
USE io_files, ONLY : prefix, iunocc
USE wvfct, ONLY : gamma_only, nbnd, npwx, npw, igk
USE lsda_mod, ONLY : lsda, nspin, current_spin, isk
USE scf, ONLY : v
USE mp_global, ONLY : me_pool, my_pool_id
USE basis, ONLY : natomwfc
USE becmod, ONLY : becp
@ -125,7 +126,7 @@ SUBROUTINE force_hub(forceh)
DO m2 = 1,ldim
DO m1 = 1,ldim
forceh(ipol,alpha) = forceh(ipol,alpha) - &
v_hub(m2,m1,is,na) * dns(m1,m2,is,na)
v%ns(m2,m1,is,na) * dns(m1,m2,is,na)
END DO
END DO
END DO

4
PW/hinit1.f90
View File

@ -18,7 +18,7 @@ SUBROUTINE hinit1()
USE gsmooth, ONLY : doublegrid
USE ldaU, ONLY : lda_plus_u
USE lsda_mod, ONLY : nspin
USE scf, ONLY : vrs, vltot, vr
USE scf, ONLY : vrs, vltot, v, kedtau
USE vlocal, ONLY : strf
USE control_flags, ONLY : pot_order
USE realus, ONLY : tqr, qpointlist
@ -53,7 +53,7 @@ SUBROUTINE hinit1()
!
! ... define the total local potential (external+scf)
!
CALL set_vrs( vrs, vltot, vr, nrxx, nspin, doublegrid )
CALL set_vrs( vrs, vltot, v%of_r, kedtau, v%kin_r, nrxx, nspin, doublegrid )
!
IF ( tqr ) CALL qpointlist()
!

6
PW/newd.f90
View File

@ -31,7 +31,7 @@ SUBROUTINE newd_g()
g, gg, ngm, gstart, ig1, ig2, ig3, &
eigts1, eigts2, eigts3, nl
USE lsda_mod, ONLY : nspin
USE scf, ONLY : vr, vltot
USE scf, ONLY : v, vltot
USE uspp, ONLY : deeq, dvan, deeq_nc, dvan_so, okvan, indv
USE uspp_param, ONLY : upf, lmaxq, nh, nhm
USE wvfct, ONLY : gamma_only
@ -123,11 +123,11 @@ SUBROUTINE newd_g()
!
IF ( nspin0 == 4 .AND. is /= 1 ) THEN
!
psic(:) = vr(:,is)
psic(:) = v%of_r(:,is)
!
ELSE
!
psic(:) = vltot(:) + vr(:,is)
psic(:) = vltot(:) + v%of_r(:,is)
!
END IF
!

14
PW/potinit.f90
View File

@ -36,12 +36,12 @@ SUBROUTINE potinit()
USE gsmooth, ONLY : doublegrid
USE control_flags, ONLY : lscf
USE scf, ONLY : rho, rho_core, rhog_core, &
vltot, vr, vrs
USE funct, ONLY : dft_is_meta
vltot, v, vrs, kedtau
USE funct, ONLY : dft_is_meta
USE wavefunctions_module, ONLY : psic
USE ener, ONLY : ehart, etxc, vtxc
USE ldaU, ONLY : niter_with_fixed_ns
USE ldaU, ONLY : lda_plus_u, Hubbard_lmax, v_hub, eth
USE ldaU, ONLY : lda_plus_u, Hubbard_lmax, eth
USE noncollin_module, ONLY : noncolin, report
USE io_files, ONLY : tmp_dir, prefix, iunocc, input_drho
USE spin_orb, ONLY : domag
@ -143,13 +143,13 @@ SUBROUTINE potinit()
IF ( nspin > 1 ) CALL errore &
( 'potinit', 'spin polarization not allowed in drho', 1 )
!
CALL read_rho ( vr, 1, input_drho )
CALL read_rho ( v%of_r, 1, input_drho )
!
WRITE( UNIT = stdout, &
FMT = '(/5X,"a scf correction to at. rho is read from",A)' ) &
TRIM( input_drho )
!
rho%of_r = rho%of_r + vr
rho%of_r = rho%of_r + v%of_r
!
END IF
!
@ -214,11 +214,11 @@ SUBROUTINE potinit()
! ... compute the potential and store it in vr
!
CALL v_of_rho( rho, rho_core, rhog_core, &
ehart, etxc, vtxc, eth, etotefield, charge, vr, v_hub )
ehart, etxc, vtxc, eth, etotefield, charge, v )
!
! ... define the total local potential (external+scf)
!
CALL set_vrs( vrs, vltot, vr, nrxx, nspin, doublegrid )
CALL set_vrs( vrs, vltot, v%of_r, kedtau, v%kin_r, nrxx, nspin, doublegrid )
!
! ... write on output the parameters used in the lda+U calculation
!

7
PW/pwcom.f90
View File

@ -289,8 +289,7 @@ MODULE vlocal
COMPLEX(DP), ALLOCATABLE :: &
strf(:,:) ! the structure factor
REAL(DP), ALLOCATABLE :: &
vloc(:,:), &! the local potential for each atom type
vnew(:,:) ! V_out - V_in, needed in scf force correction
vloc(:,:) ! the local potential for each atom type
!
END MODULE vlocal
!
@ -463,8 +462,8 @@ MODULE ldaU
!
COMPLEX(DP), ALLOCATABLE :: &
swfcatom(:,:) ! orthogonalized atomic wfcs
REAL(DP), ALLOCATABLE :: &
v_hub(:,:,:,:) ! the hubbard contribution to the potential
! REAL(DP), ALLOCATABLE :: &
! v_hub(:,:,:,:) ! the hubbard contribution to the potential
REAL(DP) :: &
eth, &! the Hubbard contribution to the energy
Hubbard_U(ntypx), &! the Hubbard U

6
PW/read_file.f90
View File

@ -31,7 +31,7 @@ SUBROUTINE read_file()
nl, gstart
USE gsmooth, ONLY : ngms, nls, nrx1s, nr1s, nr2s, nr3s
USE spin_orb, ONLY : so, lspinorb
USE scf, ONLY : rho, rho_core, rhog_core, vr
USE scf, ONLY : rho, rho_core, rhog_core, v
USE wavefunctions_module, ONLY : psic
USE vlocal, ONLY : strf
USE io_files, ONLY : tmp_dir, prefix, iunpun, nwordwfc, iunwfc
@ -44,7 +44,7 @@ SUBROUTINE read_file()
USE uspp, ONLY : okvan
USE paw_variables, ONLY : okpaw
USE paw_init, ONLY : paw_init_onecenter, allocate_paw_internals
USE ldaU, ONLY : v_hub, eth
USE ldaU, ONLY : eth
!
IMPLICIT NONE
!
@ -213,7 +213,7 @@ SUBROUTINE read_file()
! ... recalculate the potential
!
CALL v_of_rho( rho, rho_core, rhog_core, &
ehart, etxc, vtxc, eth, etotefield, charge, vr, v_hub )
ehart, etxc, vtxc, eth, etotefield, charge, v )
!
! ... reads the wavefunctions and writes them in 'distributed' form
! ... to unit iunwfc (for compatibility)

6
PW/realus.f90
View File

@ -464,7 +464,7 @@ MODULE realus
USE cell_base, ONLY : omega
USE gvect, ONLY : nr1, nr2, nr3, nrxx
USE lsda_mod, ONLY : nspin
USE scf, ONLY : vr, vltot
USE scf, ONLY : v, vltot
USE uspp, ONLY : okvan, deeq, deeq_nc, dvan, dvan_so
USE uspp_param, ONLY : upf, nh, nhm
USE noncollin_module, ONLY : noncolin
@ -527,9 +527,9 @@ MODULE realus
DO is = 1, nspin0
!
IF ( nspin0 == 4 .AND. is /= 1 ) THEN
aux(:) = vr(:,is)
aux(:) = v%of_r(:,is)
ELSE
aux(:) = vltot(:) + vr(:,is)
aux(:) = vltot(:) + v%of_r(:,is)
END IF
!
iqs = 0

7
PW/restart_in_ions.f90
View File

@ -18,10 +18,9 @@ subroutine restart_in_ions (iter, ik_, dr2)
nrxx
USE klist, ONLY: nks
USE lsda_mod, ONLY: nspin
USE scf, ONLY : rho, rho_core, rhog_core
USE scf, ONLY : rho, rho_core, rhog_core, v, vnew
USE ldaU, ONLY : eth
USE control_flags, ONLY: restart, tr2, ethr
USE vlocal, ONLY: vnew
USE wvfct, ONLY: nbnd, et
USE wavefunctions_module, ONLY : evc, psic
implicit none
@ -54,7 +53,7 @@ subroutine restart_in_ions (iter, ik_, dr2)
read (iunres, err=10, end=10) ( (et(ibnd,ik), ibnd=1,nbnd), ik=1,nks)
read (iunres, err=10, end=10) etot, tr2
! vnew = V(in)-V(out) is needed in the scf correction term to forces
read (iunres, err=10, end=10) vnew
read (iunres, err=10, end=10) vnew%of_r
close (unit = iunres, status = 'keep')
WRITE( stdout, '(5x,"Calculation restarted from IONS ",i3)')
!
@ -82,7 +81,7 @@ subroutine restart_in_ions (iter, ik_, dr2)
! recalculate etxc, vtxc, ehart, needed by stress calculation
!
CALL v_of_rho( rho, rho_core, rhog_core, &
ehart, etxc, vtxc, eth, etotefield, charge, psic, psic )
ehart, etxc, vtxc, eth, etotefield, charge, v )
!
! restart procedure completed
!

4
PW/save_in_electrons.f90
View File

@ -14,7 +14,7 @@ subroutine save_in_electrons (iter, dr2)
USE klist, ONLY: nks
USE control_flags, ONLY: io_level, conv_elec, tr2, ethr
USE wvfct, ONLY: nbnd, et
USE vlocal, ONLY: vnew
USE scf, ONLY: vnew
implicit none
character :: where * 20
! are we in the right place?
@ -43,7 +43,7 @@ subroutine save_in_electrons (iter, dr2)
write (iunres) ( (et (ibnd, ik), ibnd = 1, nbnd), ik = 1, nks)
write (iunres) etot, tr2
! vnew = V(in)-V(out) is needed in the scf correction term to forces
write (iunres) vnew
write (iunres) vnew%of_r
else
!
! save iteration number

12
PW/scf_mod.f90
View File

@ -39,19 +39,17 @@ MODULE scf
REAL(DP), ALLOCATABLE :: ns(:,:,:,:)! the LDA+U occupation matrix
END TYPE mix_type
type (scf_type) :: rho
type (scf_type) :: rho ! the charge density and its other components
type (scf_type) :: v ! the scf potential
type (scf_type) :: vnew ! used to correct the forces
REAL(DP) :: v_of_0 ! vltot(G=0)
REAL(DP), ALLOCATABLE :: &
vr(:,:), &! the Hartree + xc potential in real space
vltot(:), &! the local potential in real space
vrs(:,:), &! the total pot. in real space (smooth grig)
rho_core(:), &! the core charge in real space
! tauk(:,:), &! kinetic energy density in real space (dense grid)
kedtau(:,:), &! position dependent kinetic energy enhancement factor
! used in META-GGA in real space (smooth grid)
kedtaur(:,:) ! position dependent kinetic energy enhancement factor
! used in META-GGA in real space (dense grid)
kedtau(:,:) ! position dependent kinetic energy enhancement factor
COMPLEX(DP), ALLOCATABLE :: &
rhog_core(:) ! the core charge in reciprocal space

8
PW/set_vrs.f90
View File

@ -6,19 +6,22 @@
! or http://www.gnu.org/copyleft/gpl.txt .
!
!--------------------------------------------------------------------
subroutine set_vrs (vrs, vltot, vr, nrxx, nspin, doublegrid)
subroutine set_vrs (vrs, vltot, vr, kedtau, kedtaur,nrxx, nspin, doublegrid)
!--------------------------------------------------------------------
! set the total local potential vrs on the smooth mesh to be used in
! h_psi, adding the (spin dependent) scf (H+xc) part and the sum of
! all the local pseudopotential contributions.
!
USE kinds
USE funct, only : dft_is_meta
USE gsmooth, only : nrxxs
implicit none
integer :: nspin, nrxx
! input: number of spin components: 1 if lda, 2 if lsd, 4 if noncolinear
! input: the fft grid dimension
real(DP) :: vrs (nrxx, nspin), vltot (nrxx), vr (nrxx, nspin)
real(DP) :: vrs (nrxx, nspin), vltot (nrxx), vr (nrxx, nspin), &
kedtau(nrxxs,nspin), kedtaur(nrxx,nspin)
! output: total local potential on the smooth grid
! vrs=vltot+vr
! input: the total local pseudopotential
@ -44,6 +47,7 @@ subroutine set_vrs (vrs, vltot, vr, nrxx, nspin, doublegrid)
! ... and interpolate it on the smooth mesh if necessary
!
if (doublegrid) call interpolate (vrs (1, is), vrs (1, is), - 1)
if (dft_is_meta()) call interpolate(kedtaur(1,is),kedtau(1,is),-1)
enddo
return

6
PW/stres_hub.f90
View File

@ -20,8 +20,8 @@ SUBROUTINE stres_hub ( sigmah )
USE ions_base, ONLY : nat, ityp
USE cell_base, ONLY : omega, at, bg
USE ldaU, ONLY : hubbard_lmax, hubbard_l, hubbard_u, &
hubbard_alpha, U_projection, v_hub
USE scf, ONLY :
hubbard_alpha, U_projection
USE scf, ONLY : v
USE lsda_mod, ONLY : nspin
USE symme, ONLY : s, nsym
USE io_files, ONLY : prefix, iunocc
@ -91,7 +91,7 @@ SUBROUTINE stres_hub ( sigmah )
DO m2 = 1, 2 * Hubbard_l(nt) + 1
DO m1 = 1, 2 * Hubbard_l(nt) + 1
sigmah(ipol,jpol) = sigmah(ipol,jpol) - omin1 * &
v_hub(m2,m1,is,na) * dns(m1,m2,is,na)
v%ns(m2,m1,is,na) * dns(m1,m2,is,na)
END DO
END DO
END DO

6
PW/update_pot.f90
View File

@ -192,8 +192,8 @@ SUBROUTINE extrapolate_charge( rho_extr )
USE gvect, ONLY : nrxx, ngm, g, gg, gstart, nr1, nr2, nr3, &
nl, eigts1, eigts2, eigts3, nrx1, nrx2, nrx3
USE lsda_mod, ONLY : lsda, nspin
USE scf, ONLY : rho, rho_core, rhog_core, vr
USE ldaU, ONLY : eth, v_hub
USE scf, ONLY : rho, rho_core, rhog_core, v
USE ldaU, ONLY : eth
USE wavefunctions_module, ONLY : psic
USE control_flags, ONLY : alpha0, beta0
USE ener, ONLY : ehart, etxc, vtxc
@ -380,7 +380,7 @@ SUBROUTINE extrapolate_charge( rho_extr )
END DO
!
CALL v_of_rho( rho, rho_core, rhog_core, &
ehart, etxc, vtxc, eth, etotefield, charge, vr, v_hub )
ehart, etxc, vtxc, eth, etotefield, charge, v )
!
IF ( ABS( charge - nelec ) / charge > 1.D-7 ) THEN
!

45
PW/v_of_rho.f90
View File

@ -7,7 +7,7 @@
!
!----------------------------------------------------------------------------
SUBROUTINE v_of_rho( rho, rho_core, rhog_core, &
ehart, etxc, vtxc, eth, etotefield, charge, v, v_hub )
ehart, etxc, vtxc, eth, etotefield, charge, v )
!----------------------------------------------------------------------------
!
! ... This routine computes the Hartree and Exchange and Correlation
@ -28,18 +28,19 @@ SUBROUTINE v_of_rho( rho, rho_core, rhog_core, &
IMPLICIT NONE
!
TYPE(scf_type), INTENT(IN) :: rho ! the valence charge
TYPE(scf_type), INTENT(INOUT) :: v ! the scf (Hxc) potential
!!!!!!!!!!!!!!!!! NB: NOTE that in F90 derived data type must be INOUT and
!!!!!!!!!!!!!!!!! not just OUT because otherwise their allocatable or pointer
!!!!!!!!!!!!!!!!! components are NOT defined !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
REAL(DP), INTENT(IN) :: rho_core(nrxx)
! input: the core charge
COMPLEX(DP), INTENT(IN) :: rhog_core(ngm)
! input: the core charge in reciprocal space
REAL(DP), INTENT(OUT) :: vtxc, etxc, ehart, eth, charge, etotefield, &
v(nrxx,nspin), &
v_hub(2*Hubbard_lmax+1,2*Hubbard_lmax+1,nspin,nat)
REAL(DP), INTENT(OUT) :: vtxc, etxc, ehart, eth, charge, etotefield
! output: the integral V_xc * rho
! output: the E_xc energy
! output: the hartree energy
! output: the integral of the charge
! output: the H+xc_up potential
!
INTEGER :: is
!
@ -48,26 +49,26 @@ SUBROUTINE v_of_rho( rho, rho_core, rhog_core, &
! ... calculate exchange-correlation potential
!
if (dft_is_meta()) then
call v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v )
call v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v%of_r, v%kin_r )
else
CALL v_xc( rho, rho_core, rhog_core, etxc, vtxc, v )
CALL v_xc( rho, rho_core, rhog_core, etxc, vtxc, v%of_r )
endif
!
! ... add a magnetic field (if any)
!
CALL add_bfield( v, rho%of_r )
CALL add_bfield( v%of_r, rho%of_r )
!
! ... calculate hartree potential
!
CALL v_h( rho%of_g, ehart, charge, v )
CALL v_h( rho%of_g, ehart, charge, v%of_r )
!
if (lda_plus_u) call v_Hubbard(rho%ns,v_hub,eth)
if (lda_plus_u) call v_Hubbard(rho%ns,v%ns,eth)
!
! ... add an electric field
!
DO is = 1, nspin
!
CALL add_efield( rho%of_r, v(1,is), etotefield, 0 )
CALL add_efield( rho%of_r, v%of_r(1,is), etotefield, 0 )
!
END DO
!
@ -77,7 +78,7 @@ SUBROUTINE v_of_rho( rho, rho_core, rhog_core, &
!
END SUBROUTINE v_of_rho
!
SUBROUTINE v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v )
SUBROUTINE v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v, kedtaur )
!----------------------------------------------------------------------------
!
! ... Exchange-Correlation potential Vxc(r) from n(r)
@ -99,7 +100,7 @@ SUBROUTINE v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v )
! the core charge
COMPLEX(DP), INTENT(IN) :: rhog_core(ngm)
! input: the core charge in reciprocal space
REAL(DP), INTENT(OUT) :: v(nrxx,nspin), vtxc, etxc
REAL(DP), INTENT(OUT) :: v(nrxx,nspin), kedtaur(nrxx,nspin), vtxc, etxc
! V_xc potential
! integral V_xc * rho
! E_xc energy
@ -130,7 +131,7 @@ SUBROUTINE v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v )
rhoneg = 0.D0
!
! IF (get_igcx()==7.AND.get_igcc()==6) THEN
call v_xc_tpss( rho, rho_core, rhog_core, etxc, vtxc, v )
call v_xc_tpss( rho, rho_core, rhog_core, etxc, vtxc, v, kedtaur )
! ELSE
! CALL errore('v_xc_meta','wrong igcx and/or igcc',1)
! ENDIF
@ -138,7 +139,7 @@ SUBROUTINE v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v )
RETURN
END SUBROUTINE v_xc_meta
!
SUBROUTINE v_xc_tpss( rho, rho_core, rhog_core, etxc, vtxc, v )
SUBROUTINE v_xc_tpss( rho, rho_core, rhog_core, etxc, vtxc, v, kedtaur )
! ===================
!--------------------------------------------------------------------
! use gvecp, only: ng => ngm
@ -146,7 +147,7 @@ SUBROUTINE v_xc_tpss( rho, rho_core, rhog_core, etxc, vtxc, v )
USE gvect, ONLY : nrxx, nrx1,nrx2,nrx3,nr1,nr2,nr3, &
g,nl,ngm
USE gsmooth, ONLY : nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, nrxxs
USE scf, ONLY : scf_type, kedtau, kedtaur
USE scf, ONLY : scf_type
USE lsda_mod, ONLY : nspin
USE cell_base, ONLY : omega, alat
USE constants, ONLY : e2
@ -156,7 +157,7 @@ SUBROUTINE v_xc_tpss( rho, rho_core, rhog_core, etxc, vtxc, v )
TYPE (scf_type), INTENT(IN) :: rho
REAL(DP),INTENT(IN) :: rho_core(nrxx)
COMPLEX(DP),INTENT(IN) :: rhog_core(ngm)
REAL(DP),INTENT(OUT) :: etxc, vtxc, v(nrxx,nspin)
REAL(DP),INTENT(OUT) :: etxc, vtxc, v(nrxx,nspin), kedtaur(nrxx,nspin)
! integer nspin , nnr
! real(8) grho(nnr,3,nspin), rho(nnr,nspin),kedtau(nnr,nspin)
! output: excrho: exc * rho ; E_xc = \int excrho(r) d_r
@ -269,13 +270,6 @@ SUBROUTINE v_xc_tpss( rho, rho_core, rhog_core, etxc, vtxc, v )
ENDIF
ENDDO
!
IF (nspin==1) then
CALL interpolate(kedtaur(1,1),kedtau(1,1),-1)
else
CALL interpolate(kedtaur(1,1),kedtau(1,1),-1)
CALL interpolate(kedtaur(1,2),kedtau(1,2),-1)
endif
!
ALLOCATE( dh( nrxx ) )
!
! ... second term of the gradient correction :
@ -637,8 +631,7 @@ SUBROUTINE v_hubbard(ns, v_hub, eth)
REAL(DP), INTENT(OUT) :: eth
INTEGER :: is, na, nt, m1, m2
!
! Now the contribution to the total energy is computed. The corrections
! needed to obtain a variational expression are already included
! Now the contribution to the total energy is computed.
!
eth = 0.d0
v_hub(:,:,:,:) = 0.d0

7
PW/vhpsi.f90
View File

@ -16,8 +16,9 @@ subroutine vhpsi (ldap, np, mp, psip, hpsi)
!
USE kinds, ONLY : DP
USE ldaU, ONLY : Hubbard_lmax, Hubbard_l, HUbbard_U, Hubbard_alpha, &
v_hub, swfcatom
swfcatom
USE lsda_mod, ONLY : nspin, current_spin
USE scf, ONLY : v
USE ions_base, ONLY : nat, ntyp => nsp, ityp
USE basis, ONLY : natomwfc
USE wvfct, ONLY : gamma_only
@ -72,8 +73,8 @@ subroutine vhpsi (ldap, np, mp, psip, hpsi)
do m1 = 1, 2 * Hubbard_l(nt) + 1
temp = 0.d0
do m2 = 1, 2 * Hubbard_l(nt) + 1
temp = temp + v_hub ( m1, m2, current_spin, na) * &
proj (offset(na)+m2, ibnd)
temp = temp + v%ns( m1, m2, current_spin, na) * &
proj (offset(na)+m2, ibnd)
enddo
if (gamma_only) then
call DAXPY (2*np, temp, swfcatom(1,offset(na)+m1), 1, &

7
PW/vhpsi_nc.f90
View File

@ -16,7 +16,8 @@ subroutine vhpsi_nc (ldap, np, mp, psip, hpsi)
!
USE kinds, ONLY: DP
USE ldaU, ONLY: Hubbard_lmax, Hubbard_l, Hubbard_U, Hubbard_alpha, &
v_hub, swfcatom
swfcatom
USE scf, ONLY: v
USE lsda_mod, ONLY: nspin, current_spin
USE ions_base, ONLY : nat, ntyp => nsp, ityp
USE basis, ONLY: natomwfc
@ -62,8 +63,8 @@ subroutine vhpsi_nc (ldap, np, mp, psip, hpsi)
do m1 = 1, 2 * Hubbard_l(nt) + 1
temp = 0.d0
do m2 = 1, 2 * Hubbard_l(nt) + 1
temp = temp + v_hub ( m1, m2, current_spin, na) * &
proj (offset(na)+m2, ibnd)
temp = temp + v%ns( m1, m2, current_spin, na) * &
proj (offset(na)+m2, ibnd)
enddo
call ZAXPY (np, temp, swfcatom(1,offset(na)+m1), 1, &
hpsi(1,1,ibnd), 1)

8
PWCOND/poten.f90
View File

@ -15,7 +15,7 @@ SUBROUTINE poten(vppot,nrz,z)
! local potential in each slab.
!
USE pwcom
USE scf, only : vltot, vr
USE scf, only : vltot, v
USE noncollin_module, ONLY : noncolin, npol
USE cond
USE mp, ONLY : mp_bcast
@ -101,12 +101,12 @@ vppot = 0.d0
DO ispin=1,nspin_eff
IF (noncolin) THEN
IF (ispin==1) THEN
auxr(:) = vltot(:)+vr(:,1)
auxr(:) = vltot(:)+v%of_r(:,1)
ELSE
auxr(:) = vr(:,ispin)
auxr(:) = v%of_r(:,ispin)
ENDIF
ELSE
auxr(:) = vltot(:) + vr(:,iofspin)
auxr(:) = vltot(:) + v%of_r(:,iofspin)
ENDIF
!
! To collect the potential from different CPUs

28
VdW/eff_pot.f90
View File

@ -17,7 +17,7 @@ subroutine eff_pot (rho, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, nl, &
! USE wvfct, ONLY : g2kin, wg, nbndx, et, nbnd, npwx, &
! igk, npw
USE uspp, ONLY : nkb
USE scf, ONLY : vr, vltot, vrs, rho_core
USE scf, ONLY : v, vltot, vrs, rho_core
USE gsmooth, ONLY : nls, nlsm, nr1s, nr2s, nr3s, nrx1s,&
nrx2s, nrx3s, nrxxs, doublegrid
USE eff_v, ONLY : rho_fft, veff
@ -37,7 +37,7 @@ subroutine eff_pot (rho, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, nl, &
!
! output
!
real(kind=DP), allocatable :: v (:,:), rho_in(:,:), k_gamma(:)
real(kind=DP), allocatable :: vv (:,:), rho_in(:,:), k_gamma(:)
!
! local variables
!
@ -56,7 +56,7 @@ subroutine eff_pot (rho, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, nl, &
call start_clock('eff_pot')
!
allocate (aux(2,nrxx), aux1(2,ngm), psi(2,nrxx), psi_smooth(2,nrxx),S(nrxx) )
allocate ( v(nrxx, nspin) )
allocate ( vv(nrxx, nspin) )
allocate ( rho_in(nrxx, nspin) )
!
tpiba2 = (fpi / 2.d0 / alat) **2
@ -149,12 +149,12 @@ subroutine eff_pot (rho, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, nl, &
is = 1
if (.false.) then
do ir = 1, nrxx
v(ir,is) = -aux(1,ir)
vv(ir,is) = -aux(1,ir)
end do
else
do ir = 1, nrxx
if (abs(psi_smooth(1,ir)) > eps ) then
v(ir,is) = -aux(1,ir) / psi_smooth(1,ir)
vv(ir,is) = -aux(1,ir) / psi_smooth(1,ir)
else
avg1 = avg1 - aux(1,ir)
avg2 = avg2 + psi_smooth(1,ir)
@ -169,7 +169,7 @@ subroutine eff_pot (rho, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, nl, &
if (nnn > 0 ) then
do ir = 1, nrxx
if (abs(psi_smooth(1,ir)) <= eps ) then
v(ir,is) = avg1 / avg2
vv(ir,is) = avg1 / avg2
end if
end do
end if
@ -179,8 +179,8 @@ subroutine eff_pot (rho, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, nl, &
!
vstart=20000
do ir = 1, nrxx
vrs(ir,1) = vr(ir,1) + vltot(ir)
v(ir,is) = erf(abs(psi_smooth(1,ir)*dble(vstart)))*v(ir,is) + &
vrs(ir,1) = v%of_r(ir,1) + vltot(ir)
vv(ir,is) = erf(abs(psi_smooth(1,ir)*dble(vstart)))*vv(ir,is) + &
(1.d0-erf( abs( psi_smooth(1,ir)*dble(vstart) ) ))*&
vrs(ir,is)
@ -188,7 +188,7 @@ vstart=20000
!
! check the quality of trial potential
!
CALL check_v_eff(v(1,1), charge)
CALL check_v_eff(vv(1,1), charge)
WRITE( stdout, '(/,10x,"Charge difference due to V_eff (initial) ",f10.8,/)' ) charge
!
! iterative procedure of V_eff if necessary
@ -204,7 +204,7 @@ vstart=20000
! whose the GS solution should be the square root
! of the charge density
!
call check_v_eff(v, charge)
call check_v_eff(vv, charge)
nite = nite + 1
!#ifdef __PARA
! call ireduce(1, nite)
@ -226,9 +226,9 @@ vstart=20000
! set the optmized eff. potential to veff
!
100 continue
veff(:,:) = v(:,:)
veff(:,:) = vv(:,:)
!
deallocate ( v, rho_in )
deallocate ( vv, rho_in )
deallocate (aux,aux1,psi,psi_smooth,S)
!
call stop_clock('eff_pot')
@ -255,7 +255,7 @@ vstart=20000
!CALL flush_unit( stdout )
s2 = 0.d0
do ir = 1, nrxx
S(ir) = psi_smooth(1,ir) * ( aux(1,ir) + v(ir,1)*psi_smooth(1,ir) )
S(ir) = psi_smooth(1,ir) * ( aux(1,ir) + vv(ir,1)*psi_smooth(1,ir) )
s2 = s2 + S(ir)**2
enddo
#ifdef __PARA
@ -314,7 +314,7 @@ vstart=20000
! Update V-eff
!
do ir = 1, nrxx
v(ir,1) = v(ir,1) + alp*S(ir) + beta
vv(ir,1)= vv(ir,1) + alp*S(ir) + beta
S(ir) = S(ir) * (1.d0 + alp*psi_smooth(1,ir)**2) + &
beta*psi_smooth(1,ir)**2
enddo

2
examples/example27/run_example
View File

@ -117,7 +117,7 @@ ATOMIC_POSITIONS bohr
C 0.109484E+01 0.137081E+01 -0.496954E+00
EOF
$ECHO " running the cp.x SCF calculation...\c"
$CP_COMMAND < c4h6.cp.metaGGA.in > c4h6.cp.metaGGA.out
#$CP_COMMAND < c4h6.cp.metaGGA.in > c4h6.cp.metaGGA.out
check_failure $?
$ECHO " done"