PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:34:46 +01:00
Code Releases Activity

new routine to fix GTpp

Browse Source
This commit is contained in:
Pierre-Francois Loos 2024-09-12 20:10:02 +02:00
parent 6290634d87
commit 8ec7d3d170
8 changed files with 218 additions and 196 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

105
src/GT/RG0T0pp.f90
View File

@ -41,7 +41,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
! Local variables ! Local variables
logical :: print_T = .true. logical :: print_T = .true.
double precision :: lambda logical :: plot_self = .false.
integer :: isp_T integer :: isp_T
integer :: iblock integer :: iblock
@ -81,10 +81,6 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
write(*,*)'*********************************' write(*,*)'*********************************'
write(*,*) write(*,*)
! Initialization
lambda = 1d0
! TDA for T ! TDA for T
if(TDA_T) then if(TDA_T) then
@ -125,33 +121,13 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,lambda,eHF,ERI,Cpp) if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,lambda,eHF,ERI,Dpp) call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp)
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,lambda,ERI,Bpp) call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp)
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T)) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T))
deallocate(Bpp,Cpp,Dpp)
!print*, 'LAPACK:'
!print*, Om2s
!print*, Om1s
!n_states = nOOs + 5 deallocate(Bpp,Cpp,Dpp)
!n_states_diag = n_states + 4
!allocate(Om(nOOs+nVVs), R(nOOs+nVVs,n_states_diag))
!allocate(supp_data_dbl(1), supp_data_int(1))
!supp_data_int(1) = 0
!supp_data_dbl(1) = 0.d0
!call ppLR_davidson(iblock, TDA_T, nC, nO, nR, nOrb, nOOs, nVVs, &
! 1.d0, & ! lambda
! eHF(1), &
! 0.d0, & ! eF
! ERI(1,1,1,1), &
! supp_data_int(1), 1, &
! supp_data_dbl(1), 1, &
! Om(1), R(1,1), n_states, n_states_diag, "RPA")
!deallocate(supp_data_dbl, supp_data_int)
!deallocate(Om, R)
!stop
if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-beta)',nVVs,Om1s) if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-beta)',nVVs,Om1s)
if(print_T) call print_excitation_energies('ppRPA@RHF','2h (alpha-beta)',nOOs,Om2s) if(print_T) call print_excitation_energies('ppRPA@RHF','2h (alpha-beta)',nOOs,Om2s)
@ -168,34 +144,13 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,lambda,eHF,ERI,Cpp) if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,lambda,eHF,ERI,Dpp) call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp)
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,lambda,ERI,Bpp) call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp)
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T)) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
deallocate(Bpp,Cpp,Dpp) deallocate(Bpp,Cpp,Dpp)
!print*, 'LAPACK:'
!print*, Om2t
!print*, Om1t
!n_states = nOOt + 5
!n_states_diag = n_states + 4
!allocate(Om(nOOt+nVVt), R(nOOt+nVVt,n_states_diag))
!allocate(supp_data_dbl(1), supp_data_int(1))
!supp_data_int(1) = 0
!supp_data_dbl(1) = 0.d0
!call ppLR_davidson(iblock, TDA_T, nC, nO, nR, nOrb, nOOt, nVVt, &
! 1.d0, & ! lambda
! eHF(1), &
! 0.d0, & ! eF
! ERI(1,1,1,1), &
! supp_data_int(1), 1, &
! supp_data_dbl(1), 1, &
! Om(1), R(1,1), n_states, n_states_diag, "RPA")
!deallocate(Om, R)
!deallocate(supp_data_dbl)
!stop
if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-alpha)',nVVt,Om1t) if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-alpha)',nVVt,Om1t)
if(print_T) call print_excitation_energies('ppRPA@RHF','2h (alpha-alpha)',nOOt,Om2t) if(print_T) call print_excitation_energies('ppRPA@RHF','2h (alpha-alpha)',nOOt,Om2t)
@ -249,8 +204,8 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
end if end if
! call RGTpp_plot_self_energy(nOrb,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,eGT,Om1s,rho1s,Om2s,rho2s, & if(plot_self) call RGTpp_plot_self_energy(nOrb,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,eGT,Om1s,rho1s,Om2s,rho2s, &
! Om1t,rho1t,Om2t,rho2t) Om1t,rho1t,Om2t,rho2t)
!---------------------------------------------- !----------------------------------------------
! Dump results ! Dump results
@ -264,9 +219,9 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,lambda,eGT,ERI,Cpp) if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,lambda,eGT,ERI,Dpp) call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,lambda,ERI,Bpp) call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T)) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T))
@ -278,9 +233,9 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,lambda,eGT,ERI,Cpp) if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,lambda,eGT,ERI,Dpp) call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp)
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,lambda,ERI,Bpp) call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp)
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T)) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
@ -295,17 +250,10 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
if(dophBSE) then if(dophBSE) then
call RGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & call RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, & Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, &
ERI,dipole_int,eHF,eGT,EcBSE) ERI,dipole_int,eHF,eGT,EcBSE)
if(exchange_kernel) then
EcBSE(1) = 0.5d0*EcBSE(1)
EcBSE(2) = 1.5d0*EcBSE(1)
end if
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp correlation energy (singlet) =',EcBSE(1),' au' write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp correlation energy (singlet) =',EcBSE(1),' au'
@ -319,29 +267,10 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
if(doACFDT) then if(doACFDT) then
write(*,*) '--------------------------------------------------------'
write(*,*) 'Adiabatic connection version of phBSE correlation energy'
write(*,*) '--------------------------------------------------------'
write(*,*)
if(doXBS) then
write(*,*) '*** scaled screening version (XBS) ***'
write(*,*)
end if
call RGTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS, & call RGTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS, &
nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, &
Om2t,X2t,Y2t,rho1t,rho2t,ERI,eHF,eGT,EcBSE) Om2t,X2t,Y2t,rho1t,rho2t,ERI,eHF,eGT,EcBSE)
if(exchange_kernel) then
EcBSE(1) = 0.5d0*EcBSE(1)
EcBSE(2) = 1.5d0*EcBSE(2)
end if
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp@RHF correlation energy (singlet) = ',EcBSE(1),' au' write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp@RHF correlation energy (singlet) = ',EcBSE(1),' au'

