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

Offsets of atomic wavefunctions used for LDA+U projections is now computed only once in setup and stored, 

instead of recomputing them several times (they do not change during the run).
GS


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@5997 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
sclauzer 2009-10-07 13:11:59 +00:00
parent 35f89b13ff
commit 53385d4cc5
7 changed files with 117 additions and 104 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
PW/clean_pw.f90
View File

@ -29,7 +29,7 @@ SUBROUTINE clean_pw( lflag )
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
USE ldaU, ONLY : lda_plus_u, oatwfc, swfcatom
USE extfield, ONLY : forcefield
USE fft_base, ONLY : dfftp, dffts
USE stick_base, ONLY : sticks_deallocate
@ -140,7 +140,13 @@ SUBROUTINE clean_pw( lflag )
! ... arrays allocated in allocate_wfc.f90 ( and never deallocated )
!
IF ( ALLOCATED( evc ) ) DEALLOCATE( evc )
IF ( ALLOCATED( swfcatom ) ) DEALLOCATE( swfcatom )
!
! ... arrays allocated for LDA+U calculations
!
IF ( lda_plus_u ) THEN
IF ( ALLOCATED( oatwfc ) ) DEALLOCATE( oatwfc )
IF ( ALLOCATED( swfcatom ) ) DEALLOCATE( swfcatom )
ENDIF
!
! ... fft structures allocated in data_structure.f90
!

32
PW/force_hub.f90
View File

@ -20,7 +20,7 @@ SUBROUTINE force_hub(forceh)
USE cell_base, ONLY : at, bg
USE ldaU, ONLY : hubbard_lmax, hubbard_l, hubbard_u, &
hubbard_alpha, U_projection, &
swfcatom
swfcatom, oatwfc
USE symme, ONLY : s, nsym, irt
USE io_files, ONLY : prefix, iunocc
USE wvfct, ONLY : nbnd, npwx, npw, igk
@ -32,7 +32,6 @@ SUBROUTINE force_hub(forceh)
USE basis, ONLY : natomwfc
USE becmod, ONLY : bec_type, becp, calbec, allocate_bec_type, deallocate_bec_type
USE uspp, ONLY : nkb, vkb
USE uspp_param, ONLY : upf
USE wavefunctions_module, ONLY : evc
USE klist, ONLY : nks, xk, ngk
USE io_files, ONLY : iunigk, nwordwfc, iunwfc, &
@ -47,39 +46,22 @@ SUBROUTINE force_hub(forceh)
! spsi(npwx,nbnd)
REAL (DP), ALLOCATABLE :: dns(:,:,:,:)
! dns(ldim,ldim,nspin,nat) ! the derivative of the atomic occupations
INTEGER, ALLOCATABLE :: offset(:)
! offset(nat) : offset of d electrons of atom d in the natomwfc ordering
COMPLEX (DP) :: c_one, c_zero
INTEGER :: alpha, na, nt, is, m1, m2, ipol, ldim, l, n, ik
INTEGER :: counter
INTEGER :: alpha, na, nt, is, m1, m2, ipol, ldim, ik
IF (U_projection .NE. "atomic") CALL errore("force_hub", &
" forces for this U_projection_type not implemented",1)
ldim= 2 * Hubbard_lmax + 1
ALLOCATE ( dns(ldim,ldim,nspin,nat), offset(nat), spsi(npwx,nbnd) )
ALLOCATE ( dns(ldim,ldim,nspin,nat), spsi(npwx,nbnd) )
call allocate_bec_type ( nkb, nbnd, becp)
call allocate_bec_type ( natomwfc, nbnd, proj )
forceh(:,:) = 0.d0
counter = 0
DO na=1,nat
offset(na) = 0
nt=ityp(na)
DO n=1,upf(nt)%nwfc
IF (upf(nt)%oc(n) >= 0.d0) THEN
l=upf(nt)%lchi(n)
IF (l == Hubbard_l(nt)) offset(na) = counter
counter = counter + 2 * l + 1
END IF
END DO
END DO
IF (counter /= natomwfc) &
CALL errore('new_ns','Internal error: nstart<>counter',1)
! Offset of atomic wavefunctions initialized in setup and stored in oatwfc
!
! we start a loop on k points
!
@ -106,9 +88,9 @@ SUBROUTINE force_hub(forceh)
DO ipol = 1,3
DO alpha = 1,nat ! the displaced atom
IF ( gamma_only ) THEN
CALL dndtau_gamma(ldim,offset,proj%r,swfcatom,spsi,alpha,ipol,dns)
CALL dndtau_gamma(ldim,oatwfc,proj%r,swfcatom,spsi,alpha,ipol,dns)
ELSE
CALL dndtau_k (ldim,offset,proj%k,swfcatom,spsi,alpha,ipol,ik,dns)
CALL dndtau_k (ldim,oatwfc,proj%k,swfcatom,spsi,alpha,ipol,ik,dns)
ENDIF
DO na = 1,nat ! the Hubbard atom
nt = ityp(na)
@ -131,7 +113,7 @@ SUBROUTINE force_hub(forceh)
CALL mp_sum( forceh, inter_pool_comm )
#endif
DEALLOCATE(dns, offset, spsi)
DEALLOCATE(dns, spsi)
call deallocate_bec_type (proj)
call deallocate_bec_type (becp)

