PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-31 16:45:49 +01:00
Code Releases Activity

TDA in pp-RPA and GT

Browse Source
This commit is contained in:
Pierre-Francois Loos 2021-10-18 22:05:26 +02:00
parent 6f11f0909f
commit 22aeccb222
9 changed files with 86 additions and 109 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
input/methods
View File

@ -13,9 +13,9 @@
# G0F2* evGF2* qsGF2* G0F3 evGF3
F F F F F
# G0W0* evGW* qsGW* ufG0W0 ufGW
T F F F F
F F F F F
# G0T0 evGT qsGT
F F T
T F F
# MCMP2
F
# * unrestricted version available

12
src/LR/linear_response.f90
View File

@ -7,9 +7,17 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E
! Input variables
logical,intent(in) :: dRPA,TDA,BSE
logical,intent(in) :: dRPA
logical,intent(in) :: TDA
logical,intent(in) :: BSE
double precision,intent(in) :: eta
integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS
integer,intent(in) :: ispin
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: lambda
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)

126
src/LR/linear_response_pp.f90
View File

@ -1,4 +1,4 @@
subroutine linear_response_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV, &
subroutine linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV, &
e,ERI,Omega1,X1,Y1,Omega2,X2,Y2,EcRPA)
! Compute the p-p channel of the linear response: see Scuseria et al. JCP 139, 104113 (2013)
@ -8,7 +8,13 @@ subroutine linear_response_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV, &
! Input variables
integer,intent(in) :: ispin,nBas,nC,nO,nV,nR
integer,intent(in) :: ispin
logical,intent(in) :: TDA
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nOO
integer,intent(in) :: nVV
double precision,intent(in) :: e(nBas)
@ -16,7 +22,6 @@ subroutine linear_response_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV, &
! Local variables
logical :: dump_matrices = .false.
integer :: ab,cd,ij,kl
integer :: p,q,r,s
double precision :: trace_matrix
@ -43,92 +48,55 @@ subroutine linear_response_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV, &
allocate(B(nVV,nOO),C(nVV,nVV),D(nOO,nOO),M(nOO+nVV,nOO+nVV),Z(nOO+nVV,nOO+nVV),Omega(nOO+nVV))
!-------------------------------------------------!
! Solve the p-p eigenproblem !
!-------------------------------------------------!
! !
! | C -B | | X1 X2 | | w1 0 | | X1 X2 | !
! | | | | = | | | | !
! | Bt -D | | Y1 Y2 | | 0 w2 | | Y1 Y2 | !
! !
!-------------------------------------------------!
! Build B, C and D matrices for the pp channel
call linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B)
call linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C)
call linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D)
!------------------------------------------------------------------------
! Solve the p-p eigenproblem
!------------------------------------------------------------------------
!
! | C -B | | X1 X2 | | w1 0 | | X1 X2 |
! | | | | = | | | |
! | Bt -D | | Y1 Y2 | | 0 w2 | | Y1 Y2 |
!
if(TDA) then
! Diagonal blocks
M( 1:nVV , 1:nVV) = + C(1:nVV,1:nVV)
M(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - D(1:nOO,1:nOO)
! Off-diagonal blocks
M( 1:nVV ,nVV+1:nOO+nVV) = - B(1:nVV,1:nOO)
M(nVV+1:nOO+nVV, 1:nVV) = + transpose(B(1:nVV,1:nOO))
! Dump ppRPA matrices
if(dump_matrices) then
open(unit=42,file='B.dat')
open(unit=43,file='C.dat')
open(unit=44,file='D.dat')
open(unit=45,file='ERI.dat')
open(unit=46,file='eps.dat')
do ab=1,nVV
do ij=1,nOO
if(abs(B(ab,ij)) > 1d-15) write(42,*) ab,ij,B(ab,ij)
end do
end do
do ab=1,nVV
do cd=1,nVV
if(abs(C(ab,cd)) > 1d-15) write(43,*) ab,cd,C(ab,cd)
end do
end do
do ij=1,nOO
do kl=1,nOO
if(abs(D(ij,kl)) > 1d-15) write(44,*) ij,kl,D(ij,kl)
end do
end do
do p=1,nBas
write(46,*) p,e(p)
do q=1,nBas
do r=1,nBas
do s=1,nBas
if(abs(ERI(p,q,r,s)) > 1d-15) write(45,*) p,q,r,s,ERI(p,q,r,s)
end do
end do
end do
end do
X1(:,:) = +C(:,:)
Y1(:,:) = 0d0
call diagonalize_matrix(nVV,X1,Omega1)
close(42)
close(43)
close(44)
close(45)
close(46)
X2(:,:) = 0d0
Y2(:,:) = -D(:,:)
call diagonalize_matrix(nOO,Y2,Omega2)
else
call linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B)
! Diagonal blocks
M( 1:nVV , 1:nVV) = + C(1:nVV,1:nVV)
M(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - D(1:nOO,1:nOO)
! Off-diagonal blocks
M( 1:nVV ,nVV+1:nOO+nVV) = - B(1:nVV,1:nOO)
M(nVV+1:nOO+nVV, 1:nVV) = + transpose(B(1:nVV,1:nOO))
! Diagonalize the p-h matrix
if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Omega,Z)
! Split the various quantities in p-p and h-h parts
call sort_ppRPA(nOO,nVV,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),Y2(:,:))
end if
! Diagonalize the p-h matrix
if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Omega,Z)
! allocate(order(nOO+nVV))
! call quick_sort(Omega(:),order(:),nOO+nVV)
! call matout(nOO+nVV,1,Omega(:)*HaToeV)
! Split the various quantities in p-p and h-h parts
call sort_ppRPA(nOO,nVV,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),Y2(:,:))
! call matout(32,1,(Omega1(:) - Omega1(1))*HaToeV)
! Compute the RPA correlation energy
EcRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) )

