PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-04-25 01:34:57 +02:00
Code Releases Activity

starting work on G Parquet

Browse Source
This commit is contained in:
Antoine Marie 2025-03-22 17:20:56 +01:00
parent 45fc92b764
commit 36685fcc73
10 changed files with 314 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/AOtoMO/AOtoMO_ERI_RHF.f90
View File

@ -32,6 +32,7 @@ subroutine AOtoMO_ERI_RHF(nBas,nOrb,c,ERI_AO,ERI_MO)
, ERI_AO(1,1,1,1), nBas, c(1,1), nBas &
, 0.d0, a2(1,1,1,1), nBas*nBas*nBas)
call dgemm( 'T', 'N', nBas*nBas*nOrb, nOrb, nBas, 1.d0 &
, a2(1,1,1,1), nBas, c(1,1), nBas &
, 0.d0, a1(1,1,1,1), nBas*nBas*nOrb)
@ -50,5 +51,5 @@ subroutine AOtoMO_ERI_RHF(nBas,nOrb,c,ERI_AO,ERI_MO)
, 0.d0, ERI_MO(1,1,1,1), nOrb*nOrb*nOrb)
deallocate(a2)
end subroutine

5
src/GW/RG0W0.f90
View File

@ -1,5 +1,5 @@
subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF,eGW_out)
! Perform G0W0 calculation
@ -61,6 +61,7 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
double precision,allocatable :: eGWlin(:)
double precision,allocatable :: eGW(:)
double precision,intent(inout):: eGW_out(nOrb)
! Output variables
@ -171,6 +172,8 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
call print_RG0W0(nOrb,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM)
eGW_out(:) = eGW(:)
!---------------------------!
! Perform phBSE calculation !
!---------------------------!

7
src/GW/RGW.f90
View File

@ -1,7 +1,7 @@
subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_diis,doACFDT, &
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF, &
S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF,eGW)
! Restricted GW module
@ -70,6 +70,9 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_dii
logical :: doccG0W0,doccGW
double precision,intent(inout) :: eGW(nOrb)
!------------------------------------------------------------------------
! Perform G0W0 calculation
!------------------------------------------------------------------------
@ -78,7 +81,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_dii
call wall_time(start_GW)
call RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF,eGW)
call wall_time(end_GW)
t_GW = end_GW - start_GW

2
src/LR/print_excitation_energies.f90
View File

@ -14,7 +14,7 @@ subroutine print_excitation_energies(method,manifold,nS,Om)
! Local variables
integer,parameter :: maxS = 10
integer,parameter :: maxS = 25
integer :: m
write(*,*)

281
src/Parquet/GParquet.f90 Normal file
View File

