PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-09-27 12:00:50 +02:00
Code Releases Activity

fix bug in sort_ppRPA

Browse Source
This commit is contained in:
Pierre-Francois Loos 2024-09-12 17:14:04 +02:00
parent ff7cff0963
commit 6290634d87
8 changed files with 54 additions and 53 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/GT/RG0T0pp.f90
View File

@ -40,7 +40,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
! Local variables ! Local variables
logical :: print_T = .false. logical :: print_T = .true.
double precision :: lambda double precision :: lambda
integer :: isp_T integer :: isp_T
@ -173,6 +173,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,lambda,ERI,Bpp) if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,lambda,ERI,Bpp)
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T)) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
deallocate(Bpp,Cpp,Dpp) deallocate(Bpp,Cpp,Dpp)
!print*, 'LAPACK:' !print*, 'LAPACK:'
!print*, Om2t !print*, Om2t

2
src/GT/RGTeh_QP_graph.f90
View File

@ -42,7 +42,7 @@ subroutine RGTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eOld,e
! Run Newton's algorithm to find the root ! Run Newton's algorithm to find the root
write(*,*)'-----------------------------------------------------' write(*,*)'-----------------------------------------------------'
write(*,'(A5,1X,A3,1X,A15,1X,A15,1X,A10)') 'Orb.','It.','e_GTehlin (eV)','e_GTehlin (eV)','Z' write(*,'(A5,1X,A3,1X,A15,1X,A15,1X,A10)') 'Orb.','It.','e_GTehlin (eV)','e_GTeh (eV)','Z'
write(*,*)'-----------------------------------------------------' write(*,*)'-----------------------------------------------------'
do p=nC+1,nBas-nR do p=nC+1,nBas-nR

2
src/GT/RGTpp_QP_graph.f90
View File

@ -45,7 +45,7 @@ subroutine RGTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1
! Run Newton's algorithm to find the root ! Run Newton's algorithm to find the root
write(*,*)'-----------------------------------------------------' write(*,*)'-----------------------------------------------------'
write(*,'(A5,1X,A3,1X,A15,1X,A15,1X,A10)') 'Orb.','It.','e_GTpplin (eV)','e_GTpplin (eV)','Z' write(*,'(A5,1X,A3,1X,A15,1X,A15,1X,A10)') 'Orb.','It.','e_GTpplin (eV)','e_GTpp (eV)','Z'
write(*,*)'-----------------------------------------------------' write(*,*)'-----------------------------------------------------'
do p=nC+1,nBas-nR do p=nC+1,nBas-nR

1
src/GW/RGW_ppBSE.f90
View File

@ -76,7 +76,6 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
double precision,intent(out) :: EcBSE(nspin) double precision,intent(out) :: EcBSE(nspin)
!--------------------------------- !---------------------------------
! Compute (singlet) RPA screening ! Compute (singlet) RPA screening
!--------------------------------- !---------------------------------

30
src/GW/RGW_ppBSE_dynamic_perturbation.f90
View File

