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

Various improvements to EXX by Ivan (see release-notes)

This commit is contained in:
Paolo Giannozzi 2017-12-21 17:45:37 +01:00
parent 7f7608d669
commit ed63ae87d5
11 changed files with 783 additions and 342 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
Doc/release-notes
View File

@ -1,5 +1,18 @@
New in development version:
* EXX: ACE can be projected on arbitrary number of orbitals
Slight change in SCF convergence with ACE (dexx) properly
accounting for difference between 2*<new|V(old)|new>
and <new|V(old)|new> + <old|V(new)|old> (see "fock3" term)
SCDM and localized exchange is now compatible with ecutfock
Localized exchange can deal with virtual orbitals (empty bands)
(Ivan Carnimeo)
Problems fixed in development version:
* Phonons with images wasn't working in v.6.2.1 with old PP files
containing "&" in the header section (Dec.20, 2017)
* PWscf in "manypw" mode wasn't working due to a bad IF (Dec.12, 2017)
Incompatible changes in development version:

3
Modules/input_parameters.f90
View File

@ -442,6 +442,7 @@ MODULE input_parameters
REAL(DP) :: ecutfock = -1.d0
! variables used in Lin Lin's ACE and SCDM
REAL(DP) :: localization_thr = 0.0_dp, scdmden=0.10d0, scdmgrd=0.20d0
INTEGER :: n_proj = 0
LOGICAL :: scdm=.FALSE.
LOGICAL :: ace=.TRUE.
@ -629,7 +630,7 @@ MODULE input_parameters
edir, emaxpos, eopreg, eamp, smearing, starting_ns_eigenvalue, &
U_projection_type, input_dft, la2F, assume_isolated, &
nqx1, nqx2, nqx3, ecutfock, localization_thr, scdm, ace, &
scdmden, scdmgrd, &
scdmden, scdmgrd, n_proj, &
exxdiv_treatment, x_gamma_extrapolation, yukawa, ecutvcut, &
exx_fraction, screening_parameter, ref_alat, &
noncolin, lspinorb, starting_spin_angle, lambda, angle1, angle2, &

2
Modules/read_namelists.f90
View File

@ -220,6 +220,7 @@ MODULE read_namelists_module
! ... EXX
!
ace=.TRUE.
n_proj = 0
localization_thr = 0.0_dp
scdm=.FALSE.
scdmden=0.10d0
@ -805,6 +806,7 @@ MODULE read_namelists_module
CALL mp_bcast( scdm, ionode_id, intra_image_comm )
CALL mp_bcast( scdmden, ionode_id, intra_image_comm )
CALL mp_bcast( scdmgrd, ionode_id, intra_image_comm )
CALL mp_bcast( n_proj, ionode_id, intra_image_comm )
CALL mp_bcast( nqx1, ionode_id, intra_image_comm )
CALL mp_bcast( nqx2, ionode_id, intra_image_comm )
CALL mp_bcast( nqx3, ionode_id, intra_image_comm )

11
PW/src/electrons.f90
View File

@ -41,7 +41,7 @@ SUBROUTINE electrons()
lambda, report
USE uspp, ONLY : okvan
USE exx, ONLY : aceinit,exxinit, exxenergy2, exxenergy, exxbuff, &
fock0, fock1, fock2, dexx, use_ace, local_thr
fock0, fock1, fock2, fock3, dexx, use_ace, local_thr
USE funct, ONLY : dft_is_hybrid, exx_is_active
USE control_flags, ONLY : adapt_thr, tr2_init, tr2_multi, gamma_only
!
@ -90,6 +90,7 @@ SUBROUTINE electrons()
IF (dft_is_hybrid() .AND. adapt_thr ) tr2= tr2_init
fock0 = 0.D0
fock1 = 0.D0
fock3 = 0.D0
IF (.NOT. exx_is_active () ) fock2 = 0.D0
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -214,7 +215,7 @@ SUBROUTINE electrons()
!
CALL exxinit(DoLoc)
IF( DoLoc ) CALL localize_orbitals( )
IF (use_ace) CALL aceinit ( )
IF (use_ace) CALL aceinit ( fock3 )
!
! fock2 is the exchange energy calculated for orbitals at step n,
! using orbitals at step n in the expression of exchange
@ -231,7 +232,11 @@ SUBROUTINE electrons()
! there is some numerical problem. One such cause could be that
! the treatment of the divergence in exact exchange has failed.
!
dexx = fock1 - 0.5D0*(fock0+fock2)
IF (use_ace) THEN
dexx = 0.5D0 * ((fock1-fock0)+(fock3-fock2))
ELSE
dexx = fock1 - 0.5D0*(fock0+fock2)
END IF
!
IF ( dexx < 0.0_dp ) THEN
IF( Doloc ) THEN

207
PW/src/exx.f90
View File

