From 14318a217851fd85c370e971688a2e81b4a74764 Mon Sep 17 00:00:00 2001 From: Antoine MARIE Date: Mon, 17 Mar 2025 14:56:56 +0100 Subject: [PATCH] update on Parquet --- input/options.default | 2 + src/LR/ppRLR_C.f90 | 4 +- src/LR/ppRLR_D.f90 | 4 +- src/Parquet/README.md | 22 + src/Parquet/RParquet.f90 | 511 ++++++++++++++++++++ src/Parquet/R_eh_singlet_Gam.f90 | 155 ++++++ src/Parquet/R_eh_triplet_Gam.f90 | 155 ++++++ src/Parquet/R_irred_Parquet_self_energy.f90 | 94 ++++ src/Parquet/R_pp_singlet_Gam.f90 | 225 +++++++++ src/Parquet/R_pp_triplet_Gam.f90 | 217 +++++++++ src/Parquet/R_screened_integrals.f90 | 293 +++++++++++ src/QuAcK/GQuAcK.f90 | 4 +- src/QuAcK/QuAcK.f90 | 10 +- src/QuAcK/RQuAcK.f90 | 14 +- src/QuAcK/read_options.f90 | 16 +- 15 files changed, 1714 insertions(+), 12 deletions(-) create mode 100644 src/Parquet/README.md create mode 100644 src/Parquet/RParquet.f90 create mode 100644 src/Parquet/R_eh_singlet_Gam.f90 create mode 100644 src/Parquet/R_eh_triplet_Gam.f90 create mode 100644 src/Parquet/R_irred_Parquet_self_energy.f90 create mode 100644 src/Parquet/R_pp_singlet_Gam.f90 create mode 100644 src/Parquet/R_pp_triplet_Gam.f90 create mode 100644 src/Parquet/R_screened_integrals.f90 diff --git a/input/options.default b/input/options.default index 1a1be58..2e47df4 100644 --- a/input/options.default +++ b/input/options.default @@ -18,3 +18,5 @@ F F F F T # HFB: temperature sigma chem_pot_HF restart_HFB 0.05 1.00 T F +# Parquet: max_it_macro conv_one_body max_it_micro conv_two_body + 1 0.00001 1 0.00001 diff --git a/src/LR/ppRLR_C.f90 b/src/LR/ppRLR_C.f90 index 658d6c5..e9aec8a 100644 --- a/src/LR/ppRLR_C.f90 +++ b/src/LR/ppRLR_C.f90 @@ -32,8 +32,8 @@ subroutine ppRLR_C(ispin,nOrb,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp) ! Define the chemical potential - eF = e(nO) + e(nO+1) -! eF = 0d0 +! eF = e(nO) + e(nO+1) + eF = 0d0 ! Build C matrix for the singlet manifold diff --git a/src/LR/ppRLR_D.f90 b/src/LR/ppRLR_D.f90 index 6878bf6..bcaca6d 100644 --- a/src/LR/ppRLR_D.f90 +++ b/src/LR/ppRLR_D.f90 @@ -30,8 +30,8 @@ subroutine ppRLR_D(ispin,nOrb,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp) ! Define the chemical potential - eF = e(nO) + e(nO+1) -! eF = 0d0 +! eF = e(nO) + e(nO+1) + eF = 0d0 ! Build the D matrix for the singlet manifold diff --git a/src/Parquet/README.md b/src/Parquet/README.md new file mode 100644 index 0000000..bec4abd --- /dev/null +++ b/src/Parquet/README.md @@ -0,0 +1,22 @@ +# Overview of the Parquet implementation + +## Parameters controling the run + +The parameters provided by the user are: +- `max_it_macro` and `max_it_micro` which set the maximum number of iterations of the macro (one-body) and micro (two-body) self-consistent cycles. +- `conv_one_body` and `conv_two_body` which set the convergence threshold of the macro (one-body) and micro (two-body) self-consistent cycles. +- +- + +The hard-coded parameters are: +- `linearize` which control whether the quasiparticle equation will be linearized or not. Note that the Newton-Raphson has not been implemented yet. +- `TDA` which control whether the Tamm-Dancoff approximation is enforced for the BSE problems or not. +- +- + +## Files and their routines +`RParquet.f90` is the main file for the restricted Parquet calculation, it is called by `RQuack.f90`. The main task of this file is to control the self-consistent cycles. + +## TODO list + +- [ ] Write the TODO list diff --git a/src/Parquet/RParquet.f90 b/src/Parquet/RParquet.f90 new file mode 100644 index 0000000..cdb1f07 --- /dev/null +++ b/src/Parquet/RParquet.f90 @@ -0,0 +1,511 @@ +subroutine RParquet(max_it_macro,conv_one_body,max_it_micro,conv_two_body,nOrb,nC,nO,nV,nR,nS,eHF,ERI) + +! Spatial orbital Parquet implementation + 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_macro,max_it_micro + double precision,intent(in) :: conv_one_body,conv_two_body + 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 :: n_it_macro,n_it_micro + double precision :: err_one_body,err_two_body + double precision :: err_eh_sing,err_eh_trip + double precision :: err_hh_sing,err_hh_trip + double precision :: err_ee_sing,err_ee_trip + double precision :: start_t, end_t, t + + integer :: nOOs,nOOt + integer :: nVVs,nVVt + + double precision :: EcRPA + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: 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 :: X1s(:,:),X1t(:,:) + double precision,allocatable :: Y1s(:,:),Y1t(:,:) + double precision,allocatable :: ee_sing_Om(:), old_ee_sing_Om(:) + double precision,allocatable :: ee_trip_Om(:), old_ee_trip_Om(:) + double precision,allocatable :: X2s(:,:),X2t(:,:) + double precision,allocatable :: Y2s(:,:),Y2t(:,:) + double precision,allocatable :: hh_sing_Om(:), old_hh_sing_Om(:) + double precision,allocatable :: hh_trip_Om(:), old_hh_trip_Om(:) + + double precision,allocatable :: eh_sing_rho(:,:,:),eh_trip_rho(:,:,:) + double precision,allocatable :: ee_sing_rho(:,:,:),hh_sing_rho(:,:,:) + double precision,allocatable :: ee_trip_rho(:,:,:),hh_trip_rho(:,:,:) + + double precision,allocatable :: eh_sing_Gam_A(:,:),eh_sing_Gam_B(:,:) + double precision,allocatable :: eh_trip_Gam_A(:,:),eh_trip_Gam_B(:,:) + double precision,allocatable :: pp_sing_Gam_B(:,:),pp_sing_Gam_C(:,:),pp_sing_Gam_D(:,:) + double precision,allocatable :: pp_trip_Gam_B(:,:),pp_trip_Gam_C(:,:),pp_trip_Gam_D(:,:) + double precision,allocatable :: eh_sing_Gam(:,:,:,:),eh_trip_Gam(:,:,:,:) + double precision,allocatable :: pp_sing_Gam(:,:,:,:),pp_trip_Gam(:,:,:,:) + + double precision,allocatable :: eParquetlin(:),eParquet(:),old_eParquet(:) + double precision,allocatable :: SigC(:) + double precision,allocatable :: Z(:) + double precision :: EcGM + +! Output variables + nOOs = nO*(nO + 1)/2 + nVVs = nV*(nV + 1)/2 + nOOt = nO*(nO - 1)/2 + nVVt = nV*(nV - 1)/2 + + allocate(eParquet(nOrb),old_eParquet(nOrb)) + + write(*,*) + write(*,*)'**********************************' + write(*,*)'* Restricted Parquet Calculation *' + write(*,*)'**********************************' + write(*,*) + +! Print parameters + write(*,*)'Parameters for this run:' + write(*,*)'Maximum number of macro iteration:', max_it_macro + write(*,*)'Convergence threshold for one-body energies:', conv_one_body + write(*,*)'Maximum number of micro iteration:', max_it_micro + write(*,*)'Convergence threshold for two-body energies:', conv_two_body + + if (linearize) then + write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' + write(*,*) + else + write(*,*) ' *** Quasiparticle energies obtained by root search *** ' + write(*,*) + endif + +! Initialization + n_it_macro = 1 + err_one_body = 1d0 + n_it_micro = 1 + err_two_body = 1d0 + old_eParquet(:) = eHF(:) + + write(*,*) + write(*,*)'************ Solving initial linear-response problems ************' + write(*,*)'------------------------------------------------------------------' + + !------------------- + ! Density channel + !------------------- + allocate(Aph(nS,nS),Bph(nS,nS),eh_sing_Om(nS),sing_XpY(nS,nS),sing_XmY(nS,nS),old_eh_sing_Om(nS)) + Aph(:,:) = 0d0 + Bph(:,:) = 0d0 + call wall_time(start_t) + if(.not.TDA) call phRLR_B(1,.false.,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) + call phRLR_A(1,.false.,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) + call phRLR(TDA,nS,Aph,Bph,EcRPA,eh_sing_Om,sing_XpY,sing_XmY) + call wall_time(end_t) + t = end_t - start_t + write(*,*) + write(*,'(A51,1X,F9.3,A8)') 'Total wall time for initial singlet phRPA problem =',t,' seconds' + !if (print_phLR) call print_excitation_energies('phRPA@RHF','singlet',nS,eh_sing_Om) + call print_excitation_energies('phRPA@RHF','singlet',nS,eh_sing_Om) + + deallocate(Aph,Bph) + + !------------------- + ! Magnetic channel + !------------------- + allocate(Aph(nS,nS),Bph(nS,nS),eh_trip_Om(nS),trip_XpY(nS,nS),trip_XmY(nS,nS),old_eh_trip_Om(nS)) + Aph(:,:) = 0d0 + Bph(:,:) = 0d0 + call wall_time(start_t) + if(.not.TDA) call phRLR_B(2,.false.,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) + call phRLR_A(2,.false.,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) + call phRLR(TDA,nS,Aph,Bph,EcRPA,eh_trip_Om,trip_XpY,trip_XmY) + call wall_time(end_t) + t = end_t - start_t + write(*,*) + write(*,'(A51,1X,F9.3,A8)') 'Total wall time for initial triplet phRPA problem =',t,' seconds' + !if (print_phLR) call print_excitation_energies('phRPA@RHF','triplet',nS,eh_trip_Om) + call print_excitation_energies('phRPA@RHF','triplet',nS,eh_trip_Om) + + deallocate(Aph,Bph) + + !------------------- + ! Singlet channel + !------------------- + allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs), & + ee_sing_Om(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs),old_ee_sing_Om(nVVs), & + hh_sing_Om(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs),old_hh_sing_Om(nOOs)) + Bpp(:,:)=0d0 + Cpp(:,:)=0d0 + Dpp(:,:)=0d0 + call wall_time(start_t) + if(.not.TDA) call ppRLR_B(1,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) + call ppRLR_C(1,nOrb,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp) + call ppRLR_D(1,nOrb,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp) + call ppRLR(TDA,nOOs,nVVs,Bpp,Cpp,Dpp,ee_sing_Om,X1s,Y1s,hh_sing_Om,X2s,Y2s,EcRPA) + call wall_time(end_t) + t = end_t - start_t + write(*,*) + write(*,'(A51,1X,F9.3,A8)') 'Total wall time for initial singlet ppRPA problem =',t,' seconds' + if (print_ppLR) call print_excitation_energies('ppRPA@RHF','2p (singlets)',nVVs,ee_sing_Om) + if (print_ppLR) call print_excitation_energies('ppRPA@RHF','2h (singlets)',nOOs,hh_sing_Om) + + deallocate(Bpp,Cpp,Dpp) + + !------------------- + ! Triplet channel + !------------------- + allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt), & + ee_trip_Om(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt),old_ee_trip_Om(nVVt), & + hh_trip_Om(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt),old_hh_trip_Om(nOOt)) + Bpp(:,:)=0d0 + Cpp(:,:)=0d0 + Dpp(:,:)=0d0 + call wall_time(start_t) + if(.not.TDA) call ppRLR_B(2,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp) + call ppRLR_C(2,nOrb,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp) + call ppRLR_D(2,nOrb,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp) + call ppRLR(TDA,nOOt,nVVt,Bpp,Cpp,Dpp,ee_trip_Om,X1t,Y1t,hh_trip_Om,X2t,Y2t,EcRPA) + call wall_time(end_t) + t = end_t - start_t + write(*,*) + write(*,'(A51,1X,F9.3,A8)') 'Total wall time for initial triplet ppRPA problem =',t,' seconds' + write(*,*) + if(print_ppLR) call print_excitation_energies('ppRPA@RHF','2p (triplets)',nVVt,ee_trip_Om) + if(print_ppLR) call print_excitation_energies('ppRPA@RHF','2h (triplets)',nOOt,hh_trip_Om) + + deallocate(Bpp,Cpp,Dpp) + + !------------------- + ! Updating + !------------------- + old_eh_sing_Om(:) = eh_sing_Om(:) + old_eh_trip_Om(:) = eh_trip_Om(:) + old_ee_sing_Om(:) = ee_sing_Om(:) + old_hh_sing_Om(:) = hh_sing_Om(:) + old_ee_trip_Om(:) = ee_trip_Om(:) + old_hh_trip_Om(:) = hh_trip_Om(:) + + deallocate(eh_sing_Om,eh_trip_Om,ee_sing_Om,hh_sing_Om,ee_trip_Om,hh_trip_Om) + + allocate(eh_sing_rho(nOrb,nOrb,nS)) + call R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,sing_XpY,eh_sing_rho) + deallocate(sing_XpY,sing_XmY) + + allocate(eh_trip_rho(nOrb,nOrb,nS)) + call R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,trip_XpY,eh_trip_rho) + deallocate(trip_XpY,trip_XmY) + + allocate(ee_sing_rho(nOrb,nOrb,nVVs),hh_sing_rho(nOrb,nOrb,nOOs)) + allocate(pp_sing_Gam(nOrb,nOrb,nOrb,nOrb)) + pp_sing_Gam(:,:,:,:) = 0d0 + call R_pp_singlet_screened_integral(nOrb,nC,nO,nV,nR,nOOs,nVVs,ERI,pp_sing_Gam,X1s,Y1s,ee_sing_rho,X2s,Y2s,hh_sing_rho) + deallocate(X1s,Y1s,X2s,Y2s) + deallocate(pp_sing_Gam) + + allocate(ee_trip_rho(nOrb,nOrb,nVVt),hh_trip_rho(nOrb,nOrb,nOOt)) + allocate(pp_trip_Gam(nOrb,nOrb,nOrb,nOrb)) + pp_trip_Gam(:,:,:,:) = 0d0 + call R_pp_triplet_screened_integral(nOrb,nC,nO,nV,nR,nOOt,nVVt,ERI,pp_trip_Gam,X1t,Y1t,ee_trip_rho,X2t,Y2t,hh_trip_rho) + deallocate(X1t,Y1t,X2t,Y2t) + deallocate(pp_trip_Gam) + +!------------------------------------------------------------------------ +! Loop on one-body energies +!------------------------------------------------------------------------ + + do while(err_one_body > conv_one_body .and. n_it_macro <= max_it_macro) + + write(*,*) + write(*,*)'************ Macro iteration number ',n_it_macro,' ************' + write(*,*)'---------------------------------------------------------------' + write(*,*) + + + do while(err_two_body > conv_two_body .and. n_it_micro <= max_it_micro) + + !TODO add some timers everywhere + write(*,*) + write(*,*)' Micro iteration number ',n_it_micro + write(*,*)' -----------------------------------' + write(*,*) + + !------------------- + ! Density channel + !------------------- + write(*,*)'Diagonalizing singlet ehBSE:' + write(*,*)'----------------------------' + allocate(Aph(nS,nS),Bph(nS,nS),eh_sing_Om(nS),sing_XpY(nS,nS),sing_XmY(nS,nS),eh_sing_Gam_A(nS,nS),eh_sing_Gam_B(nS,nS)) + Aph(:,:) = 0d0 + Bph(:,:) = 0d0 + call wall_time(start_t) + if(.not.TDA) call phRLR_B(1,.false.,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) + call phRLR_A(1,.false.,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) + + call R_eh_singlet_Gamma_A(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + old_eh_sing_Om,eh_sing_rho,old_eh_trip_Om,eh_trip_rho, & + old_ee_sing_Om,ee_sing_rho,old_ee_trip_Om,ee_trip_rho, & + old_hh_sing_Om,hh_sing_rho,old_hh_trip_Om,hh_trip_rho, eh_sing_Gam_A) + call R_eh_singlet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + old_eh_sing_Om,eh_sing_rho,old_eh_trip_Om,eh_trip_rho, & + old_ee_sing_Om,ee_sing_rho,old_ee_trip_Om,ee_trip_rho, & + old_hh_sing_Om,hh_sing_rho,old_hh_trip_Om,hh_trip_rho, eh_sing_Gam_B) + + Aph(:,:) = Aph(:,:) + eh_sing_Gam_A(:,:) + Bph(:,:) = Bph(:,:) + eh_sing_Gam_B(:,:) + call phRLR(TDA,nS,Aph,Bph,EcRPA,eh_sing_Om,sing_XpY,sing_XmY) + call wall_time(end_t) + t = end_t - start_t + write(*,'(A52,1X,F9.3,A8)') 'Total wall time for singlet phBSE problem =',t,' seconds' + write(*,*) + if (print_phLR) call print_excitation_energies('phBSE@Parquet','singlet',nS,eh_sing_Om) + err_eh_sing = maxval(abs(old_eh_sing_Om - eh_sing_Om)) + deallocate(Aph,Bph,eh_sing_Gam_A,eh_sing_Gam_B) + + !------------------- + ! Magnetic channel + !------------------- + write(*,*)'Diagonalizing triplet ehBSE:' + write(*,*)'----------------------------' + allocate(Aph(nS,nS),Bph(nS,nS),eh_trip_Om(nS),trip_XpY(nS,nS),trip_XmY(nS,nS),eh_trip_Gam_A(nS,nS),eh_trip_Gam_B(nS,nS)) + Aph(:,:) = 0d0 + Bph(:,:) = 0d0 + call wall_time(start_t) + if(.not.TDA) call phRLR_B(2,.false.,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) + call phRLR_A(2,.false.,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) + + call R_eh_triplet_Gamma_A(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + old_eh_sing_Om,eh_sing_rho,old_eh_trip_Om,eh_trip_rho, & + old_ee_sing_Om,ee_sing_rho,old_ee_trip_Om,ee_trip_rho, & + old_hh_sing_Om,hh_sing_rho,old_hh_trip_Om,hh_trip_rho, eh_trip_Gam_A) + call R_eh_triplet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + old_eh_sing_Om,eh_sing_rho,old_eh_trip_Om,eh_trip_rho, & + old_ee_sing_Om,ee_sing_rho,old_ee_trip_Om,ee_trip_rho, & + old_hh_sing_Om,hh_sing_rho,old_hh_trip_Om,hh_trip_rho, eh_trip_Gam_B) + + Aph(:,:) = Aph(:,:) + eh_trip_Gam_A(:,:) + Bph(:,:) = Bph(:,:) + eh_trip_Gam_B(:,:) + call phRLR(TDA,nS,Aph,Bph,EcRPA,eh_trip_Om,trip_XpY,trip_XmY) + call wall_time(end_t) + t = end_t - start_t + write(*,'(A52,1X,F9.3,A8)') 'Total wall time for triplet phBSE problem =',t,' seconds' + write(*,*) + if (print_phLR) call print_excitation_energies('phBSE@Parquet','triplet',nS,eh_trip_Om) + err_eh_trip = maxval(abs(old_eh_trip_Om - eh_trip_Om)) + deallocate(Aph,Bph,eh_trip_Gam_A,eh_trip_Gam_B) + + !------------------- + ! Singlet channel + !------------------- + write(*,*)'Diagonalizing singlet ppBSE:' + write(*,*)'----------------------------' + allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs), & + ee_sing_Om(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & + hh_sing_Om(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & + pp_sing_Gam_B(nVVs,nOOs),pp_sing_Gam_C(nVVs,nVVs),pp_sing_Gam_D(nOOs,nOOs)) + Bpp(:,:) = 0d0 + Cpp(:,:) = 0d0 + Dpp(:,:) = 0d0 + call wall_time(start_t) + if(.not.TDA) call ppRLR_B(1,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) + call ppRLR_C(1,nOrb,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp) + call ppRLR_D(1,nOrb,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp) + + call R_pp_singlet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,& + old_eh_sing_Om,eh_sing_rho,old_eh_trip_Om,eh_trip_rho, pp_sing_Gam_B) + call R_pp_singlet_Gamma_C(nOrb,nC,nO,nV,nR,nS,nVVs,old_eh_sing_Om,eh_sing_rho,old_eh_trip_Om,eh_trip_rho, pp_sing_Gam_C) + call R_pp_singlet_Gamma_D(nOrb,nC,nO,nV,nR,nS,nOOs,old_eh_sing_Om,eh_sing_rho,old_eh_trip_Om,eh_trip_rho, pp_sing_Gam_D) + + Bpp(:,:) = Bpp(:,:) + pp_sing_Gam_B(:,:) + Cpp(:,:) = Cpp(:,:) + pp_sing_Gam_C(:,:) + Dpp(:,:) = Dpp(:,:) + pp_sing_Gam_D(:,:) + + call ppRLR(TDA,nOOs,nVVs,Bpp,Cpp,Dpp,ee_sing_Om,X1s,Y1s,hh_sing_Om,X2s,Y2s,EcRPA) + call wall_time(end_t) + t = end_t - start_t + write(*,'(A52,1X,F9.3,A8)') 'Total wall time for singlet ppBSE problem =',t,' seconds' + write(*,*) + if (print_ppLR) call print_excitation_energies('ppBSE@Parquet','2p (singlets)',nVVs,ee_sing_Om) + if (print_ppLR) call print_excitation_energies('ppBSE@Parquet','2h (singlets)',nOOs,hh_sing_Om) + err_ee_sing = maxval(abs(old_ee_sing_Om - ee_sing_Om)) + err_hh_sing = maxval(abs(old_hh_sing_Om - hh_sing_Om)) + deallocate(Bpp,Cpp,Dpp,pp_sing_Gam_B,pp_sing_Gam_C,pp_sing_Gam_D) + + !------------------- + ! Triplet channel + !------------------- + write(*,*)'Diagonalizing triplet ppBSE:' + write(*,*)'----------------------------' + allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt), & + ee_trip_Om(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & + hh_trip_Om(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & + pp_trip_Gam_B(nVVt,nOOt),pp_trip_Gam_C(nVVt,nVVt),pp_trip_Gam_D(nOOt,nOOt)) + Bpp(:,:)=0d0 + Cpp(:,:)=0d0 + Dpp(:,:)=0d0 + call wall_time(start_t) + if(.not.TDA) call ppRLR_B(2,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp) + call ppRLR_C(2,nOrb,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp) + call ppRLR_D(2,nOrb,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp) + + call R_pp_triplet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOt,nVVt,& + old_eh_sing_Om,eh_sing_rho,old_eh_trip_Om,eh_trip_rho, pp_trip_Gam_B) + call R_pp_triplet_Gamma_C(nOrb,nC,nO,nV,nR,nS,nVVt,old_eh_sing_Om,eh_sing_rho,old_eh_trip_Om,eh_trip_rho, pp_trip_Gam_C) + call R_pp_triplet_Gamma_D(nOrb,nC,nO,nV,nR,nS,nOOt,old_eh_sing_Om,eh_sing_rho,old_eh_trip_Om,eh_trip_rho, pp_trip_Gam_D) + + Bpp(:,:) = Bpp(:,:) + pp_trip_Gam_B(:,:) + Cpp(:,:) = Cpp(:,:) + pp_trip_Gam_C(:,:) + Dpp(:,:) = Dpp(:,:) + pp_trip_Gam_D(:,:) + + call ppRLR(TDA,nOOt,nVVt,Bpp,Cpp,Dpp,ee_trip_Om,X1t,Y1t,hh_trip_Om,X2t,Y2t,EcRPA) + call wall_time(end_t) + t = end_t - start_t + write(*,'(A52,1X,F9.3,A8)') 'Total wall time for triplet ppBSE problem =',t,' seconds' + write(*,*) + if(print_ppLR) call print_excitation_energies('ppBSE@Parquet','2p (triplets)',nVVt,ee_trip_Om) + if(print_ppLR) call print_excitation_energies('ppBSE@Parquet','2h (triplets)',nOOt,hh_trip_Om) + err_ee_trip = maxval(abs(old_ee_trip_Om - ee_trip_Om)) + err_hh_trip = maxval(abs(old_hh_trip_Om - hh_trip_Om)) + deallocate(Bpp,Cpp,Dpp,pp_trip_Gam_B,pp_trip_Gam_C,pp_trip_Gam_D) + + write(*,*) + write(*,*)'Error for density channel:', err_eh_sing + write(*,*)'Error for magnetic channel:',err_eh_trip + write(*,*)'Error for singlet channel:', max(err_ee_sing,err_hh_sing) + write(*,*)'Error for triplet channel:', max(err_ee_trip,err_hh_trip) + write(*,*) + + !------------------- + ! Updating + !------------------- + old_eh_sing_Om(:) = eh_sing_Om(:) + old_eh_trip_Om(:) = eh_trip_Om(:) + old_ee_sing_Om(:) = ee_sing_Om(:) + old_hh_sing_Om(:) = hh_sing_Om(:) + old_ee_trip_Om(:) = ee_trip_Om(:) + old_hh_trip_Om(:) = hh_trip_Om(:) + + allocate(pp_trip_Gam(nOrb,nOrb,nOrb,nOrb)) + call wall_time(start_t) + call R_pp_triplet_Gamma(nOrb,nC,nR,nS,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho,pp_trip_Gam) + call wall_time(end_t) + t = end_t - start_t + write(*,'(A52,1X,F9.3,A8)') 'Total wall time for building pp triplet Gamma =',t,' seconds' + write(*,*) + allocate(pp_sing_Gam(nOrb,nOrb,nOrb,nOrb)) + call wall_time(start_t) + call R_pp_singlet_Gamma(nOrb,nC,nR,nS,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho,pp_sing_Gam) + call wall_time(end_t) + t = end_t - start_t + write(*,'(A52,1X,F9.3,A8)') 'Total wall time for building pp singlet Gamma =',t,' seconds' + write(*,*) + + deallocate(eh_sing_Om,eh_trip_Om,ee_sing_Om,hh_sing_Om,ee_trip_Om,hh_trip_Om) + deallocate(eh_sing_rho,eh_trip_rho,ee_sing_rho,ee_trip_rho,hh_sing_rho,hh_trip_rho) + + allocate(eh_sing_rho(nOrb,nOrb,nS)) + call R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,sing_XpY,eh_sing_rho) + deallocate(sing_XpY,sing_XmY) + + allocate(eh_trip_rho(nOrb,nOrb,nS)) + call R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,trip_XpY,eh_trip_rho) + deallocate(trip_XpY,trip_XmY) + + allocate(ee_sing_rho(nOrb,nOrb,nVVs),hh_sing_rho(nOrb,nOrb,nOOs)) + call R_pp_singlet_screened_integral(nOrb,nC,nO,nV,nR,nOOs,nVVs,ERI,pp_sing_Gam,X1s,Y1s,ee_sing_rho,X2s,Y2s,hh_sing_rho) + deallocate(X1s,Y1s,X2s,Y2s,pp_sing_Gam) + + allocate(ee_trip_rho(nOrb,nOrb,nVVt),hh_trip_rho(nOrb,nOrb,nOOt)) + call R_pp_triplet_screened_integral(nOrb,nC,nO,nV,nR,nOOt,nVVt,ERI,pp_trip_Gam,X1t,Y1t,ee_trip_rho,X2t,Y2t,hh_trip_rho) + deallocate(X1t,Y1t,X2t,Y2t,pp_trip_Gam) + + ! Convergence criteria + err_two_body = max(err_eh_sing,err_eh_trip,err_ee_sing,err_ee_trip,err_hh_sing,err_hh_trip) + + n_it_micro = n_it_micro + 1 + end do + !------------------------------------------------------------------------ + ! End main loop + !------------------------------------------------------------------------ + + ! Did it actually converge? + if(n_it_micro == max_it_micro+1) 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','singlet',nS,old_eh_sing_Om) + call print_excitation_energies('phBSE@Parquet','triplet',nS,old_eh_trip_Om) + call print_excitation_energies('ppBSE@Parquet','2p (singlets)',nVVs,old_ee_sing_Om) + call print_excitation_energies('ppBSE@Parquet','2h (singlets)',nOOs,old_hh_sing_Om) + call print_excitation_energies('ppBSE@Parquet','2p (triplets)',nVVt,old_ee_trip_Om) + call print_excitation_energies('ppBSE@Parquet','2h (triplets)',nOOt,old_hh_trip_Om) + end if + + allocate(eParquetlin(nOrb),Z(nOrb),SigC(nOrb)) + write(*,*) 'Building self-energy' + call wall_time(start_t) + call R_irred_Parquet_self_energy(nOrb,nC,nO,nV,nR,old_eParquet,EcGM,SigC,Z) + call wall_time(end_t) + t = end_t - start_t + write(*,'(A52,1X,F9.3,A8)') 'Total wall time for building self energy =',t,' seconds' + write(*,*) + + eParquetlin(:) = eHF(:) !+ Z(:)*SigC(:) + ! Solve the quasi-particle equation + if(linearize) then + eParquet(:) = eParquetlin(:) + else + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Newton-Raphson for Dyson equation not implemented ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + stop + end if + deallocate(eParquetlin,Z,SigC) + + err_one_body = maxval(abs(old_eParquet - eParquet)) + old_eParquet(:) = eParquet(:) + + n_it_macro = n_it_macro + 1 + end do ! End the macro loop + + ! Did it actually converge? + if(n_it_macro == max_it_macro+1) 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 diff --git a/src/Parquet/R_eh_singlet_Gam.f90 b/src/Parquet/R_eh_singlet_Gam.f90 new file mode 100644 index 0000000..546bb5a --- /dev/null +++ b/src/Parquet/R_eh_singlet_Gam.f90 @@ -0,0 +1,155 @@ +subroutine R_eh_singlet_Gamma_A(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho, & + ee_sing_Om,ee_sing_rho,ee_trip_Om,ee_trip_rho, & + hh_sing_Om,hh_sing_rho,hh_trip_Om,hh_trip_rho, eh_sing_Gam_A) + + +! Compute irreducible vertex in the triplet pp channel + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt + double precision,intent(in) :: eh_sing_Om(nS) + double precision,intent(in) :: eh_sing_rho(nOrb,nOrb,nS) + double precision,intent(in) :: eh_trip_Om(nS) + double precision,intent(in) :: eh_trip_rho(nOrb,nOrb,nS) + double precision,intent(in) :: ee_sing_Om(nVVs) + double precision,intent(in) :: ee_sing_rho(nOrb,nOrb,nVVs) + double precision,intent(in) :: ee_trip_Om(nVVt) + double precision,intent(in) :: ee_trip_rho(nOrb,nOrb,nVVt) + double precision,intent(in) :: hh_sing_Om(nOOs) + double precision,intent(in) :: hh_sing_rho(nOrb,nOrb,nOOs) + double precision,intent(in) :: hh_trip_Om(nVVs) + double precision,intent(in) :: hh_trip_rho(nOrb,nOrb,nVVs) + +! Local variables + integer :: i,a,j,b + integer :: ia,jb + integer :: n + double precision,external :: Kronecker_delta + +! Output variables + double precision, intent(out) :: eh_sing_Gam_A(nS,nS) + +! Initialization + eh_sing_Gam_A(:,:) = 0d0 + + ia = 0 + do i=nC+1,nO + do a=nO+1,nOrb-nR + ia = ia + 1 + jb = 0 + do j=nC+1,nO + do b=nO+1,norb-nR + jb = jb + 1 + + do n=1,nS + eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) & + + eh_sing_rho(a,b,n)*eh_sing_rho(j,i,n)/eh_sing_Om(n) & + + 3d0 * eh_trip_rho(a,b,n)*eh_trip_rho(j,i,n)/eh_trip_Om(n) + end do + + do n=1,nVVs + eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) & + + ee_sing_rho(a,j,n)*ee_sing_rho(i,b,n)/ee_sing_Om(n) + end do + + do n=1,nOOs + eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) & + - hh_sing_rho(a,j,n)*hh_sing_rho(i,b,n)/hh_sing_Om(n) + end do + + do n=1,nVVt + eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) & + + 3d0 * ee_trip_rho(a,j,n)*ee_trip_rho(i,b,n)/ee_trip_Om(n) + end do + + do n=1,nOOt + eh_sing_Gam_A(ia,jb) = eh_sing_Gam_A(ia,jb) & + - 3d0 * hh_trip_rho(a,j,n)*hh_trip_rho(i,b,n)/hh_trip_Om(n) + end do + + enddo + enddo + enddo + enddo + +end subroutine R_eh_singlet_Gamma_A + +subroutine R_eh_singlet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho, & + ee_sing_Om,ee_sing_rho,ee_trip_Om,ee_trip_rho, & + hh_sing_Om,hh_sing_rho,hh_trip_Om,hh_trip_rho, eh_sing_Gam_B) + + +! Compute irreducible vertex in the triplet pp channel + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt + double precision,intent(in) :: eh_sing_Om(nS) + double precision,intent(in) :: eh_sing_rho(nOrb,nOrb,nS) + double precision,intent(in) :: eh_trip_Om(nS) + double precision,intent(in) :: eh_trip_rho(nOrb,nOrb,nS) + double precision,intent(in) :: ee_sing_Om(nVVs) + double precision,intent(in) :: ee_sing_rho(nOrb,nOrb,nVVs) + double precision,intent(in) :: ee_trip_Om(nVVt) + double precision,intent(in) :: ee_trip_rho(nOrb,nOrb,nVVt) + double precision,intent(in) :: hh_sing_Om(nOOs) + double precision,intent(in) :: hh_sing_rho(nOrb,nOrb,nOOs) + double precision,intent(in) :: hh_trip_Om(nVVs) + double precision,intent(in) :: hh_trip_rho(nOrb,nOrb,nVVs) + +! Local variables + integer :: i,a,j,b + integer :: ia,jb + integer :: n + double precision,external :: Kronecker_delta + +! Output variables + double precision, intent(out) :: eh_sing_Gam_B(nS,nS) + +! Initialization + eh_sing_Gam_B(:,:) = 0d0 + + ia = 0 + do i=nC+1,nO + do a=nO+1,nOrb-nR + ia = ia + 1 + jb = 0 + do j=nC+1,nO + do b=nO+1,norb-nR + jb = jb + 1 + + do n=1,nS + eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) & + + eh_sing_rho(a,j,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) & + + 3d0 * eh_trip_rho(a,j,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n) + end do + + do n=1,nVVs + eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) & + + ee_sing_rho(a,b,n)*ee_sing_rho(i,j,n)/ee_sing_Om(n) + end do + + do n=1,nOOs + eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) & + - hh_sing_rho(a,b,n)*hh_sing_rho(i,j,n)/hh_sing_Om(n) + end do + + do n=1,nVVt + eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) & + + 3d0 * ee_trip_rho(a,b,n)*ee_trip_rho(i,j,n)/ee_trip_Om(n) + end do + + do n=1,nOOt + eh_sing_Gam_B(ia,jb) = eh_sing_Gam_B(ia,jb) & + - 3d0 * hh_trip_rho(a,b,n)*hh_trip_rho(i,j,n)/hh_trip_Om(n) + end do + + enddo + enddo + enddo + enddo + +end subroutine R_eh_singlet_Gamma_B diff --git a/src/Parquet/R_eh_triplet_Gam.f90 b/src/Parquet/R_eh_triplet_Gam.f90 new file mode 100644 index 0000000..eb019c8 --- /dev/null +++ b/src/Parquet/R_eh_triplet_Gam.f90 @@ -0,0 +1,155 @@ +subroutine R_eh_triplet_Gamma_A(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho, & + ee_sing_Om,ee_sing_rho,ee_trip_Om,ee_trip_rho, & + hh_sing_Om,hh_sing_rho,hh_trip_Om,hh_trip_rho, eh_trip_Gam_A) + + +! Compute irreducible vertex in the triplet pp channel + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt + double precision,intent(in) :: eh_sing_Om(nS) + double precision,intent(in) :: eh_sing_rho(nOrb,nOrb,nS) + double precision,intent(in) :: eh_trip_Om(nS) + double precision,intent(in) :: eh_trip_rho(nOrb,nOrb,nS) + double precision,intent(in) :: ee_sing_Om(nVVs) + double precision,intent(in) :: ee_sing_rho(nOrb,nOrb,nVVs) + double precision,intent(in) :: ee_trip_Om(nVVt) + double precision,intent(in) :: ee_trip_rho(nOrb,nOrb,nVVt) + double precision,intent(in) :: hh_sing_Om(nOOs) + double precision,intent(in) :: hh_sing_rho(nOrb,nOrb,nOOs) + double precision,intent(in) :: hh_trip_Om(nVVs) + double precision,intent(in) :: hh_trip_rho(nOrb,nOrb,nVVs) + +! Local variables + integer :: i,a,j,b + integer :: ia,jb + integer :: n + double precision,external :: Kronecker_delta + +! Output variables + double precision, intent(out) :: eh_trip_Gam_A(nS,nS) + +! Initialization + eh_trip_Gam_A(:,:) = 0d0 + + ia = 0 + do i=nC+1,nO + do a=nO+1,nOrb-nR + ia = ia + 1 + jb = 0 + do j=nC+1,nO + do b=nO+1,norb-nR + jb = jb + 1 + + do n=1,nS + eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) & + + eh_sing_rho(a,b,n)*eh_sing_rho(j,i,n)/eh_sing_Om(n) & + - eh_trip_rho(a,b,n)*eh_trip_rho(j,i,n)/eh_trip_Om(n) + end do + + do n=1,nVVs + eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) & + - ee_sing_rho(a,j,n)*ee_sing_rho(i,b,n)/ee_sing_Om(n) + end do + + do n=1,nOOs + eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) & + + hh_sing_rho(a,j,n)*hh_sing_rho(i,b,n)/hh_sing_Om(n) + end do + + do n=1,nVVt + eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) & + + ee_trip_rho(a,j,n)*ee_trip_rho(i,b,n)/ee_trip_Om(n) + end do + + do n=1,nOOt + eh_trip_Gam_A(ia,jb) = eh_trip_Gam_A(ia,jb) & + - hh_trip_rho(a,j,n)*hh_trip_rho(i,b,n)/hh_trip_Om(n) + end do + + enddo + enddo + enddo + enddo + +end subroutine R_eh_triplet_Gamma_A + +subroutine R_eh_triplet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho, & + ee_sing_Om,ee_sing_rho,ee_trip_Om,ee_trip_rho, & + hh_sing_Om,hh_sing_rho,hh_trip_Om,hh_trip_rho, eh_trip_Gam_B) + + +! Compute irreducible vertex in the triplet pp channel + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt + double precision,intent(in) :: eh_sing_Om(nS) + double precision,intent(in) :: eh_sing_rho(nOrb,nOrb,nS) + double precision,intent(in) :: eh_trip_Om(nS) + double precision,intent(in) :: eh_trip_rho(nOrb,nOrb,nS) + double precision,intent(in) :: ee_sing_Om(nVVs) + double precision,intent(in) :: ee_sing_rho(nOrb,nOrb,nVVs) + double precision,intent(in) :: ee_trip_Om(nVVt) + double precision,intent(in) :: ee_trip_rho(nOrb,nOrb,nVVt) + double precision,intent(in) :: hh_sing_Om(nOOs) + double precision,intent(in) :: hh_sing_rho(nOrb,nOrb,nOOs) + double precision,intent(in) :: hh_trip_Om(nVVs) + double precision,intent(in) :: hh_trip_rho(nOrb,nOrb,nVVs) + +! Local variables + integer :: i,a,j,b + integer :: ia,jb + integer :: n + double precision,external :: Kronecker_delta + +! Output variables + double precision, intent(out) :: eh_trip_Gam_B(nS,nS) + +! Initialization + eh_trip_Gam_B(:,:) = 0d0 + + ia = 0 + do i=nC+1,nO + do a=nO+1,nOrb-nR + ia = ia + 1 + jb = 0 + do j=nC+1,nO + do b=nO+1,norb-nR + jb = jb + 1 + + do n=1,nS + eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) & + + eh_sing_rho(a,j,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) & + - eh_trip_rho(a,j,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n) + end do + + do n=1,nVVs + eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) & + - ee_sing_rho(a,b,n)*ee_sing_rho(i,j,n)/ee_sing_Om(n) + end do + + do n=1,nOOs + eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) & + + hh_sing_rho(a,b,n)*hh_sing_rho(i,j,n)/hh_sing_Om(n) + end do + + do n=1,nVVt + eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) & + + ee_trip_rho(a,b,n)*ee_trip_rho(i,j,n)/ee_trip_Om(n) + end do + + do n=1,nOOt + eh_trip_Gam_B(ia,jb) = eh_trip_Gam_B(ia,jb) & + - hh_trip_rho(a,b,n)*hh_trip_rho(i,j,n)/hh_trip_Om(n) + end do + + enddo + enddo + enddo + enddo + +end subroutine R_eh_triplet_Gamma_B diff --git a/src/Parquet/R_irred_Parquet_self_energy.f90 b/src/Parquet/R_irred_Parquet_self_energy.f90 new file mode 100644 index 0000000..e385ff2 --- /dev/null +++ b/src/Parquet/R_irred_Parquet_self_energy.f90 @@ -0,0 +1,94 @@ +subroutine R_irred_Parquet_self_energy(nOrb,nC,nO,nV,nR,e,EcGM,SigC,Z) + +! Compute correlation part of the self-energy with only irreducible vertices contribution + implicit none + include 'parameters.h' + +! Input variables + integer,intent(in) :: nOrb,nC,nO,nV,nR + double precision,intent(in) :: e(nOrb) +! Local variables + integer :: p,i,j,a,b + double precision :: D2p1h,D2h1p +! Output variables + double precision,intent(out) :: EcGM + double precision,intent(out) :: SigC(nOrb) + double precision,intent(out) :: Z(nOrb) + +!----------------------------! +! Static Parquet self-energy ! +!----------------------------! + SigC(:) = 0d0 + ! 2h1p part of the correlation self-energy + do p=nC+1,nOrb-nR + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nOrb-nR + + D2h1p = e(p) + e(a) - e(i) - e(j) + SigC(p) = SigC(p) !+ 2d0*rho(p,i,m)**2*(1d0-exp(-2d0*s*Dpim*Dpim))/Dpim + + end do + end do + end do + end do + ! 2p1h part of the correlation self-energy + do p=nC+1,nOrb-nR + do i=nC+1,nO + do a=nO+1,nOrb-nR + do b=nO+1,nOrb-nR + + D2p1h = e(p) + e(i) - e(a) - e(b) + SigC(p) = SigC(p) !+ 2d0*rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam + + end do + end do + end do + end do +!------------------------! +! Renormalization factor ! +!------------------------! + Z(:) = 0d0 + ! 2h1p part of the renormlization factor + do p=nC+1,nOrb-nR + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nOrb-nR + + D2h1p = e(p) + e(a) - e(i) - e(j) + Z(p) = Z(p) + + end do + end do + end do + end do + ! 2p1h part of the renormlization factor + do p=nC+1,nOrb-nR + do i=nC+1,nO + do a=nO+1,nOrb-nR + do b=nO+1,nOrb-nR + + D2p1h = e(p) + e(i) - e(a) - e(b) + Z(p) = Z(p) + + end do + end do + end do + end do + + Z(:) = 1d0/(1d0 - Z(:)) + +!-------------------------------------! +! Galitskii-Migdal correlation energy ! +!-------------------------------------! + + EcGM = 0d0 + ! do i=nC+1,nO + ! do a=nO+1,nOrb-nR + ! do m=1,nS + + ! end do + ! end do + ! end do + +end subroutine R_irred_Parquet_self_energy diff --git a/src/Parquet/R_pp_singlet_Gam.f90 b/src/Parquet/R_pp_singlet_Gam.f90 new file mode 100644 index 0000000..9bc0e7b --- /dev/null +++ b/src/Parquet/R_pp_singlet_Gam.f90 @@ -0,0 +1,225 @@ +subroutine R_pp_singlet_Gamma(nOrb,nC,nR,nS,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho,pp_sing_Gam) + +! Compute irreducible vertex in the triplet pp channel + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nR,nS + double precision,intent(in) :: eh_sing_Om(nS) + double precision,intent(in) :: eh_sing_rho(nOrb,nOrb,nS) + double precision,intent(in) :: eh_trip_Om(nS) + double precision,intent(in) :: eh_trip_rho(nOrb,nOrb,nS) + +! Local variables + integer :: p,q,r,s + integer :: n + double precision,external :: Kronecker_delta + +! Output variables + double precision, intent(out) :: pp_sing_Gam(nOrb,nOrb,nOrb,nOrb) + +! Initialization + pp_sing_Gam(:,:,:,:) = 0d0 + +! !$OMP PARALLEL DEFAULT(NONE) & +! !$OMP PRIVATE(i, j, ij, k, l, kl, n) & +! !$OMP SHARED(nC, nOrb, nO, nS, pp_sing_Gam_D, eh_sing_rho, eh_sing_Om, eh_trip_rho, eh_trip_Om) +! !$OMP DO COLLAPSE(2) + + do s = nC+1, nOrb-nR + do r = nC+1, nOrb-nR + do q = nC+1, nOrb-nR + do p = nC+1, nOrb-nR + + do n=1,nS + pp_sing_Gam(p,q,r,s) = pp_sing_Gam(p,q,r,s) & + - eh_sing_rho(p,r,n)*eh_sing_rho(q,s,n)/eh_sing_Om(n) & + + 3d0 * eh_trip_rho(p,r,n)*eh_trip_rho(q,s,n)/eh_trip_Om(n) & + - eh_sing_rho(p,s,n)*eh_sing_rho(q,r,n)/eh_sing_Om(n) & + + 3d0 * eh_trip_rho(p,s,n)*eh_trip_rho(q,r,n)/eh_trip_Om(n) + end do + + pp_sing_Gam(p,q,r,s) = pp_sing_Gam(p,q,r,s)/sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(r,s))) + + end do + end do + end do + end do +! !$OMP END DO +! !$OMP END PARALLEL + +end subroutine R_pp_singlet_Gamma + +subroutine R_pp_singlet_Gamma_D(nOrb,nC,nO,nV,nR,nS,nOOs,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho, pp_sing_Gam_D) + +! Compute irreducible vertex in the triplet pp channel + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nOOs + double precision,intent(in) :: eh_sing_Om(nS) + double precision,intent(in) :: eh_sing_rho(nOrb,nOrb,nS) + double precision,intent(in) :: eh_trip_Om(nS) + double precision,intent(in) :: eh_trip_rho(nOrb,nOrb,nS) + +! Local variables + integer :: i,j,k,l + integer :: ij,kl + integer :: n + double precision,external :: Kronecker_delta + +! Output variables + double precision, intent(out) :: pp_sing_Gam_D(nOOs,nOOs) + +! Initialization + pp_sing_Gam_D(:,:) = 0d0 + +! !$OMP PARALLEL DEFAULT(NONE) & +! !$OMP PRIVATE(i, j, ij, k, l, kl, n) & +! !$OMP SHARED(nC, nOrb, nO, nS, pp_sing_Gam_D, eh_sing_rho, eh_sing_Om, eh_trip_rho, eh_trip_Om) +! !$OMP DO COLLAPSE(2) + + ij = 0 + do i=nC+1,nO + do j=i,nO + ij = ij + 1 + + kl = 0 + do k=nC+1,nO + do l=k,nO + kl = kl +1 + + do n=1,nS + pp_sing_Gam_D(ij,kl) = pp_sing_Gam_D(ij,kl) & + - eh_sing_rho(i,k,n)*eh_sing_rho(j,l,n)/eh_sing_Om(n) & + + 3d0 * eh_trip_rho(i,k,n)*eh_trip_rho(j,l,n)/eh_trip_Om(n) & + - eh_sing_rho(i,l,n)*eh_sing_rho(j,k,n)/eh_sing_Om(n) & + + 3d0 * eh_trip_rho(i,l,n)*eh_trip_rho(j,k,n)/eh_trip_Om(n) + end do + + pp_sing_Gam_D(ij,kl) = pp_sing_Gam_D(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) + + end do + end do + end do + end do +! !$OMP END DO +! !$OMP END PARALLEL + +end subroutine + +subroutine R_pp_singlet_Gamma_C(nOrb,nC,nO,nV,nR,nS,nVVs,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho, pp_sing_Gam_C) + +! Compute irreducible vertex in the triplet pp channel + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nVVs + double precision,intent(in) :: eh_sing_Om(nS) + double precision,intent(in) :: eh_sing_rho(nOrb,nOrb,nS) + double precision,intent(in) :: eh_trip_Om(nS) + double precision,intent(in) :: eh_trip_rho(nOrb,nOrb,nS) + +! Local variables + integer :: a,b,c,d + integer :: ab,cd + integer :: n + double precision,external :: Kronecker_delta + +! Output variables + double precision, intent(out) :: pp_sing_Gam_C(nVVs,nVVs) + +! Initialization + pp_sing_Gam_C(:,:) = 0d0 + +! !$OMP PARALLEL DEFAULT(NONE) & +! !$OMP PRIVATE(a, b, ab, c, d, cd, n) & +! !$OMP SHARED(nC, nOrb, nO, nS, pp_sing_Gam_C, eh_sing_rho, eh_sing_Om, eh_trip_rho, eh_trip_Om) +! !$OMP DO COLLAPSE(2) + + ab = 0 + do a=nO+1,nOrb - nR + do b=a,nOrb - nR + ab = ab + 1 + + cd = 0 + do c=nO+1,nOrb - nR + do d=c,nOrb - nR + cd = cd +1 + + do n=1,nS + pp_sing_Gam_C(ab,cd) = pp_sing_Gam_C(ab,cd) & + - eh_sing_rho(a,c,n)*eh_sing_rho(b,d,n)/eh_sing_Om(n) & + + 3d0 * eh_trip_rho(a,c,n)*eh_trip_rho(b,d,n)/eh_trip_Om(n) & + - eh_sing_rho(a,d,n)*eh_sing_rho(b,c,n)/eh_sing_Om(n) & + + 3d0 * eh_trip_rho(a,d,n)*eh_trip_rho(b,c,n)/eh_trip_Om(n) + end do + + pp_sing_Gam_C(ab,cd) = pp_sing_Gam_C(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) + + end do + end do + end do + end do +! !$OMP END DO +! !$OMP END PARALLEL + +end subroutine + +subroutine R_pp_singlet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho,pp_sing_Gam_B) + +! Compute irreducible vertex in the triplet pp channel + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nOOs,nVVs + double precision,intent(in) :: eh_sing_Om(nS) + double precision,intent(in) :: eh_sing_rho(nOrb,nOrb,nS) + double precision,intent(in) :: eh_trip_Om(nS) + double precision,intent(in) :: eh_trip_rho(nOrb,nOrb,nS) + +! Local variables + integer :: a,b,i,j + integer :: ab,ij + integer :: n + double precision,external :: Kronecker_delta + +! Output variables + double precision, intent(out) :: pp_sing_Gam_B(nVVs,nOOs) + +! Initialization + pp_sing_Gam_B(:,:) = 0d0 + +! !$OMP PARALLEL DEFAULT(NONE) & +! !$OMP PRIVATE(a, b, ab, i, j, ij, n) & +! !$OMP SHARED(nC, nOrb, nO, nS, pp_sing_Gam_B, eh_sing_rho, eh_sing_Om, eh_trip_rho, eh_trip_Om) +! !$OMP DO COLLAPSE(2) + + ab = 0 + do a=nO+1,nOrb - nR + do b=a,nOrb - nR + ab = ab + 1 + + ij = 0 + do i=nC+1,nO + do j=i,nO + ij = ij +1 + + do n=1,nS + pp_sing_Gam_B(ab,ij) = pp_sing_Gam_B(ab,ij) & + - eh_sing_rho(a,i,n)*eh_sing_rho(b,j,n)/eh_sing_Om(n) & + + 3d0 * eh_trip_rho(a,i,n)*eh_trip_rho(b,j,n)/eh_trip_Om(n) & + - eh_sing_rho(a,j,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) & + + 3d0 * eh_trip_rho(a,j,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n) + end do + + pp_sing_Gam_B(ab,ij) = pp_sing_Gam_B(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) + + end do + end do + end do + end do +! !$OMP END DO +! !$OMP END PARALLEL + +end subroutine diff --git a/src/Parquet/R_pp_triplet_Gam.f90 b/src/Parquet/R_pp_triplet_Gam.f90 new file mode 100644 index 0000000..7c3efbb --- /dev/null +++ b/src/Parquet/R_pp_triplet_Gam.f90 @@ -0,0 +1,217 @@ +subroutine R_pp_triplet_Gamma(nOrb,nC,nR,nS,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho, pp_trip_Gam) + +! Compute irreducible vertex in the triplet pp channel + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nR,nS + double precision,intent(in) :: eh_sing_Om(nS) + double precision,intent(in) :: eh_sing_rho(nOrb,nOrb,nS) + double precision,intent(in) :: eh_trip_Om(nS) + double precision,intent(in) :: eh_trip_rho(nOrb,nOrb,nS) + +! Local variables + integer :: p,q,r,s + integer :: n + double precision,external :: Kronecker_delta + +! Output variables + double precision, intent(out) :: pp_trip_Gam(nOrb,nOrb,nOrb,nOrb) + +! Initialization + pp_trip_Gam(:,:,:,:) = 0d0 + +! !$OMP PARALLEL DEFAULT(NONE) & +! !$OMP PRIVATE(i, j, ij, k, l, kl, n) & +! !$OMP SHARED(nC, nOrb, nO, nS, pp_trip_Gam_D, eh_sing_rho, eh_sing_Om, eh_trip_rho, eh_trip_Om) +! !$OMP DO COLLAPSE(2) + + do s = nC+1, nOrb-nR + do r = nC+1, nOrb-nR + do q = nC+1, nOrb-nR + do p = nC+1, nOrb-nR + + do n=1,nS + pp_trip_Gam(p,q,r,s) = pp_trip_Gam(p,q,r,s) & + - eh_sing_rho(p,r,n)*eh_sing_rho(q,s,n)/eh_sing_Om(n) & + - eh_trip_rho(p,r,n)*eh_trip_rho(q,s,n)/eh_trip_Om(n) & + + eh_sing_rho(p,s,n)*eh_sing_rho(q,r,n)/eh_sing_Om(n) & + + eh_trip_rho(p,s,n)*eh_trip_rho(q,r,n)/eh_trip_Om(n) + end do + + end do + end do + end do + end do +! !$OMP END DO +! !$OMP END PARALLEL + +end subroutine + +subroutine R_pp_triplet_Gamma_D(nOrb,nC,nO,nV,nR,nS,nOOt,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho, pp_trip_Gam_D) + +! Compute irreducible vertex in the triplet pp channel + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nOOt + double precision,intent(in) :: eh_sing_Om(nS) + double precision,intent(in) :: eh_sing_rho(nOrb,nOrb,nS) + double precision,intent(in) :: eh_trip_Om(nS) + double precision,intent(in) :: eh_trip_rho(nOrb,nOrb,nS) + +! Local variables + integer :: i,j,k,l + integer :: ij,kl + integer :: n + double precision,external :: Kronecker_delta + +! Output variables + double precision, intent(out) :: pp_trip_Gam_D(nOOt,nOOt) + +! Initialization + pp_trip_Gam_D(:,:) = 0d0 + +! !$OMP PARALLEL DEFAULT(NONE) & +! !$OMP PRIVATE(i, j, ij, k, l, kl, n) & +! !$OMP SHARED(nC, nOrb, nO, nS, pp_trip_Gam_D, eh_sing_rho, eh_sing_Om, eh_trip_rho, eh_trip_Om) +! !$OMP DO COLLAPSE(2) + + ij = 0 + do i=nC+1,nO + do j=i+1,nO + ij = ij + 1 + + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl +1 + + do n=1,nS + pp_trip_Gam_D(ij,kl) = pp_trip_Gam_D(ij,kl) & + - eh_sing_rho(i,k,n)*eh_sing_rho(j,l,n)/eh_sing_Om(n) & + - eh_trip_rho(i,k,n)*eh_trip_rho(j,l,n)/eh_trip_Om(n) & + + eh_sing_rho(i,l,n)*eh_sing_rho(j,k,n)/eh_sing_Om(n) & + + eh_trip_rho(i,l,n)*eh_trip_rho(j,k,n)/eh_trip_Om(n) + end do + + end do + end do + end do + end do +! !$OMP END DO +! !$OMP END PARALLEL + +end subroutine + +subroutine R_pp_triplet_Gamma_C(nOrb,nC,nO,nV,nR,nS,nVVt,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho, pp_trip_Gam_C) + +! Compute irreducible vertex in the triplet pp channel + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nVVt + double precision,intent(in) :: eh_sing_Om(nS) + double precision,intent(in) :: eh_sing_rho(nOrb,nOrb,nS) + double precision,intent(in) :: eh_trip_Om(nS) + double precision,intent(in) :: eh_trip_rho(nOrb,nOrb,nS) + +! Local variables + integer :: a,b,c,d + integer :: ab,cd + integer :: n + double precision,external :: Kronecker_delta + +! Output variables + double precision, intent(out) :: pp_trip_Gam_C(nVVt,nVVt) + +! Initialization + pp_trip_Gam_C(:,:) = 0d0 + +! !$OMP PARALLEL DEFAULT(NONE) & +! !$OMP PRIVATE(a, b, ab, c, d, cd, n) & +! !$OMP SHARED(nC, nOrb, nO, nS, pp_trip_Gam_C, eh_sing_rho, eh_sing_Om, eh_trip_rho, eh_trip_Om) +! !$OMP DO COLLAPSE(2) + + ab = 0 + do a=nO+1,nOrb - nR + do b=a+1,nOrb - nR + ab = ab + 1 + + cd = 0 + do c=nO+1,nOrb - nR + do d=c+1,nOrb - nR + cd = cd +1 + + do n=1,nS + pp_trip_Gam_C(ab,cd) = pp_trip_Gam_C(ab,cd) & + - eh_sing_rho(a,c,n)*eh_sing_rho(b,d,n)/eh_sing_Om(n) & + - eh_trip_rho(a,c,n)*eh_trip_rho(b,d,n)/eh_trip_Om(n) & + + eh_sing_rho(a,d,n)*eh_sing_rho(b,c,n)/eh_sing_Om(n) & + + eh_trip_rho(a,d,n)*eh_trip_rho(b,c,n)/eh_trip_Om(n) + end do + + end do + end do + end do + end do +! !$OMP END DO +! !$OMP END PARALLEL + +end subroutine + +subroutine R_pp_triplet_Gamma_B(nOrb,nC,nO,nV,nR,nS,nOOt,nVVt,eh_sing_Om,eh_sing_rho,eh_trip_Om,eh_trip_rho,pp_trip_Gam_B) + +! Compute irreducible vertex in the triplet pp channel + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nO,nV,nR,nS,nOOt,nVVt + double precision,intent(in) :: eh_sing_Om(nS) + double precision,intent(in) :: eh_sing_rho(nOrb,nOrb,nS) + double precision,intent(in) :: eh_trip_Om(nS) + double precision,intent(in) :: eh_trip_rho(nOrb,nOrb,nS) + +! Local variables + integer :: a,b,i,j + integer :: ab,ij + integer :: n + double precision,external :: Kronecker_delta + +! Output variables + double precision, intent(out) :: pp_trip_Gam_B(nVVt,nOOt) + +! Initialization + pp_trip_Gam_B(:,:) = 0d0 + +! !$OMP PARALLEL DEFAULT(NONE) & +! !$OMP PRIVATE(a, b, ab, i, j, ij, n) & +! !$OMP SHARED(nC, nOrb, nO, nS, pp_trip_Gam_B, eh_sing_rho, eh_sing_Om, eh_trip_rho, eh_trip_Om) +! !$OMP DO COLLAPSE(2) + + ab = 0 + do a=nO+1,nOrb - nR + do b=a+1,nOrb - nR + ab = ab + 1 + + ij = 0 + do i=nC+1,nO + do j=i+1,nO + ij = ij +1 + + do n=1,nS + pp_trip_Gam_B(ab,ij) = pp_trip_Gam_B(ab,ij) & + - eh_sing_rho(a,i,n)*eh_sing_rho(b,j,n)/eh_sing_Om(n) & + - eh_trip_rho(a,i,n)*eh_trip_rho(b,j,n)/eh_trip_Om(n) & + + eh_sing_rho(a,j,n)*eh_sing_rho(b,i,n)/eh_sing_Om(n) & + + eh_trip_rho(a,j,n)*eh_trip_rho(b,i,n)/eh_trip_Om(n) + end do + + end do + end do + end do + end do +! !$OMP END DO +! !$OMP END PARALLEL + +end subroutine diff --git a/src/Parquet/R_screened_integrals.f90 b/src/Parquet/R_screened_integrals.f90 new file mode 100644 index 0000000..97c24c6 --- /dev/null +++ b/src/Parquet/R_screened_integrals.f90 @@ -0,0 +1,293 @@ +subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,XpY,rho) + +! Compute excitation densities + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nO,nR,nS + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: XpY(nS,nS) + +! Local variables + integer :: ia,jb,p,q,j,b + +! Output variables + double precision,intent(out) :: rho(nOrb,nOrb,nS) + + rho(:,:,:) = 0d0 + !$OMP PARALLEL & + !$OMP SHARED(nC,nOrb,nR,nO,nS,rho,ERI,XpY) & + !$OMP PRIVATE(q,p,jb,ia) & + !$OMP DEFAULT(NONE) + !$OMP DO + do q=nC+1,nOrb-nR + do p=nC+1,nOrb-nR + jb = 0 + do j=nC+1,nO + do b=nO+1,nOrb-nR + jb = jb + 1 + do ia=1,nS + rho(p,q,ia) = rho(p,q,ia) + (2d0*ERI(p,j,q,b) - ERI(p,j,b,q))*XpY(ia,jb) + end do + end do + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + +end subroutine R_eh_singlet_screened_integral + +subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,XpY,rho) + +! Compute excitation densities + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nO,nR,nS + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: XpY(nS,nS) + +! Local variables + integer :: ia,jb,p,q,j,b + +! Output variables + double precision,intent(out) :: rho(nOrb,nOrb,nS) + + rho(:,:,:) = 0d0 + !$OMP PARALLEL & + !$OMP SHARED(nC,nOrb,nR,nO,nS,rho,ERI,XpY) & + !$OMP PRIVATE(q,p,jb,ia) & + !$OMP DEFAULT(NONE) + !$OMP DO + do q=nC+1,nOrb-nR + do p=nC+1,nOrb-nR + jb = 0 + do j=nC+1,nO + do b=nO+1,nOrb-nR + jb = jb + 1 + do ia=1,nS + rho(p,q,ia) = rho(p,q,ia) - ERI(p,j,b,q)*XpY(ia,jb) + end do + end do + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + +end subroutine R_eh_triplet_screened_integral + + +subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_sing_Gam,X1,Y1,rho1,X2,Y2,rho2) + +! Compute excitation densities in the singlet pp channel + + implicit none + +! Input variables + + + integer,intent(in) :: nOrb,nC,nO,nV,nR + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: pp_sing_Gam(nOrb,nOrb,nOrb,nOrb) + integer,intent(in) :: nOO + integer,intent(in) :: nVV + double precision,intent(in) :: X1(nVV,nVV) + double precision,intent(in) :: Y1(nOO,nVV) + double precision,intent(in) :: X2(nVV,nOO) + double precision,intent(in) :: Y2(nOO,nOO) + +! Local variables + + integer :: i,j,k,l + integer :: a,b,c,d + integer :: p,q + integer :: ab,cd,ij,kl + double precision,external :: Kronecker_delta + +! Output variables + + double precision,intent(out) :: rho1(nOrb,nOrb,nVV) + double precision,intent(out) :: rho2(nOrb,nOrb,nOO) + + integer :: dim_1, dim_2 + +! Initialization + + rho1(:,:,:) = 0d0 + rho2(:,:,:) = 0d0 + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(p, q, a, b, ab, c, d, cd, i, j, ij, k, l, kl) & + !$OMP SHARED(nC, nOrb, nR, nO, rho1, rho2, ERI, pp_sing_Gam, X1, Y1, X2, Y2) + !$OMP DO COLLAPSE(2) + do q=nC+1,nOrb-nR + do p=nC+1,nOrb-nR + + ab = 0 + do a=nO+1,nOrb-nR + do b=a,nOrb-nR + ab = ab + 1 + + cd = 0 + do c=nO+1,nOrb-nR + do d=c,nOrb-nR + cd = cd + 1 + + rho1(p,q,ab) = rho1(p,q,ab) & + + (ERI(p,q,c,d) + ERI(p,q,d,c) + 0d0*pp_sing_Gam(p,q,c,d))*X1(cd,ab)/ & + sqrt(1d0 + Kronecker_delta(c,d)) + end do + end do + + kl = 0 + do k=nC+1,nO + do l=k,nO + kl = kl + 1 + rho1(p,q,ab) = rho1(p,q,ab) & + + (ERI(p,q,k,l) + ERI(p,q,l,k) + 0d0*pp_sing_Gam(p,q,k,l))*Y1(kl,ab)/ & + sqrt(1d0 + Kronecker_delta(k,l)) + end do + end do + + end do + end do + + ij = 0 + do i=nC+1,nO + do j=i,nO + ij = ij + 1 + cd = 0 + do c=nO+1,nOrb-nR + do d=c,nOrb-nR + cd = cd + 1 + rho2(p,q,ij) = rho2(p,q,ij) & + + (ERI(p,q,c,d) + ERI(p,q,d,c) + 0d0*pp_sing_Gam(p,q,c,d))*X2(cd,ij)/ & + sqrt(1d0 + Kronecker_delta(c,d)) + end do + end do + + kl = 0 + do k=nC+1,nO + do l=k,nO + kl = kl + 1 + rho2(p,q,ij) = rho2(p,q,ij) & + + (ERI(p,q,k,l) + ERI(p,q,l,k) + 0d0*pp_sing_Gam(p,q,k,l))*Y2(kl,ij)/ & + sqrt(1d0 + Kronecker_delta(k,l)) + end do + end do + end do + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + +end subroutine R_pp_singlet_screened_integral + + + + +subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,pp_trip_Gam,X1,Y1,rho1,X2,Y2,rho2) + +! Compute excitation densities in the triplet pp channel + implicit none + +! Input variables + integer,intent(in) :: nOrb,nC,nO,nV,nR + integer,intent(in) :: nOO + integer,intent(in) :: nVV + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: pp_trip_Gam(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: X1(nVV,nVV) + double precision,intent(in) :: Y1(nOO,nVV) + double precision,intent(in) :: X2(nVV,nOO) + double precision,intent(in) :: Y2(nOO,nOO) + +! Local variables + integer :: i,j,k,l + integer :: a,b,c,d + integer :: p,q + integer :: ab,cd,ij,kl + double precision,external :: Kronecker_delta + +! Output variables + double precision,intent(out) :: rho1(nOrb,nOrb,nVV) + double precision,intent(out) :: rho2(nOrb,nOrb,nOO) + integer :: dim_1, dim_2 + +! Initialization + rho1(:,:,:) = 0d0 + rho2(:,:,:) = 0d0 + + dim_1 = (nOrb - nO) * (nOrb - nO - 1) / 2 + dim_2 = nO * (nO - 1) / 2 + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(p, q, a, b, ab, c, d, cd, i, j, ij, k, l, kl) & + !$OMP SHARED(nC, nOrb, nR, nO, rho1, rho2, ERI, pp_trip_Gam, X1, Y1, X2, Y2) + !$OMP DO COLLAPSE(2) + do q = nC+1, nOrb-nR + do p = nC+1, nOrb-nR + ab = 0 + + do a = nO+1, nOrb-nR + do b = a+1, nOrb-nR + ab = ab + 1 + + cd = 0 + do c = nO+1, nOrb-nR + do d = c+1, nOrb-nR + cd = cd + 1 + + rho1(p,q,ab) = rho1(p,q,ab) & + + (ERI(p,q,c,d) - ERI(p,q,d,c) + 0d0*pp_trip_Gam(p,q,c,d))*X1(cd,ab) + end do ! d + end do ! c + + kl = 0 + do k = nC+1, nO + do l = k+1, nO + + kl = kl + 1 + + rho1(p,q,ab) = rho1(p,q,ab) & + + (ERI(p,q,k,l) - ERI(p,q,l,k) + 0d0*pp_trip_Gam(p,q,k,l))*Y1(kl,ab) + end do ! l + end do ! k + end do ! b + end do ! a + + ij = 0 + do i = nC+1, nO + do j = i+1, nO + ij = ij + 1 + + cd = 0 + do c = nO+1, nOrb-nR + do d = c+1, nOrb-nR + cd = cd + 1 + + rho2(p,q,ij) = rho2(p,q,ij) & + + (ERI(p,q,c,d) - ERI(p,q,d,c) + 0d0*pp_trip_Gam(p,q,c,d))*X2(cd,ij) + end do ! d + end do ! c + + kl = 0 + do k = nC+1, nO + do l = k+1, nO + kl = kl + 1 + + rho2(p,q,ij) = rho2(p,q,ij) & + + (ERI(p,q,k,l) - ERI(p,q,l,k) + 0d0*pp_trip_Gam(p,q,k,l))*Y2(kl,ij) + end do ! l + end do ! k + end do ! j + end do ! i + end do ! p + end do ! q + !$OMP END DO + !$OMP END PARALLEL + +end subroutine R_pp_triplet_screened_integral diff --git a/src/QuAcK/GQuAcK.f90 b/src/QuAcK/GQuAcK.f90 index 0ee9e95..30fa842 100644 --- a/src/QuAcK/GQuAcK.f90 +++ b/src/QuAcK/GQuAcK.f90 @@ -341,7 +341,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 wall_time(end_Parquet) t_Parquet = end_Parquet - start_Parquet diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 50d3a64..a9477ca 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -78,6 +78,10 @@ program QuAcK logical :: restart_hfb double precision :: temperature,sigma + + integer :: max_it_macro,max_it_micro + double precision :: conv_one_body,conv_two_body + character(len=256) :: working_dir ! Check if the right number of arguments is provided @@ -141,7 +145,8 @@ program QuAcK maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, & doACFDT,exchange_kernel,doXBS, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA, & - temperature,sigma,chem_pot_hf,restart_hfb) + temperature,sigma,chem_pot_hf,restart_hfb, & + max_it_macro,conv_one_body,max_it_micro,conv_two_body) !------------------! ! Hardware ! @@ -257,7 +262,8 @@ program QuAcK guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA, & 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) + dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS, & + max_it_macro,conv_one_body,max_it_micro,conv_two_body) endif endif diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index b63f124..63dc813 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -1,13 +1,14 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, & doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW, & - doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh,doParquet, & + doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh,doParquet, & nNuc,nBas,nOrb,nC,nO,nV,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,singlet,triplet,TDA, & 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) + dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS, & + max_it_macro,conv_one_body,max_it_micro,conv_two_body) ! Restricted branch of QuAcK @@ -81,6 +82,9 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2, 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,doCI,doRPA,doGF,doGW,doGT @@ -370,11 +374,13 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2, if(doParquet) then call wall_time(start_Parquet) - + call RParquet(max_it_macro,conv_one_body,max_it_micro,conv_two_body, & + nOrb,nC,nO,nV,nR,nS, & + eHF,ERI_MO) call wall_time(end_Parquet) t_Parquet = end_Parquet - start_Parquet - write(*,'(A65,1X,F9.3,A8)') 'Total wall time for Parquet module = ',t_Parquet,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for Parquet module = ', t_Parquet, ' seconds' write(*,*) end if diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 65b9354..27ca5ee 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -8,7 +8,8 @@ subroutine read_options(working_dir, maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, & doACFDT,exchange_kernel,doXBS, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA, & - temperature,sigma,chem_pot_hf,restart_hfb) + temperature,sigma,chem_pot_hf,restart_hfb, & + max_it_macro,conv_one_body,max_it_micro,conv_two_body) ! Read desired methods @@ -78,6 +79,9 @@ subroutine read_options(working_dir, double precision,intent(out) :: temperature double precision,intent(out) :: sigma + integer,intent(out) :: max_it_macro,max_it_micro + double precision,intent(out) :: conv_one_body,conv_two_body + ! Local variables character(len=1) :: ans1,ans2,ans3,ans4,ans5 @@ -235,7 +239,17 @@ subroutine read_options(working_dir, if(ans1 == 'T') chem_pot_hf = .true. if(ans2 == 'T') restart_hfb = .true. + + ! Options for Parquet module + max_it_macro = 1 + conv_one_body = 0.01 + max_it_micro = 1 + conv_two_body = 0.01 + + read(1,*) + read(1,*) max_it_macro,conv_one_body,max_it_micro,conv_two_body + endif ! Close file with options