33
PW/new_ns.f90
View File

@ -20,7 +20,7 @@ SUBROUTINE new_ns(ns)
USE ions_base, ONLY : nat, ityp
USE basis, ONLY : natomwfc
USE klist, ONLY : nks, ngk
USE ldaU, ONLY : Hubbard_lmax, Hubbard_l, &
USE ldaU, ONLY : Hubbard_lmax, Hubbard_l, oatwfc, &
Hubbard_U, Hubbard_alpha, swfcatom
USE symme, ONLY : d1, d2, d3
USE lsda_mod, ONLY : lsda, current_spin, nspin, isk
@ -31,7 +31,6 @@ SUBROUTINE new_ns(ns)
USE gvect, ONLY : gstart
USE io_files, ONLY : iunigk, nwordwfc, iunwfc, nwordatwfc, iunsat
USE buffers, ONLY : get_buffer
USE uspp_param, ONLY : upf
USE mp_global, ONLY : intra_pool_comm, inter_pool_comm
USE mp, ONLY : mp_sum
@ -42,13 +41,11 @@ SUBROUTINE new_ns(ns)
!
REAL(DP) :: ns(2*Hubbard_lmax+1,2*Hubbard_lmax+1,nspin,nat)
INTEGER :: ik, ibnd, is, i, na, nb, nt, isym, n, counter, m1, m2, &
m0, m00, l, ldim
INTEGER, ALLOCATABLE :: offset (:)
INTEGER :: ik, ibnd, is, i, na, nb, nt, isym, m1, m2, &
m0, m00, ldim
! counter on k points
! " " bands
! " " spins
! offset of d electrons of atom d
! in the natomwfc ordering
REAL(DP) , ALLOCATABLE :: nr (:,:,:,:)
REAL(DP) :: t0, scnds
@ -62,24 +59,12 @@ SUBROUTINE new_ns(ns)
t0 = scnds ()
ldim = 2 * Hubbard_lmax + 1
ALLOCATE( offset(nat), proj(natomwfc,nbnd), nr(ldim,ldim,nspin,nat) )
ALLOCATE( proj(natomwfc,nbnd), nr(ldim,ldim,nspin,nat) )
!
! D_Sl for l=1, l=2 and l=3 are already initialized, for l=0 D_S0 is 1
!
counter = 0
DO na = 1, nat
nt = ityp (na)
DO n = 1, upf(nt)%nwfc
IF (upf(nt)%oc (n) >= 0.d0) THEN
l = upf(nt)%lchi (n)
IF (l == Hubbard_l(nt)) offset (na) = counter
counter = counter + 2 * l + 1
ENDIF
ENDDO
ENDDO
IF (counter.NE.natomwfc) CALL errore ('new_ns', 'nstart<>counter', 1)
! Offset of atomic wavefunctions initialized in setup and stored in oatwfc
!
nr (:,:,:,:) = 0.d0
ns (:,:,:,:) = 0.d0
!
@ -127,8 +112,8 @@ SUBROUTINE new_ns(ns)
DO m2 = m1, 2 * Hubbard_l(nt) + 1
DO ibnd = 1, nbnd
nr(m1,m2,current_spin,na) = nr(m1,m2,current_spin,na) + &
wg(ibnd,ik) * DBLE( proj(offset(na)+m2,ibnd) * &
CONJG(proj(offset(na)+m1,ibnd)) )
wg(ibnd,ik) * DBLE( proj(oatwfc(na)+m2,ibnd) * &
CONJG(proj(oatwfc(na)+m1,ibnd)) )
ENDDO
ENDDO
ENDDO
@ -219,7 +204,7 @@ SUBROUTINE new_ns(ns)
ENDIF
ENDDO
DEALLOCATE ( offset, proj, nr )
DEALLOCATE ( proj, nr )
RETURN

