From 77fea9a9bd271e51bca0e7dfa1948b61e7e468c7 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 13 Sep 2024 09:33:43 +0200 Subject: [PATCH] redefine ppLR matrix --- src/LR/ppGLR.f90 | 127 +++++++++++++++++++++-------------------- src/LR/ppLR.f90 | 59 ++++++++++--------- src/LR/ppULR.f90 | 3 +- src/RPA/sort_ppRPA.f90 | 50 +++++++++------- 4 files changed, 126 insertions(+), 113 deletions(-) diff --git a/src/LR/ppGLR.f90 b/src/LR/ppGLR.f90 index 536015b..362fcd7 100644 --- a/src/LR/ppGLR.f90 +++ b/src/LR/ppGLR.f90 @@ -1,4 +1,4 @@ -subroutine ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) +subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) ! ! Solve the pp-RPA linear eigenvalue problem @@ -28,98 +28,103 @@ subroutine ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) implicit none include 'parameters.h' - logical, intent(in) :: TDA - integer, intent(in) :: nOO, nVV - double precision, intent(in) :: Bpp(nVV,nOO), Cpp(nVV,nVV), Dpp(nOO,nOO) - double precision, intent(out) :: Om1(nVV), X1(nVV,nVV), Y1(nOO,nVV) - double precision, intent(out) :: Om2(nOO), X2(nVV,nOO), Y2(nOO,nOO) - double precision, intent(out) :: EcRPA + logical,intent(in) :: TDA + integer,intent(in) :: nOO,nVV + double precision,intent(in) :: Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO) + double precision,intent(out) :: Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV) + double precision,intent(out) :: Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO) + double precision,intent(out) :: EcRPA - logical :: imp_bio, verbose - integer :: i, j, N - double precision :: EcRPA1, EcRPA2 - double precision :: thr_d, thr_nd, thr_deg - double precision,allocatable :: M(:,:), Z(:,:), Om(:) + logical :: imp_bio,verbose + integer :: i,j + integer :: Npp + double precision :: EcRPA1,EcRPA2 + double precision :: thr_d,thr_nd,thr_deg + double precision,allocatable :: M(:,:),Z(:,:),Om(:) - double precision, external :: trace_matrix + double precision,external :: trace_matrix + nPP = nOO + nVV + allocate(M(nPP,nPP),Z(nPP,nPP),Om(nPP)) - N = nOO + nVV - - allocate(M(N,N), Z(N,N), Om(N)) + ! Hermitian case for TDA if(TDA) then X1(:,:) = +Cpp(:,:) Y1(:,:) = 0d0 - if(nVV > 0) call diagonalize_matrix(nVV, X1, Om1) + if(nVV > 0) call diagonalize_matrix(nVV,X1,Om1) X2(:,:) = 0d0 Y2(:,:) = -Dpp(:,:) - if(nOO > 0) call diagonalize_matrix(nOO, Y2, Om2) + if(nOO > 0) call diagonalize_matrix(nOO,Y2,Om2) else ! Diagonal blocks - M( 1:nVV , 1:nVV) = + Cpp(1:nVV,1:nVV) - M(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - Dpp(1:nOO,1:nOO) + + M( 1:nVV, 1:nVV) = + Cpp(1:nVV,1:nVV) + M(nVV+1:nPP,nVV+1:nPP) = - Dpp(1:nOO,1:nOO) ! Off-diagonal blocks - 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 + M( 1:nVV,nVV+1:nPP) = + Bpp(1:nVV,1:nOO) + M(nVV+1:nPP, 1:nVV) = - transpose(Bpp(1:nVV,1:nOO)) - ! 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) +! if((nOO .eq. 0) .or. (nVV .eq. 0)) then - else + ! Diagonalize the pp matrix - 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 + if(nPP > 0) call diagonalize_general_matrix(nPP,M,Om,Z) - endif + ! Split the various quantities in ee and hh parts + + call sort_ppRPA(nOO,nVV,nPP,Om,Z,Om1,X1,Y1,Om2,X2,Y2) + +! 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 + +! endif end if ! Compute the RPA correlation energy - EcRPA = 0.5d0 * (sum(Om1) - sum(Om2) - trace_matrix(nVV, Cpp) - trace_matrix(nOO, Dpp)) - EcRPA1 = +sum(Om1) - trace_matrix(nVV, Cpp) - EcRPA2 = -sum(Om2) - trace_matrix(nOO, Dpp) + + EcRPA = 0.5d0 * (sum(Om1) - sum(Om2) - trace_matrix(nVV,Cpp) - trace_matrix(nOO,Dpp)) + EcRPA1 = +sum(Om1) - trace_matrix(nVV,Cpp) + EcRPA2 = -sum(Om2) - trace_matrix(nOO,Dpp) if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) then print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!' endif - deallocate(M, Z, Om) + deallocate(M,Z,Om) end subroutine - - diff --git a/src/LR/ppLR.f90 b/src/LR/ppLR.f90 index 6d62908..362fcd7 100644 --- a/src/LR/ppLR.f90 +++ b/src/LR/ppLR.f90 @@ -28,54 +28,59 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) implicit none include 'parameters.h' - logical, intent(in) :: TDA - integer, intent(in) :: nOO, nVV - double precision, intent(in) :: Bpp(nVV,nOO), Cpp(nVV,nVV), Dpp(nOO,nOO) - double precision, intent(out) :: Om1(nVV), X1(nVV,nVV), Y1(nOO,nVV) - double precision, intent(out) :: Om2(nOO), X2(nVV,nOO), Y2(nOO,nOO) - double precision, intent(out) :: EcRPA + logical,intent(in) :: TDA + integer,intent(in) :: nOO,nVV + double precision,intent(in) :: Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO) + double precision,intent(out) :: Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV) + double precision,intent(out) :: Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO) + double precision,intent(out) :: EcRPA - logical :: imp_bio, verbose - integer :: i, j, N - double precision :: EcRPA1, EcRPA2 - double precision :: thr_d, thr_nd, thr_deg - double precision,allocatable :: M(:,:), Z(:,:), Om(:) + logical :: imp_bio,verbose + integer :: i,j + integer :: Npp + double precision :: EcRPA1,EcRPA2 + double precision :: thr_d,thr_nd,thr_deg + double precision,allocatable :: M(:,:),Z(:,:),Om(:) - double precision, external :: trace_matrix + double precision,external :: trace_matrix - N = nOO + nVV + nPP = nOO + nVV - allocate(M(N,N), Z(N,N), Om(N)) + allocate(M(nPP,nPP),Z(nPP,nPP),Om(nPP)) + + ! Hermitian case for TDA if(TDA) then X1(:,:) = +Cpp(:,:) Y1(:,:) = 0d0 - if(nVV > 0) call diagonalize_matrix(nVV, X1, Om1) + if(nVV > 0) call diagonalize_matrix(nVV,X1,Om1) X2(:,:) = 0d0 Y2(:,:) = -Dpp(:,:) - if(nOO > 0) call diagonalize_matrix(nOO, Y2, Om2) + if(nOO > 0) call diagonalize_matrix(nOO,Y2,Om2) else ! Diagonal blocks - M( 1:nVV , 1:nVV) = + Cpp(1:nVV,1:nVV) - M(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - Dpp(1:nOO,1:nOO) + + M( 1:nVV, 1:nVV) = + Cpp(1:nVV,1:nVV) + M(nVV+1:nPP,nVV+1:nPP) = - Dpp(1:nOO,1:nOO) ! Off-diagonal blocks - 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( 1:nVV,nVV+1:nPP) = + Bpp(1:nVV,1:nOO) + M(nVV+1:nPP, 1:nVV) = - transpose(Bpp(1:nVV,1:nOO)) ! if((nOO .eq. 0) .or. (nVV .eq. 0)) then ! Diagonalize the pp matrix - if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Om,Z) + if(nPP > 0) call diagonalize_general_matrix(nPP,M,Om,Z) - ! Split the various quantities in p-p and h-h parts + ! Split the various quantities in ee and hh parts - call sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) + call sort_ppRPA(nOO,nVV,nPP,Om,Z,Om1,X1,Y1,Om2,X2,Y2) ! else @@ -112,9 +117,9 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) ! Compute the RPA correlation energy - EcRPA = 0.5d0 * (sum(Om1) - sum(Om2) - trace_matrix(nVV, Cpp) - trace_matrix(nOO, Dpp)) - EcRPA1 = +sum(Om1) - trace_matrix(nVV, Cpp) - EcRPA2 = -sum(Om2) - trace_matrix(nOO, Dpp) + EcRPA = 0.5d0 * (sum(Om1) - sum(Om2) - trace_matrix(nVV,Cpp) - trace_matrix(nOO,Dpp)) + EcRPA1 = +sum(Om1) - trace_matrix(nVV,Cpp) + EcRPA2 = -sum(Om2) - trace_matrix(nOO,Dpp) if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) then print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!' @@ -123,5 +128,3 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) deallocate(M,Z,Om) end subroutine - - diff --git a/src/LR/ppULR.f90 b/src/LR/ppULR.f90 index 49afb79..f9e40af 100644 --- a/src/LR/ppULR.f90 +++ b/src/LR/ppULR.f90 @@ -81,8 +81,7 @@ subroutine ppULR(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,nHaa,nHab,nHbb,nH ! Split the various quantities in p-p and h-h parts - call sort_ppRPA(nHt,nPt,Om(:),Z(:,:),Om1(:),X1(:,:),Y1(:,:),Om2(:),X2(:,:),& - Y2(:,:)) + call sort_ppRPA(nHt,nPt,nHt+nPt,Om,Z,Om1,X1,Y1,Om2,X2,Y2) ! Compute the RPA correlation energy diff --git a/src/RPA/sort_ppRPA.f90 b/src/RPA/sort_ppRPA.f90 index 24c6f59..1c1c80a 100644 --- a/src/RPA/sort_ppRPA.f90 +++ b/src/RPA/sort_ppRPA.f90 @@ -1,4 +1,4 @@ -subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) +subroutine sort_ppRPA(nOO,nVV,nPP,Om,Z,Om1,X1,Y1,Om2,X2,Y2) ! Compute the metric matrix for pp-RPA @@ -9,14 +9,13 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) integer,intent(in) :: nOO integer,intent(in) :: nVV - double precision,intent(in) :: Om(nOO+nVV) - double precision,intent(in) :: Z(nOO+nVV,nOO+nVV) + integer,intent(in) :: nPP + double precision,intent(in) :: Om(nPP) + double precision,intent(in) :: Z(nPP,nPP) ! Local variables integer :: pq,ab,ij -! integer :: deg1,ab_start,ab_end -! integer :: deg2,ij_start,ij_end double precision,allocatable :: M(:,:) double precision,allocatable :: Z1(:,:) double precision,allocatable :: Z2(:,:) @@ -42,7 +41,7 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) ! Memory allocation - allocate(M(nOO+nVV,nOO+nVV),Z1(nOO+nVV,nVV),Z2(nOO+nVV,nOO),order1(nVV),order2(nOO)) + allocate(M(nPP,nPP),Z1(nPP,nVV),Z2(nPP,nOO),order1(nVV),order2(nOO)) ! Initializatiom @@ -71,19 +70,19 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) ab = 0 ij = 0 - do pq=1,nOO+nVV + do pq=1,nPP if(Om(pq) > 0d0) then ab = ab + 1 Om1(ab) = Om(pq) - Z1(1:nOO+nVV,ab) = Z(1:nOO+nVV,pq) + Z1(1:nPP,ab) = Z(1:nPP,pq) else ij = ij + 1 Om2(ij) = Om(pq) - Z2(1:nOO+nVV,ij) = Z(1:nOO+nVV,pq) + Z2(1:nPP,ij) = Z(1:nPP,pq) end if @@ -99,7 +98,7 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) end do call quick_sort(Om1,order1,nVV) - call set_order(Z1,order1,nOO+nVV,nVV) + call set_order(Z1,order1,nPP,nVV) end if @@ -110,17 +109,17 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) end do call quick_sort(Om2,order2,nOO) - call set_order(Z2,order2,nOO+nVV,nOO) + call set_order(Z2,order2,nPP,nOO) end if allocate(S1(nVV,nVV),S2(nOO,nOO),O1(nVV,nVV),O2(nOO,nOO)) - allocate(tmp1(nOO+nVV,nVV),tmp2(nOO+nVV,nOO)) + allocate(tmp1(nPP,nVV),tmp2(nPP,nOO)) - if(nVV > 0) call dgemm ('N', 'N', nOO+nVV, nVV, nOO+nVV, 1d0, M, nOO+nVV, Z1, nOO+nVV, 0d0, tmp1, nOO+nVV) - if(nVV > 0) call dgemm ('T', 'N', nVV , nVV, nOO+nVV, 1d0, Z1, nOO+nVV, tmp1, nOO+nVV, 0d0, S1, nVV) - if(nOO > 0) call dgemm ('N', 'N', nOO+nVV, nOO, nOO+nVV, 1d0, M, nOO+nVV, -1d0*Z2, nOO+nVV, 0d0, tmp2, nOO+nVV) - if(nOO > 0) call dgemm ('T', 'N', nOO , nOO, nOO+nVV, 1d0, Z2, nOO+nVV, tmp2, nOO+nVV, 0d0, S2, nOO) + if(nVV > 0) call dgemm ('N','N',nPP,nVV,nPP,1d0, M,nPP,Z1,nPP,0d0,tmp1,nPP) + if(nVV > 0) call dgemm ('T','N',nVV,nVV,nPP,1d0,Z1,nPP,tmp1,nPP,0d0,S1,nVV) + if(nOO > 0) call dgemm ('N','N',nPP,nOO,nPP,1d0, M,nPP,-1d0*Z2,nPP,0d0,tmp2, nPP) + if(nOO > 0) call dgemm ('T','N',nOO,nOO,nPP,1d0,Z2,nPP,tmp2,nPP,0d0,S2,nOO) ! S1 = + matmul(transpose(Z1),matmul(M,Z1)) ! S2 = - matmul(transpose(Z2),matmul(M,Z2)) @@ -128,9 +127,9 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) if(nVV > 0) call orthogonalize_matrix(1,nVV,S1,O1) if(nOO > 0) call orthogonalize_matrix(1,nOO,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',nPP,nVV,nVV,1d0,Z1,nPP,O1,nVV,0d0,tmp1,nPP) Z1 = tmp1 - if(nOO > 0) call dgemm ('N', 'N', nOO+nVV,nOO,nOO, 1d0, Z2, nOO+nVV, O2, nOO,0d0, tmp2, nOO+nVV) + if(nOO > 0) call dgemm ('N','N',nPP,nOO,nOO,1d0,Z2,nPP,O2,nOO,0d0,tmp2,nPP) Z2 = tmp2 ! Z1 = matmul(Z1,O1) @@ -138,11 +137,11 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) ! Define submatrices X1, Y1, X2, & Y2 - X1(1:nVV,1:nVV) = + Z1( 1: nVV,1:nVV) - Y1(1:nOO,1:nVV) = - Z1(nVV+1:nOO+nVV,1:nVV) + X1(1:nVV,1:nVV) = Z1( 1: nVV,1:nVV) + Y1(1:nOO,1:nVV) = Z1(nVV+1:nPP,1:nVV) - X2(1:nVV,1:nOO) = + Z2( 1: nVV,1:nOO) - Y2(1:nOO,1:nOO) = - Z2(nVV+1:nOO+nVV,1:nOO) + X2(1:nVV,1:nOO) = Z2( 1: nVV,1:nOO) + Y2(1:nOO,1:nOO) = Z2(nVV+1:nPP,1:nOO) ! call matout(nVV,nVV,X1) ! call matout(nOO,nVV,Y1) @@ -150,4 +149,11 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) ! call matout(nVV,nOO,X2) ! call matout(nOO,nOO,Y2) +! Check orthonormality + +! call matout(nVV,nVV,matmul(transpose(X1),X1) - matmul(transpose(Y1),Y1)) +! call matout(nOO,nOO,matmul(transpose(X2),X2) - matmul(transpose(Y2),Y2)) +! call matout(nVV,nOO,matmul(transpose(X1),X2) - matmul(transpose(Y1),Y2)) +! call matout(nOO,nVV,matmul(transpose(X2),X1) - matmul(transpose(Y2),Y1)) + end subroutine