PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-04-23 16:54:54 +02:00
Code Releases Activity

update on Parquet

Browse Source
This commit is contained in:
Antoine Marie 2025-03-17 14:56:56 +01:00
parent 10e8b9c885
commit 14318a2178
15 changed files with 1714 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
input/options.default
View File

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

4
src/LR/ppRLR_C.f90
View File

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

4
src/LR/ppRLR_D.f90
View File

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

22
src/Parquet/README.md Normal file
View File

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

511
src/Parquet/RParquet.f90 Normal file
View File

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

155
src/Parquet/R_eh_singlet_Gam.f90 Normal file
View File

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

155
src/Parquet/R_eh_triplet_Gam.f90 Normal file
View File

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

94
src/Parquet/R_irred_Parquet_self_energy.f90 Normal file
View File

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

225
src/Parquet/R_pp_singlet_Gam.f90 Normal file
View File

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

217
src/Parquet/R_pp_triplet_Gam.f90 Normal file
View File

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

293
src/Parquet/R_screened_integrals.f90 Normal file
View File

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

4
src/QuAcK/GQuAcK.f90
View File

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

10
src/QuAcK/QuAcK.f90
View File

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

14
src/QuAcK/RQuAcK.f90
View File

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

16
src/QuAcK/read_options.f90
View File

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