2
PW/pwcom.f90
View File

@ -437,6 +437,8 @@ MODULE ldaU
conv_ns ! .TRUE. if ns are converged
CHARACTER(LEN=30) :: & ! 'atomic', 'ortho-atomic', 'file'
U_projection ! specifies how input coordinates are given
INTEGER, ALLOCATABLE :: &
oatwfc(:) ! offset of atomic wfcs used for projections
!
END MODULE ldaU
!

71
PW/setup.f90
View File

@ -64,7 +64,7 @@ SUBROUTINE setup()
USE uspp_param, ONLY : upf
USE uspp, ONLY : okvan
USE ldaU, ONLY : lda_plus_u, Hubbard_U, &
Hubbard_l, Hubbard_alpha, Hubbard_lmax
Hubbard_l, Hubbard_alpha, Hubbard_lmax, oatwfc
USE bp, ONLY : gdir, lberry, nppstr, lelfield, nx_el, nppstr_3d,l3dstring, efield
USE fixed_occ, ONLY : f_inp, tfixed_occ
USE funct, ONLY : set_dft_from_name
@ -710,6 +710,9 @@ SUBROUTINE setup()
CALL errore( 'setup', &
& 'lda_plus_u calculation but Hubbard_l not set', 1 )
!
ALLOCATE ( oatwfc(nat) )
CALL offset_atom_wfc ( nat, ntyp, ityp, natomwfc, Hubbard_l, oatwfc )
!
ELSE
!
Hubbard_lmax = 0
@ -777,6 +780,72 @@ FUNCTION n_atom_wfc( nat, ityp )
RETURN
!
END FUNCTION n_atom_wfc
!
!----------------------------------------------------------------------------
SUBROUTINE offset_atom_wfc( nat, ntyp, ityp, ntotwfc, Hubbard_l, offset )
!----------------------------------------------------------------------------
!
! For each Hubbard atom, compute the index of the projector in the
! list of atomic wavefunctions
!
USE uspp_param, ONLY : upf
USE noncollin_module, ONLY : noncolin
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: nat, ntyp, ityp(nat), ntotwfc, Hubbard_l(ntyp)
!
INTEGER, INTENT(OUT) :: offset(nat)
!
INTEGER :: counter, na, nt, n
!
!
counter = 0
offset(:) = -99
!
!
DO na = 1, nat
!
nt = ityp(na)
!
DO n = 1, upf(nt)%nwfc
!
IF ( upf(nt)%oc(n) >= 0.D0 ) THEN
!
IF ( noncolin ) THEN
! N.B.: presently LDA+U not yet implemented for noncolin
!
IF ( upf(nt)%has_so ) THEN
!
counter = counter + 2 * upf(nt)%lchi(n)
!
IF ( ABS( upf(nt)%jchi(n)-upf(nt)%lchi(n) - 0.5D0 ) < 1.D-6 ) &
counter = counter + 2
!
ELSE
!
counter = counter + 2 * ( 2 * upf(nt)%lchi(n) + 1 )
!
END IF
!
ELSE
!
IF ( upf(nt)%lchi(n) == Hubbard_l(nt) ) offset(na) = counter
!
counter = counter + 2 * upf(nt)%lchi(n) + 1
!
END IF
END IF
END DO
END DO
!
IF (counter.NE.ntotwfc) &
CALL errore ('offset_atom_wfc', 'wrong number of wavefunctions', 1)
!
RETURN
!
END SUBROUTINE offset_atom_wfc
SUBROUTINE check_para_diag( nelec )
!