16
src/GT/RGTpp_phACFTD.f90
View File

@ -85,6 +85,22 @@ subroutine RGTpp_phACFDT(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
allocate(Aph(nS,nS),Bph(nS,nS),TAs(nS,nS),TBs(nS,nS),TAt(nS,nS),TBt(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS)) allocate(Aph(nS,nS),Bph(nS,nS),TAs(nS,nS),TBs(nS,nS),TAt(nS,nS),TBt(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS))
allocate(Ec(nAC,nspin)) allocate(Ec(nAC,nspin))
! Hello World
write(*,*) '-------------------------------------------------------------'
write(*,*) ' Adiabatic connection version of BSE@GTpp correlation energy '
write(*,*) '-------------------------------------------------------------'
write(*,*)
! eXtended BSE
if(doXBS) then
write(*,*) '*** scaled screening version (XBS) ***'
write(*,*)
end if
! Antisymmetrized kernel version ! Antisymmetrized kernel version
if(exchange_kernel) then if(exchange_kernel) then

10
src/GT/RGTpp_phBSE.f90
View File

@ -1,4 +1,4 @@
subroutine RGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, &
Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,rho1ab,rho2ab,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, & Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,rho1ab,rho2ab,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, &
ERI,dipole_int,eT,eGT,EcBSE) ERI,dipole_int,eT,eGT,EcBSE)
@ -9,6 +9,7 @@ subroutine RGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
! Input variables ! Input variables
logical,intent(in) :: exchange_kernel
logical,intent(in) :: TDA_T logical,intent(in) :: TDA_T
logical,intent(in) :: TDA logical,intent(in) :: TDA
logical,intent(in) :: dBSE logical,intent(in) :: dBSE
@ -189,4 +190,11 @@ subroutine RGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
dipole_int,OmBSE,XpY_BSE,XmY_BSE,TAab,TAaa) dipole_int,OmBSE,XpY_BSE,XmY_BSE,TAab,TAaa)
end if end if
if(exchange_kernel) then
EcBSE(1) = 0.5d0*EcBSE(1)
EcBSE(2) = 1.5d0*EcBSE(1)
end if
end subroutine end subroutine

