PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-05-06 07:05:33 +02:00
Code Releases Activity

fixed complex qsGF

Browse Source
This commit is contained in:
Loris Burth 2025-05-05 13:43:20 +02:00
parent b7489c7bd4
commit 0808dc4be8
5 changed files with 104 additions and 105 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
src/GF/complex_RGF.f90
View File

@ -1,6 +1,6 @@
subroutine complex_RGF(dotest,docG0F2,doevGF2,doqsGF2,maxSCF, &
thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize, &
eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF, &
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, &
CAP_AO,CAP_MO)
@ -27,7 +27,7 @@ subroutine complex_RGF(dotest,docG0F2,doevGF2,doqsGF2,maxSCF,
logical,intent(in) :: triplet
logical,intent(in) :: linearize
double precision,intent(in) :: eta
logical,intent(in) :: regularize
logical,intent(in) :: doSRG
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
@ -71,7 +71,7 @@ subroutine complex_RGF(dotest,docG0F2,doevGF2,doqsGF2,maxSCF,
call wall_time(start_GF)
call complex_cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS, &
linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS, &
ENuc,ERHF,ERI_MO,CAP_MO,dipole_int_MO,eHF)
call wall_time(end_GF)
@ -84,7 +84,7 @@ subroutine complex_RGF(dotest,docG0F2,doevGF2,doqsGF2,maxSCF,
call wall_time(start_GF)
call complex_evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,triplet, &
linearize,eta,regularize,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)
call wall_time(end_GF)
t_GF = end_GF - start_GF
@ -97,7 +97,7 @@ subroutine complex_RGF(dotest,docG0F2,doevGF2,doqsGF2,maxSCF,
call wall_time(start_GF)
call complex_qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA, &
dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc, &
dBSE,dTDA,singlet,triplet,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, &
CAP_AO,CAP_MO)

10
src/GF/complex_cRG0F2.f90
View File

@ -1,4 +1,4 @@
subroutine complex_cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, &
subroutine complex_cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,doSRG, &
nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,CAP,dipole_int,eHF)
! Perform a one-shot second-order Green function calculation
@ -19,7 +19,7 @@ subroutine complex_cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,l
logical,intent(in) :: triplet
logical,intent(in) :: linearize
double precision,intent(in) :: eta
logical,intent(in) :: regularize
logical,intent(in) :: doSRG
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
@ -66,10 +66,10 @@ subroutine complex_cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,l
Re_eGFlin(nOrb),Im_eGFlin(nOrb), Re_eGF(nOrb),Im_eGF(nOrb),Re_eHF(nOrb),Im_eHF(nOrb))
Re_eHF(:) = real(eHF(:))
Im_eHF(:) = aimag(eHF(:))
flow = 100d0
flow = 500d0
! Frequency-dependent second-order contribution
if(regularize) then
if(doSRG) then
call complex_cRGF2_SRG_self_energy_diag(flow,eta,nOrb,nC,nO,nV,nR,Re_eHF,Im_eHF,ERI,Re_SigC,Im_SigC,Re_Z,Im_Z)
else
call complex_cRGF2_self_energy_diag(eta,nOrb,nC,nO,nV,nR,Re_eHF,Im_eHF,ERI,Re_SigC,Im_SigC,Re_Z,Im_Z)
@ -88,7 +88,7 @@ subroutine complex_cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,l
write(*,*) ' *** Quasiparticle energies obtained by root search *** '
write(*,*)
call complex_cRGF2_QP_graph(flow,regularize,eta,nOrb,nC,nO,nV,nR,Re_eHF,Im_eHF,ERI,Re_eGFlin,Im_eGFlin,&
call complex_cRGF2_QP_graph(doSRG,eta,flow,nOrb,nC,nO,nV,nR,Re_eHF,Im_eHF,ERI,Re_eGFlin,Im_eGFlin,&
Re_eHF,Im_eHF,Re_eGF,Im_eGF,Re_Z,Im_Z)
end if

6
src/GF/complex_cRGF2_QP_graph.f90
View File

@ -1,4 +1,4 @@
subroutine complex_cRGF2_QP_graph(flow,reg,eta,nBas,nC,nO,nV,nR,Re_eHF,Im_eHF,&
subroutine complex_cRGF2_QP_graph(doSRG,eta,flow,nBas,nC,nO,nV,nR,Re_eHF,Im_eHF,&
ERI,Re_eGFlin,Im_eGFlin,Re_eOld,Im_eold,Re_eGF,Im_eGF,Re_Z,Im_Z)
! Compute the graphical solution of the complex GF2 QP equation
@ -10,7 +10,7 @@ subroutine complex_cRGF2_QP_graph(flow,reg,eta,nBas,nC,nO,nV,nR,Re_eHF,Im_eHF,&
double precision,intent(in) :: eta
double precision,intent(in) :: flow
logical, intent(in) :: reg
logical, intent(in) :: doSRG
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
@ -62,7 +62,7 @@ subroutine complex_cRGF2_QP_graph(flow,reg,eta,nBas,nC,nO,nV,nR,Re_eHF,Im_eHF,&
do while (abs(cmplx(Re_f,Im_f,kind=8)) > thresh .and. nIt < maxIt)
nIt = nIt + 1
if(reg) then
if(doSRG) then
call complex_cRGF2_SRG_SigC_dSigC(flow,p,eta,nBas,nC,nO,nV,nR,Re_w,Im_w,Re_eOld,Im_eOld,ERI,&
Re_SigC,Im_SigC,Re_dSigC,Im_dSigC)
else

10
src/GF/complex_evRGF2.f90
View File

@ -1,5 +1,5 @@
subroutine complex_evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,triplet, &
linearize,eta,regularize,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)
! Perform eigenvalue self-consistent second-order Green function calculation
@ -22,7 +22,7 @@ subroutine complex_evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max
logical,intent(in) :: triplet
logical,intent(in) :: linearize
double precision,intent(in) :: eta
logical,intent(in) :: regularize
logical,intent(in) :: doSRG
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
@ -104,7 +104,7 @@ subroutine complex_evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max
! Frequency-dependent second-order contribution
if(regularize) then
if(doSRG) then
call complex_cRGF2_SRG_self_energy_diag(flow,eta,nOrb,nC,nO,nV,nR,Re_eGF,Im_eGF,ERI,Re_SigC,Im_SigC,Re_Z,Im_Z)
@ -125,7 +125,7 @@ subroutine complex_evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max
write(*,*) ' *** Quasiparticle energies obtained by root search *** '
write(*,*)
call complex_cRGF2_QP_graph(flow,regularize,eta,nOrb,nC,nO,nV,nR,Re_eHF,Im_eHF,&
call complex_cRGF2_QP_graph(doSRG,eta,flow,nOrb,nC,nO,nV,nR,Re_eHF,Im_eHF,&
ERI,Re_eOld,Im_eOld,Re_eOld,Im_eOld,Re_eGF,Im_eGF,Re_Z,Im_Z)
eGF = cmplx(Re_eGF,Im_eGF,kind=8)
end if
@ -134,7 +134,7 @@ subroutine complex_evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max
! Print results
!call RMP2(.false.,regularize,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
!call RMP2(.false.,doSRG,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
call print_complex_evRGF2(nOrb,nO,nSCF,Conv,Re_eHF,Im_eHF,ENuc,ERHF,Re_SigC,Im_SigC,Re_Z,Im_Z,Re_eGF,Im_eGF)
! DIIS extrapolation

171
src/GF/complex_qsRGF2.f90
View File

@ -1,5 +1,5 @@
subroutine complex_qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA, &
dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc, &
dBSE,dTDA,singlet,triplet,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, &
CAP_AO,CAP_MO)
@ -24,7 +24,7 @@ subroutine complex_qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA, &
logical,intent(in) :: singlet
logical,intent(in) :: triplet
double precision,intent(in) :: eta
logical,intent(in) :: regularize
logical,intent(in) :: doSRG
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
@ -68,12 +68,11 @@ subroutine complex_qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA, &
complex*16 :: Ec
complex*16 :: EcBSE(nspin)
complex*16,allocatable :: error_diis(:,:)
complex*16,allocatable :: err_diis(:,:)
complex*16,allocatable :: F_diis(:,:)
complex*16,allocatable :: c(:,:)
complex*16,allocatable :: cp(:,:)
complex*16,allocatable :: eGF(:)
complex*16,allocatable :: eOld(:)
complex*16,allocatable :: P(:,:)
complex*16,allocatable :: F(:,:)
complex*16,allocatable :: Fp(:,:)
@ -82,7 +81,7 @@ subroutine complex_qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA, &
complex*16,allocatable :: SigC(:,:)
complex*16,allocatable :: SigCp(:,:)
complex*16,allocatable :: Z(:)
complex*16,allocatable :: error(:,:)
complex*16,allocatable :: err(:,:)
! Hello world
@ -93,6 +92,17 @@ subroutine complex_qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA, &
write(*,*)'********************************'
write(*,*)
! SRG regularization
flow = 500d0
if(doSRG) then
write(*,*) '*** SRG regularized qsGF2 scheme ***'
write(*,*)
end if
! Warning
write(*,*) '!! ERIs in MO basis will be overwritten in qsGF2 !!'
@ -101,6 +111,7 @@ subroutine complex_qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA, &
! Stuff
nBas_Sq = nBas*nBas
! TDA
if(TDA) then
@ -111,8 +122,6 @@ subroutine complex_qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA, &
! Memory allocation
allocate(eGF(nOrb))
allocate(eOld(nOrb))
allocate(c(nBas,nOrb))
allocate(cp(nOrb,nOrb))
@ -122,14 +131,14 @@ subroutine complex_qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA, &
allocate(F(nBas,nBas))
allocate(J(nBas,nBas))
allocate(K(nBas,nBas))
allocate(error(nBas,nBas))
allocate(err(nBas,nBas))
allocate(Z(nOrb))
allocate(SigC(nOrb,nOrb))
allocate(SigCp(nBas,nBas))
allocate(error_diis(nBas_Sq,max_diis))
allocate(err_diis(nBas_Sq,max_diis))
allocate(F_diis(nBas_Sq,max_diis))
! Initialization
@ -139,13 +148,13 @@ subroutine complex_qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA, &
ispin = 1
Conv = 1d0
P(:,:) = PHF(:,:)
eOld(:) = eHF(:)
eGF(:) = eHF(:)
c(:,:) = cHF(:,:)
F_diis(:,:) = 0d0
error_diis(:,:) = 0d0
err_diis(:,:) = 0d0
rcond = 0d0
flow = 500d0
!------------------------------------------------------------------------
! Main loop
@ -171,7 +180,7 @@ subroutine complex_qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA, &
! Compute self-energy and renormalization factor
if(regularize) then
if(doSRG) then
call complex_cRGF2_SRG_self_energy(flow,eta, nOrb, nC, nO, nV, nR, eGF, ERI_MO, SigC, Z)
@ -191,52 +200,15 @@ subroutine complex_qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA, &
F(:,:) = cmplx(Hc(:,:),CAP_AO(:,:),kind=8) + J(:,:) + 0.5d0*K(:,:) + SigCp(:,:)
if(nBas .ne. nOrb) then
call complex_complex_AOtoMO(nBas, nOrb, c(1,1), F(1,1), Fp(1,1))
call complex_MOtoAO(nBas, nOrb, S(1,1), c(1,1), Fp(1,1), F(1,1))
call complex_complex_AOtoMO(nBas, nOrb, c, F, Fp)
call complex_MOtoAO(nBas, nOrb, S, c, Fp, F)
endif
! Compute commutator and convergence criteria
error = matmul(F, matmul(P, S)) - matmul(matmul(S, P), F)
err = matmul(F, matmul(P, S)) - matmul(matmul(S, P), F)
! DIIS extrapolation
n_diis = min(n_diis+1, max_diis)
if(abs(rcond) > 1d-7) then
call complex_DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,error_diis,F_diis,error,F)
else
n_diis = 0
end if
! Diagonalize Hamiltonian in AO basis
if(nBas .eq. nOrb) then
Fp = matmul(transpose(X), matmul(F, X))
cp(:,:) = Fp(:,:)
call complex_diagonalize_matrix(nOrb, cp, eGF)
call complex_orthogonalize_matrix(nBas,cp)
c = matmul(X, cp)
else
Fp = matmul(transpose(c), matmul(F, c))
cp(:,:) = Fp(:,:)
call complex_diagonalize_matrix(nOrb, cp, eGF)
call complex_orthogonalize_matrix(nBas,cp)
c = matmul(c, cp)
endif
! Compute new density matrix in the AO basis
P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO)))
! Save quasiparticles energy for next cycle
Conv = maxval(abs(eGF - eOld))
eOld(:) = eGF(:)
!------------------------------------------------------------------------
! Compute total energy
!------------------------------------------------------------------------
Conv = maxval(abs(err))
! Kinetic energy
@ -258,14 +230,41 @@ subroutine complex_qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA, &
Ex = 0.25d0*complex_trace_matrix(nBas, matmul(P, K))
! Correlation energy
!call RMP2(.false., regularize, nOrb, nC, nO, nV, nR, ERI_MO, ENuc, EqsGF2, eGF, Ec)
! Total energy
EqsGF2 = ET + EV + EJ + Ex + Ec
EqsGF2 = ET + EV + EJ + Ex + Ec + EW
! DIIS extrapolation
if(max_diis>1) then
n_diis = min(n_diis+1,max_diis)
call complex_DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,F)
end if
! Diagonalize Hamiltonian in AO basis
if(nBas == nOrb) then
Fp = matmul(transpose(X), matmul(F, X))
cp(:,:) = Fp(:,:)
call complex_diagonalize_matrix(nOrb, cp, eGF)
call complex_orthogonalize_matrix(nBas,cp)
c = matmul(X, cp)
else
Fp = matmul(transpose(c), matmul(F, c))
cp(:,:) = Fp(:,:)
call complex_diagonalize_matrix(nOrb, cp, eGF)
call complex_orthogonalize_matrix(nBas,cp)
c = matmul(c, cp)
endif
call complex_complex_AOtoMO(nBas,nOrb,c,SigCp,SigC)
! Compute new density matrix in the AO basis
P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO)))
!------------------------------------------------------------------------
! Print results
@ -289,14 +288,14 @@ subroutine complex_qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA, &
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, error, error_diis, F_diis)
deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, err, err_diis, F_diis)
stop
end if
! Deallocate memory
deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, error, error_diis, F_diis)
deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, err, err_diis, F_diis)
!! Perform phBSE@GF2 calculation
!
@ -318,31 +317,31 @@ subroutine complex_qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA, &
! Perform ppBSE@GF2 calculation
if(doppBSE) then
call RGF2_ppBSE(TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, &
nC, nO, nV, nR, ERI_MO, dipole_int_MO, eGF, EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGF2 correlation energy (singlet) =',EcBSE(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGF2 correlation energy (triplet) =',3d0*EcBSE(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGF2 correlation energy =',EcBSE(1) + 3d0*EcBSE(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGF2 total energy =',ENuc + EqsGF2 + EcBSE(1) + 3d0*EcBSE(2),' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end if
! Testing zone
if(dotest) then
call dump_test_value('R','qsGF2 correlation energy',Ec)
call dump_test_value('R','qsGF2 HOMO energy',eGF(nO))
call dump_test_value('R','qsGF2 LUMO energy',eGF(nO+1))
end if
!
! if(doppBSE) then
!
! call RGF2_ppBSE(TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, &
! nC, nO, nV, nR, ERI_MO, dipole_int_MO, eGF, EcBSE)
!
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGF2 correlation energy (singlet) =',EcBSE(1),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGF2 correlation energy (triplet) =',3d0*EcBSE(2),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGF2 correlation energy =',EcBSE(1) + 3d0*EcBSE(2),' au'
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGF2 total energy =',ENuc + EqsGF2 + EcBSE(1) + 3d0*EcBSE(2),' au'
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
!
! end if
!
!! Testing zone
!
! if(dotest) then
!
! call dump_test_value('R','qsGF2 correlation energy',Ec)
! call dump_test_value('R','qsGF2 HOMO energy',eGF(nO))
! call dump_test_value('R','qsGF2 LUMO energy',eGF(nO+1))
!
! end if
end subroutine