From 6290634d8783bf23b183fb06d6b5319a7a5c5056 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 12 Sep 2024 17:14:04 +0200 Subject: [PATCH] fix bug in sort_ppRPA --- src/GT/RG0T0pp.f90 | 5 +- src/GT/RGTeh_QP_graph.f90 | 2 +- src/GT/RGTpp_QP_graph.f90 | 2 +- src/GW/RGW_ppBSE.f90 | 1 - src/GW/RGW_ppBSE_dynamic_perturbation.f90 | 32 ++++++------- src/LR/ppLR.f90 | 58 +++++++++++------------ src/RPA/sort_ppRPA.f90 | 5 +- src/utils/orthogonalization_matrix.f90 | 2 +- 8 files changed, 54 insertions(+), 53 deletions(-) diff --git a/src/GT/RG0T0pp.f90 b/src/GT/RG0T0pp.f90 index 67dca13..474ea95 100644 --- a/src/GT/RG0T0pp.f90 +++ b/src/GT/RG0T0pp.f90 @@ -40,7 +40,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d ! Local variables - logical :: print_T = .false. + logical :: print_T = .true. double precision :: lambda 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) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T)) + deallocate(Bpp,Cpp,Dpp) !print*, 'LAPACK:' !print*, Om2t @@ -244,7 +245,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d write(*,*) call RGTpp_QP_graph(eta,nOrb,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, & - Om1t,rho1t,Om2t,rho2t,eGTlin,eHF,eGT,Z) + Om1t,rho1t,Om2t,rho2t,eGTlin,eHF,eGT,Z) end if diff --git a/src/GT/RGTeh_QP_graph.f90 b/src/GT/RGTeh_QP_graph.f90 index 2bfc163..f9e6bce 100644 --- a/src/GT/RGTeh_QP_graph.f90 +++ b/src/GT/RGTeh_QP_graph.f90 @@ -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 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(*,*)'-----------------------------------------------------' do p=nC+1,nBas-nR diff --git a/src/GT/RGTpp_QP_graph.f90 b/src/GT/RGTpp_QP_graph.f90 index 29d3387..692196c 100644 --- a/src/GT/RGTpp_QP_graph.f90 +++ b/src/GT/RGTpp_QP_graph.f90 @@ -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 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(*,*)'-----------------------------------------------------' do p=nC+1,nBas-nR diff --git a/src/GW/RGW_ppBSE.f90 b/src/GW/RGW_ppBSE.f90 index a823c2f..808caa9 100644 --- a/src/GW/RGW_ppBSE.f90 +++ b/src/GW/RGW_ppBSE.f90 @@ -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) - !--------------------------------- ! Compute (singlet) RPA screening !--------------------------------- diff --git a/src/GW/RGW_ppBSE_dynamic_perturbation.f90 b/src/GW/RGW_ppBSE_dynamic_perturbation.f90 index d7422a5..255cd85 100644 --- a/src/GW/RGW_ppBSE_dynamic_perturbation.f90 +++ b/src/GW/RGW_ppBSE_dynamic_perturbation.f90 @@ -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) ! 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 logical,intent(in) :: dTDA double precision,intent(in) :: eta - integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO 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) :: nVV - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: eW(nBas) - double precision,intent(in) :: eGW(nBas) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: eW(nOrb) + double precision,intent(in) :: eGW(nOrb) + double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart) 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) :: X1(nVV,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 - 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))) Om1_dyn(ab) = + dot_product(X1(:,ab),matmul(KC_dyn - KC_sta,X1(:,ab))) 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_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_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_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,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,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))) & + 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 - 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))) Om2_dyn(kl) = - dot_product(Y2(:,ij),matmul(KD_dyn - KD_sta,Y2(:,ij))) 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_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_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_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,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,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))) & + dot_product(Y2(:,ij),matmul(ZD_dyn,Y2(:,ij))) diff --git a/src/LR/ppLR.f90 b/src/LR/ppLR.f90 index 498e508..ff2469a 100644 --- a/src/LR/ppLR.f90 +++ b/src/LR/ppLR.f90 @@ -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(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 if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV, M, Om, Z) ! Split the various quantities in p-p and h-h parts 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_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 - imp_bio = .True. ! impose bi-orthogonality - verbose = .False. - 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 +! 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_deg = 1d-8 ! to determine if two eigenvectors are degenerate or not +! imp_bio = .True. ! impose bi-orthogonality +! verbose = .False. +! 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 - endif +! endif end if diff --git a/src/RPA/sort_ppRPA.f90 b/src/RPA/sort_ppRPA.f90 index 8518ecb..85bed67 100644 --- a/src/RPA/sort_ppRPA.f90 +++ b/src/RPA/sort_ppRPA.f90 @@ -17,6 +17,7 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) integer :: pq,ab,ij ! integer :: deg1,ab_start,ab_end ! integer :: deg2,ij_start,ij_end + integer :: OO,VV double precision,allocatable :: M(:,:) double precision,allocatable :: Z1(:,:) 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)) ! S2 = - matmul(transpose(Z2),matmul(M,Z2)) - if(nVV > 0) call orthogonalization_matrix(nVV,S1,O1) - if(nOO > 0) call orthogonalization_matrix(nOO,S2,O2) + if(nVV > 0) call orthogonalization_matrix(nVV,VV,S1,O1) + 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) Z1 = tmp1 diff --git a/src/utils/orthogonalization_matrix.f90 b/src/utils/orthogonalization_matrix.f90 index fd9e59c..66b3bb7 100644 --- a/src/utils/orthogonalization_matrix.f90 +++ b/src/utils/orthogonalization_matrix.f90 @@ -19,7 +19,7 @@ subroutine orthogonalization_matrix(nBas,nOrb,S,X) ! Output variables - integer :: nOrb + integer,intent(out) :: nOrb double precision,intent(out) :: X(nBas,nBas) debug = .false.