mirror of https://gitlab.com/QEF/q-e.git
Initial version of the localization algorithm for EXX by Ivan Carnimeo
git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@13780 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
parent
baf00d25a9
commit
0e3ec3449a
|
@ -177,7 +177,7 @@
|
|||
CALL exx_grid_init()
|
||||
CALL exx_div_check()
|
||||
call stop_exx()
|
||||
call exxinit
|
||||
call exxinit(.false.)
|
||||
call start_exx()
|
||||
current_k= 1
|
||||
!the following is very important
|
||||
|
@ -432,7 +432,7 @@
|
|||
!if EXX is one calculates stuff for Fock operator
|
||||
if(dft_is_hybrid()) then
|
||||
!NOT_TO_BE_INCLUDED_START
|
||||
call exxinit
|
||||
call exxinit(.false.)
|
||||
current_k= 1
|
||||
!NOT_TO_BE_INCLUDED_END
|
||||
endif
|
||||
|
|
|
@ -440,8 +440,9 @@ MODULE input_parameters
|
|||
REAL(DP) :: conv_thr_multi = 0.1_DP
|
||||
REAL(DP) :: ecutfock = -1.d0
|
||||
! variables used in Lin Lin's ACE and SCDM
|
||||
REAL(DP) :: localization_thr = 0.0_dp
|
||||
LOGICAL :: scdm=.FALSE., ace=.TRUE.
|
||||
REAL(DP) :: localization_thr = 0.0_dp, scdmden=0.10d0, scdmgrd=0.20d0
|
||||
LOGICAL :: scdm=.FALSE.
|
||||
LOGICAL :: ace=.TRUE.
|
||||
|
||||
! parameters for external electric field
|
||||
INTEGER :: edir = 0
|
||||
|
@ -616,6 +617,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, &
|
||||
exxdiv_treatment, x_gamma_extrapolation, yukawa, ecutvcut, &
|
||||
exx_fraction, screening_parameter, ref_alat, &
|
||||
noncolin, lspinorb, starting_spin_angle, lambda, angle1, angle2, &
|
||||
|
|
|
@ -219,9 +219,11 @@ MODULE read_namelists_module
|
|||
!
|
||||
! ... EXX
|
||||
!
|
||||
ace=.TRUE.
|
||||
localization_thr = 0.0_dp
|
||||
scdm=.FALSE.
|
||||
ace=.TRUE.
|
||||
scdmden=0.10d0
|
||||
scdmgrd=0.20d0
|
||||
!
|
||||
! ... electric fields
|
||||
!
|
||||
|
@ -798,9 +800,11 @@ MODULE read_namelists_module
|
|||
|
||||
! ... EXX
|
||||
|
||||
CALL mp_bcast( ace, ionode_id, intra_image_comm )
|
||||
CALL mp_bcast( localization_thr, ionode_id, intra_image_comm )
|
||||
CALL mp_bcast( scdm, ionode_id, intra_image_comm )
|
||||
CALL mp_bcast( ace, 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( nqx1, ionode_id, intra_image_comm )
|
||||
CALL mp_bcast( nqx2, ionode_id, intra_image_comm )
|
||||
CALL mp_bcast( nqx3, ionode_id, intra_image_comm )
|
||||
|
|
|
@ -185,6 +185,7 @@ s_psi.o \
|
|||
save_in_cbands.o \
|
||||
save_in_electrons.o \
|
||||
scale_h.o \
|
||||
loc_scdm.o \
|
||||
scf_mod.o \
|
||||
set_kplusq.o \
|
||||
set_kup_and_kdw.o \
|
||||
|
|
|
@ -40,8 +40,8 @@ SUBROUTINE electrons()
|
|||
USE noncollin_module, ONLY : noncolin, magtot_nc, i_cons, bfield, &
|
||||
lambda, report
|
||||
USE uspp, ONLY : okvan
|
||||
USE exx, ONLY : exxinit, exxenergy2, exxenergy, exxbuff, &
|
||||
fock0, fock1, fock2, dexx, use_ace
|
||||
USE exx, ONLY : aceinit,exxinit, exxenergy2, exxenergy, exxbuff, &
|
||||
fock0, fock1, fock2, 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
|
||||
!
|
||||
|
@ -49,6 +49,7 @@ SUBROUTINE electrons()
|
|||
USE paw_onecenter, ONLY : PAW_potential
|
||||
USE paw_symmetry, ONLY : PAW_symmetrize_ddd
|
||||
USE ions_base, ONLY : nat
|
||||
USE loc_scdm, ONLY : use_scdm, localize_orbitals
|
||||
!
|
||||
!
|
||||
IMPLICIT NONE
|
||||
|
@ -70,8 +71,10 @@ SUBROUTINE electrons()
|
|||
! when using adaptive thresholds.
|
||||
LOGICAL :: first, exst
|
||||
REAL(DP) :: etot_cmp_paw(nat,2,2)
|
||||
LOGICAL :: DoLoc
|
||||
!
|
||||
!
|
||||
DoLoc = local_thr.gt.0.0d0
|
||||
exxen = 0.0d0
|
||||
iter = 0
|
||||
first = .true.
|
||||
|
@ -103,7 +106,7 @@ SUBROUTINE electrons()
|
|||
ELSE IF ( iter < 0 .OR. iter > niter ) THEN
|
||||
iter = 0
|
||||
ELSE
|
||||
READ (iunres, *) fock0, fock1, fock2
|
||||
READ (iunres, *) exxen, fock0, fock1, fock2
|
||||
! FIXME: et and wg should be read from xml file
|
||||
READ (iunres, *) (wg(1:nbnd,ik),ik=1,nks)
|
||||
READ (iunres, *) (et(1:nbnd,ik),ik=1,nks)
|
||||
|
@ -111,11 +114,13 @@ SUBROUTINE electrons()
|
|||
! ... if restarting here, exx was already active
|
||||
! ... initialize stuff for exx
|
||||
first = .false.
|
||||
CALL exxinit()
|
||||
CALL exxinit(DoLoc)
|
||||
IF( DoLoc ) CALL localize_orbitals( )
|
||||
! FIXME: ugly hack, overwrites exxbuffer from exxinit
|
||||
CALL seqopn (iunres, 'restart_exx', 'unformatted', exst)
|
||||
IF (exst) READ (iunres, iostat=ios) exxbuff
|
||||
IF (ios /= 0) WRITE(stdout,'(5x,"Error in EXX restart!")')
|
||||
IF (use_ace) CALL aceinit ( )
|
||||
!
|
||||
CALL v_of_rho( rho, rho_core, rhog_core, &
|
||||
ehart, etxc, vtxc, eth, etotefield, charge, v)
|
||||
|
@ -150,7 +155,7 @@ SUBROUTINE electrons()
|
|||
& i6)') iter
|
||||
CALL seqopn (iunres, 'restart_e', 'formatted', exst)
|
||||
WRITE (iunres, *) iter-1, tr2, dexx
|
||||
WRITE (iunres, *) fock0, fock1, fock2
|
||||
WRITE (iunres, *) exxen, fock0, fock1, fock2
|
||||
WRITE (iunres, *) (wg(1:nbnd,ik),ik=1,nks)
|
||||
WRITE (iunres, *) (et(1:nbnd,ik),ik=1,nks)
|
||||
CLOSE (unit=iunres, status='keep')
|
||||
|
@ -172,14 +177,15 @@ SUBROUTINE electrons()
|
|||
! Activate exact exchange, set orbitals used in its calculation,
|
||||
! then calculate exchange energy (will be useful at next step)
|
||||
!
|
||||
CALL exxinit()
|
||||
CALL exxinit(DoLoc)
|
||||
IF( DoLoc ) CALL localize_orbitals( )
|
||||
IF ( use_ace) THEN
|
||||
CALL aceinit ( )
|
||||
fock2 = exxenergyace()
|
||||
ELSE
|
||||
fock2 = exxenergy2()
|
||||
ENDIF
|
||||
exxen = 0.50d0*fock2
|
||||
write (6,*) 'fock energy ',exxen
|
||||
etot = etot - etxc
|
||||
!
|
||||
! Recalculate potential because XC functional has changed,
|
||||
|
@ -206,7 +212,9 @@ SUBROUTINE electrons()
|
|||
!
|
||||
! Set new orbitals for the calculation of the exchange term
|
||||
!
|
||||
CALL exxinit()
|
||||
CALL exxinit(DoLoc)
|
||||
IF( DoLoc ) CALL localize_orbitals( )
|
||||
IF (use_ace) CALL aceinit ( )
|
||||
!
|
||||
! fock2 is the exchange energy calculated for orbitals at step n,
|
||||
! using orbitals at step n in the expression of exchange
|
||||
|
@ -224,13 +232,16 @@ SUBROUTINE electrons()
|
|||
! the treatment of the divergence in exact exchange has failed.
|
||||
!
|
||||
dexx = fock1 - 0.5D0*(fock0+fock2)
|
||||
IF ( dexx < 0d0 ) THEN
|
||||
! WRITE(stdout,'(5x,a,1e12.3)') "BEWARE: negative dexx:", dexx
|
||||
! dexx = ABS(dexx)
|
||||
WRITE( stdout, * ) "dexx:", dexx
|
||||
CALL errore( 'electrons', 'dexx is negative! &
|
||||
& Check that exxdiv_treatment is appropriate for the system,&
|
||||
& or ecutfock may be too low', 1 )
|
||||
!
|
||||
IF ( dexx < 0.0_dp ) THEN
|
||||
IF( Doloc ) THEN
|
||||
WRITE(stdout,'(5x,a,1e12.3)') "BEWARE: negative dexx:", dexx
|
||||
dexx = ABS ( dexx )
|
||||
ELSE
|
||||
CALL errore( 'electrons', 'dexx is negative! &
|
||||
& Check that exxdiv_treatment is appropriate for the system,&
|
||||
& or ecutfock may be too low', 1 )
|
||||
ENDIF
|
||||
ENDIF
|
||||
!
|
||||
! remove the estimate exchange energy exxen used in the inner SCF
|
||||
|
@ -273,7 +284,7 @@ SUBROUTINE electrons()
|
|||
conv_elec=.FALSE.
|
||||
CALL seqopn (iunres, 'restart_e', 'formatted', exst)
|
||||
WRITE (iunres, *) iter, tr2, dexx
|
||||
WRITE (iunres, *) fock0, fock1, fock2
|
||||
WRITE (iunres, *) exxen, fock0, fock1, fock2
|
||||
! FIXME: et and wg are written to xml file
|
||||
WRITE (iunres, *) (wg(1:nbnd,ik),ik=1,nks)
|
||||
WRITE (iunres, *) (et(1:nbnd,ik),ik=1,nks)
|
||||
|
@ -372,6 +383,9 @@ SUBROUTINE electrons_scf ( printout, exxen )
|
|||
USE wrappers, ONLY : memstat
|
||||
!
|
||||
USE plugin_variables, ONLY : plugin_etot
|
||||
USE exx, ONLY : exxinit
|
||||
USE loc_scdm, ONLY : use_scdm, localize_orbitals
|
||||
USE funct, ONLY : dft_is_hybrid, stop_exx
|
||||
!
|
||||
IMPLICIT NONE
|
||||
!
|
||||
|
@ -407,6 +421,7 @@ SUBROUTINE electrons_scf ( printout, exxen )
|
|||
!
|
||||
REAL(DP), EXTERNAL :: ewald, get_clock
|
||||
REAL(DP) :: etot_cmp_paw(nat,2,2)
|
||||
!
|
||||
!
|
||||
iter = 0
|
||||
dr2 = 0.0_dp
|
||||
|
|
385
PW/src/exx.f90
385
PW/src/exx.f90
|
@ -32,6 +32,7 @@ MODULE exx
|
|||
COMPLEX(DP), ALLOCATABLE :: psi_exx(:,:), hpsi_exx(:,:)
|
||||
COMPLEX(DP), ALLOCATABLE :: evc_exx(:,:), psic_exx(:)
|
||||
INTEGER :: lda_original, n_original
|
||||
integer :: ibnd_buff_start, ibnd_buff_end
|
||||
!
|
||||
! general purpose vars
|
||||
!
|
||||
|
@ -56,19 +57,19 @@ MODULE exx
|
|||
INTEGER :: ibnd_end = 0 ! ending band index used in bgrp parallelization
|
||||
|
||||
COMPLEX(DP), ALLOCATABLE :: exxbuff(:,:,:)
|
||||
! temporary buffer for wfc storage
|
||||
! temporary (complex) buffer for wfc storage
|
||||
REAL(DP), ALLOCATABLE :: locbuff(:,:,:)
|
||||
! temporary (real) buffer for wfc storage
|
||||
!
|
||||
LOGICAL :: use_ace=.true. ! true: Use Lin Lin's ACE method
|
||||
! false: do not use ACE, use old algorithm instead
|
||||
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
|
||||
INTEGER :: nbndproj
|
||||
LOGICAL :: domat
|
||||
REAL(DP):: local_thr ! threshold for Lin Lin's SCDM localized orbitals:
|
||||
! discard contribution to V_x if overlap between
|
||||
! localized orbitals is smaller than "local_thr"
|
||||
!
|
||||
LOGICAL :: use_scdm=.false. ! if .true. enable Lin Lin's SCDM localization
|
||||
! currently implemented only within ACE formalism
|
||||
REAL(DP):: local_thr ! threshold for Lin Lin's SCDM localized orbitals:
|
||||
! discard contribution to V_x if overlap between
|
||||
! localized orbitals is smaller than "local_thr"
|
||||
!
|
||||
#if defined(__USE_INTEL_HBM_DIRECTIVES)
|
||||
!DIR$ ATTRIBUTES FASTMEM :: exxbuff
|
||||
|
@ -294,7 +295,7 @@ MODULE exx
|
|||
IF ( allocated(x_occupation) ) DEALLOCATE(x_occupation)
|
||||
IF ( allocated(xkq_collect) ) DEALLOCATE(xkq_collect)
|
||||
IF ( allocated(exxbuff) ) DEALLOCATE(exxbuff)
|
||||
!civn
|
||||
IF ( allocated(locbuff) ) DEALLOCATE(locbuff)
|
||||
IF ( allocated(xi) ) DEALLOCATE(xi)
|
||||
!
|
||||
IF(allocated(becxx)) THEN
|
||||
|
@ -752,20 +753,26 @@ MODULE exx
|
|||
exxalfa = get_exx_fraction()
|
||||
CALL start_exx()
|
||||
CALL weights()
|
||||
CALL exxinit()
|
||||
IF(local_thr.gt.0.0d0) Call errore('exx_restart','SCDM with restart NYI',1)
|
||||
CALL exxinit(.false.)
|
||||
IF ( use_ace) CALL aceinit ( )
|
||||
fock0 = exxenergy2()
|
||||
!
|
||||
RETURN
|
||||
!------------------------------------------------------------------------
|
||||
END SUBROUTINE exx_restart
|
||||
!------------------------------------------------------------------------
|
||||
!
|
||||
!------------------------------------------------------------------------
|
||||
SUBROUTINE exxinit()
|
||||
SUBROUTINE exxinit(DoLoc)
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
! This SUBROUTINE is run before the first H_psi() of each iteration.
|
||||
! It saves the wavefunctions for the right density matrix, in real space
|
||||
!
|
||||
! DoLoc = .true. ... Real Array exxbuff(ir, nbnd, nkqs)
|
||||
! .false. ... Complex Array exxbuff(ir, nbnd/2, nkqs)
|
||||
!
|
||||
USE wavefunctions_module, ONLY : evc, psic
|
||||
USE io_files, ONLY : nwordwfc, iunwfc_exx
|
||||
USE buffers, ONLY : get_buffer
|
||||
|
@ -791,7 +798,7 @@ MODULE exx
|
|||
IMPLICIT NONE
|
||||
INTEGER :: ik,ibnd, i, j, k, ir, isym, ikq, ig
|
||||
INTEGER :: h_ibnd
|
||||
INTEGER :: ibnd_loop_start, ibnd_buff_start, ibnd_buff_end
|
||||
INTEGER :: ibnd_loop_start
|
||||
INTEGER :: ipol, jpol
|
||||
REAL(dp), ALLOCATABLE :: occ(:,:)
|
||||
COMPLEX(DP),ALLOCATABLE :: temppsic(:)
|
||||
|
@ -811,8 +818,11 @@ MODULE exx
|
|||
INTEGER, EXTERNAL :: global_kpoint_index
|
||||
INTEGER :: ibnd_start_new, ibnd_end_new
|
||||
INTEGER :: ibnd_exx, evc_offset
|
||||
LOGICAL :: DoLoc
|
||||
!
|
||||
CALL start_clock ('exxinit')
|
||||
write(stdout,'(A,L)') 'exxinit: DoLoc=',DoLoc
|
||||
IF(DoLoc.and..not.gamma_only) Call errore('exxinit','SCDM with K-points NYI',1)
|
||||
!
|
||||
CALL transform_evc_to_exx(2)
|
||||
!
|
||||
|
@ -900,23 +910,38 @@ MODULE exx
|
|||
ibnd_buff_end = ibnd_end_new
|
||||
ENDIF
|
||||
!
|
||||
IF (.not. allocated(exxbuff)) THEN
|
||||
IF (gamma_only) THEN
|
||||
ALLOCATE( exxbuff(nrxxs*npol, ibnd_buff_start:ibnd_buff_end, nks))
|
||||
ELSE
|
||||
ALLOCATE( exxbuff(nrxxs*npol, ibnd_buff_start:ibnd_buff_end, nkqs))
|
||||
END IF
|
||||
IF( DoLoc) then
|
||||
IF (.not. allocated(locbuff)) ALLOCATE( locbuff(nrxxs*npol, nbnd, nks))
|
||||
ELSE
|
||||
IF (.not. allocated(exxbuff)) THEN
|
||||
IF (gamma_only) THEN
|
||||
ALLOCATE( exxbuff(nrxxs*npol, ibnd_buff_start:ibnd_buff_end, nks))
|
||||
ELSE
|
||||
ALLOCATE( exxbuff(nrxxs*npol, ibnd_buff_start:ibnd_buff_end, nkqs))
|
||||
END IF
|
||||
END IF
|
||||
END IF
|
||||
|
||||
!assign buffer
|
||||
!$omp parallel do collapse(3) default(shared) firstprivate(npol,nrxxs,nkqs,ibnd_buff_start,ibnd_buff_end) private(ir,ibnd,ikq,ipol)
|
||||
DO ikq=1,SIZE(exxbuff,3)
|
||||
DO ibnd=ibnd_buff_start,ibnd_buff_end
|
||||
DO ir=1,nrxxs*npol
|
||||
exxbuff(ir,ibnd,ikq)=(0.0_DP,0.0_DP)
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDDO
|
||||
IF(DoLoc) then
|
||||
DO ikq=1,SIZE(locbuff,3)
|
||||
DO ibnd=1, x_nbnd_occ
|
||||
DO ir=1,nrxxs*npol
|
||||
locbuff(ir,ibnd,ikq)=(0.0_DP,0.0_DP)
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDDO
|
||||
ELSE
|
||||
DO ikq=1,SIZE(exxbuff,3)
|
||||
DO ibnd=ibnd_buff_start,ibnd_buff_end
|
||||
DO ir=1,nrxxs*npol
|
||||
exxbuff(ir,ibnd,ikq)=(0.0_DP,0.0_DP)
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDDO
|
||||
END IF
|
||||
|
||||
!
|
||||
! This is parallelized over pools. Each pool computes only its k-points
|
||||
!
|
||||
|
@ -973,8 +998,12 @@ MODULE exx
|
|||
ENDIF
|
||||
|
||||
CALL invfft ('CustomWave', psic_exx, exx_fft%dfftt)
|
||||
|
||||
exxbuff(1:nrxxs,h_ibnd,ik)=psic_exx(1:nrxxs)
|
||||
IF(DoLoc) then
|
||||
locbuff(1:nrxxs,ibnd-ibnd_loop_start+evc_offset+1,ik)=Dble( psic_exx(1:nrxxs) )
|
||||
locbuff(1:nrxxs,ibnd-ibnd_loop_start+evc_offset+2,ik)=Aimag( psic_exx(1:nrxxs) )
|
||||
ELSE
|
||||
exxbuff(1:nrxxs,h_ibnd,ik)=psic_exx(1:nrxxs)
|
||||
END IF
|
||||
|
||||
ENDDO
|
||||
!
|
||||
|
@ -1163,8 +1192,6 @@ MODULE exx
|
|||
!
|
||||
CALL change_data_structure(.FALSE.)
|
||||
!
|
||||
IF ( use_ace) CALL aceinit ( )
|
||||
!
|
||||
CALL stop_clock ('exxinit')
|
||||
!
|
||||
!-----------------------------------------------------------------------
|
||||
|
@ -5205,11 +5232,6 @@ END SUBROUTINE compute_becpsi
|
|||
END SUBROUTINE exxbuff_comm_gamma
|
||||
!-----------------------------------------------------------------------
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
SUBROUTINE matprt(label,n,m,A)
|
||||
IMPLICIT NONE
|
||||
|
@ -5224,6 +5246,7 @@ IMPLICIT NONE
|
|||
DO i = 1,n
|
||||
WRITE(stdout,frmt) A(i,:)
|
||||
ENDDO
|
||||
|
||||
END SUBROUTINE
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
SUBROUTINE errinfo(routine,message,INFO)
|
||||
|
@ -5287,35 +5310,50 @@ SUBROUTINE aceinit( )
|
|||
END SUBROUTINE aceinit
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
SUBROUTINE aceinit_gamma(nnpw,nbnd,phi,xitmp,becpsi,exxe)
|
||||
USE becmod, ONLY : bec_type
|
||||
USE wvfct, ONLY : current_k
|
||||
USE becmod, ONLY : bec_type
|
||||
USE wvfct, ONLY : current_k
|
||||
USE mp, ONLY : mp_stop
|
||||
!
|
||||
! compute xi(npw,nbndproj) for the ACE method
|
||||
!
|
||||
IMPLICIT NONE
|
||||
INTEGER :: nnpw,nbnd
|
||||
INTEGER :: nnpw,nbnd, nrxxs
|
||||
COMPLEX(DP) :: phi(nnpw,nbnd)
|
||||
real(DP), ALLOCATABLE :: mexx(:,:)
|
||||
COMPLEX(DP) :: xitmp(nnpw,nbndproj)
|
||||
INTEGER :: i
|
||||
real(DP) :: exxe
|
||||
real(DP), PARAMETER :: Zero=0.0d0, One=1.0d0, Two=2.0d0, Pt5=0.50d0
|
||||
TYPE(bec_type), INTENT(in) :: becpsi
|
||||
logical :: domat0
|
||||
|
||||
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)
|
||||
|
||||
nrxxs= exx_fft%dfftt%nnr * npol
|
||||
|
||||
ALLOCATE( mexx(nbndproj,nbndproj) )
|
||||
xitmp = (Zero,Zero)
|
||||
mexx = Zero
|
||||
! |xi> = Vx[phi]|phi>
|
||||
CALL vexx(nnpw, nnpw, nbndproj, phi, xitmp, becpsi)
|
||||
! mexx = <phi|Vx[phi]|phi>
|
||||
CALL matcalc('exact',.true.,.false.,nnpw,nbndproj,nbndproj,phi,xitmp,mexx,exxe)
|
||||
! |xi> = -One * Vx[phi]|phi> * rmexx^T
|
||||
CALL aceupdate(nbndproj,nnpw,xitmp,mexx)
|
||||
!
|
||||
IF( local_thr.gt.0.0d0 ) then
|
||||
! CALL measure_localization(locbuff(1,1,1), nrxxs, nbndproj, mexx)
|
||||
CALL measure_localization_G(locbuff(1,1,1), nrxxs, nbndproj, mexx)
|
||||
CALL vexx_loc(nnpw, nbndproj, xitmp, mexx)
|
||||
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)
|
||||
! mexx = <phi|Vx[phi]|phi>
|
||||
CALL matcalc('exact',.true.,.false.,nnpw,nbndproj,nbndproj,phi,xitmp,mexx,exxe)
|
||||
! |xi> = -One * Vx[phi]|phi> * rmexx^T
|
||||
CALL aceupdate(nbndproj,nnpw,xitmp,mexx)
|
||||
END IF
|
||||
DEALLOCATE( mexx )
|
||||
|
||||
CALL stop_clock( 'aceinit' )
|
||||
|
@ -5410,7 +5448,7 @@ IMPLICIT NONE
|
|||
DO i = 1,n
|
||||
ee = ee + wg(i,current_k)*mat(i,i)
|
||||
ENDDO
|
||||
!WRITE(stdout,'(A,f16.8,A)') string//label, ee, ' Ry'
|
||||
WRITE(stdout,'(A,f16.8,A)') string//label, ee, ' Ry'
|
||||
ENDIF
|
||||
|
||||
CALL stop_clock('matcalc')
|
||||
|
@ -5577,9 +5615,6 @@ IMPLICIT NONE
|
|||
CALL stop_clock('aceupdate')
|
||||
|
||||
END SUBROUTINE
|
||||
|
||||
|
||||
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
SUBROUTINE vexxace_k(nnpw,nbnd,phi,exxe,vphi)
|
||||
USE becmod, ONLY : calbec
|
||||
|
@ -5633,6 +5668,258 @@ IMPLICIT NONE
|
|||
|
||||
END SUBROUTINE
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE vexx_loc(npw, nbnd, hpsi, mexx)
|
||||
USE noncollin_module, ONLY : npol
|
||||
USE cell_base, ONLY : omega, alat
|
||||
USE wvfct, ONLY : current_k
|
||||
USE klist, ONLY : xk, nks, nkstot
|
||||
USE fft_interfaces, ONLY : fwfft, invfft
|
||||
USE mp, ONLY : mp_stop, mp_barrier, mp_sum
|
||||
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>
|
||||
! mexx in input contains the localization integrals, in output the exchange matrix
|
||||
!
|
||||
integer :: npw, nbnd, nrxxs, npairs
|
||||
integer :: ig, ir, ik, ikq, iq, ibnd, jbnd, kbnd, NQR
|
||||
INTEGER :: current_ik
|
||||
real(DP) :: ovpairs(2), mexx(nbnd,nbnd), exxe
|
||||
complex(DP) :: hpsi(npw,nbnd)
|
||||
COMPLEX(DP),ALLOCATABLE :: rhoc(:), vc(:), RESULT(:,:)
|
||||
REAL(DP),ALLOCATABLE :: fac(:)
|
||||
REAL(DP) :: xkp(3), xkq(3)
|
||||
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)
|
||||
|
||||
CALL start_clock('vexxloc')
|
||||
|
||||
write(stdout,'(A,f24.12)') 'local_thr =', local_thr
|
||||
nrxxs= exx_fft%dfftt%nnr
|
||||
NQR = nrxxs*npol
|
||||
|
||||
! exchange projected onto localized orbitals (nbnd=occ only)
|
||||
ALLOCATE( fac(exx_fft%ngmt) )
|
||||
ALLOCATE( rhoc(nrxxs), vc(NQR) )
|
||||
ALLOCATE( RESULT(nrxxs,nbnd) )
|
||||
|
||||
current_ik = global_kpoint_index ( nkstot, current_k )
|
||||
xkp = xk(:,current_k)
|
||||
|
||||
vc=(0.0d0, 0.0d0)
|
||||
npairs = 0
|
||||
ovpairs = 0.0d0
|
||||
|
||||
DO iq=1,nqs
|
||||
ikq = index_xkq(current_ik,iq)
|
||||
ik = index_xk(ikq)
|
||||
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 kbnd = 1, ibnd-1
|
||||
ovpairs(2) = ovpairs(2) + mexx(ibnd,kbnd)
|
||||
IF(mexx(ibnd,kbnd).gt.local_thr) then
|
||||
ovpairs(1) = ovpairs(1) + mexx(ibnd,kbnd)
|
||||
DO ir = 1, NQR
|
||||
rhoc(ir) = locbuff(ir,ibnd,ikq) * locbuff(ir,kbnd,ikq) / omega
|
||||
ENDDO
|
||||
npairs = npairs + 1
|
||||
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,kbnd) = RESULT(ir,kbnd) + locbuff(ir,ibnd,nkqs) * vc(ir)
|
||||
ENDDO
|
||||
DO ir = 1, NQR
|
||||
RESULT(ir,ibnd) = RESULT(ir,ibnd) + locbuff(ir,kbnd,nkqs) * vc(ir)
|
||||
ENDDO
|
||||
! write(stdout,'(2I4,2f12.6,I10,2f12.6)') ibnd, kbnd, mexx(ibnd,kbnd), ovpairs(1)
|
||||
end if
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
DO jbnd = 1, nbnd
|
||||
CALL fwfft( 'CustomWave' , RESULT(:,jbnd), exx_fft%dfftt )
|
||||
DO ig = 1, npw
|
||||
hpsi(ig,jbnd) = hpsi(ig,jbnd) - exxalfa*RESULT(exx_fft%nlt(ig),jbnd)
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
ENDDO
|
||||
DEALLOCATE( fac, vc )
|
||||
DEALLOCATE( RESULT )
|
||||
|
||||
! Localized functions to G-space and exchange matrix onto localized functions
|
||||
allocate( RESULT(npw,nbnd) )
|
||||
RESULT = (0.0d0,0.0d0)
|
||||
DO jbnd = 1, nbnd
|
||||
rhoc(:) = dble(locbuff(:,jbnd,nkqs)) + (0.0d0,1.0d0)*0.0d0
|
||||
CALL fwfft( 'CustomWave' , rhoc, exx_fft%dfftt )
|
||||
DO ig = 1, npw
|
||||
RESULT(ig,jbnd) = rhoc(exx_fft%nlt(ig))
|
||||
ENDDO
|
||||
ENDDO
|
||||
deallocate ( rhoc )
|
||||
CALL matcalc('M1-',.true.,.false.,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, ' % ', float(npairs)/float((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
|
||||
|
||||
CALL stop_clock('vexxloc')
|
||||
|
||||
END SUBROUTINE vexx_loc
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
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/float(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 = float(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/float(nat)
|
||||
if(abs(loc_diag-float(NBands)).gt.epss) Call errore('measure_localization','Orthonormality broken',1)
|
||||
|
||||
call stop_clock('measure')
|
||||
|
||||
END SUBROUTINE
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
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.,.false.,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 = float(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/float(nat)
|
||||
! if(abs(loc_diag-float(NBands)).gt.epss) Call errore('measure_localization','Orthonormality broken',1)
|
||||
|
||||
call stop_clock('measure')
|
||||
|
||||
END SUBROUTINE
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
END MODULE exx
|
||||
!-----------------------------------------------------------------------
|
||||
|
|
|
@ -139,7 +139,8 @@ SUBROUTINE iosys()
|
|||
yukawa_ => yukawa, &
|
||||
ecutvcut_ => ecutvcut, &
|
||||
ecutfock_ => ecutfock, &
|
||||
local_thr, use_scdm, use_ace
|
||||
use_ace, local_thr
|
||||
USE loc_scdm, ONLY : use_scdm, scdm_den, scdm_grd
|
||||
!
|
||||
USE lsda_mod, ONLY : nspin_ => nspin, &
|
||||
starting_magnetization_ => starting_magnetization, &
|
||||
|
@ -211,7 +212,6 @@ SUBROUTINE iosys()
|
|||
USE read_pseudo_mod, ONLY : readpp
|
||||
|
||||
USE qmmm, ONLY : qmmm_config
|
||||
|
||||
!
|
||||
! ... CONTROL namelist
|
||||
!
|
||||
|
@ -242,7 +242,8 @@ SUBROUTINE iosys()
|
|||
x_gamma_extrapolation, nqx1, nqx2, nqx3, &
|
||||
exxdiv_treatment, yukawa, ecutvcut, &
|
||||
exx_fraction, screening_parameter, ecutfock, &
|
||||
gau_parameter, localization_thr, scdm, ace, &
|
||||
gau_parameter, localization_thr, scdm, ace, &
|
||||
scdmden, scdmgrd, &
|
||||
edir, emaxpos, eopreg, eamp, noncolin, lambda, &
|
||||
angle1, angle2, constrained_magnetization, &
|
||||
B_field, fixed_magnetization, report, lspinorb,&
|
||||
|
@ -295,9 +296,8 @@ SUBROUTINE iosys()
|
|||
! ... CARDS
|
||||
!
|
||||
USE input_parameters, ONLY : k_points, xk, wk, nk1, nk2, nk3, &
|
||||
k1, k2, k3, nkstot
|
||||
k1, k2, k3, nkstot
|
||||
USE input_parameters, ONLY : nconstr_inp, trd_ht, rd_ht, cell_units
|
||||
USE input_parameters, ONLY : deallocate_input_parameters
|
||||
!
|
||||
USE constraints_module, ONLY : init_constraint
|
||||
USE read_namelists_module, ONLY : read_namelists, sm_not_set
|
||||
|
@ -306,6 +306,7 @@ SUBROUTINE iosys()
|
|||
USE tsvdw_module, ONLY : vdw_isolated, vdw_econv_thr
|
||||
USE us, ONLY : spline_ps_ => spline_ps
|
||||
!
|
||||
USE input_parameters, ONLY : deallocate_input_parameters
|
||||
USE wyckoff, ONLY : nattot, sup_spacegroup
|
||||
USE qexsd_module, ONLY : qexsd_input_obj
|
||||
USE qes_types_module, ONLY: input_type
|
||||
|
@ -1536,13 +1537,11 @@ SUBROUTINE iosys()
|
|||
exxdiv_treatment_ = trim(exxdiv_treatment)
|
||||
yukawa_ = yukawa
|
||||
ecutvcut_ = ecutvcut
|
||||
use_ace = ace
|
||||
local_thr = localization_thr
|
||||
use_scdm = scdm
|
||||
use_ace = ace
|
||||
IF ( .NOT.scdm .AND. localization_thr > 0.0_dp ) &
|
||||
CALL infomsg ('iosys', 'localization threshold needs SCDM')
|
||||
IF ( scdm .AND..NOT.ace ) &
|
||||
CALL errore ('iosys', 'SCDM implemented only with ACE',1)
|
||||
scdm_den = scdmden
|
||||
scdm_grd = scdmgrd
|
||||
!
|
||||
IF(ecutfock <= 0.0_DP) THEN
|
||||
! default case
|
||||
|
|
|
@ -0,0 +1,261 @@
|
|||
! Copyright (C) 2005-2015 Quantum ESPRESSO group
|
||||
! This file is distributed under the terms of the
|
||||
! GNU General Public License. See the file `License'
|
||||
! in the root directory of the present distribution,
|
||||
! or http://www.gnu.org/copyleft/gpl.txt .
|
||||
!
|
||||
!--------------------------------------
|
||||
MODULE loc_scdm
|
||||
!--------------------------------------
|
||||
!
|
||||
! Variables and subroutines related to the calculation of
|
||||
! SCDM localization of molecular orbitals
|
||||
! Implements ACE: Lin Lin, J. Chem. Theory Comput. 2016, 12, 2242
|
||||
!
|
||||
USE kinds, ONLY : DP
|
||||
USE io_global, ONLY : stdout
|
||||
USE exx, ONLY : exx_fft, x_nbnd_occ, locbuff, nkqs
|
||||
|
||||
IMPLICIT NONE
|
||||
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
|
||||
|
||||
CONTAINS
|
||||
!
|
||||
!------------------------------------------------------------------------
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
SUBROUTINE SCDM_PGG(psi, NQR, nbnd_eff)
|
||||
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)
|
||||
!
|
||||
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
|
||||
|
||||
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
|
||||
|
||||
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
|
||||
|
||||
#if defined (__MPI)
|
||||
! ir_end = MIN(nnr,dfftp%nr1x*dfftp%nr2x*dfftp%npp(me_bgrp+1))
|
||||
ir_end = dfftp%nr1x*dfftp%my_nr2p*dfftp%my_nr3p
|
||||
#else
|
||||
ir_end = 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 / float(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 / float(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
|
||||
npt = npt + 1
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
nptot = npt
|
||||
cpu_npt(me_bgrp) = npt
|
||||
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(*,'(2(A,I8))') ' npt = ', npt, ' procID= ', me_bgrp
|
||||
write(stdout,*) ' reduced matrix, allocate: ', nptot
|
||||
write(stdout,*) ' cpu_npt = ', cpu_npt(:)
|
||||
|
||||
|
||||
! find the map of the index
|
||||
allocate( small(nbnd_eff,nptot), list(nptot) )
|
||||
small = 0.0d0
|
||||
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
|
||||
n = n + 1
|
||||
ncpu_start = sum(cpu_npt(0:me_bgrp-1))
|
||||
icpu = ncpu_start+n
|
||||
small(:,icpu) = psi(ir,:)
|
||||
list(icpu) = ir
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
! call matprt('small',nbnd_eff,nptot,small)
|
||||
call mp_sum(small,intra_bgrp_comm)
|
||||
! call matprt('small',nbnd_eff,nptot,small)
|
||||
call mp_sum(list,intra_bgrp_comm)
|
||||
|
||||
lwork = 4*nptot
|
||||
allocate( pivot(nptot), tau(nptot), work(lwork) )
|
||||
tau = 0.0d0
|
||||
work = 0.0d0
|
||||
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))
|
||||
! if(pivot(i).le.ncpu_end.and.pivot(i).ge.ncpu_start+1) &
|
||||
! write(*,'(A,3I9,2f16.8)') 'pivoting: ', me_bgrp, i, pivot(i), den(j), grad
|
||||
end do
|
||||
deallocate( den, grad_den)
|
||||
deallocate( tau, work )
|
||||
deallocate( 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
|
||||
! mat(:,i) = QRbuff(:,list(pivot(i)))
|
||||
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 matprt('Q',nbnd_eff,nbnd_eff,mat)
|
||||
call mp_sum(mat,intra_bgrp_comm)
|
||||
! call matprt('Q',nbnd_eff,nbnd_eff,mat)
|
||||
! deallocate( QRbuff )
|
||||
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)),:)
|
||||
! mat(i,:) = QRbuff(list(pivot(i)),:)
|
||||
end do
|
||||
deallocate( cpu_npt )
|
||||
! call matprt('Q2',nbnd_eff,nbnd_eff,mat)
|
||||
call mp_sum(mat,intra_bgrp_comm)
|
||||
! call matprt('Q2',nbnd_eff,nbnd_eff,mat)
|
||||
|
||||
! 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 invchol(n,A)
|
||||
USE exx, ONLY : errinfo
|
||||
IMPLICIT NONE
|
||||
!
|
||||
! given a matrix A, returns the inverse of the Cholesky decomposition of A
|
||||
! for real matrices
|
||||
!
|
||||
real(DP) :: A(n,n)
|
||||
integer :: n, INFO
|
||||
|
||||
INFO = -1
|
||||
CALL DPOTRF( 'L', n, A, n, INFO )
|
||||
CALL errinfo('DPOTRF','Cholesky failed in invchol.',INFO)
|
||||
INFO = -1
|
||||
CALL DTRTRI( 'L', 'N', n, A, n, INFO )
|
||||
CALL errinfo('DTRTRI','inversion failed in invchol.',INFO)
|
||||
|
||||
END SUBROUTINE
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
SUBROUTINE localize_orbitals( )
|
||||
!
|
||||
! Driver for SCDM orbital localization
|
||||
!
|
||||
USE noncollin_module, ONLY : npol
|
||||
USE exx, ONLY : locbuff, measure_localization, measure_localization_G
|
||||
!
|
||||
implicit none
|
||||
integer :: NQR
|
||||
|
||||
NQR = exx_fft%dfftt%nnr * npol
|
||||
|
||||
! CALL measure_localization(locbuff(1,1,1), NQR, x_nbnd_occ)
|
||||
CALL measure_localization_G(locbuff(1,1,1), NQR, x_nbnd_occ)
|
||||
CALL SCDM_PGG(locbuff(1,1,1), NQR, x_nbnd_occ)
|
||||
! CALL measure_localization(locbuff(1,1,1), NQR, x_nbnd_occ)
|
||||
CALL measure_localization_G(locbuff(1,1,1), NQR, x_nbnd_occ)
|
||||
|
||||
END SUBROUTINE localize_orbitals
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
END MODULE loc_scdm
|
||||
!-----------------------------------------------------------------------
|
|
@ -412,6 +412,7 @@ electrons.o : extfield.o
|
|||
electrons.o : exx.o
|
||||
electrons.o : io_rho_xml.o
|
||||
electrons.o : ldaU.o
|
||||
electrons.o : loc_scdm.o
|
||||
electrons.o : newd.o
|
||||
electrons.o : paw_onecenter.o
|
||||
electrons.o : paw_symmetry.o
|
||||
|
@ -885,6 +886,7 @@ init_vloc.o : pwcom.o
|
|||
input.o : ../../KS_Solvers/CG/constants.o
|
||||
input.o : ../../Modules/bfgs_module.o
|
||||
input.o : ../../Modules/cell_base.o
|
||||
input.o : ../../Modules/check_stop.o
|
||||
input.o : ../../Modules/constraints_module.o
|
||||
input.o : ../../Modules/control_flags.o
|
||||
input.o : ../../Modules/fcp_variables.o
|
||||
|
@ -921,6 +923,7 @@ input.o : esm.o
|
|||
input.o : extfield.o
|
||||
input.o : exx.o
|
||||
input.o : ldaU.o
|
||||
input.o : loc_scdm.o
|
||||
input.o : martyna_tuckerman.o
|
||||
input.o : pwcom.o
|
||||
input.o : realus.o
|
||||
|
@ -967,6 +970,19 @@ ldaU.o : ../../Modules/ions_base.o
|
|||
ldaU.o : ../../Modules/kind.o
|
||||
ldaU.o : ../../Modules/parameters.o
|
||||
ldaU.o : atomic_wfc_mod.o
|
||||
loc_scdm.o : ../../KS_Solvers/CG/constants.o
|
||||
loc_scdm.o : ../../Modules/cell_base.o
|
||||
loc_scdm.o : ../../Modules/fft_base.o
|
||||
loc_scdm.o : ../../Modules/io_global.o
|
||||
loc_scdm.o : ../../Modules/ions_base.o
|
||||
loc_scdm.o : ../../Modules/kind.o
|
||||
loc_scdm.o : ../../Modules/mp_bands.o
|
||||
loc_scdm.o : ../../Modules/noncol.o
|
||||
loc_scdm.o : ../../Modules/xc_vdW_DF.o
|
||||
loc_scdm.o : ../../UtilXlib/mp.o
|
||||
loc_scdm.o : exx.o
|
||||
loc_scdm.o : pwcom.o
|
||||
loc_scdm.o : scf_mod.o
|
||||
make_pointlists.o : ../../Modules/cell_base.o
|
||||
make_pointlists.o : ../../Modules/fft_base.o
|
||||
make_pointlists.o : ../../Modules/io_global.o
|
||||
|
|
|
@ -165,6 +165,8 @@ SUBROUTINE print_clock_pw()
|
|||
CALL print_clock( 'fft_scatt_tg' )
|
||||
CALL print_clock( 'ALLTOALL' )
|
||||
#endif
|
||||
CALL print_clock( 'localization' )
|
||||
CALL print_clock( 'measure' )
|
||||
!
|
||||
IF ( lda_plus_U ) THEN
|
||||
WRITE( stdout, '(/,5X,"Hubbard U routines")' )
|
||||
|
@ -179,10 +181,10 @@ SUBROUTINE print_clock_pw()
|
|||
CALL print_clock( 'exx_grid' )
|
||||
CALL print_clock( 'exxinit' )
|
||||
CALL print_clock( 'vexx' )
|
||||
!civn
|
||||
CALL print_clock( 'matcalc' )
|
||||
CALL print_clock( 'aceupdate' )
|
||||
CALL print_clock( 'vexxace' )
|
||||
CALL print_clock( 'vexxloc' )
|
||||
CALL print_clock( 'aceinit' )
|
||||
CALL print_clock( 'exxenergy' )
|
||||
IF( okvan) THEN
|
||||
|
|
Loading…
Reference in New Issue