From 1e0ed28f67a118131955143dfca2612251ef0495 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 25 Mar 2025 10:31:17 +0100 Subject: [PATCH] DIIS in parquet - working on it --- src/GW/GG0W0.f90 | 6 ++- src/GW/GGW.f90 | 8 +++- src/GW/RG0W0.f90 | 4 +- src/GW/RGW.f90 | 3 +- src/Parquet/GParquet.f90 | 98 ++++++++++++++++++++++++++++++++-------- src/QuAcK/GQuAcK.f90 | 10 ++-- 6 files changed, 103 insertions(+), 26 deletions(-) diff --git a/src/GW/GG0W0.f90 b/src/GW/GG0W0.f90 index 780ec2e..a58989c 100644 --- a/src/GW/GG0W0.f90 +++ b/src/GW/GG0W0.f90 @@ -1,5 +1,5 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & - linearize,eta,doSRG,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) + linearize,eta,doSRG,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF,eGW_out) ! Perform G0W0 calculation implicit none @@ -58,6 +58,8 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! Output variables + double precision,intent(out) :: eGW_out(nBas) + ! Hello world write(*,*) @@ -157,6 +159,8 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA call print_GG0W0(nBas,nO,eHF,ENuc,EGHF,SigC,Z,eGW,EcRPA,EcGM) + eGW_out(:) = eGW(:) + ! Deallocate memory deallocate(SigC,Z,Om,XpY,XmY,rho) diff --git a/src/GW/GGW.f90 b/src/GW/GGW.f90 index adacd18..b4ab3a2 100644 --- a/src/GW/GGW.f90 +++ b/src/GW/GGW.f90 @@ -1,6 +1,6 @@ subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE, & TDA_W,TDA,dBSE,dTDA,linearize,eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,V,Hc, & - ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,eHF) + ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,eHF,eGW) ! GW module @@ -63,6 +63,10 @@ subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchan double precision :: start_GW ,end_GW ,t_GW +! Output variables + + double precision,intent(out) :: eGW(nBas2) + !------------------------------------------------------------------------ ! Perform G0W0 calculatiom !------------------------------------------------------------------------ @@ -71,7 +75,7 @@ subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchan call wall_time(start_GW) call GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & - linearize,eta,doSRG,nBas2,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) + linearize,eta,doSRG,nBas2,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF,eGW) call wall_time(end_GW) t_GW = end_GW - start_GW diff --git a/src/GW/RG0W0.f90 b/src/GW/RG0W0.f90 index 1be8512..813322a 100644 --- a/src/GW/RG0W0.f90 +++ b/src/GW/RG0W0.f90 @@ -61,8 +61,10 @@ 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 + + double precision,intent(out) :: eGW_out(nOrb) ! Output variables diff --git a/src/GW/RGW.f90 b/src/GW/RGW.f90 index 23eadc6..0c683d6 100644 --- a/src/GW/RGW.f90 +++ b/src/GW/RGW.f90 @@ -70,8 +70,9 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_dii logical :: doccG0W0,doccGW +! Output variables - double precision,intent(inout) :: eGW(nOrb) + double precision,intent(out) :: eGW(nOrb) !------------------------------------------------------------------------ ! Perform G0W0 calculation diff --git a/src/Parquet/GParquet.f90 b/src/Parquet/GParquet.f90 index 2e411ff..c02269e 100644 --- a/src/Parquet/GParquet.f90 +++ b/src/Parquet/GParquet.f90 @@ -52,15 +52,38 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, double precision,allocatable :: SigC(:) double precision,allocatable :: Z(:) double precision :: EcGM - + + integer :: max_diis,n_diis + double precision :: rcond + double precision,allocatable :: err_diis(:,:) + double precision,allocatable :: Om_diis(:,:) + double precision,allocatable :: err(:) + double precision,allocatable :: Om(:) + ! Output variables ! None +! Useful parameters + nOO = nO*(nO - 1)/2 nVV = nV*(nV - 1)/2 allocate(eQP(nOrb),eOld(nOrb)) +! DIIS parameters + + max_diis = 10 + n_diis = 0 + rcond = 0d0 + + allocate(err_diis(nS+nOO+nVV,max_diis),Om_diis(nS+nOO+nVV,max_diis)) + allocate(err(nS+nOO+nVV),Om(nS+nOO+nVV)) + + err_diis(:,:) = 0d0 + Om_diis(:,:) = 0d0 + +! Start + write(*,*) write(*,*)'***********************************' write(*,*)'* Generalized Parquet Calculation *' @@ -108,9 +131,9 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, ee_rho(:,:,:) = 0d0 hh_rho(:,:,:) = 0d0 - old_eh_Om(:) = 1d0 - old_ee_Om(:) = 1d0 - old_hh_Om(:) = 1d0 + old_eh_Om(:) = 0d0 + old_ee_Om(:) = 0d0 + old_hh_Om(:) = 0d0 !-----------------------------------------! ! Main loop for one-body self-consistency ! @@ -152,26 +175,42 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, ! Build eh effective interaction write(*,*) 'Computing eh effective interaction...' - call wall_time(start_t) - call G_eh_Gamma(nOrb,nC,nO,nV,nR,nS,nOO,nVV, & - old_eh_Om,eh_rho,old_ee_Om,ee_rho,old_hh_Om,hh_rho, & - eh_Gam) - call wall_time(end_t) - t = end_t - start_t + if(n_it_2b == 1) then + + eh_Gam(:,:,:,:) = 0d0 - write(*,'(A50,1X,F9.3,A8)') 'Wall time for eh Gamma =',t,' seconds' - write(*,*) + else + + call wall_time(start_t) + call G_eh_Gamma(nOrb,nC,nO,nV,nR,nS,nOO,nVV, & + old_eh_Om,eh_rho,old_ee_Om,ee_rho,old_hh_Om,hh_rho, & + eh_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(*,*) + + end if ! Build singlet pp effective interaction write(*,*) 'Computing pp effective interaction...' - call wall_time(start_t) - call G_pp_Gamma(nOrb,nC,nO,nV,nR,nS,old_eh_Om,eh_rho,pp_Gam) - call wall_time(end_t) - t = end_t - start_t + if(n_it_2b == 1) then - write(*,'(A50,1X,F9.3,A8)') 'Wall time for pp Gamma =',t,' seconds' - write(*,*) + pp_Gam(:,:,:,:) = 0d0 + + else + + call wall_time(start_t) + call G_pp_Gamma(nOrb,nC,nO,nV,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(*,*) + + end if !-----------------! ! eh channel ! @@ -292,6 +331,29 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, write(*,*) '----------------------------------------' write(*,*) + !--------------------! + ! DIIS extrapolation ! + !--------------------! + + err( 1:nS ) = eh_Om(:) - old_eh_Om(:) + err( nS+1:nS+nVV ) = ee_Om(:) - old_ee_Om(:) + err(nVV+nS+1:nS+nVV+nOO) = hh_Om(:) - old_hh_Om(:) + + Om( 1:nS ) = eh_Om(:) + Om( nS+1:nS+nVV ) = ee_Om(:) + Om(nVV+nS+1:nS+nVV+nOO) = hh_Om(:) + + if(max_diis > 1) then + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nS+nOO+nVV,nS+nOO+nVV,n_diis,err_diis,Om_diis,err,Om) + + end if + + eh_Om(:) = Om( 1:nS ) + ee_Om(:) = Om( nS+1:nS+nVV ) + hh_Om(:) = Om(nVV+nS+1:nS+nVV+nOO) + !----------! ! Updating ! !----------! diff --git a/src/QuAcK/GQuAcK.f90 b/src/QuAcK/GQuAcK.f90 index c166b9d..b6d7534 100644 --- a/src/QuAcK/GQuAcK.f90 +++ b/src/QuAcK/GQuAcK.f90 @@ -106,6 +106,8 @@ subroutine GQuAcK(working_dir,dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dop integer :: nBas2 integer :: nS + double precision,allocatable :: eGW(:) + write(*,*) write(*,*) '*******************************' write(*,*) '* Generalized Branch of QuAcK *' @@ -121,6 +123,7 @@ subroutine GQuAcK(working_dir,dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dop allocate(cHF(nBas2,nBas2),eHF(nBas2),PHF(nBas2,nBas2),FHF(nBas2,nBas2), & dipole_int_MO(nBas2,nBas2,ncart),ERI_MO(nBas2,nBas2,nBas2,nBas2)) + allocate(eGW(nBas2)) allocate(ERI_AO(nBas,nBas,nBas,nBas)) call wall_time(start_int) @@ -307,9 +310,10 @@ subroutine GQuAcK(working_dir,dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dop if(doGW) then call wall_time(start_GW) - call GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT,exchange_kernel,doXBS, & - dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc, & - nBas,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + call GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT,exchange_kernel,doXBS, & + dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc, & + nBas,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,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