fix bug in TDDFPT/src/lr_alloc_init

fix pseudo_hermitian code
This commit is contained in:
Oscar Baseggio 2021-01-15 17:13:49 +01:00 committed by Baseggio
parent f85ba9f344
commit d7ca8fa110
3 changed files with 12 additions and 8 deletions

View File

@ -67,7 +67,7 @@ lr_dav_debug.o\
plugin_tddfpt_potential.o\
linear_solvers.o\
orthogonalize_omega.o\
lr_sternheimer.o
lr_sternheimer.o\
lr_apply_liouvillian_magnons.o\
lr_dvpsi_magnons.o\
lr_Opsi_magnons.o\

View File

@ -66,13 +66,9 @@ subroutine lanczos_nonhermitian_c(j, npwx_npol, nbnd_occ, nksq, qj_r, Aqj_r, &
!
CALL zscal(size_evc,cmplx(1.0d0/beta,0.0d0,kind=dp),qj_r(1,1,1,1),1)
CALL zscal(size_evc,cmplx(1.0d0/beta,0.0d0,kind=dp),Aqj_r(1,1,1,1),1)
! qj_r(:,:,:,:) = qj_r(:,:,:,:) / cmplx(beta,0.0d0,kind=dp)
! Aqj_r(:,:,:,:) = Aqj_r(:,:,:,:) / cmplx(beta,0.0d0,kind=dp)
!
CALL zscal(size_evc,1.0d0/conjg(gamma),qj_l(1,1,1,1),1)
CALL zscal(size_evc,1.0d0/conjg(gamma),Aqj_l(1,1,1,1),1)
! qj_l(:,:,:,:) = qj_l(:,:,:,:) / conjg(gamma)
! Aqj_l(:,:,:,:) = Aqj_l(:,:,:,:) / conjg(gamma)
!
! computes alpha
!

View File

@ -38,6 +38,8 @@ subroutine lanczos_pseudohermitian_c(j, npwx_npol, nbnd_occ, nksq, qj_r, Aqj_r,
!! upper coefficient of the tridiagonal matrix
COMPLEX(kind=dp), INTENT(OUT) :: zeta(n_ipol)
!! (u,q_j) products
! dummy variable
COMPLEX(kind=dp) :: Aqj_l(npwx_npol,nbnd_occ,nksq, 2)
!
COMPLEX(kind=dp),EXTERNAL :: lr_dot_magnons
!
@ -60,13 +62,19 @@ subroutine lanczos_pseudohermitian_c(j, npwx_npol, nbnd_occ, nksq, qj_r, Aqj_r,
!
CALL zscal(size_evc,cmplx(1.0d0/beta,0.0d0,kind=dp),qj_r(1,1,1,1),1)
CALL zscal(size_evc,cmplx(1.0d0/beta,0.0d0,kind=dp),Aqj_r(1,1,1,1),1)
! qj_r(:,:,:,:) = qj_r(:,:,:,:) / cmplx(beta,0.0d0,kind=dp)
! Aqj_r(:,:,:,:) = Aqj_r(:,:,:,:) / cmplx(beta,0.0d0,kind=dp)
!
!
! computes alpha
!
alpha = (0.0d0, 0.0d0)
Aqj_l = Aqj_r
Aqj_l(:,:,:,2) = -Aqj_l(:,:,:,2)
!
alpha = lr_dot_magnons(Aqj_l, Aqj_r)
!
IF (force_zero_alpha) THEN
alpha = (0.0d0,0.0d0)
ENDIF
!
! Calculation of zeta coefficients.
! See Eq.(35) in Malcioglu et al., Comput. Phys. Commun. 182, 1744 (2011).
@ -85,7 +93,7 @@ subroutine lanczos_pseudohermitian_c(j, npwx_npol, nbnd_occ, nksq, qj_r, Aqj_r,
! Renormalization will be done in the begining of the next iteration.
! In the non-Hermitian case, similar operation needs to be done also for p(i).
!
! CALL zaxpy(size_evc,-alpha,qj_r(1,1,1,1),1,Aqj_r(1,1,1,1),1)
CALL zaxpy(size_evc,-alpha,qj_r(1,1,1,1),1,Aqj_r(1,1,1,1),1)
CALL zaxpy(size_evc,-gamma,qjold_r(1,1,1,1),1,Aqj_r(1,1,1,1),1)
!
! X. Ge: Throw away q(i-1), and make q(i+1) to be the current vector,