PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:34:46 +01:00
Code Releases Activity

redefine ppLR matrix

Browse Source
This commit is contained in:
Pierre-Francois Loos 2024-09-13 09:33:43 +02:00
parent ac03c0e109
commit 77fea9a9bd
4 changed files with 126 additions and 113 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

123
src/LR/ppGLR.f90
View File

@ -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 ! 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 implicit none
include 'parameters.h' include 'parameters.h'
logical, intent(in) :: TDA logical,intent(in) :: TDA
integer, intent(in) :: nOO, nVV integer,intent(in) :: nOO,nVV
double precision, intent(in) :: Bpp(nVV,nOO), Cpp(nVV,nVV), Dpp(nOO,nOO) 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) :: Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV)
double precision, intent(out) :: Om2(nOO), X2(nVV,nOO), Y2(nOO,nOO) double precision,intent(out) :: Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO)
double precision, intent(out) :: EcRPA double precision,intent(out) :: EcRPA
logical :: imp_bio, verbose logical :: imp_bio,verbose
integer :: i, j, N integer :: i,j
double precision :: EcRPA1, EcRPA2 integer :: Npp
double precision :: thr_d, thr_nd, thr_deg double precision :: EcRPA1,EcRPA2
double precision,allocatable :: M(:,:), Z(:,:), Om(:) 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 ! Hermitian case for TDA
allocate(M(N,N), Z(N,N), Om(N))
if(TDA) then if(TDA) then
X1(:,:) = +Cpp(:,:) X1(:,:) = +Cpp(:,:)
Y1(:,:) = 0d0 Y1(:,:) = 0d0
if(nVV > 0) call diagonalize_matrix(nVV, X1, Om1) if(nVV > 0) call diagonalize_matrix(nVV,X1,Om1)
X2(:,:) = 0d0 X2(:,:) = 0d0
Y2(:,:) = -Dpp(:,:) Y2(:,:) = -Dpp(:,:)
if(nOO > 0) call diagonalize_matrix(nOO, Y2, Om2) if(nOO > 0) call diagonalize_matrix(nOO,Y2,Om2)
else else
! Diagonal blocks ! 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 ! 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 .eq. 0) .or. (nVV .eq. 0)) then
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 ! Diagonalize the pp matrix
thr_d = 1d-6 ! to determine if diagonal elements of L.T x R are close enouph to 1 if(nPP > 0) call diagonalize_general_matrix(nPP,M,Om,Z)
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 ! Split the various quantities in ee and hh parts
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 call sort_ppRPA(nOO,nVV,nPP,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
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 ! 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 end if
! Compute the RPA correlation energy ! 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) EcRPA = 0.5d0 * (sum(Om1) - sum(Om2) - trace_matrix(nVV,Cpp) - trace_matrix(nOO,Dpp))
EcRPA2 = -sum(Om2) - 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 if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) then
print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!' print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!'
endif endif
deallocate(M, Z, Om) deallocate(M,Z,Om)
end subroutine end subroutine

59
src/LR/ppLR.f90
View File

@ -28,54 +28,59 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
implicit none implicit none
include 'parameters.h' include 'parameters.h'
logical, intent(in) :: TDA logical,intent(in) :: TDA
integer, intent(in) :: nOO, nVV integer,intent(in) :: nOO,nVV
double precision, intent(in) :: Bpp(nVV,nOO), Cpp(nVV,nVV), Dpp(nOO,nOO) 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) :: Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV)
double precision, intent(out) :: Om2(nOO), X2(nVV,nOO), Y2(nOO,nOO) double precision,intent(out) :: Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO)
double precision, intent(out) :: EcRPA double precision,intent(out) :: EcRPA
logical :: imp_bio, verbose logical :: imp_bio,verbose
integer :: i, j, N integer :: i,j
double precision :: EcRPA1, EcRPA2 integer :: Npp
double precision :: thr_d, thr_nd, thr_deg double precision :: EcRPA1,EcRPA2
double precision,allocatable :: M(:,:), Z(:,:), Om(:) 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 if(TDA) then
X1(:,:) = +Cpp(:,:) X1(:,:) = +Cpp(:,:)
Y1(:,:) = 0d0 Y1(:,:) = 0d0
if(nVV > 0) call diagonalize_matrix(nVV, X1, Om1) if(nVV > 0) call diagonalize_matrix(nVV,X1,Om1)
X2(:,:) = 0d0 X2(:,:) = 0d0
Y2(:,:) = -Dpp(:,:) Y2(:,:) = -Dpp(:,:)
if(nOO > 0) call diagonalize_matrix(nOO, Y2, Om2) if(nOO > 0) call diagonalize_matrix(nOO,Y2,Om2)
else else
! Diagonal blocks ! 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 ! 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 ! if((nOO .eq. 0) .or. (nVV .eq. 0)) then
! Diagonalize the pp matrix ! 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 ! 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 ! Compute the RPA correlation energy
EcRPA = 0.5d0 * (sum(Om1) - sum(Om2) - trace_matrix(nVV, Cpp) - 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) EcRPA1 = +sum(Om1) - trace_matrix(nVV,Cpp)
EcRPA2 = -sum(Om2) - trace_matrix(nOO, Dpp) EcRPA2 = -sum(Om2) - trace_matrix(nOO,Dpp)
if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) then if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) then
print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!' 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) deallocate(M,Z,Om)
end subroutine end subroutine