@ -66,6 +66,7 @@ MODULE exx
LOGICAL :: use_ace ! true: Use Lin Lin's ACE method
! false: do not use ACE, use old algorithm instead
COMPLEX(DP), ALLOCATABLE :: xi(:,:,:) ! ACE projectors
COMPLEX(DP), ALLOCATABLE :: evc0(:,:) ! old wfc (G-space) needed to compute fock3
INTEGER :: nbndproj
LOGICAL :: domat
REAL(DP):: local_thr ! threshold for Lin Lin's SCDM localized orbitals:
@ -133,9 +134,10 @@ MODULE exx
!
! energy related variables
!
REAL(DP) :: fock0 = 0.0_DP, & ! sum <phi|Vx(phi)|phi>
fock1 = 0.0_DP, & ! sum <psi|vx(phi)|psi>
fock2 = 0.0_DP, & ! sum <psi|vx(psi)|psi>
REAL(DP) :: fock0 = 0.0_DP, & ! sum <old|Vx(old)|old>
fock1 = 0.0_DP, & ! sum <new|vx(old)|new>
fock2 = 0.0_DP, & ! sum <new|vx(new)|new>
fock3 = 0.0_DP, & ! sum <old|vx(new)|old>
dexx = 0.0_DP ! fock1 - 0.5*(fock2+fock0)
!
! custom fft grids
@ -325,6 +327,7 @@ MODULE exx
IF ( allocated(locbuff) ) DEALLOCATE(locbuff)
IF ( allocated(locmat) ) DEALLOCATE(locmat)
IF ( allocated(xi) ) DEALLOCATE(xi)
IF ( allocated(evc0) ) DEALLOCATE(evc0)
!
IF(allocated(becxx)) THEN
DO ikq = 1, nkqs
@ -968,7 +971,7 @@ MODULE exx
!
IF( DoLoc) then
IF (.not. allocated(locbuff)) ALLOCATE( locbuff(nrxxs*npol, nbnd, nks))
IF (.not. allocated(locmat)) ALLOCATE( locmat(x_nbnd_occ, x_nbnd_occ))
IF (.not. allocated(locmat)) ALLOCATE( locmat(nbnd, nbnd))
ELSE
IF (.not. allocated(exxbuff)) THEN
IF (gamma_only) THEN
@ -5287,7 +5290,7 @@ END SUBROUTINE compute_becpsi
END SUBROUTINE exxbuff_comm_gamma
!-----------------------------------------------------------------------
SUBROUTINE aceinit( )
SUBROUTINE aceinit( exex )
!
USE wvfct, ONLY : nbnd, npwx, current_k
USE klist, ONLY : nks, xk, ngk, igk_k
@ -5307,8 +5310,13 @@ SUBROUTINE aceinit( )
REAL (DP) :: ee, eexx
INTEGER :: ik, npw
TYPE(bec_type) :: becpsi
REAL (DP), OPTIONAL, INTENT(OUT) :: exex
!
nbndproj = nbnd
IF(nbndproj<x_nbnd_occ.or.nbndproj>nbnd) THEN
write(stdout,'(3(A,I4))') ' occ = ', x_nbnd_occ, ' proj = ', nbndproj, ' tot = ', nbnd
CALL errore('aceinit','n_proj must be between occ and tot.',1)
END IF
IF (.not. allocated(xi)) ALLOCATE( xi(npwx*npol,nbndproj,nks) )
IF ( okvan ) CALL allocate_bec_type( nkb, nbnd, becpsi)
eexx = 0.0d0
@ -5331,13 +5339,14 @@ SUBROUTINE aceinit( )
ENDDO
CALL mp_sum( eexx, inter_pool_comm)
! WRITE(stdout,'(/,5X,"ACE energy",f15.8)') eexx
IF(present(exex)) exex = eexx
IF ( okvan ) CALL deallocate_bec_type(becpsi)
domat = .false.
END SUBROUTINE aceinit
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE aceinit_gamma(nnpw,nbnd,phi,xitmp,becpsi,exxe)
USE becmod, ONLY : bec_type
USE wvfct, ONLY : current_k
USE wvfct, ONLY : current_k, npwx
USE mp, ONLY : mp_stop
!
! compute xi(npw,nbndproj) for the ACE method
@ -5354,8 +5363,10 @@ IMPLICIT NONE
CALL start_clock( 'aceinit' )
IF(nbndproj>nbnd) CALL errore('aceinit','nbndproj greater than nbnd.',1)
IF(nbndproj<=0) CALL errore('aceinit','nbndproj le 0.',1)
IF (.not. allocated(evc0)) then
ALLOCATE( evc0(npwx*npol,nbndproj) )
evc0 = (Zero,Zero)
END IF
nrxxs= exx_fft%dfftt%nnr * npol
@ -5365,11 +5376,8 @@ IMPLICIT NONE
!
IF( local_thr.gt.0.0d0 ) then
CALL vexx_loc(nnpw, nbndproj, xitmp, mexx)
Call MatSymm('S','L',mexx,nbndproj)
CALL aceupdate(nbndproj,nnpw,xitmp,mexx)
domat0 = domat
domat = .true.
CALL vexxace_gamma(nnpw,nbndproj,phi,exxe)
domat = domat0
ELSE
! |xi> = Vx[phi]|phi>
CALL vexx(nnpw, nnpw, nbndproj, phi, xitmp, becpsi)
@ -5380,6 +5388,12 @@ IMPLICIT NONE
END IF
DEALLOCATE( mexx )
domat0 = domat
domat = .true.
CALL vexxace_gamma(nnpw,nbndproj,evc0,exxe)
evc0 = phi
domat = domat0
CALL stop_clock( 'aceinit' )
END SUBROUTINE aceinit_gamma
@ -5590,11 +5604,11 @@ USE mp_bands, ONLY : intra_bgrp_comm, me_bgrp, nproc_bgrp
implicit none
!
! Exact exchange with SCDM orbitals.
! - Vx|phi> = Vx|psi> <psi|Vx|psi>^(-1) <psi|Vx|phi>
! Vx|phi> = Vx|psi> <psi|Vx|psi>^(-1) <psi|Vx|phi>
! locmat contains localization integrals
! mexx contains in output the exchange matrix
!
integer :: npw, nbnd, nrxxs, npairs
integer :: npw, nbnd, nrxxs, npairs, ntot
integer :: ig, ir, ik, ikq, iq, ibnd, jbnd, kbnd, NQR
INTEGER :: current_ik
real(DP) :: ovpairs(2), mexx(nbnd,nbnd), exxe
@ -5605,19 +5619,16 @@ implicit none
INTEGER, EXTERNAL :: global_kpoint_index
write(stdout,'(A)') '---------------------------------'
write(stdout,'(A)') 'Exact-exchange with SCDM orbitals'
write(stdout,'(A)') '---------------------------------'
if(nbnd.gt.x_nbnd_occ) CALL errore('vexx_loc', 'nbnd > x_nbnd_occ not debugged yet',1)
write(stdout,'(5X,A)') ' '
write(stdout,'(5X,A)') 'Exact-exchange with localized orbitals'
CALL start_clock('vexxloc')
write(stdout,'(A,f24.12)') 'local_thr =', local_thr
write(stdout,'(7X,A,f24.12)') 'local_thr =', local_thr
nrxxs= exx_fft%dfftt%nnr
NQR = nrxxs*npol
! exchange projected onto localized orbitals (nbnd=occ only)
! exchange projected onto localized orbitals
ALLOCATE( fac(exx_fft%ngmt) )
ALLOCATE( rhoc(nrxxs), vc(NQR) )
ALLOCATE( RESULT(nrxxs,nbnd) )
@ -5635,25 +5646,28 @@ implicit none
xkq = xkq_collect(:,ikq)
CALL g2_convolution(exx_fft%ngmt, exx_fft%gt, xkp, xkq, fac)
RESULT = (0.0d0, 0.0d0)
DO ibnd = 1,x_nbnd_occ
! write(*,'(I1,A,2I1,A)') ibnd, '(',ibnd,ibnd,')'
DO ir = 1, NQR
rhoc(ir) = locbuff(ir,ibnd,ikq) * locbuff(ir,ibnd,ikq) / omega
ENDDO
CALL fwfft ('Custom', rhoc, exx_fft%dfftt)
vc=(0.0d0, 0.0d0)
DO ig = 1, exx_fft%ngmt
vc(exx_fft%nlt(ig)) = fac(ig) * rhoc(exx_fft%nlt(ig))
vc(exx_fft%nltm(ig)) = fac(ig) * rhoc(exx_fft%nltm(ig))
ENDDO
CALL invfft ('Custom', vc, exx_fft%dfftt)
DO ir = 1, NQR
RESULT(ir,ibnd) = RESULT(ir,ibnd) + locbuff(ir,ibnd,nkqs) * vc(ir)
ENDDO
DO ibnd = 1, nbnd
IF(x_occupation(ibnd,ikq).gt.0.0d0) THEN
DO ir = 1, NQR
rhoc(ir) = locbuff(ir,ibnd,ikq) * locbuff(ir,ibnd,ikq) / omega
ENDDO
CALL fwfft ('Custom', rhoc, exx_fft%dfftt)
vc=(0.0d0, 0.0d0)
DO ig = 1, exx_fft%ngmt
vc(exx_fft%nlt(ig)) = fac(ig) * rhoc(exx_fft%nlt(ig))
vc(exx_fft%nltm(ig)) = fac(ig) * rhoc(exx_fft%nltm(ig))
ENDDO
CALL invfft ('Custom', vc, exx_fft%dfftt)
DO ir = 1, NQR
RESULT(ir,ibnd) = RESULT(ir,ibnd) + locbuff(ir,ibnd,nkqs) * vc(ir)
ENDDO
END IF
DO kbnd = 1, ibnd-1
ovpairs(2) = ovpairs(2) + locmat(ibnd,kbnd)
IF(locmat(ibnd,kbnd).gt.local_thr) then
IF((locmat(ibnd,kbnd).gt.local_thr).and. &
((x_occupation(ibnd,ikq).gt.0.0d0).or.(x_occupation(kbnd,ikq).gt.0.0d0))) then
! write(stdout,'(2I4,3f12.6,A)') ibnd, kbnd, x_occupation(ibnd,ikq), x_occupation(kbnd,ikq), locmat(ibnd,kbnd), ' IN '
ovpairs(1) = ovpairs(1) + locmat(ibnd,kbnd)
DO ir = 1, NQR
rhoc(ir) = locbuff(ir,ibnd,ikq) * locbuff(ir,kbnd,ikq) / omega
@ -5667,13 +5681,14 @@ implicit none
ENDDO
CALL invfft ('Custom', vc, exx_fft%dfftt)
DO ir = 1, NQR
RESULT(ir,kbnd) = RESULT(ir,kbnd) + locbuff(ir,ibnd,nkqs) * vc(ir)
RESULT(ir,kbnd) = RESULT(ir,kbnd) + x_occupation(ibnd,ikq) * locbuff(ir,ibnd,nkqs) * vc(ir)
ENDDO
DO ir = 1, NQR
RESULT(ir,ibnd) = RESULT(ir,ibnd) + locbuff(ir,kbnd,nkqs) * vc(ir)
RESULT(ir,ibnd) = RESULT(ir,ibnd) + x_occupation(kbnd,ikq) * locbuff(ir,kbnd,nkqs) * vc(ir)
ENDDO
! write(stdout,'(2I4,2f12.6,I10,2f12.6)') ibnd, kbnd, mexx(ibnd,kbnd), ovpairs(1)
end if
! ELSE
! write(stdout,'(2I4,3f12.6,A)') ibnd, kbnd, x_occupation(ibnd,ikq), x_occupation(kbnd,ikq), locmat(ibnd,kbnd), ' OUT '
END IF
ENDDO
ENDDO
@ -5702,14 +5717,114 @@ implicit none
CALL matcalc('M1-',.true.,0,npw,nbnd,nbnd,RESULT,hpsi,mexx,exxe)
deallocate( RESULT )
write(stdout,'(2(A,I12),A,f12.2)') 'Pairs(full): ', (nbnd*nbnd-nbnd)/2, &
' Pairs(red): ', npairs, ' % ', dble(npairs)/dble((nbnd*nbnd-nbnd)/2 )*100.0d0
write(stdout,'(3(A,f12.6))') 'OvPairs(included): ', ovpairs(1), &
' OvPairs(tot): ', ovpairs(2), ' OvPairs(%): ', ovpairs(1)/ovpairs(2)*100.0d0
ntot = x_nbnd_occ * (x_nbnd_occ-1)/2 + x_nbnd_occ * (nbnd-x_nbnd_occ)
write(stdout,'(7X,2(A,I12),A,f12.2)') ' Pairs(full): ', ntot, ' Pairs(included): ', npairs, ' Pairs(%): ', dble(npairs)/dble(ntot)*100.0d0
write(stdout,'(7X,3(A,f12.6))') 'OvPairs(full): ', ovpairs(2), ' OvPairs(included): ', ovpairs(1), ' OvPairs(%): ', ovpairs(1)/ovpairs(2)*100.0d0
CALL stop_clock('vexxloc')
END SUBROUTINE vexx_loc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE compute_density(DoPrint,Shift,CenterPBC,SpreadPBC,Overlap,PsiI,PsiJ,NQR,ibnd,jbnd,&
RRhoIJ,CRhoIJ)
USE constants, ONLY : pi, bohr_radius_angs
USE cell_base, ONLY : alat, omega
USE mp, ONLY : mp_sum
USE mp_bands, ONLY : intra_bgrp_comm
!
! Manipulate density: get pair density, center, spread, absolute overlap
! DoPrint : ... whether to print or not the quantities
! Shift : ... .false. refer the centers to the cell -L/2 ... +L/2
! .true. shift the centers to the cell 0 ... L (presumably the one given in input)
!
IMPLICIT NONE
INTEGER, INTENT(IN) :: NQR, ibnd, jbnd
REAL(DP), INTENT(IN) :: PsiI(NQR), PsiJ(NQR)
LOGICAL, INTENT(IN) :: DoPrint
LOGICAL, INTENT(IN) :: Shift ! .false. Centers with respect to the minimum image cell convention
! .true. Centers shifted to the input cell
REAL(DP), INTENT(OUT) :: Overlap, CenterPBC(3), SpreadPBC(3)
REAL(DP), INTENT(OUT), OPTIONAL :: RRhoIJ(NQR)
COMPLEX(DP), INTENT(OUT), OPTIONAL :: CRhoIJ(NQR)
REAL(DP) :: lenx, leny, lenz, vol, rbuff, TotSpread
INTEGER :: nr1, nr2, nr3
INTEGER :: ir, i, j, k , idx, j0, k0
COMPLEX(DP) :: cbuff(3)
REAL(DP), PARAMETER :: Zero=0.0d0, One=1.0d0, Two=2.0d0
nr1 = exx_fft%dfftt%nr1x
nr2 = exx_fft%dfftt%nr2x
nr3 = exx_fft%dfftt%nr3x
lenx = alat/dble(nr1)
leny = alat/dble(nr2)
lenz = alat/dble(nr3)
vol = lenx*leny*lenz
CenterPBC = Zero
SpreadPBC = Zero
Overlap = Zero
cbuff = (Zero, Zero)
rbuff = Zero
j0 = exx_fft%dfftt%my_i0r2p ; k0 = exx_fft%dfftt%my_i0r3p
DO ir = 1, exx_fft%dfftt%nr1x*exx_fft%dfftt%my_nr2p*exx_fft%dfftt%my_nr3p
!
! ... three dimensional indexes
!
idx = ir -1
k = idx / (exx_fft%dfftt%nr1x*exx_fft%dfftt%my_nr2p)
idx = idx - (exx_fft%dfftt%nr1x*exx_fft%dfftt%my_nr2p)*k
k = k + k0
IF ( k .GE. exx_fft%dfftt%nr3 ) CYCLE
j = idx / exx_fft%dfftt%nr1x
idx = idx - exx_fft%dfftt%nr1x * j
j = j + j0
IF ( j .GE. exx_fft%dfftt%nr2 ) CYCLE
i = idx
IF ( i .GE. exx_fft%dfftt%nr1 ) CYCLE
!
rbuff = PsiI(ir) * PsiJ(ir) / omega
IF(present(RRhoIJ)) RRhoIJ(ir) = rbuff
IF(present(CRhoIJ)) CRhoIJ(ir) = (One,Zero) * rbuff + (Zero, One) * Zero
Overlap = Overlap + abs(rbuff)*vol
cbuff(1) = cbuff(1) + rbuff*exp((Zero,One)*Two*pi*DBLE(i)/DBLE(nr1))*vol
cbuff(2) = cbuff(2) + rbuff*exp((Zero,One)*Two*pi*DBLE(j)/DBLE(nr2))*vol
cbuff(3) = cbuff(3) + rbuff*exp((Zero,One)*Two*pi*DBLE(k)/DBLE(nr3))*vol
ENDDO
call mp_sum(cbuff,intra_bgrp_comm)
call mp_sum(Overlap,intra_bgrp_comm)
CenterPBC(1) = alat/Two/pi*aimag(log(cbuff(1)))
CenterPBC(2) = alat/Two/pi*aimag(log(cbuff(2)))
CenterPBC(3) = alat/Two/pi*aimag(log(cbuff(3)))
IF(Shift) then
if(CenterPBC(1).lt.Zero) CenterPBC(1) = CenterPBC(1) + alat
if(CenterPBC(2).lt.Zero) CenterPBC(2) = CenterPBC(2) + alat
if(CenterPBC(3).lt.Zero) CenterPBC(3) = CenterPBC(3) + alat
END IF
rbuff = dble(cbuff(1))**2 + aimag(cbuff(1))**2
SpreadPBC(1) = -(alat/Two/pi)**2 * dlog(rbuff)
rbuff = dble(cbuff(2))**2 + aimag(cbuff(2))**2
SpreadPBC(2) = -(alat/Two/pi)**2 * dlog(rbuff)
rbuff = dble(cbuff(3))**2 + aimag(cbuff(3))**2
SpreadPBC(3) = -(alat/Two/pi)**2 * dlog(rbuff)
TotSpread = (SpreadPBC(1) + SpreadPBC(2) + SpreadPBC(3))*bohr_radius_angs**2
IF(DoPrint) then
write(stdout,'(A,2I4)') 'MOs: ', ibnd, jbnd
write(stdout,'(A,10f12.6)') 'Absolute Overlap: ', Overlap
write(stdout,'(A,10f12.6)') 'Center(PBC)[A]: ', CenterPBC(1)*bohr_radius_angs, CenterPBC(2)*bohr_radius_angs, CenterPBC(3)*bohr_radius_angs
write(stdout,'(A,10f12.6)') 'Spread [A**2]: ', SpreadPBC(1)*bohr_radius_angs**2, SpreadPBC(2)*bohr_radius_angs**2, SpreadPBC(3)*bohr_radius_angs**2
write(stdout,'(A,10f12.6)') 'Total Spread [A**2]: ', TotSpread
END IF
IF(TotSpread.lt.Zero) Call errore('compute_density','Negative spread found',1)
END SUBROUTINE compute_density
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
END MODULE exx
!-----------------------------------------------------------------------

