10
1
mirror of https://github.com/pfloos/quack synced 2025-04-01 06:21:37 +02:00

DIIS in parquet - working on it

This commit is contained in:
Pierre-Francois Loos 2025-03-25 10:31:17 +01:00
parent 86c42b7001
commit 1e0ed28f67
6 changed files with 103 additions and 26 deletions

View File

@ -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)

View File

@ -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

View File

@ -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

View File

@ -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

View File

@ -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 !
!----------!

View File

@ -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