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

Bugfixes for AIX machines and intel compilers 

git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@6593 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
obm 2010-04-03 01:27:40 +00:00
parent 0a0321b82a
commit 0d4d58b532
8 changed files with 399 additions and 146 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
TDDFPT/Doc/TODO Normal file
View File

@ -0,0 +1,9 @@
1) Possible optimizations
1.1 revc0 each other band is not used in gamma point case
2) Feature enhancements
2.1 Exact exchange
2.2 PAW
2.3 Spin
3) Known bugs
3.1 Precision is not achieved in xl compiler

3
TDDFPT/src/lr_apply_liouvillian.f90
View File

@ -45,7 +45,8 @@ subroutine lr_apply_liouvillian( evc1, evc1_new, sevc1_new, interaction )
! !
implicit none implicit none
! !
complex(kind=dp) :: evc1(npwx,nbnd,nks), evc1_new(npwx,nbnd,nks), sevc1_new(npwx,nbnd,nks) complex(kind=dp),intent(in) :: evc1(npwx,nbnd,nks)
complex(kind=dp),intent(out) :: evc1_new(npwx,nbnd,nks), sevc1_new(npwx,nbnd,nks)
logical, intent(in) :: interaction logical, intent(in) :: interaction
! !
! Local variables ! Local variables

91
TDDFPT/src/lr_dvpsi_e.f90
View File

@ -52,12 +52,15 @@ subroutine lr_dvpsi_e(ik,ipol,dvpsi)
USE lr_variables, ONLY : lr_verbosity, evc0 USE lr_variables, ONLY : lr_verbosity, evc0
USE io_global, ONLY : stdout USE io_global, ONLY : stdout
!DEBUG
USE lr_variables, ONLY: check_all_bands_gamma, check_density_gamma,check_vector_gamma,check_vector_f
!
!
implicit none implicit none
! !
integer, intent(IN) :: ipol, ik integer, intent(IN) :: ipol, ik
! !
complex(kind=dp) :: dvpsi(npwx,nbnd) complex(kind=dp),intent(out) :: dvpsi(npwx,nbnd)
real(kind=dp) :: atnorm real(kind=dp) :: atnorm
! !
@ -112,9 +115,7 @@ subroutine lr_dvpsi_e(ik,ipol,dvpsi)
!OBM!!!! eprec is also calculated on the fly for each k point !OBM!!!! eprec is also calculated on the fly for each k point
allocate(eprec(nbnd)) allocate(eprec(nbnd))
!OBM!!!! !OBM!!!!
evc(:,:)=evc0(:,:,ik) evc(:,:)=evc0(:,:,ik)
!print *, "evc", evc(1:3,1)
!print *,"npw_k(",ik,")=",npw_k(ik)
if (nkb > 0) then if (nkb > 0) then
allocate (dvkb (npwx, nkb), dvkb1(npwx, nkb)) allocate (dvkb (npwx, nkb), dvkb1(npwx, nkb))
dvkb (:,:) = (0.d0, 0.d0) dvkb (:,:) = (0.d0, 0.d0)
@ -125,7 +126,12 @@ subroutine lr_dvpsi_e(ik,ipol,dvpsi)
gk (2, ig) = (xk (2, ik) + g (2, igk (ig) ) ) * tpiba gk (2, ig) = (xk (2, ik) + g (2, igk (ig) ) ) * tpiba
gk (3, ig) = (xk (3, ik) + g (3, igk (ig) ) ) * tpiba gk (3, ig) = (xk (3, ik) + g (3, igk (ig) ) ) * tpiba
g2kin (ig) = gk (1, ig) **2 + gk (2, ig) **2 + gk (3, ig) **2 g2kin (ig) = gk (1, ig) **2 + gk (2, ig) **2 + gk (3, ig) **2
enddo enddo
if (lr_verbosity > 10 ) then
write(stdout,'("lr_dvpsi_e g2kin:",F15.8)') SUM(g2kin(:))
endif
! !
! this is the kinetic contribution to [H,x]: -2i (k+G)_ipol * psi ! this is the kinetic contribution to [H,x]: -2i (k+G)_ipol * psi
! !
@ -153,7 +159,12 @@ subroutine lr_dvpsi_e(ik,ipol,dvpsi)
! ! ! !
! enddo ! enddo
! print *, "lr_dvpsi_e d0psi kinetic contribution", obm_debug ! print *, "lr_dvpsi_e d0psi kinetic contribution", obm_debug
if (lr_verbosity > 10 ) then
write(stdout,'("lr_dvpsi_e d0psi kinetic contribution:")')
do ibnd=1,nbnd
call check_vector_gamma(d0psi(:,ibnd))
enddo
endif
!!obm_debug !!obm_debug
@ -168,7 +179,38 @@ subroutine lr_dvpsi_e(ik,ipol,dvpsi)
! and this is the contribution from nonlocal pseudopotentials ! and this is the contribution from nonlocal pseudopotentials
! !
call gen_us_dj (ik, dvkb) call gen_us_dj (ik, dvkb)
call gen_us_dy (ik, at (1, ipol), dvkb1) call gen_us_dy (ik, at (1, ipol), dvkb1)
if (lr_verbosity > 10 ) then
write(stdout,'("lr_dvpsi_e dvkb:")')
jkb=0
do nt = 1, ntyp
do na = 1, nat
if (nt == ityp (na)) then
do ikb = 1, nh (nt)
jkb = jkb + 1
call check_vector_f(dvkb(:,jkb))
enddo
endif
enddo
enddo
endif
if (lr_verbosity > 10 ) then
write(stdout,'("lr_dvpsi_e dvkb1:")')
jkb=0
do nt = 1, ntyp
do na = 1, nat
if (nt == ityp (na)) then
do ikb = 1, nh (nt)
jkb = jkb + 1
call check_vector_f(dvkb1(:,jkb))
enddo
endif
enddo
enddo
endif
do ig = 1, npw_k(ik) do ig = 1, npw_k(ik)
if (g2kin (ig) < 1.0d-10) then if (g2kin (ig) < 1.0d-10) then
gk (1, ig) = 0.d0 gk (1, ig) = 0.d0
@ -201,9 +243,23 @@ subroutine lr_dvpsi_e(ik,ipol,dvpsi)
deallocate (gk) deallocate (gk)
!OBM!!!be careful, from bwalker, why?!!!!! !OBM!!!be careful, from bwalker, why?!!!!!
work(:,:)=(0.0d0,1.0d0)*work(:,:) work(:,:)=(0.0d0,1.0d0)*work(:,:)
!OBM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !OBM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!print *, "nonlocal contribution"!, d0psi(1:3,1) if (lr_verbosity > 10 ) then
!CALL lr_normalise( d0psi(:,:), anorm) write(stdout,'("lr_dvpsi_e non-local contribution:")')
jkb=0
do nt = 1, ntyp
do na = 1, nat
if (nt == ityp (na)) then
do ikb = 1, nh (nt)
jkb = jkb + 1
call check_vector_gamma(work(:,jkb))
enddo
endif
enddo
enddo
endif
if(gamma_only) then if(gamma_only) then
call lr_dvpsi_e_gamma() call lr_dvpsi_e_gamma()
else if (noncolin) then else if (noncolin) then
@ -219,7 +275,18 @@ subroutine lr_dvpsi_e(ik,ipol,dvpsi)
! ! ! !
! enddo ! enddo
! print *, "lr_dvpsi_e norm of dvpsi being returned=", obm_debug ! print *, "lr_dvpsi_e norm of dvpsi being returned=", obm_debug
! !
if (lr_verbosity > 10 ) then
write(stdout,'("lr_dvpsi_e d0psi after lr_dvpsi_e case specific calc:")')
do ibnd=1,nbnd
call check_vector_gamma(d0psi(:,ibnd))
enddo
write(stdout,'("lr_dvpsi_e dvpsii after lr_dvpsi_e case specific calc:")')
do ibnd=1,nbnd
call check_vector_gamma(dvpsi(:,ibnd))
enddo
endif
!!obm_debug !!obm_debug
IF (nkb > 0) THEN IF (nkb > 0) THEN