35
src/GT/evRGTpp.f90
View File

@ -222,12 +222,12 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
! DIIS extrapolation ! DIIS extrapolation
if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis) n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nOrb,nOrb,n_diis,error_diis,e_diis,eGT(:)-eOld(:),eGT(:)) call DIIS_extrapolation(rcond,nOrb,nOrb,n_diis,error_diis,e_diis,eGT(:)-eOld(:),eGT(:))
! Reset DIIS if required end if
if(abs(rcond) < 1d-15) n_diis = 0
! Save quasiparticles energy for next cycle ! Save quasiparticles energy for next cycle
@ -246,17 +246,10 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
if(BSE) then if(BSE) then
call RGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & call RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, & Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, &
ERI,dipole_int,eGT,eGT,EcBSE) ERI,dipole_int,eGT,eGT,EcBSE)
if(exchange_kernel) then
EcRPA(1) = 0.5d0*EcRPA(1)
EcRPA(2) = 1.5d0*EcRPA(1)
end if
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGTpp correlation energy (singlet) =',EcBSE(1) write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGTpp correlation energy (singlet) =',EcBSE(1)
@ -270,29 +263,10 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
if(doACFDT) then if(doACFDT) then
write(*,*) '------------------------------------------------------'
write(*,*) 'Adiabatic connection version of BSE correlation energy'
write(*,*) '------------------------------------------------------'
write(*,*)
if(doXBS) then
write(*,*) '*** scaled screening version (XBS) ***'
write(*,*)
end if
call RGTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS, & call RGTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS, &
nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, &
Om2t,X2t,Y2t,rho1t,rho2t,ERI,eGT,eGT,EcBSE) Om2t,X2t,Y2t,rho1t,rho2t,ERI,eGT,eGT,EcBSE)
if(exchange_kernel) then
EcBSE(1) = 0.5d0*EcBSE(1)
EcBSE(2) = 1.5d0*EcBSE(2)
end if
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'AC@phBSE@evGTpp correlation energy (singlet) =',EcBSE(1) write(*,'(2X,A50,F20.10)') 'AC@phBSE@evGTpp correlation energy (singlet) =',EcBSE(1)
@ -306,7 +280,6 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
end if end if
! Testing zone ! Testing zone
if(dotest) then if(dotest) then

93
src/GT/qsRGTpp.f90
View File

