Commit ce47a54e authored by Marco Govoni's avatar Marco Govoni
Browse files

Added apply_operator

parent c3ad10a7
...@@ -48,7 +48,6 @@ MODULE dfpt_module ...@@ -48,7 +48,6 @@ MODULE dfpt_module
USE io_push, ONLY : io_push_title USE io_push, ONLY : io_push_title
USE mp_world, ONLY : mpime,world_comm USE mp_world, ONLY : mpime,world_comm
USE types_bz_grid, ONLY : k_grid, q_grid, compute_phase USE types_bz_grid, ONLY : k_grid, q_grid, compute_phase
USE types_coulomb, ONLY : pot3D
! !
IMPLICIT NONE IMPLICIT NONE
! !
...@@ -68,7 +67,7 @@ MODULE dfpt_module ...@@ -68,7 +67,7 @@ MODULE dfpt_module
INTEGER :: npwkq INTEGER :: npwkq
! !
REAL(DP) :: g0(3) REAL(DP) :: g0(3)
REAL(DP) :: anorm, prod REAL(DP) :: anorm
REAL(DP), ALLOCATABLE :: eprec(:) REAL(DP), ALLOCATABLE :: eprec(:)
! !
COMPLEX(DP), ALLOCATABLE :: dvpsi(:,:),dpsi(:,:) COMPLEX(DP), ALLOCATABLE :: dvpsi(:,:),dpsi(:,:)
...@@ -171,14 +170,11 @@ MODULE dfpt_module ...@@ -171,14 +170,11 @@ MODULE dfpt_module
! !
DO ipert = 1, m DO ipert = 1, m
! !
ALLOCATE( aux_g(npwqx) ) ALLOCATE( aux_g(npwqx) ); aux_g = 0._DP
ALLOCATE( aux_r(dffts%nnr) ) ALLOCATE( aux_r(dffts%nnr) ); aux_r = 0._DP
!
aux_g=0._DP
aux_r=0._DP
! !
DO ig = 1, npwq ! perturbation acts only on body DO CONCURRENT (ig = 1:npwq)
aux_g(ig) = dvg(ig,ipert) * pot3D%sqvc(ig) aux_g(ig) = dvg(ig,ipert)
ENDDO ENDDO
! !
! ... inverse Fourier transform of the perturbation: (q+)G ---> R ! ... inverse Fourier transform of the perturbation: (q+)G ---> R
...@@ -200,7 +196,7 @@ MODULE dfpt_module ...@@ -200,7 +196,7 @@ MODULE dfpt_module
DO ibnd=1,nbndval-MOD(nbndval,2),2 DO ibnd=1,nbndval-MOD(nbndval,2),2
! !
CALL double_invfft_gamma(dffts,npw,npwx,evc(1,ibnd),evc(1,ibnd+1),psic,'Wave') CALL double_invfft_gamma(dffts,npw,npwx,evc(1,ibnd),evc(1,ibnd+1),psic,'Wave')
DO ir=1,dffts%nnr DO CONCURRENT (ir=1:dffts%nnr)
psic(ir) = psic(ir) * REAL(aux_r(ir),KIND=DP) psic(ir) = psic(ir) * REAL(aux_r(ir),KIND=DP)
ENDDO ENDDO
CALL double_fwfft_gamma(dffts,npw,npwx,psic,dvpsi(1,ibnd),dvpsi(1,ibnd+1),'Wave') CALL double_fwfft_gamma(dffts,npw,npwx,psic,dvpsi(1,ibnd),dvpsi(1,ibnd+1),'Wave')
...@@ -212,7 +208,7 @@ MODULE dfpt_module ...@@ -212,7 +208,7 @@ MODULE dfpt_module
ibnd=nbndval ibnd=nbndval
! !
CALL single_invfft_gamma(dffts,npw,npwx,evc(1,ibnd),psic,'Wave') CALL single_invfft_gamma(dffts,npw,npwx,evc(1,ibnd),psic,'Wave')
DO ir=1,dffts%nnr DO CONCURRENT (ir=1:dffts%nnr)
psic(ir) = CMPLX( REAL(psic(ir),KIND=DP) * REAL(aux_r(ir),KIND=DP), 0._DP, KIND=DP) psic(ir) = CMPLX( REAL(psic(ir),KIND=DP) * REAL(aux_r(ir),KIND=DP), 0._DP, KIND=DP)
ENDDO ENDDO
CALL single_fwfft_gamma(dffts,npw,npwx,psic,dvpsi(1,ibnd),'Wave') CALL single_fwfft_gamma(dffts,npw,npwx,psic,dvpsi(1,ibnd),'Wave')
...@@ -230,7 +226,7 @@ MODULE dfpt_module ...@@ -230,7 +226,7 @@ MODULE dfpt_module
! ... construct right-hand-side term of Sternheimer equation: ! ... construct right-hand-side term of Sternheimer equation:
! ... product of wavefunction at [k-q], phase and perturbation in real space ! ... product of wavefunction at [k-q], phase and perturbation in real space
! !
DO ir = 1, dffts%nnr DO CONCURRENT (ir = 1:dffts%nnr)
psic(ir) = psic(ir) * phase(ir) * aux_r(ir) psic(ir) = psic(ir) * phase(ir) * aux_r(ir)
ENDDO ENDDO
! !
...@@ -248,7 +244,7 @@ MODULE dfpt_module ...@@ -248,7 +244,7 @@ MODULE dfpt_module
! !
CALL single_invfft_k(dffts,npwkq,npwx,evckmq(npwx+1,ibnd),psic,'Wave',igk_k(1,ikqs)) CALL single_invfft_k(dffts,npwkq,npwx,evckmq(npwx+1,ibnd),psic,'Wave',igk_k(1,ikqs))
! !
DO ir = 1, dffts%nnr DO CONCURRENT (ir = 1:dffts%nnr)
psic(ir) = psic(ir) * phase(ir) * aux_r(ir) psic(ir) = psic(ir) * phase(ir) * aux_r(ir)
ENDDO ENDDO
! !
...@@ -292,9 +288,8 @@ MODULE dfpt_module ...@@ -292,9 +288,8 @@ MODULE dfpt_module
DO ibnd=1,nbndval DO ibnd=1,nbndval
! !
CALL double_invfft_gamma(dffts,npw,npwx,evc(1,ibnd),dpsi(1,ibnd),psic,'Wave') CALL double_invfft_gamma(dffts,npw,npwx,evc(1,ibnd),dpsi(1,ibnd),psic,'Wave')
DO ir=1,dffts%nnr DO CONCURRENT (ir=1:dffts%nnr)
prod = REAL( psic(ir),KIND=DP) * DIMAG( psic(ir)) aux_r(ir) = aux_r(ir) + CMPLX( REAL( psic(ir),KIND=DP) * DIMAG( psic(ir)) , 0.0_DP, KIND=DP)
aux_r(ir) = aux_r(ir) + CMPLX( prod, 0.0_DP, KIND=DP)
ENDDO ENDDO
! !
ENDDO ENDDO
...@@ -313,7 +308,7 @@ MODULE dfpt_module ...@@ -313,7 +308,7 @@ MODULE dfpt_module
! !
CALL single_invfft_k(dffts,npw,npwx,dpsi(1,ibnd),dpsic,'Wave',igk_k(1,iks)) CALL single_invfft_k(dffts,npw,npwx,dpsi(1,ibnd),dpsic,'Wave',igk_k(1,iks))
! !
DO ir = 1, dffts%nnr DO CONCURRENT (ir = 1: dffts%nnr)
aux_r(ir) = aux_r(ir) + CONJG( psic(ir) * phase(ir) ) * dpsic(ir) aux_r(ir) = aux_r(ir) + CONJG( psic(ir) * phase(ir) ) * dpsic(ir)
ENDDO ENDDO
! !
...@@ -326,7 +321,7 @@ MODULE dfpt_module ...@@ -326,7 +321,7 @@ MODULE dfpt_module
! !
CALL single_invfft_k(dffts,npw,npwx,dpsi(npwx+1,ibnd),dpsic,'Wave',igk_k(1,iks)) CALL single_invfft_k(dffts,npw,npwx,dpsi(npwx+1,ibnd),dpsic,'Wave',igk_k(1,iks))
! !
DO ir = 1, dffts%nnr DO CONCURRENT (ir = 1: dffts%nnr)
aux_r(ir) = aux_r(ir) + CONJG( psic(ir) * phase(ir) ) * dpsic(ir) aux_r(ir) = aux_r(ir) + CONJG( psic(ir) * phase(ir) ) * dpsic(ir)
ENDDO ENDDO
! !
...@@ -347,8 +342,8 @@ MODULE dfpt_module ...@@ -347,8 +342,8 @@ MODULE dfpt_module
CALL single_fwfft_k(dffts,npwq,npwqx,aux_r,aux_g,'Wave',igq_q(1,iq)) CALL single_fwfft_k(dffts,npwq,npwqx,aux_r,aux_g,'Wave',igq_q(1,iq))
ENDIF ENDIF
! !
DO ig = 1, npwq ! pert acts only on body DO CONCURRENT( ig = 1: npwq) ! pert acts only on body
dng(ig,ipert) = dng(ig,ipert) + 2._DP * k_grid%weight(iks) * aux_g(ig) * pot3D%sqvc(ig) / omega dng(ig,ipert) = dng(ig,ipert) + 2._DP * k_grid%weight(iks) * aux_g(ig) / omega
ENDDO ENDDO
! !
DEALLOCATE( aux_g ) DEALLOCATE( aux_g )
......
...@@ -26,6 +26,7 @@ wstat_setup.o \ ...@@ -26,6 +26,7 @@ wstat_setup.o \
wstat_restart.o \ wstat_restart.o \
wstat_memory_report.o \ wstat_memory_report.o \
wstat_tools.o \ wstat_tools.o \
apply_operator.o \
davidson_diago.o \ davidson_diago.o \
wstat.o wstat.o
......
!
! Copyright (C) 2015-2017 M. Govoni
! 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 .
!
! This file is part of WEST.
!
! Contributors to this file:
! Marco Govoni
!
!-------------------------------------------------------------------------
SUBROUTINE apply_operator (m,dvg,dng,tr2,iq)
!-----------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE mp, ONLY : mp_barrier
USE westcom, ONLY : npwqx,npwq
USE mp_world, ONLY : world_comm
USE types_coulomb, ONLY : pot3D
USE dfpt_module, ONLY : dfpt
!
IMPLICIT NONE
!
! I/O
!
INTEGER, INTENT(IN), OPTIONAL :: iq
REAL(DP),INTENT(IN), OPTIONAL :: tr2
INTEGER, INTENT(IN) :: m
COMPLEX(DP), INTENT(IN) :: dvg(npwqx,m)
COMPLEX(DP), INTENT(OUT) :: dng(npwqx,m)
!
! Workspace
!
INTEGER :: ipert, ig
INTEGER :: iq_
!
REAL(DP) :: tr2_
!
COMPLEX(DP), ALLOCATABLE ::aux_g(:,:)
!
LOGICAL :: l_outsource
!
CALL mp_barrier( world_comm )
!
l_outsource = .FALSE.
!
IF(PRESENT(iq)) THEN
iq_ = iq
ELSE
iq_ = 1
ENDIF
IF(PRESENT(tr2)) THEN
tr2_ = tr2
ELSE
tr2_ = 1.d-8
ENDIF
!
dng=0._DP
!
ALLOCATE( aux_g(npwqx,m) ); aux_g=0._DP
!
DO CONCURRENT (ipert = 1:m, ig = 1:npwq)
aux_g(ig,ipert) = dvg(ig,ipert) * pot3D%sqvc(ig) ! perturbation acts only on body
ENDDO
!
IF( l_outsource ) THEN
CALL calc_outsourced(m,aux_g,dng)
ELSE
CALL dfpt(m,aux_g,dng,tr2_,iq_)
ENDIF
!
DO CONCURRENT (ipert = 1:m, ig = 1:npwq)
dng(ig,ipert) = dng(ig,ipert) * pot3D%sqvc(ig) ! perturbation acts only on body
ENDDO
!
CALL mp_barrier( world_comm )
!
END SUBROUTINE
!
!
!-------------------------------------------------------------------------
SUBROUTINE calc_outsourced (m,dvg,dng)
!-----------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE westcom, ONLY : npwqx
!
IMPLICIT NONE
!
! I/O
!
INTEGER, INTENT(IN) :: m
COMPLEX(DP), INTENT(IN) :: dvg(npwqx,m)
COMPLEX(DP), INTENT(OUT) :: dng(npwqx,m)
!
dng = dvg ! placeholder
!
END SUBROUTINE
...@@ -53,7 +53,7 @@ SUBROUTINE davidson_diago_gamma ( ) ...@@ -53,7 +53,7 @@ SUBROUTINE davidson_diago_gamma ( )
USE wstat_tools, ONLY : diagox,serial_diagox,build_hr,symm_hr_distr,redistribute_vr_distr,& USE wstat_tools, ONLY : diagox,serial_diagox,build_hr,symm_hr_distr,redistribute_vr_distr,&
& update_with_vr_distr,refresh_with_vr_distr & update_with_vr_distr,refresh_with_vr_distr
USE types_coulomb, ONLY : pot3D USE types_coulomb, ONLY : pot3D
USE dfpt_module USE dfpt_module, ONLY : dfpt
! !
IMPLICIT NONE IMPLICIT NONE
! !
...@@ -473,7 +473,7 @@ SUBROUTINE davidson_diago_k ( ) ...@@ -473,7 +473,7 @@ SUBROUTINE davidson_diago_k ( )
& update_with_vr_distr,refresh_with_vr_distr & update_with_vr_distr,refresh_with_vr_distr
USE types_bz_grid, ONLY : q_grid USE types_bz_grid, ONLY : q_grid
USE types_coulomb, ONLY : pot3D USE types_coulomb, ONLY : pot3D
USE dfpt_module USE dfpt_module, ONLY : dfpt
! !
IMPLICIT NONE IMPLICIT NONE
! !
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment