From 36685fcc734d3c2f24b7ed377c23cbf0bc865eba Mon Sep 17 00:00:00 2001 From: Antoine MARIE Date: Sat, 22 Mar 2025 17:20:56 +0100 Subject: [PATCH] starting work on G Parquet --- src/AOtoMO/AOtoMO_ERI_RHF.f90 | 3 +- src/GW/RG0W0.f90 | 5 +- src/GW/RGW.f90 | 7 +- src/LR/print_excitation_energies.f90 | 2 +- src/Parquet/GParquet.f90 | 281 +++++++++++++++++++++++++++ src/Parquet/RParquet.f90 | 13 +- src/Parquet/R_screened_integrals.f90 | 1 + src/QuAcK/GQuAcK.f90 | 12 +- src/QuAcK/QuAcK.f90 | 6 +- src/QuAcK/RQuAcK.f90 | 5 +- 10 files changed, 314 insertions(+), 21 deletions(-) create mode 100644 src/Parquet/GParquet.f90 diff --git a/src/AOtoMO/AOtoMO_ERI_RHF.f90 b/src/AOtoMO/AOtoMO_ERI_RHF.f90 index a248ce1..ac57590 100644 --- a/src/AOtoMO/AOtoMO_ERI_RHF.f90 +++ b/src/AOtoMO/AOtoMO_ERI_RHF.f90 @@ -32,6 +32,7 @@ subroutine AOtoMO_ERI_RHF(nBas,nOrb,c,ERI_AO,ERI_MO) , ERI_AO(1,1,1,1), nBas, c(1,1), nBas & , 0.d0, a2(1,1,1,1), nBas*nBas*nBas) + call dgemm( 'T', 'N', nBas*nBas*nOrb, nOrb, nBas, 1.d0 & , a2(1,1,1,1), nBas, c(1,1), nBas & , 0.d0, a1(1,1,1,1), nBas*nBas*nOrb) @@ -50,5 +51,5 @@ subroutine AOtoMO_ERI_RHF(nBas,nOrb,c,ERI_AO,ERI_MO) , 0.d0, ERI_MO(1,1,1,1), nOrb*nOrb*nOrb) deallocate(a2) - + end subroutine diff --git a/src/GW/RG0W0.f90 b/src/GW/RG0W0.f90 index 5362af0..1be8512 100644 --- a/src/GW/RG0W0.f90 +++ b/src/GW/RG0W0.f90 @@ -1,5 +1,5 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & - linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) + linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF,eGW_out) ! Perform G0W0 calculation @@ -61,6 +61,7 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA double precision,allocatable :: eGWlin(:) double precision,allocatable :: eGW(:) + double precision,intent(inout):: eGW_out(nOrb) ! Output variables @@ -171,6 +172,8 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA call print_RG0W0(nOrb,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM) + eGW_out(:) = eGW(:) + !---------------------------! ! Perform phBSE calculation ! !---------------------------! diff --git a/src/GW/RGW.f90 b/src/GW/RGW.f90 index aac941d..23eadc6 100644 --- a/src/GW/RGW.f90 +++ b/src/GW/RGW.f90 @@ -1,7 +1,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_diis,doACFDT, & exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, & linearize,eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF, & - S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF,eGW) ! Restricted GW module @@ -70,6 +70,9 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_dii logical :: doccG0W0,doccGW + + double precision,intent(inout) :: eGW(nOrb) + !------------------------------------------------------------------------ ! Perform G0W0 calculation !------------------------------------------------------------------------ @@ -78,7 +81,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_dii call wall_time(start_GW) call RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & - linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) + linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF,eGW) call wall_time(end_GW) t_GW = end_GW - start_GW diff --git a/src/LR/print_excitation_energies.f90 b/src/LR/print_excitation_energies.f90 index 8d1448c..b909f21 100644 --- a/src/LR/print_excitation_energies.f90 +++ b/src/LR/print_excitation_energies.f90 @@ -14,7 +14,7 @@ subroutine print_excitation_energies(method,manifold,nS,Om) ! Local variables - integer,parameter :: maxS = 10 + integer,parameter :: maxS = 25 integer :: m write(*,*) diff --git a/src/Parquet/GParquet.f90 b/src/Parquet/GParquet.f90 new file mode 100644 index 0000000..c49f8b7 --- /dev/null +++ b/src/Parquet/GParquet.f90 @@ -0,0 +1,281 @@ +subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,ERI) + +! Parquet approximation based on restricted orbitals + + implicit none + include 'parameters.h' + +! Hard-coded parameters + + logical :: linearize = .true. + logical :: TDA = .true. + logical :: print_phLR = .false. + logical :: print_ppLR = .false. + +! Input variables + + integer,intent(in) :: max_it_1b,max_it_2b + double precision,intent(in) :: conv_1b,conv_2b + integer,intent(in) :: nOrb,nC,nO,nV,nR,nS + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + +! Local variables + + integer :: ispin + + integer :: n_it_1b,n_it_2b + double precision :: err_1b,err_2b + double precision :: err_eh, err_hh, err_ee + double precision :: start_t, end_t, t + double precision :: start_1b, end_1b, t_1b + double precision :: start_2b, end_2b, t_2b + + integer :: nOO,nVV + + double precision :: EcRPA + double precision,allocatable :: Aph(:,:), Bph(:,:) + double precision,allocatable :: XpY(:,:),XmY(:,:) + double precision,allocatable :: eh_Om(:), old_eh_Om(:) + + double precision,allocatable :: Bpp(:,:), Cpp(:,:), Dpp(:,:) + double precision,allocatable :: X1(:,:),Y1(:,:) + double precision,allocatable :: ee_Om(:), old_ee_Om(:) + double precision,allocatable :: X2(:,:),Y2(:,:) + double precision,allocatable :: hh_Om(:), old_hh_Om(:) + + double precision,allocatable :: eh_rho(:,:,:), ee_rho(:,:,:), hh_rho(:,:,:) + + double precision,allocatable :: eh_Gam_A(:,:),eh_Gam_B(:,:) + double precision,allocatable :: pp_Gam_B(:,:),pp_Gam_C(:,:),pp_Gam_D(:,:) + double precision,allocatable :: eh_Gam(:,:,:,:),pp_Gam(:,:,:,:) + + double precision,allocatable :: eQPlin(:),eQP(:),eOld(:) + double precision,allocatable :: SigC(:) + double precision,allocatable :: Z(:) + double precision :: EcGM + +! Output variables +! None + + nOO = nO*(nO - 1)/2 + nVV = nV*(nV - 1)/2 + + allocate(eQP(nOrb),eOld(nOrb)) + + write(*,*) + write(*,*)'***********************************' + write(*,*)'* Generalized Parquet Calculation *' + write(*,*)'***********************************' + write(*,*) + + ! Print parameters + + write(*,*)'---------------------------------------------------------------' + write(*,*)' Parquet parameters for one-body and two-body self-consistency ' + write(*,*)'---------------------------------------------------------------' + write(*,'(1X,A50,1X,I5)') 'Maximum number for one-body self-consistency:', max_it_1b + write(*,'(1X,A50,1X,E10.5)') 'Convergence threshold for one-body energies:', conv_1b + write(*,*)'---------------------------------------------------------------' + write(*,'(1X,A50,1X,I5)') 'Maximum number for two-body self-consistency:', max_it_2b + write(*,'(1X,A50,1X,E10.5)') 'Convergence threshold for two-body energies:', conv_2b + write(*,*)'---------------------------------------------------------------' + write(*,*) + + if(linearize) then + write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' + write(*,*) + else + write(*,*) ' *** Quasiparticle energies obtained by root search *** ' + write(*,*) + endif + + ! Memory allocation + + allocate(old_eh_Om(nS),old_ee_Om(nVV),old_hh_Om(nOO)) + allocate(eh_rho(nOrb,nOrb,nS),ee_rho(nOrb,nOrb,nVV),hh_rho(nOrb,nOrb,nOO)) + +! Initialization + + n_it_1b = 0 + err_1b = 1d0 + + n_it_2b = 0 + err_2b = 1d0 + + eQP(:) = eHF(:) + eOld(:) = eHF(:) + + eh_rho(:,:,:) = 0d0 + ee_rho(:,:,:) = 0d0 + hh_rho(:,:,:) = 0d0 + + old_eh_Om(:) = 1d0 + old_ee_Om(:) = 1d0 + old_hh_Om(:) = 1d0 + + !-----------------------------------------! + ! Main loop for one-body self-consistency ! + !-----------------------------------------! + + do while(err_1b > conv_1b .and. n_it_1b < max_it_1b) + + n_it_1b = n_it_1b + 1 + call wall_time(start_1b) + + write(*,*) + write(*,*)'=====================================' + write(*,'(1X,A30,1X,I4)') 'One-body iteration #',n_it_1b + write(*,*)'=====================================' + write(*,*) + + !-----------------------------------------! + ! Main loop for two-body self-consistency ! + !-----------------------------------------! + + do while(err_2b > conv_2b .and. n_it_2b < max_it_2b) + + n_it_2b = n_it_2b + 1 + call wall_time(start_2b) + + write(*,*)' ***********************************' + write(*,'(1X,A30,1X,I4)') 'Two-body iteration #',n_it_2b + write(*,*)' ***********************************' + write(*,*) + + !--------------------------------! + ! Compute effective interactions ! + !--------------------------------! + + ! Memory allocation + allocate(eh_Gam(nOrb,nOrb,nOrb,nOrb)) + allocate(pp_Gam(nOrb,nOrb,nOrb,nOrb)) + + ! Build eh effective interaction + write(*,*) 'Computing eh effective interaction...' + + call wall_time(start_t) + !call R_eh_Gamma(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + ! old_eh_Om,eh_rho,old_ee_Om,ee_rho,old_hh_Om,hh_rho, & + ! eh_trip_Gam) + call wall_time(end_t) + t = end_t - start_t + + write(*,'(A50,1X,F9.3,A8)') 'Wall time for eh Gamma =',t,' seconds' + write(*,*) + + ! Build singlet pp effective interaction + write(*,*) 'Computing pp effective interaction...' + + call wall_time(start_t) + !call R_pp_Gamma(nOrb,nC,nR,nS,old_eh_Om,eh_rho,pp_Gam) + call wall_time(end_t) + t = end_t - start_t + + write(*,'(A50,1X,F9.3,A8)') 'Wall time for pp Gamma =',t,' seconds' + write(*,*) + + + + call wall_time(end_2b) + t_2b = end_2b - start_2b + write(*,'(A50,1X,F9.3,A8)') 'Wall time for two-body iteration =',t_2b,' seconds' + write(*,*) + + end do + !---------------------------------------------! + ! End main loop for two-body self-consistency ! + !---------------------------------------------! + + ! Did it actually converge? + + if(n_it_2b == max_it_2b) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Two-body convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + stop + + else + + write(*,*) + write(*,*)'****************************************************' + write(*,*)' Two-body convergence success ' + write(*,*)'****************************************************' + write(*,*) + + call print_excitation_energies('phBSE@Parquet','1h1p',nS,old_eh_Om) + call print_excitation_energies('ppBSE@Parquet','2p',nVV,old_ee_Om) + call print_excitation_energies('ppBSE@Parquet','2h',nOO,old_hh_Om) + + end if + + allocate(eQPlin(nOrb),Z(nOrb),SigC(nOrb)) + + write(*,*) 'Building self-energy' + + call wall_time(start_t) + !call G_irred_Parquet_self_energy(nOrb,nC,nO,nV,nR,eOld,EcGM,SigC,Z) + call wall_time(end_t) + t = end_t - start_t + write(*,'(A50,1X,F9.3,A8)') 'Wall time for self energy =',t,' seconds' + write(*,*) + + eQPlin(:) = eHF(:) !+ Z(:)*SigC(:) + + ! Solve the quasi-particle equation + + if(linearize) then + + eQP(:) = eQPlin(:) + + else + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Newton-Raphson for Dyson equation not implemented ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + stop + + end if + + deallocate(eQPlin,Z,SigC) + + ! Check one-body converge + + err_1b = maxval(abs(eOld - eQP)) + eOld(:) = eQP(:) + + call wall_time(end_1b) + t_1b = end_1b - start_1b + write(*,'(A50,1X,F9.3,A8)') 'Wall time for one-body iteration =',t_1b,' seconds' + + end do + !---------------------------------------------! + ! End main loop for one-body self-consistency ! + !---------------------------------------------! + + ! Did it actually converge? + if(n_it_1b == max_it_1b) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' One-body convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + stop + + else + + write(*,*) + write(*,*)'****************************************************' + write(*,*)' One-body convergence success ' + write(*,*)'****************************************************' + write(*,*) + + end if + +end subroutine GParquet diff --git a/src/Parquet/RParquet.f90 b/src/Parquet/RParquet.f90 index d24d8d6..a146ee0 100644 --- a/src/Parquet/RParquet.f90 +++ b/src/Parquet/RParquet.f90 @@ -37,16 +37,13 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, integer :: nVVs,nVVt double precision :: EcRPA - double precision,allocatable :: Aph(:,:) - double precision,allocatable :: Bph(:,:) + double precision,allocatable :: Aph(:,:), Bph(:,:) double precision,allocatable :: sing_XpY(:,:),trip_XpY(:,:) double precision,allocatable :: sing_XmY(:,:),trip_XmY(:,:) double precision,allocatable :: eh_sing_Om(:), old_eh_sing_Om(:) double precision,allocatable :: eh_trip_Om(:), old_eh_trip_Om(:) - double precision,allocatable :: Bpp(:,:) - double precision,allocatable :: Cpp(:,:) - double precision,allocatable :: Dpp(:,:) + double precision,allocatable :: Bpp(:,:), Cpp(:,:), Dpp(:,:) double precision,allocatable :: X1s(:,:),X1t(:,:) double precision,allocatable :: Y1s(:,:),Y1t(:,:) double precision,allocatable :: ee_sing_Om(:), old_ee_sing_Om(:) @@ -73,7 +70,8 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, double precision :: EcGM ! Output variables - +! None + nOOs = nO*(nO + 1)/2 nVVs = nV*(nV + 1)/2 nOOt = nO*(nO - 1)/2 @@ -166,7 +164,6 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, n_it_2b = n_it_2b + 1 call wall_time(start_2b) - !TODO add some timers everywhere write(*,*)' ***********************************' write(*,'(1X,A30,1X,I4)') 'Two-body iteration #',n_it_2b write(*,*)' ***********************************' @@ -659,4 +656,4 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, end if -end subroutine +end subroutine RParquet diff --git a/src/Parquet/R_screened_integrals.f90 b/src/Parquet/R_screened_integrals.f90 index 2011dd8..de2e8f7 100644 --- a/src/Parquet/R_screened_integrals.f90 +++ b/src/Parquet/R_screened_integrals.f90 @@ -31,6 +31,7 @@ subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Gam,XpY,r rho(p,q,ia) = rho(p,q,ia) & + (2d0*ERI(p,j,q,b) - ERI(p,j,b,q))*XpY(ia,jb) & + 1d0*eh_sing_Gam(p,j,q,b)*XpY(ia,jb) + end do end do end do diff --git a/src/QuAcK/GQuAcK.f90 b/src/QuAcK/GQuAcK.f90 index 30fa842..7dcc441 100644 --- a/src/QuAcK/GQuAcK.f90 +++ b/src/QuAcK/GQuAcK.f90 @@ -7,7 +7,8 @@ subroutine GQuAcK(working_dir,dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dop TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, & maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, & - dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS) + dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS, & + max_it_macro,conv_one_body,max_it_micro,conv_two_body) implicit none include 'parameters.h' @@ -74,6 +75,9 @@ subroutine GQuAcK(working_dir,dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dop logical,intent(in) :: dophBSE,dophBSE2,doppBSE,dBSE,dTDA logical,intent(in) :: doACFDT,exchange_kernel,doXBS + integer,intent(in) :: max_it_macro,max_it_micro + double precision,intent(in) :: conv_one_body,conv_two_body + ! Local variables logical :: doMP,doCC,doRPA,doGF,doGW,doGT @@ -341,9 +345,9 @@ subroutine GQuAcK(working_dir,dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dop if(doParquet) then call wall_time(start_Parquet) - write(*,'(A65,1X,F9.3,A8)') 'The Parquet method is not implemented in spin-orbital yet :(' - write(*,'(A65,1X,F9.3,A8)') 'Try running the RHF version!' - write(*,*) + call GParquet(max_it_macro,conv_one_body,max_it_micro,conv_two_body, & + nBas,nC,nO,nV,nR,nS, & + eHF,ERI_MO) call wall_time(end_Parquet) t_Parquet = end_Parquet - start_Parquet diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index a9477ca..03abb0d 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -289,14 +289,14 @@ program QuAcK if(doGQuAcK) & call GQuAcK(working_dir,doGtest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, & - doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2,doG0T0pp,doevGTpp,doqsGTpp,doParquet, & + doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2,doG0T0pp,doevGTpp,doqsGTpp,doParquet, & nNuc,nBas,sum(nC),sum(nO),sum(nV),sum(nR),ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO, & maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, & maxSCF_CC,max_diis_CC,thresh_CC,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, & maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, & - dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS) - + dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS, & + max_it_macro,conv_one_body,max_it_micro,conv_two_body) !--------------------------! ! Bogoliubov QuAcK branch ! diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index 63dc813..2bf23f7 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -112,6 +112,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2, double precision,allocatable :: ERI_MO(:,:,:,:) integer :: ixyz integer :: nS + double precision,allocatable :: eGW(:) write(*,*) write(*,*) '******************************' @@ -130,6 +131,8 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2, allocate(dipole_int_MO(nOrb,nOrb,ncart)) allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb)) + allocate(eGW(nOrb)) + allocate(ERI_AO(nBas,nBas,nBas,nBas)) call wall_time(start_int) call read_2e_integrals(working_dir,nBas,ERI_AO) @@ -337,7 +340,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2, call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF_GW,thresh_GW,max_diis_GW, & doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, & lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T, & - V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF,eGW) call wall_time(end_GW) t_GW = end_GW - start_GW