44
PW/stres_hub.f90
View File

@ -133,11 +133,10 @@ SUBROUTINE dndepsilon ( dns,ldim,ipol,jpol )
USE control_flags, ONLY : gamma_only
USE klist, ONLY : nks, xk, ngk
USE ldaU, ONLY : swfcatom, Hubbard_l, &
Hubbard_U, Hubbard_alpha
Hubbard_U, Hubbard_alpha, oatwfc
USE lsda_mod, ONLY : lsda, nspin, current_spin, isk
USE wvfct, ONLY : nbnd, npwx, npw, igk, wg
USE uspp, ONLY : nkb, vkb
USE uspp_param, ONLY : upf
USE becmod, ONLY : bec_type, becp, calbec, &
allocate_bec_type, deallocate_bec_type
USE io_files, ONLY : iunigk, nwordwfc, iunwfc, &
@ -158,38 +157,23 @@ SUBROUTINE dndepsilon ( dns,ldim,ipol,jpol )
INTEGER :: ik, & ! counter on k points
ibnd, & ! " " bands
is, & ! " " spins
i, na, nt, n, counter, m1, m2, l
na, nt, m1, m2
INTEGER, ALLOCATABLE :: offset(:)
! offset(nat) ! offset of d electrons of atom d in the natomwfc ordering
COMPLEX (DP), ALLOCATABLE :: spsi(:,:)
type (bec_type) :: proj, dproj
! COMPLEX (DP), ALLOCATABLE :: dproj(:,:)
! REAL (DP), ALLOCATABLE :: drproj(:,:)
!
!
ALLOCATE (offset(nat), spsi(npwx,nbnd) )
ALLOCATE ( spsi(npwx,nbnd) )
call allocate_bec_type( natomwfc,nbnd, proj)
call allocate_bec_type ( natomwfc,nbnd, dproj )
call allocate_bec_type ( nkb,nbnd, becp )
!
! D_Sl for l=1 and l=2 are already initialized, for l=0 D_S0 is 1
!
counter = 0
DO na=1,nat
offset(na) = 0
nt=ityp(na)
DO n=1,upf(nt)%nwfc
IF (upf(nt)%oc(n) >= 0.d0) THEN
l=upf(nt)%lchi(n)
IF (l.EQ.Hubbard_l(nt)) offset(na) = counter
counter = counter + 2 * l + 1
END IF
END DO
END DO
IF(counter.NE.natomwfc) CALL errore('new_ns','nstart<>counter',1)
! Offset of atomic wavefunctions initialized in setup and stored in oatwfc
dns(:,:,:,:) = 0.d0
!
! we start a loop on k points
@ -230,19 +214,19 @@ SUBROUTINE dndepsilon ( dns,ldim,ipol,jpol )
DO ibnd = 1,nbnd
dns(m1,m2,current_spin,na) = &
dns(m1,m2,current_spin,na) + wg(ibnd,ik) *&
(proj%r(offset(na)+m1,ibnd) * &
dproj%r(offset(na)+m2,ibnd) + &
dproj%r(offset(na)+m1,ibnd) * &
proj%r(offset(na)+m2,ibnd))
(proj%r(oatwfc(na)+m1,ibnd) * &
dproj%r(oatwfc(na)+m2,ibnd) + &
dproj%r(oatwfc(na)+m1,ibnd) * &
proj%r(oatwfc(na)+m2,ibnd))
END DO
ELSE
DO ibnd = 1,nbnd
dns(m1,m2,current_spin,na) = &
dns(m1,m2,current_spin,na) + wg(ibnd,ik) *&
DBLE(proj%k(offset(na)+m1,ibnd) * &
CONJG(dproj%k(offset(na)+m2,ibnd) ) + &
dproj%k(offset(na)+m1,ibnd)* &
CONJG(proj%k(offset(na)+m2,ibnd) ) )
DBLE(proj%k(oatwfc(na)+m1,ibnd) * &
CONJG(dproj%k(oatwfc(na)+m2,ibnd) ) + &
dproj%k(oatwfc(na)+m1,ibnd)* &
CONJG(proj%k(oatwfc(na)+m2,ibnd) ) )
END DO
END IF
END DO
@ -273,7 +257,7 @@ SUBROUTINE dndepsilon ( dns,ldim,ipol,jpol )
END DO
END DO
DEALLOCATE (offset, spsi)
DEALLOCATE ( spsi )
call deallocate_bec_type (proj)
call deallocate_bec_type (dproj)
call deallocate_bec_type (becp)