85
TDDFPT/src/lr_lanczos.f90
View File

@ -38,7 +38,8 @@ contains
v_loc_psir, s_psir_gamma, igk_k,npw_k, real_space_debug v_loc_psir, s_psir_gamma, igk_k,npw_k, real_space_debug
USE lr_variables, ONLY : lr_verbosity, charge_response USE lr_variables, ONLY : lr_verbosity, charge_response
use charg_resp, only : w_T_beta_store,w_T,lr_calc_F use charg_resp, only : w_T_beta_store,w_T,lr_calc_F
!Debugging
USE lr_variables, ONLY: check_all_bands_gamma, check_density_gamma,check_vector_gamma
! !
implicit none implicit none
! !
@ -51,7 +52,7 @@ contains
! !
! Local variables ! Local variables
! !
real(kind=dp) :: alpha, beta, gamma real(kind=dp) :: alpha, beta, gamma, temp
! !
complex(kind=dp) :: zeta complex(kind=dp) :: zeta
! !
@ -61,14 +62,14 @@ contains
! !
If (lr_verbosity > 5) THEN If (lr_verbosity > 5) THEN
WRITE(stdout,'("<lr_lanczos_one_step>")') WRITE(stdout,'("<lr_lanczos_one_step>")')
endif
If (lr_verbosity > 10) THEN
print *, "Real space = ", real_space print *, "Real space = ", real_space
print *, "Real space debug ", real_space_debug print *, "Real space debug ", real_space_debug
print *, "TQR = ", tqr print *, "TQR = ", tqr
endif endif
! !
call start_clock('one_step') call start_clock('one_step')
!print *, "norm of d0psi", lr_dot(d0psi(1,1,1,1),d0psi(1,1,1,1))
pol_index=1 pol_index=1
if ( n_ipol /= 1 ) pol_index=LR_polarization if ( n_ipol /= 1 ) pol_index=LR_polarization
! !
@ -117,7 +118,28 @@ contains
evc1_new(:,:,:,2)=evc1_new(:,:,:,1) evc1_new(:,:,:,2)=evc1_new(:,:,:,1)
sevc1_new(:,:,:,2)=sevc1_new(:,:,:,1) sevc1_new(:,:,:,2)=sevc1_new(:,:,:,1)
endif endif
LR_iteration=1 LR_iteration=1
!Debugging
if (lr_verbosity >10) then
!write(stdout,'("evc1_new(1), first step")')
!do ibnd=1,nbnd
! call check_vector_gamma(evc1_new(:,ibnd,1,1))
!enddo
!write(stdout,'("sevc1_new(1), first step")')
!do ibnd=1,nbnd
! call check_vector_gamma(sevc1_new(:,ibnd,1,1))
!enddo
!write(stdout,'("evc1_new(2), first step")')
!do ibnd=1,nbnd
! call check_vector_gamma(evc1_new(:,ibnd,1,2))
!enddo
!write(stdout,'("sevc1_new(2), first step")')
!do ibnd=1,nbnd
! call check_vector_gamma(sevc1_new(:,ibnd,1,2))
!enddo
temp=dble(lr_dot(evc1_new(1,1,1,1),sevc1_new(1,1,1,2)))
write(stdout,'("<evc1_new(1)|sevc1_new(2)> first step",E15.8)') temp
endif
endif endif
! !
! The lanczos algorithm starts here ! The lanczos algorithm starts here
@ -125,12 +147,6 @@ contains
! !
! Left and right vectors are orthogonalised wrto ground state wf ! Left and right vectors are orthogonalised wrto ground state wf
!print *, "norm of evc1,1 before lr_ortho", lr_dot(evc1_new(1,1,1,1),evc1_new(1,1,1,1))
!print *, "norm of evc1,2 before lr_ortho", lr_dot(evc1_new(1,1,1,2),evc1_new(1,1,1,2))
!print *, "norm of evc0 before lr_ortho", lr_dot(evc0(1,1,1),evc0(1,1,1))
!print *, "norm of sevc0 before lr_ortho", lr_dot(sevc0(1,1,1),sevc0(1,1,1))
!print *,"nks=", nks
!OBM: Notice that here "orthogonalization" is not strictly the true word, as the norm of the vectors change !OBM: Notice that here "orthogonalization" is not strictly the true word, as the norm of the vectors change
!This is due to how the uspp scheme is implemented, the beta are evc1_new(left).sevc1_new(right), that is, !This is due to how the uspp scheme is implemented, the beta are evc1_new(left).sevc1_new(right), that is,
!a mixing of two vectors, thus the resultant vector from belov should be devoid from S, which affects the norm !a mixing of two vectors, thus the resultant vector from belov should be devoid from S, which affects the norm
@ -139,10 +155,12 @@ contains
call lr_ortho(evc1_new(:,:,ik,1), evc0(:,:,ik), ik, ik, sevc0(:,:,ik),.true.) call lr_ortho(evc1_new(:,:,ik,1), evc0(:,:,ik), ik, ik, sevc0(:,:,ik),.true.)
call lr_ortho(evc1_new(:,:,ik,2), evc0(:,:,ik), ik, ik, sevc0(:,:,ik),.true.) call lr_ortho(evc1_new(:,:,ik,2), evc0(:,:,ik), ik, ik, sevc0(:,:,ik),.true.)
enddo enddo
if (lr_verbosity >10) then
!print *, "norm of evc1,1 after lr_ortho", lr_dot(evc1_new(1,1,1,1),evc1_new(1,1,1,1)) temp=dble(lr_dot(evc1_new(1,1,1,1),sevc0(1,1,1)))
!print *, "norm of evc1,2 after lr_ortho", lr_dot(evc1_new(1,1,1,1),evc1_new(1,1,1,2)) write(stdout,'("<evc1_new(1)|sevc0>",E15.8)') temp
temp=dble(lr_dot(evc1_new(1,1,1,2),sevc0(1,1,1)))
write(stdout,'("<evc1_new(2)|sevc0>",E15.8)') temp
endif
! !
! By construction <p|Lq>=0 should be 0, forcing this both conserves resources and increases stability ! By construction <p|Lq>=0 should be 0, forcing this both conserves resources and increases stability
@ -215,7 +233,6 @@ contains
end if end if
!print *, "norm of sevc1,1 after spsi", lr_dot(sevc1_new(1,1,1,1),sevc1_new(1,1,1,1)) !print *, "norm of sevc1,1 after spsi", lr_dot(sevc1_new(1,1,1,1),sevc1_new(1,1,1,1))
!print *, "norm of sevc1,2 after spsi", lr_dot(sevc1_new(1,1,1,1),sevc1_new(1,1,1,2)) !print *, "norm of sevc1,2 after spsi", lr_dot(sevc1_new(1,1,1,1),sevc1_new(1,1,1,2))
!Resume the LR !Resume the LR
! !
! Orthogonality requirement as proposed by Y. Saad beta=sqrt(|qdash.pdash|) gamma=sign(qdash.pdash)*beta ! Orthogonality requirement as proposed by Y. Saad beta=sqrt(|qdash.pdash|) gamma=sign(qdash.pdash)*beta
@ -276,6 +293,18 @@ contains
evc1_new(:,:,:,:)=(0.0d0,0.0d0) evc1_new(:,:,:,:)=(0.0d0,0.0d0)
sevc1_new(:,:,:,:)=(0.0d0,0.0d0) sevc1_new(:,:,:,:)=(0.0d0,0.0d0)
! !
if (lr_verbosity >10) then
write(stdout,'("evc1(1), rotate")')
do ibnd=1,nbnd
call check_vector_gamma(evc1(:,ibnd,1,1))
enddo
write(stdout,'("evc1(2), rotate")')
do ibnd=1,nbnd
call check_vector_gamma(evc1(:,ibnd,1,2))
enddo
endif
! !
! !
if(.not.ltammd) then if(.not.ltammd) then
@ -294,7 +323,17 @@ contains
call zcopy(size_evc,evc1_new(:,:,:,1),1,evc1_new(:,:,:,2),1) !evc1_new(,1) = evc1_new(,2) call zcopy(size_evc,evc1_new(:,:,:,1),1,evc1_new(:,:,:,2),1) !evc1_new(,1) = evc1_new(,2)
call zcopy(size_evc,evc1_new(:,:,:,1),1,evc1_new(:,:,:,2),1) !evc1_new(,1) = evc1_new(,2) call zcopy(size_evc,evc1_new(:,:,:,1),1,evc1_new(:,:,:,2),1) !evc1_new(,1) = evc1_new(,2)
end if end if
if (lr_verbosity >10) then
write(stdout,'("evc1(1), apply L")')
do ibnd=1,nbnd
call check_vector_gamma(evc1_new(:,ibnd,1,1))
enddo
write(stdout,'("evc1(2), apply L")')
do ibnd=1,nbnd
call check_vector_gamma(evc1_new(:,ibnd,1,2))
enddo
endif
! !
! qdash(i+1)=f(q(i))-gamma*q(i-1) ! qdash(i+1)=f(q(i))-gamma*q(i-1)
! pdash(i+1)=f(p(i))-beta*p(i-1) ! pdash(i+1)=f(p(i))-beta*p(i-1)
@ -302,7 +341,17 @@ contains
! !
!OBM BLAS !OBM BLAS
call zaxpy(size_evc,-cmplx(gamma,0.0d0,kind=dp),evc1_old(:,:,:,1),1,evc1_new(:,:,:,1),1) call zaxpy(size_evc,-cmplx(gamma,0.0d0,kind=dp),evc1_old(:,:,:,1),1,evc1_new(:,:,:,1),1)
call zaxpy(size_evc,-cmplx(beta,0.0d0,kind=dp),evc1_old(:,:,:,2),1,evc1_new(:,:,:,2),1) call zaxpy(size_evc,-cmplx(beta,0.0d0,kind=dp),evc1_old(:,:,:,2),1,evc1_new(:,:,:,2),1)
if (lr_verbosity >10) then
write(stdout,'("evc1(1), final")')
do ibnd=1,nbnd
call check_vector_gamma(evc1_new(:,ibnd,1,1))
enddo
write(stdout,'("evc1(2), final")')
do ibnd=1,nbnd
call check_vector_gamma(evc1_new(:,ibnd,1,2))
enddo
endif
! !
! Writing files for restart ! Writing files for restart
! !

