From e8706afc54112047e7f6cbc19f7397be9cb57f03 Mon Sep 17 00:00:00 2001 From: paulatto Date: Fri, 9 Nov 2007 12:45:43 +0000 Subject: [PATCH] Added routine to symmetrize becsum when PAW calculation is being performed. It still needs some test. As a collateral effect d1,d2 and d3 variables in modules symme (pwcom.f90) are now marked as TARGET. File make.depend updated accordingly. git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@4412 c92efa57-630b-4861-b058-cf58834340f0 --- PW/grid_paw_routines.f90 | 6 +- PW/make.depend | 4 + PW/print_clock_pw.f90 | 1 + PW/pwcom.f90 | 3 +- PW/rad_paw_routines.f90 | 175 +++++++++++++++++++++++++++++++++++++-- PW/sum_band.f90 | 9 ++ 6 files changed, 189 insertions(+), 9 deletions(-) diff --git a/PW/grid_paw_routines.f90 b/PW/grid_paw_routines.f90 index 576a6527a..42437caf6 100644 --- a/PW/grid_paw_routines.f90 +++ b/PW/grid_paw_routines.f90 @@ -39,6 +39,7 @@ SUBROUTINE atomic_becsum() USE ions_base, ONLY : nat, ityp USE lsda_mod, ONLY : nspin, starting_magnetization USE grid_paw_variables, ONLY : tpawp, okpaw + USE rad_paw_routines, ONLY : PAW_symmetrize IMPLICIT NONE INTEGER :: ispin, na, nt, ijh, ih, jh, nb, mb ! @@ -68,7 +69,8 @@ SUBROUTINE atomic_becsum() END IF ijh = ijh + 1 ! - jh_loop: DO jh = ( ih + 1 ), nh(nt) + jh_loop: & + DO jh = ( ih + 1 ), nh(nt) !mb = indv(jh,nt) DO ispin = 1, nspin becsum(ijh,na,ispin) = 0._DP @@ -84,6 +86,8 @@ SUBROUTINE atomic_becsum() PRINT '(1f20.10)', becsum(:,1,1) #endif + CALL PAW_symmetrize(becsum) + END SUBROUTINE atomic_becsum ! Analogous to PW/allocate_nlpot.f90 diff --git a/PW/make.depend b/PW/make.depend index 72577c6bd..81ffde1c0 100644 --- a/PW/make.depend +++ b/PW/make.depend @@ -520,6 +520,7 @@ grid_paw_routines.o : ../Modules/uspp.o grid_paw_routines.o : ../Modules/wavefunctions.o grid_paw_routines.o : noncol.o grid_paw_routines.o : pwcom.o +grid_paw_routines.o : rad_paw_routines.o grid_paw_routines.o : scf_mod.o gweights.o : ../Modules/kind.o h_1psi.o : ../Modules/kind.o @@ -878,6 +879,7 @@ qvan2.o : ../Modules/uspp.o qvan2.o : pwcom.o rad_paw_routines.o : ../Modules/atom.o rad_paw_routines.o : ../Modules/constants.o +rad_paw_routines.o : ../Modules/control_flags.o rad_paw_routines.o : ../Modules/functionals.o rad_paw_routines.o : ../Modules/grid_paw_variables.o rad_paw_routines.o : ../Modules/ions_base.o @@ -1213,6 +1215,7 @@ struct_fact.o : ../Modules/kind.o sum_band.o : ../Modules/cell_base.o sum_band.o : ../Modules/control_flags.o sum_band.o : ../Modules/functionals.o +sum_band.o : ../Modules/grid_paw_variables.o sum_band.o : ../Modules/io_files.o sum_band.o : ../Modules/ions_base.o sum_band.o : ../Modules/kind.o @@ -1223,6 +1226,7 @@ sum_band.o : ../Modules/wavefunctions.o sum_band.o : buffers.o sum_band.o : noncol.o sum_band.o : pwcom.o +sum_band.o : rad_paw_routines.o sum_band.o : realus.o sum_band.o : scf_mod.o sumkg.o : ../Modules/kind.o diff --git a/PW/print_clock_pw.f90 b/PW/print_clock_pw.f90 index d9413c217..b1f54a0c4 100644 --- a/PW/print_clock_pw.f90 +++ b/PW/print_clock_pw.f90 @@ -200,6 +200,7 @@ SUBROUTINE print_clock_pw() CALL print_clock ('PAW_ddot') CALL print_clock ('PAW_rad_init') CALL print_clock ('PAW_energy') + CALL print_clock ('PAW_symme') ! second level routines: CALL print_clock ('PAW_rho_lm') CALL print_clock ('PAW_h_pot') diff --git a/PW/pwcom.f90 b/PW/pwcom.f90 index 60d603d6e..8bfee438c 100644 --- a/PW/pwcom.f90 +++ b/PW/pwcom.f90 @@ -209,11 +209,10 @@ MODULE symme irt(:,:) ! symmetric atom for each atom and sym.op. LOGICAL :: & invsym ! if .TRUE. the system has inversion symmetry - REAL(DP) :: & + REAL(DP),TARGET :: & d1(3,3,48), &! matrices for rotating spherical d2(5,5,48), &! harmonics (d1 for l=1, ...) d3(7,7,48) ! - ! END MODULE symme diff --git a/PW/rad_paw_routines.f90 b/PW/rad_paw_routines.f90 index b182d8751..269ce9959 100644 --- a/PW/rad_paw_routines.f90 +++ b/PW/rad_paw_routines.f90 @@ -54,12 +54,13 @@ MODULE rad_paw_routines IMPLICIT NONE ! entry points: - PUBLIC :: PAW_potential ! prepare paw potential and store it, - ! also computes energy if required - PUBLIC :: PAW_integrate ! computes \int v(r) \times n(r) dr - PUBLIC :: PAW_ddot ! error estimate for mix_rho - PUBLIC :: PAW_init ! initialize - PUBLIC :: PAW_newd ! computes descreening coefficients + PUBLIC :: PAW_potential ! prepare paw potential and store it, + ! also computes energy if required + PUBLIC :: PAW_integrate ! computes \int v(r) \times n(r) dr + PUBLIC :: PAW_ddot ! error estimate for mix_rho + PUBLIC :: PAW_init ! initialize + PUBLIC :: PAW_newd ! computes descreening coefficients + PUBLIC :: PAW_symmetrize ! symmetrize becsums ! PRIVATE SAVE @@ -233,6 +234,167 @@ SUBROUTINE PAW_potential(becsum, energy, e_cmp) END SUBROUTINE PAW_potential +SUBROUTINE PAW_symmetrize(becsum) + USE kinds, ONLY : DP + USE lsda_mod, ONLY : nspin + USE uspp_param, ONLY : nhm, lmaxq + USE ions_base, ONLY : nat, ityp + USE symme, ONLY : nsym, irt, d1, d2, d3 + USE uspp, ONLY : nhtolm,nhtol + USE uspp_param, ONLY : nh + USE grid_paw_variables, ONLY : tpawp + USE wvfct, ONLY : gamma_only + USE control_flags, ONLY : nosym + + REAL(DP), INTENT(INOUT) :: becsum(nhm*(nhm+1)/2,nat,nspin)! cross band occupations + + REAL(DP) :: becsym(nhm*(nhm+1)/2,nat,nspin)! symmetrized becsum +! REAL(DP) :: bectwo(nhm*(nhm+1)/2,nat,nspin)! twice symmetrized becsum + + INTEGER :: is, na, nt ! counters on spin, atom, atom-type + INTEGER :: ma ! atom symmetric to na + INTEGER :: ih,jh, ijh ! counters for augmentation channels + INTEGER :: lm_i, lm_j, &! angular momentums of non-symmetrized becsum + l_i, l_j, m_i, m_j + INTEGER :: p_i, p_j ! projector indexes of ih and jh (tricky!) + INTEGER :: m_o, m_u ! counters for sums on m + INTEGER :: oh, uh, ouh ! auxiliary indexes corresponding to m_o and m_u + INTEGER :: isym ! counter for symmetry operation + CHARACTER(len=8) :: do_symme + + ! The following mess is necessary because the symmetrization operation + ! in LDA+U code is simpler than in PAW, so the required quantities are + ! represented in a simple but not general way. + ! I will fix this when everything works. + REAL(DP), TARGET :: d0(1,1,48) + TYPE symmetryzation_tensor + REAL(DP),POINTER :: d(:,:,:) + END TYPE symmetryzation_tensor + TYPE(symmetryzation_tensor) :: D(0:3) + d0(1,1,:) = 1._dp + D(0)%d => d0 ! d0(1,1,48) + D(1)%d => d1 ! d1(3,3,48) + D(2)%d => d2 ! d2(5,5,48) + D(3)%d => d3 ! d3(7,7,48) + +! => lm = l**2 + m +! => ih = lm + (l+proj)**2 <-- if proj starts from zero! +! = lm + proj**2 + 2*l*proj +! = m + l**2 + proj**2 + 2*l*proj +! ^^^ +! Known ih and m_i I can compute the index oh of a different m = m_o but +! the same augmentation channel (l_i = l_o, proj_i = proj_o): +! +! oh = ih - m_i + m_o +#define __DEBUG_PAW_SYM +#ifdef __DEBUG_PAW_SYM + nt = 1 + na = 1 + is = 1 + DO ih = 1, nh(nt) + DO jh = 1, nh(nt) + IF (jh > ih) THEN + ijh = nh(nt)*(ih-1) - ih*(ih-1)/2 + jh + ELSE + ijh = nh(nt)*(jh-1) - jh*(jh-1)/2 + ih + ENDIF + write(*,"(1f10.5)", advance='no') becsum(ijh,na,is) + ENDDO + write(*,*) + ENDDO +#endif + + IF( gamma_only .or. nosym ) RETURN + + CALL start_clock('PAW_symme') + + becsym(:,:,:) = 0._dp + + DO is = 1, nspin + DO na = 1, nat + nt = ityp(na) + ! No need to symmetrize non-PAW atoms + IF ( .not. tpawp(nt) ) CYCLE + ! + DO ih = 1, nh(nt) + DO jh = ih, nh(nt) ! note: jh >= ih + ijh = nh(nt)*(ih-1) - ih*(ih-1)/2 + jh + ! + lm_i = nhtolm(ih,nt) + lm_j = nhtolm(jh,nt) + ! + l_i = nhtol(ih,nt) + l_j = nhtol(jh,nt) + ! + m_i = lm_i - l_i**2 + m_j = lm_j - l_j**2 + ! + DO isym = 1,nsym + ma = irt(isym,na) + DO m_o = 1, 2*l_i +1 + DO m_u = 1, 2*l_j +1 + oh = ih - m_i + m_o + uh = jh - m_j + m_u + IF ( oh >= uh ) THEN + ouh = nh(nt)*(uh-1) - uh*(uh-1)/2 + oh + ELSE + ouh = nh(nt)*(oh-1) - oh*(oh-1)/2 + uh + ENDIF + ! + becsym(ijh, na, is) = becsym(ijh, na, is) & + + D(l_i)%d(m_o,m_i, isym) * D(l_j)%d(m_u,m_j, isym) & + * becsum(ouh, ma, is) / nsym + ENDDO ! m_o + ENDDO ! m_u + ENDDO ! isym + ENDDO ! ih + ENDDO ! jh + ENDDO ! nat + ENDDO ! nspin + +#define __DEBUG_PAW_SYM +#ifdef __DEBUG_PAW_SYM + nt = 1 + na = 1 + is = 1 + write(*,*) + DO ih = 1, nh(nt) + DO jh = 1, nh(nt) + IF (jh > ih) THEN + ijh = nh(nt)*(ih-1) - ih*(ih-1)/2 + jh + ELSE + ijh = nh(nt)*(jh-1) - jh*(jh-1)/2 + ih + ENDIF + write(*,"(1f10.5)", advance='no') becsym(ijh,na,is) + ENDDO + write(*,*) + ENDDO + write(*,*) + DO ih = 1, nh(nt) + DO jh = 1, nh(nt) + IF (jh > ih) THEN + ijh = nh(nt)*(ih-1) - ih*(ih-1)/2 + jh + ELSE + ijh = nh(nt)*(jh-1) - jh*(jh-1)/2 + ih + ENDIF + write(*,"(1f10.5)", advance='no') becsym(ijh,na,is)-becsum(ijh,na,is) + ENDDO + write(*,*) + ENDDO + write(*,*) "________________________________________________" +#endif + + ! Apply symmetrization: + CALL get_environment_variable('DO_SYMME', do_symme) + IF (trim(do_symme)=='.false.') THEN + write(*,*) "**** skipping PAW symmetrization!" + ELSE + becsum(:,:,:) = becsym(:,:,:) + ENDIF + + CALL stop_clock('PAW_symme') + +END SUBROUTINE PAW_symmetrize !___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! @@ -1349,6 +1511,7 @@ SUBROUTINE PAW_rho_lm(i, becsum, pfunc, rho_lm, aug) ! This subroutine computes the angular momentum components of rho ! using the following formula: ! rho(\vec{r}) = \sum_{LM} Y_{LM} \sum_{i,j} (\hat{r}) a_{LM}^{(lm)_i(lm)_j} becsum_ij pfunc_ij(r) + ! where a_{LM}^{(lm)_i(lm)_j} are the Clebsh-Gordan coefficients. ! ! actually different angular momentum components are stored separately: ! rho^{LM}(\vec{r}) = \sum_{i,j} (\hat{r}) a_{LM}^{(lm)_i(lm)_j} becsum_ij pfunc_ij(r) diff --git a/PW/sum_band.f90 b/PW/sum_band.f90 index 7a6ef76b9..763c9a0ba 100644 --- a/PW/sum_band.f90 +++ b/PW/sum_band.f90 @@ -43,6 +43,8 @@ SUBROUTINE sum_band() root_image, npool, my_pool_id USE mp, ONLY : mp_bcast USE funct, ONLY : dft_is_meta + USE rad_paw_routines, ONLY : PAW_symmetrize + USE grid_paw_variables, ONLY : okpaw ! IMPLICIT NONE ! @@ -110,6 +112,7 @@ SUBROUTINE sum_band() ! END IF ! + ! ! ... If a double grid is used, interpolate onto the fine grid ! IF ( doublegrid ) THEN @@ -186,6 +189,12 @@ SUBROUTINE sum_band() END IF ! #endif + ! ... Needed for PAW: becsum has to be symmetrized so that they reflect a real integral + ! in k-space, not only on the irreducible zone. For USPP there is no need to do this as + ! becsums are only used to compute the density, which is symmetrized later. + ! + IF ( okpaw ) CALL PAW_symmetrize(becsum) + ! if (dft_is_meta() ) deallocate (kplusg) ! ! ... synchronize rho%of_g to the calculated rho%of_r