3
src/LR/ppULR.f90
View File

@ -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 ! Split the various quantities in p-p and h-h parts
call sort_ppRPA(nHt,nPt,Om(:),Z(:,:),Om1(:),X1(:,:),Y1(:,:),Om2(:),X2(:,:),& call sort_ppRPA(nHt,nPt,nHt+nPt,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
Y2(:,:))
! Compute the RPA correlation energy ! Compute the RPA correlation energy

50
src/RPA/sort_ppRPA.f90
View File

@ -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 ! 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) :: nOO
integer,intent(in) :: nVV integer,intent(in) :: nVV
double precision,intent(in) :: Om(nOO+nVV) integer,intent(in) :: nPP
double precision,intent(in) :: Z(nOO+nVV,nOO+nVV) double precision,intent(in) :: Om(nPP)
double precision,intent(in) :: Z(nPP,nPP)
! Local variables ! Local variables
integer :: pq,ab,ij integer :: pq,ab,ij
! integer :: deg1,ab_start,ab_end
! integer :: deg2,ij_start,ij_end
double precision,allocatable :: M(:,:) double precision,allocatable :: M(:,:)
double precision,allocatable :: Z1(:,:) double precision,allocatable :: Z1(:,:)
double precision,allocatable :: Z2(:,:) double precision,allocatable :: Z2(:,:)
@ -42,7 +41,7 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
! Memory allocation ! 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 ! Initializatiom
@ -71,19 +70,19 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
ab = 0 ab = 0
ij = 0 ij = 0
do pq=1,nOO+nVV do pq=1,nPP
if(Om(pq) > 0d0) then if(Om(pq) > 0d0) then
ab = ab + 1 ab = ab + 1
Om1(ab) = Om(pq) Om1(ab) = Om(pq)
Z1(1:nOO+nVV,ab) = Z(1:nOO+nVV,pq) Z1(1:nPP,ab) = Z(1:nPP,pq)
else else
ij = ij + 1 ij = ij + 1
Om2(ij) = Om(pq) Om2(ij) = Om(pq)
Z2(1:nOO+nVV,ij) = Z(1:nOO+nVV,pq) Z2(1:nPP,ij) = Z(1:nPP,pq)
end if end if
@ -99,7 +98,7 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
end do end do
call quick_sort(Om1,order1,nVV) call quick_sort(Om1,order1,nVV)
call set_order(Z1,order1,nOO+nVV,nVV) call set_order(Z1,order1,nPP,nVV)
end if end if
@ -110,17 +109,17 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
end do end do
call quick_sort(Om2,order2,nOO) call quick_sort(Om2,order2,nOO)
call set_order(Z2,order2,nOO+nVV,nOO) call set_order(Z2,order2,nPP,nOO)
end if end if
allocate(S1(nVV,nVV),S2(nOO,nOO),O1(nVV,nVV),O2(nOO,nOO)) 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 ('N','N',nPP,nVV,nPP,1d0, M,nPP,Z1,nPP,0d0,tmp1,nPP)
if(nVV > 0) call dgemm ('T', 'N', nVV , nVV, nOO+nVV, 1d0, Z1, nOO+nVV, tmp1, nOO+nVV, 0d0, S1, nVV) 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', nOO+nVV, nOO, nOO+nVV, 1d0, M, nOO+nVV, -1d0*Z2, nOO+nVV, 0d0, tmp2, nOO+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, nOO+nVV, 1d0, Z2, nOO+nVV, tmp2, nOO+nVV, 0d0, S2, nOO) 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)) ! S1 = + matmul(transpose(Z1),matmul(M,Z1))
! S2 = - matmul(transpose(Z2),matmul(M,Z2)) ! 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(nVV > 0) call orthogonalize_matrix(1,nVV,S1,O1)
if(nOO > 0) call orthogonalize_matrix(1,nOO,S2,O2) 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 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 Z2 = tmp2
! Z1 = matmul(Z1,O1) ! 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 ! Define submatrices X1, Y1, X2, & Y2
X1(1:nVV,1:nVV) = + Z1( 1: nVV,1:nVV) X1(1:nVV,1:nVV) = Z1( 1: nVV,1:nVV)
Y1(1:nOO,1:nVV) = - Z1(nVV+1:nOO+nVV,1:nVV) Y1(1:nOO,1:nVV) = Z1(nVV+1:nPP,1:nVV)
X2(1:nVV,1:nOO) = + Z2( 1: nVV,1:nOO) X2(1:nVV,1:nOO) = Z2( 1: nVV,1:nOO)
Y2(1:nOO,1:nOO) = - Z2(nVV+1:nOO+nVV,1:nOO) Y2(1:nOO,1:nOO) = Z2(nVV+1:nPP,1:nOO)
! call matout(nVV,nVV,X1) ! call matout(nVV,nVV,X1)
! call matout(nOO,nVV,Y1) ! 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(nVV,nOO,X2)
! call matout(nOO,nOO,Y2) ! 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 end subroutine