33
TDDFPT/src/lr_main.f90
View File

@ -34,14 +34,15 @@ PROGRAM lr_main
USE control_ph, ONLY : nbnd_occ USE control_ph, ONLY : nbnd_occ
use wvfct, only : nbnd use wvfct, only : nbnd
!Debugging
USE lr_variables, ONLY: check_all_bands_gamma, check_density_gamma,check_vector_gamma
! !
IMPLICIT NONE IMPLICIT NONE
! !
! Local variables ! Local variables
! !
!INTEGER :: iter, ip, op !INTEGER :: iter, ip, op
INTEGER :: ip,pol_index,ibnd_occ,ibnd_virt INTEGER :: ip,pol_index,ibnd_occ,ibnd_virt,ibnd
INTEGER :: iter_restart,iteration INTEGER :: iter_restart,iteration
LOGICAL :: rflag LOGICAL :: rflag
CHARACTER (len=9) :: code = 'TDDFPT' CHARACTER (len=9) :: code = 'TDDFPT'
@ -89,8 +90,6 @@ PROGRAM lr_main
endif endif
!OBM_DEBUG !OBM_DEBUG
! !
!if (restart) call test_restart()
!
! !
!Initialisation of degauss/openshell related stuff !Initialisation of degauss/openshell related stuff
! !
@ -138,7 +137,7 @@ PROGRAM lr_main
CALL lr_read_d0psi() CALL lr_read_d0psi()
else else
CALL lr_solve_e() CALL lr_solve_e()
endif endif
! !
! Set up initial stuff for derivatives ! Set up initial stuff for derivatives
! !
@ -200,10 +199,28 @@ PROGRAM lr_main
! !
write(stdout,'(/5x,"Starting Lanczos loop",1x,i8)') LR_polarization write(stdout,'(/5x,"Starting Lanczos loop",1x,i8)') LR_polarization
! !
END IF END IF
! if (lr_verbosity >10) then
CALL sd0psi() write(stdout,'("d0psi")')
do ibnd=1,nbnd
call check_vector_gamma(d0psi(:,ibnd,1,pol_index))
enddo
endif
!
CALL sd0psi() !after this d0psi is Sd0psi
! !
if (lr_verbosity >10) then
write(stdout,'("initial evc1")')
do ibnd=1,nbnd
call check_vector_gamma(evc1(:,ibnd,1,1))
enddo
write(stdout,'("initial sd0psi")')
do ibnd=1,nbnd
call check_vector_gamma(d0psi(:,ibnd,1,pol_index))
enddo
endif
! !
lancz_loop1 : DO iteration = iter_restart, itermax lancz_loop1 : DO iteration = iter_restart, itermax
! !