@ -1,9 +1,6 @@
subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA, &
! --- dBSE,dTDA,singlet,triplet,eta,regularize,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)
subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, TDA_T, TDA, &
dBSE, dTDA, singlet, triplet, eta, regularize, 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)
! Perform a quasiparticle self-consistent GT calculation ! Perform a quasiparticle self-consistent GT calculation
@ -34,7 +31,8 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d
double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc double precision,intent(in) :: ENuc
integer,intent(in) :: nBas, nOrb integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nC,nO,nV,nR,nS integer,intent(in) :: nC,nO,nV,nR,nS
double precision,intent(in) :: ERHF double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nOrb) double precision,intent(in) :: eHF(nOrb)
@ -154,10 +152,10 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d
allocate(error_diis(nBas_Sq,max_diis)) allocate(error_diis(nBas_Sq,max_diis))
allocate(F_diis(nBas_Sq,max_diis)) allocate(F_diis(nBas_Sq,max_diis))
allocate(Om1s(nVVs), X1s(nVVs,nVVs), Y1s(nOOs,nVVs), rho1s(nOrb,nOrb,nVVs)) allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs),rho1s(nOrb,nOrb,nVVs))
allocate(Om2s(nOOs), X2s(nVVs,nOOs), Y2s(nOOs,nOOs), rho2s(nOrb,nOrb,nOOs)) allocate(Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs),rho2s(nOrb,nOrb,nOOs))
allocate(Om1t(nVVt), X1t(nVVt,nVVt), Y1t(nOOt,nVVt), rho1t(nOrb,nOrb,nVVt)) allocate(Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt),rho1t(nOrb,nOrb,nVVt))
allocate(Om2t(nOOt), X2t(nVVt,nOOt), Y2t(nOOt,nOOt), rho2t(nOrb,nOrb,nOOt)) allocate(Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt),rho2t(nOrb,nOrb,nOOt))
! Initialization ! Initialization
@ -193,7 +191,7 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d
! AO to MO transformation of two-electron integrals ! AO to MO transformation of two-electron integrals
call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO) call AOtoMO_ERI_RHF(nBas,nOrb,c,ERI_AO,ERI_MO)
! Compute linear response ! Compute linear response
@ -202,9 +200,9 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp)
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp) call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp) call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp)
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp)
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
@ -215,9 +213,9 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp)
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp) call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp) call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp)
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp)
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
@ -246,14 +244,14 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d
Sig = 0.5d0*(Sig + transpose(Sig)) Sig = 0.5d0*(Sig + transpose(Sig))
call MOtoAO(nBas, nOrb, S, c, Sig, Sigp) call MOtoAO(nBas,nOrb,S,c,Sig,Sigp)
! Solve the quasi-particle equation ! Solve the quasi-particle equation
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + Sigp(:,:) F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + Sigp(:,:)
if(nBas .ne. nOrb) then if(nBas .ne. nOrb) then
call AOtoMO(nBas, nOrb, c(1,1), F(1,1), Fp(1,1)) call AOtoMO(nBas,nOrb,c(1,1),F(1,1),Fp(1,1))
call MOtoAO(nBas, nOrb, S(1,1), c(1,1), Fp(1,1), F(1,1)) call MOtoAO(nBas,nOrb,S(1,1),c(1,1),Fp(1,1),F(1,1))
endif endif
! Compute commutator and convergence criteria ! Compute commutator and convergence criteria
@ -272,20 +270,20 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d
! Diagonalize Hamiltonian in AO basis ! Diagonalize Hamiltonian in AO basis
if(nBas .eq. nOrb) then if(nBas .eq. nOrb) then
Fp = matmul(transpose(X), matmul(F, X)) Fp = matmul(transpose(X),matmul(F,X))
cp(:,:) = Fp(:,:) cp(:,:) = Fp(:,:)
call diagonalize_matrix(nOrb, cp, eGT) call diagonalize_matrix(nOrb,cp,eGT)
c = matmul(X, cp) c = matmul(X,cp)
else else
Fp = matmul(transpose(c), matmul(F, c)) Fp = matmul(transpose(c),matmul(F,c))
cp(:,:) = Fp(:,:) cp(:,:) = Fp(:,:)
call diagonalize_matrix(nOrb, cp, eGT) call diagonalize_matrix(nOrb,cp,eGT)
c = matmul(c, cp) c = matmul(c,cp)
endif endif
! Compute new density matrix in the AO basis ! Compute new density matrix in the AO basis
P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO))) P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
! Save quasiparticles energy for next cycle ! Save quasiparticles energy for next cycle
@ -319,9 +317,8 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d
! Print results ! Print results
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
call print_qsRGTpp(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, & call print_qsRGTpp(nBas,nOrb,nO,nSCF,Conv,thresh,eHF,eGT,c,Sig,Z, &
eGT, c, Sig, Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, & ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole)
EqsGT, dipole)
end do end do
!------------------------------------------------------------------------ !------------------------------------------------------------------------
@ -338,34 +335,27 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*) write(*,*)
deallocate(c, cp, P, F, Fp, J, K, Sig, Sigp, Z, error, error_diis, F_diis) deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Z,error,error_diis,F_diis)
deallocate(Om1s, X1s, Y1s, rho1s) deallocate(Om1s,X1s,Y1s,rho1s)
deallocate(Om2s, X2s, Y2s, rho2s) deallocate(Om2s,X2s,Y2s,rho2s)
deallocate(Om1t, X1t, Y1t, rho1t) deallocate(Om1t,X1t,Y1t,rho1t)
deallocate(Om2t, X2t, Y2t, rho2t) deallocate(Om2t,X2t,Y2t,rho2t)
stop stop
end if end if
! Deallocate memory ! Deallocate memory
deallocate(c, cp, P, F, Fp, J, K, Sig, Sigp, Z, error, error_diis, F_diis) deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Z,error,error_diis,F_diis)
! Perform BSE calculation ! Perform BSE calculation
if(dophBSE) then if(dophBSE) then
call RGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & call RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, & Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, &
ERI_MO,dipole_int_MO,eGT,eGT,EcBSE) ERI_MO,dipole_int_MO,eGT,eGT,EcBSE)
if(exchange_kernel) then
EcBSE(1) = 0.5d0*EcBSE(1)
EcBSE(2) = 1.5d0*EcBSE(2)
end if
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGTpp correlation energy (singlet) =',EcBSE(1) write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGTpp correlation energy (singlet) =',EcBSE(1)
@ -379,18 +369,6 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d
if(doACFDT) then if(doACFDT) then
write(*,*) '------------------------------------------------------'
write(*,*) 'Adiabatic connection version of BSE correlation energy'
write(*,*) '------------------------------------------------------'
write(*,*)
if(doXBS) then
write(*,*) '*** scaled screening version (XBS) ***'
write(*,*)
end if
call RGTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS, & call RGTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS, &
nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, &
Om2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcBSE) Om2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcBSE)
@ -408,7 +386,6 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d
end if end if
! Testing zone ! Testing zone
if(dotest) then if(dotest) then
@ -419,9 +396,9 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d
end if end if
deallocate(Om1s, X1s, Y1s, rho1s) deallocate(Om1s,X1s,Y1s,rho1s)
deallocate(Om2s, X2s, Y2s, rho2s) deallocate(Om2s,X2s,Y2s,rho2s)
deallocate(Om1t, X1t, Y1t, rho1t) deallocate(Om1t,X1t,Y1t,rho1t)
deallocate(Om2t, X2t, Y2t, rho2t) deallocate(Om2t,X2t,Y2t,rho2t)
end subroutine end subroutine