@ -0,0 +1,281 @@
subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,ERI)
! Parquet approximation based on restricted orbitals
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_1b,max_it_2b
double precision,intent(in) :: conv_1b,conv_2b
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 :: ispin
integer :: n_it_1b,n_it_2b
double precision :: err_1b,err_2b
double precision :: err_eh, err_hh, err_ee
double precision :: start_t, end_t, t
double precision :: start_1b, end_1b, t_1b
double precision :: start_2b, end_2b, t_2b
integer :: nOO,nVV
double precision :: EcRPA
double precision,allocatable :: Aph(:,:), Bph(:,:)
double precision,allocatable :: XpY(:,:),XmY(:,:)
double precision,allocatable :: eh_Om(:), old_eh_Om(:)
double precision,allocatable :: Bpp(:,:), Cpp(:,:), Dpp(:,:)
double precision,allocatable :: X1(:,:),Y1(:,:)
double precision,allocatable :: ee_Om(:), old_ee_Om(:)
double precision,allocatable :: X2(:,:),Y2(:,:)
double precision,allocatable :: hh_Om(:), old_hh_Om(:)
double precision,allocatable :: eh_rho(:,:,:), ee_rho(:,:,:), hh_rho(:,:,:)
double precision,allocatable :: eh_Gam_A(:,:),eh_Gam_B(:,:)
double precision,allocatable :: pp_Gam_B(:,:),pp_Gam_C(:,:),pp_Gam_D(:,:)
double precision,allocatable :: eh_Gam(:,:,:,:),pp_Gam(:,:,:,:)
double precision,allocatable :: eQPlin(:),eQP(:),eOld(:)
double precision,allocatable :: SigC(:)
double precision,allocatable :: Z(:)
double precision :: EcGM
! Output variables
! None
nOO = nO*(nO - 1)/2
nVV = nV*(nV - 1)/2
allocate(eQP(nOrb),eOld(nOrb))
write(*,*)
write(*,*)'***********************************'
write(*,*)'* Generalized Parquet Calculation *'
write(*,*)'***********************************'
write(*,*)
! Print parameters
write(*,*)'---------------------------------------------------------------'
write(*,*)' Parquet parameters for one-body and two-body self-consistency '
write(*,*)'---------------------------------------------------------------'
write(*,'(1X,A50,1X,I5)') 'Maximum number for one-body self-consistency:', max_it_1b
write(*,'(1X,A50,1X,E10.5)') 'Convergence threshold for one-body energies:', conv_1b
write(*,*)'---------------------------------------------------------------'
write(*,'(1X,A50,1X,I5)') 'Maximum number for two-body self-consistency:', max_it_2b
write(*,'(1X,A50,1X,E10.5)') 'Convergence threshold for two-body energies:', conv_2b
write(*,*)'---------------------------------------------------------------'
write(*,*)
if(linearize) then
write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
write(*,*)
else
write(*,*) ' *** Quasiparticle energies obtained by root search *** '
write(*,*)
endif
! Memory allocation
allocate(old_eh_Om(nS),old_ee_Om(nVV),old_hh_Om(nOO))
allocate(eh_rho(nOrb,nOrb,nS),ee_rho(nOrb,nOrb,nVV),hh_rho(nOrb,nOrb,nOO))
! Initialization
n_it_1b = 0
err_1b = 1d0
n_it_2b = 0
err_2b = 1d0
eQP(:) = eHF(:)
eOld(:) = eHF(:)
eh_rho(:,:,:) = 0d0
ee_rho(:,:,:) = 0d0
hh_rho(:,:,:) = 0d0
old_eh_Om(:) = 1d0
old_ee_Om(:) = 1d0
old_hh_Om(:) = 1d0
!-----------------------------------------!
! Main loop for one-body self-consistency !
!-----------------------------------------!
do while(err_1b > conv_1b .and. n_it_1b < max_it_1b)
n_it_1b = n_it_1b + 1
call wall_time(start_1b)
write(*,*)
write(*,*)'====================================='
write(*,'(1X,A30,1X,I4)') 'One-body iteration #',n_it_1b
write(*,*)'====================================='
write(*,*)
!-----------------------------------------!
! Main loop for two-body self-consistency !
!-----------------------------------------!
do while(err_2b > conv_2b .and. n_it_2b < max_it_2b)
n_it_2b = n_it_2b + 1
call wall_time(start_2b)
write(*,*)' ***********************************'
write(*,'(1X,A30,1X,I4)') 'Two-body iteration #',n_it_2b
write(*,*)' ***********************************'
write(*,*)
!--------------------------------!
! Compute effective interactions !
!--------------------------------!
! Memory allocation
allocate(eh_Gam(nOrb,nOrb,nOrb,nOrb))
allocate(pp_Gam(nOrb,nOrb,nOrb,nOrb))
! Build eh effective interaction
write(*,*) 'Computing eh effective interaction...'
call wall_time(start_t)
!call R_eh_Gamma(nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
! old_eh_Om,eh_rho,old_ee_Om,ee_rho,old_hh_Om,hh_rho, &
! eh_trip_Gam)
call wall_time(end_t)
t = end_t - start_t
write(*,'(A50,1X,F9.3,A8)') 'Wall time for eh Gamma =',t,' seconds'
write(*,*)
! Build singlet pp effective interaction
write(*,*) 'Computing pp effective interaction...'
call wall_time(start_t)
!call R_pp_Gamma(nOrb,nC,nR,nS,old_eh_Om,eh_rho,pp_Gam)
call wall_time(end_t)
t = end_t - start_t
write(*,'(A50,1X,F9.3,A8)') 'Wall time for pp Gamma =',t,' seconds'
write(*,*)
call wall_time(end_2b)
t_2b = end_2b - start_2b
write(*,'(A50,1X,F9.3,A8)') 'Wall time for two-body iteration =',t_2b,' seconds'
write(*,*)
end do
!---------------------------------------------!
! End main loop for two-body self-consistency !
!---------------------------------------------!
! Did it actually converge?
if(n_it_2b == max_it_2b) 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','1h1p',nS,old_eh_Om)
call print_excitation_energies('ppBSE@Parquet','2p',nVV,old_ee_Om)
call print_excitation_energies('ppBSE@Parquet','2h',nOO,old_hh_Om)
end if
allocate(eQPlin(nOrb),Z(nOrb),SigC(nOrb))
write(*,*) 'Building self-energy'
call wall_time(start_t)
!call G_irred_Parquet_self_energy(nOrb,nC,nO,nV,nR,eOld,EcGM,SigC,Z)
call wall_time(end_t)
t = end_t - start_t
write(*,'(A50,1X,F9.3,A8)') 'Wall time for self energy =',t,' seconds'
write(*,*)
eQPlin(:) = eHF(:) !+ Z(:)*SigC(:)
! Solve the quasi-particle equation
if(linearize) then
eQP(:) = eQPlin(:)
else
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Newton-Raphson for Dyson equation not implemented '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
stop
end if
deallocate(eQPlin,Z,SigC)
! Check one-body converge
err_1b = maxval(abs(eOld - eQP))
eOld(:) = eQP(:)
call wall_time(end_1b)
t_1b = end_1b - start_1b
write(*,'(A50,1X,F9.3,A8)') 'Wall time for one-body iteration =',t_1b,' seconds'
end do
!---------------------------------------------!
! End main loop for one-body self-consistency !
!---------------------------------------------!
! Did it actually converge?
if(n_it_1b == max_it_1b) 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 GParquet