169
TDDFPT/src/lr_ortho.f90
View File

@ -46,7 +46,7 @@ COMPLEX(DP), INTENT(INOUT) :: dvpsi(npwx*npol,nbnd)
COMPLEX(DP), INTENT(IN) :: sevc(npwx*npol,nbnd) ! work space allocated by COMPLEX(DP), INTENT(IN) :: sevc(npwx*npol,nbnd) ! work space allocated by
! the calling routine (was called dpsi) ! the calling routine (was called dpsi)
!real(kind=dp), intent(IN) :: lr_alpha_pv !This is calculated manually in tddfpt !real(kind=dp), intent(IN) :: lr_alpha_pv !This is calculated manually in tddfpt
logical, optional :: inverse !if .true. |dvspi> = |dvpsi> - |evq><sevc|dvpsi> instead of |dvspi> = |dvpsi> - |sevc><evq|dvpsi> logical, intent(in):: inverse !if .true. |dvspi> = |dvpsi> - |evq><sevc|dvpsi> instead of |dvspi> = |dvpsi> - |sevc><evq|dvpsi>
logical:: inverse_mode logical:: inverse_mode
@ -57,11 +57,11 @@ CALL start_clock ('lr_ortho')
If (lr_verbosity > 5) WRITE(stdout,'("<lr_ortho>")') If (lr_verbosity > 5) WRITE(stdout,'("<lr_ortho>")')
! !
if (.not. present(inverse)) then !if (.not. present(inverse)) then
inverse_mode=.false. ! inverse_mode=.false.
else !else
inverse_mode=inverse inverse_mode=inverse
endif !endif
if (gamma_only) then if (gamma_only) then
! !
call lr_ortho_gamma() call lr_ortho_gamma()
@ -205,116 +205,73 @@ contains
! !
if (lgauss) then if (lgauss) then
call errore ('lr_ortho', "degauss with gamma point algorithms",1) call errore ('lr_ortho', "degauss with gamma point algorithms",1)
!
! metallic case
!
ps = (0.d0, 0.d0)
if (inverse_mode) then
CALL ZGEMM( 'C', 'N', nbnd, nbnd_occ (ikk), npw_k(ikk), (1.d0,0.d0), &
sevc, npwx, dvpsi, npwx, (0.d0,0.d0), ps, nbnd )
else
CALL ZGEMM( 'C', 'N', nbnd, nbnd_occ (ikk), npw_k(ikk), (1.d0,0.d0), &
evq, npwx, dvpsi, npwx, (0.d0,0.d0), ps, nbnd )
endif
!
DO ibnd = 1, nbnd_occ (ikk)
wg1 = wgauss ((ef-et(ibnd,ikk)) / degauss, ngauss)
w0g = w0gauss((ef-et(ibnd,ikk)) / degauss, ngauss) / degauss
DO jbnd = 1, nbnd
wgp = wgauss ( (ef - et (jbnd, ikq) ) / degauss, ngauss)
deltae = et (jbnd, ikq) - et (ibnd, ikk)
theta = wgauss (deltae / degauss, 0)
wwg = wg1 * (1.d0 - theta) + wgp * theta
IF (jbnd <= nbnd_occ (ikq) ) THEN
IF (abs (deltae) > 1.0d-5) THEN
wwg = wwg + alpha_pv * theta * (wgp - wg1) / deltae
ELSE
!
! if the two energies are too close takes the limit
! of the 0/0 ratio
!
wwg = wwg - alpha_pv * theta * w0g
ENDIF
ENDIF
!
ps(jbnd,ibnd) = wwg * ps(jbnd,ibnd)
!
ENDDO
call DSCAL (2*npw_k(ikk), wg1, dvpsi(1,ibnd), 1)
END DO
nbnd_eff=nbnd
ELSE ELSE
! !
! insulators ! insulators
! ps = <evq|dvpsi> ! ps = <evq|dvpsi>
! in old version it was ps = <S evc0|sv> ! in old version it was ps = <S evc0|sv>
ps = 0.d0 ps = 0.d0
if (inverse_mode) then if (inverse_mode) then
CALL DGEMM( 'C', 'N', nbnd, nbnd ,2*npw_k(1), & CALL DGEMM( 'C', 'N', nbnd, nbnd ,2*npw_k(1), &
2.d0, sevc, 2*npwx, dvpsi, 2*npwx, & 2.d0, sevc, 2*npwx, dvpsi, 2*npwx, &
0.d0, ps, nbnd ) 0.d0, ps, nbnd )
else !ps = 2*<sevc|dvpsi>
CALL DGEMM( 'C', 'N', nbnd, nbnd ,2*npw_k(1), & else
2.d0, evq, 2*npwx, dvpsi, 2*npwx, & CALL DGEMM( 'C', 'N', nbnd, nbnd ,2*npw_k(1), &
0.d0, ps, nbnd ) 2.d0, evq, 2*npwx, dvpsi, 2*npwx, &
endif 0.d0, ps, nbnd )
nbnd_eff=nbnd !ps = 2*<evq|dvpsi>
if (gstart == 2) then endif
if (inverse_mode) then nbnd_eff=nbnd
CALL DGER( nbnd, nbnd, -1.D0, sevc, 2*npwx, dvpsi, 2*npwx, ps, nbnd ) if (gstart == 2) then
else if (inverse_mode) then
CALL DGER( nbnd, nbnd, -1.D0, evq, 2*npwx, dvpsi, 2*npwx, ps, nbnd ) CALL DGER( nbnd, nbnd, -1.D0, sevc, 2*npwx, dvpsi, 2*npwx, ps, nbnd )
endif !PS = PS - sevc*dvpsi
endif else
END IF CALL DGER( nbnd, nbnd, -1.D0, evq, 2*npwx, dvpsi, 2*npwx, ps, nbnd )
!PS = PS - evc*dvpsi
endif
endif
END IF
#ifdef __PARA #ifdef __PARA
call mp_sum(ps(:,1:nbnd_eff),intra_pool_comm) call mp_sum(ps(:,:),intra_pool_comm)
#endif #endif
! in the original dpsi was used as a storage for sevc, since in ! in the original dpsi was used as a storage for sevc, since in
! tddfpt we have it stored in memory as sevc0 this part is obsolote ! tddfpt we have it stored in memory as sevc0 this part is obsolote
!! !!
!! dpsi is used as work space to store S|evc> !! dpsi is used as work space to store S|evc>
!! !!
!IF (noncolin) THEN !IF (noncolin) THEN
! IF (okvan) CALL calbec ( npw_k(ikk), vkb, evq, becp_nc, nbnd_eff ) ! IF (okvan) CALL calbec ( npw_k(ikk), vkb, evq, becp_nc, nbnd_eff )
!ELSE !ELSE
! IF (okvan) CALL calbec ( npwq, vkb, evq, becp, nbnd_eff) ! IF (okvan) CALL calbec ( npwq, vkb, evq, becp, nbnd_eff)
!ENDIF !ENDIF
!CALL s_psi (npwx, npwq, nbnd_eff, evq, dpsi) !CALL s_psi (npwx, npwq, nbnd_eff, evq, dpsi)
ps_c = cmplx(ps, 0.d0, dp) ps_c = cmplx(ps, 0.d0, dp)
!! !!
!! |dvspi> = -(|dvpsi> - S|evq><evq|dvpsi>) !! |dvspi> = -(|dvpsi> - S|evq><evq|dvpsi>)
!! !!
!OBM!!! changed to |dvspi> = |dvpsi> - |sevc><evq|dvpsi> !OBM!!! changed to |dvspi> = |dvpsi> - |sevc><evq|dvpsi>
if (lgauss) then if (lgauss) then
! !errore ?
! metallic case ELSE
! !
if (inverse_mode) then ! Insulators: note that nbnd_occ(ikk)=nbnd_occ(ikq) in an insulator
CALL ZGEMM( 'N', 'N', npw_k(1), nbnd_occ(1), nbnd, & !
(1.d0,0.d0), evq, npwx, ps_c, nbnd, (-1.0d0,0.d0), & if (inverse_mode) then
dvpsi, npwx ) CALL ZGEMM( 'N', 'N', npw_k(1), nbnd, nbnd, &
else (-1.d0,0.d0), evq, npwx, ps_c, nbnd, (1.0d0,0.d0), &
CALL ZGEMM( 'N', 'N', npw_k(1), nbnd_occ(1), nbnd, & dvpsi, npwx )
(1.d0,0.d0), sevc, npwx, ps_c, nbnd, (-1.0d0,0.d0), & !dvpsi=dvpsi-|evq><sevc|dvpsi>
dvpsi, npwx ) else
endif CALL ZGEMM( 'N', 'N', npw_k(1), nbnd, nbnd, &
ELSE (-1.d0,0.d0), sevc, npwx, ps_c, nbnd, (1.0d0,0.d0), &
! dvpsi, npwx )
! Insulators: note that nbnd_occ(ikk)=nbnd_occ(ikq) in an insulator !dvpsi=dvpsi-|sevc><evq|dvpsi>
! endif
if (inverse_mode) then ENDIF
CALL ZGEMM( 'N', 'N', npw_k(1), nbnd, nbnd, &
(-1.d0,0.d0), evq, npwx, ps_c, nbnd, (1.0d0,0.d0), &
dvpsi, npwx )
else
CALL ZGEMM( 'N', 'N', npw_k(1), nbnd, nbnd, &
(-1.d0,0.d0), sevc, npwx, ps_c, nbnd, (1.0d0,0.d0), &
dvpsi, npwx )
endif
ENDIF
DEALLOCATE(ps) DEALLOCATE(ps)
DEALLOCATE(ps_c) DEALLOCATE(ps_c)
end subroutine lr_ortho_gamma end subroutine lr_ortho_gamma

