diff --git a/GWW/pw4gww/produce_wannier_gamma.f90 b/GWW/pw4gww/produce_wannier_gamma.f90 index d14392996..ba7b86b17 100644 --- a/GWW/pw4gww/produce_wannier_gamma.f90 +++ b/GWW/pw4gww/produce_wannier_gamma.f90 @@ -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 diff --git a/Modules/input_parameters.f90 b/Modules/input_parameters.f90 index 23bff3883..a32b676fa 100644 --- a/Modules/input_parameters.f90 +++ b/Modules/input_parameters.f90 @@ -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, & diff --git a/Modules/read_namelists.f90 b/Modules/read_namelists.f90 index b2dc8f4db..fa20e3123 100644 --- a/Modules/read_namelists.f90 +++ b/Modules/read_namelists.f90 @@ -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 ) diff --git a/PW/src/Makefile b/PW/src/Makefile index ce82765c5..2dcc53434 100644 --- a/PW/src/Makefile +++ b/PW/src/Makefile @@ -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 \ diff --git a/PW/src/electrons.f90 b/PW/src/electrons.f90 index aed6ad11b..07bac9056 100644 --- a/PW/src/electrons.f90 +++ b/PW/src/electrons.f90 @@ -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 diff --git a/PW/src/exx.f90 b/PW/src/exx.f90 index b4910ac16..4c510256b 100644 --- a/PW/src/exx.f90 +++ b/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 = - 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 = + 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> ^(-1) +! 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 !----------------------------------------------------------------------- diff --git a/PW/src/input.f90 b/PW/src/input.f90 index 5b4aa55ab..3c4312533 100644 --- a/PW/src/input.f90 +++ b/PW/src/input.f90 @@ -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 diff --git a/PW/src/loc_scdm.f90 b/PW/src/loc_scdm.f90 new file mode 100644 index 000000000..9bdce743e --- /dev/null +++ b/PW/src/loc_scdm.f90 @@ -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 +!----------------------------------------------------------------------- diff --git a/PW/src/make.depend b/PW/src/make.depend index 869a912f0..477b14267 100644 --- a/PW/src/make.depend +++ b/PW/src/make.depend @@ -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 diff --git a/PW/src/print_clock_pw.f90 b/PW/src/print_clock_pw.f90 index 150fd9f92..8d8a52d25 100644 --- a/PW/src/print_clock_pw.f90 +++ b/PW/src/print_clock_pw.f90 @@ -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