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

One more bug fix in the paw+exx code (the atomic's index of projector was computed incorrectly when the crystal indeces where given in a certain order). 

Now everything seems to work quite reliably, a few calculation with PAW on molecules where successful, however:

there is still a nasty bug that cause wrong results (usually in the form of S matrix not positive, or wrong charge) when ALL this coditions are met:
1. PAW and/or US
2. using Gygi-Baldereschi divergence threatment with x_gamma_extrapolation
3. K_POINTS gamma

I hope to fix it soonish
LP



git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@13139 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
paulatto 2016-11-09 14:54:16 +00:00
parent 9489e1ff28
commit 01b6e1f73e
4 changed files with 93 additions and 72 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
PW/src/exx.f90
View File

@ -674,7 +674,7 @@ MODULE exx
USE uspp, ONLY : nkb, vkb, okvan
USE us_exx, ONLY : rotate_becxx
USE paw_variables, ONLY : okpaw
USE paw_exx, ONLY : PAW_init_keeq
USE paw_exx, ONLY : PAW_init_fock_kernel
USE mp_pools, ONLY : me_pool, my_pool_id, root_pool, nproc_pool, &
inter_pool_comm, my_pool_id, intra_pool_comm
USE mp_orthopools, ONLY : intra_orthopool_comm
@ -928,7 +928,7 @@ MODULE exx
!
! Initialize 4-wavefunctions one-center Fock integrals
! \int \psi_a(r)\phi_a(r)\phi_b(r')\psi_b(r')/|r-r'|
IF(okpaw) CALL PAW_init_keeq()
IF(okpaw) CALL PAW_init_fock_kernel()
!
#if defined(__EXX_ACE)
CALL aceinit ( )

153
PW/src/paw_exx.f90
View File

@ -12,11 +12,11 @@
MODULE paw_exx
!=----------------------------------------------------------------------------=!
USE kinds, ONLY : DP
TYPE paw_keeq_type
TYPE paw_fockrnl_type
REAL(DP),POINTER :: k(:,:,:,:)
END TYPE paw_keeq_type
TYPE(paw_keeq_type),ALLOCATABLE :: ke(:)
LOGICAL,PRIVATE :: paw_has_init_keeq = .false.
END TYPE paw_fockrnl_type
TYPE(paw_fockrnl_type),ALLOCATABLE :: ke(:)
LOGICAL,PRIVATE :: paw_has_init_paw_fockrnl = .false.
CONTAINS
!
@ -46,8 +46,8 @@ MODULE paw_exx
!
!RETURN
!
IF(.not.paw_has_init_keeq) &
CALL errore("PAW_newdxx", "you have to initialize paw keeq before", 1)
IF(.not.paw_has_init_paw_fockrnl) &
CALL errore("PAW_newdxx", "you have to initialize paw paw_fockrnl before", 1)
!
CALL start_clock( 'PAW_newdxx' )
!
@ -112,8 +112,8 @@ MODULE paw_exx
INTEGER :: np, na
INTEGER :: ih, jh, oh, uh
INTEGER :: ikb, jkb, okb, ukb, ijkb0
IF(.not.paw_has_init_keeq) &
CALL errore("PAW_xx_energy", "you have to initialize paw keeq before", 1)
IF(.not.paw_has_init_paw_fockrnl) &
CALL errore("PAW_xx_energy", "you have to initialize paw paw_fockrnl before", 1)
!
CALL start_clock("PAW_xx_nrg")
PAW_xx_energy = 0._dp
@ -160,114 +160,117 @@ MODULE paw_exx
!=----------------------------------------------------------------------------=!
!
!=----------------------------------------------------------------------------=!
SUBROUTINE PAW_init_keeq()
SUBROUTINE PAW_init_fock_kernel()
!=----------------------------------------------------------------------------=!
! Driver to compute the 2-electron 4-wavefunctions integrals
! Driver to compute the 2-electron 4-wavefunctions integrals that constitue
! the kernel of the Fock operator
USE kinds, ONLY : DP
USE ions_base, ONLY : ntyp => nsp
USE uspp_param, ONLY : nh
IMPLICIT NONE
INTEGER :: ns, ih,jh,oh,uh
INTEGER :: ns !, ih,jh,oh,uh
REAL(DP),ALLOCATABLE :: k_ae(:,:,:,:), k_ps(:,:,:,:)
IF(paw_has_init_keeq) RETURN !CALL errore("PAW_init_keeq", "already init paw keeq", 1)
paw_has_init_keeq = .true.
IF(paw_has_init_paw_fockrnl) RETURN !CALL errore("PAW_init_fock_kernel", "already init paw paw_fockrnl", 1)
paw_has_init_paw_fockrnl = .true.
!
! We have one matrix for the all electron and one for the pseudo part for each atomic specie
ALLOCATE(ke(ntyp))
CALL allocate_keeq(ntyp, nh, ke)
CALL allocate_paw_fockrnl(ntyp, nh, ke)
DO ns = 1,ntyp
!
ALLOCATE(k_ae(nh(ns),nh(ns),nh(ns),nh(ns)))
CALL PAW_keeq('AE', ns, k_ae)
CALL PAW_fock_onecenter('AE', ns, k_ae)
!
ALLOCATE(k_ps(nh(ns),nh(ns),nh(ns),nh(ns)))
CALL PAW_keeq('PS', ns, k_ps)
CALL PAW_fock_onecenter('PS', ns, k_ps)
!
ke(ns)%k = k_ae - k_ps
! Symmetrize wrt the on-site wavefunctions indexes as the hartree kernel is not
! perfectly symmetrical: the asymmetry accumulates and causes S matrix to be non-positive
! definite (especially with many k-points)
DO ih = 1, nh(ns)
DO jh = 1, nh(ns)
DO oh = 1, nh(ns)
DO uh = 1, nh(ns)
!
ke(ns)%k(ih,jh,oh,uh) = 0.25_dp * ( &
k_ae(ih,jh,oh,uh)-k_ps(ih,jh,oh,uh) &
+ k_ae(oh,uh,ih,jh)-k_ps(oh,uh,ih,jh) &
+ k_ae(jh,ih,uh,oh)-k_ps(jh,ih,uh,oh) &
+ k_ae(uh,oh,jh,ih)-k_ps(uh,oh,jh,ih) )
!
ENDDO
ENDDO
ENDDO
ENDDO
!
! DO ih = 1, nh(ns)
! DO jh = 1, nh(ns)
! DO oh = 1, nh(ns)
! DO uh = 1, nh(ns)
! !
! ke(ns)%k(ih,jh,oh,uh) = 0.25_dp * ( &
! k_ae(ih,jh,oh,uh)-k_ps(ih,jh,oh,uh) &
! + k_ae(oh,uh,ih,jh)-k_ps(oh,uh,ih,jh) &
! + k_ae(jh,ih,uh,oh)-k_ps(jh,ih,uh,oh) &
! + k_ae(uh,oh,jh,ih)-k_ps(uh,oh,jh,ih) )
! !
! ENDDO
! ENDDO
! ENDDO
! ENDDO
! !
DEALLOCATE(k_ae, k_ps)
!
ENDDO
!=----------------------------------------------------------------------------=!
END SUBROUTINE PAW_init_keeq
END SUBROUTINE PAW_init_fock_kernel
!=----------------------------------------------------------------------------=!
!
!=----------------------------------------------------------------------------=!
SUBROUTINE PAW_destroy_keeq()
SUBROUTINE PAW_clean_fock_kernel()
!=----------------------------------------------------------------------------=!
! ke_ae and ke_ps for later use
USE ions_base, ONLY : ityp, ntyp => nsp
IMPLICIT NONE
IF(.not.paw_has_init_keeq) RETURN !CALL errore("PAW_destroy_keeq", "nothing to destroy :(", 1)
paw_has_init_keeq = .false.
IF(.not.paw_has_init_paw_fockrnl) RETURN
paw_has_init_paw_fockrnl = .false.
!
! We have one matrix for the all electron and one for the pseudo part for each atomic specie
CALL deallocate_keeq(ntyp, ke)
CALL deallocate_paw_fockrnl(ntyp, ke)
DEALLOCATE(ke)
!
!=----------------------------------------------------------------------------=!
END SUBROUTINE PAW_destroy_keeq
END SUBROUTINE PAW_clean_fock_kernel
!=----------------------------------------------------------------------------=!
!
!=----------------------------------------------------------------------------=!
SUBROUTINE allocate_keeq(ntp, nh, keeq)
SUBROUTINE allocate_paw_fockrnl(ntp, nh, paw_fockrnl)
!=----------------------------------------------------------------------------=!
IMPLICIT NONE
INTEGER,INTENT(in) :: ntp
INTEGER,INTENT(in) :: nh(ntp)
TYPE(paw_keeq_type),INTENT(inout) :: keeq(ntp)
TYPE(paw_fockrnl_type),INTENT(inout) :: paw_fockrnl(ntp)
INTEGER :: i
!
DO i = 1,ntp
ALLOCATE(keeq(i)%k(nh(i),nh(i),nh(i),nh(i)))
ALLOCATE(paw_fockrnl(i)%k(nh(i),nh(i),nh(i),nh(i)))
ENDDO
RETURN
!=----------------------------------------------------------------------------=!
END SUBROUTINE allocate_keeq
END SUBROUTINE allocate_paw_fockrnl
!=----------------------------------------------------------------------------=!
!=----------------------------------------------------------------------------=!
SUBROUTINE deallocate_keeq(ntp, keeq)
SUBROUTINE deallocate_paw_fockrnl(ntp, paw_fockrnl)
!=----------------------------------------------------------------------------=!
IMPLICIT NONE
INTEGER,INTENT(in) :: ntp
TYPE(paw_keeq_type),INTENT(inout) :: keeq(ntp)
TYPE(paw_fockrnl_type),INTENT(inout) :: paw_fockrnl(ntp)
INTEGER :: i
!
DO i = 1,ntp
DEALLOCATE(keeq(i)%k)
DEALLOCATE(paw_fockrnl(i)%k)
ENDDO
!
RETURN
!=----------------------------------------------------------------------------=!
END SUBROUTINE deallocate_keeq
END SUBROUTINE deallocate_paw_fockrnl
!=----------------------------------------------------------------------------=!
!
!=----------------------------------------------------------------------------=!
SUBROUTINE PAW_keeq(what, np, keeq)
SUBROUTINE PAW_fock_onecenter(what, np, paw_fockrnl)
!=----------------------------------------------------------------------------=!
! Compute the 2-electron 4-wavefunctions integrals and i.e. the exchange integral
! between two one-center wavefunctions. Includes augmentation in the pseudo case.
! Store it in global variables ke_ae and ke_ps for later use.
! Compute the 2-electron 4-wavefunctions integrals i.e. the exchange integral
! for one atomic species (either pseudo or all-electron)
! Includes augmentation in the pseudo case.
USE constants, ONLY : e2
USE atom, ONLY : g => rgrid
USE ions_base, ONLY : nat, ityp, ntyp => nsp
@ -277,7 +280,7 @@ MODULE paw_exx
USE lsda_mod, ONLY : nspin
IMPLICIT NONE
INTEGER, INTENT(in) :: np ! atomic type
REAL(DP),INTENT(out) :: keeq(nh(np),nh(np),nh(np),nh(np))
REAL(DP),INTENT(out) :: paw_fockrnl(nh(np),nh(np),nh(np),nh(np))
CHARACTER(len=2),INTENT(in) :: what ! "AE"= all-electron or "PS"=pseudo
!
TYPE(paw_info) :: i
@ -289,15 +292,15 @@ MODULE paw_exx
REAL(DP) :: kexx, e
IF(what/="AE" .and. what /="PS") &
CALL errore("PAW_keeq", "can only do all-electron or pseudo", 1)
CALL errore("PAW_fock_onecenter", "can only do all-electron or pseudo", 1)
! Only wavefunctions on the same atom exchange, and the result only depends on the atom type
IF (.not.upf(np)%tpawp) THEN
keeq = 0._dp
paw_fockrnl = 0._dp
RETURN
ENDIF
!
CALL start_clock('PAW_keeq')
CALL start_clock('PAW_fock_onecenter')
!
i%a = -1 ! atom's index (UNUSED HERE)
i%t = np ! type of atom
@ -313,6 +316,8 @@ MODULE paw_exx
ALLOCATE(aux(i%m))
!
DO ih = 1, nh(i%t)
! It would be better to do jh = ih, nh(i%t), but it is a bit complicated
! inside rho_lm_ij we always set ih as the highest of the two indexes
DO jh = 1, nh(i%t)
!
rho_lm_ij = 0._dp
@ -333,17 +338,18 @@ MODULE paw_exx
CALL PAW_rho_lm_ij(i, oh, uh, upf(i%t)%paw%pfunc, rho_lm_ou)
ENDIF
!
! Now I have rho_ij and rho_ou, I have to compute the 4-wfcs hartree kernel keeq=K_ijou
! Now I have rho_ij and rho_ou, I have to compute the 4-wfcs hartree kernel paw_fockrnl=K_ijou
kexx = 0._dp
DO lm = 1, i%l**2
DO k = 1, i%m
aux(k) = v_lm(k,lm) * rho_lm_ou(k,lm)
ENDDO
CALL simpson(upf(i%t)%kkbeta, aux, g(i%t)%rab, e)
!CALL simpson(upf(i%t)%kkbeta, aux, g(i%t)%rab, e)
CALL simpson(i%m, aux, g(i%t)%rab, e)
kexx = kexx + e
!
ENDDO
keeq(ih,jh,oh,uh) = e2*kexx ! = K_ijlk : Eq. 33 Ref. 1
paw_fockrnl(ih,jh,oh,uh) = e2*kexx ! = K_ijlk : Eq. 33 Ref. 1
!
ENDDO
ENDDO
@ -352,14 +358,14 @@ MODULE paw_exx
!
DEALLOCATE(aux, v_lm, rho_lm_ij, rho_lm_ou)
!
CALL stop_clock('PAW_keeq')
CALL stop_clock('PAW_fock_onecenter')
!
!=----------------------------------------------------------------------------=!
END SUBROUTINE PAW_keeq
END SUBROUTINE PAW_fock_onecenter
!=----------------------------------------------------------------------------=!
!
!=----------------------------------------------------------------------------=!
SUBROUTINE PAW_rho_lm_ij(i, ih, jh, pfunc, rho_lm, aug)
SUBROUTINE PAW_rho_lm_ij(i, ih_, jh_, pfunc, rho_lm, aug)
!=----------------------------------------------------------------------------=!
! Computes the fake two-wavefunctions density i.e. phi_i(r)phi_j(r)^*,
! Represent it as spherical harmonics. Details: this is a generalized version of PAW_rho_lm.
@ -369,22 +375,37 @@ MODULE paw_exx
USE paw_variables, ONLY : paw_info
IMPLICIT NONE
TYPE(paw_info), INTENT(IN) :: i ! atom's minimal info
INTEGER,INTENT(in) :: ih, jh
INTEGER,INTENT(in) :: ih_, jh_
INTEGER :: ih, jh
REAL(DP), INTENT(IN) :: pfunc(i%m,i%b,i%b) ! psi_i * psi_j
REAL(DP), INTENT(OUT) :: rho_lm(i%m,i%l**2) ! AE charge density on rad. grid
REAL(DP), OPTIONAL,INTENT(IN) :: &
aug(i%m,i%b*(i%b+1)/2,0:2*upf(i%t)%lmax) ! augmentation functions (only for PS part)
aug(i%m,(i%b*(i%b+1))/2,0:2*upf(i%t)%lmax) ! augmentation functions (only for PS part)
INTEGER :: nb_, mb_
INTEGER :: nb, mb, &
nmb, & ! composite "triangular" index for pfunc nmb = 1,nh*(nh+1)/2
lm,lp,l ! counters for angular momentum lm = l**2+m
REAL(DP) :: pref
! initialize density
rho_lm(:,:) = 0._dp
! loop on all pfunc for this kind of pseudo
nb = indv(ih,i%t)
mb = indv(jh,i%t)
nmb = mb * (mb-1)/2 + nb ! mb has to be .ge. nb
! nb has to be less or equal mb for nmb to be compute correctly
! the matrix is symmetric, hence the order does not matter
nb_ = indv(ih_,i%t)
mb_ = indv(jh_,i%t)
IF(mb_>=nb_)THEN
ih = ih_
jh = jh_
nb = nb_
mb = mb_
ELSE
ih = jh_
jh = ih_
nb = mb_
mb = nb_
ENDIF
nmb = (mb*(mb-1))/2 + nb ! mb has to be .ge. nb
!
angular_momentum: &
DO lp = 1, lpx (nhtolm(jh,i%t), nhtolm(ih,i%t)) !lmaxq**2

4
PW/src/paw_onecenter.f90
View File

@ -1142,7 +1142,7 @@ SUBROUTINE PAW_rho_lm(i, becsum, pfunc, rho_lm, aug)
REAL(DP), INTENT(IN) :: pfunc(i%m,i%b,i%b) ! psi_i * psi_j
REAL(DP), INTENT(OUT) :: rho_lm(i%m,i%l**2,nspin_mag) ! AE charge density on rad. grid
REAL(DP), OPTIONAL,INTENT(IN) :: &
aug(i%m,i%b*(i%b+1)/2,0:2*upf(i%t)%lmax) ! augmentation functions (only for PS part)
aug(i%m,(i%b*(i%b+1))/2,0:2*upf(i%t)%lmax) ! augmentation functions (only for PS part)
REAL(DP) :: pref ! workspace (ap*becsum)
@ -1177,7 +1177,7 @@ SUBROUTINE PAW_rho_lm(i, becsum, pfunc, rho_lm, aug)
ijh = ijh+1
nb = indv(ih,i%t)
mb = indv(jh,i%t)
nmb = mb * (mb-1)/2 + nb ! mb has to be .ge. nb
nmb = (mb*(mb-1))/2 + nb ! mb has to be .ge. nb
!write(*,'(99i4)') nb,mb,nmb
IF (ABS(becsum(ijh,i%a,ispin)) < eps12) CYCLE
!

4
PW/src/us_exx.f90
View File

@ -348,8 +348,8 @@ MODULE us_exx
CALL errore('newdxx_g', 'need gamma tricks for this flag: '//flag, 2 )
IF ( gamma_only .AND. add_complex ) &
CALL errore('newdxx_g', 'gamma trick not good for this flag: '//flag, 3 )
IF ( ( add_complex .AND. .NOT. PRESENT(becphi_c) ) .OR. &
( add_real .AND. .NOT. PRESENT(becphi_r) ) .OR. &
IF ( ( add_complex .AND..NOT. PRESENT(becphi_c) ) .OR. &
( add_real .AND..NOT. PRESENT(becphi_r) ) .OR. &
( add_imaginary.AND..NOT. PRESENT(becphi_r) ) ) &
CALL errore('newdxx_g', 'called with incorrect arguments', 2 )
!