25
TDDFPT/src/lr_read_wf.f90
View File

@ -73,6 +73,8 @@ subroutine normal_read()
! !
!The usual way of reading wavefunctions !The usual way of reading wavefunctions
! !
USE lr_variables, ONLY: check_all_bands_gamma, check_density_gamma,&
check_vector_gamma
IMPLICIT NONE IMPLICIT NONE
nwordwfc = 2 * nbnd * npwx nwordwfc = 2 * nbnd * npwx
@ -287,6 +289,14 @@ subroutine normal_read()
! !
!print * , "evc0 ",evc0(1:3,1,1) !print * , "evc0 ",evc0(1:3,1,1)
! !
if (lr_verbosity >10) then
call check_all_bands_gamma(evc0(:,:,1),sevc0(:,:,1),nbnd,nbnd)
write(stdout,'("evc0")')
do ibnd=1,nbnd
call check_vector_gamma(evc0(:,ibnd,1))
enddo
call check_density_gamma(revc0(:,:,1),nbnd)
endif
! !
!OBM!!! debug--- !OBM!!! debug---
!CALL lr_normalise( evc0(:,:,1), obm_debug) !CALL lr_normalise( evc0(:,:,1), obm_debug)
@ -306,6 +316,8 @@ subroutine virt_read()
! !
USE control_ph, ONLY : nbnd_occ USE control_ph, ONLY : nbnd_occ
use gvect, only : nrxx use gvect, only : nrxx
USE lr_variables, ONLY: check_all_bands_gamma, check_density_gamma,&
check_vector_gamma
IMPLICIT NONE IMPLICIT NONE
complex(kind=dp), allocatable :: evc_all(:,:,:) complex(kind=dp), allocatable :: evc_all(:,:,:)
complex(kind=dp), allocatable :: sevc_all(:,:,:) complex(kind=dp), allocatable :: sevc_all(:,:,:)
@ -423,7 +435,7 @@ use gvect, only : nrxx
! Inverse fourier transform of evc_all ! Inverse fourier transform of evc_all
! !
! !
revc0=(0.0d0,0.0d0) revc_all=(0.0d0,0.0d0)
! !
if ( gamma_only ) then if ( gamma_only ) then
! !
@ -481,8 +493,11 @@ use gvect, only : nrxx
! !
nbnd=nbnd_occ(1) nbnd=nbnd_occ(1)
! !
evc0=(0.0d0,0.0d0)
evc0(:,:,:)=evc_all(:,1:nbnd,:) evc0(:,:,:)=evc_all(:,1:nbnd,:)
sevc0=(0.0d0,0.0d0)
sevc0(:,:,:)=sevc_all(:,1:nbnd,:) sevc0(:,:,:)=sevc_all(:,1:nbnd,:)
revc0=(0.0d0,0.0d0)
revc0(:,:,:)=revc_all(:,1:nbnd,:) revc0(:,:,:)=revc_all(:,1:nbnd,:)
if (nkb>0) then if (nkb>0) then
if (gamma_only) then if (gamma_only) then
@ -506,6 +521,14 @@ use gvect, only : nrxx
endif endif
endif endif
endif endif
if (lr_verbosity >10) then
call check_all_bands_gamma(evc_all(:,:,1),sevc_all(:,:,1),nbnd_total,nbnd)
call check_density_gamma(revc_all(:,:,1),nbnd)
write(stdout,'("evc0")')
do ibnd=1,nbnd
call check_vector_gamma(evc0(:,ibnd,1))
enddo
endif
if (nkb>0) then if (nkb>0) then
if (gamma_only) then if (gamma_only) then
deallocate(becp1_all) deallocate(becp1_all)