64
PW/src/gradcorr.f90
View File

@ -419,6 +419,70 @@ SUBROUTINE gradient( nrxx, a, ngm, g, nl, ga )
RETURN
!
END SUBROUTINE gradient
!----------------------------------------------------------------------------
SUBROUTINE exx_gradient( nrxx, a, ngm, g, nl, ga )
!----------------------------------------------------------------------------
!
! ... Calculates ga = \grad a in R-space (a is also in R-space, exx grid)
!
USE constants, ONLY : tpi
USE cell_base, ONLY : tpiba
USE kinds, ONLY : DP
USE control_flags, ONLY : gamma_only
USE exx, ONLY : exx_fft
USE fft_interfaces,ONLY : fwfft, invfft
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: nrxx
INTEGER, INTENT(IN) :: ngm, nl(ngm)
REAL(DP), INTENT(IN) :: a(nrxx), g(3,ngm)
REAL(DP), INTENT(OUT) :: ga(3,nrxx)
!
INTEGER :: ipol
COMPLEX(DP), ALLOCATABLE :: aux(:), gaux(:)
!
!
ALLOCATE( aux( nrxx ) )
ALLOCATE( gaux( nrxx ) )
!
aux = CMPLX( a(:), 0.D0 ,kind=DP)
!
! ... bring a(r) to G-space, a(G) ...
!
CALL fwfft ('Custom', aux, exx_fft%dfftt)
!
! ... multiply by (iG) to get (\grad_ipol a)(G) ...
!
DO ipol = 1, 3
!
gaux(:) = CMPLX(0.d0,0.d0, kind=dp)
!
gaux(nl(:)) = g(ipol,:) * &
CMPLX( -AIMAG( aux(nl(:)) ), REAL( aux(nl(:)) ) ,kind=DP)
!
IF ( gamma_only ) THEN
!
gaux(exx_fft%nltm(:)) = CMPLX( REAL( gaux(nl(:)) ), -AIMAG( gaux(nl(:)) ) ,kind=DP)
!
END IF
!
! ... bring back to R-space, (\grad_ipol a)(r) ...
!
CALL invfft ('Custom', gaux, exx_fft%dfftt)
!
! ...and add the factor 2\pi/a missing in the definition of G
!
ga(ipol,:) = tpiba * DBLE( gaux(:) )
!
END DO
!
DEALLOCATE( gaux )
DEALLOCATE( aux )
!
RETURN
!
END SUBROUTINE exx_gradient
!
!----------------------------------------------------------------------------
SUBROUTINE grad_dot( nrxx, a, ngm, g, nl, alat, da )