13
src/Parquet/RParquet.f90
View File

@ -37,16 +37,13 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
integer :: nVVs,nVVt
double precision :: EcRPA
double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:)
double precision,allocatable :: Aph(:,:), 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 :: Bpp(:,:), Cpp(:,:), Dpp(:,:)
double precision,allocatable :: X1s(:,:),X1t(:,:)
double precision,allocatable :: Y1s(:,:),Y1t(:,:)
double precision,allocatable :: ee_sing_Om(:), old_ee_sing_Om(:)
@ -73,7 +70,8 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
double precision :: EcGM
! Output variables
! None
nOOs = nO*(nO + 1)/2
nVVs = nV*(nV + 1)/2
nOOt = nO*(nO - 1)/2
@ -166,7 +164,6 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
n_it_2b = n_it_2b + 1
call wall_time(start_2b)
!TODO add some timers everywhere
write(*,*)' ***********************************'
write(*,'(1X,A30,1X,I4)') 'Two-body iteration #',n_it_2b
write(*,*)' ***********************************'
@ -659,4 +656,4 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
end if
end subroutine
end subroutine RParquet

1
src/Parquet/R_screened_integrals.f90
View File

@ -31,6 +31,7 @@ subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Gam,XpY,r
rho(p,q,ia) = rho(p,q,ia) &
+ (2d0*ERI(p,j,q,b) - ERI(p,j,b,q))*XpY(ia,jb) &
+ 1d0*eh_sing_Gam(p,j,q,b)*XpY(ia,jb)
end do
end do
end do

12
src/QuAcK/GQuAcK.f90
View File

@ -7,7 +7,8 @@ subroutine GQuAcK(working_dir,dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dop
TDA,maxSCF_GF,max_diis_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)
implicit none
include 'parameters.h'
@ -74,6 +75,9 @@ subroutine GQuAcK(working_dir,dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dop
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,doRPA,doGF,doGW,doGT
@ -341,9 +345,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 GParquet(max_it_macro,conv_one_body,max_it_micro,conv_two_body, &
nBas,nC,nO,nV,nR,nS, &
eHF,ERI_MO)
call wall_time(end_Parquet)
t_Parquet = end_Parquet - start_Parquet

6
src/QuAcK/QuAcK.f90
View File

@ -289,14 +289,14 @@ program QuAcK
if(doGQuAcK) &
call GQuAcK(working_dir,doGtest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2,doG0T0pp,doevGTpp,doqsGTpp,doParquet, &
doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2,doG0T0pp,doevGTpp,doqsGTpp,doParquet, &
nNuc,nBas,sum(nC),sum(nO),sum(nV),sum(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,TDA,maxSCF_GF,max_diis_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)
!--------------------------!
! Bogoliubov QuAcK branch !

5
src/QuAcK/RQuAcK.f90
View File

@ -112,6 +112,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,
double precision,allocatable :: ERI_MO(:,:,:,:)
integer :: ixyz
integer :: nS
double precision,allocatable :: eGW(:)
write(*,*)
write(*,*) '******************************'
@ -130,6 +131,8 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,
allocate(dipole_int_MO(nOrb,nOrb,ncart))
allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb))
allocate(eGW(nOrb))
allocate(ERI_AO(nBas,nBas,nBas,nBas))
call wall_time(start_int)
call read_2e_integrals(working_dir,nBas,ERI_AO)
@ -337,7 +340,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,
call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF_GW,thresh_GW,max_diis_GW, &
doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T, &
V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF,eGW)
call wall_time(end_GW)
t_GW = end_GW - start_GW