130
TDDFPT/src/lr_variables.f90
View File

@ -170,5 +170,135 @@ module lr_variables
! !
!real(kind=dp) :: broadening !Broadening !real(kind=dp) :: broadening !Broadening
integer :: plot_type ! Dumps rho as: 1=xyzd 2=xsf 3=cube integer :: plot_type ! Dumps rho as: 1=xyzd 2=xsf 3=cube
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Debugging subroutines
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
contains
!----------------------------------------------------------------------------
SUBROUTINE check_vector_gamma (x)
! Checks the inner product for a given vector, and its imaginary and real
! component
! input, evc
! output : screen output
USE mp_global, ONLY : inter_pool_comm, intra_pool_comm
use mp, only : mp_sum
use realus, only : npw_k
use gvect, only : gstart
use io_global, only : stdout
IMPLICIT NONE
!input/output
complex(kind=dp),intent(in) :: x(:)
!local
real(kind=dp) :: temp_gamma
real(kind=dp), external :: DDOT
temp_gamma = 2.D0*DDOT(2*npw_k(1),x(:),1,x(:),1)
if (gstart==2) temp_gamma = temp_gamma - dble(x(1))*dble(x(1))
#ifdef __PARA
call mp_sum(temp_gamma, intra_pool_comm)
#endif
write(stdout,'("<x> = ",E15.8)') temp_gamma
END SUBROUTINE check_vector_gamma
!----------------------------------------------------------------------------
SUBROUTINE check_vector_f (x)
! Checks the inner product for a given vector, and its imaginary and real
! component
! input, evc
! output : screen output
USE mp_global, ONLY : inter_pool_comm, intra_pool_comm
use mp, only : mp_sum
use realus, only : npw_k
use gvect, only : gstart
use io_global, only : stdout
IMPLICIT NONE
!input/output
complex(kind=dp),intent(in) :: x(:)
!local
complex(kind=dp) :: temp_f
complex(kind=dp), external :: ZDOTC
temp_f = ZDOTC(npw_k(1),x(:),1,x(:),1)
#ifdef __PARA
call mp_sum(temp_f, intra_pool_comm)
#endif
write(stdout,'("<x> = ",2E15.8,1X)') temp_f
END SUBROUTINE check_vector_f
!----------------------------------------------------------------------------
SUBROUTINE check_all_bands_gamma (x,sx,nbnd1,nbnd2)
! Checks all bands of given KS states for orthoganilty
! input, evc and sevc
! output : screen output
USE mp_global, ONLY : inter_pool_comm, intra_pool_comm
use mp, only : mp_sum
use realus, only : npw_k
use io_global, only : stdout
use gvect, only : gstart
IMPLICIT NONE
!input/output
integer,intent(in) :: nbnd1,nbnd2 !Total number of bands for x and sx
complex(kind=dp),intent(in) :: x(:,:), sx(:,:)
!local
integer :: ibnd, jbnd
real(kind=dp) :: temp_gamma
real(kind=dp), external :: DDOT
do ibnd=1,nbnd1
do jbnd=ibnd,nbnd2
!
temp_gamma = 2.D0*DDOT(2*npw_k(1),x(:,ibnd),1,sx(:,jbnd),1)
if (gstart==2) temp_gamma = temp_gamma - dble(x(1,ibnd))*dble(sx(1,jbnd))
#ifdef __PARA
call mp_sum(temp_gamma, intra_pool_comm)
#endif
write(stdout,'("<x,",I02,"|S|x,",I02,"> =",E15.8)') ibnd,jbnd,temp_gamma
enddo
enddo
END SUBROUTINE check_all_bands_gamma
!----------------------------------------------------------------------------
SUBROUTINE check_density_gamma (rx,nbnd)
! Checks the contirbution of a given function transformed into real space
! input, revc
! output : screen output
USE mp_global, ONLY : inter_pool_comm, intra_pool_comm
use mp, only : mp_sum
use realus, only : npw_k
use wvfct, only : wg
use gvect, only : nrxx
use io_global, only : stdout
use cell_base, only : omega
IMPLICIT NONE
!input/output
integer,intent(in) :: nbnd !Total number of bands for x and sx
complex(kind=dp),intent(in) :: rx(:,:)
!local
integer :: ibnd
real(kind=dp) :: temp_gamma,w1,w2
do ibnd=1,nbnd,2
w1=wg(ibnd,1)/omega
!
if(ibnd<nbnd) then
w2=wg(ibnd+1,1)/omega
else
w2=w1
endif
temp_gamma=SUM(w1*DBLE(rx(1:nrxx,ibnd))*DBLE(rx(1:nrxx,ibnd))&
+w2*aimag(rx(1:nrxx,ibnd))*aimag(rx(1:nrxx,ibnd)))
#ifdef __PARA
call mp_sum(temp_gamma, intra_pool_comm)
#endif
write(stdout,'("Contribution of bands ",I02," and ",I02," to total density",E15.8)') ibnd,ibnd+1,temp_gamma
enddo
!
!
END SUBROUTINE check_density_gamma
!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
end module lr_variables end module lr_variables
!---------------------------------------------------------------------------- !----------------------------------------------------------------------------