29
PW/vhpsi.f90
View File

@ -16,13 +16,12 @@ subroutine vhpsi (ldap, np, mps, psip, hpsi)
USE kinds, ONLY : DP
USE becmod, ONLY : bec_type, calbec, allocate_bec_type, deallocate_bec_type
USE ldaU, ONLY : Hubbard_lmax, Hubbard_l, HUbbard_U, Hubbard_alpha, &
swfcatom
swfcatom, oatwfc
USE lsda_mod, ONLY : nspin, current_spin
USE scf, ONLY : v
USE ions_base, ONLY : nat, ntyp => nsp, ityp
USE basis, ONLY : natomwfc
USE gvect, ONLY : gstart
USE uspp_param,ONLY : upf
USE control_flags, ONLY : gamma_only
USE mp_global, ONLY: intra_pool_comm
USE mp, ONLY: mp_sum
@ -33,26 +32,12 @@ subroutine vhpsi (ldap, np, mps, psip, hpsi)
complex(DP), intent(in) :: psip (ldap, mps)
complex(DP), intent(inout) :: hpsi (ldap, mps)
!
integer :: ibnd, i, na, nt, n, counter, m1, m2, l
integer, allocatable :: offset (:)
! offset of localized electrons of atom na in the natomwfc ordering
integer :: ibnd, na, nt, m1, m2
complex(DP) :: temp
type (bec_type) :: proj
!
allocate ( offset(nat) )
counter = 0
do na = 1, nat
nt = ityp (na)
do n = 1, upf(nt)%nwfc
if (upf(nt)%oc (n) >= 0.d0) then
l = upf(nt)%lchi (n)
if (l.eq.Hubbard_l(nt)) offset (na) = counter
counter = counter + 2 * l + 1
endif
enddo
enddo
! Offset of atomic wavefunctions initialized in setup and stored in oatwfc
!
if (counter /= natomwfc) call errore ('vhpsi', 'nstart<>counter', 1)
call allocate_bec_type ( natomwfc,mps, proj )
CALL calbec (np, swfcatom, psip, proj)
do ibnd = 1, mps
@ -64,16 +49,16 @@ subroutine vhpsi (ldap, np, mps, psip, hpsi)
if (gamma_only) then
do m2 = 1, 2 * Hubbard_l(nt) + 1
temp = temp + v%ns( m1, m2, current_spin, na) * &
proj%r(offset(na)+m2, ibnd)
proj%r(oatwfc(na)+m2, ibnd)
enddo
call daxpy (2*np, temp, swfcatom(1,offset(na)+m1), 1, &
call daxpy (2*np, temp, swfcatom(1,oatwfc(na)+m1), 1, &
hpsi(1,ibnd), 1)
else
do m2 = 1, 2 * Hubbard_l(nt) + 1
temp = temp + v%ns( m1, m2, current_spin, na) * &
proj%k(offset(na)+m2, ibnd)
proj%k(oatwfc(na)+m2, ibnd)
enddo
call zaxpy (np, temp, swfcatom(1,offset(na)+m1), 1, &
call zaxpy (np, temp, swfcatom(1,oatwfc(na)+m1), 1, &
hpsi(1,ibnd), 1)
endif
enddo