mirror of
https://github.com/pfloos/quack
synced 2025-05-06 15:14:55 +02:00
First go for fully complex RG0W0
This commit is contained in:
Loris Burth
parent
cdf6f5ecf2
commit
9b663eeca6
85
src/GW/complex_RGW.f90
Normal file
85
src/GW/complex_RGW.f90
Normal file
@ -0,0 +1,85 @@
|
||||
subroutine complex_RGW(dotest,docG0W0,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,CAP_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
|
||||
|
||||
! Restricted GW module
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
logical,intent(in) :: dotest
|
||||
|
||||
logical,intent(in) :: docG0W0
|
||||
|
||||
integer,intent(in) :: maxSCF
|
||||
integer,intent(in) :: max_diis
|
||||
double precision,intent(in) :: thresh
|
||||
logical,intent(in) :: doACFDT
|
||||
logical,intent(in) :: exchange_kernel
|
||||
logical,intent(in) :: doXBS
|
||||
logical,intent(in) :: dophBSE
|
||||
logical,intent(in) :: dophBSE2
|
||||
logical,intent(in) :: TDA_W
|
||||
logical,intent(in) :: TDA
|
||||
logical,intent(in) :: dBSE
|
||||
logical,intent(in) :: dTDA
|
||||
logical,intent(in) :: doppBSE
|
||||
logical,intent(in) :: singlet
|
||||
logical,intent(in) :: triplet
|
||||
logical,intent(in) :: linearize
|
||||
double precision,intent(in) :: eta
|
||||
logical,intent(in) :: doSRG
|
||||
|
||||
integer,intent(in) :: nNuc
|
||||
double precision,intent(in) :: ZNuc(nNuc)
|
||||
double precision,intent(in) :: rNuc(nNuc,ncart)
|
||||
double precision,intent(in) :: ENuc
|
||||
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nOrb
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nS
|
||||
|
||||
complex*16,intent(in) :: ERHF
|
||||
complex*16,intent(in) :: eHF(nOrb)
|
||||
complex*16,intent(in) :: cHF(nBas,nOrb)
|
||||
complex*16,intent(in) :: PHF(nBas,nBas)
|
||||
double precision,intent(in) :: S(nBas,nBas)
|
||||
double precision,intent(in) :: T(nBas,nBas)
|
||||
double precision,intent(in) :: V(nBas,nBas)
|
||||
double precision,intent(in) :: Hc(nBas,nBas)
|
||||
double precision,intent(in) :: X(nBas,nOrb)
|
||||
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
|
||||
complex*16,intent(in) :: ERI_MO(nOrb,nOrb,nOrb,nOrb)
|
||||
complex*16,intent(in) :: CAP_MO(nOrb,nOrb)
|
||||
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
|
||||
complex*16,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision :: start_GW ,end_GW ,t_GW
|
||||
|
||||
logical :: doccG0W0,doccGW
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Perform cG0W0 calculation
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(docG0W0) then
|
||||
call wall_time(start_GW)
|
||||
call complex_cRG0W0(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,CAP_MO,dipole_int_MO,eHF)
|
||||
call wall_time(end_GW)
|
||||
|
||||
t_GW = end_GW - start_GW
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for complex G0W0 = ',t_GW,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
end subroutine
|
||||
107
src/GW/complex_RGW_QP_graph.f90
Normal file
107
src/GW/complex_RGW_QP_graph.f90
Normal file
@ -0,0 +1,107 @@
|
||||
subroutine complex_RGW_QP_graph(doSRG,eta,flow,nBas,nC,nO,nV,nR,nS,Re_eHF,Im_eHF,&
|
||||
Om,rho,Re_eGWlin,Im_eGWlin, &
|
||||
Re_eOld,Im_eOld,Re_eGW,Im_eGW,Re_Z,Im_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) :: nS
|
||||
|
||||
logical,intent(in) :: doSRG
|
||||
double precision,intent(in) :: eta
|
||||
double precision,intent(in) :: flow
|
||||
double precision,intent(in) :: Re_eHF(nBas)
|
||||
double precision,intent(in) :: Im_eHF(nBas)
|
||||
complex*16,intent(in) :: Om(nS)
|
||||
complex*16,intent(in) :: rho(nBas,nBas,nS)
|
||||
|
||||
double precision,intent(in) :: Re_eGWlin(nBas)
|
||||
double precision,intent(in) :: Im_eGWlin(nBas)
|
||||
double precision,external :: cRGW_Re_SigC,cRGW_Re_dSigC
|
||||
double precision,external :: cRGW_Im_SigC,cRGW_Im_dSigC
|
||||
double precision,intent(in) :: Re_eOld(nBas)
|
||||
double precision,intent(in) :: Im_eOld(nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: p
|
||||
integer :: nIt
|
||||
integer,parameter :: maxIt = 64
|
||||
double precision,parameter :: thresh = 1d-6
|
||||
double precision :: Re_SigC,Re_dSigC
|
||||
double precision :: Im_SigC,Im_dSigC
|
||||
double precision :: Re_f,Im_f,Re_df,Im_df
|
||||
double precision :: Re_w
|
||||
double precision :: Im_w
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: Re_eGW(nBas),Im_eGW(nBas)
|
||||
double precision,intent(out) :: Re_Z(nBas),Im_Z(nBas)
|
||||
|
||||
! Run Newton's algorithm to find the root
|
||||
|
||||
write(*,*)'-----------------------------------------------------'
|
||||
write(*,'(A5,1X,A3,1X,A16,1X,A16,1X,A10)') 'Orb.','It.','Re(e_GWlin) (eV)','Re(e_GW (eV))','Re(Z)'
|
||||
write(*,'(A5,1X,A3,1X,A16,1X,A16,1X,A10)') 'Orb.','It.','Im(e_GWlin) (eV)','Im(e_GW (eV))','Im(Z)'
|
||||
write(*,*)'-----------------------------------------------------'
|
||||
|
||||
do p=nC+1,nBas-nR
|
||||
|
||||
Re_w = Re_eGWlin(p)
|
||||
Im_w = Im_eGWlin(p)
|
||||
nIt = 0
|
||||
Re_f = 1d0
|
||||
Im_f = 1d0
|
||||
|
||||
do while (sqrt(Re_f**2+Im_f**2) > thresh .and. nIt < maxIt)
|
||||
|
||||
nIt = nIt + 1
|
||||
|
||||
|
||||
! Re_SigC = cRGW_Re_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_eOld,Im_eold,Om,rho)
|
||||
! Im_SigC = cRGW_Im_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_eOld,Im_eold,Om,rho)
|
||||
! Re_dSigC = cRGW_Re_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_eOld,Im_eold,Om,rho)
|
||||
! Im_dSigC = cRGW_Im_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_eOld,Im_eold,Om,rho)
|
||||
!
|
||||
|
||||
Re_f = Re_w - Re_eHF(p) - Re_SigC
|
||||
Im_f = Im_w - Im_eHF(p) - Im_SigC
|
||||
Re_df = (1d0 - Re_dSigC)/((1d0 - Re_dSigC)**2 + Im_dSigC**2)
|
||||
Im_df = Im_dSigC/((1d0 - Re_dSigC)**2 + Im_dSigC**2)
|
||||
Re_w = Re_w - Re_df*Re_f + Im_df*Im_f
|
||||
Im_w = Im_w - Re_f*Im_df - Re_df*Im_f
|
||||
|
||||
end do
|
||||
|
||||
if(nIt == maxIt) then
|
||||
|
||||
Re_eGW(p) = Re_eGWlin(p)
|
||||
write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6,1X,A12)') p,nIt,Re_eGWlin(p)*HaToeV,Re_eGW(p)*HaToeV,Re_Z(p),'Cvg Failed!'
|
||||
|
||||
else
|
||||
|
||||
Re_eGW(p) = Re_w
|
||||
Im_eGW(p) = Im_w
|
||||
Re_Z(p) = Re_df
|
||||
Im_Z(p) = Im_df
|
||||
|
||||
write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6)') p,nIt,Re_eGWlin(p)*HaToeV,Re_eGW(p)*HaToeV,Re_Z(p)
|
||||
write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6)') p,nIt,Im_eGWlin(p)*HaToeV,Im_eGW(p)*HaToeV,Im_Z(p)
|
||||
|
||||
end if
|
||||
|
||||
write(*,*)'-----------------------------------------------------'
|
||||
end do
|
||||
write(*,*)
|
||||
|
||||
end subroutine
|
||||
46
src/GW/complex_RGW_excitation_density.f90
Normal file
46
src/GW/complex_RGW_excitation_density.f90
Normal file
@ -0,0 +1,46 @@
|
||||
subroutine complex_RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
|
||||
|
||||
! Compute excitation densities
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nOrb
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nS
|
||||
complex*16,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
|
||||
complex*16,intent(in) :: XpY(nS,nS)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: ia,jb,p,q,j,b
|
||||
|
||||
! Output variables
|
||||
|
||||
complex*16,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,q,b)*XpY(ia,jb)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
end subroutine
|
||||
97
src/GW/complex_RGW_self_energy_diag.f90
Normal file
97
src/GW/complex_RGW_self_energy_diag.f90
Normal file
@ -0,0 +1,97 @@
|
||||
subroutine complex_RGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho,EcGM,Re_Sig,Im_Sig,Re_Z,Im_Z)
|
||||
|
||||
! Compute diagonal of the correlation part of the self-energy and the renormalization factor
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
double precision,intent(in) :: eta
|
||||
integer,intent(in) :: nBas
|
||||
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) :: Re_e(nBas)
|
||||
double precision,intent(in) :: Im_e(nBas)
|
||||
complex*16,intent(in) :: Om(nS)
|
||||
complex*16,intent(in) :: rho(nBas,nBas,nS)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,a,p,m
|
||||
double precision :: eps
|
||||
complex*16 :: num
|
||||
double precision :: eta_tilde
|
||||
double precision,allocatable :: Re_DS(:)
|
||||
double precision,allocatable :: Im_DS(:)
|
||||
complex*16 :: tmp
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: Re_Sig(nBas)
|
||||
double precision,intent(out) :: Im_Sig(nBas)
|
||||
double precision,intent(out) :: Re_Z(nBas)
|
||||
double precision,intent(out) :: Im_Z(nBas)
|
||||
complex*16,intent(out) :: EcGM
|
||||
|
||||
! Initialize
|
||||
allocate(Re_DS(nBas),Im_DS(nBas))
|
||||
Re_Sig(:) = 0d0
|
||||
Im_Sig(:) = 0d0
|
||||
Re_DS(:) = 0d0
|
||||
Im_DS(:) = 0d0
|
||||
|
||||
!----------------!
|
||||
! GW self-energy !
|
||||
!----------------!
|
||||
|
||||
! Occupied part of the correlation self-energy
|
||||
do p=nC+1,nBas-nR
|
||||
do i=nC+1,nO
|
||||
do m=1,nS
|
||||
eps = Re_e(p) - Re_e(i) + real(Om(m))
|
||||
eta_tilde = eta - Im_e(p) + Im_e(i) - aimag(Om(m))
|
||||
num = 2d0*rho(p,i,m)**2
|
||||
tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),eta_tilde/(eps**2+eta_tilde**2),kind=8)
|
||||
Re_Sig(p) = Re_Sig(p) + real(tmp)
|
||||
Im_Sig(p) = Im_Sig(p) + aimag(tmp)
|
||||
tmp = num*cmplx(-(eps**2-eta_tilde**2)/(eps**2 + eta_tilde**2)**2,&
|
||||
-2*eta_tilde*eps/(eps**2 + eta_tilde**2)**2,kind=8)
|
||||
Re_DS(p) = Re_DS(p) + real(tmp)
|
||||
Im_DS(p) = Im_DS(p) + aimag(tmp)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Virtual part of the correlation self-energy
|
||||
|
||||
do p=nC+1,nBas-nR
|
||||
do a=nO+1,nBas-nR
|
||||
do m=1,nS
|
||||
|
||||
eps = Re_e(p) - Re_e(a) - real(Om(m))
|
||||
eta_tilde = eta + Im_e(p) - Im_e(a) - aimag(Om(m))
|
||||
num = 2d0*rho(p,a,m)**2
|
||||
tmp = num*cmplx(eps/(eps**2 + eta_tilde**2)**2,&
|
||||
-eta_tilde/(eps**2 + eta_tilde**2)**2,kind=8)
|
||||
Re_Sig(p) = Re_Sig(p) + real(tmp)
|
||||
Im_Sig(p) = Im_Sig(p) + aimag(tmp)
|
||||
tmp = num*cmplx(-(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2,&
|
||||
2*eta_tilde*eps/eps/(eps**2 + eta_tilde**2)**2,kind=8)
|
||||
Re_DS(p) = Re_DS(p) + real(tmp)
|
||||
Im_DS(p) = Im_DS(p) + aimag(tmp)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
|
||||
! Compute renormalization factor from derivative
|
||||
Re_Z(:) = (1d0-Re_DS(:))/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
|
||||
Im_Z(:) = Im_DS(:)/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
|
||||
deallocate(Re_DS,Im_DS)
|
||||
end subroutine
|
||||
@ -35,11 +35,11 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nS
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: ERHF
|
||||
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
|
||||
double precision,intent(in) :: CAP(nOrb,nOrb)
|
||||
double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart)
|
||||
double precision,intent(in) :: eHF(nOrb)
|
||||
complex*16,intent(in) :: ERHF
|
||||
complex*16,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
|
||||
complex*16,intent(in) :: CAP(nOrb,nOrb)
|
||||
complex*16,intent(in) :: dipole_int(nOrb,nOrb,ncart)
|
||||
complex*16,intent(in) :: eHF(nOrb)
|
||||
|
||||
! Local variables
|
||||
|
||||
@ -49,26 +49,29 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,
|
||||
integer :: isp_W
|
||||
integer :: p
|
||||
double precision :: flow
|
||||
double precision :: EcRPA
|
||||
double precision :: EcBSE(nspin)
|
||||
double precision :: EcGM
|
||||
double precision,allocatable :: Aph(:,:)
|
||||
double precision,allocatable :: Bph(:,:)
|
||||
complex*16 :: EcRPA
|
||||
complex*16 :: EcBSE(nspin)
|
||||
complex*16 :: EcGM
|
||||
complex*16,allocatable :: Aph(:,:)
|
||||
complex*16,allocatable :: Bph(:,:)
|
||||
double precision,allocatable :: Re_SigC(:)
|
||||
double precision,allocatable :: Im_SigC(:)
|
||||
double precision,allocatable :: Re_Z(:)
|
||||
double precision,allocatable :: Im_Z(:)
|
||||
double precision,allocatable :: Om(:)
|
||||
double precision,allocatable :: XpY(:,:)
|
||||
double precision,allocatable :: XmY(:,:)
|
||||
double precision,allocatable :: rho(:,:,:)
|
||||
complex*16,allocatable :: Om(:)
|
||||
double precision,allocatable :: Re_Om(:)
|
||||
double precision,allocatable :: Im_Om(:)
|
||||
complex*16,allocatable :: XpY(:,:)
|
||||
complex*16,allocatable :: XmY(:,:)
|
||||
complex*16,allocatable :: rho(:,:,:)
|
||||
|
||||
|
||||
double precision,allocatable :: Re_eGWlin(:)
|
||||
double precision, allocatable :: Im_eGWlin(:)
|
||||
double precision,allocatable :: Re_eGW(:)
|
||||
double precision,allocatable :: Im_eGW(:)
|
||||
double precision, allocatable :: e_cap(:)
|
||||
double precision, allocatable :: Re_eHF(:)
|
||||
double precision, allocatable :: Im_eHF(:)
|
||||
|
||||
! Hello world
|
||||
|
||||
@ -103,10 +106,9 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,
|
||||
! Memory allocation
|
||||
|
||||
allocate(Aph(nS,nS),Bph(nS,nS),Re_SigC(nOrb),Im_SigC(nOrb),Re_Z(nOrb),Im_Z(nOrb),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS), &
|
||||
Re_eGW(nOrb),Im_eGW(nOrb),Re_eGWlin(nOrb),Im_eGWlin(nOrb),e_cap(nOrb))
|
||||
do p = 1, nOrb
|
||||
e_cap(p) = CAP(p,p)
|
||||
end do
|
||||
Re_eGW(nOrb),Im_eGW(nOrb),Re_eGWlin(nOrb),Im_eGWlin(nOrb),Re_eHF(nOrb),Im_eHF(nOrb),Re_Om(nS),Im_Om(nS))
|
||||
Re_eHF(:) = real(eHF(:))
|
||||
Im_eHF(:) = aimag(eHF(:))
|
||||
!-------------------!
|
||||
! Compute screening !
|
||||
!-------------------!
|
||||
@ -115,6 +117,9 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,
|
||||
if(.not.TDA_W) call complex_phRLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
|
||||
|
||||
call complex_phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
|
||||
Re_Om(:) = real(Om(:))
|
||||
Im_Om(:) = aimag(Om(:))
|
||||
call complex_vecout(nS,Om)
|
||||
|
||||
!if(print_W) call print_excitation_energies('phRPA@RHF','singlet',nS,Om)
|
||||
|
||||
@ -123,22 +128,19 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,
|
||||
!--------------------------!
|
||||
|
||||
call complex_RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
|
||||
|
||||
write(*,*) rho(1,1,1)
|
||||
!------------------------!
|
||||
! Compute GW self-energy !
|
||||
!------------------------!
|
||||
|
||||
!!!! STOPPED HERE PROCEED WITH IMPLEMENTING COMPLEX_CRGW_SELF_ENERGY_DIAG
|
||||
!!!!
|
||||
call complex_cRGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,Re_SigC,Im_SigC,Re_Z,Im_Z,e_cap)
|
||||
|
||||
call complex_RGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_eHF,Im_eHF,Om,rho,EcGM,Re_SigC,Im_SigC,Re_Z,Im_Z)
|
||||
|
||||
!-----------------------------------!
|
||||
! Solve the quasi-particle equation !
|
||||
!-----------------------------------!
|
||||
|
||||
! Linearized or graphical solution?
|
||||
Re_eGWlin(:) = eHF(:) + Re_Z(:)*Re_SigC(:) - Im_Z(:)*Im_SigC(:)
|
||||
Im_eGWlin(:) = e_cap(:) + Re_Z(:)*Im_SigC(:) + Im_Z(:)*Re_SigC(:)
|
||||
Re_eGWlin(:) = Re_eHF(:) + Re_Z(:)*Re_SigC(:) - Im_Z(:)*Im_SigC(:)
|
||||
Im_eGWlin(:) = Im_eHF(:) + Re_Z(:)*Im_SigC(:) + Im_Z(:)*Re_SigC(:)
|
||||
|
||||
if(linearize) then
|
||||
|
||||
@ -152,94 +154,96 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,
|
||||
|
||||
write(*,*) ' *** Quasiparticle energies obtained by root search *** '
|
||||
write(*,*)
|
||||
|
||||
call cRGW_QP_graph(doSRG,eta,flow,nOrb,nC,nO,nV,nR,nS,eHF,e_cap,Om,rho,Re_eGWlin,Im_eGWlin,eHF,e_cap,Re_eGW,Im_eGW,Re_Z,Im_Z)
|
||||
write(*,*) "ROOT SEARCH NOT IMPLEMENTED YET"
|
||||
! call complex_RGW_QP_graph(doSRG,eta,flow,nOrb,nC,nO,nV,nR,nS, &
|
||||
! Re_eHF,Im_eHF,Om,rho,Re_eGWlin,Im_eGWlin,Re_eHF,Im_eHF, &
|
||||
! Re_eGW,Im_eGW,Re_Z,Im_Z)
|
||||
end if
|
||||
|
||||
! Plot self-energy, renormalization factor, and spectral function
|
||||
|
||||
if(plot_self) call RGW_plot_self_energy(nOrb,eta,nC,nO,nV,nR,nS,eHF,eHF,Om,rho)
|
||||
|
||||
! Cumulant expansion
|
||||
|
||||
! call RGWC(dotest,eta,nOrb,nC,nO,nV,nR,nS,Om,rho,eHF,eHF,eGW,Z)
|
||||
|
||||
! Compute the RPA correlation energy
|
||||
|
||||
call phRLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,Re_eGW,ERI,Aph)
|
||||
if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
|
||||
|
||||
call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
|
||||
|
||||
!--------------!
|
||||
! Dump results !
|
||||
!--------------!
|
||||
|
||||
call print_cRG0W0(nOrb,nO,eHF,ENuc,ERHF,Re_SigC,Im_SigC,Re_Z,Im_Z,Re_eGW,Im_eGW,EcRPA,EcGM,CAP)
|
||||
!---------------------------!
|
||||
! Perform phBSE calculation !
|
||||
!---------------------------!
|
||||
!
|
||||
! if(dophBSE) then
|
||||
! if(plot_self) call RGW_plot_self_energy(nOrb,eta,nC,nO,nV,nR,nS,eHF,eHF,Om,rho)
|
||||
!
|
||||
! call RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta, &
|
||||
! nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,Re_eGW,EcBSE)
|
||||
!! Cumulant expansion
|
||||
!
|
||||
! write(*,*)
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au'
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,*)
|
||||
!! call RGWC(dotest,eta,nOrb,nC,nO,nV,nR,nS,Om,rho,eHF,eHF,eGW,Z)
|
||||
!
|
||||
! ! Compute the BSE correlation energy via the adiabatic connection fluctuation dissipation theorem
|
||||
!! Compute the RPA correlation energy
|
||||
!
|
||||
! if(doACFDT) then
|
||||
! call phRLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,Re_eGW,ERI,Aph)
|
||||
! if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
|
||||
!
|
||||
! call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,eHF,Re_eGW,EcBSE)
|
||||
! call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
|
||||
!
|
||||
! write(*,*)
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au'
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,*)
|
||||
!
|
||||
! end if
|
||||
!
|
||||
! end if
|
||||
!!--------------!
|
||||
!! Dump results !
|
||||
!!--------------!
|
||||
!
|
||||
call print_complex_cRG0W0(nOrb,nO,Re_eHF,Im_eHF,ENuc,ERHF,Re_SigC,Im_SigC,Re_Z,Im_Z,Re_eGW,Im_eGW,EcRPA,EcGM)
|
||||
!!---------------------------!
|
||||
!! Perform ppBSE calculation !
|
||||
!! Perform phBSE calculation !
|
||||
!!---------------------------!
|
||||
!
|
||||
! if(doppBSE) then
|
||||
!
|
||||
! call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,Re_eGW,EcBSE)
|
||||
!
|
||||
! write(*,*)
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au'
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,*)
|
||||
!
|
||||
! end if
|
||||
!
|
||||
!! Testing zone
|
||||
!
|
||||
! if(dotest) then
|
||||
!
|
||||
! call dump_test_value('R','G0W0 correlation energy',EcRPA)
|
||||
! call dump_test_value('R','G0W0 HOMO energy',Re_eGW(nO))
|
||||
! call dump_test_value('R','G0W0 LUMO energy',Re_eGW(nO+1))
|
||||
!
|
||||
! end if
|
||||
!
|
||||
!!
|
||||
!! if(dophBSE) then
|
||||
!!
|
||||
!! call RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta, &
|
||||
!! nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,Re_eGW,EcBSE)
|
||||
!!
|
||||
!! write(*,*)
|
||||
!! write(*,*)'-------------------------------------------------------------------------------'
|
||||
!! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au'
|
||||
!! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au'
|
||||
!! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au'
|
||||
!! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au'
|
||||
!! write(*,*)'-------------------------------------------------------------------------------'
|
||||
!! write(*,*)
|
||||
!!
|
||||
!! ! Compute the BSE correlation energy via the adiabatic connection fluctuation dissipation theorem
|
||||
!!
|
||||
!! if(doACFDT) then
|
||||
!!
|
||||
!! call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,eHF,Re_eGW,EcBSE)
|
||||
!!
|
||||
!! write(*,*)
|
||||
!! write(*,*)'-------------------------------------------------------------------------------'
|
||||
!! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au'
|
||||
!! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au'
|
||||
!! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au'
|
||||
!! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au'
|
||||
!! write(*,*)'-------------------------------------------------------------------------------'
|
||||
!! write(*,*)
|
||||
!!
|
||||
!! end if
|
||||
!!
|
||||
!! end if
|
||||
!!
|
||||
!!!---------------------------!
|
||||
!!! Perform ppBSE calculation !
|
||||
!!!---------------------------!
|
||||
!!
|
||||
!! if(doppBSE) then
|
||||
!!
|
||||
!! call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,Re_eGW,EcBSE)
|
||||
!!
|
||||
!! write(*,*)
|
||||
!! write(*,*)'-------------------------------------------------------------------------------'
|
||||
!! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au'
|
||||
!! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au'
|
||||
!! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au'
|
||||
!! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au'
|
||||
!! write(*,*)'-------------------------------------------------------------------------------'
|
||||
!! write(*,*)
|
||||
!!
|
||||
!! end if
|
||||
!!
|
||||
!!! Testing zone
|
||||
!!
|
||||
!! if(dotest) then
|
||||
!!
|
||||
!! call dump_test_value('R','G0W0 correlation energy',EcRPA)
|
||||
!! call dump_test_value('R','G0W0 HOMO energy',Re_eGW(nO))
|
||||
!! call dump_test_value('R','G0W0 LUMO energy',Re_eGW(nO+1))
|
||||
!!
|
||||
!! end if
|
||||
!!
|
||||
end subroutine
|
||||
|
||||
@ -1,80 +0,0 @@
|
||||
subroutine RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
|
||||
|
||||
! Compute excitation densities
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nOrb
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nS
|
||||
complex*16,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
|
||||
complex*16,intent(in) :: XpY(nS,nS)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: ia,jb,p,q,j,b
|
||||
complex*16, allocatable :: tmp(:,:,:)
|
||||
|
||||
! Output variables
|
||||
|
||||
complex*16,intent(out) :: rho(nOrb,nOrb,nS)
|
||||
|
||||
if(nOrb .lt. 256) then
|
||||
|
||||
allocate(tmp(nOrb,nOrb,nS))
|
||||
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(p, q, j, b, jb) &
|
||||
!$OMP SHARED(nOrb, nC, nO, nR, ERI, tmp)
|
||||
!$OMP DO COLLAPSE(2)
|
||||
do p = 1, nOrb
|
||||
do q = 1, nOrb
|
||||
jb = 0
|
||||
do j = nC+1, nO
|
||||
do b = nO+1, nOrb-nR
|
||||
jb = jb + 1
|
||||
tmp(p,q,jb) = ERI(p,j,q,b)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call zgemm("N", "T", nOrb*nOrb, nS, nS, 1.d0, &
|
||||
tmp(1,1,1), nOrb*nOrb, XpY(1,1), nS, &
|
||||
0.d0, rho(1,1,1), nOrb*nOrb)
|
||||
|
||||
deallocate(tmp)
|
||||
|
||||
else
|
||||
|
||||
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,q,b)*XpY(ia,jb)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
endif
|
||||
|
||||
end subroutine
|
||||
66
src/GW/print_complex_cRG0W0.f90
Normal file
66
src/GW/print_complex_cRG0W0.f90
Normal file
@ -0,0 +1,66 @@
|
||||
subroutine print_complex_cRG0W0(nBas,nO,Re_eHF,Im_eHF,ENuc,ERHF,Re_SigC,Im_SigC,Re_Z,Im_Z,Re_eGW,Im_eGW,EcRPA,EcGM)
|
||||
|
||||
! Print one-electron energies and other stuff for G0W0
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
integer,intent(in) :: nBas,nO
|
||||
double precision,intent(in) :: ENuc
|
||||
complex*16,intent(in) :: ERHF
|
||||
complex*16,intent(in) :: EcRPA
|
||||
complex*16,intent(in) :: EcGM
|
||||
double precision,intent(in) :: Re_eHF(nBas)
|
||||
double precision,intent(in) :: Im_eHF(nBas)
|
||||
double precision,intent(in) :: Re_SigC(nBas)
|
||||
double precision,intent(in) :: Im_SigC(nBas)
|
||||
double precision,intent(in) :: Re_Z(nBas)
|
||||
double precision,intent(in) :: Im_Z(nBas)
|
||||
double precision,intent(in) :: Re_eGW(nBas)
|
||||
double precision,intent(in) :: Im_eGW(nBas)
|
||||
|
||||
integer :: p,index_homo,index_lumo
|
||||
double precision :: Re_eHOMO,Re_eLUMO,Im_eHOMO,Im_eLUMO,Re_Gap,Im_Gap
|
||||
|
||||
! HOMO and LUMO
|
||||
|
||||
index_homo = maxloc(Re_eGW(1:nO),1)
|
||||
Re_eHOMO = Re_eGW(index_homo)
|
||||
Im_eHOMO = Im_eGW(index_homo)
|
||||
index_lumo = minloc(Re_eGW(nO+1:nBas),1) + nO
|
||||
Re_eLUMO = Re_eGW(index_lumo)
|
||||
Im_eLUMO = Im_eGW(index_lumo)
|
||||
Re_Gap = Re_eLUMO-Re_eHOMO
|
||||
Im_Gap = Im_eLUMO-Im_eHOMO
|
||||
|
||||
! Dump results
|
||||
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)' cG0W0@RHF calculation '
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
|
||||
'|','#','|','Re(e_HF) (eV)','|','Re(Sig_GW) (eV)','|','Re(Z)','|','Re(e_GW) (eV)','|'
|
||||
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
|
||||
'|','#','|','Im(e_HF) (eV)','|','Im(Sig_GW) (eV)','|','Im(Z)','|','Im(e_GW) (eV)','|'
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
do p=1,nBas
|
||||
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
|
||||
'|',p,'|',Re_eHF(p)*HaToeV,'|',Re_SigC(p)*HaToeV,'|',Re_Z(p),'|',Re_eGW(p)*HaToeV,'|'
|
||||
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
|
||||
'|',p,'|',Im_eHF(p)*HaToeV,'|',Im_SigC(p)*HaToeV,'|',Im_Z(p),'|',Im_eGW(p)*HaToeV,'|'
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
end do
|
||||
write(*,*)
|
||||
write(*,*)'---------------------------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A60,F15.6,A5,F15.6,A3)') 'cG0W0@RHF HOMO energy = ',Re_eHOMO*HaToeV,' + i*',Im_eHOMO*HaToeV,' eV'
|
||||
write(*,'(2X,A60,F15.6,A5,F15.6,A3)') 'cG0W0@RHF LUMO energy = ',Re_eLUMO*HaToeV,' + i*',Im_eLUMO*HaToeV,' eV'
|
||||
write(*,'(2X,A60,F15.6,A5,F15.6,A3)') 'cG0W0@RHF HOMO-LUMO gap = ',Re_Gap*HaToeV,' + i*',Im_Gap*HaToeV,' eV'
|
||||
write(*,*)'---------------------------------------------------------------------------------------------------'
|
||||
! write(*,'(2X,A60,F15.6,A3)') 'phRPA@cG0W0@RHF total energy = ',ENuc + ERHF + EcRPA,' au'
|
||||
! write(*,'(2X,A60,F15.6,A3)') 'phRPA@cG0W0@RHF correlation energy = ',EcRPA,' au'
|
||||
! write(*,'(2X,A60,F15.6,A3)') ' GM@cG0W0@RHF total energy = ',ENuc + ERHF + EcGM,' au'
|
||||
! write(*,'(2X,A60,F15.6,A3)') ' GM@cG0W0@RHF correlation energy = ',EcGM,' au'
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,*)
|
||||
!
|
||||
end subroutine
|
||||
@ -1,4 +1,4 @@
|
||||
subroutine phRLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
|
||||
subroutine complex_phRLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
|
||||
|
||||
! Compute linear response
|
||||
|
||||
@ -14,7 +14,7 @@ subroutine phRLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
|
||||
|
||||
! Local variables
|
||||
|
||||
complex*16 :: complex_trace_matrix
|
||||
complex*16,external :: complex_trace_matrix
|
||||
complex*16 :: t1, t2
|
||||
complex*16,allocatable :: RPA_matrix(:,:)
|
||||
complex*16,allocatable :: Z(:,:)
|
||||
@ -42,20 +42,22 @@ subroutine phRLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
|
||||
else
|
||||
|
||||
allocate(RPA_matrix(2*nS,2*nS),OmOmminus(2*nS))
|
||||
|
||||
call complex_diagonalize_matrix(2*nS,RPA_matrix,Om)
|
||||
RPA_matrix(1:nS,1:nS) = Aph(:,:)
|
||||
RPA_matrix(1:nS,nS+1:2*nS) = Bph(:,:)
|
||||
RPA_matrix(nS+1:2*nS,1:nS) = -Bph(:,:)
|
||||
RPA_matrix(nS+1:2*nS,nS+1:2*nS) = -Aph(:,:)
|
||||
call complex_diagonalize_matrix_without_sort(2*nS,RPA_matrix,OmOmminus)
|
||||
call complex_sort_eigenvalues_RPA(2*nS,OmOmminus,RPA_matrix)
|
||||
call complex_normalize_RPA(nS,RPA_matrix)
|
||||
Om(:) = OmOmminus(1:nS)
|
||||
call complex_vecout(nS,Om)
|
||||
call complex_vecout(nS,Om(nS+1:2*nS))
|
||||
if (maxval(abs(Om(:)+Om(nS+1:2*nS))) > 1e10) &
|
||||
if(maxval(abs(OmOmminus(1:nS)+OmOmminus(nS+1:2*nS))) > 1e-12) &
|
||||
call print_warning('We dont find a Om and -Om structure as solution of the RPA. There might be a problem somewhere.')
|
||||
if(minval(abs(Om(:))) < 0d0) &
|
||||
call print_warning('You may have instabilities in linear response: A-B is not positive definite!!')
|
||||
XpY(:,:) = RPA_matrix(1:nS,1:nS) + RPA_matrix(nS+1:2*nS,nS+1:2*nS)
|
||||
XmY(:,:) = RPA_matrix(1:nS,1:nS) - RPA_matrix(nS+1:2*nS,nS+1:2*nS)
|
||||
|
||||
deallocate(RPA_matrix)
|
||||
|
||||
XpY(:,:) = RPA_matrix(1:nS,1:nS) + RPA_matrix(1:nS,nS+1:2*nS)
|
||||
XmY(:,:) = RPA_matrix(1:nS,1:nS) - RPA_matrix(1:nS,nS+1:2*nS)
|
||||
call complex_matout(nS,nS,XpY)
|
||||
deallocate(RPA_matrix,OmOmminus)
|
||||
end if
|
||||
|
||||
! Compute the RPA correlation energy
|
||||
|
||||
@ -363,7 +363,7 @@ doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3
|
||||
call wall_time(start_GF)
|
||||
call RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2,renorm_GF,maxSCF_GF, &
|
||||
thresh_GF,max_diis_GF,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,lin_GF, &
|
||||
eta_GF,reg_GF,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF, &
|
||||
eta_GF,reg_GF,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,complex_ERHF, &
|
||||
S,X,T,V,Hc,ERI_AO,ERI_MO,CAP_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
|
||||
call wall_time(end_GF)
|
||||
|
||||
@ -396,7 +396,7 @@ doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3
|
||||
!-----------!
|
||||
|
||||
doGW = doG0W0 .or. doevGW .or. doqsGW .or. doufG0W0 .or. doufGW .or. docG0W0
|
||||
if(doGW) then
|
||||
if(doGW .and. .not. docRHF) then
|
||||
|
||||
call wall_time(start_GW)
|
||||
call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,docG0W0,maxSCF_GW,thresh_GW,max_diis_GW, &
|
||||
@ -416,8 +416,18 @@ doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3
|
||||
! complex GW module !
|
||||
!-------------------!
|
||||
|
||||
! IMPLEMENT LATER TO TREAT EVERYTHING COMPLEX, i.e. start from complex HF orbitals
|
||||
if(doGW .and. docRHF) then
|
||||
call wall_time(start_GW)
|
||||
call complex_RGW(dotest,docG0W0,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,complex_ERI_MO,complex_CAP_MO,dipole_int_AO,complex_dipole_int_MO,complex_PHF,complex_cHF,complex_eHF)
|
||||
call wall_time(end_GW)
|
||||
|
||||
t_GW = end_GW - start_GW
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GW = ',t_GW,' seconds'
|
||||
write(*,*)
|
||||
end if
|
||||
!-----------------!
|
||||
! T-matrix module !
|
||||
!-----------------!
|
||||
|
||||
@ -19,6 +19,53 @@ subroutine complex_orthogonalize_matrix(N, vectors)
|
||||
vectors = matmul(vectors,transpose(Linv))
|
||||
deallocate(L,Linv)
|
||||
end subroutine
|
||||
subroutine complex_orthonormalize(N,vectors,A)
|
||||
|
||||
! Orthonormalize vectors matrix, such that vectors^T A vectors = Identity
|
||||
! A and vectors are assumed quadratic NxN matrices
|
||||
|
||||
! Input variables
|
||||
implicit none
|
||||
integer, intent(in) :: N
|
||||
complex*16, intent(inout) :: vectors(N, N)
|
||||
complex*16, intent(inout) :: A(N, N)
|
||||
|
||||
! Local variables
|
||||
integer :: i, j
|
||||
complex*16,allocatable :: L(:,:),Linv(:,:)
|
||||
complex*16 :: proj
|
||||
complex*16 :: norm
|
||||
|
||||
! Copy the input matrix to a temporary matrix
|
||||
allocate(L(N,N),Linv(N,N))
|
||||
L = matmul(matmul(transpose(vectors),A),vectors)
|
||||
call complex_cholesky_decomp(N,L)
|
||||
call complex_inverse_matrix(N,L,Linv)
|
||||
vectors = matmul(vectors,transpose(Linv))
|
||||
deallocate(L,Linv)
|
||||
end subroutine
|
||||
subroutine complex_normalize_RPA(nS,XYYX)
|
||||
|
||||
! Orthonormalize vectors matrix, such that RPA^T (1 0; 0 -1) RPA = Identity
|
||||
|
||||
! Input variables
|
||||
implicit none
|
||||
integer, intent(in) :: nS
|
||||
complex*16, intent(inout) :: XYYX(2*nS, 2*nS)
|
||||
|
||||
! Local variables
|
||||
integer :: i
|
||||
complex*16,allocatable :: A(:,:)
|
||||
|
||||
allocate(A(2*nS,2*nS))
|
||||
A(:,:) = (0d0,0d0)
|
||||
do i=1,nS
|
||||
A(i,i) = 1
|
||||
A(i+nS,i+nS) = -1
|
||||
end do
|
||||
call complex_orthonormalize(2*nS,XYYX,A)
|
||||
deallocate(A)
|
||||
end subroutine
|
||||
subroutine complex_gram_schmidt(N, vectors)
|
||||
|
||||
! Input variables
|
||||
|
||||
36
src/utils/complex_sort_eigenvalues_RPA.f90
Normal file
36
src/utils/complex_sort_eigenvalues_RPA.f90
Normal file
@ -0,0 +1,36 @@
|
||||
subroutine complex_sort_eigenvalues_RPA(n, evals, evecs)
|
||||
|
||||
! Sorts the eigenvecs and eigenvals like in the rpa problem, i.e. first the positive ones in ascending order than the negative ones
|
||||
! in descending order
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: n
|
||||
complex*16, intent(inout) :: evals(n)
|
||||
complex*16, intent(inout) :: evecs(n, n)
|
||||
integer :: i, j, k
|
||||
complex*16 :: temp_val
|
||||
complex*16 :: temp_vec(n)
|
||||
|
||||
! bubble sort (or any sorting algorithm can be applied)
|
||||
do i = 1, n-1
|
||||
do j = i+1, n
|
||||
! sorting condition: first ascending for positives, then by absolute value for negatives
|
||||
if ((real(evals(i)) > 0 .and. real(evals(j)) > 0 .and. real(evals(i)) > real(evals(j))) .or. &
|
||||
(real(evals(i)) < 0 .and. real(evals(j)) < 0 .and. abs(real(evals(i))) > abs(real(evals(j)))) .or. &
|
||||
(real(evals(i)) < 0 .and. real(evals(j)) > 0)) then
|
||||
|
||||
! swap eigenvalues
|
||||
temp_val = evals(i)
|
||||
evals(i) = evals(j)
|
||||
evals(j) = temp_val
|
||||
|
||||
! swap corresponding eigenvectors
|
||||
temp_vec = evecs(:, i)
|
||||
evecs(:, i) = evecs(:, j)
|
||||
evecs(:, j) = temp_vec
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
|
||||
end subroutine
|
||||
|
||||
@ -170,6 +170,47 @@ subroutine complex_diagonalize_matrix(N,A,e)
|
||||
|
||||
end subroutine
|
||||
|
||||
subroutine complex_diagonalize_matrix_without_sort(N,A,e)
|
||||
|
||||
! Diagonalize a general complex matrix
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: N
|
||||
complex*16,intent(inout) :: A(N,N)
|
||||
complex*16,intent(out) :: e(N)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: lwork,info
|
||||
double precision,allocatable :: rwork(:)
|
||||
complex*16,allocatable :: work(:)
|
||||
complex*16,allocatable :: VL(:,:)
|
||||
complex*16,allocatable :: VR(:,:)
|
||||
|
||||
! Memory allocation
|
||||
allocate(work(1),rwork(2*N),VL(1,1),VR(N,N))
|
||||
lwork = -1
|
||||
call zgeev('N','V',N,A,N,e,VL,1,VR,N,work,lwork,rwork,info)
|
||||
lwork = max(1,int(real(work(1))))
|
||||
|
||||
deallocate(work)
|
||||
allocate(work(lwork))
|
||||
|
||||
call zgeev('N','V',N,A,N,e,VL,N,VR,N,work,lwork,rwork,info)
|
||||
call complex_sort_eigenvalues(N,e,VR)
|
||||
|
||||
|
||||
deallocate(work)
|
||||
A = VR
|
||||
|
||||
if(info /= 0) then
|
||||
print*,'Problem in diagonalize_matrix (zgeev)!!'
|
||||
end if
|
||||
|
||||
end subroutine
|
||||
|
||||
|
||||
subroutine svd(N,A,U,D,Vt)
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user