4
src/MBPT/Bethe_Salpeter_Tmatrix.f90
View File

@ -83,7 +83,7 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
ispin = 1
iblock = 3
call linear_response_pp(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,eT,ERI, &
call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,eT,ERI, &
Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
@ -98,7 +98,7 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
ispin = 2
iblock = 4
call linear_response_pp(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,eT,ERI, &
call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,eT,ERI, &
Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)

8
src/MBPT/G0T0.f90
View File

@ -103,7 +103,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
! Compute linear response
call linear_response_pp(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,ERI_MO, &
call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,ERI_MO, &
Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
! EcRPA(ispin) = 1d0*EcRPA(ispin)
@ -120,7 +120,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
! Compute linear response
call linear_response_pp(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,ERI_MO, &
call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,ERI_MO, &
Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
! EcRPA(ispin) = 2d0*EcRPA(ispin)
@ -183,11 +183,11 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
ispin = 1
iblock = 3
call linear_response_pp(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,eG0T0,ERI_MO, &
call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,eG0T0,ERI_MO, &
Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
ispin = 2
iblock = 4
call linear_response_pp(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,eG0T0,ERI_MO, &
call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,eG0T0,ERI_MO, &
Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
EcRPA(1) = EcRPA(1) - EcRPA(2)
EcRPA(2) = 3d0*EcRPA(2)

8
src/MBPT/evGT.f90
View File

@ -132,7 +132,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
! Compute linear response
call linear_response_pp(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,ERI_MO, &
call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,ERI_MO, &
Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
!----------------------------------------------
@ -144,7 +144,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
! Compute linear response
call linear_response_pp(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,ERI_MO, &
call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,ERI_MO, &
Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
!----------------------------------------------
@ -222,11 +222,11 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
ispin = 1
iblock = 3
call linear_response_pp(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,ERI_MO, &
call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,ERI_MO, &
Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
ispin = 2
iblock = 4
call linear_response_pp(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,ERI_MO, &
call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,ERI_MO, &
Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
EcRPA(1) = EcRPA(1) - EcRPA(2)
EcRPA(2) = 3d0*EcRPA(2)

8
src/MBPT/qsGT.f90
View File

@ -188,13 +188,13 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T
ispin = 1
iblock = 3
call linear_response_pp(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,ERI_MO, &
call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,ERI_MO, &
Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
ispin = 2
iblock = 4
call linear_response_pp(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,ERI_MO, &
call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,ERI_MO, &
Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
! Compute correlation part of the self-energy
@ -305,11 +305,11 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T
ispin = 1
iblock = 3
call linear_response_pp(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,ERI_MO, &
call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,ERI_MO, &
Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
ispin = 2
iblock = 4
call linear_response_pp(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,ERI_MO, &
call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,ERI_MO, &
Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
EcRPA(1) = EcRPA(1) - EcRPA(2)
EcRPA(2) = 3d0*EcRPA(2)

2
src/QuAcK/QuAcK.f90
View File

@ -792,7 +792,7 @@ program QuAcK
if(doppRPA) then
call cpu_time(start_ppRPA)
call ppRPA(singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF)
call ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_ppRPA)
t_ppRPA = end_ppRPA - start_ppRPA

23
src/RPA/ppRPA.f90
View File

@ -1,4 +1,4 @@
subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e)
subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e)
! Perform pp-RPA calculation
@ -7,8 +7,9 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER
! Input variables
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
logical,intent(in) :: TDA
logical,intent(in) :: singlet
logical,intent(in) :: triplet
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
@ -47,7 +48,7 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER
! Singlet manifold
if(singlet_manifold) then
if(singlet) then
ispin = 1
@ -61,9 +62,9 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER
allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), &
Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin))
call linear_response_pp(ispin,.false.,.false.,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, &
Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), &
Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), &
call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, &
Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), &
Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), &
Ec_ppRPA(ispin))
call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin))
@ -75,7 +76,7 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER
! Triplet manifold
if(triplet_manifold) then
if(triplet) then
ispin = 2
@ -90,9 +91,9 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER
Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin))
call linear_response_pp(ispin,.false.,.false.,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, &
Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), &
Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), &
call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, &
Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), &
Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), &
Ec_ppRPA(ispin))
call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin))