From 608a62dd6ada17dc91c90011ce00c3eb8830b8c9 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 17 Jul 2023 23:05:49 +0200 Subject: [PATCH] remove construction of matrices in ppRPA --- src/GT/G0T0pp.f90 | 44 ++++++++++++++++++++++++++++++++++++++++---- src/GT/evGTpp.f90 | 23 +++++++++++++++++++++-- src/GT/qsGTpp.f90 | 25 ++++++++++++++++++++++--- src/LR/ppLR.f90 | 4 +--- 4 files changed, 84 insertions(+), 12 deletions(-) diff --git a/src/GT/G0T0pp.f90 b/src/GT/G0T0pp.f90 index a57ccb0..c5c149d 100644 --- a/src/GT/G0T0pp.f90 +++ b/src/GT/G0T0pp.f90 @@ -51,6 +51,9 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp double precision :: EcAC(nspin) double precision :: EcppBSE(nspin) double precision :: EcGM + double precision,allocatable :: Bpp(:,:) + double precision,allocatable :: Cpp(:,:) + double precision,allocatable :: Dpp(:,:) double precision,allocatable :: Om1ab(:),Om1aa(:) double precision,allocatable :: X1ab(:,:),X1aa(:,:) double precision,allocatable :: Y1ab(:,:),Y1aa(:,:) @@ -102,7 +105,15 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp ! Compute linear response - call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,eHF,ERI_MO,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin)) + allocate(Bpp(nVVab,nOOab),Cpp(nVVab,nVVab),Dpp(nOOab,nOOab)) + + if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,ERI_MO,Bpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,eHF,ERI_MO,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,eHF,ERI_MO,Dpp) + + call ppLR(TDA_T,nOOab,nVVab,Bpp,Cpp,Dpp,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin)) + + deallocate(Bpp,Cpp,Dpp) call print_excitation('pp-RPA (N+2)',iblock,nVVab,Om1ab(:)) call print_excitation('pp-RPA (N-2)',iblock,nOOab,Om2ab(:)) @@ -116,7 +127,15 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp ! Compute linear response - call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,eHF,ERI_MO,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin)) + allocate(Bpp(nVVaa,nOOaa),Cpp(nVVaa,nVVaa),Dpp(nOOaa,nOOaa)) + + if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,ERI_MO,Bpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,eHF,ERI_MO,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,eHF,ERI_MO,Dpp) + + call ppLR(TDA_T,nOOaa,nVVaa,Bpp,Cpp,Dpp,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin)) + + deallocate(Bpp,Cpp,Dpp) call print_excitation('pp-RPA (N+2)',iblock,nVVaa,Om1aa(:)) call print_excitation('pp-RPA (N-2)',iblock,nOOaa,Om2aa(:)) @@ -198,11 +217,28 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp ispin = 1 iblock = 3 - call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,eGT,ERI_MO,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin)) + allocate(Bpp(nVVab,nOOab),Cpp(nVVab,nVVab),Dpp(nOOab,nOOab)) + + if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,ERI_MO,Bpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,eGT,ERI_MO,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,eGT,ERI_MO,Dpp) + + call ppLR(TDA_T,nOOab,nVVab,Bpp,Cpp,Dpp,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin)) + + deallocate(Bpp,Cpp,Dpp) + ispin = 2 iblock = 4 - call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,eGT,ERI_MO,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin)) + allocate(Bpp(nVVaa,nOOaa),Cpp(nVVaa,nVVaa),Dpp(nOOaa,nOOaa)) + + if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,ERI_MO,Bpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,eGT,ERI_MO,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,eGT,ERI_MO,Dpp) + + call ppLR(TDA_T,nOOaa,nVVaa,Bpp,Cpp,Dpp,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin)) + + deallocate(Bpp,Cpp,Dpp) EcRPA(1) = EcRPA(1) - EcRPA(2) EcRPA(2) = 3d0*EcRPA(2) diff --git a/src/GT/evGTpp.f90 b/src/GT/evGTpp.f90 index e05e59e..f821c3e 100644 --- a/src/GT/evGTpp.f90 +++ b/src/GT/evGTpp.f90 @@ -60,6 +60,9 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T double precision,allocatable :: e_diis(:,:) double precision,allocatable :: eGT(:) double precision,allocatable :: eOld(:) + double precision,allocatable :: Bpp(:,:) + double precision,allocatable :: Cpp(:,:) + double precision,allocatable :: Dpp(:,:) double precision,allocatable :: Om1s(:),Om1t(:) double precision,allocatable :: X1s(:,:),X1t(:,:) double precision,allocatable :: Y1s(:,:),Y1t(:,:) @@ -132,7 +135,15 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T ! Compute linear response - call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,ERI_MO,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) + allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) + + if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eHF,ERI_MO,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eHF,ERI_MO,Dpp) + + call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) + + deallocate(Bpp,Cpp,Dpp) !---------------------------------------------- ! alpha-alpha block @@ -143,7 +154,15 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T ! Compute linear response - call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) + allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) + + if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eHF,ERI_MO,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eHF,ERI_MO,Dpp) + + call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) + + deallocate(Bpp,Cpp,Dpp) EcRPA(1) = EcRPA(1) - EcRPA(2) EcRPA(2) = 3d0*EcRPA(2) diff --git a/src/GT/qsGTpp.f90 b/src/GT/qsGTpp.f90 index 776e8a8..0f8da08 100644 --- a/src/GT/qsGTpp.f90 +++ b/src/GT/qsGTpp.f90 @@ -76,6 +76,9 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T double precision,allocatable :: cp(:,:) double precision,allocatable :: eGT(:) double precision,allocatable :: eOld(:) + double precision,allocatable :: Bpp(:,:) + double precision,allocatable :: Cpp(:,:) + double precision,allocatable :: Dpp(:,:) double precision,allocatable :: Om1s(:),Om1t(:) double precision,allocatable :: X1s(:,:),X1t(:,:) double precision,allocatable :: Y1s(:,:),Y1t(:,:) @@ -188,12 +191,28 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T ispin = 1 iblock = 3 - call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,ERI_MO,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) + allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) + + if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,ERI_MO,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,ERI_MO,Dpp) + + call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) + + deallocate(Bpp,Cpp,Dpp) ispin = 2 iblock = 4 - call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) + allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) + + if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO,Dpp) + + call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) + + deallocate(Bpp,Cpp,Dpp) EcRPA(1) = EcRPA(1) - EcRPA(2) EcRPA(2) = 3d0*EcRPA(2) @@ -201,7 +220,7 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T ! Compute correlation part of the self-energy EcGM = 0d0 - Sig(:,:) = 0d0 + Sig(:,:) = 0d0 Z(:) = 0d0 iblock = 3 diff --git a/src/LR/ppLR.f90 b/src/LR/ppLR.f90 index 64abe6c..e63dbda 100644 --- a/src/LR/ppLR.f90 +++ b/src/LR/ppLR.f90 @@ -1,4 +1,4 @@ -subroutine ppLR(TDA,nBas,nOO,nVV,Om1,X1,Y1,Om2,X2,Y2,EcRPA) +subroutine ppLR(TDA,nOO,nVV,B,C,D,Om1,X1,Y1,Om2,X2,Y2,EcRPA) ! Compute the p-p channel of the linear response: see Scuseria et al. JCP 139, 104113 (2013) @@ -61,8 +61,6 @@ subroutine ppLR(TDA,nBas,nOO,nVV,Om1,X1,Y1,Om2,X2,Y2,EcRPA) else - call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B) - ! Diagonal blocks M( 1:nVV , 1:nVV) = + C(1:nVV,1:nVV)