mirror of
https://github.com/pfloos/quack
synced 2024-12-22 04:14:26 +01:00
GG0T0pp
This commit is contained in:
parent
d5a396200e
commit
b6652240d3
255
src/GT/GG0T0pp.f90
Normal file
255
src/GT/GG0T0pp.f90
Normal file
@ -0,0 +1,255 @@
|
||||
subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE, &
|
||||
linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
|
||||
|
||||
! Perform one-shot calculation with a T-matrix self-energy (G0T0)
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
logical,intent(in) :: dotest
|
||||
|
||||
logical,intent(in) :: doACFDT
|
||||
logical,intent(in) :: exchange_kernel
|
||||
logical,intent(in) :: doXBS
|
||||
logical,intent(in) :: dophBSE
|
||||
logical,intent(in) :: doppBSE
|
||||
logical,intent(in) :: TDA_T
|
||||
logical,intent(in) :: TDA
|
||||
logical,intent(in) :: dBSE
|
||||
logical,intent(in) :: dTDA
|
||||
logical,intent(in) :: linearize
|
||||
double precision,intent(in) :: eta
|
||||
logical,intent(in) :: regularize
|
||||
|
||||
integer,intent(in) :: nOrb
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nS
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: EGHF
|
||||
double precision,intent(in) :: eHF(nOrb)
|
||||
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
|
||||
double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart)
|
||||
|
||||
! Local variables
|
||||
|
||||
logical :: print_T = .false.
|
||||
|
||||
integer :: nOO
|
||||
integer :: nVV
|
||||
double precision :: EcRPA
|
||||
double precision :: EcBSE
|
||||
double precision :: EcGM
|
||||
double precision,allocatable :: Bpp(:,:)
|
||||
double precision,allocatable :: Cpp(:,:)
|
||||
double precision,allocatable :: Dpp(:,:)
|
||||
double precision,allocatable :: Om1(:)
|
||||
double precision,allocatable :: X1(:,:)
|
||||
double precision,allocatable :: Y1(:,:)
|
||||
double precision,allocatable :: rho1(:,:,:)
|
||||
double precision,allocatable :: Om2(:)
|
||||
double precision,allocatable :: X2(:,:)
|
||||
double precision,allocatable :: Y2(:,:)
|
||||
double precision,allocatable :: rho2(:,:,:)
|
||||
double precision,allocatable :: Sig(:)
|
||||
double precision,allocatable :: Z(:)
|
||||
double precision,allocatable :: eGT(:)
|
||||
double precision,allocatable :: eGTlin(:)
|
||||
double precision, allocatable :: Om(:), R(:,:)
|
||||
|
||||
|
||||
! Output variables
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'**********************************'
|
||||
write(*,*)'* Generalized G0T0pp Calculation *'
|
||||
write(*,*)'**********************************'
|
||||
write(*,*)
|
||||
|
||||
|
||||
! TDA for T
|
||||
|
||||
if(TDA_T) then
|
||||
write(*,*) 'Tamm-Dancoff approximation activated for pp T-matrix!'
|
||||
write(*,*)
|
||||
end if
|
||||
|
||||
! TDA
|
||||
|
||||
if(TDA) then
|
||||
write(*,*) 'Tamm-Dancoff approximation activated!'
|
||||
write(*,*)
|
||||
end if
|
||||
|
||||
! Dimensions of the pp-RPA linear reponse matrices
|
||||
|
||||
nOO = nO*(nO - 1)/2
|
||||
nVV = nV*(nV - 1)/2
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), &
|
||||
rho1(nOrb,nOrb,nVV),rho2(nOrb,nOrb,nOO),Sig(nOrb),Z(nOrb),eGT(nOrb),eGTlin(nOrb))
|
||||
|
||||
!----------------------------------------------
|
||||
! Compute T-matrix
|
||||
!----------------------------------------------
|
||||
|
||||
! Compute linear response
|
||||
|
||||
allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO))
|
||||
|
||||
call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eHF,ERI,Cpp)
|
||||
call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eHF,ERI,Dpp)
|
||||
if(.not.TDA_T) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
|
||||
|
||||
call ppGLR(TDA_T,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
|
||||
|
||||
deallocate(Bpp,Cpp,Dpp)
|
||||
|
||||
if(print_T) call print_excitation_energies('ppRPA@GHF','2p',nVV,Om1)
|
||||
if(print_T) call print_excitation_energies('ppRPA@FHF','2h',nOO,Om2)
|
||||
|
||||
!----------------------------------------------
|
||||
! Compute excitation densities
|
||||
!----------------------------------------------
|
||||
|
||||
call GGTpp_excitation_density(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2)
|
||||
|
||||
!----------------------------------------------
|
||||
! Compute T-matrix version of the self-energy
|
||||
!----------------------------------------------
|
||||
|
||||
if(regularize) call GTpp_regularization(nOrb,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2)
|
||||
|
||||
call GGTpp_self_energy_diag(eta,nOrb,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2,EcGM,Sig,Z)
|
||||
|
||||
!----------------------------------------------
|
||||
! Solve the quasi-particle equation
|
||||
!----------------------------------------------
|
||||
|
||||
eGTlin(:) = eHF(:) + Z(:)*Sig(:)
|
||||
|
||||
if(linearize) then
|
||||
|
||||
write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
|
||||
write(*,*)
|
||||
|
||||
eGT(:) = eGTlin(:)
|
||||
|
||||
else
|
||||
|
||||
write(*,*) ' *** Quasiparticle energies obtained by root search *** '
|
||||
write(*,*)
|
||||
|
||||
call GGTpp_QP_graph(eta,nOrb,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2,eGTlin,eHF,eGT,Z)
|
||||
|
||||
end if
|
||||
|
||||
! call GGTpp_plot_self_energy(nOrb,nC,nO,nV,nR,nOO,nVV,eHF,eGT,Om1,rho1,Om2,rho2)
|
||||
|
||||
!----------------------------------------------
|
||||
! Dump results
|
||||
!----------------------------------------------
|
||||
|
||||
! Compute the ppRPA correlation energy
|
||||
|
||||
allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO))
|
||||
|
||||
call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eGT,ERI,Cpp)
|
||||
call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eGT,ERI,Dpp)
|
||||
if(.not.TDA_T) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
|
||||
|
||||
call ppGLR(TDA_T,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
|
||||
|
||||
deallocate(Bpp,Cpp,Dpp)
|
||||
|
||||
call print_GG0T0pp(nOrb,nO,eHF,ENuc,EGHF,Sig,Z,eGT,EcGM,EcRPA)
|
||||
|
||||
! Perform BSE calculation
|
||||
|
||||
! if(dophBSE) then
|
||||
|
||||
! call GGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV, &
|
||||
! Om1,X1,Y1,Om2,X2,Y2,rho1,rho2,ERI,dipole_int,eHF,eGT,EcBSE)
|
||||
|
||||
! if(exchange_kernel) then
|
||||
|
||||
! EcBSE = 0.5d0*EcBSE
|
||||
|
||||
! end if
|
||||
|
||||
! write(*,*)
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp@GHF correlation energy = ',EcBSE,' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp@GHF total energy = ',ENuc + EGHF + EcBSE,' au'
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,*)
|
||||
|
||||
! Compute the BSE correlation energy via the adiabatic connection
|
||||
|
||||
! 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,eta,nOrb,nC,nO,nV,nR,nS, &
|
||||
! nOO,nVV,Om1,X1,Y1,Om2,X2,Y2,rho1,rho2,ERI,eHF,eGT,EcBSE)
|
||||
|
||||
! if(exchange_kernel) then
|
||||
|
||||
! EcBSE = 0.5d0*EcBSE
|
||||
|
||||
! end if
|
||||
|
||||
! write(*,*)
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp@GHF correlation energy = ',EcBSE,' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp@GHF total energy = ',ENuc + EGHF + EcBSE,' au'
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,*)
|
||||
|
||||
! end if
|
||||
|
||||
! end if
|
||||
|
||||
! if(doppBSE) then
|
||||
|
||||
! call GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV, &
|
||||
! Om1,X1,Y1,Om2,X2,Y2,rho1,rho2,ERI,dipole_int,eHF,eGT,EcBSE)
|
||||
|
||||
! write(*,*)
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp@GHF correlation energy = ',EcBSE,' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp@GHF total energy = ',ENuc + EGHF + EcBSE,' au'
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,*)
|
||||
|
||||
! end if
|
||||
|
||||
! Testing zone
|
||||
|
||||
if(dotest) then
|
||||
|
||||
call dump_test_value('G','G0T0pp correlation energy',EcRPA)
|
||||
call dump_test_value('G','G0T0pp HOMO energy',eGT(nO))
|
||||
call dump_test_value('G','G0T0pp LUMO energy',eGT(nO+1))
|
||||
|
||||
end if
|
||||
|
||||
end subroutine
|
||||
87
src/GT/GGTpp_QP_graph.f90
Normal file
87
src/GT/GGTpp_QP_graph.f90
Normal file
@ -0,0 +1,87 @@
|
||||
subroutine GGTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2,eGTlin,eOld,eGT,Z)
|
||||
|
||||
! Compute the graphical solution of the QP equation
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nOO
|
||||
integer,intent(in) :: nVV
|
||||
|
||||
double precision,intent(in) :: eta
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: Om1(nVV)
|
||||
double precision,intent(in) :: rho1(nBas,nBas,nVV)
|
||||
double precision,intent(in) :: Om2(nOO)
|
||||
double precision,intent(in) :: rho2(nBas,nBas,nOO)
|
||||
|
||||
double precision,intent(in) :: eGTlin(nBas)
|
||||
double precision,intent(in) :: eOld(nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: p
|
||||
integer :: nIt
|
||||
integer,parameter :: maxIt = 64
|
||||
double precision,parameter :: thresh = 1d-6
|
||||
double precision,external :: GGTpp_SigC,GGTpp_dSigC
|
||||
double precision :: SigC,dSigC
|
||||
double precision :: f,df
|
||||
double precision :: w
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: eGT(nBas)
|
||||
double precision,intent(out) :: Z(nBas)
|
||||
|
||||
! Run Newton's algorithm to find the root
|
||||
|
||||
write(*,*)'-----------------------------------------------------'
|
||||
write(*,'(A5,1X,A3,1X,A15,1X,A15,1X,A10)') 'Orb.','It.','e_GTpplin (eV)','e_GTpplin (eV)','Z'
|
||||
write(*,*)'-----------------------------------------------------'
|
||||
|
||||
do p=nC+1,nBas-nR
|
||||
|
||||
w = eGTlin(p)
|
||||
nIt = 0
|
||||
f = 1d0
|
||||
|
||||
do while (abs(f) > thresh .and. nIt < maxIt)
|
||||
|
||||
nIt = nIt + 1
|
||||
|
||||
SigC = GGTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,eOld,Om1,rho1,Om2,rho2)
|
||||
dSigC = GGTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,eOld,Om1,rho1,Om2,rho2)
|
||||
f = w - eHF(p) - SigC
|
||||
df = 1d0/(1d0 - dSigC)
|
||||
w = w - df*f
|
||||
|
||||
end do
|
||||
|
||||
if(nIt == maxIt) then
|
||||
|
||||
eGT(p) = eGTlin(p)
|
||||
write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6,1X,A12)') p,nIt,eGTlin(p)*HaToeV,eGT(p)*HaToeV,Z(p),'Cvg Failed!'
|
||||
|
||||
else
|
||||
|
||||
eGT(p) = w
|
||||
Z(p) = df
|
||||
|
||||
write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6)') p,nIt,eGTlin(p)*HaToeV,eGT(p)*HaToeV,Z(p)
|
||||
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
write(*,*)'-----------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
end subroutine
|
||||
62
src/GT/GGTpp_SigC.f90
Normal file
62
src/GT/GGTpp_SigC.f90
Normal file
@ -0,0 +1,62 @@
|
||||
double precision function GGTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2)
|
||||
|
||||
! Compute diagonal of the correlation part of the self-energy
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: p
|
||||
double precision,intent(in) :: w
|
||||
double precision,intent(in) :: eta
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nOO
|
||||
integer,intent(in) :: nVV
|
||||
|
||||
double precision,intent(in) :: e(nBas)
|
||||
double precision,intent(in) :: Om1(nVV)
|
||||
double precision,intent(in) :: rho1(nBas,nBas,nVV)
|
||||
double precision,intent(in) :: Om2(nOO)
|
||||
double precision,intent(in) :: rho2(nBas,nBas,nOO)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,a,cd,kl
|
||||
double precision :: eps
|
||||
|
||||
! Initialize
|
||||
|
||||
GGTpp_SigC = 0d0
|
||||
|
||||
!-------------------------------------------!
|
||||
! Occupied part of the T-matrix self-energy !
|
||||
!-------------------------------------------!
|
||||
|
||||
do i=nC+1,nO
|
||||
|
||||
do cd=1,nVV
|
||||
eps = w + e(i) - Om1(cd)
|
||||
GGTpp_SigC = GGTpp_SigC + rho1(p,i,cd)**2*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
|
||||
end do
|
||||
|
||||
!------------------------------------------!
|
||||
! Virtual part of the T-matrix self-energy !
|
||||
!------------------------------------------!
|
||||
|
||||
do a=nO+1,nBas-nR
|
||||
|
||||
do kl=1,nOO
|
||||
eps = w + e(a) - Om2(kl)
|
||||
GGTpp_SigC = GGTpp_SigC + rho2(p,a,kl)**2*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
|
||||
end do
|
||||
|
||||
end function
|
||||
62
src/GT/GGTpp_dSigC.f90
Normal file
62
src/GT/GGTpp_dSigC.f90
Normal file
@ -0,0 +1,62 @@
|
||||
double precision function GGTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2)
|
||||
|
||||
! Compute diagonal of the correlation part of the self-energy
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: p
|
||||
double precision,intent(in) :: w
|
||||
double precision,intent(in) :: eta
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nOO
|
||||
integer,intent(in) :: nVV
|
||||
|
||||
double precision,intent(in) :: e(nBas)
|
||||
double precision,intent(in) :: Om1(nVV)
|
||||
double precision,intent(in) :: rho1(nBas,nBas,nVV)
|
||||
double precision,intent(in) :: Om2(nOO)
|
||||
double precision,intent(in) :: rho2(nBas,nBas,nOO)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,a,cd,kl
|
||||
double precision :: eps
|
||||
|
||||
! Initialize
|
||||
|
||||
GGTpp_dSigC = 0d0
|
||||
|
||||
!-------------------------------------------!
|
||||
! Occupied part of the T-matrix self-energy !
|
||||
!-------------------------------------------!
|
||||
|
||||
do i=nC+1,nO
|
||||
|
||||
do cd=1,nVV
|
||||
eps = w + e(i) - Om1(cd)
|
||||
GGTpp_dSigC = GGTpp_dSigC - rho1(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
end do
|
||||
|
||||
end do
|
||||
|
||||
!------------------------------------------!
|
||||
! Virtual part of the T-matrix self-energy !
|
||||
!------------------------------------------!
|
||||
|
||||
do a=nO+1,nBas-nR
|
||||
|
||||
do kl=1,nOO
|
||||
eps = w + e(a) - Om2(kl)
|
||||
GGTpp_dSigC = GGTpp_dSigC - rho2(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
end do
|
||||
|
||||
end do
|
||||
|
||||
end function
|
||||
180
src/GT/GGTpp_excitation_density.f90
Normal file
180
src/GT/GGTpp_excitation_density.f90
Normal file
@ -0,0 +1,180 @@
|
||||
subroutine GGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2)
|
||||
|
||||
! Compute excitation densities for T-matrix self-energy
|
||||
|
||||
! TODO
|
||||
! debug DGEMM for nC != 0
|
||||
! and nR != 0
|
||||
|
||||
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: ispin
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
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(nBas,nBas,nVV)
|
||||
double precision,intent(out) :: rho2(nBas,nBas,nOO)
|
||||
|
||||
integer :: dim_1, dim_2
|
||||
double precision, allocatable :: ERI_1(:,:,:)
|
||||
double precision, allocatable :: ERI_2(:,:,:)
|
||||
|
||||
! Initialization
|
||||
|
||||
rho1(:,:,:) = 0d0
|
||||
rho2(:,:,:) = 0d0
|
||||
|
||||
dim_1 = (nBas - nO) * (nBas - nO - 1) / 2
|
||||
dim_2 = nO * (nO - 1) / 2
|
||||
|
||||
if((dim_1 .eq. 0) .or. (dim_2 .eq. 0)) then
|
||||
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(p, q, a, b, ab, c, d, cd, i, j, ij, k, l, kl) &
|
||||
!$OMP SHARED(nC, nBas, nR, nO, rho1, rho2, ERI, X1, Y1, X2, Y2)
|
||||
!$OMP DO COLLAPSE(2)
|
||||
do q = nC+1, nBas-nR
|
||||
do p = nC+1, nBas-nR
|
||||
|
||||
ab = 0
|
||||
|
||||
do a = nO+1, nBas-nR
|
||||
do b = a+1, nBas-nR
|
||||
|
||||
ab = ab + 1
|
||||
|
||||
cd = 0
|
||||
do c = nO+1, nBas-nR
|
||||
do d = c+1, nBas-nR
|
||||
|
||||
cd = cd + 1
|
||||
|
||||
rho1(p,q,ab) = rho1(p,q,ab) &
|
||||
+ (ERI(p,q,c,d) - ERI(p,q,d,c))*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))*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, nBas-nR
|
||||
do d = c+1, nBas-nR
|
||||
|
||||
cd = cd + 1
|
||||
|
||||
rho2(p,q,ij) = rho2(p,q,ij) &
|
||||
+ (ERI(p,q,c,d) - ERI(p,q,d,c))*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))*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
|
||||
|
||||
else
|
||||
|
||||
allocate(ERI_1(nBas,nBas,dim_1), ERI_2(nBas,nBas,dim_2))
|
||||
ERI_1 = 0.d0
|
||||
ERI_2 = 0.d0
|
||||
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(p, q, c, d, cd, k, l, kl) &
|
||||
!$OMP SHARED(nC, nBas, nR, nO, ERI_1, ERI_2, ERI)
|
||||
!$OMP DO COLLAPSE(2)
|
||||
do q = nC+1, nBas-nR
|
||||
do p = nC+1, nBas-nR
|
||||
cd = 0
|
||||
do c = nO+1, nBas-nR
|
||||
do d = c+1, nBas-nR
|
||||
cd = cd + 1
|
||||
ERI_1(p,q,cd) = ERI(p,q,c,d) - ERI(p,q,d,c)
|
||||
enddo
|
||||
enddo
|
||||
kl = 0
|
||||
do k = nC+1, nO
|
||||
do l = k+1, nO
|
||||
kl = kl + 1
|
||||
ERI_2(p,q,kl) = ERI(p,q,k,l) - ERI(p,q,l,k)
|
||||
end do
|
||||
end do
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call dgemm("N", "N", nBas*nBas, dim_1, dim_1, 1.d0, &
|
||||
ERI_1(1,1,1), nBas*nBas, X1(1,1), dim_1, &
|
||||
0.d0, rho1(1,1,1), nBas*nBas)
|
||||
|
||||
call dgemm("N", "N", nBas*nBas, dim_1, dim_2, 1.d0, &
|
||||
ERI_2(1,1,1), nBas*nBas, Y1(1,1), dim_2, &
|
||||
1.d0, rho1(1,1,1), nBas*nBas)
|
||||
|
||||
call dgemm("N", "N", nBas*nBas, dim_2, dim_1, 1.d0, &
|
||||
ERI_1(1,1,1), nBas*nBas, X2(1,1), dim_1, &
|
||||
0.d0, rho2(1,1,1), nBas*nBas)
|
||||
|
||||
call dgemm("N", "N", nBas*nBas, dim_2, dim_2, 1.d0, &
|
||||
ERI_2(1,1,1), nBas*nBas, Y2(1,1), dim_2, &
|
||||
1.d0, rho2(1,1,1), nBas*nBas)
|
||||
|
||||
deallocate(ERI_1, ERI_2)
|
||||
|
||||
endif
|
||||
|
||||
end subroutine
|
||||
105
src/GT/GGTpp_self_energy_diag.f90
Normal file
105
src/GT/GGTpp_self_energy_diag.f90
Normal file
@ -0,0 +1,105 @@
|
||||
subroutine GGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2,EcGM,Sig,Z)
|
||||
|
||||
! Compute diagonal of the correlation part of the T-matrix self-energy
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
double precision,intent(in) :: eta
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nOO
|
||||
integer,intent(in) :: nVV
|
||||
double precision,intent(in) :: e(nBas)
|
||||
double precision,intent(in) :: Om1(nVV)
|
||||
double precision,intent(in) :: rho1(nBas,nBas,nVV)
|
||||
double precision,intent(in) :: Om2(nOO)
|
||||
double precision,intent(in) :: rho2(nBas,nBas,nOO)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,j,a,b,p,cd,kl
|
||||
double precision :: num,eps
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(inout) :: EcGM
|
||||
double precision,intent(inout) :: Sig(nBas)
|
||||
double precision,intent(inout) :: Z(nBas)
|
||||
|
||||
! Initialization
|
||||
|
||||
Sig(:) = 0d0
|
||||
Z(:) = 0d0
|
||||
EcGM = 0d0
|
||||
|
||||
!--------------------------------------!
|
||||
! Occupied part of the Tpp self-energy !
|
||||
!--------------------------------------!
|
||||
|
||||
do p=nC+1,nBas-nR
|
||||
do i=nC+1,nO
|
||||
|
||||
do cd=1,nVV
|
||||
eps = e(p) + e(i) - Om1(cd)
|
||||
num = rho1(p,i,cd)**2
|
||||
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
|
||||
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
!------------------------------------------!
|
||||
! Virtual part of the T-matrix self-energy !
|
||||
!------------------------------------------!
|
||||
|
||||
do p=nC+1,nBas-nR
|
||||
do a=nO+1,nBas-nR
|
||||
|
||||
do kl=1,nOO
|
||||
eps = e(p) + e(a) - Om2(kl)
|
||||
num = rho2(p,a,kl)**2
|
||||
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
|
||||
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
!-------------------------------------!
|
||||
! Galitskii-Migdal correlation energy !
|
||||
!-------------------------------------!
|
||||
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
|
||||
do cd=1,nVV
|
||||
eps = e(i) + e(j) - Om1(cd)
|
||||
num = rho1(i,j,cd)**2
|
||||
EcGM = EcGM + num*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
do a=nO+1,nBas-nR
|
||||
do b=nO+1,nBas-nR
|
||||
|
||||
do kl=1,nOO
|
||||
eps = e(a) + e(b) - Om2(kl)
|
||||
num = rho2(a,b,kl)**2
|
||||
EcGM = EcGM - num*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
Z(:) = 1d0/(1d0 - Z(:))
|
||||
|
||||
end subroutine
|
||||
@ -207,8 +207,8 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0W0 correlation energy =',EcBSE,' au'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0W0 total energy =',ENuc + EGHF + EcBSE,' au'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@GHF correlation energy = ',EcBSE,' au'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@GHF total energy = ',ENuc + EGHF + EcBSE,' au'
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user