11
PW/src/input.f90
View File

@ -139,7 +139,7 @@ SUBROUTINE iosys()
yukawa_ => yukawa, &
ecutvcut_ => ecutvcut, &
ecutfock_ => ecutfock, &
use_ace, local_thr
use_ace, nbndproj, local_thr
USE loc_scdm, ONLY : use_scdm, scdm_den, scdm_grd
!
USE lsda_mod, ONLY : nspin_ => nspin, &
@ -245,7 +245,7 @@ SUBROUTINE iosys()
exxdiv_treatment, yukawa, ecutvcut, &
exx_fraction, screening_parameter, ecutfock, &
gau_parameter, localization_thr, scdm, ace, &
scdmden, scdmgrd, &
scdmden, scdmgrd, n_proj, &
edir, emaxpos, eopreg, eamp, noncolin, lambda, &
angle1, angle2, constrained_magnetization, &
B_field, fixed_magnetization, report, lspinorb,&
@ -1601,16 +1601,17 @@ SUBROUTINE iosys()
yukawa_ = yukawa
ecutvcut_ = ecutvcut
use_ace = ace
nbndproj = n_proj
local_thr = localization_thr
use_scdm = scdm
scdm_den = scdmden
scdm_grd = scdmgrd
IF ( local_thr > 0.0_dp .AND. .NOT. gamma_only) &
CALL errore('input','localization for k-points not yet implemented',1)
CALL errore('input','localization for k-points not implemented',1)
IF ( local_thr > 0.0_dp .AND. .NOT. use_ace ) &
CALL errore('input','localization without ACE not yet implemented',1)
CALL errore('input','localization without ACE not implemented',1)
IF ( local_thr > 0.0_dp .AND. nspin > 1 ) &
CALL errore('input','spin-polarized localization not yet implemented',1)
CALL errore('input','spin-polarized localization not implemented',1)
IF ( use_scdm ) CALL errore('input','use_scdm not yet implemented',1)
!
IF(ecutfock <= 0.0_DP) THEN

80
PW/src/interpolate.f90
View File

@ -166,3 +166,83 @@ subroutine cinterpolate (v, vs, iflag)
call stop_clock ('interpolate')
return
end subroutine cinterpolate
!
subroutine exx_interpolate (v, vs, iflag)
!
! This subroutine interpolates :
! vs on the exx mesh to v on the density mesh (iflag>0)
! vs is unchanged on output
! v on the density mesh to vs on the exx mesh (iflag<=0)
! v is unchanged on output
! V and Vs are real and in real space . V and Vs may coincide
!
USE kinds, ONLY: DP
USE gvect, ONLY: nl, nlm, g
USE control_flags, ONLY: gamma_only
USE fft_base, ONLY : dfftp
USE exx, ONLY : exx_fft
USE fft_interfaces,ONLY : fwfft, invfft
!
implicit none
real(DP) :: v (dfftp%nnr), vs (exx_fft%dfftt%nnr)
! function on density mesh
! function on exx mesh
complex(DP), allocatable :: aux (:), auxs (:)
! work array on density mesh
! work array on exx mesh
integer :: iflag
! gives the direction of the interpolation
integer :: ig, ir
call start_clock ('interpolate')
if (iflag <= 0) then
!
! from density to exx
!
allocate (aux( dfftp%nnr))
allocate (auxs(exx_fft%dfftt%nnr))
aux (:) = (1.0d0,0.0d0) * v (:)
CALL fwfft ('Dense', aux, dfftp)
auxs (:) = (0.d0, 0.d0)
do ig = 1, exx_fft%ngmt
auxs (exx_fft%nlt(ig)) = aux(nl(ig))
enddo
if (gamma_only) then
do ig = 1, exx_fft%ngmt
auxs(exx_fft%nltm(ig) ) = aux (nlm(ig) )
enddo
end if
CALL invfft ('Custom', auxs, exx_fft%dfftt)
vs (:) = real(auxs (:))
deallocate (auxs)
deallocate (aux)
else
!
! from exx to density
!
allocate (aux( dfftp%nnr))
allocate (auxs(exx_fft%dfftt%nnr))
auxs (:) = vs (:)
CALL fwfft ('Custom', auxs, exx_fft%dfftt)
aux (:) = (0.d0, 0.d0)
do ig = 1, exx_fft%ngmt
aux (nl (ig) ) = auxs (exx_fft%nlt (ig) )
enddo
if (gamma_only) then
do ig = 1, exx_fft%ngmt
aux (nlm(ig) ) = auxs (exx_fft%nltm(ig) )
enddo
end if
CALL invfft ('Dense', aux, dfftp)
v (:) = aux (:)
deallocate (auxs)
deallocate (aux)
endif
call stop_clock ('interpolate')
return
end subroutine exx_interpolate

649
PW/src/loc_scdm.f90
View File