@ -1,4 +1,4 @@
subroutine RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int, & subroutine RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int, &
OmRPA,rho_RPA,Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta) OmRPA,rho_RPA,Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta)
! Compute dynamical effects via perturbation theory for BSE ! Compute dynamical effects via perturbation theory for BSE
@ -11,7 +11,7 @@ subroutine RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO
integer,intent(in) :: ispin integer,intent(in) :: ispin
logical,intent(in) :: dTDA logical,intent(in) :: dTDA
double precision,intent(in) :: eta double precision,intent(in) :: eta
integer,intent(in) :: nBas integer,intent(in) :: nOrb
integer,intent(in) :: nC integer,intent(in) :: nC
integer,intent(in) :: nO integer,intent(in) :: nO
integer,intent(in) :: nV integer,intent(in) :: nV
@ -20,12 +20,12 @@ subroutine RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO
integer,intent(in) :: nOO integer,intent(in) :: nOO
integer,intent(in) :: nVV integer,intent(in) :: nVV
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eW(nBas) double precision,intent(in) :: eW(nOrb)
double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: eGW(nOrb)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart)
double precision,intent(in) :: OmRPA(nS) double precision,intent(in) :: OmRPA(nS)
double precision,intent(in) :: rho_RPA(nBas,nBas,nS) double precision,intent(in) :: rho_RPA(nOrb,nOrb,nS)
double precision,intent(in) :: Om1(nVV) double precision,intent(in) :: Om1(nVV)
double precision,intent(in) :: X1(nVV,nVV) double precision,intent(in) :: X1(nVV,nVV)
double precision,intent(in) :: Y1(nOO,nVV) double precision,intent(in) :: Y1(nOO,nVV)
@ -76,16 +76,16 @@ subroutine RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO
if(dTDA) then if(dTDA) then
call RGW_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,eGW,OmRPA,rho_RPA,Om1(ab),KC_dyn,ZC_dyn) call RGW_ppBSE_dynamic_kernel_C(ispin,eta,nOrb,nC,nO,nV,nR,nS,nVV,1d0,eGW,OmRPA,rho_RPA,Om1(ab),KC_dyn,ZC_dyn)
Z1_dyn(ab) = + dot_product(X1(:,ab),matmul(ZC_dyn,X1(:,ab))) Z1_dyn(ab) = + dot_product(X1(:,ab),matmul(ZC_dyn,X1(:,ab)))
Om1_dyn(ab) = + dot_product(X1(:,ab),matmul(KC_dyn - KC_sta,X1(:,ab))) Om1_dyn(ab) = + dot_product(X1(:,ab),matmul(KC_dyn - KC_sta,X1(:,ab)))
else else
call RGW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGW,OmRPA,rho_RPA,KB_dyn) call RGW_ppBSE_dynamic_kernel_B(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGW,OmRPA,rho_RPA,KB_dyn)
call RGW_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,eGW,OmRPA,rho_RPA,Om1(ab),KC_dyn,ZC_dyn) call RGW_ppBSE_dynamic_kernel_C(ispin,eta,nOrb,nC,nO,nV,nR,nS,nVV,1d0,eGW,OmRPA,rho_RPA,Om1(ab),KC_dyn,ZC_dyn)
call RGW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,eGW,OmRPA,rho_RPA,Om1(ab),KD_dyn,ZD_dyn) call RGW_ppBSE_dynamic_kernel_D(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,1d0,eGW,OmRPA,rho_RPA,Om1(ab),KD_dyn,ZD_dyn)
Z1_dyn(ab) = dot_product(X1(:,ab),matmul(ZC_dyn,X1(:,ab))) & Z1_dyn(ab) = dot_product(X1(:,ab),matmul(ZC_dyn,X1(:,ab))) &
+ dot_product(Y1(:,ab),matmul(ZD_dyn,Y1(:,ab))) + dot_product(Y1(:,ab),matmul(ZD_dyn,Y1(:,ab)))
@ -119,16 +119,16 @@ subroutine RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO
if(dTDA) then if(dTDA) then
call RGW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,eGW,OmRPA,rho_RPA,Om2(ij),KD_dyn,ZD_dyn) call RGW_ppBSE_dynamic_kernel_D(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,1d0,eGW,OmRPA,rho_RPA,Om2(ij),KD_dyn,ZD_dyn)
Z2_dyn(kl) = + dot_product(Y2(:,ij),matmul(ZD_dyn,Y2(:,ij))) Z2_dyn(kl) = + dot_product(Y2(:,ij),matmul(ZD_dyn,Y2(:,ij)))
Om2_dyn(kl) = - dot_product(Y2(:,ij),matmul(KD_dyn - KD_sta,Y2(:,ij))) Om2_dyn(kl) = - dot_product(Y2(:,ij),matmul(KD_dyn - KD_sta,Y2(:,ij)))
else else
call RGW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGW,OmRPA,rho_RPA,KB_dyn) call RGW_ppBSE_dynamic_kernel_B(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGW,OmRPA,rho_RPA,KB_dyn)
call RGW_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,eGW,OmRPA,rho_RPA,Om2(ij),KC_dyn,ZC_dyn) call RGW_ppBSE_dynamic_kernel_C(ispin,eta,nOrb,nC,nO,nV,nR,nS,nVV,1d0,eGW,OmRPA,rho_RPA,Om2(ij),KC_dyn,ZC_dyn)
call RGW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,eGW,OmRPA,rho_RPA,Om2(ij),KD_dyn,ZD_dyn) call RGW_ppBSE_dynamic_kernel_D(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,1d0,eGW,OmRPA,rho_RPA,Om2(ij),KD_dyn,ZD_dyn)
Z2_dyn(kl) = dot_product(X2(:,ij),matmul(ZC_dyn,X2(:,ij))) & Z2_dyn(kl) = dot_product(X2(:,ij),matmul(ZC_dyn,X2(:,ij))) &
+ dot_product(Y2(:,ij),matmul(ZD_dyn,Y2(:,ij))) + dot_product(Y2(:,ij),matmul(ZD_dyn,Y2(:,ij)))

58
src/LR/ppLR.f90
View File