10
src/LR/ppLR.f90
View File

@ -43,7 +43,6 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
double precision, external :: trace_matrix double precision, external :: trace_matrix
N = nOO + nVV N = nOO + nVV
allocate(M(N,N), Z(N,N), Om(N)) allocate(M(N,N), Z(N,N), Om(N))
@ -70,10 +69,13 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
! if((nOO .eq. 0) .or. (nVV .eq. 0)) then ! if((nOO .eq. 0) .or. (nVV .eq. 0)) then
! Diagonalize the p-p matrix ! Diagonalize the pp matrix
if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV, M, Om, Z)
if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Om,Z)
! Split the various quantities in p-p and h-h parts ! Split the various quantities in p-p and h-h parts
call sort_ppRPA(nOO, nVV, Om, Z, Om1, X1, Y1, Om2, X2, Y2)
call sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
! else ! else

5
src/RPA/sort_ppRPA.f90
View File

@ -17,7 +17,6 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
integer :: pq,ab,ij integer :: pq,ab,ij
! integer :: deg1,ab_start,ab_end ! integer :: deg1,ab_start,ab_end
! integer :: deg2,ij_start,ij_end ! integer :: deg2,ij_start,ij_end
integer :: OO,VV
double precision,allocatable :: M(:,:) double precision,allocatable :: M(:,:)
double precision,allocatable :: Z1(:,:) double precision,allocatable :: Z1(:,:)
double precision,allocatable :: Z2(:,:) double precision,allocatable :: Z2(:,:)
@ -211,8 +210,8 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
! S1 = + matmul(transpose(Z1),matmul(M,Z1)) ! S1 = + matmul(transpose(Z1),matmul(M,Z1))
! S2 = - matmul(transpose(Z2),matmul(M,Z2)) ! S2 = - matmul(transpose(Z2),matmul(M,Z2))
if(nVV > 0) call orthogonalization_matrix(nVV,VV,S1,O1) if(nVV > 0) call orthogonalize_matrix(1,nVV,S1,O1)
if(nOO > 0) call orthogonalization_matrix(nOO,OO,S2,O2) if(nOO > 0) call orthogonalize_matrix(1,nOO,S2,O2)
if(nVV > 0) call dgemm ('N', 'N', nOO+nVV,nVV,nVV, 1d0, Z1, nOO+nVV, O1, nVV,0d0, tmp1, nOO+nVV) if(nVV > 0) call dgemm ('N', 'N', nOO+nVV,nVV,nVV, 1d0, Z1, nOO+nVV, O1, nVV,0d0, tmp1, nOO+nVV)
Z1 = tmp1 Z1 = tmp1