@ -21,86 +21,347 @@ MODULE loc_scdm
SAVE
LOGICAL :: use_scdm=.false. ! if .true. enable Lin Lin's SCDM localization
! currently implemented only within ACE formalism
LOGICAL :: scdm_dipole=.false.
REAL(DP) :: scdm_den, scdm_grd
REAL(DP), PARAMETER :: Zero=0.0d0, One=1.0d0, Two=2.0d0, Three=2.0d0
CONTAINS
!
!------------------------------------------------------------------------
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE SCDM_PGG(psi, NQR, nbnd_eff)
SUBROUTINE localize_orbitals( )
!
! Driver for SCDM orbital localization
!
USE noncollin_module, ONLY : npol
USE wvfct, ONLY : nbnd
USE control_flags, ONLY : gamma_only
!
implicit none
integer :: NGrid
character(len=1) :: HowTo
if(.not.gamma_only) CALL errore('localize_orbitals', 'k-points NYI.',1)
NGrid = exx_fft%dfftt%nnr * npol
HowTo = 'G' ! How to compute the absolute overlap integrals
locmat = One
CALL measure_localization(HowTo,locbuff(1,1,1),NGrid,x_nbnd_occ,nkqs,locmat(1:x_nbnd_occ,1:x_nbnd_occ))
CALL SCDM_PGG(locbuff(1,1,1), NGrid, x_nbnd_occ)
locmat = One
CALL measure_localization(HowTo,locbuff(1,1,1),NGrid,x_nbnd_occ,nkqs,locmat(1:x_nbnd_occ,1:x_nbnd_occ))
END SUBROUTINE localize_orbitals
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE measure_localization(CFlag, orbt, NGrid, NBands, NKK, MatLoc)
USE cell_base, ONLY : alat, omega, at, bg
USE exx, ONLY : compute_density
USE constants, ONLY : bohr_radius_angs
implicit none
!
! Analyze various localization criteria for the wavefunctions in orbt
! Calculation of absolute overlap:
! CFlag = 'R' real space integral (exact but slow)
! 'G' via FFT (fast but less accurate)
!
INTEGER :: NGrid, NBands, jbnd, kbnd, NKK
REAL(DP) :: orbt(NGrid, NBands,NKK)
REAL(DP) :: loc_diag, loc_off, tmp, DistMax
REAL(DP) :: RDist(3), SpreadPBC(3), TotSpread
REAL(DP), OPTIONAL :: MatLoc(NBands,NBands)
REAL(DP), ALLOCATABLE :: CenterPBC(:,:), Mat(:,:)
CHARACTER(LEN=1) :: CFlag
REAL(DP), PARAMETER :: epss=0.0010d0
ALLOCATE( Mat(NBands,NBands), CenterPBC(3,NBands) )
IF(CFlag.eq.'R') then
Call AbsOvR(orbt, NGrid, NBands, NKK, Mat)
ELSEIF(CFlag.eq.'G') then
call AbsOvG(orbt, NGrid, NBands, NKK, Mat)
ELSE
call errore('measure_localization','Wrong CFlag',1)
END IF
loc_diag = Zero
loc_off = Zero
TotSpread = Zero
DistMax = Zero
DO jbnd = 1, NBands
loc_diag = loc_diag + Mat(jbnd,jbnd)
call compute_density(.false.,.false.,CenterPBC(1,jbnd), SpreadPBC, tmp, orbt(1,jbnd,1), orbt(1,jbnd,1), NGrid, jbnd, jbnd)
TotSpread = TotSpread + SpreadPBC(1) + SpreadPBC(2) + SpreadPBC(3)
DO kbnd = 1, jbnd - 1
loc_off = loc_off + Mat(jbnd,kbnd)
RDist(1) = (CenterPBC(1,jbnd) - CenterPBC(1,kbnd))/alat
RDist(2) = (CenterPBC(2,jbnd) - CenterPBC(2,kbnd))/alat
RDist(3) = (CenterPBC(3,jbnd) - CenterPBC(3,kbnd))/alat
CALL cryst_to_cart( 1, RDist, bg, -1 )
RDist(:) = RDist(:) - ANINT( RDist(:) )
CALL cryst_to_cart( 1, RDist, at, 1 )
tmp = alat *bohr_radius_angs* sqrt( RDist(1)**2 + RDist(2)**2 + RDist(3)**2 )
IF(DistMax.lt.tmp) DistMax = tmp
ENDDO
ENDDO
write(stdout,'(7X,A,f12.6,A)') 'Max Dist [A] = ', bohr_radius_angs*alat*sqrt(Three)/Two, ' (sqrt(3)*L/2)'
write(stdout,'(7X,A,f12.6)') 'Max Dist Found [A] =', DistMax
write(stdout,'(7X,A,f12.6,I3)') 'Total Charge =', loc_diag
write(stdout,'(7X,A,f12.6,I3)') 'Total Abs. Overlap =', loc_off
write(stdout,'(7X,A,f12.6,I3)') 'Total Spread [A**2] =', TotSpread * bohr_radius_angs**2
write(stdout,'(7X,A,f12.6,I3)') 'Aver. Spread [A**2] =', TotSpread * bohr_radius_angs**2/ dble(NBands)
IF(present(MatLoc)) MAtLoc = Mat
DEALLOCATE( CenterPBC, Mat )
END SUBROUTINE measure_localization
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE AbsOvG(orbt, NGrid, NBands, NKK, Mat)
USE fft_interfaces, ONLY : fwfft
USE wvfct, ONLY : npwx
implicit none
!
! Compute the Absolute Overlap in G-space
! (cutoff might not be accurate for the moduli of the wavefunctions)
!
INTEGER :: NGrid, NBands, jbnd, ig, NKK
REAL(DP) :: orbt(NGrid, NBands,NKK), Mat(NBands,NBands), tmp
COMPLEX(DP), ALLOCATABLE :: buffer(:), Gorbt(:,:)
call start_clock('measure')
write(stdout,'(5X,A)') ' '
write(stdout,'(5X,A)') 'Absolute Overlap calculated in G-space'
! Localized functions to G-space and Overlap matrix onto localized functions
allocate( buffer(NGrid), Gorbt(npwx,NBands) )
Mat = Zero
buffer = (Zero,Zero)
Gorbt = (Zero,Zero)
DO jbnd = 1, NBands
buffer(:) = abs(dble(orbt(:,jbnd,NKK))) + (Zero,One)*Zero
CALL fwfft( 'CustomWave' , buffer, exx_fft%dfftt )
DO ig = 1, npwx
Gorbt(ig,jbnd) = buffer(exx_fft%nlt(ig))
ENDDO
ENDDO
CALL matcalc('Coeff-',.false.,0,npwx,NBands,NBands,Gorbt,Gorbt,Mat,tmp)
deallocate ( buffer, Gorbt )
call stop_clock('measure')
END SUBROUTINE AbsOvG
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE AbsOvR(orbt, NGrid, NBands, NKK, Mat)
USE mp, ONLY : mp_sum
USE mp_bands, ONLY : intra_bgrp_comm
implicit none
!
! Compute the Absolute Overlap in R-space
! (Exact but slow)
!
INTEGER :: NGrid, NBands, nxxs, ir, jbnd, kbnd, NKK
REAL(DP) :: orbt(NGrid, NBands, NKK), Mat(NBands,NBands)
REAL(DP) :: cost, tmp
call start_clock('measure')
write(stdout,'(5X,A)') ' '
write(stdout,'(5X,A)') 'Absolute Overlap calculated in R-space'
nxxs = exx_fft%dfftt%nr1x *exx_fft%dfftt%nr2x *exx_fft%dfftt%nr3x
cost = One/dble(nxxs)
Mat = Zero
DO jbnd = 1, NBands
Mat(jbnd,jbnd) = Mat(jbnd,jbnd) + cost * sum( abs(orbt(:,jbnd,NKK)) * abs(orbt(:,jbnd,NKK)))
DO kbnd = 1, jbnd - 1
tmp = cost * sum( abs(orbt(:,jbnd,NKK)) * abs(orbt(:,kbnd,NKK)) )
Mat(jbnd,kbnd) = Mat(jbnd,kbnd) + tmp
Mat(kbnd,jbnd) = Mat(kbnd,jbnd) + tmp
ENDDO
ENDDO
call mp_sum(mat,intra_bgrp_comm)
call stop_clock('measure')
END SUBROUTINE AbsOvR
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE SCDM_PGG(psi, NGrid, NBands)
USE mp_bands, ONLY : nproc_bgrp
!
! density matrix localization (I/O in psi)
!
IMPLICIT NONE
INTEGER, INTENT(IN) :: NGrid, NBands
REAL(DP), INTENT(INOUT) :: psi(NGrid,NBands)
REAL(DP), ALLOCATABLE :: QRbuff(:,:), mat(:,:)
INTEGER, ALLOCATABLE :: pivot(:), list(:), cpu_npt(:)
REAL(DP), ALLOCATABLE :: den(:), grad_den(:,:)
INTEGER :: nptot
REAL(DP) :: ThrDen, ThrGrd
call start_clock('localization')
write(stdout,'(5X,A)') ' '
write(stdout,'(5X,A)') 'SCDM localization with prescreening'
allocate( den(exx_fft%dfftt%nnr), grad_den(3, exx_fft%dfftt%nnr) )
Call scdm_thresholds( den, grad_den, ThrDen, ThrGrd )
allocate( cpu_npt(0:nproc_bgrp-1) )
Call scdm_points( den, grad_den, ThrDen, ThrGrd, cpu_npt, nptot )
allocate( list(nptot), pivot(nptot) )
Call scdm_prescreening( NGrid, NBands, psi, den, grad_den, ThrDen, ThrGrd, &
cpu_npt, nptot, list, pivot )
deallocate( den, grad_den )
! Psi(pivot(1:NBands),:) in mat
allocate( mat(NBands,NBands) )
Call scdm_fill( nptot, NGrid, NBands, cpu_npt, pivot, list, psi, Mat)
! Pc = Psi * Psi(pivot(1:NBands),:)' in QRbuff
allocate( QRbuff(NGrid, NBands) )
QRbuff = Zero
CALL DGEMM( 'N' , 'N' , NGrid, NBands, NBands, One, psi, NGrid, mat, NBands, Zero, QRbuff, NGrid)
! Orthonormalization
! Pc(pivot(1:NBands),:) in mat
Call scdm_fill( nptot, NGrid, NBands, cpu_npt, pivot, list, QRBuff, mat)
deallocate( cpu_npt )
! Cholesky(psi)^(-1) in mat
CALL invchol(NBands,mat)
Call MatSymm('U','L',mat, NBands)
! Phi = Pc * Chol^(-1) = QRbuff * mat
psi = Zero
CALL DGEMM( 'N' , 'N' , NGrid, NBands, NBands, One, QRbuff, NGrid, mat, NBands, Zero, psi, NGrid)
deallocate( QRbuff, mat, pivot, list )
write(stdout,'(7X,A)') 'SCDM-PGG done '
call stop_clock('localization')
END SUBROUTINE SCDM_PGG
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE scdm_thresholds( den, grad_den, ThrDen, ThrGrd )
USE cell_base, ONLY : omega
USE mp, ONLY : mp_stop, mp_barrier, mp_sum
USE mp_bands, ONLY : intra_bgrp_comm, me_bgrp, nproc_bgrp
USE vdW_DF, ONLY : numerical_gradient
USE fft_base, ONLY : dfftp
USE scf, ONLY : rho
USE lsda_mod, ONLY : nspin
!
! density matrix decomposition (I/O in psi)
!
USE mp, ONLY : mp_sum
USE mp_bands, ONLY : intra_bgrp_comm
IMPLICIT NONE
integer :: NQR, nbnd_eff, lwork, INFO, i, j, n, ir, jbnd, ir_end, nnr
real(DP), allocatable :: QRbuff(:,:), tau(:), work(:), mat(:,:),mat2(:,:)
integer, allocatable :: pivot(:)
real(DP) :: charge, grad, psi(NQR,nbnd_eff)
integer :: npt, nptot, icpu, ncpu_start, ncpu_end, nxxs
integer, allocatable :: list(:), cpu_npt(:)
real(DP), allocatable :: small(:,:)
real(DP), allocatable :: den(:), grad_den(:,:)
! real(DP), parameter :: ThDen = 0.10d0, ThGrad=0.200d0
REAL(DP), INTENT(OUT) :: den(exx_fft%dfftt%nnr), grad_den(3, exx_fft%dfftt%nnr)
REAL(DP), INTENT(OUT) :: ThrDen, ThrGrd
call start_clock('localization')
write(stdout,'(A)') '--------------------------------'
write(stdout,'(A)') 'Gradient-based SCDM localization'
write(stdout,'(A)') '--------------------------------'
write(stdout,'(2(A,f12.6))') ' scdm_den = ', scdm_den, ' scdm_grd = ',scdm_grd
REAL(DP), ALLOCATABLE :: temp(:)
REAL(DP) :: charge, grad, DenAve, GrdAve
INTEGER :: ir, ir_end, nxxs, nxtot
if(exx_fft%dfftt%nnr.ne.dfftp%nnr) then
write(stdout,*) exx_fft%dfftt%nnr, dfftp%nnr
call errore( 'SCDM_PGG', 'density and orbital grids not identical',1)
end if
nnr = dfftp%nnr
! interpolate density to the exx grid
allocate( temp(dfftp%nnr))
temp(:) = rho%of_r(:,1)
IF ( nspin == 2 ) temp(:) = temp(:) + rho%of_r(:,2)
Call exx_interpolate(temp, den, -1)
deallocate( temp )
#if defined (__MPI)
ir_end = dfftp%nr1x*dfftp%my_nr2p*dfftp%my_nr3p
ir_end = exx_fft%dfftt%nr1x*exx_fft%dfftt%my_nr2p*exx_fft%dfftt%my_nr3p
#else
ir_end = nnr
ir_end = exx_fft%dfftt%nnr
#endif
nxtot = exx_fft%dfftt%nr1x *exx_fft%dfftt%nr2x *exx_fft%dfftt%nr3x
nxxs = exx_fft%dfftt%nnr
charge = Zero
DenAve = Zero
do ir = 1, ir_end
charge = charge + den(ir) * omega / dble(nxtot)
DenAve = DenAve + den(ir)
end do
call mp_sum(DenAve,intra_bgrp_comm)
call mp_sum(charge,intra_bgrp_comm)
DenAve = DenAve / dble(nxtot)
write(stdout,'(7x,A,f12.6)') 'Charge = ', charge
write(stdout,'(7x,A,f12.6)') 'DenAve = ', DenAve
ThrDen = scdm_den
! gradient on the exx grid
Call exx_gradient( nxxs, den , exx_fft%ngmt, exx_fft%gt, exx_fft%nlt, grad_den )
charge = Zero
GrdAve = Zero
do ir = 1, ir_end
grad = sqrt( grad_den(1,ir)**2 + grad_den(2,ir)**2 + grad_den(3,ir)**2 )
charge = charge + grad * omega / dble(nxtot)
GrdAve = GrdAve + grad
end do
call mp_sum(GrdAve,intra_bgrp_comm)
call mp_sum(charge,intra_bgrp_comm)
GrdAve = GrdAve / dble(nxtot)
write(stdout,'(7X,A,f12.6)') 'GradTot = ', charge
write(stdout,'(7X,A,f12.6)') 'GrdAve = ', GrdAve
ThrGrd = scdm_grd
write(stdout,'(7x,2(A,f12.6))') 'scdm_den = ', scdm_den, ' scdm_grd = ',scdm_grd
write(stdout,'(7x,2(A,f12.6))') 'ThrDen = ', ThrDen, ' ThrGrd = ',ThrGrd
END SUBROUTINE scdm_thresholds
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE scdm_fill( nptot, NGrid, NBands, CPUPts, Pivot, List, Vect, Mat)
USE mp, ONLY : mp_sum
USE mp_bands, ONLY : intra_bgrp_comm, me_bgrp, nproc_bgrp
!
! Fill the matrix Mat with the elements of Vect
! mapped by CPUPts, Pivot and List
!
IMPLICIT NONE
INTEGER, INTENT(IN) :: NBands, NGrid, nptot
INTEGER, INTENT(IN) :: CPUPts(0:nproc_bgrp-1), Pivot(nptot), List(nptot)
REAL(DP), INTENT(IN) :: Vect(NGrid, NBands)
REAL(DP), INTENT(OUT) :: Mat(NBands,NBands)
INTEGER :: i, NStart, NEnd
Mat = Zero
do i = 1, NBands
NStart = sum(CPUPts(0:me_bgrp-1))
NEnd = sum(CPUPts(0:me_bgrp))
if(Pivot(i).le.NEnd.and.Pivot(i).ge.NStart+1) Mat(:,i) = Vect(List(pivot(i)),:)
end do
call mp_sum(Mat,intra_bgrp_comm)
END SUBROUTINE scdm_fill
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE scdm_points ( den, grad_den, ThrDen, ThrGrd, cpu_npt, nptot )
USE mp, ONLY : mp_sum
USE mp_bands, ONLY : intra_bgrp_comm, me_bgrp, nproc_bgrp
!
! find the relevant points for the allocation
!
IMPLICIT NONE
INTEGER, INTENT(OUT) :: cpu_npt(0:nproc_bgrp-1), nptot
REAL(DP), INTENT(IN) :: den(exx_fft%dfftt%nnr), grad_den(3, exx_fft%dfftt%nnr)
REAL(DP), INTENT(IN) :: ThrDen, ThrGrd
INTEGER :: npt, ir, ir_end
REAL(DP) :: grad
#if defined (__MPI)
ir_end = exx_fft%dfftt%nr1x*exx_fft%dfftt%my_nr2p*exx_fft%dfftt%my_nr3p
#else
ir_end = exx_fft%dfftt%nnr
#endif
nxxs = exx_fft%dfftt%nr1x *exx_fft%dfftt%nr2x *exx_fft%dfftt%nr3x
allocate( den(nnr), grad_den(nnr, 3) )
charge = 0.0d0
den(:) = rho%of_r(:,1)
IF ( nspin == 2 ) den(:) = den(:) + rho%of_r(:,2)
do ir = 1, ir_end
charge = charge + den(ir) * omega / dble(nxxs)
end do
call mp_sum(charge,intra_bgrp_comm)
write(stdout,'(A,f12.6)') ' charge = ', charge
! find numerical gradient of the density
call numerical_gradient (den, grad_den)
charge = 0.0d0
do ir = 1, ir_end
grad = sqrt( grad_den(ir,1)**2 + grad_den(ir,2)**2 + grad_den(ir,3)**2 )
charge = charge + grad * omega / dble(nxxs)
end do
call mp_sum(charge,intra_bgrp_comm)
write(stdout,'(A,f12.6)') ' grad = ', charge
! find the relevant point for the allocation
allocate( cpu_npt(0:nproc_bgrp-1) )
npt = 0
cpu_npt(:) = 0
do ir = 1, ir_end
grad = 1.0d0
if(den(ir).gt.scdm_den) then
grad = sqrt( grad_den(ir,1)**2 + grad_den(ir,2)**2 + grad_den(ir,3)**2 )
if(grad.lt.scdm_grd) then
if(den(ir).gt.ThrDen) then
grad = sqrt( grad_den(1,ir)**2 + grad_den(2,ir)**2 + grad_den(3,ir)**2 )
if(grad.lt.ThrGrd) then
npt = npt + 1
end if
end if
@ -110,247 +371,69 @@ IMPLICIT NONE
call mp_sum(nptot,intra_bgrp_comm)
if(nptot.le.0) call errore('SCDM_PGG', 'No points prescreened. Loose the thresholds', 1)
call mp_sum(cpu_npt,intra_bgrp_comm)
write(stdout,'(2(A,I8))') ' max npt = ', maxval(cpu_npt(:)), &
' min npt = ', minval(cpu_npt(:))
write(stdout,*) ' reduced matrix, allocate: ', nptot
! write(stdout,*) ' cpu_npt = ', cpu_npt(:)
write(stdout,'(7X,2(A,I8))') 'Max npt = ', maxval(cpu_npt(:)), ' Min npt = ', minval(cpu_npt(:))
write(stdout,'(7X,2(A,I10))') 'Reduced matrix, allocate: ', nptot, ' out of ', exx_fft%dfftt%nnr
END SUBROUTINE scdm_points
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE scdm_prescreening ( NGrid, NBands, psi, den, grad_den, ThrDen, ThrGrd, cpu_npt, nptot, &
list, pivot )
USE mp, ONLY : mp_sum
USE mp_bands, ONLY : intra_bgrp_comm, me_bgrp, nproc_bgrp
!
! Get List from ThrDen and ThrGrd, and Pivot from the QRCP of small
!
IMPLICIT NONE
INTEGER, INTENT(OUT) :: list(nptot), pivot(nptot)
INTEGER, INTENT(IN) :: cpu_npt(0:nproc_bgrp-1), nptot
INTEGER, INTENT(IN) :: NGrid, NBands
REAL(DP), INTENT(IN) :: psi(NGrid,NBands)
REAL(DP), INTENT(IN) :: den(exx_fft%dfftt%nnr), grad_den(3, exx_fft%dfftt%nnr)
REAL(DP), INTENT(IN) :: ThrDen, ThrGrd
! find the map of the index
allocate( small(nbnd_eff,nptot), list(nptot) )
small = 0.0d0
INTEGER :: ir, ir_end, INFO, lwork
integer :: n, npt, ncpu_start
REAL(DP) :: grad
REAL(DP), ALLOCATABLE :: small(:,:), tau(:), work(:)
#if defined (__MPI)
ir_end = exx_fft%dfftt%nr1x*exx_fft%dfftt%my_nr2p*exx_fft%dfftt%my_nr3p
#else
ir_end = exx_fft%dfftt%nnr
#endif
! find the map of the indeces
allocate( small(NBands,nptot) )
small = Zero
list = 0
n = 0
do ir = 1, ir_end
grad = 1.0d0
if(den(ir).gt.scdm_den) then
grad = sqrt( grad_den(ir,1)**2 + grad_den(ir,2)**2 + grad_den(ir,3)**2 )
if(grad.lt.scdm_grd) then
grad = One
if(den(ir).gt.ThrDen) then
grad = sqrt( grad_den(1,ir)**2 + grad_den(2,ir)**2 + grad_den(3,ir)**2 )
if(grad.lt.ThrGrd) then
n = n + 1
ncpu_start = sum(cpu_npt(0:me_bgrp-1))
icpu = ncpu_start+n
small(:,icpu) = psi(ir,:)
list(icpu) = ir
npt = ncpu_start+n
small(:,npt) = psi(ir,:)
list(npt) = ir
end if
end if
end do
call mp_sum(small,intra_bgrp_comm)
call mp_sum(list,intra_bgrp_comm)
! perform the QRCP on the small matrix and get pivot
lwork = 4*nptot
allocate( pivot(nptot), tau(nptot), work(lwork) )
tau = 0.0d0
work = 0.0d0
allocate( tau(nptot), work(lwork) )
tau = Zero
work = Zero
pivot = 0
INFO = -1
CALL DGEQP3( nbnd_eff, nptot, small, nbnd_eff, pivot, tau, work, lwork, INFO )
do i = 1, nbnd_eff
j = list(pivot(i))
grad = sqrt( grad_den(j,1)**2 + grad_den(j,2)**2 + grad_den(j,3)**2 )
ncpu_start = 0
ncpu_end = 0
ncpu_start = sum(cpu_npt(0:me_bgrp-1))
ncpu_end = sum(cpu_npt(0:me_bgrp))
end do
deallocate( den, grad_den)
deallocate( tau, work )
deallocate( small )
CALL DGEQP3( NBands, nptot, small, NBands, pivot, tau, work, lwork, INFO )
deallocate( tau, work, small )
! Psi(pivot(1:nbnd_eff),:) in mat
lwork = 3*nbnd_eff
allocate( mat(nbnd_eff,nbnd_eff), mat2(nbnd_eff,nbnd_eff), tau(nbnd_eff), work(lwork) )
mat = 0.0d0
do i = 1, nbnd_eff
ncpu_start = 0
ncpu_end = 0
ncpu_start = sum(cpu_npt(0:me_bgrp-1))
ncpu_end = sum(cpu_npt(0:me_bgrp))
if(pivot(i).le.ncpu_end.and.pivot(i).ge.ncpu_start+1) mat(:,i) = psi(list(pivot(i)),:)
end do
call mp_sum(mat,intra_bgrp_comm)
deallocate( mat2, tau, work )
! Pc = Psi * Psi(pivot(1:nbnd_eff),:)' in QRbuff
allocate( QRbuff(NQR, nbnd_eff) )
QRbuff = 0.0d0
CALL DGEMM( 'N' , 'N' , NQR, nbnd_eff, nbnd_eff, 1.0d0, psi, NQR, mat, nbnd_eff, 0.0d0, QRbuff, NQR)
! Orthonormalization
! Pc(pivot(1:nbnd_eff),:) in mat
mat = 0.0d0
do i = 1, nbnd_eff
ncpu_start = 0
ncpu_end = 0
ncpu_start = sum(cpu_npt(0:me_bgrp-1))
ncpu_end = sum(cpu_npt(0:me_bgrp))
if(pivot(i).le.ncpu_end.and.pivot(i).ge.ncpu_start+1) mat(i,:) = QRBuff(list(pivot(i)),:)
end do
deallocate( cpu_npt )
call mp_sum(mat,intra_bgrp_comm)
! Cholesky(psi)^(-1) in mat
CALL invchol(nbnd_eff,mat)
allocate( mat2(nbnd_eff,nbnd_eff) )
mat2 = 0.0d0
do i = 1, nbnd_eff
do j = i, nbnd_eff
mat2(i,j) = mat(j,i)
end do
end do
deallocate( mat )
! Phi = Pc * Chol^(-1) = QRbuff * mat
psi = 0.0d0
CALL DGEMM( 'N' , 'N' , NQR, nbnd_eff, nbnd_eff, 1.0d0, QRbuff, NQR, mat2, nbnd_eff, 0.0d0, psi, NQR)
deallocate( QRbuff, mat2, pivot, list )
write(stdout,'(A)') ' SCDM-PGG done '
call stop_clock('localization')
END SUBROUTINE SCDM_PGG
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE localize_orbitals( )
!
! Driver for SCDM orbital localization
!
USE noncollin_module, ONLY : npol
USE wvfct, ONLY : nbnd
!
implicit none
integer :: NQR
if( nbnd > x_nbnd_occ) CALL errore('localize_orbitals', 'nbnd > x_nbnd_occ allowed',1)
NQR = exx_fft%dfftt%nnr * npol
! CALL measure_localization(locbuff(1,1,1), NQR, x_nbnd_occ, locmat)
CALL measure_localization_G(locbuff(1,1,1), NQR, x_nbnd_occ, locmat)
CALL SCDM_PGG(locbuff(1,1,1), NQR, x_nbnd_occ)
! CALL measure_localization(locbuff(1,1,1), NQR, x_nbnd_occ, locmat)
CALL measure_localization_G(locbuff(1,1,1), NQR, x_nbnd_occ, locmat)
END SUBROUTINE localize_orbitals
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE measure_localization(orbt, NGrid, NBands, MatLoc)
USE kinds, ONLY : DP
USE cell_base, ONLY : omega
USE ions_base, ONLY : nat
USE mp, ONLY : mp_sum
USE mp_bands, ONLY : intra_bgrp_comm
implicit none
integer :: NGrid, NBands, nxxs, ir, jbnd, kbnd
real(DP) :: orbt(NGrid, NBands)
real(DP) :: cost, loc_diag, loc_off, tmp
real(DP), allocatable :: Mat(:,:)
real(DP), optional :: MatLoc(NBands,NBands)
real(DP), parameter :: epss=0.0010d0
call start_clock('measure')
write(stdout,'(A)') '-------------------'
write(stdout,'(A)') 'Localization matrix'
write(stdout,'(A)') '-------------------'
nxxs = exx_fft%dfftt%nr1x *exx_fft%dfftt%nr2x *exx_fft%dfftt%nr3x
cost = 1.0d0/dble(nxxs)
loc_diag = 0.0d0
loc_off = 0.0d0
allocate( Mat(NBands,NBands) )
Mat = 0.0d0
DO ir = 1, NGrid
DO jbnd = 1, NBands
Mat(jbnd,jbnd) = Mat(jbnd,jbnd) + cost * abs(orbt(ir,jbnd)) * abs(orbt(ir,jbnd))
DO kbnd = 1, jbnd - 1
tmp = cost * abs(orbt(ir,jbnd)) * abs(orbt(ir,kbnd))
Mat(jbnd,kbnd) = Mat(jbnd,kbnd) + tmp
Mat(kbnd,jbnd) = Mat(kbnd,jbnd) + tmp
ENDDO
ENDDO
ENDDO
call mp_sum(mat,intra_bgrp_comm)
DO jbnd = 1, NBands
loc_diag = loc_diag + Mat(jbnd,jbnd)
DO kbnd = 1, jbnd - 1
loc_off = loc_off + Mat(jbnd,kbnd)
ENDDO
ENDDO
IF(present(MatLoc)) MAtLoc = Mat
deallocate( Mat )
write(stdout,'(A,f12.6,I3)') ' Total Charge =', loc_diag
write(stdout,'(A,f12.6,I3)') ' Total Localization =', loc_off
tmp = dble(NBands*(NBands-1))/2.0d0
write(stdout,'(A,f12.6,I3)') ' Localization per orbital pair =', loc_off/tmp
write(stdout,'(A,f12.6,I3)') ' Localization per unit vol =', loc_off/omega
write(stdout,'(A,f12.6,I3)') ' Localization per atom =', loc_off/dble(nat)
if(abs(loc_diag-dble(NBands)).gt.epss) Call errore('measure_localization','Orthonormality broken',1)
call stop_clock('measure')
END SUBROUTINE measure_localization
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE measure_localization_G(orbt, NGrid, NBands, MatLoc)
USE kinds, ONLY : DP
USE cell_base, ONLY : omega
USE ions_base, ONLY : nat
USE mp, ONLY : mp_sum
USE mp_bands, ONLY : intra_bgrp_comm
USE fft_interfaces, ONLY : fwfft, invfft
USE wvfct, ONLY : npwx
implicit none
integer :: NGrid, NBands, nxxs, ir, jbnd, kbnd, ig
real(DP) :: orbt(NGrid, NBands)
real(DP) :: cost, loc_diag, loc_off, tmp
real(DP), allocatable :: Mat(:,:)
complex(DP), allocatable :: buffer(:), Gorbt(:,:)
real(DP), optional :: MatLoc(NBands,NBands)
real(DP), parameter :: epss=0.000010d0
call start_clock('measure')
write(stdout,'(A)') '-----------------------'
write(stdout,'(A)') 'Localization matrix (G)'
write(stdout,'(A)') '-----------------------'
! Localized functions to G-space and exchange matrix onto localized functions
allocate( buffer(NGrid), Gorbt(npwx,NBands) )
allocate( Mat(NBands,NBands) )
Mat = 0.0d0
buffer = (0.0d0,0.0d0)
Gorbt = (0.0d0,0.0d0)
DO jbnd = 1, NBands
buffer(:) = abs(dble(orbt(:,jbnd))) + (0.0d0,1.0d0)*0.0d0
CALL fwfft( 'CustomWave' , buffer, exx_fft%dfftt )
DO ig = 1, npwx
Gorbt(ig,jbnd) = buffer(exx_fft%nlt(ig))
ENDDO
ENDDO
deallocate ( buffer )
CALL matcalc('Coeff-',.false.,0,npwx,NBands,NBands,Gorbt,Gorbt,Mat,tmp)
deallocate ( Gorbt )
loc_diag = 0.0d0
loc_off = 0.0d0
DO jbnd = 1, NBands
loc_diag = loc_diag + Mat(jbnd,jbnd)
DO kbnd = 1, jbnd - 1
loc_off = loc_off + Mat(jbnd,kbnd)
ENDDO
ENDDO
IF(present(MatLoc)) MAtLoc = Mat
deallocate( Mat )
write(stdout,'(A,f12.6,I3)') ' Total Charge =', loc_diag
write(stdout,'(A,f12.6,I3)') ' Total Localization =', loc_off
tmp = dble(NBands*(NBands-1))/2.0d0
write(stdout,'(A,f12.6,I3)') ' Localization per orbital pair =', loc_off/tmp
write(stdout,'(A,f12.6,I3)') ' Localization per unit vol =', loc_off/omega
write(stdout,'(A,f12.6,I3)') ' Localization per atom =', loc_off/dble(nat)
! if(abs(loc_diag-dble(NBands)).gt.epss) Call errore('measure_localization','Orthonormality broken',1)
call stop_clock('measure')
END SUBROUTINE measure_localization_G
END SUBROUTINE scdm_prescreening
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
END MODULE loc_scdm
!-----------------------------------------------------------------------

