diff --git a/input/options.default b/input/options.default index 93be274..1dcd72d 100644 --- a/input/options.default +++ b/input/options.default @@ -18,5 +18,5 @@ F F F F T # HFB: temperature sigma chem_pot_HF restart_HFB 0.05 1.00 T F -# Parquet: TDAeh TDApp max_it_1b conv_1b max_it_2b conv_2b lin reg - F F 1 0.00001 1 0.00001 F 0.0 +# Parquet: TDAeh TDApp max_it_1b conv_1b max_it_2b conv_2b DIIS_1b DIIS_2b lin reg + T T 10 0.00001 10 0.00001 2 2 T 100.0 diff --git a/src/Parquet/GParquet.f90 b/src/Parquet/GParquet.f90 index 266c671..d8d26b9 100644 --- a/src/Parquet/GParquet.f90 +++ b/src/Parquet/GParquet.f90 @@ -1,6 +1,7 @@ -subroutine GParquet(TDAeh,TDApp,linearize,eta,ENuc,max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,EGHF,eHF,ERI) +subroutine GParquet(TDAeh,TDApp,max_diis_1b,max_diis_2b,linearize,eta,ENuc,max_it_1b,conv_1b,max_it_2b,conv_2b, & + nOrb,nC,nO,nV,nR,nS,EGHF,eHF,ERI) -! Parquet approximation based on restricted orbitals +! Parquet approximation based on spin orbitals implicit none include 'parameters.h' @@ -14,6 +15,8 @@ subroutine GParquet(TDAeh,TDApp,linearize,eta,ENuc,max_it_1b,conv_1b,max_it_2b,c logical,intent(in) :: TDAeh logical,intent(in) :: TDApp + integer,intent(in) :: max_diis_1b + integer,intent(in) :: max_diis_2b logical,intent(in) :: linearize double precision,intent(in) :: eta double precision,intent(in) :: ENuc @@ -61,9 +64,9 @@ subroutine GParquet(TDAeh,TDApp,linearize,eta,ENuc,max_it_1b,conv_1b,max_it_2b,c double precision,allocatable :: Z(:) double precision :: EcGM ! DIIS - integer :: max_diis,n_diis + integer :: n_diis_2b double precision :: rcond - double precision,allocatable :: err_diis(:,:) + double precision,allocatable :: err_diis_2b(:,:) double precision,allocatable :: Phi_diis(:,:) double precision,allocatable :: err(:) double precision,allocatable :: Phi(:) @@ -82,14 +85,12 @@ subroutine GParquet(TDAeh,TDApp,linearize,eta,ENuc,max_it_1b,conv_1b,max_it_2b,c ! DIIS parameters - max_diis = 1 - n_diis = 0 rcond = 1d0 - allocate(err_diis(2*nOrb**4,max_diis),Phi_diis(2*nOrb**4,max_diis)) + allocate(err_diis_2b(2*nOrb**4,max_diis_2b),Phi_diis(2*nOrb**4,max_diis_2b)) allocate(err(2*nOrb**4),Phi(2*nOrb**4)) - err_diis(:,:) = 0d0 + err_diis_2b(:,:) = 0d0 Phi_diis(:,:) = 0d0 ! Start @@ -109,11 +110,13 @@ subroutine GParquet(TDAeh,TDApp,linearize,eta,ENuc,max_it_1b,conv_1b,max_it_2b,c write(*,'(1X,A50,1X,E10.5)') 'Convergence threshold for one-body energies:',conv_1b write(*,'(1X,A50,1X,L5)') 'Linearization of quasiparticle equation?',conv_1b write(*,'(1X,A50,1X,E10.5)') 'Strenght of SRG regularization:',eta + write(*,'(1X,A50,1X,I5)') 'Maximum length of DIIS expansion:',max_diis_1b write(*,*)'---------------------------------------------------------------' write(*,'(1X,A50,1X,I5)') 'Maximum number of two-body iteration:',max_it_2b write(*,'(1X,A50,1X,E10.5)') 'Convergence threshold for two-body energies:',conv_2b write(*,'(1X,A50,1X,L5)') 'TDA for eh excitation energies?',TDAeh write(*,'(1X,A50,1X,L5)') 'TDA for pp excitation energies?',TDApp + write(*,'(1X,A50,1X,I5)') 'Maximum length of DIIS expansion:',max_diis_2b write(*,*)'---------------------------------------------------------------' write(*,*) @@ -382,10 +385,10 @@ subroutine GParquet(TDAeh,TDApp,linearize,eta,ENuc,max_it_1b,conv_1b,max_it_2b,c end do end do - if(max_diis > 1) then + if(max_diis_2b > 1) then - n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,2*nOrb**4,2*nOrb**4,n_diis,err_diis,Phi_diis,err,Phi) + n_diis_2b = min(n_diis_2b+1,max_diis_2b) + call DIIS_extrapolation(rcond,2*nOrb**4,2*nOrb**4,n_diis_2b,err_diis_2b,Phi_diis,err,Phi) end if diff --git a/src/Parquet/RParquet.f90 b/src/Parquet/RParquet.f90 index f519c43..310e4cc 100644 --- a/src/Parquet/RParquet.f90 +++ b/src/Parquet/RParquet.f90 @@ -1,4 +1,5 @@ -subroutine RParquet(TDAeh,TDApp,linearize,eta,ENuc,max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,ERHF,eHF,ERI) +subroutine RParquet(TDAeh,TDApp,max_diis_1b,max_diis_2b,linearize,eta,ENuc,max_it_1b,conv_1b,max_it_2b,conv_2b, & + nOrb,nC,nO,nV,nR,nS,ERHF,eHF,ERI) ! Parquet approximation based on restricted orbitals @@ -14,6 +15,8 @@ subroutine RParquet(TDAeh,TDApp,linearize,eta,ENuc,max_it_1b,conv_1b,max_it_2b,c logical,intent(in) :: TDAeh logical,intent(in) :: TDApp + integer,intent(in) :: max_diis_1b + integer,intent(in) :: max_diis_2b logical,intent(in) :: linearize double precision,intent(in) :: eta double precision,intent(in) :: ENuc diff --git a/src/QuAcK/GQuAcK.f90 b/src/QuAcK/GQuAcK.f90 index d858ac7..9555608 100644 --- a/src/QuAcK/GQuAcK.f90 +++ b/src/QuAcK/GQuAcK.f90 @@ -8,7 +8,7 @@ subroutine GQuAcK(working_dir,dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dop 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, & - TDAeh,TDApp,max_it_1b,conv_1b,max_it_2b,conv_2b,lin_parquet,reg_parquet) + TDAeh,TDApp,max_diis_1b,max_diis_2b,max_it_1b,conv_1b,max_it_2b,conv_2b,lin_parquet,reg_parquet) implicit none include 'parameters.h' @@ -77,6 +77,7 @@ subroutine GQuAcK(working_dir,dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dop integer,intent(in) :: max_it_1b,max_it_2b double precision,intent(in) :: conv_1b,conv_2b + integer,intent(in) :: max_diis_1b,max_diis_2b logical,intent(in) :: TDAeh,TDApp double precision,intent(in) :: reg_parquet logical,intent(in) :: lin_parquet @@ -352,7 +353,8 @@ subroutine GQuAcK(working_dir,dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dop if(doParquet) then call wall_time(start_Parquet) - call GParquet(TDAeh,TDApp,lin_parquet,reg_parquet,ENuc,max_it_1b,conv_1b,max_it_2b,conv_2b,nBas2,nC,nO,nV,nR,nS,EGHF,eHF,ERI_MO) + call GParquet(TDAeh,TDApp,max_diis_1b,max_diis_2b,lin_parquet,reg_parquet,ENuc,max_it_1b,conv_1b,max_it_2b,conv_2b, & + nBas2,nC,nO,nV,nR,nS,EGHF,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 63b507e..a41675a 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -80,6 +80,7 @@ program QuAcK integer :: max_it_1b,max_it_2b double precision :: conv_1b,conv_2b + integer :: max_diis_1b,max_diis_2b logical :: TDAeh,TDApp double precision :: reg_parquet logical :: lin_parquet @@ -149,7 +150,7 @@ program QuAcK doACFDT,exchange_kernel,doXBS, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA, & temperature,sigma,chem_pot_hf,restart_hfb, & - TDAeh,TDApp,max_it_1b,conv_1b,max_it_2b,conv_2b,lin_parquet,reg_parquet) + TDAeh,TDApp,max_diis_1b,max_diis_2b,max_it_1b,conv_1b,max_it_2b,conv_2b,lin_parquet,reg_parquet) !------------------! ! Hardware ! @@ -266,7 +267,7 @@ program QuAcK maxSCF_GF,max_diis_GF,renorm_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, & - TDAeh,TDApp,max_it_1b,conv_1b,max_it_2b,conv_2b,lin_parquet,reg_parquet) + TDAeh,TDApp,max_diis_1b,max_diis_2b,max_it_1b,conv_1b,max_it_2b,conv_2b,lin_parquet,reg_parquet) endif endif @@ -299,7 +300,7 @@ program QuAcK 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, & - TDAeh,TDApp,max_it_1b,conv_1b,max_it_2b,conv_2b,lin_parquet,reg_parquet) + TDAeh,TDApp,max_diis_1b,max_diis_2b,max_it_1b,conv_1b,max_it_2b,conv_2b,lin_parquet,reg_parquet) !--------------------------! ! Bogoliubov QuAcK branch ! diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index 20c5642..a47ee60 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -8,7 +8,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2, maxSCF_GF,max_diis_GF,renorm_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, & - TDAeh,TDApp,max_it_1b,conv_1b,max_it_2b,conv_2b,lin_parquet,reg_parquet) + TDAeh,TDApp,max_diis_1b,max_diis_2b,max_it_1b,conv_1b,max_it_2b,conv_2b,lin_parquet,reg_parquet) ! Restricted branch of QuAcK @@ -84,6 +84,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2, integer,intent(in) :: max_it_1b,max_it_2b double precision,intent(in) :: conv_1b,conv_2b + integer,intent(in) :: max_diis_1b,max_diis_2b logical,intent(in) :: TDAeh,TDApp double precision,intent(in) :: reg_parquet logical,intent(in) :: lin_parquet @@ -380,7 +381,8 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2, if(doParquet) then call wall_time(start_Parquet) - call RParquet(TDAeh,TDApp,lin_parquet,reg_parquet,ENuc,max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,ERHF,eHF,ERI_MO) + call RParquet(TDAeh,TDApp,max_diis_1b,max_diis_2b,lin_parquet,reg_parquet,ENuc,max_it_1b,conv_1b,max_it_2b,conv_2b, & + nOrb,nC,nO,nV,nR,nS,ERHF,eHF,ERI_MO) call wall_time(end_Parquet) t_Parquet = end_Parquet - start_Parquet diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 8fe56f6..028fa64 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -9,7 +9,7 @@ subroutine read_options(working_dir, doACFDT,exchange_kernel,doXBS, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA, & temperature,sigma,chem_pot_hf,restart_hfb, & - TDAeh,TDApp,max_it_1b,conv_1b,max_it_2b,conv_2b,lin_parquet,reg_parquet) + TDAeh,TDApp,max_diis_1b,max_diis_2b,max_it_1b,conv_1b,max_it_2b,conv_2b,lin_parquet,reg_parquet) ! Read desired methods @@ -79,11 +79,12 @@ subroutine read_options(working_dir, double precision,intent(out) :: temperature double precision,intent(out) :: sigma - integer :: max_it_1b,max_it_2b - double precision :: conv_1b,conv_2b - logical :: TDAeh,TDApp - double precision :: reg_parquet - logical :: lin_parquet + integer,intent(out) :: max_it_1b,max_it_2b + double precision,intent(out) :: conv_1b,conv_2b + integer,intent(out) :: max_diis_1b,max_diis_2b + logical,intent(out) :: TDAeh,TDApp + double precision,intent(out) :: reg_parquet + logical,intent(out) :: lin_parquet ! Local variables @@ -247,6 +248,8 @@ subroutine read_options(working_dir, TDAeh = .false. TDApp = .false. + max_diis_1b = 1 + max_diis_2b = 1 max_it_1b = 1 conv_1b = 1d-2 max_it_2b = 1 @@ -255,7 +258,7 @@ subroutine read_options(working_dir, reg_parquet = 0d0 read(1,*) - read(1,*) ans1,ans2,max_it_1b,conv_1b,max_it_2b,conv_2b,ans3,reg_parquet + read(1,*) ans1,ans2,max_it_1b,conv_1b,max_it_2b,conv_2b,max_diis_1b,max_diis_2b,ans3,reg_parquet if(ans1 == 'T') TDAeh = .true. if(ans2 == 'T') TDApp = .true. diff --git a/tests/inp/options.RHF b/tests/inp/options.RHF index daea672..07c570a 100644 --- a/tests/inp/options.RHF +++ b/tests/inp/options.RHF @@ -18,5 +18,5 @@ F F F F T # HFB: temperature sigma chem_pot_HF restart_HFB 0.05 1.00 T F -# Parquet: TDAeh TDApp max_it_1b conv_1b max_it_2b conv_2b lin reg - F F 1 0.00001 1 0.00001 F 0.0 +# Parquet: TDAeh TDApp max_it_1b conv_1b max_it_2b conv_2b DIIS_1b DIIS_2b lin reg + T T 10 0.00001 10 0.00001 2 2 T 100.0