@ -68,43 +68,43 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
M( 1:nVV ,nVV+1:nOO+nVV) = - Bpp(1:nVV,1:nOO) M( 1:nVV ,nVV+1:nOO+nVV) = - Bpp(1:nVV,1:nOO)
M(nVV+1:nOO+nVV, 1:nVV) = + transpose(Bpp(1:nVV,1:nOO)) M(nVV+1:nOO+nVV, 1:nVV) = + transpose(Bpp(1:nVV,1:nOO))
if((nOO .eq. 0) .or. (nVV .eq. 0)) then ! if((nOO .eq. 0) .or. (nVV .eq. 0)) then
! Diagonalize the p-p matrix ! Diagonalize the p-p matrix
if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV, M, Om, Z) if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV, M, Om, Z)
! Split the various quantities in p-p and h-h parts ! Split the various quantities in p-p and h-h parts
call sort_ppRPA(nOO, nVV, Om, Z, Om1, X1, Y1, Om2, X2, Y2) call sort_ppRPA(nOO, nVV, Om, Z, Om1, X1, Y1, Om2, X2, Y2)
else ! else
thr_d = 1d-6 ! to determine if diagonal elements of L.T x R are close enouph to 1 ! thr_d = 1d-6 ! to determine if diagonal elements of L.T x R are close enouph to 1
thr_nd = 1d-6 ! to determine if non-diagonal elements of L.T x R are close enouph to 1 ! thr_nd = 1d-6 ! to determine if non-diagonal elements of L.T x R are close enouph to 1
thr_deg = 1d-8 ! to determine if two eigenvectors are degenerate or not ! thr_deg = 1d-8 ! to determine if two eigenvectors are degenerate or not
imp_bio = .True. ! impose bi-orthogonality ! imp_bio = .True. ! impose bi-orthogonality
verbose = .False. ! verbose = .False.
call diagonalize_nonsym_matrix(N, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose) ! call diagonalize_nonsym_matrix(N, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose)
!
! do i = 1, nOO
! Om2(i) = Om(i)
! do j = 1, nVV
! X2(j,i) = Z(j,i)
! enddo
! do j = 1, nOO
! Y2(j,i) = Z(nVV+j,i)
! enddo
! enddo
!
! do i = 1, nVV
! Om1(i) = Om(nOO+i)
! do j = 1, nVV
! X1(j,i) = M(j,nOO+i)
! enddo
! do j = 1, nOO
! Y1(j,i) = M(nVV+j,nOO+i)
! enddo
! enddo
do i = 1, nOO ! endif
Om2(i) = Om(i)
do j = 1, nVV
X2(j,i) = Z(j,i)
enddo
do j = 1, nOO
Y2(j,i) = Z(nVV+j,i)
enddo
enddo
do i = 1, nVV
Om1(i) = Om(nOO+i)
do j = 1, nVV
X1(j,i) = M(j,nOO+i)
enddo
do j = 1, nOO
Y1(j,i) = M(nVV+j,nOO+i)
enddo
enddo
endif
end if end if

5
src/RPA/sort_ppRPA.f90
View File

@ -17,6 +17,7 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
integer :: pq,ab,ij integer :: pq,ab,ij
! integer :: deg1,ab_start,ab_end ! integer :: deg1,ab_start,ab_end
! integer :: deg2,ij_start,ij_end ! integer :: deg2,ij_start,ij_end
integer :: OO,VV
double precision,allocatable :: M(:,:) double precision,allocatable :: M(:,:)
double precision,allocatable :: Z1(:,:) double precision,allocatable :: Z1(:,:)
double precision,allocatable :: Z2(:,:) double precision,allocatable :: Z2(:,:)
@ -210,8 +211,8 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
! S1 = + matmul(transpose(Z1),matmul(M,Z1)) ! S1 = + matmul(transpose(Z1),matmul(M,Z1))
! S2 = - matmul(transpose(Z2),matmul(M,Z2)) ! S2 = - matmul(transpose(Z2),matmul(M,Z2))
if(nVV > 0) call orthogonalization_matrix(nVV,S1,O1) if(nVV > 0) call orthogonalization_matrix(nVV,VV,S1,O1)
if(nOO > 0) call orthogonalization_matrix(nOO,S2,O2) if(nOO > 0) call orthogonalization_matrix(nOO,OO,S2,O2)
if(nVV > 0) call dgemm ('N', 'N', nOO+nVV,nVV,nVV, 1d0, Z1, nOO+nVV, O1, nVV,0d0, tmp1, nOO+nVV) if(nVV > 0) call dgemm ('N', 'N', nOO+nVV,nVV,nVV, 1d0, Z1, nOO+nVV, O1, nVV,0d0, tmp1, nOO+nVV)
Z1 = tmp1 Z1 = tmp1

2
src/utils/orthogonalization_matrix.f90
View File

@ -19,7 +19,7 @@ subroutine orthogonalization_matrix(nBas,nOrb,S,X)
! Output variables ! Output variables
integer :: nOrb integer,intent(out) :: nOrb
double precision,intent(out) :: X(nBas,nBas) double precision,intent(out) :: X(nBas,nBas)
debug = .false. debug = .false.