From 22aeccb222a5ca2d90587316a025a7f2a4317b28 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 18 Oct 2021 22:05:26 +0200 Subject: [PATCH] TDA in pp-RPA and GT --- input/methods | 4 +- src/LR/linear_response.f90 | 12 ++- src/LR/linear_response_pp.f90 | 126 +++++++++++----------------- src/MBPT/Bethe_Salpeter_Tmatrix.f90 | 4 +- src/MBPT/G0T0.f90 | 8 +- src/MBPT/evGT.f90 | 8 +- src/MBPT/qsGT.f90 | 8 +- src/QuAcK/QuAcK.f90 | 2 +- src/RPA/ppRPA.f90 | 23 ++--- 9 files changed, 86 insertions(+), 109 deletions(-) diff --git a/input/methods b/input/methods index ab7063e..ddcffcb 100644 --- a/input/methods +++ b/input/methods @@ -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 diff --git a/src/LR/linear_response.f90 b/src/LR/linear_response.f90 index 831be12..ee80c72 100644 --- a/src/LR/linear_response.f90 +++ b/src/LR/linear_response.f90 @@ -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) diff --git a/src/LR/linear_response_pp.f90 b/src/LR/linear_response_pp.f90 index 399ee7b..4b35ec0 100644 --- a/src/LR/linear_response_pp.f90 +++ b/src/LR/linear_response_pp.f90 @@ -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(:,:)) ) diff --git a/src/MBPT/Bethe_Salpeter_Tmatrix.f90 b/src/MBPT/Bethe_Salpeter_Tmatrix.f90 index 63cdf1b..f16a703 100644 --- a/src/MBPT/Bethe_Salpeter_Tmatrix.f90 +++ b/src/MBPT/Bethe_Salpeter_Tmatrix.f90 @@ -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) diff --git a/src/MBPT/G0T0.f90 b/src/MBPT/G0T0.f90 index 75e52f6..125cde2 100644 --- a/src/MBPT/G0T0.f90 +++ b/src/MBPT/G0T0.f90 @@ -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) diff --git a/src/MBPT/evGT.f90 b/src/MBPT/evGT.f90 index 5c2e2e3..e657bf6 100644 --- a/src/MBPT/evGT.f90 +++ b/src/MBPT/evGT.f90 @@ -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) diff --git a/src/MBPT/qsGT.f90 b/src/MBPT/qsGT.f90 index 1f61fb6..59c43c7 100644 --- a/src/MBPT/qsGT.f90 +++ b/src/MBPT/qsGT.f90 @@ -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) diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 931b21f..946770a 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -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 diff --git a/src/RPA/ppRPA.f90 b/src/RPA/ppRPA.f90 index 3c45cc5..e5f1ea6 100644 --- a/src/RPA/ppRPA.f90 +++ b/src/RPA/ppRPA.f90 @@ -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))