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

Miscellaneous improvements in TDDFPT: 

1) Use the new routine h_prec (preconditioning matrix)
2) Use the routine g2_kin
3) Adding more documentation
4) Extension of compute_d0psi_rs to k-points version


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@12514 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
timrov 2016-06-21 12:35:43 +00:00
parent 8ab9d17e00
commit 288b3e3e9c
10 changed files with 297 additions and 305 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
TDDFPT/Examples/make.sys
View File

@ -7,7 +7,7 @@
#RUNNER = openmpi-1.4.3/bin/mpirun -np 8
#PWSCF = "../../../bin/pw.x" -ndiag 4
RUNNER =
RUNNER =
PWSCF = "../../../bin/pw.x"
TDDFPT_LANCZOS = "../../../bin/turbo_lanczos.x"
TDDFPT_DAVIDSON = "../../../bin/turbo_davidson.x"

82
TDDFPT/src/lr_apply_liouvillian.f90
View File

@ -21,13 +21,14 @@ SUBROUTINE lr_apply_liouvillian( evc1, evc1_new, interaction )
!
! interaction=.true. corresponds to eq.(32)
! interaction=.false. corresponds to eq.(33)
! in Ralph Gebauer, Brent Walker J. Chem. Phys., 127, 164106 (2007)
! in B. Walker and R. Gebauer, J. Chem. Phys., 127, 164106 (2007)
!
! Modified by Osman Baris Malcioglu in 2009
! Modified by Simone Binnie in 2012 (EXX)
! Modified by Xiaochuan Ge in 2013 (Davidson)
! Modified by Iurii Timrov in 2014 (Environ)
!
USE kinds, ONLY : DP
USE ions_base, ONLY : ityp, nat, ntyp=>nsp
USE cell_base, ONLY : tpiba2
USE fft_base, ONLY : dffts, dfftp
@ -35,7 +36,6 @@ SUBROUTINE lr_apply_liouvillian( evc1, evc1_new, interaction )
USE gvecs, ONLY : nls, nlsm
USE gvect, ONLY : nl, ngm, gstart, g, gg
USE io_global, ONLY : stdout
USE kinds, ONLY : dp
USE klist, ONLY : nks, xk, ngk, igk_k
USE lr_variables, ONLY : evc0, sevc0, revc0, rho_1, rho_1c, &
& ltammd, size_evc, no_hxc, lr_exx, &
@ -69,23 +69,22 @@ SUBROUTINE lr_apply_liouvillian( evc1, evc1_new, interaction )
!
IMPLICIT NONE
!
COMPLEX(kind=dp),INTENT(in) :: evc1(npwx*npol,nbnd,nks)
COMPLEX(kind=dp),INTENT(out) :: evc1_new(npwx*npol,nbnd,nks)
LOGICAL, INTENT(in) :: interaction
COMPLEX(DP), INTENT(IN) :: evc1(npwx*npol,nbnd,nks)
COMPLEX(DP), INTENT(OUT) :: evc1_new(npwx*npol,nbnd,nks)
LOGICAL, INTENT(IN) :: interaction
!
! Local variables
!
INTEGER :: ir, ibnd, ik, ig, ia, mbia
INTEGER :: ijkb0, na, nt, ih, jh, ikb, jkb, iqs,jqs
REAL(kind=dp), ALLOCATABLE :: dvrs(:,:), dvrss(:), &
& d_deeq(:,:,:,:), w1(:), w2(:)
COMPLEX(kind=dp), ALLOCATABLE :: dvrs_temp(:,:), spsi1(:,:), &
& dvrsc(:,:), dvrssc(:), &
& sevc1_new(:,:,:)
REAL(DP), ALLOCATABLE :: dvrs(:,:), dvrss(:), d_deeq(:,:,:,:), &
& w1(:), w2(:)
COMPLEX(DP), ALLOCATABLE :: dvrs_temp(:,:), spsi1(:,:), dvrsc(:,:), &
& dvrssc(:), sevc1_new(:,:,:)
!
! Environ related arrays
!
REAL(kind=dp), ALLOCATABLE :: &
REAL(DP), ALLOCATABLE :: &
dv_pol(:), & ! response polarization potential
dv_epsilon(:) ! response dielectric potential
!
@ -157,7 +156,7 @@ SUBROUTINE lr_apply_liouvillian( evc1, evc1_new, interaction )
!
ALLOCATE( dvrs_temp(dfftp%nnr, nspin) )
!
dvrs_temp = CMPLX( dvrs, 0.0d0, kind=dp )
dvrs_temp = CMPLX( dvrs, 0.0d0, kind=DP )
!
DEALLOCATE ( dvrs ) ! to save memory
!
@ -254,7 +253,7 @@ SUBROUTINE lr_apply_liouvillian( evc1, evc1_new, interaction )
! Here we add the two terms:
! [H(k) - E(k)] * evc1(k) + dV_HXC * revc0(k)
!
CALL zaxpy(size_evc,CMPLX(1.0d0,0.0d0,kind=dp),&
CALL zaxpy(size_evc,CMPLX(1.0d0,0.0d0,kind=DP),&
& evc1_new(:,:,:), 1, sevc1_new(:,:,:), 1)
!
ELSEIF ( interaction .and. ltammd ) THEN
@ -266,7 +265,7 @@ SUBROUTINE lr_apply_liouvillian( evc1, evc1_new, interaction )
!
! Here evc1_new contains the interaction
!
CALL zaxpy(size_evc,CMPLX(0.5d0,0.0d0,kind=dp),&
CALL zaxpy(size_evc,CMPLX(0.5d0,0.0d0,kind=DP),&
& evc1_new(:,:,:) , 1, sevc1_new(:,:,:),1)
!
ELSE
@ -279,7 +278,7 @@ SUBROUTINE lr_apply_liouvillian( evc1, evc1_new, interaction )
ENDIF
!
IF (gstart == 2 .AND. gamma_only ) sevc1_new(1,:,:) = &
& CMPLX( REAL( sevc1_new(1,:,:), dp ), 0.0d0 ,dp )
& CMPLX( REAL( sevc1_new(1,:,:), DP ), 0.0d0, DP )
!
IF (gstart==2 .and. gamma_only) THEN
DO ik=1,nks
@ -358,8 +357,8 @@ CONTAINS
!
IMPLICIT NONE
!
REAL(kind=dp), ALLOCATABLE :: becp2(:,:)
REAL(kind=dp), ALLOCATABLE :: tg_dvrss(:)
REAL(DP), ALLOCATABLE :: becp2(:,:)
REAL(DP), ALLOCATABLE :: tg_dvrss(:)
LOGICAL :: use_tg
INTEGER :: v_siz, incr, ioff
INTEGER :: ibnd_start_gamma, ibnd_end_gamma
@ -454,7 +453,7 @@ CONTAINS
!
DO ir=1, dffts%nr1x*dffts%nr2x*dffts%tg_npp( me_bgrp + 1 )
!
tg_psic(ir) = tg_revc0(ir,ibnd,1)*CMPLX(tg_dvrss(ir),0.0d0,dp)
tg_psic(ir) = tg_revc0(ir,ibnd,1)*CMPLX(tg_dvrss(ir),0.0d0,DP)
!
ENDDO
!
@ -462,7 +461,7 @@ CONTAINS
!
DO ir = 1,dffts%nnr
!
psic(ir) = revc0(ir,ibnd,1)*CMPLX(dvrss(ir),0.0d0,dp)
psic(ir) = revc0(ir,ibnd,1)*CMPLX(dvrss(ir),0.0d0,DP)
!
ENDDO
!
@ -558,14 +557,14 @@ CONTAINS
!
IF (lr_exx .AND. .NOT.interaction) CALL lr_exx_kernel_noint(evc1,evc1_new)
!
! Note: In the gamma_only case the kinetic energy g2kin has
! been already computed in lr_dvpsi_e.
! The kinetic energy g2kin was already computed when
! calling the routine lr_solve_e.
!
! Call h_psi on evc1 such that h.evc1 = sevc1_new
! Compute sevc1_new = H*evc1
!
CALL h_psi(npwx,ngk(1),nbnd,evc1(1,1,1),sevc1_new(1,1,1)) ! g2kin is needed here
CALL h_psi(npwx,ngk(1),nbnd,evc1(1,1,1),sevc1_new(1,1,1))
!
! spsi1 = s*evc1
! Compute spsi1 = S*evc1
!
IF (real_space_debug > 9 ) THEN
DO ibnd = 1,nbnd,2
@ -581,7 +580,7 @@ CONTAINS
!
DO ibnd = 1,nbnd
!
CALL zaxpy(ngk(1), CMPLX(-(et(ibnd,1)-scissor),0.0d0,dp), &
CALL zaxpy(ngk(1), CMPLX(-(et(ibnd,1)-scissor),0.0d0,DP), &
& spsi1(:,ibnd), 1, sevc1_new(:,ibnd,1), 1)
!
ENDDO
@ -601,7 +600,7 @@ SUBROUTINE lr_apply_liouvillian_k()
!
IMPLICIT NONE
!
COMPLEX(kind=dp), ALLOCATABLE :: becp2(:,:)
COMPLEX(DP), ALLOCATABLE :: becp2(:,:)
!
IF ( nkb > 0 .AND. okvan ) THEN
!
@ -711,35 +710,32 @@ SUBROUTINE lr_apply_liouvillian_k()
!
IF (lr_exx .AND. .NOT.interaction) CALL lr_exx_kernel_noint(evc1,evc1_new)
!
DO ik=1,nks
DO ik = 1, nks
!
CALL init_us_2(ngk(ik),igk_k(1,ik),xk(1,ik),vkb)
! US case: Compute the beta function vkb at point k
!
DO ig = 1, ngk(ik)
!
g2kin(ig)=((xk(1,ik)+g(1,igk_k(ig,ik)))**2 &
+(xk(2,ik)+g(2,igk_k(ig,ik)))**2 &
+(xk(3,ik)+g(3,igk_k(ig,ik)))**2)*tpiba2
!
ENDDO
CALL init_us_2(ngk(ik), igk_k(1,ik), xk(1,ik), vkb)
!
! Compute the kinetic energy g2kin: (k+G)^2
!
CALL g2_kin(ik)
!
current_k = ik
!
! Compute sevc1_new = H*evc1
!
CALL h_psi(npwx,ngk(ik),nbnd,evc1(1,1,ik),sevc1_new(1,1,ik))
!
! Compute spsi1 = S*evc1
!
CALL s_psi(npwx,ngk(ik),nbnd,evc1(1,1,ik),spsi1)
!
! Subtract the eigenvalues
! IT: One may want to add scissor
!
DO ibnd=1,nbnd
DO ibnd = 1, nbnd
!
DO ig=1,ngk(ik)
!
sevc1_new(ig,ibnd,ik)=sevc1_new(ig,ibnd,ik) &
-cmplx(et(ibnd,ik),0.0d0,dp)*spsi1(ig,ibnd)
!
ENDDO
CALL zaxpy(ngk(ik), CMPLX(-(et(ibnd,ik)-scissor),0.0d0,dp), &
& spsi1(:,ibnd), 1, sevc1_new(:,ibnd,ik), 1)
!
ENDDO
!

33
TDDFPT/src/lr_apply_liouvillian_eels.f90
View File

@ -16,21 +16,13 @@ SUBROUTINE lr_apply_liouvillian_eels ( evc1, evc1_new, interaction )
!
! Written by Iurii Timrov (2013)
!
USE ions_base, ONLY : ityp, nat, ntyp=>nsp
USE cell_base, ONLY : tpiba2
USE kinds, ONLY : DP
USE fft_base, ONLY : dfftp, dffts
USE fft_parallel, ONLY : tg_cgather
USE fft_interfaces, ONLY : fwfft, invfft
USE gvecs, ONLY : nls, nlsm, ngms, doublegrid
USE gvect, ONLY : nl, nlm, ngm, g, gg
USE io_global, ONLY : stdout
USE kinds, ONLY : dp
USE klist, ONLY : nks, xk, igk_k, ngk
USE lr_variables, ONLY : evc0, no_hxc
USE lsda_mod, ONLY : nspin, current_spin
USE wvfct, ONLY : nbnd, npwx, g2kin, et, current_k
USE gvecw, ONLY : gcutw
USE io_global, ONLY : stdout
USE klist, ONLY : xk, igk_k, ngk
USE lr_variables, ONLY : no_hxc
USE lsda_mod, ONLY : current_spin
USE wvfct, ONLY : nbnd, npwx, et, current_k
USE uspp, ONLY : vkb
USE io_files, ONLY : iunwfc, nwordwfc
USE wavefunctions_module, ONLY : evc, psic, psic_nc
@ -46,12 +38,12 @@ SUBROUTINE lr_apply_liouvillian_eels ( evc1, evc1_new, interaction )
IMPLICIT NONE
!
COMPLEX(kind=dp),INTENT(in) :: evc1(npwx*npol,nbnd,nksq)
COMPLEX(kind=dp),INTENT(out) :: evc1_new(npwx*npol,nbnd,nksq)
COMPLEX(DP), INTENT(IN) :: evc1(npwx*npol,nbnd,nksq)
COMPLEX(DP), INTENT(OUT) :: evc1_new(npwx*npol,nbnd,nksq)
LOGICAL, INTENT(IN) :: interaction
!
! Local variables
!
LOGICAL, INTENT(in) :: interaction
LOGICAL :: interaction1
INTEGER :: i,j,ir, ibnd, ig, ia, ios, is, incr, v_siz, ipol
INTEGER :: ik, &
@ -63,8 +55,7 @@ SUBROUTINE lr_apply_liouvillian_eels ( evc1, evc1_new, interaction )
& sevc1_new(:,:,:), &
! Change of the Hartree and exchange-correlation (HXC) potential
& tg_psic(:,:), tg_dvrssc(:,:)
! Task groups: wfct in R-space
! Task groups: HXC potential
! Task groups: wfct in R-space and the response HXC potential
!
CALL start_clock('lr_apply')
!
@ -266,11 +257,7 @@ SUBROUTINE lr_apply_liouvillian_eels ( evc1, evc1_new, interaction )
!
! Compute the kinetic energy g2kin: (k+q+G)^2
!
DO ig = 1, npwq
g2kin (ig) = ( (xk (1,ikq) + g (1, igk_k(ig,ikq)) ) **2 + &
(xk (2,ikq) + g (2, igk_k(ig,ikq)) ) **2 + &
(xk (3,ikq) + g (3, igk_k(ig,ikq)) ) **2 ) * tpiba2
ENDDO
CALL g2_kin(ikq)
!
IF (noncolin) THEN
IF (.NOT. ALLOCATED(psic_nc)) ALLOCATE(psic_nc(dfftp%nnr,npol))

1
TDDFPT/src/lr_calc_dens.f90
View File

@ -323,7 +323,6 @@ CONTAINS
USE io_global, ONLY : stdout
USE realus, ONLY : real_space, invfft_orbital_gamma,&
& initialisation_level,&
& fwfft_orbital_gamma,&
& calbec_rs_gamma,&
& add_vuspsir_gamma, v_loc_psir,&
& real_space_debug

4
TDDFPT/src/lr_dealloc.f90
View File

@ -64,11 +64,11 @@ SUBROUTINE lr_dealloc()
IF (allocated(dmuxc)) DEALLOCATE(dmuxc)
IF (allocated(igk_k)) DEALLOCATE(igk_k)
IF (allocated(ngk)) DEALLOCATE(ngk)
IF (allocated(ikks)) DEALLOCATE(ikks)
IF (allocated(ikqs)) DEALLOCATE(ikqs)
!
! EELS-related variables
!
IF (allocated(ikks)) DEALLOCATE(ikks)
IF (allocated(ikqs)) DEALLOCATE(ikqs)
IF (allocated(dpsi)) DEALLOCATE(dpsi)
IF (allocated(dvpsi)) DEALLOCATE(dvpsi)
IF (allocated(eigqts)) DEALLOCATE(eigqts)

127
TDDFPT/src/lr_dvpsi_e.f90
View File

@ -12,14 +12,15 @@ SUBROUTINE lr_dvpsi_e(ik,ipol,dvpsi)
! On output: dvpsi contains P_c^+ x | psi_ik > in crystal axis
! (projected on at(*,ipol) )
!
! dvpsi is COMPUTED and WRITTEN on file (vkb and evc must be set)
! OBM: ^ This is now handled elesewhere
! dvpsi is computed here (vkb and evc must be set)
! and it is written on file in the routine lr_solve_e.
!
! See J. Tobik and A. Dal Corso, JCP 120, 9934 (2004)
! See Ref.[1] : J. Tobik and A. Dal Corso, JCP 120, 9934 (2004)
! for the details of the theory implemented in this routine.
!
! Modified by Osman Baris Malcioglu (2009)
! Rebased wrt PHONON routines. S J Binnie (2011)
! Modified by Iurii Timrov (2016)
!
USE kinds, ONLY : DP
USE cell_base, ONLY : tpiba2, at
@ -52,7 +53,7 @@ SUBROUTINE lr_dvpsi_e(ik,ipol,dvpsi)
! npw: number of plane-waves at point ik
REAL(kind=dp) :: atnorm
COMPLEX(kind=dp),ALLOCATABLE :: d0psi(:,:)
REAL(DP), ALLOCATABLE :: h_diag (:,:), eprec(:)
REAL(DP), ALLOCATABLE :: h_diag (:,:)
! diagonal part of h_scf
real(DP) :: anorm
! preconditioning cut-off
@ -69,17 +70,12 @@ SUBROUTINE lr_dvpsi_e(ik,ipol,dvpsi)
!
CALL start_clock ('lr_dvpsi_e')
!
IF (lr_verbosity > 5) WRITE(stdout,'("<lr_dvpsi_e>")')
!
conv_root = .TRUE.
!
ALLOCATE(d0psi(npwx*npol,nbnd))
ALLOCATE ( d0psi(npwx*npol,nbnd) )
d0psi = (0.d0, 0.d0)
dvpsi = (0.d0, 0.d0)
!
ALLOCATE (h_diag( npwx*npol, nbnd))
h_diag = 0.d0
!
npw = ngk(ik)
!
CALL allocate_bec_type ( nkb, nbnd, becp1 )
@ -97,53 +93,35 @@ SUBROUTINE lr_dvpsi_e(ik,ipol,dvpsi)
CALL orthogonalize(d0psi, evc, ik, ik, sevc0(:,:,ik), npw, .true.)
d0psi = -d0psi
!
! Calculate the kinetic energy g2kin
! Calculate the kinetic energy g2kin: (k+G)^2
!
DO ig = 1, npw
g2kin(ig) = SUM((xk(1:3,ik) + g(1:3,igk_k(ig,ik)))**2) * tpiba2
ENDDO
CALL g2_kin(ik)
!
! d0psi contains P^+_c [H-eS,x] psi_v for the polarization direction ipol
! Now solve the linear systems (H-e_vS)*P_c(x*psi_v)=P_c^+ [H-e_vS,x]*psi_v
! Calculate the preconditioning matrix h_diag used by cgsolve_all
!
! eprec is now calculated on the fly for each k point
ALLOCATE ( h_diag(npwx*npol, nbnd) )
CALL h_prec(ik, evc, h_diag)
!
ALLOCATE(eprec(nbnd))
CALL lr_calc_eprec(eprec)
!
DO ibnd = 1, nbnd_occ (ik)
DO ig = 1, npw
h_diag (ig, ibnd) = 1.d0 / max (1.0d0, g2kin (ig) / eprec (ibnd) )
ENDDO
IF (noncolin) THEN
DO ig = 1, npw
h_diag (ig+npwx, ibnd) = 1.d0/max(1.0d0,g2kin(ig)/eprec(ibnd))
ENDDO
ENDIF
ENDDO
! d0psi contains P^+_c [H-eS,x] psi_v for the polarization direction ipol
! Now solve the linear systems (H+Q-e_vS)*P_c(x*psi_v)=P_c^+ [H-e_vS,x]*psi_v
! See Eq.(9) in Ref. [1]
!
CALL cgsolve_all (ch_psi_all, cg_psi, et (1, ik), d0psi, dvpsi, &
h_diag, npwx, npw, thresh, ik, lter, conv_root, anorm, &
nbnd_occ(ik), 1)
h_diag, npwx, npw, thresh, ik, lter, conv_root, anorm, nbnd_occ(ik), 1)
!
IF (.not.conv_root) WRITE( stdout, '(5x,"ik",i4," ibnd",i4, &
& " linter: root not converged ",e10.3)') &
& " lr_dvpsi_e: root not converged ",e10.3)') &
ik, ibnd, anorm
!
FLUSH( stdout )
!
! we have now obtained P_c x |psi>.
! In the case of USPP this quantity is needed for the Born
! effective charges, so we save it to disc
DEALLOCATE (h_diag)
!
! In the US case we obtain P_c x |psi>, but we need P_c^+ x | psi>,
! therefore we apply S again, and then subtract the additional term
! furthermore we add the term due to dipole of the augmentation charges.
! See Eq.(10) in Ref. [1]
!
IF (okvan) THEN
!
! for effective charges
!
ALLOCATE (spsi ( npwx*npol, nbnd))
CALL calbec (npw, vkb, dvpsi, becp )
CALL s_psi(npwx,npw,nbnd,dvpsi,spsi)
@ -154,7 +132,6 @@ SUBROUTINE lr_dvpsi_e(ik,ipol,dvpsi)
CALL qdipol_cryst()
CALL adddvepsi_us(becp1,becp2,ipol,ik,dvpsi)
DEALLOCATE (dpqq)
!
ENDIF
!
IF (okvan) CALL calbec ( npw, vkb, evc, becp, nbnd)
@ -164,11 +141,9 @@ SUBROUTINE lr_dvpsi_e(ik,ipol,dvpsi)
CALL orthogonalize(dvpsi, evc, ik, ik, sevc0(:,:,ik), npw, .true.)
dvpsi = -dvpsi
!
DEALLOCATE (h_diag)
DEALLOCATE (eprec)
DEALLOCATE (d0psi)
!
! OBM: Addendum to PH dvpsi
! US case: apply the S^{-1} operator
!
IF (okvan) THEN
ALLOCATE (spsi ( npwx*npol, nbnd))
@ -181,13 +156,9 @@ SUBROUTINE lr_dvpsi_e(ik,ipol,dvpsi)
! Here we include the correct normalization
! for Lanczos initial wfcs
!
atnorm = dsqrt(at(1,ipol)**2+at(2,ipol)**2+at(3,ipol)**2)
atnorm = DSQRT(at(1,ipol)**2 + at(2,ipol)**2 + at(3,ipol)**2)
!
dvpsi(:,:) = dvpsi(:,:)/atnorm
!
! nrec = (ipol - 1)*nksq + ik
! call davcio(dvpsi, lrebar, iuebar, nrec, 1)
! this_pcxpsi_is_on_file(ik,ipol) = .true.
dvpsi(:,:) = dvpsi(:,:) / atnorm
!
CALL deallocate_bec_type ( becp1 )
IF (nkb > 0) CALL deallocate_bec_type ( becp2 )
@ -196,62 +167,6 @@ SUBROUTINE lr_dvpsi_e(ik,ipol,dvpsi)
!
RETURN
!
CONTAINS
SUBROUTINE lr_calc_eprec(eprec)
USE kinds, ONLY : DP
USE gvect, ONLY : gstart
USE wvfct, ONLY : npwx, nbnd, g2kin
USE wavefunctions_module, ONLY : evc
USE klist, ONLY : ngk
USE mp, ONLY : mp_sum
USE mp_global, ONLY : intra_bgrp_comm
IMPLICIT NONE
REAL(KIND=DP), INTENT(INOUT) :: eprec(nbnd)
COMPLEX(KIND=DP), ALLOCATABLE :: work (:,:)
REAL(KIND=DP), EXTERNAL :: ddot
COMPLEX(KIND=DP), EXTERNAL :: ZDOTC
! the scalar products
!
ALLOCATE (work(npwx,nbnd))
!
DO ibnd=1,nbnd
!
work = 0.d0
!
DO ig = 1, ngk(ik)
work(ig,1) = g2kin(ig)*evc(ig,ibnd)
ENDDO
!
IF (gamma_only) THEN
!
eprec(ibnd) = 2.0d0*DDOT(2*npw,evc(1,ibnd),1,work,1)
!
IF (gstart==2) THEN
eprec(ibnd) = eprec(ibnd)-DBLE(evc(1,ibnd))*DBLE(work(1,ibnd))
ENDIF
!
eprec(ibnd) = 1.35d0*eprec(ibnd)
!
ELSE
eprec(ibnd) = 1.35d0*ZDOTC(ngk(ik),evc(1,ibnd),1,work,1)
ENDIF
!
ENDDO
!
#ifdef __MPI
CALL mp_sum(eprec, intra_bgrp_comm)
#endif
!
DEALLOCATE(work)
!
RETURN
!
END SUBROUTINE lr_calc_eprec
END SUBROUTINE lr_dvpsi_e

8
TDDFPT/src/lr_init_nfo.f90
View File

@ -80,12 +80,12 @@ SUBROUTINE lr_init_nfo()
!
nksq = nks
!
! Note: in this case the array ikqs is
! needed only in order to use the routine
! ch_psi_all
! Note: in this case the arrays ikks and ikqs
! are needed only in h_prec and ch_psi_all.
!
ALLOCATE(ikqs(nksq))
ALLOCATE(ikks(nksq), ikqs(nksq))
DO ik = 1, nksq
ikks(ik) = ik
ikqs(ik) = ik
ENDDO
!

29
TDDFPT/src/lr_restart.f90
View File

@ -18,15 +18,14 @@ SUBROUTINE lr_restart(iter_restart,rflag)
USE io_global, ONLY : stdout, ionode_id
USE control_flags, ONLY : gamma_only
USE klist, ONLY : nks, xk, ngk, igk_k
USE cell_base, ONLY : tpiba2
USE gvect, ONLY : g
USE io_files, ONLY : tmp_dir, prefix, diropn, wfc_dir
USE lr_variables, ONLY : itermax, evc1, evc1_old, &
restart, nwordrestart, iunrestart,project,nbnd_total,F, &
bgz_suffix, beta_store, gamma_store, zeta_store, norm0, &
lr_verbosity, charge_response, LR_polarization, n_ipol, eels, sum_rule
lr_verbosity, charge_response, LR_polarization, n_ipol, &
eels, sum_rule
USE charg_resp, ONLY : resonance_condition, rho_1_tot,rho_1_tot_im
USE wvfct, ONLY : nbnd, g2kin, npwx
USE wvfct, ONLY : nbnd
USE becmod, ONLY : bec_type, becp, calbec
USE uspp, ONLY : vkb
USE io_global, ONLY : ionode
@ -37,9 +36,10 @@ SUBROUTINE lr_restart(iter_restart,rflag)
IMPLICIT NONE
!
INTEGER, INTENT(OUT) :: iter_restart
LOGICAL, INTENT(OUT) :: rflag
!
CHARACTER(len=6), EXTERNAL :: int_to_char
INTEGER, INTENT(out) :: iter_restart
LOGICAL, INTENT(out) :: rflag
!
! local variables
!
@ -61,20 +61,15 @@ SUBROUTINE lr_restart(iter_restart,rflag)
!
rflag = .false.
!
! Optical case: restarting kintic-energy and ultrasoft
! Optical case: recompute the kintic-energy g2kin and
! beta functions vkb (needed only in the US case).
! Note, this is done only in the gamma_only case,
! because in the k-points version all is recomputed
! on-the-fly for every k point.
!
IF (gamma_only) THEN
!
DO ig=1,npwx
!
g2kin(ig)=((xk(1,1)+g(1,igk_k(ig,1)))**2 &
+ (xk(2,1)+g(2,igk_k(ig,1)))**2 &
+ (xk(3,1)+g(3,igk_k(ig,1)))**2)*tpiba2
!
ENDDO
!
CALL g2_kin(1)
CALL init_us_2(ngk(1),igk_k(:,1),xk(:,1),vkb)
!
ENDIF
!
! Reading Lanczos coefficients

274
TDDFPT/src/lr_solve_e.f90
View File

@ -50,16 +50,11 @@ SUBROUTINE lr_solve_e
real (kind=dp) :: anorm
CHARACTER(len=256) :: tmp_dir_saved
!
IF (lr_verbosity > 5) WRITE(stdout,'("<lr_solve_e>")')
!
CALL start_clock ('lr_solve_e')
!
IF( lr_verbosity > 1 ) &
WRITE(stdout,'(5X,"lr_solve_e: alpha_pv=",1X,e12.5)') alpha_pv
!
IF (eels) THEN
!
! EELS
! EELS case
!
DO ik = 1, nksq
CALL lr_dvpsi_eels(ik, d0psi(:,:,ik,1), d0psi2(:,:,ik,1))
@ -69,50 +64,67 @@ SUBROUTINE lr_solve_e
!
! Optical case
!
! TDDFPT+hybrids: Current implementation of the commutator [H,r]
! (commutator_Hx_psi) does not contain the contribution [V_EXX,r],
! hence the relative intensities of the peaks in the spectrum
! will be wrong. One should use d0psi_rs=.true. in this case
! (see the documentation and the explanation in CPC 185, 2080 (2014)).
!
!
IF (lr_exx .AND. .NOT.d0psi_rs) WRITE( stdout, '(/5X,"WARNING!!!", &
& " TDDFPT + hybrids will give wrong intensities.", &
& /5x,"Set d0psi_rs=.true. in order to have correct intensities.")' )
!
DO ik = 1, nks
IF (.NOT.d0psi_rs) THEN
!
current_k = ik
! Compute dipole in the G space using the commutator [H,r]
!
IF ( lsda ) current_spin = isk(ik)
! TDDFPT+hybrids: Current implementation of the commutator [H,r]
! (commutator_Hx_psi) does not contain the contribution [V_EXX,r],
! hence the relative intensities of the peaks in the spectrum
! will be wrong. One should use d0psi_rs=.true. in this case
! (see the documentation and the explanation in CPC 185, 2080 (2014)).
!
IF (lr_exx) WRITE( stdout, '(/5X,"WARNING!!!", &
& " TDDFPT + hybrids will give wrong intensities.", &
& /5x,"Set d0psi_rs=.true. in order to have correct intensities.")' )
!
evc(:,:) = evc0(:,:,ik)
DO ik = 1, nks
!
current_k = ik
!
IF ( lsda ) current_spin = isk(ik)
!
evc(:,:) = evc0(:,:,ik)
!
! US case: calculate beta-functions vkb.
!
CALL init_us_2(ngk(ik), igk_k(:,ik), xk(:,ik), vkb)
!
! Compute d0psi = P_c^+ r psi_k
!
IF ( n_ipol==3 ) THEN
DO ip=1,3
CALL lr_dvpsi_e(ik,ip,d0psi(:,:,ik,ip))
ENDDO
ELSEIF ( n_ipol==1 ) THEN
CALL lr_dvpsi_e(ik,LR_polarization,d0psi(:,:,ik,1))
ENDIF
!
ENDDO
!
! Ultrasoft case: calculate beta-functions vkb.
ELSE
!
CALL init_us_2(ngk(ik), igk_k(:,ik), xk(:,ik), vkb)
! Compute dipole in the R space. This option can be used
! only for finite systems (e.g. molecules).
!
! Computes/reads P_c^+ x psi_kpoint into d0psi array
CALL compute_d0psi_rs(n_ipol)
!
IF ( n_ipol==3 ) THEN
DO ip=1,3
CALL lr_dvpsi_e(ik,ip,d0psi(:,:,ik,ip))
ENDDO
ELSEIF ( n_ipol==1 ) THEN
CALL lr_dvpsi_e(ik,LR_polarization,d0psi(:,:,ik,1))
! In the gamma_only case we compute the kinetic energy
! and the beta functions vkb here just once, and then
! they are used throughout the calculation. In the
! k-points version these quantities are always recomputed.
!
IF (gamma_only) THEN
CALL g2_kin(1)
CALL init_us_2(ngk(1), igk_k(:,1), xk(:,1), vkb)
ENDIF
!
ENDDO
!
ENDIF
!
ENDIF
!
IF (gstart == 2 .and. gamma_only) d0psi(1,:,:,:) = cmplx(dble(d0psi(1,:,:,:)),0.0d0,dp)
!
! X. Ge: compute d0psi in real-space, will overwrite d0psi calculated before.
!
IF (.NOT. eels) THEN
IF (d0psi_rs) CALL compute_d0psi_rs(n_ipol)
ENDIF
IF (gstart==2 .and. gamma_only) d0psi(1,:,:,:) = &
& CMPLX(DBLE(d0psi(1,:,:,:)),0.0d0,dp)
!
! Writing of d0psi to the file.
! This is a parallel writing, done in wfc_dir
@ -123,8 +135,10 @@ SUBROUTINE lr_solve_e
!
DO ip = 1, n_ipol
!
IF (n_ipol==1) CALL diropn ( iund0psi, 'd0psi.'//trim(int_to_char(LR_polarization)), nwordd0psi, exst)
IF (n_ipol==3) CALL diropn ( iund0psi, 'd0psi.'//trim(int_to_char(ip)), nwordd0psi, exst)
IF (n_ipol==1) CALL diropn ( iund0psi, 'd0psi.'// &
& trim(int_to_char(LR_polarization)), nwordd0psi, exst)
IF (n_ipol==3) CALL diropn ( iund0psi, 'd0psi.'// &
& trim(int_to_char(ip)), nwordd0psi, exst)
!
CALL davcio(d0psi(1,1,1,ip),nwordd0psi,iund0psi,1,1)
!
@ -136,7 +150,8 @@ SUBROUTINE lr_solve_e
!
IF (eels) THEN
!
CALL diropn ( iund0psi, 'd0psi2.'//trim(int_to_char(LR_polarization)), nwordd0psi, exst)
CALL diropn ( iund0psi, 'd0psi2.'// &
& trim(int_to_char(LR_polarization)), nwordd0psi, exst)
CALL davcio(d0psi2(1,1,1,1),nwordd0psi,iund0psi,1,1)
CLOSE( UNIT = iund0psi)
!
@ -152,9 +167,10 @@ CONTAINS
SUBROUTINE compute_d0psi_rs( n_ipol )
!
! ... original code from routine compute_dipole
! ... Original code from routine compute_dipole
! ... modified to calculate the d0psi in the real space
! ... by Xiaochuan Ge, Oct, 2013
! ... Modified by I. Timrov, June 2016
!
USE kinds, ONLY : DP
USE cell_base, ONLY : at, bg, alat, omega
@ -167,39 +183,39 @@ SUBROUTINE compute_d0psi_rs( n_ipol )
USE lr_variables, ONLY : evc0,sevc0,d0psi,lshift_d0psi
USE wavefunctions_module, ONLY : psic
USE uspp, ONLY : okvan
USE realus, ONLY : fwfft_orbital_gamma,invfft_orbital_gamma
USE realus, ONLY : fwfft_orbital_gamma, invfft_orbital_gamma, &
fwfft_orbital_k, invfft_orbital_k
!
IMPLICIT NONE
!
! ... Define variables
complex(dp) :: wfck(npwx,1)
! ... Local variables
REAL(DP),allocatable :: r(:,:)
complex(dp),allocatable :: psic_temp(:)
INTEGER :: i, j, k, ip, ir, ir_end, index0,ib,n_ipol
!
REAL(DP), ALLOCATABLE :: r(:,:)
COMPLEX(DP), ALLOCATABLE :: psic_temp(:), evc_temp(:,:)
INTEGER :: i, j, k, ip, ir, ir_end, index0, ibnb, n_ipol
REAL(DP) :: inv_nr1, inv_nr2, inv_nr3
!
IF (.not. ALLOCATED(psic)) ALLOCATE(psic(dfftp%nnr))
WRITE(stdout,'(/,5X,"Calculation of the dipole in real space")')
!
IF (.NOT. ALLOCATED(psic)) ALLOCATE(psic(dfftp%nnr))
ALLOCATE(evc_temp(npwx,nbnd))
ALLOCATE(psic_temp(dfftp%nnr))
ALLOCATE(r(dfftp%nnr,n_ipol))
!
! ... Initialization
WRITE(stdout,'(/,5x,"Calculating d0psi in the real space.")')
evc_temp(:,:) = (0.0d0, 0.0d0)
psic_temp(:) = (0.0d0, 0.0d0)
r(:,:) = 0.0d0
!
IF (okvan) THEN
!
WRITE(stdout,'(10x,"At this moment d0psi_rs is not available for USPP !!!",//)')
WRITE(stdout,'(10x,"Real space dipole + USPP is not supported",/)')
#ifdef __MPI
call mp_barrier(intra_bgrp_comm)
! call mp_stop(100)
CALL mp_barrier(intra_bgrp_comm)
#endif
stop
endif
STOP
!
ENDIF
!
IF (nks > 1) CALL errore( 'compute_d0psi_rs', 'nks>1 is not supported', 1 )
!
! Calculate r
! Calculate r
!
inv_nr1 = 1.D0 / DBLE( dfftp%nr1 )
inv_nr2 = 1.D0 / DBLE( dfftp%nr2 )
@ -231,43 +247,99 @@ SUBROUTINE compute_d0psi_rs( n_ipol )
!
ENDDO
!
If (lshift_d0psi) call shift_d0psi(r,n_ipol)
! Shift the dipole to the center of cell
!
DO ib = 1, nbnd
IF (lshift_d0psi) CALL shift_d0psi(r,n_ipol)
!
! Calculate the product r * psi(r)
!
IF (gamma_only) THEN
!
wfck(:,1) = evc0(:,ib,1)
! gamma_only version
! Note: two bands are treated at the same time.
!
CALL invfft_orbital_gamma(wfck(:,:),1,1)
!
psic_temp(:) = psic(:)
!
DO ip = 1, n_ipol
!
! Apply the dipole operator
!
DO ir = 1, dfftp%nnr
psic(ir) = r(ir,ip) * alat * psic_temp(ir)
ENDDO
!
! Convert to G-space
!
CALL fwfft_orbital_gamma(wfck(:,:),1,1)
!
d0psi(:,ib,1,ip) = wfck(:,1)
!
DO ibnd = 1, nbnd, 2
!
! FFT to R-space: evc_temp -> psic
!
CALL invfft_orbital_gamma(evc0(:,:,1), ibnd, nbnd)
!
psic_temp(:) = psic(:)
!
DO ip = 1, n_ipol
!
! Multiple in R space the dipole and
! the ground state wavefunction
!
DO ir = 1, dfftp%nnr
psic(ir) = r(ir,ip) * alat * psic_temp(ir)
ENDDO
!
! FFT to G space: psic -> evc_temp
!
CALL fwfft_orbital_gamma(evc_temp, ibnd, nbnd)
!
d0psi(:,ibnd, 1,ip) = evc_temp(:,ibnd)
d0psi(:,ibnd+1,1,ip) = evc_temp(:,ibnd+1)
!
ENDDO
!
ENDDO
!
ENDDO
ELSE
!
! k-points version
! Note: one band is treated at a time
!
DO ik = 1, nks
!
DO ibnd = 1, nbnd
!
! FFT to R-space: evc0 -> psic
!
CALL invfft_orbital_k(evc0(:,:,ik), ibnd, nbnd, ik)
!
psic_temp(:) = psic(:)
!
DO ip = 1, n_ipol
!
! Multiple in R space the dipole and
! the ground state wavefunction
!
DO ir = 1, dfftp%nnr
psic(ir) = r(ir,ip) * alat * psic_temp(ir)
ENDDO
!
! FFT to G space: psic -> evc_temp
!
CALL fwfft_orbital_k(evc_temp, ibnd, nbnd, ik)
!
d0psi(:,ibnd,ik,ip) = evc_temp(:,ibnd)
!
ENDDO
!
ENDDO
!
ENDDO
!
ENDIF
!
! Orthogonalized batch orbitals to occupied minifold
! Orthogonalization: apply P_c to d0psi
! Output: d0psi = P_c r * psi(r)
!
DO ip = 1, n_ipol
CALL orthogonalize(d0psi(:,:,1,ip), evc0(:,:,1), 1, 1, sevc0(:,:,1), ngk(1), .true.)
d0psi(:,:,1,ip) = -d0psi(:,:,1,ip)
DO ik = 1, nks
!
CALL orthogonalize(d0psi(:,:,ik,ip), evc0(:,:,ik), &
& ik, ik, sevc0(:,:,ik), ngk(ik), .true.)
d0psi(:,:,ik,ip) = -d0psi(:,:,ik,ip)
!
ENDDO
ENDDO
!
DEALLOCATE(r)
DEALLOCATE(evc_temp)
DEALLOCATE(psic_temp)
DEALLOCATE(r)
!
RETURN
!
@ -283,37 +355,43 @@ SUBROUTINE shift_d0psi( r, n_ipol )
use ions_base, only : nat, tau
USE io_global, ONLY : stdout
use cell_base, only : alat, at
!
IMPLICIT NONE
REAL(dp),intent(inout) :: r(dfftp%nnr,n_ipol)
integer,intent(in) :: n_ipol
real(dp) :: mmin(3), mmax(3), center(3), origin(3), check_cell
integer :: ip, iatm, ir, ip2
WRITE(stdout,'(/,5X,"Wavefunctions are shifted to the center of cell for the calculation of d0psi. ")')
!
WRITE(stdout,'(/,5X,"Dipole is shifted to the center of cell", &
& " for the calculation of d0psi. ")')
!
check_cell = 0.d0
do ip = 1, 3
do ip2 = 1,3
if(.not. ip .eq. ip2) check_cell = check_cell + at(ip,ip2)**2
enddo
enddo
if(check_cell .gt. 1.d-5) call errore('lr_read_wfc'," I am not sure that this type of super cell is supported now. ",1)
!
! XG: I am not sure that this type of super cell is supported now.
!
if (check_cell .gt. 1.d-5) call errore('lr_read_wfc', &
& "This type of the supercell is not supported",1)
!
mmin(:) = 2000.d0
mmax(:)= -2000.d0
!
do ip = 1, n_ipol
do iatm = 1, nat
mmin(ip) = min(mmin(ip), tau(ip,iatm))
mmax(ip) = max(mmax(ip), tau(ip,iatm))
enddo
enddo
!
center(:)= 0.5d0*(mmin(:)+mmax(:))
do ip = 1, n_ipol
origin(ip)= center(ip)-0.5d0*at(ip,ip)
enddo
!
do ir = 1, dfftp%nnr
r(ir,:)= r(ir,:) - origin(:)
do ip = 1, n_ipol
@ -321,9 +399,9 @@ SUBROUTINE shift_d0psi( r, n_ipol )
if(r(ir,ip) .gt. at(ip,ip)) r(ir,ip)=r(ir,ip)-at(ip,ip)
enddo
enddo
!
RETURN
!
END SUBROUTINE shift_d0psi
END SUBROUTINE lr_solve_e

42
TDDFPT/src/make.depend
View File

@ -64,20 +64,14 @@ lr_apply_liouvillian.o : ../../PW/src/realus.o
lr_apply_liouvillian.o : ../../PW/src/scf_mod.o
lr_apply_liouvillian.o : lr_exx_kernel.o
lr_apply_liouvillian.o : lr_variables.o
lr_apply_liouvillian_eels.o : ../../FFTXlib/fft_interfaces.o
lr_apply_liouvillian_eels.o : ../../FFTXlib/fft_parallel.o
lr_apply_liouvillian_eels.o : ../../LR_Modules/dv_of_drho.o
lr_apply_liouvillian_eels.o : ../../LR_Modules/lrcom.o
lr_apply_liouvillian_eels.o : ../../Modules/cell_base.o
lr_apply_liouvillian_eels.o : ../../Modules/fft_base.o
lr_apply_liouvillian_eels.o : ../../Modules/gvecw.o
lr_apply_liouvillian_eels.o : ../../Modules/io_files.o
lr_apply_liouvillian_eels.o : ../../Modules/io_global.o
lr_apply_liouvillian_eels.o : ../../Modules/ions_base.o
lr_apply_liouvillian_eels.o : ../../Modules/kind.o
lr_apply_liouvillian_eels.o : ../../Modules/mp_bands.o
lr_apply_liouvillian_eels.o : ../../Modules/noncol.o
lr_apply_liouvillian_eels.o : ../../Modules/recvec.o
lr_apply_liouvillian_eels.o : ../../Modules/uspp.o
lr_apply_liouvillian_eels.o : ../../Modules/wavefunctions.o
lr_apply_liouvillian_eels.o : ../../PW/src/buffers.o
@ -249,8 +243,6 @@ lr_dvpsi_e.o : ../../Modules/control_flags.o
lr_dvpsi_e.o : ../../Modules/io_global.o
lr_dvpsi_e.o : ../../Modules/ions_base.o
lr_dvpsi_e.o : ../../Modules/kind.o
lr_dvpsi_e.o : ../../Modules/mp.o
lr_dvpsi_e.o : ../../Modules/mp_global.o
lr_dvpsi_e.o : ../../Modules/noncol.o
lr_dvpsi_e.o : ../../Modules/recvec.o
lr_dvpsi_e.o : ../../Modules/uspp.o
@ -441,7 +433,6 @@ lr_readin.o : lr_charg_resp.o
lr_readin.o : lr_dav_variables.o
lr_readin.o : lr_variables.o
lr_restart.o : ../../Modules/becmod.o
lr_restart.o : ../../Modules/cell_base.o
lr_restart.o : ../../Modules/control_flags.o
lr_restart.o : ../../Modules/fft_base.o
lr_restart.o : ../../Modules/io_files.o
@ -450,7 +441,6 @@ lr_restart.o : ../../Modules/kind.o
lr_restart.o : ../../Modules/mp.o
lr_restart.o : ../../Modules/mp_world.o
lr_restart.o : ../../Modules/noncol.o
lr_restart.o : ../../Modules/recvec.o
lr_restart.o : ../../Modules/uspp.o
lr_restart.o : ../../PW/src/pwcom.o
lr_restart.o : lr_charg_resp.o
@ -505,6 +495,38 @@ lr_solve_e.o : ../../Modules/wavefunctions.o
lr_solve_e.o : ../../PW/src/pwcom.o
lr_solve_e.o : ../../PW/src/realus.o
lr_solve_e.o : lr_variables.o
lr_solve_e_cp.o : ../../LR_Modules/lrcom.o
lr_solve_e_cp.o : ../../Modules/cell_base.o
lr_solve_e_cp.o : ../../Modules/control_flags.o
lr_solve_e_cp.o : ../../Modules/fft_base.o
lr_solve_e_cp.o : ../../Modules/io_files.o
lr_solve_e_cp.o : ../../Modules/io_global.o
lr_solve_e_cp.o : ../../Modules/ions_base.o
lr_solve_e_cp.o : ../../Modules/kind.o
lr_solve_e_cp.o : ../../Modules/mp.o
lr_solve_e_cp.o : ../../Modules/mp_global.o
lr_solve_e_cp.o : ../../Modules/recvec.o
lr_solve_e_cp.o : ../../Modules/uspp.o
lr_solve_e_cp.o : ../../Modules/wavefunctions.o
lr_solve_e_cp.o : ../../PW/src/pwcom.o
lr_solve_e_cp.o : ../../PW/src/realus.o
lr_solve_e_cp.o : lr_variables.o
lr_solve_e_cp2.o : ../../LR_Modules/lrcom.o
lr_solve_e_cp2.o : ../../Modules/cell_base.o
lr_solve_e_cp2.o : ../../Modules/control_flags.o
lr_solve_e_cp2.o : ../../Modules/fft_base.o
lr_solve_e_cp2.o : ../../Modules/io_files.o
lr_solve_e_cp2.o : ../../Modules/io_global.o
lr_solve_e_cp2.o : ../../Modules/ions_base.o
lr_solve_e_cp2.o : ../../Modules/kind.o
lr_solve_e_cp2.o : ../../Modules/mp.o
lr_solve_e_cp2.o : ../../Modules/mp_global.o
lr_solve_e_cp2.o : ../../Modules/recvec.o
lr_solve_e_cp2.o : ../../Modules/uspp.o
lr_solve_e_cp2.o : ../../Modules/wavefunctions.o
lr_solve_e_cp2.o : ../../PW/src/pwcom.o
lr_solve_e_cp2.o : ../../PW/src/realus.o
lr_solve_e_cp2.o : lr_variables.o
lr_summary.o : ../../LR_Modules/lrcom.o
lr_summary.o : ../../Modules/cell_base.o
lr_summary.o : ../../Modules/fft_base.o