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

General clean-up of TDDFT routines. Some work courtesy of Iurii Timrov 

git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@9439 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
sbinnie 2012-09-20 09:54:04 +00:00
parent 287793385e
commit 41246be616
4 changed files with 935 additions and 854 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

218
TDDFPT/src/lr_apply_liouvillian.f90
View File

@ -1,10 +1,14 @@
!-----------------------------------------------------------------------
SUBROUTINE lr_apply_liouvillian( evc1, evc1_new, sevc1_new, interaction )
!---------------------------------------------------------------------
! ... applies linear response operator to response wavefunctions
! OBM: or to be more exact this function is responsible for calculating L.q(i) and (L^T).p(i)
! where q is evc1 and return partial qdash(i+1) and pdash(i+1) in evc1_new . Ultrasoft additions are
! handled here...
! Applies linear response operator to response wavefunctions
! (H - E)*psi(k+q) + HXC
!
! Or to be more exact this routine is responsible for calculating
! L.q(i) and (L^T).p(i)
! Where q is evc1
! It returns partial qdash(i+1) and pdash(i+1) in evc1_new.
!
! interaction=.true. corresponds to eq.(32)
! interaction=.false. corresponds to eq.(33)
! in Ralph Gebauer, Brent Walker J. Chem. Phys., 127, 164106 (2007)
@ -16,13 +20,14 @@ SUBROUTINE lr_apply_liouvillian( evc1, evc1_new, sevc1_new, interaction )
USE cell_base, ONLY : tpiba2
USE fft_base, ONLY : dffts
USE fft_interfaces, ONLY : fwfft
USE gvecs, ONLY : nls, nlsm
USE gvecs, ONLY : nls, nlsm
USE gvect, ONLY : nl, ngm, gstart, g, gg
USE fft_base, ONLY : dfftp
USE io_global, ONLY : stdout
USE kinds, ONLY : dp
USE klist, ONLY : nks, xk
USE lr_variables, ONLY : evc0, revc0, rho_1, lr_verbosity, ltammd, size_evc, no_hxc
USE lr_variables, ONLY : evc0, revc0, rho_1, lr_verbosity,&
& ltammd, size_evc, no_hxc
USE realus, ONLY : igk_k,npw_k
USE lsda_mod, ONLY : nspin
USE uspp, ONLY : vkb, nkb, okvan
@ -30,91 +35,128 @@ SUBROUTINE lr_apply_liouvillian( evc1, evc1_new, sevc1_new, interaction )
USE wavefunctions_module, ONLY : psic
USE wvfct, ONLY : nbnd, npwx, igk, g2kin, et
USE control_flags, ONLY : gamma_only
USE realus, ONLY : real_space, fft_orbital_gamma, initialisation_level, &
bfft_orbital_gamma, calbec_rs_gamma, add_vuspsir_gamma, &
v_loc_psir, s_psir_gamma, real_space_debug, &
betasave, box_beta, maxbox_beta,newq_r
USE lr_variables, ONLY : lr_verbosity
USE io_global, ONLY : stdout
USE dfunct, ONLY : newq
USE control_flags, ONLY : tqr
USE realus, ONLY : real_space, fft_orbital_gamma,&
& initialisation_level,&
& bfft_orbital_gamma,&
& calbec_rs_gamma,&
& add_vuspsir_gamma, v_loc_psir,&
& s_psir_gamma, real_space_debug,&
& betasave, box_beta, maxbox_beta,&
& newq_r
USE lr_variables, ONLY : lr_verbosity
USE io_global, ONLY : stdout
USE dfunct, ONLY : newq
USE control_flags, ONLY : tqr
!
!
IMPLICIT NONE
!
COMPLEX(kind=dp),INTENT(in) :: evc1(npwx,nbnd,nks)
COMPLEX(kind=dp),INTENT(out) :: evc1_new(npwx,nbnd,nks), sevc1_new(npwx,nbnd,nks)
COMPLEX(kind=dp),INTENT(out) :: evc1_new(npwx,nbnd,nks),&
& sevc1_new(npwx,nbnd,nks)
! output : sevc1_new = S * evc1_new
!
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(:)
COMPLEX(kind=dp), ALLOCATABLE :: dvrs_temp(:,:) !OBM This waste of memory was already there in lr_dv_of_drho
real(kind=dp), ALLOCATABLE :: d_deeq(:,:,:,:)
!
REAL(kind=dp), ALLOCATABLE :: dvrs(:,:), dvrss(:)
REAL(kind=dp), ALLOCATABLE :: d_deeq(:,:,:,:)
!
COMPLEX(kind=dp), ALLOCATABLE :: dvrs_temp(:,:)
COMPLEX(kind=dp), ALLOCATABLE :: spsi1(:,:)
!
COMPLEX(kind=dp) :: fp, fm
REAL(DP), ALLOCATABLE, DIMENSION(:) :: w1, w2
!
REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: w1, w2
!
!
IF (lr_verbosity > 5) THEN
WRITE(stdout,'("<lr_apply_liouvillian>")')
ENDIF
!
CALL start_clock('lr_apply')
IF (interaction) CALL start_clock('lr_apply_int')
!
IF (interaction) CALL start_clock('lr_apply_int')
IF (.not.interaction) CALL start_clock('lr_apply_no')
!
ALLOCATE( d_deeq(nhm, nhm, nat, nspin) )
d_deeq(:,:,:,:)=0.0d0
!
ALLOCATE( spsi1(npwx, nbnd) )
spsi1(:,:)=(0.0d0,0.0d0)
!
IF( interaction ) THEN !If true, the full L is calculated
IF( interaction ) THEN
!
! Calculate the full L
!
ALLOCATE( dvrs(dfftp%nnr, nspin) )
ALLOCATE( dvrss(dffts%nnr) )
dvrs(:,:)=0.0d0
dvrss(:)=0.0d0
!
! Calculation of the charge density response
!
CALL lr_calc_dens( evc1, .false. )
!
! Given the change of the charge density
! we calculate the change of the Hartree and XC potential and
! put it in dvrss()
!
IF (no_hxc) THEN
!OBM no_hxc controls the hartree excange correlation addition, if true, they are not added
dvrs(:,1)=0.0d0
CALL interpolate (dvrs(:,1),dvrss,-1)
!
! With no_hxc=.true. we recover the independent electrion
! approximation, so we zero the interation.
!
dvrs(:,1)=0.0d0
CALL interpolate (dvrs(:,1),dvrss,-1)
!
ELSE
dvrs(:,1)=rho_1(:,1)
!
!call lr_dv_of_drho(dvrs)
ALLOCATE( dvrs_temp(dfftp%nnr, nspin) )
dvrs_temp = CMPLX( dvrs, 0.0d0, kind=dp )
! OBM: This memory copy was hidden in lr_dv_of_drho, can it be avoided?
DEALLOCATE ( dvrs )
CALL dv_of_drho(0,dvrs_temp,.false.)
ALLOCATE ( dvrs(dfftp%nnr, nspin) ) !SJB Worth getting rid of this memory bottle neck for the moment.
dvrs=dble(dvrs_temp)
DEALLOCATE(dvrs_temp)
!
IF ( okvan ) THEN
IF ( tqr ) THEN
CALL newq_r(dvrs,d_deeq,.true.)
ELSE
ALLOCATE( psic(dfftp%nnr) )
psic(:)=(0.0d0,0.0d0)
CALL newq(dvrs,d_deeq,.true.)
DEALLOCATE( psic )
!
dvrs(:,1)=rho_1(:,1)
!
! In the gamma_only case dvrs is real, but dv_of_drho expects
! a complex array on input, hence this temporary variable.
!
ALLOCATE( dvrs_temp(dfftp%nnr, nspin) )
dvrs_temp = CMPLX( dvrs, 0.0d0, kind=dp )
!
DEALLOCATE ( dvrs )
!
CALL dv_of_drho(0,dvrs_temp,.FALSE.)
!
ALLOCATE ( dvrs(dfftp%nnr, nspin) )
dvrs=DBLE(dvrs_temp)
DEALLOCATE(dvrs_temp)
!
IF ( okvan ) THEN
IF ( tqr ) THEN
CALL newq_r(dvrs,d_deeq,.TRUE.)
ELSE
ALLOCATE( psic(dfftp%nnr) )
psic(:)=(0.0d0,0.0d0)
!
CALL newq(dvrs,d_deeq,.TRUE.)
!
DEALLOCATE( psic )
ENDIF
ENDIF
ENDIF
CALL add_paw_to_deeq(d_deeq)
!
CALL interpolate (dvrs(:,1),dvrss,-1)
CALL add_paw_to_deeq(d_deeq)
!
! Put the nteraction on the smooth grid.
!
CALL interpolate (dvrs(:,1),dvrss,-1)
ENDIF
!
ENDIF
!
! Make sure the psic workspace is availible
!
ALLOCATE ( psic (dfftp%nnr) )
!
IF( gamma_only ) THEN
!
CALL lr_apply_liouvillian_gamma()
@ -124,44 +166,48 @@ SUBROUTINE lr_apply_liouvillian( evc1, evc1_new, sevc1_new, interaction )
CALL lr_apply_liouvillian_k()
!
ENDIF
!
DEALLOCATE ( psic )
!
IF ( interaction .and. (.not.ltammd) ) THEN
!
! Normal interaction
!
WRITE(stdout,'(5X,"lr_apply_liouvillian: applying interaction: normal")')
WRITE(stdout,'(5X,"lr_apply_liouvillian: &
&applying interaction: normal")')
!
! Here evc1_new contains the interaction
! Here we add the two terms:
! [H(k) - E(k)] * evc1(k) + P_c(k) [dV_{HXC} * revc0(k)]
!
!OBM, blas
!sevc1_new=sevc1_new+(1.0d0,0.0d0)*evc1_new
CALL zaxpy(size_evc,cmplx(1.0d0,0.0d0,kind=dp),evc1_new(:,:,:),1,sevc1_new(:,:,:),1)
CALL zaxpy(size_evc,CMPLX(1.0d0,0.0d0,kind=dp),&
&evc1_new(:,:,:), 1, sevc1_new(:,:,:), 1)
!
!
ELSEIF ( interaction .and. ltammd ) THEN
!
! Tamm-dancoff interaction
!
WRITE(stdout,'(5X,"lr_apply_liouvillian: applying interaction: tamm-dancoff")')
WRITE(stdout,'(5X,"lr_apply_liouvillian:&
& applying interaction: tamm-dancoff")')
!
! Here evc1_new contains the interaction
!
!OBM, blas
!sevc1_new=sevc1_new+(0.50d0,0.0d0)*evc1_new
CALL zaxpy(size_evc,cmplx(0.5d0,0.0d0,kind=dp),evc1_new(:,:,:),1,sevc1_new(:,:,:),1)
CALL zaxpy(size_evc,CMPLX(0.5d0,0.0d0,kind=dp),&
&evc1_new(:,:,:) , 1 , sevc1_new(:,:,:),1)
!
!
ELSE
!
! Non interacting
!
WRITE(stdout,'(5X,"lr_apply_liouvillian: not applying interaction")')
WRITE(stdout,'(5X,"lr_apply_liouvillian:&
& not applying interaction")')
!
ENDIF
!
IF (gstart == 2 .and. gamma_only ) sevc1_new(1,:,:)=cmplx(real(sevc1_new(1,:,:),dp),0.0d0,dp)
! (OBM: Why there is this check?)
IF (gstart == 2 .AND. gamma_only ) sevc1_new(1,:,:) = &
&CMPLX( REAL( sevc1_new(1,:,:), dp ), 0.0d0 ,dp )
!
IF(gstart==2 .and. gamma_only) THEN
!
DO ik=1,nks
@ -183,9 +229,19 @@ SUBROUTINE lr_apply_liouvillian( evc1, evc1_new, sevc1_new, interaction )
!
ENDIF
!
! Here we apply the S^{-1} operator.
! See equations after Eq.(47) of B. Walker et al., J. Chem. Phys.
! 127, 164106 (2007).
!
! S^{-1} ( H(k)*evc1(k) - E(k) * S evc1(k) )
! or
! S^{-1} ( H(k)*evc1(k) - E(k) * S evc1(k) + P_c(k) [dV_{HXC} *
! revc0(k)] )
!
DO ik=1,nks
!
CALL sm1_psi(.false.,ik,npwx,npw_k(ik),nbnd,sevc1_new(1,1,ik),evc1_new(1,1,ik))
CALL sm1_psi(.FALSE., ik, npwx, npw_k(ik), nbnd, &
&sevc1_new(1,1,ik), evc1_new(1,1,ik))
!
ENDDO
!
@ -194,10 +250,9 @@ SUBROUTINE lr_apply_liouvillian( evc1, evc1_new, sevc1_new, interaction )
DEALLOCATE(d_deeq)
DEALLOCATE(spsi1)
!
9000 FORMAT(/5x,'lr_apply_liouvillian: ibnd=',1X,i2,1X,'sevc1_new(G=0)[',i1,']',2(1x,e12.5)/)
!
CALL stop_clock('lr_apply')
IF (interaction) CALL stop_clock('lr_apply_int')
IF (interaction) CALL stop_clock('lr_apply_int')
IF (.not.interaction) CALL stop_clock('lr_apply_no')
!
RETURN
@ -206,7 +261,6 @@ CONTAINS
!
SUBROUTINE lr_apply_liouvillian_gamma()
!
!use becmod, only : bec_type,becp
USE lr_variables, ONLY : becp1
!
real(kind=dp), ALLOCATABLE :: becp2(:,:)
@ -271,19 +325,17 @@ CONTAINS
!
DO ibnd=1,nbnd,2
!
! Product with the potential vrs = (vltot+vr)
! revc0 is on smooth grid. psic is used upto smooth grid
! Product with the potential vrs = (vltot+vr)
! revc0 is on smooth grid. psic is used up to smooth grid
!
DO ir=1,dffts%nnr
!
psic(ir)=revc0(ir,ibnd,1)*cmplx(dvrss(ir),0.0d0,dp)
!
ENDDO
!
!print *,"1"
IF (real_space_debug > 7 .and. okvan .and. nkb > 0) THEN
!THE REAL SPACE PART (modified from s_psi)
!print *, "lr_apply_liouvillian:Experimental interaction part not using vkb"
!fac = sqrt(omega)
!
ijkb0 = 0
@ -297,7 +349,6 @@ CONTAINS
IF ( ityp(ia) == nt ) THEN
!
mbia = maxbox_beta(ia)
!print *, "mbia=",mbia
ALLOCATE( w1(nh(nt)), w2(nh(nt)) )
w1 = 0.D0
w2 = 0.D0
@ -308,7 +359,8 @@ CONTAINS
!
jkb = ijkb0 + jh
w1(ih) = w1(ih) + becp2(jkb, ibnd)
IF ( ibnd+1 <= nbnd ) w2(ih) = w2(ih) + becp2(jkb, ibnd+1)
IF ( ibnd+1 <= nbnd ) w2(ih) = w2(ih) + &
&becp2(jkb, ibnd+1)
!
ENDDO
!
@ -323,7 +375,10 @@ CONTAINS
DO ir = 1, mbia
!
iqs = jqs + ir
psic( box_beta(ir,ia) ) = psic( box_beta(ir,ia) ) + betasave(ia,ih,ir)*cmplx( w1(ih), w2(ih) )
psic( box_beta(ir,ia) ) = &
&psic( box_beta(ir,ia) ) + &
&betasave(ia,ih,ir)*&
&CMPLX( w1(ih), w2(ih) )
!
ENDDO
!
@ -340,12 +395,11 @@ CONTAINS
ENDDO
ENDIF
!print *,"2"
!
! Back to reciprocal space This part is equivalent to bfft_orbital_gamma
! Back to reciprocal space
!
CALL bfft_orbital_gamma (evc1_new(:,:,1), ibnd, nbnd,.false.)
!print *,"3"
!
ENDDO
!
!
@ -353,9 +407,9 @@ CONTAINS
!The non real_space part
CALL dgemm( 'N', 'N', 2*npw_k(1), nbnd, nkb, 1.d0, vkb, &
2*npwx, becp2, nkb, 1.d0, evc1_new, 2*npwx )
!print *, "lr_apply_liouvillian:interaction part using vkb"
!
ENDIF
!
CALL stop_clock('interaction')
!
ENDIF
@ -380,7 +434,8 @@ CONTAINS
!
DO ibnd=1,nbnd
!
CALL zaxpy(npw_k(1), cmplx(-et(ibnd,1),0.0d0,dp), spsi1(:,ibnd), 1, sevc1_new(:,ibnd,1), 1)
CALL zaxpy(npw_k(1), CMPLX(-et(ibnd,1),0.0d0,dp), &
&spsi1(:,ibnd), 1, sevc1_new(:,ibnd,1), 1)
!
ENDDO
!
@ -392,12 +447,11 @@ CONTAINS
!
SUBROUTINE lr_apply_liouvillian_k()
!
!use becmod, only : becp
USE lr_variables, ONLY : becp1_c
!
COMPLEX(kind=dp), ALLOCATABLE :: becp2(:,:)
!
IF( nkb > 0 .and. okvan ) THEN
IF( nkb > 0 .AND. okvan ) THEN
!
ALLOCATE(becp2(nkb,nbnd))
becp2(:,:)=(0.0d0,0.0d0)

489
TDDFPT/src/lr_calc_dens.f90
View File

@ -7,153 +7,158 @@ SUBROUTINE lr_calc_dens( evc1, response_calc )
!
! Modified by Osman Baris Malcioglu in 2009
!
! Input : evc1 (qdash etc) Output: rho_1 (=2*sum_v (revc0_v(r) . revc1_v(r,w) v:valance state index, r denotes a transformation t
! in case of uspps becsum is also calculated here
! in case of charge response calculation, the rho_tot is calculated here
! Input : evc1 (qdash etc)
! Output: rho_1 (=2*sum_v (revc0_v(r) . revc1_v(r,w) )
! where v:valance state index, r denotes a transformation t
!
USE ions_base, ONLY : ityp,nat,ntyp=>nsp
USE cell_base, ONLY : omega
USE ener, ONLY : ef
USE gvecs, ONLY : nls,nlsm,doublegrid
USE fft_base, ONLY : dffts, dfftp
USE fft_interfaces, ONLY : invfft
USE io_global, ONLY : stdout
USE kinds, ONLY : dp
USE klist, ONLY : nks,xk,wk
USE lr_variables, ONLY : evc0,revc0,rho_1,lr_verbosity, &
charge_response, itermax,&
cube_save, rho_1_tot,rho_1_tot_im, &
LR_iteration, LR_polarization, &
project,evc0_virt,F,nbnd_total,n_ipol, becp1_virt
USE lsda_mod, ONLY : current_spin, isk
USE wavefunctions_module, ONLY : psic
USE wvfct, ONLY : nbnd,et,wg,npwx,npw
USE control_flags, ONLY : gamma_only
USE uspp, ONLY : vkb,nkb,okvan,qq,becsum
USE uspp_param, ONLY : upf, nh
USE io_global, ONLY : ionode, stdout
USE io_files, ONLY : tmp_dir, prefix
USE mp, ONLY : mp_sum
USE mp_global, ONLY : inter_pool_comm, intra_pool_comm,nproc
USE realus, ONLY : igk_k,npw_k,addusdens_r
USE charg_resp, ONLY : w_T, lr_dump_rho_tot_cube,&
lr_dump_rho_tot_xyzd, &
lr_dump_rho_tot_xcrys,&
resonance_condition,epsil
USE noncollin_module, ONLY : nspin_mag
USE control_flags, ONLY : tqr
USE becmod, ONLY : becp
! In case of US PP, becsum is also calculated here
! In case of charge response calculation, the rho_tot is calculated here
!
USE ions_base, ONLY : ityp, nat, ntyp=>nsp
USE cell_base, ONLY : omega
USE ener, ONLY : ef
USE gvecs, ONLY : nls, nlsm, doublegrid
USE fft_base, ONLY : dffts, dfftp
USE fft_interfaces, ONLY : invfft
USE io_global, ONLY : stdout
USE kinds, ONLY : dp
USE klist, ONLY : nks,xk,wk
USE lr_variables, ONLY : evc0,revc0,rho_1,lr_verbosity,&
& charge_response, itermax,&
& cube_save, rho_1_tot&
&,rho_1_tot_im, LR_iteration,&
& LR_polarization, project,&
& evc0_virt, F,nbnd_total,&
& n_ipol, becp1_virt
USE lsda_mod, ONLY : current_spin, isk
USE wavefunctions_module, ONLY : psic
USE wvfct, ONLY : nbnd, et, wg, npwx, npw
USE control_flags, ONLY : gamma_only
USE uspp, ONLY : vkb, nkb, okvan, qq, becsum
USE uspp_param, ONLY : upf, nh
USE io_global, ONLY : ionode, stdout
USE io_files, ONLY : tmp_dir, prefix
USE mp, ONLY : mp_sum
USE mp_global, ONLY : inter_pool_comm, intra_pool_comm,&
& nproc
USE realus, ONLY : igk_k,npw_k, addusdens_r
USE charg_resp, ONLY : w_T, lr_dump_rho_tot_cube,&
& lr_dump_rho_tot_xyzd, &
& lr_dump_rho_tot_xcrys,&
& resonance_condition, epsil
USE noncollin_module, ONLY : nspin_mag
USE control_flags, ONLY : tqr
USE becmod, ONLY : becp
USE constants, ONLY : eps12
!
IMPLICIT NONE
!
CHARACTER(len=6), EXTERNAL :: int_to_char
CHARACTER(len=6), EXTERNAL :: int_to_char
!
COMPLEX(kind=dp), INTENT(in) :: evc1(npwx,nbnd,nks)
LOGICAL, INTENT(in) :: response_calc
LOGICAL, INTENT(in) :: response_calc
!
! functions
real(kind=dp) :: ddot
REAL(kind=dp) :: ddot
!
! local variables
INTEGER :: ir,ik,ibnd,jbnd,ig,ijkb0,np,na,ijh,ih,jh,ikb,jkb,ispin
INTEGER :: i, j, k, l
real(kind=dp) :: w1,w2,scal
real(kind=dp) :: rho_sum!,weight
real(kind=dp), ALLOCATABLE :: rho_sum_resp_x(:),rho_sum_resp_y(:),rho_sum_resp_z(:) ! These are temporary buffers for response cha
!complex(kind=dp), external :: ZDOTC
!complex(kind=dp), allocatable :: spsi(:,:)
INTEGER :: ir, ik, ibnd, jbnd, ig, ijkb0, np, na
INTEGER :: ijh,ih,jh,ikb,jkb ,ispin
INTEGER :: i, j, k, l
REAL(kind=dp) :: w1, w2, scal
REAL(kind=dp) :: rho_sum
!
! These are temporary buffers for the response
REAL(kind=dp), ALLOCATABLE :: rho_sum_resp_x(:), rho_sum_resp_y(:),&
& rho_sum_resp_z(:)
!
CHARACTER(len=256) :: tempfile, filename
!
!OBM DEBUG
COMPLEX(kind=dp),EXTERNAL :: lr_dot
!
IF (lr_verbosity > 5) THEN
WRITE(stdout,'("<lr_calc_dens>")')
WRITE(stdout,'("<lr_calc_dens>")')
ENDIF
!
CALL start_clock('lr_calc_dens')
!
!allocate(spsi(npwx,nbnd))
!spsi(:,:)=(0.0d0,0.0d0)
!
ALLOCATE( psic(dfftp%nnr) )
psic(:)=(0.0d0,0.0d0)
rho_1(:,:)=0.0d0
psic(:) = (0.0d0,0.0d0)
rho_1(:,:) = 0.0d0
!
IF(gamma_only) THEN
CALL lr_calc_dens_gamma()
ELSE
CALL lr_calc_dens_k()
ENDIF
!print *, "rho_1 after lr_calc_dens calculates",SUM(rho_1)
!print *, "norm of evc1 after lr_calc_dens calculates", lr_dot(evc1(1,1,1),evc1(1,1,1))
!
! ... If a double grid is used, interpolate onto the fine grid
! If a double grid is used, interpolate onto the fine grid
!
IF ( doublegrid ) CALL interpolate(rho_1,rho_1,1)
!
! ... Here we add the Ultrasoft contribution to the charge
! Here we add the Ultrasoft contribution to the charge density
! response.
!
!IF ( okvan ) CALL lr_addusdens(rho_1)
!print *, "rho_1 before addusdens",SUM(rho_1)
!call start_clock('lrcd_usdens') !TQR makes a huge gain here
IF(okvan) THEN
IF (tqr) THEN
CALL addusdens_r(rho_1,.false.)
ELSE
CALL addusdens(rho_1)
ENDIF
!
IF (tqr) THEN
CALL addusdens_r(rho_1,.FALSE.)
ELSE
CALL addusdens(rho_1)
ENDIF
ENDIF
DEALLOCATE ( psic )
!call stop_clock('lrcd_usdens')
!
!print *, "rho_1 after addusdens",SUM(rho_1)
! The psic workspace can present a memory bottleneck
DEALLOCATE ( psic )
!
#ifdef __MPI
!call poolreduce(nrxx,rho_1)
CALL mp_sum(rho_1, inter_pool_comm)
#endif
!
! check response charge density sums to 0
!call start_clock('lrcd_sp') !Minimal lag, no need to improve
IF (lr_verbosity > 0) THEN
DO ispin = 1, nspin_mag
rho_sum=0.0d0
rho_sum=sum(rho_1(:,ispin))
!
#ifdef __MPI
CALL mp_sum(rho_sum, intra_pool_comm )
#endif
!
rho_sum=rho_sum*omega/(dfftp%nr1*dfftp%nr2*dfftp%nr3)
!
IF(abs(rho_sum)>1.0d-12) THEN
IF (tqr) THEN
WRITE(stdout,'(5X, "lr_calc_dens: Charge drift due to real space implementation = " ,1X,e12.5)')&
rho_sum
!seems useless
!rho_sum=rho_sum/(1.0D0*nrxxs)
!rho_1(:,ispin)=rho_1(:,ispin)-rho_sum
ELSE
WRITE(stdout,'(5X,"lr_calc_dens: ****** response charge density does not sum to zero")')
!
WRITE(stdout,'(5X,"lr_calc_dens: ****** response charge density =",1X,e12.5)')&
rho_sum
!
WRITE(stdout,'(5X,"lr_calc_dens: ****** response charge density, US part =",1X,e12.5)')&
scal
! call errore(' lr_calc_dens ','Linear response charge density '// &
! & 'does not sum to zero',1)
ENDIF
ENDIF
ENDDO
!
ENDIF
IF (charge_response == 2 .and. LR_iteration /=0) THEN
IF (lr_verbosity > 0) THEN
!
DO ispin = 1, nspin_mag
!
rho_sum=0.0d0
rho_sum=SUM(rho_1(:,ispin))
!
#ifdef __MPI
CALL mp_sum(rho_sum, intra_pool_comm )
#endif
!
rho_sum = rho_sum * omega / (dfftp%nr1*dfftp%nr2*dfftp%nr3)
!
IF (ABS(rho_sum) > eps12) THEN
!
IF (tqr) THEN
!
WRITE(stdout,'(5X, "lr_calc_dens: Charge drift due to &
&real space implementation = " ,1X,e12.5)') rho_sum
!
ELSE
!
WRITE(stdout,'(5X,"lr_calc_dens: ****** response &
&charge density does not sum to zero")')
!
WRITE(stdout,'(5X,"lr_calc_dens: ****** response &
&charge density =",1X,e12.5)') rho_sum
!
WRITE(stdout,'(5X,"lr_calc_dens: ****** response &
&charge density, US part =",1X,e12.5)') scal
!
ENDIF
!
ENDIF
!
ENDDO
!
ENDIF
!
IF (charge_response == 2 .AND. LR_iteration /=0) THEN
!
ALLOCATE( rho_sum_resp_x( dfftp%nr1 ) )
ALLOCATE( rho_sum_resp_y( dfftp%nr2 ) )
@ -181,41 +186,45 @@ ENDIF
CALL mp_sum(rho_sum_resp_z, intra_pool_comm)
IF (ionode) THEN
#endif
WRITE(stdout,'(5X,"Dumping plane sums of densities for iteration ",I4)') LR_iteration
!
filename = trim(prefix) // ".density_x"
tempfile = trim(tmp_dir) // trim(filename)
!
OPEN (158, file = tempfile, form = 'formatted', status = 'unknown', position = 'append')
!
DO i=1,dfftp%nr1
WRITE(158,*) rho_sum_resp_x(i)
ENDDO
!
CLOSE(158)
!
filename = trim(prefix) // ".density_y"
tempfile = trim(tmp_dir) // trim(filename)
!
OPEN (158, file = tempfile, form = 'formatted', status = 'unknown', position = 'append')
!
DO i=1,dfftp%nr2
WRITE(158,*) rho_sum_resp_y(i)
ENDDO
!
CLOSE(158)
!
filename = trim(prefix) // ".density_z"
tempfile = trim(tmp_dir) // trim(filename)
!
OPEN (158, file = tempfile, form = 'formatted', status = 'unknown', position = 'append')
!
DO i=1,dfftp%nr3
WRITE(158,*) rho_sum_resp_z(i)
ENDDO
!
CLOSE(158)
!
WRITE(stdout,'(5X,"Dumping plane sums of densities for &
&iteration ",I4)') LR_iteration
!
filename = TRIM(prefix) // ".density_x"
tempfile = TRIM(tmp_dir) // TRIM(filename)
!
OPEN (158, file = tempfile, form = 'formatted', status = &
&'unknown', position = 'append')
!
DO i=1,dfftp%nr1
WRITE(158,*) rho_sum_resp_x(i)
ENDDO
!
CLOSE(158)
!
filename = TRIM(prefix) // ".density_y"
tempfile = TRIM(tmp_dir) // TRIM(filename)
!
OPEN (158, file = tempfile, form = 'formatted', status = &
&'unknown', position = 'append')
!
DO i=1,dfftp%nr2
WRITE(158,*) rho_sum_resp_y(i)
ENDDO
!
CLOSE(158)
!
filename = TRIM(prefix) // ".density_z"
tempfile = TRIM(tmp_dir) // TRIM(filename)
!
OPEN (158, file = tempfile, form = 'formatted', status = &
&'unknown', position = 'append')
!
DO i=1,dfftp%nr3
WRITE(158,*) rho_sum_resp_z(i)
ENDDO
!
CLOSE(158)
!
#ifdef __MPI
ENDIF
#endif
@ -225,67 +234,78 @@ ENDIF
DEALLOCATE( rho_sum_resp_z )
!
ENDIF
IF (charge_response == 1 .and. response_calc) THEN
IF (LR_iteration < itermax) WRITE(stdout,'(5x,"Calculating total response charge density")')
IF (charge_response == 1 .AND. response_calc) THEN
IF (LR_iteration < itermax) WRITE(stdout,'(5x,"Calculating total &
&response charge density")')
! the charge response, it is actually equivalent to an element of
! V^T . phi_v where V^T is the is the transpose of the Krylov subspace generated
! by the Lanczos algorithm. The total charge density can be written
! as
! V^T . phi_v where V^T is the is the transpose of the Krylov
! subspace generated by the Lanczos algorithm. The total charge
! density can be written as,
!
! \sum_(lanczos iterations) (V^T.phi_v) . w_T
! Where w_T is the corresponding eigenvector from the solution of
!
! Where w_T is the corresponding eigenvector from the solution of,
!
! (w-L)e_1 = w_T
!
! notice that rho_1 is already reduced across pools above, so no parallelization is necessary
! notice that rho_1 is already reduced across pools above, so no
! parallelization is necessary
!
! the lr_calc_dens corresponds to q of x only in even iterations
!
!print *,"1"
!print *,"weight",(-1.0d0*AIMAG(w_T(LR_iteration)))
IF (resonance_condition) THEN
!singular matrix, the broadening term dominates, phi' has strong imaginary component
! Using BLAS here probably results in cmplx(rho_1(:,1),0.0d0,dp) being copied into a new array
! due to the call being to an F77 funtion.
!CALL zaxpy(dfftp%nnr, w_T(LR_iteration),cmplx(rho_1(:,1),0.0d0,dp),1,rho_1_tot_im(:,1),1) !spin not implemented
!
! Singular matrix, the broadening term dominates, phi' has
! strong imaginary component
!
! Using BLAS here would result in cmplx(rho_1(:,1),0.0d0 ,dp)
! being copied into a NEW array due to the call being to an
! F77 funtion.
!
rho_1_tot_im(1:dfftp%nnr,:) = rho_1_tot_im(1:dfftp%nnr,:) &
& + w_T(LR_iteration) * cmplx(rho_1(1:dfftp%nnr,:),0.0d0,dp)
!
ELSE
!not at resonance, the imaginary part is neglected ,these are the non-absorbing oscillations
!CALL daxpy(dfftp%nnr, dble(w_T(LR_iteration)),rho_1(:,1),1,rho_1_tot(:,1),1) !spin not implemented
!
! Not at resonance.
! The imaginary part is neglected, these are the non-absorbing
! oscillations
!
rho_1_tot(1:dfftp%nnr,:) = rho_1_tot(1:dfftp%nnr,:) &
& + dble( w_T(LR_iteration) ) * rho_1(1:dfftp%nnr,:)
ENDIF
!
ENDIF
!
!deallocate(spsi)
!
CALL stop_clock('lr_calc_dens')
!
!call stop_clock('lrcd_sp')
RETURN
!
!
ENDIF
!
!
CALL stop_clock('lr_calc_dens')
RETURN
!
CONTAINS
!
SUBROUTINE lr_calc_dens_gamma
!
USE becmod, ONLY : bec_type, becp, calbec
USE lr_variables, ONLY : becp1 !,real_space
!use real_beta, only : ccalbecr_gamma, fft_orbital_gamma
USE io_global, ONLY : stdout
USE realus, ONLY : real_space, fft_orbital_gamma, initialisation_level, &
bfft_orbital_gamma, calbec_rs_gamma, add_vuspsir_gamma, v_loc_psir,&
real_space_debug
! Gamma_only case.
!
USE becmod, ONLY : bec_type, becp, calbec
USE lr_variables, ONLY : becp1
USE io_global, ONLY : stdout
USE realus, ONLY : real_space, fft_orbital_gamma,&
& initialisation_level,&
& bfft_orbital_gamma,&
& calbec_rs_gamma,&
& add_vuspsir_gamma, v_loc_psir,&
& real_space_debug
!
DO ibnd=1,nbnd,2
!
! FFT: evc1 -> psic
!
CALL fft_orbital_gamma(evc1(:,:,1),ibnd,nbnd)
!
! Set weights of the two real bands now in psic
!
w1=wg(ibnd,1)/omega
!
IF(ibnd<nbnd) THEN
@ -293,53 +313,64 @@ CONTAINS
ELSE
w2=w1
ENDIF
!call start_clock('lrcd-lp1')
! OBM:
!
! (n'(r,w)=2*sum_v (psi_v(r) . q_v(r,w))
! where psi are the ground state valance orbitals
! and q_v are the standart batch representation (rotated)
! and q_v are the standard batch representation (rotated)
! response orbitals
! Here, since the ith iteration is the best approximate we have
! for the most dominant eigenvalues/vectors, an estimate for the response
! charge density can be calculated. This is in no way the final
! response charge density.
! the loop is over real space points.
! Here, since the ith iteration is the best approximation we
! have for the most dominant eigenvalues/vectors, an estimate
! for the response charge density can be calculated. This is
! in no way the final response charge density.
!
! The loop is over real space points.
!
DO ir=1,dffts%nnr
rho_1(ir,1)=rho_1(ir,1) &
+2.0d0*(w1*real(revc0(ir,ibnd,1),dp)*real(psic(ir),dp)&
+w2*aimag(revc0(ir,ibnd,1))*aimag(psic(ir)))
ENDDO
!
!call stop_clock('lrcd-lp1')
! OBM - psic now contains the response functions at
! real space, eagerly putting all the real space stuff at this point.
! notice that betapointlist() is called in lr_readin at the very start
IF ( real_space_debug > 6 .and. okvan) THEN
! The rbecp term
CALL calbec_rs_gamma(ibnd,nbnd,becp%r)
! OBM - psic now contains the response functions in real space.
! Eagerly putting all the real space stuff at this point.
!
! Notice that betapointlist() is called in lr_readin at the
! very start
!
IF ( real_space_debug > 6 .AND. okvan) THEN
! The rbecp term
CALL calbec_rs_gamma(ibnd,nbnd,becp%r)
!
ENDIF
!
! End of real space stuff
!
ENDDO
!
! ... If we have a US pseudopotential we compute here the becsum term
! This corresponds to the right hand side of the formula (36) in Ultrasoft paper
! be careful about calling lr_calc_dens, as it modifies this globally
!call start_clock('lrcd-us')
! If we have a US pseudopotential we compute here the becsum
! term.
! This corresponds to the right hand side of the formula (36) in
! the ultrasoft paper.
!
! Be careful about calling lr_calc_dens, as it modifies this
! globally.
!
IF ( okvan ) THEN
!
scal = 0.0d0
becsum(:,:,:) = 0.0d0
!
IF ( real_space_debug <= 6) THEN !in real space, the value is calculated above
!call pw_gemm('Y',nkb,nbnd,npw_k(1),vkb,npwx,evc1,npwx,rbecp,nkb)
IF ( real_space_debug <= 6) THEN
! In real space, the value is calculated above
CALL calbec(npw_k(1), vkb, evc1(:,:,1), becp)
!
ENDIF
!
CALL start_clock( 'becsum' )
!
DO ibnd = 1, nbnd
scal = 0.0d0
!
scal = 0.0d0
w1 = wg(ibnd,1)
ijkb0 = 0
!
@ -359,8 +390,11 @@ CONTAINS
!
becsum(ijh,na,current_spin) = &
becsum(ijh,na,current_spin) + &
2.d0 * w1 * becp%r(ikb,ibnd) * becp1(ikb,ibnd)
scal = scal + qq(ih,ih,np) *1.d0 * becp%r(ikb,ibnd) * becp1(ikb,ibnd)
2.d0 * w1 * becp%r(ikb,ibnd) *&
& becp1(ikb,ibnd)
!
scal = scal + qq(ih,ih,np) *1.d0 *&
& becp%r(ikb,ibnd) * becp1(ikb,ibnd)
!
ijh = ijh + 1
!
@ -370,10 +404,14 @@ CONTAINS
!
becsum(ijh,na,current_spin) = &
becsum(ijh,na,current_spin) + &
w1 * 2.D0 * (becp1(ikb,ibnd) * becp%r(jkb,ibnd) + &
w1 * 2.D0 * (becp1(ikb,ibnd) * &
&becp%r(jkb,ibnd) + &
becp1(jkb,ibnd) * becp%r(ikb,ibnd))
scal = scal + qq(ih,jh,np) *1.d0 * (becp%r(ikb,ibnd) * becp1(jkb,ibnd)+&
becp%r(jkb,ibnd) * becp1(ikb,ibnd))
!
scal = scal + qq(ih,jh,np) * 1.d0 *&
& (becp%r(ikb,ibnd) * &
&becp1(jkb, ibnd) + &
&becp%r(jkb,ibnd) * becp1(ikb,ibnd))
!
ijh = ijh + 1
!
@ -399,20 +437,16 @@ CONTAINS
!
ENDDO
!
! OBM debug
!write(stdout,'(5X,"lr_calc_dens: ibnd,scal=",1X,i3,1X,e12.5)')&
! ibnd,scal
ENDDO
!
CALL stop_clock( 'becsum' )
!
ENDIF
!call stop_clock('lrcd-us')
!
RETURN
!
END SUBROUTINE lr_calc_dens_gamma
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
SUBROUTINE lr_calc_dens_k
!
USE becmod, ONLY : bec_type, becp, calbec
@ -478,8 +512,13 @@ CONTAINS
!
becsum(ijh,na,current_spin) = &
becsum(ijh,na,current_spin) + &
2.d0 * w1 * real(conjg(becp%k(ikb,ibnd)) * becp1_c(ikb,ibnd,ik))
scal = scal + qq(ih,ih,np) *1.d0 * real(conjg(becp%k(ikb,ibnd)) * becp1_c(ikb,ibnd,ik))
2.d0 * w1 * &
&DBLE(CONJG(becp%k(ikb,ibnd)) *&
& becp1_c(ikb,ibnd,ik))
!
scal = scal + qq(ih,ih,np) * 1.d0 *&
& DBLE(CONJG(becp%k(ikb,ibnd)) *&
& becp1_c(ikb,ibnd,ik))
!
ijh = ijh + 1
!
@ -489,10 +528,17 @@ CONTAINS
!
becsum(ijh,na,current_spin) = &
becsum(ijh,na,current_spin) + &
w1 * 2.d0 * real(conjg(becp1_c(ikb,ibnd,ik)) * becp%k(jkb,ibnd) + &
becp1_c(jkb,ibnd,ik) * conjg(becp%k(ikb,ibnd)))
scal = scal + qq(ih,jh,np) *1.d0 * real(conjg(becp%k(ikb,ibnd)) * becp1_c(jkb,ibnd,ik)+&
becp%k(jkb,ibnd) * conjg(becp1_c(ikb,ibnd,ik)))
w1 * 2.d0 * DBLE(&
& CONJG(becp1_c(ikb,ibnd,ik)) *&
& becp%k(jkb,ibnd) + &
becp1_c(jkb,ibnd,ik) *&
& CONJG(becp%k(ikb,ibnd)))
!
scal = scal + qq(ih,jh,np) * 1.d0 * &
& DBLE(CONJG(becp%k(ikb,ibnd)) *&
& becp1_c(jkb,ibnd,ik)+&
& becp%k(jkb,ibnd) * &
& CONJG(becp1_c(ikb,ibnd,ik)))
!
ijh = ijh + 1
!
@ -518,9 +564,8 @@ CONTAINS
!
ENDDO
!
! write(stdout,'(5X,"lr_calc_dens: ibnd,scal=",1X,i3,1X,e12.5)')&
! ibnd,scal
ENDDO
!
CALL stop_clock( 'becsum' )
!
ENDDO
@ -530,7 +575,7 @@ CONTAINS
RETURN
!
END SUBROUTINE lr_calc_dens_k
!-------------------------------------------------------------------------------
!-----------------------------------------------------------------------
!--------------------------------------------------------------------
END SUBROUTINE lr_calc_dens
!-----------------------------------------------------------------------
!--------------------------------------------------------------------

93
TDDFPT/src/lr_normalise.f90
View File

@ -6,18 +6,18 @@
SUBROUTINE lr_normalise(evc1,norm)
!
!
USE kinds, ONLY : dp
USE gvect, ONLY : gstart
USE cell_base, ONLY : omega
USE io_global, ONLY : stdout
USE kinds, ONLY : dp
USE klist, ONLY : nks,xk
USE lsda_mod, ONLY : nspin
USE lr_variables, ONLY : lanc_norm
USE realus, ONLY : igk_k,npw_k
USE uspp, ONLY : vkb,nkb,okvan
USE wvfct, ONLY : nbnd,npwx,npw,wg
USE realus, ONLY : igk_k, npw_k
USE uspp, ONLY : vkb, nkb, okvan
USE wvfct, ONLY : nbnd, npwx, npw, wg
USE control_flags, ONLY : gamma_only
USE lr_variables, ONLY : lr_verbosity
USE lr_variables, ONLY : lr_verbosity
!
IMPLICIT NONE
!
@ -46,74 +46,72 @@ SUBROUTINE lr_normalise(evc1,norm)
!
CONTAINS
!
!--------------------------------------------------------------------
SUBROUTINE lr_normalise_gamma()
!
USE becmod, ONLY : bec_type, becp,calbec
!use lr_variables, only : real_space
!use real_beta, only : ccalbecr_gamma,s_psir,fft_orbital_gamma,bfft_orbital_gamma
USE realus, ONLY : real_space, fft_orbital_gamma, initialisation_level, &
bfft_orbital_gamma, calbec_rs_gamma, add_vuspsir_gamma, &
v_loc_psir, s_psir_gamma, real_space_debug
USE becmod, ONLY : bec_type, becp,calbec
USE realus, ONLY : real_space, fft_orbital_gamma,&
& initialisation_level,&
& bfft_orbital_gamma, calbec_rs_gamma,&
& add_vuspsir_gamma, v_loc_psir,&
& s_psir_gamma, real_space_debug
!
!
!
IMPLICIT NONE
!
real(kind=dp) :: prod
REAL(kind=dp) :: prod
COMPLEX(kind=dp), EXTERNAL :: lr_dot
INTEGER :: ibnd,ig
!
prod=0.0d0
!
IF ( nkb > 0 ) THEN
!
IF (real_space_debug>6) THEN
! real space & nkb > 0
!
DO ibnd=1,nbnd,2
CALL fft_orbital_gamma(evc1(:,:,1),ibnd,nbnd)
CALL calbec_rs_gamma(ibnd,nbnd,becp%r)
CALL s_psir_gamma(ibnd,nbnd)
CALL bfft_orbital_gamma(spsi(:,:,1),ibnd,nbnd)
ENDDO
!
!
ELSE
!
!the non real_space & nkb > 0 case
!
CALL calbec(npw_k(1),vkb,evc1(:,:,1),becp)
!call pw_gemm('Y',nkb,nbnd,npw_k(1),vkb,npwx,evc1(1,1,1),npwx,rbecp,nkb)
!
CALL s_psi(npwx,npw_k(1),nbnd,evc1(1,1,1),spsi)
!
!
ENDIF
IF (real_space_debug>6) THEN
! real space & nkb > 0
!
DO ibnd=1,nbnd,2
CALL fft_orbital_gamma(evc1(:,:,1),ibnd,nbnd)
CALL calbec_rs_gamma(ibnd,nbnd,becp%r)
CALL s_psir_gamma(ibnd,nbnd)
CALL bfft_orbital_gamma(spsi(:,:,1),ibnd,nbnd)
ENDDO
!
!
ELSE
!
!the non real_space & nkb > 0 case
!
CALL calbec(npw_k(1),vkb,evc1(:,:,1),becp)
!
CALL s_psi(npwx,npw_k(1),nbnd,evc1(1,1,1),spsi)
!
!
ENDIF
ELSE
! The nkb == 0 part
! JUST array copying
! The nkb == 0 part
! JUST array copying
CALL s_psi(npwx,npw_k(1),nbnd,evc1(1,1,1),spsi)
!
!
!
!
!
!
ENDIF
!The below two lines are the replicated part in real space implementation
!call calbec(npw_k(1),vkb,evc1(:,:,1),rbecp)
!call s_psi(npwx,npw_k(1),nbnd,evc1(1,1,1),spsi)
!
prod=dble( lr_dot( evc1(1,1,1),spsi(1,1,1) ) )
prod=1.0d0/sqrt(abs(prod))
!
evc1(:,:,1)=cmplx(prod,0.0d0,dp)*evc1(:,:,1)
!
WRITE(stdout,'(5X,"Norm of initial Lanczos vectors=",1x,f21.15)') 1.0d0/prod
WRITE(stdout,'(5X,"Norm of initial Lanczos vectors=",1x,f21.15)')&
& 1.0d0/prod
!
lanc_norm=1.d0/prod**2/omega
norm=1.0d0/prod
!
RETURN
END SUBROUTINE lr_normalise_gamma
!
!--------------------------------------------------------------------
SUBROUTINE lr_normalise_k()
!
USE becmod, ONLY : becp,calbec
@ -129,7 +127,6 @@ CONTAINS
!
CALL init_us_2(npw_k(ik),igk_k(1,ik),xk(1,ik),vkb)
!
!call ccalbec(nkb,npwx,npw_k(ik),nbnd,becp,vkb,evc1(1,1,ik))
CALL calbec(npw_k(ik),vkb,evc1(:,:,ik),becp)
!
ENDIF
@ -143,7 +140,9 @@ CONTAINS
!
evc1(:,:,:)=cmplx(prod,0.0d0,dp)*evc1(:,:,:)
!
WRITE(stdout,'(5X,"Norm of initial Lanczos vectors=",1x,f21.15)') 1.0d0/prod
WRITE(stdout,'(5X,"Norm of initial Lanczos vectors=",1x,f21.15)')&
& 1.0d0/prod
!
lanc_norm=1.d0/prod**2/omega
norm=1.0d0/prod
!

989
TDDFPT/src/lr_read_wf.f90
View File

File diff suppressed because it is too large Load Diff