118
src/utils/orthogonalize_matrix.f90 Normal file
View File

@ -0,0 +1,118 @@
subroutine orthogonalize_matrix(ortho_type,nBas,S,X)
! Compute the orthogonalization matrix X
implicit none
! Input variables
integer,intent(in) :: nBas,ortho_type
double precision,intent(in) :: S(nBas,nBas)
! Local variables
logical :: debug
double precision,allocatable :: UVec(:,:),Uval(:)
double precision,parameter :: thresh = 1d-6
integer :: i
! Output variables
double precision,intent(out) :: X(nBas,nBas)
debug = .false.
! Type of orthogonalization ortho_type
!
! 1 = Lowdin
! 2 = Canonical
! 3 = SVD
!
allocate(Uvec(nBas,nBas),Uval(nBas))
if(ortho_type == 1) then
! write(*,*)
! write(*,*) ' Lowdin orthogonalization'
! write(*,*)
Uvec = S
call diagonalize_matrix(nBas,Uvec,Uval)
do i=1,nBas
if(Uval(i) < thresh) then
write(*,*) 'Eigenvalue',i,' is very small in Lowdin orthogonalization = ',Uval(i)
endif
Uval(i) = 1d0/sqrt(Uval(i))
enddo
call ADAt(nBas,Uvec,Uval,X)
elseif(ortho_type == 2) then
! write(*,*)
! write(*,*) 'Canonical orthogonalization'
! write(*,*)
Uvec = S
call diagonalize_matrix(nBas,Uvec,Uval)
do i=1,nBas
if(Uval(i) > thresh) then
Uval(i) = 1d0/sqrt(Uval(i))
else
write(*,*) ' Eigenvalue',i,'too small for canonical orthogonalization'
endif
enddo
call AD(nBas,Uvec,Uval)
X = Uvec
elseif(ortho_type == 3) then
! write(*,*)
! write(*,*) ' SVD-based orthogonalization NYI'
! write(*,*)
! Uvec = S
! call diagonalize_matrix(nBas,Uvec,Uval)
! do i=1,nBas
! if(Uval(i) > thresh) then
! Uval(i) = 1d0/sqrt(Uval(i))
! else
! write(*,*) 'Eigenvalue',i,'too small for canonical orthogonalization'
! endif
! enddo
!
! call AD(nBas,Uvec,Uval)
! X = Uvec
endif
! Print results
if(debug) then
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Orthogonalization matrix'
write(*,'(A28)') '----------------------'
call matout(nBas,nBas,X)
write(*,*)
endif
end subroutine