3
PW/src/setup.f90
View File

@ -89,7 +89,7 @@ SUBROUTINE setup()
USE qes_libs_module, ONLY : qes_reset_output, qes_reset_parallel_info, qes_reset_general_info
USE qes_types_module, ONLY : output_type, parallel_info_type, general_info_type
#endif
USE exx, ONLY : ecutfock, exx_grid_init, exx_mp_init, exx_div_check
USE exx, ONLY : ecutfock, exx_grid_init, exx_mp_init, exx_div_check, nbndproj
USE funct, ONLY : dft_is_meta, dft_is_hybrid, dft_is_gradient
USE paw_variables, ONLY : okpaw
USE fcp_variables, ONLY : lfcpopt, lfcpdyn
@ -402,6 +402,7 @@ SUBROUTINE setup()
! ... set the max number of bands used in iterative diagonalization
!
nbndx = nbnd
IF(nbndproj.eq.0) nbndproj = nbnd
IF ( isolve == 0 ) nbndx = david * nbnd
!
#if defined(__MPI)

82
PW/src/utils.f90
View File

@ -36,15 +36,16 @@ SUBROUTINE matcalc (label, DoE, PrtMat, ninner, n, m, U, V, mat, ee)
mat = 0.0_dp
CALL calbec(ninner, U, V, mat, m)
IF( PrtMat .ge.2 ) CALL matprt(string//label,n,m,mat)
IF(DoE) THEN
IF(n/=m) CALL errore('matcalc','no trace for rectangular matrix.',1)
IF( PrtMat > 1 ) CALL matprt(string//label,n,m,mat)
string = 'E-'
ee = 0.0_dp
DO i = 1,n
ee = ee + wg(i,current_k)*mat(i,i)
ENDDO
IF ( PrtMat > 0 ) WRITE(stdout,'(A,f16.8,A)') string//label, ee, ' Ry'
IF ( PrtMat .ge. 1 ) WRITE(stdout,'(A,f16.8,A)') string//label, ee, ' Ry'
ENDIF
CALL stop_clock('matcalc')
@ -100,7 +101,7 @@ SUBROUTINE matcalc_k (label, DoE, PrtMat, ik, ninner, n, m, U, V, mat, ee)
END SUBROUTINE matcalc_k
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE matprt(label,n,m,A)
SUBROUTINE MatPrt(label,n,m,A)
USE kinds, ONLY : dp
USE io_global,ONLY : stdout
IMPLICIT NONE
@ -197,3 +198,78 @@ SUBROUTINE errinfo(routine,message,INFO)
END SUBROUTINE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE MatSymm( MShape, How, Mat, n )
USE kinds, ONLY : dp
USE io_global,ONLY : stdout
!
! Symmetrize the (square) matrix Mat
! How = 'U' ... copying the upper block into the lower block
! 'L' ... copying the lower block into the upper block
! 'S' ... averaging
!
! MShape = 'U' ... return the Upper Triangular (Zeros in Lower)
! 'L' ... return the Lower Triangular (Zeros in Upper)
! 'S' ... return the Square symmetric matrix
!
IMPLICIT NONE
INTEGER :: n, i, j
REAL(DP) :: Mat(n,n)
REAL(DP), ALLOCATABLE :: MatT(:,:)
REAL(DP), PARAMETER :: Zero=0.0d0, Two=2.0d0
CHARACTER(LEN=1) :: How, MShape
ALLOCATE( MatT(n,n) )
! Properly fill the lower triangular of MatT
MatT = Zero
IF(How.eq.'L') then ! use lower
do i = 1, n
MatT(i,i) = Mat(i,i)
do j = i+1, n
MatT(j,i) = Mat(j,i)
end do
end do
ELSE IF( How.eq.'U' ) then ! use upper
do i = 1, n
MatT(i,i) = Mat(i,i)
do j = i+1, n
MatT(j,i) = Mat(i,j)
end do
end do
ELSE IF( How.eq.'S' ) then ! use average
do i = 1, n
MatT(i,i) = Mat(i,i)
do j = i+1, n
MatT(j,i) = (Mat(i,j) + Mat(j,i)) / Two
end do
end do
ELSE
Call errore('MatSymm','Wrong How in MatSymm.',1)
END IF
! Properly copy the results in Mat
Mat = Zero
IF(MShape.eq.'L') then ! return lower
Mat=MatT
ELSE IF(MShape.eq.'U') then ! return upper
do i = 1, n
Mat(i,i) = MatT(i,i)
do j = i+1, n
Mat(i,j) = MatT(j,i)
end do
end do
ELSE IF(MShape.eq.'S') then ! return square
Mat=MatT
do i = 1, n
do j = i+1, n
Mat(i,j) = MatT(j,i)
end do
end do
ELSE
Call errore('MatSymm','Wrong MShape in MatSymm.',1)
END IF
DEALLOCATE( MatT )
END SUBROUTINE MatSymm
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!