mirror of
https://github.com/pfloos/quack
synced 2024-12-22 20:34:46 +01:00
evGGF2 and evGGW
This commit is contained in:
Pierre-Francois Loos
parent
a60aa8d8ee
commit
4c9841c5cd
@ -11,7 +11,7 @@
|
|||||||
# phRPA phRPAx crRPA ppRPA
|
# phRPA phRPAx crRPA ppRPA
|
||||||
F F F F
|
F F F F
|
||||||
# G0F2 evGF2 qsGF2 G0F3 evGF3
|
# G0F2 evGF2 qsGF2 G0F3 evGF3
|
||||||
F F F F F
|
F T F F F
|
||||||
# G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW
|
# G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW
|
||||||
F F F F F F
|
F F F F F F
|
||||||
# G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
|
# G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
|
||||||
|
|||||||
@ -15,4 +15,4 @@
|
|||||||
# ACFDT: AC Kx XBS
|
# ACFDT: AC Kx XBS
|
||||||
F F T
|
F F T
|
||||||
# BSE: phBSE phBSE2 ppBSE dBSE dTDA
|
# BSE: phBSE phBSE2 ppBSE dBSE dTDA
|
||||||
T F F F T
|
F F F F T
|
||||||
|
|||||||
@ -82,7 +82,7 @@ subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regu
|
|||||||
|
|
||||||
! Print results
|
! Print results
|
||||||
|
|
||||||
call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,eGF,Ec)
|
call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
|
||||||
call print_G0F2(nBas,nO,eHF,SigC,eGF,Z,ENuc,ERHF,Ec)
|
call print_G0F2(nBas,nO,eHF,SigC,eGF,Z,ENuc,ERHF,Ec)
|
||||||
|
|
||||||
! Perform BSE2 calculation
|
! Perform BSE2 calculation
|
||||||
@ -96,7 +96,7 @@ subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regu
|
|||||||
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (singlet) =',EcBSE(1)
|
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (singlet) =',EcBSE(1)
|
||||||
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (triplet) =',EcBSE(2)
|
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (triplet) =',EcBSE(2)
|
||||||
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy =',sum(EcBSE(:))
|
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy =',sum(EcBSE(:))
|
||||||
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 total energy =',ENuc + EHF + sum(EcBSE(:))
|
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 total energy =',ENuc + ERHF + sum(EcBSE(:))
|
||||||
write(*,*)'-------------------------------------------------------------------------------'
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
write(*,*)
|
write(*,*)
|
||||||
|
|
||||||
|
|||||||
@ -31,7 +31,7 @@ subroutine GG0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, &
|
|||||||
! Local variables
|
! Local variables
|
||||||
|
|
||||||
double precision :: Ec
|
double precision :: Ec
|
||||||
double precision :: EcBSE(nspin)
|
double precision :: EcBSE
|
||||||
double precision,allocatable :: eGFlin(:)
|
double precision,allocatable :: eGFlin(:)
|
||||||
double precision,allocatable :: eGF(:)
|
double precision,allocatable :: eGF(:)
|
||||||
double precision,allocatable :: SigC(:)
|
double precision,allocatable :: SigC(:)
|
||||||
@ -85,33 +85,29 @@ subroutine GG0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, &
|
|||||||
|
|
||||||
! Perform BSE2 calculation
|
! Perform BSE2 calculation
|
||||||
|
|
||||||
! if(dophBSE) then
|
if(dophBSE) then
|
||||||
!
|
|
||||||
! call GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
|
|
||||||
|
|
||||||
! write(*,*)
|
call GGF2_phBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
|
||||||
! write(*,*)'-------------------------------------------------------------------------------'
|
|
||||||
! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (singlet) =',EcBSE(1)
|
|
||||||
! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (triplet) =',EcBSE(2)
|
|
||||||
! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy =',sum(EcBSE(:))
|
|
||||||
! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 total energy =',ENuc + EHF + sum(EcBSE(:))
|
|
||||||
! write(*,*)'-------------------------------------------------------------------------------'
|
|
||||||
! write(*,*)
|
|
||||||
|
|
||||||
! end if
|
write(*,*)
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@GG0F2 correlation energy =',EcBSE,' au'
|
||||||
|
write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@GG0F2 total energy =',ENuc + EHF + EcBSE,' au'
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
! Perform ppBSE2 calculation
|
! Perform ppBSE2 calculation
|
||||||
|
|
||||||
! if(doppBSE) then
|
! if(doppBSE) then
|
||||||
!
|
!
|
||||||
! call GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
|
! call GGF2_ppBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
|
||||||
|
|
||||||
! write(*,*)
|
! write(*,*)
|
||||||
! write(*,*)'-------------------------------------------------------------------------------'
|
! write(*,*)'-------------------------------------------------------------------------------'
|
||||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy (singlet) =',EcBSE(1),' au'
|
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0F2 correlation energy =',EcBSE,' au'
|
||||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy (triplet) =',3d0*EcBSE(2),' au'
|
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0F2 total energy =',ENuc + ERHF + EcBSE,' au'
|
||||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy =',EcBSE(1) + 3d0*EcBSE(2),' au'
|
|
||||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 total energy =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2),' au'
|
|
||||||
! write(*,*)'-------------------------------------------------------------------------------'
|
! write(*,*)'-------------------------------------------------------------------------------'
|
||||||
! write(*,*)
|
! write(*,*)
|
||||||
|
|
||||||
|
|||||||
@ -78,9 +78,8 @@ subroutine GGF(doG0F2,doevGF2,doqsGF2,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA
|
|||||||
if(doevGF2) then
|
if(doevGF2) then
|
||||||
|
|
||||||
call wall_time(start_GF)
|
call wall_time(start_GF)
|
||||||
! call evGGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, &
|
call evGGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,linearize,eta,regularize, &
|
||||||
! linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF, &
|
nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
|
||||||
! ERI,dipole_int,epsHF)
|
|
||||||
call wall_time(end_GF)
|
call wall_time(end_GF)
|
||||||
|
|
||||||
t_GF = end_GF - start_GF
|
t_GF = end_GF - start_GF
|
||||||
|
|||||||
71
src/GF/GGF2_phBSE2.f90
Normal file
71
src/GF/GGF2_phBSE2.f90
Normal file
@ -0,0 +1,71 @@
|
|||||||
|
subroutine GGF2_phBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
|
||||||
|
|
||||||
|
! Compute the second-order Bethe-Salpeter excitation energies
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
include 'parameters.h'
|
||||||
|
|
||||||
|
! Input variables
|
||||||
|
|
||||||
|
logical,intent(in) :: TDA
|
||||||
|
logical,intent(in) :: dBSE
|
||||||
|
logical,intent(in) :: dTDA
|
||||||
|
|
||||||
|
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) :: nS
|
||||||
|
double precision,intent(in) :: eGF(nBas)
|
||||||
|
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||||
|
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
|
||||||
|
logical :: dRPA = .false.
|
||||||
|
integer :: ispin
|
||||||
|
double precision,allocatable :: OmBSE(:)
|
||||||
|
double precision,allocatable :: XpY(:,:)
|
||||||
|
double precision,allocatable :: XmY(:,:)
|
||||||
|
double precision,allocatable :: A_sta(:,:)
|
||||||
|
double precision,allocatable :: B_sta(:,:)
|
||||||
|
double precision,allocatable :: KA_sta(:,:)
|
||||||
|
double precision,allocatable :: KB_sta(:,:)
|
||||||
|
|
||||||
|
! Output variables
|
||||||
|
|
||||||
|
double precision,intent(out) :: EcBSE
|
||||||
|
|
||||||
|
! Memory allocation
|
||||||
|
|
||||||
|
allocate(OmBSE(nS),XpY(nS,nS),XmY(nS,nS),A_sta(nS,nS),KA_sta(nS,nS))
|
||||||
|
allocate(B_sta(nS,nS),KB_sta(nS,nS))
|
||||||
|
|
||||||
|
ispin = 3
|
||||||
|
EcBSE = 0d0
|
||||||
|
|
||||||
|
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,A_sta)
|
||||||
|
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta)
|
||||||
|
|
||||||
|
! Compute static kernel
|
||||||
|
|
||||||
|
call GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KA_sta)
|
||||||
|
if(.not.TDA) call GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_sta)
|
||||||
|
|
||||||
|
A_sta(:,:) = A_sta(:,:) + KA_sta(:,:)
|
||||||
|
if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB_sta(:,:)
|
||||||
|
|
||||||
|
! Compute phBSE2@GF2 excitation energies
|
||||||
|
|
||||||
|
call phLR(TDA,nS,A_sta,B_sta,EcBSE,OmBSE,XpY,XmY)
|
||||||
|
call print_excitation_energies('phBSE2@GGF2',ispin,nS,OmBSE)
|
||||||
|
call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY)
|
||||||
|
|
||||||
|
! Compute dynamic correction for BSE via perturbation theory
|
||||||
|
|
||||||
|
! if(dBSE) &
|
||||||
|
! call GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY)
|
||||||
|
|
||||||
|
end subroutine
|
||||||
162
src/GF/GGF2_phBSE2_static_kernel_A.f90
Normal file
162
src/GF/GGF2_phBSE2_static_kernel_A.f90
Normal file
@ -0,0 +1,162 @@
|
|||||||
|
subroutine GGF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,KA_sta)
|
||||||
|
|
||||||
|
! Compute the resonant part of the static BSE2 matrix
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
include 'parameters.h'
|
||||||
|
|
||||||
|
! Input variables
|
||||||
|
|
||||||
|
integer,intent(in) :: ispin
|
||||||
|
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
|
||||||
|
double precision,intent(in) :: eta
|
||||||
|
double precision,intent(in) :: lambda
|
||||||
|
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||||
|
double precision,intent(in) :: eGF(nBas)
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
|
||||||
|
double precision :: dem,num
|
||||||
|
integer :: i,j,k,l
|
||||||
|
integer :: a,b,c,d
|
||||||
|
integer :: ia,jb
|
||||||
|
|
||||||
|
! Output variables
|
||||||
|
|
||||||
|
double precision,intent(out) :: KA_sta(nS,nS)
|
||||||
|
|
||||||
|
! Initialization
|
||||||
|
|
||||||
|
KA_sta(:,:) = 0d0
|
||||||
|
|
||||||
|
! Second-order correlation kernel for the block A of the singlet manifold
|
||||||
|
|
||||||
|
if(ispin == 1) then
|
||||||
|
|
||||||
|
jb = 0
|
||||||
|
!$omp parallel do default(private) shared(KA_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR)
|
||||||
|
! do j=nC+1,nO
|
||||||
|
! do b=nO+1,nBas-nR
|
||||||
|
! jb = (b-nO) + (j-1)*(nBas-nO)
|
||||||
|
!
|
||||||
|
! ia = 0
|
||||||
|
! do i=nC+1,nO
|
||||||
|
! do a=nO+1,nBas-nR
|
||||||
|
! ia = (a-nO) + (i-1)*(nBas-nO)
|
||||||
|
!
|
||||||
|
!
|
||||||
|
! do k=nC+1,nO
|
||||||
|
! do c=nO+1,nBas-nR
|
||||||
|
!
|
||||||
|
! dem = - (eGF(c) - eGF(k))
|
||||||
|
! num = 2d0*ERI(j,k,i,c)*ERI(a,c,b,k) - ERI(j,k,i,c)*ERI(a,c,k,b) &
|
||||||
|
! - ERI(j,k,c,i)*ERI(a,c,b,k) + 2d0*ERI(j,k,c,i)*ERI(a,c,k,b)
|
||||||
|
|
||||||
|
! KA_sta(ia,jb) = KA_sta(ia,jb) - num*dem/(dem**2 + eta**2)
|
||||||
|
!
|
||||||
|
! dem = + (eGF(c) - eGF(k))
|
||||||
|
! num = 2d0*ERI(j,c,i,k)*ERI(a,k,b,c) - ERI(j,c,i,k)*ERI(a,k,c,b) &
|
||||||
|
! - ERI(j,c,k,i)*ERI(a,k,b,c) + 2d0*ERI(j,c,k,i)*ERI(a,k,c,b)
|
||||||
|
|
||||||
|
! KA_sta(ia,jb) = KA_sta(ia,jb) + num*dem/(dem**2 + eta**2)
|
||||||
|
!
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
|
||||||
|
! do c=nO+1,nBas-nR
|
||||||
|
! do d=nO+1,nBas-nR
|
||||||
|
!
|
||||||
|
! dem = - (eGF(c) + eGF(d))
|
||||||
|
! num = 2d0*ERI(a,j,c,d)*ERI(c,d,i,b) - ERI(a,j,c,d)*ERI(c,d,b,i) &
|
||||||
|
! - ERI(a,j,d,c)*ERI(c,d,i,b) + 2d0*ERI(a,j,d,c)*ERI(c,d,b,i)
|
||||||
|
|
||||||
|
! KA_sta(ia,jb) = KA_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2)
|
||||||
|
!
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
|
||||||
|
! do k=nC+1,nO
|
||||||
|
! do l=nC+1,nO
|
||||||
|
|
||||||
|
! dem = - (eGF(k) + eGF(l))
|
||||||
|
! num = 2d0*ERI(a,j,k,l)*ERI(k,l,i,b) - ERI(a,j,k,l)*ERI(k,l,b,i) &
|
||||||
|
! - ERI(a,j,l,k)*ERI(k,l,i,b) + 2d0*ERI(a,j,l,k)*ERI(k,l,b,i)
|
||||||
|
|
||||||
|
! KA_sta(ia,jb) = KA_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2)
|
||||||
|
!
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
!
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
!$omp end parallel do
|
||||||
|
|
||||||
|
! end if
|
||||||
|
|
||||||
|
! Second-order correlation kernel for the block A of the triplet manifold
|
||||||
|
|
||||||
|
! if(ispin == 2) then
|
||||||
|
|
||||||
|
! jb = 0
|
||||||
|
!$omp parallel do default(private) shared(KA_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR)
|
||||||
|
! do j=nC+1,nO
|
||||||
|
! do b=nO+1,nBas-nR
|
||||||
|
! jb = (b-nO) + (j-1)*(nBas-nO)
|
||||||
|
!
|
||||||
|
! ia = 0
|
||||||
|
! do i=nC+1,nO
|
||||||
|
! do a=nO+1,nBas-nR
|
||||||
|
! ia = (a-nO) + (i-1)*(nBas-nO)
|
||||||
|
!
|
||||||
|
! do k=nC+1,nO
|
||||||
|
! do c=nO+1,nBas-nR
|
||||||
|
!
|
||||||
|
! dem = - (eGF(c) - eGF(k))
|
||||||
|
! num = 2d0*ERI(j,k,i,c)*ERI(a,c,b,k) - ERI(j,k,i,c)*ERI(a,c,k,b) - ERI(j,k,c,i)*ERI(a,c,b,k)
|
||||||
|
|
||||||
|
! KA_sta(ia,jb) = KA_sta(ia,jb) - num*dem/(dem**2 + eta**2)
|
||||||
|
!
|
||||||
|
! dem = + (eGF(c) - eGF(k))
|
||||||
|
! num = 2d0*ERI(j,c,i,k)*ERI(a,k,b,c) - ERI(j,c,i,k)*ERI(a,k,c,b) - ERI(j,c,k,i)*ERI(a,k,b,c)
|
||||||
|
|
||||||
|
! KA_sta(ia,jb) = KA_sta(ia,jb) + num*dem/(dem**2 + eta**2)
|
||||||
|
!
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
|
||||||
|
! do c=nO+1,nBas-nR
|
||||||
|
! do d=nO+1,nBas-nR
|
||||||
|
!
|
||||||
|
! dem = - (eGF(c) + eGF(d))
|
||||||
|
! num = ERI(a,j,c,d)*ERI(c,d,b,i) + ERI(a,j,d,c)*ERI(c,d,i,b)
|
||||||
|
|
||||||
|
! KA_sta(ia,jb) = KA_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2)
|
||||||
|
!
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
|
||||||
|
! do k=nC+1,nO
|
||||||
|
! do l=nC+1,nO
|
||||||
|
|
||||||
|
! dem = - (eGF(k) + eGF(l))
|
||||||
|
! num = ERI(a,j,k,l)*ERI(k,l,b,i) + ERI(a,j,l,k)*ERI(k,l,i,b)
|
||||||
|
|
||||||
|
! KA_sta(ia,jb) = KA_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2)
|
||||||
|
!
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
!
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
!$omp end parallel do
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
|
||||||
|
end subroutine
|
||||||
161
src/GF/GGF2_phBSE2_static_kernel_B.f90
Normal file
161
src/GF/GGF2_phBSE2_static_kernel_B.f90
Normal file
@ -0,0 +1,161 @@
|
|||||||
|
subroutine GGF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,KB_sta)
|
||||||
|
|
||||||
|
! Compute the anti-resonant part of the static BSE2 matrix
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
include 'parameters.h'
|
||||||
|
|
||||||
|
! Input variables
|
||||||
|
|
||||||
|
integer,intent(in) :: ispin
|
||||||
|
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
|
||||||
|
double precision,intent(in) :: eta
|
||||||
|
double precision,intent(in) :: lambda
|
||||||
|
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||||
|
double precision,intent(in) :: eGF(nBas)
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
|
||||||
|
double precision :: dem,num
|
||||||
|
integer :: i,j,k,l
|
||||||
|
integer :: a,b,c,d
|
||||||
|
integer :: ia,jb
|
||||||
|
|
||||||
|
! Output variables
|
||||||
|
|
||||||
|
double precision,intent(out) :: KB_sta(nS,nS)
|
||||||
|
|
||||||
|
! Initialization
|
||||||
|
|
||||||
|
KB_sta(:,:) = 0d0
|
||||||
|
|
||||||
|
! Second-order correlation kernel for the block A of the singlet manifold
|
||||||
|
|
||||||
|
if(ispin == 1) then
|
||||||
|
|
||||||
|
jb = 0
|
||||||
|
!$omp parallel do default(private) shared(KB_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR)
|
||||||
|
! do j=nC+1,nO
|
||||||
|
! do b=nO+1,nBas-nR
|
||||||
|
! jb = (b-nO) + (j-1)*(nBas-nO)
|
||||||
|
|
||||||
|
! ia = 0
|
||||||
|
! do i=nC+1,nO
|
||||||
|
! do a=nO+1,nBas-nR
|
||||||
|
! ia = (a-nO) + (i-1)*(nBas-nO)
|
||||||
|
!
|
||||||
|
! do k=nC+1,nO
|
||||||
|
! do c=nO+1,nBas-nR
|
||||||
|
!
|
||||||
|
! dem = + eGF(k) - eGF(c)
|
||||||
|
! num = 2d0*ERI(b,k,i,c)*ERI(a,c,j,k) - ERI(b,k,i,c)*ERI(a,c,k,j) &
|
||||||
|
! - ERI(b,k,c,i)*ERI(a,c,j,k) + 2d0*ERI(b,k,c,i)*ERI(a,c,k,j)
|
||||||
|
|
||||||
|
! KB_sta(ia,jb) = KB_sta(ia,jb) - num*dem/(dem**2 + eta**2)
|
||||||
|
!
|
||||||
|
! dem = - eGF(c) + eGF(k)
|
||||||
|
! num = 2d0*ERI(b,c,i,k)*ERI(a,k,j,c) - ERI(b,c,i,k)*ERI(a,k,c,j) &
|
||||||
|
! - ERI(b,c,k,i)*ERI(a,k,j,c) + 2d0*ERI(b,c,k,i)*ERI(a,k,c,j)
|
||||||
|
|
||||||
|
! KB_sta(ia,jb) = KB_sta(ia,jb) - num*dem/(dem**2 + eta**2)
|
||||||
|
!
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
|
||||||
|
! do c=nO+1,nBas-nR
|
||||||
|
! do d=nO+1,nBas-nR
|
||||||
|
!
|
||||||
|
! dem = - eGF(c) - eGF(d)
|
||||||
|
! num = 2d0*ERI(a,b,c,d)*ERI(c,d,i,j) - ERI(a,b,c,d)*ERI(c,d,j,i) &
|
||||||
|
! - ERI(a,b,d,c)*ERI(c,d,i,j) + 2d0*ERI(a,b,d,c)*ERI(c,d,j,i)
|
||||||
|
|
||||||
|
! KB_sta(ia,jb) = KB_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2)
|
||||||
|
!
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
|
||||||
|
! do k=nC+1,nO
|
||||||
|
! do l=nC+1,nO
|
||||||
|
|
||||||
|
! dem = + eGF(k) + eGF(l)
|
||||||
|
! num = 2d0*ERI(a,b,k,l)*ERI(k,l,i,j) - ERI(a,b,k,l)*ERI(k,l,j,i) &
|
||||||
|
! - ERI(a,b,l,k)*ERI(k,l,i,j) + 2d0*ERI(a,b,l,k)*ERI(k,l,j,i)
|
||||||
|
|
||||||
|
! KB_sta(ia,jb) = KB_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2)
|
||||||
|
!
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
!
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
!$omp end parallel do
|
||||||
|
|
||||||
|
! end if
|
||||||
|
|
||||||
|
! Second-order correlation kernel for the block A of the triplet manifold
|
||||||
|
|
||||||
|
! if(ispin == 2) then
|
||||||
|
|
||||||
|
! jb = 0
|
||||||
|
!$omp parallel do default(private) shared(KB_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR)
|
||||||
|
! do j=nC+1,nO
|
||||||
|
! do b=nO+1,nBas-nR
|
||||||
|
! jb = (b-nO) + (j-1)*(nBas-nO)
|
||||||
|
!
|
||||||
|
! ia = 0
|
||||||
|
! do i=nC+1,nO
|
||||||
|
! do a=nO+1,nBas-nR
|
||||||
|
! ia = (a-nO) + (i-1)*(nBas-nO)
|
||||||
|
!
|
||||||
|
! do k=nC+1,nO
|
||||||
|
! do c=nO+1,nBas-nR
|
||||||
|
!
|
||||||
|
! dem = + eGF(k) - eGF(c)
|
||||||
|
! num = 2d0*ERI(b,k,i,c)*ERI(a,c,j,k) - ERI(b,k,i,c)*ERI(a,c,k,j) - ERI(b,k,c,i)*ERI(a,c,j,k)
|
||||||
|
|
||||||
|
! KB_sta(ia,jb) = KB_sta(ia,jb) - num*dem/(dem**2 + eta**2)
|
||||||
|
!
|
||||||
|
! dem = - eGF(c) + eGF(k)
|
||||||
|
! num = 2d0*ERI(b,c,i,k)*ERI(a,k,j,c) - ERI(b,c,i,k)*ERI(a,k,c,j) - ERI(b,c,k,i)*ERI(a,k,j,c)
|
||||||
|
|
||||||
|
! KB_sta(ia,jb) = KB_sta(ia,jb) - num*dem/(dem**2 + eta**2)
|
||||||
|
!
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
|
||||||
|
! do c=nO+1,nBas-nR
|
||||||
|
! do d=nO+1,nBas-nR
|
||||||
|
!
|
||||||
|
! dem = - eGF(c) - eGF(d)
|
||||||
|
! num = ERI(a,b,c,d)*ERI(c,d,j,i) + ERI(a,b,d,c)*ERI(c,d,i,j)
|
||||||
|
|
||||||
|
! KB_sta(ia,jb) = KB_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2)
|
||||||
|
!
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
|
||||||
|
! do k=nC+1,nO
|
||||||
|
! do l=nC+1,nO
|
||||||
|
|
||||||
|
! dem = + eGF(k) + eGF(l)
|
||||||
|
! num = ERI(a,b,k,l)*ERI(k,l,j,i) + ERI(a,b,l,k)*ERI(k,l,i,j)
|
||||||
|
|
||||||
|
! KB_sta(ia,jb) = KB_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2)
|
||||||
|
!
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
!
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
!$omp end parallel do
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
|
||||||
|
end subroutine
|
||||||
@ -107,7 +107,7 @@ subroutine evGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,tr
|
|||||||
|
|
||||||
! Print results
|
! Print results
|
||||||
|
|
||||||
call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,eGF,Ec)
|
call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
|
||||||
call print_evGF2(nBas,nO,nSCF,Conv,eHF,SigC,Z,eGF,ENuc,ERHF,Ec)
|
call print_evGF2(nBas,nO,nSCF,Conv,eHF,SigC,Z,eGF,ENuc,ERHF,Ec)
|
||||||
|
|
||||||
! DIIS extrapolation
|
! DIIS extrapolation
|
||||||
@ -153,7 +153,7 @@ subroutine evGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,tr
|
|||||||
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGF2 correlation energy (singlet) =',EcBSE(1)
|
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGF2 correlation energy (singlet) =',EcBSE(1)
|
||||||
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGF2 correlation energy (triplet) =',EcBSE(2)
|
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGF2 correlation energy (triplet) =',EcBSE(2)
|
||||||
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGF2 correlation energy =',sum(EcBSE(:))
|
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGF2 correlation energy =',sum(EcBSE(:))
|
||||||
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGF2 total energy =',ENuc + EHF + sum(EcBSE(:))
|
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGF2 total energy =',ENuc + ERHF + sum(EcBSE(:))
|
||||||
write(*,*)'-------------------------------------------------------------------------------'
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
write(*,*)
|
write(*,*)
|
||||||
|
|
||||||
|
|||||||
173
src/GF/evGGF2.f90
Normal file
173
src/GF/evGGF2.f90
Normal file
@ -0,0 +1,173 @@
|
|||||||
|
subroutine evGGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, &
|
||||||
|
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
|
||||||
|
|
||||||
|
! Perform eigenvalue self-consistent second-order Green function calculation
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
include 'parameters.h'
|
||||||
|
|
||||||
|
! Input variables
|
||||||
|
|
||||||
|
logical,intent(in) :: dophBSE
|
||||||
|
logical,intent(in) :: doppBSE
|
||||||
|
logical,intent(in) :: TDA
|
||||||
|
logical,intent(in) :: dBSE
|
||||||
|
logical,intent(in) :: dTDA
|
||||||
|
integer,intent(in) :: maxSCF
|
||||||
|
double precision,intent(in) :: thresh
|
||||||
|
integer,intent(in) :: max_diis
|
||||||
|
logical,intent(in) :: linearize
|
||||||
|
double precision,intent(in) :: eta
|
||||||
|
logical,intent(in) :: regularize
|
||||||
|
|
||||||
|
integer,intent(in) :: nBas
|
||||||
|
integer,intent(in) :: nO
|
||||||
|
integer,intent(in) :: nC
|
||||||
|
integer,intent(in) :: nV
|
||||||
|
integer,intent(in) :: nR
|
||||||
|
integer,intent(in) :: nS
|
||||||
|
double precision,intent(in) :: ENuc
|
||||||
|
double precision,intent(in) :: ERHF
|
||||||
|
double precision,intent(in) :: eHF(nBas)
|
||||||
|
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||||
|
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
|
||||||
|
integer :: nSCF
|
||||||
|
integer :: n_diis
|
||||||
|
double precision :: Ec
|
||||||
|
double precision :: EcBSE
|
||||||
|
double precision :: Conv
|
||||||
|
double precision :: rcond
|
||||||
|
double precision,allocatable :: eGF(:)
|
||||||
|
double precision,allocatable :: eOld(:)
|
||||||
|
double precision,allocatable :: SigC(:)
|
||||||
|
double precision,allocatable :: Z(:)
|
||||||
|
double precision,allocatable :: error_diis(:,:)
|
||||||
|
double precision,allocatable :: e_diis(:,:)
|
||||||
|
|
||||||
|
! Hello world
|
||||||
|
|
||||||
|
write(*,*)
|
||||||
|
write(*,*)'************************************************'
|
||||||
|
write(*,*)'| Second-order Green function calculation |'
|
||||||
|
write(*,*)'************************************************'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
! Memory allocation
|
||||||
|
|
||||||
|
allocate(SigC(nBas),Z(nBas),eGF(nBas),eOld(nBas),error_diis(nBas,max_diis),e_diis(nBas,max_diis))
|
||||||
|
|
||||||
|
! Initialization
|
||||||
|
|
||||||
|
Conv = 1d0
|
||||||
|
nSCF = 0
|
||||||
|
n_diis = 0
|
||||||
|
e_diis(:,:) = 0d0
|
||||||
|
error_diis(:,:) = 0d0
|
||||||
|
eGF(:) = eHF(:)
|
||||||
|
eOld(:) = eHF(:)
|
||||||
|
rcond = 0d0
|
||||||
|
|
||||||
|
!------------------------------------------------------------------------
|
||||||
|
! Main SCF loop
|
||||||
|
!------------------------------------------------------------------------
|
||||||
|
|
||||||
|
do while(Conv > thresh .and. nSCF < maxSCF)
|
||||||
|
|
||||||
|
! Frequency-dependent second-order contribution
|
||||||
|
|
||||||
|
if(regularize) then
|
||||||
|
|
||||||
|
! call GGF2_reg_self_energy_diag(eta,nBas,nC,nO,nV,nR,eGF,ERI,SigC,Z)
|
||||||
|
|
||||||
|
else
|
||||||
|
|
||||||
|
call GGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,eGF,ERI,SigC,Z)
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
! Solve the quasi-particle equation
|
||||||
|
|
||||||
|
if(linearize) then
|
||||||
|
|
||||||
|
else
|
||||||
|
|
||||||
|
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
call GGF2_QP_graph(eta,nBas,nC,nO,nV,nR,eHF,ERI,eOld,eOld,eGF,Z)
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
Conv = maxval(abs(eGF - eOld))
|
||||||
|
|
||||||
|
! Print results
|
||||||
|
|
||||||
|
call GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
|
||||||
|
call print_evGF2(nBas,nO,nSCF,Conv,eHF,SigC,Z,eGF,ENuc,ERHF,Ec)
|
||||||
|
|
||||||
|
! DIIS extrapolation
|
||||||
|
|
||||||
|
n_diis = min(n_diis+1,max_diis)
|
||||||
|
call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGF-eOld,eGF)
|
||||||
|
|
||||||
|
if(abs(rcond) < 1d-15) n_diis = 0
|
||||||
|
|
||||||
|
eOld(:) = eGF(:)
|
||||||
|
|
||||||
|
! Increment
|
||||||
|
|
||||||
|
nSCF = nSCF + 1
|
||||||
|
|
||||||
|
end do
|
||||||
|
!------------------------------------------------------------------------
|
||||||
|
! End main SCF loop
|
||||||
|
!------------------------------------------------------------------------
|
||||||
|
|
||||||
|
! Did it actually converge?
|
||||||
|
|
||||||
|
if(nSCF == maxSCF+1) then
|
||||||
|
|
||||||
|
write(*,*)
|
||||||
|
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||||
|
write(*,*)' Convergence failed '
|
||||||
|
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
stop
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
! Perform BSE2 calculation
|
||||||
|
|
||||||
|
if(dophBSE) then
|
||||||
|
|
||||||
|
call GGF2_phBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
|
||||||
|
|
||||||
|
write(*,*)
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@evGGF2 correlation energy =',EcBSE,' au'
|
||||||
|
write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@evGGF2 total energy =',ENuc + ERHF + EcBSE,' au'
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
! Perform ppBSE2 calculation
|
||||||
|
|
||||||
|
! if(doppBSE) then
|
||||||
|
|
||||||
|
! call GGF2_ppBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
|
||||||
|
|
||||||
|
! write(*,*)
|
||||||
|
! write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@evGGF2 correlation energy =',EcBSE),' au'
|
||||||
|
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@evGGF2 total energy =',ENuc + ERHF + EcBSE,' au'
|
||||||
|
! write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
! write(*,*)
|
||||||
|
|
||||||
|
! end if
|
||||||
|
|
||||||
|
end subroutine
|
||||||
@ -202,7 +202,7 @@ subroutine GG0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,d
|
|||||||
|
|
||||||
! end if
|
! end if
|
||||||
|
|
||||||
! end if
|
end if
|
||||||
|
|
||||||
! if(doppBSE) then
|
! if(doppBSE) then
|
||||||
|
|
||||||
@ -217,6 +217,6 @@ subroutine GG0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,d
|
|||||||
! write(*,*)'-------------------------------------------------------------------------------'
|
! write(*,*)'-------------------------------------------------------------------------------'
|
||||||
! write(*,*)
|
! write(*,*)
|
||||||
|
|
||||||
end if
|
! end if
|
||||||
|
|
||||||
end subroutine
|
end subroutine
|
||||||
|
|||||||
@ -88,8 +88,8 @@ subroutine GGW(doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT, &
|
|||||||
if(doevGW) then
|
if(doevGW) then
|
||||||
|
|
||||||
call wall_time(start_GW)
|
call wall_time(start_GW)
|
||||||
! call evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
|
call evGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
|
||||||
! linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
|
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
|
||||||
call wall_time(end_GW)
|
call wall_time(end_GW)
|
||||||
|
|
||||||
t_GW = end_GW - start_GW
|
t_GW = end_GW - start_GW
|
||||||
|
|||||||
267
src/GW/evGGW.f90
Normal file
267
src/GW/evGGW.f90
Normal file
@ -0,0 +1,267 @@
|
|||||||
|
subroutine evGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
|
||||||
|
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
|
||||||
|
|
||||||
|
! Perform self-consistent eigenvalue-only GW calculation
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
include 'parameters.h'
|
||||||
|
|
||||||
|
! Input variables
|
||||||
|
|
||||||
|
integer,intent(in) :: maxSCF
|
||||||
|
integer,intent(in) :: max_diis
|
||||||
|
double precision,intent(in) :: thresh
|
||||||
|
double precision,intent(in) :: ENuc
|
||||||
|
double precision,intent(in) :: ERHF
|
||||||
|
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) :: linearize
|
||||||
|
double precision,intent(in) :: eta
|
||||||
|
logical,intent(in) :: regularize
|
||||||
|
|
||||||
|
integer,intent(in) :: nBas
|
||||||
|
integer,intent(in) :: nC
|
||||||
|
integer,intent(in) :: nO
|
||||||
|
integer,intent(in) :: nV
|
||||||
|
integer,intent(in) :: nR
|
||||||
|
integer,intent(in) :: nS
|
||||||
|
double precision,intent(in) :: eHF(nBas)
|
||||||
|
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||||
|
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
|
||||||
|
logical :: linear_mixing
|
||||||
|
logical :: dRPA = .true.
|
||||||
|
integer :: ispin
|
||||||
|
integer :: nSCF
|
||||||
|
integer :: n_diis
|
||||||
|
double precision :: rcond
|
||||||
|
double precision :: Conv
|
||||||
|
double precision :: EcRPA
|
||||||
|
double precision :: EcBSE
|
||||||
|
double precision :: EcGM
|
||||||
|
double precision :: alpha
|
||||||
|
double precision,allocatable :: Aph(:,:)
|
||||||
|
double precision,allocatable :: Bph(:,:)
|
||||||
|
double precision,allocatable :: error_diis(:,:)
|
||||||
|
double precision,allocatable :: e_diis(:,:)
|
||||||
|
double precision,allocatable :: eGW(:)
|
||||||
|
double precision,allocatable :: eOld(:)
|
||||||
|
double precision,allocatable :: Z(:)
|
||||||
|
double precision,allocatable :: SigC(:)
|
||||||
|
double precision,allocatable :: Om(:)
|
||||||
|
double precision,allocatable :: XpY(:,:)
|
||||||
|
double precision,allocatable :: XmY(:,:)
|
||||||
|
double precision,allocatable :: rho(:,:,:)
|
||||||
|
|
||||||
|
! Hello world
|
||||||
|
|
||||||
|
write(*,*)
|
||||||
|
write(*,*)'************************************************'
|
||||||
|
write(*,*)'| Self-consistent evGW calculation |'
|
||||||
|
write(*,*)'************************************************'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
! TDA for W
|
||||||
|
|
||||||
|
if(TDA_W) then
|
||||||
|
write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
|
||||||
|
write(*,*)
|
||||||
|
end if
|
||||||
|
|
||||||
|
! TDA
|
||||||
|
|
||||||
|
if(TDA) then
|
||||||
|
write(*,*) 'Tamm-Dancoff approximation activated!'
|
||||||
|
write(*,*)
|
||||||
|
end if
|
||||||
|
|
||||||
|
! Linear mixing
|
||||||
|
|
||||||
|
linear_mixing = .false.
|
||||||
|
alpha = 0.2d0
|
||||||
|
|
||||||
|
! Memory allocation
|
||||||
|
|
||||||
|
allocate(Aph(nS,nS),Bph(nS,nS),eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas), &
|
||||||
|
Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis))
|
||||||
|
|
||||||
|
! Initialization
|
||||||
|
|
||||||
|
nSCF = 0
|
||||||
|
ispin = 3
|
||||||
|
n_diis = 0
|
||||||
|
Conv = 1d0
|
||||||
|
e_diis(:,:) = 0d0
|
||||||
|
error_diis(:,:) = 0d0
|
||||||
|
eGW(:) = eHF(:)
|
||||||
|
eOld(:) = eGW(:)
|
||||||
|
Z(:) = 1d0
|
||||||
|
rcond = 0d0
|
||||||
|
|
||||||
|
!------------------------------------------------------------------------
|
||||||
|
! Main loop
|
||||||
|
!------------------------------------------------------------------------
|
||||||
|
|
||||||
|
do while(Conv > thresh .and. nSCF <= maxSCF)
|
||||||
|
|
||||||
|
! Compute screening
|
||||||
|
|
||||||
|
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
|
||||||
|
if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
|
||||||
|
|
||||||
|
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
|
||||||
|
|
||||||
|
! Compute spectral weights
|
||||||
|
|
||||||
|
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho)
|
||||||
|
|
||||||
|
! Compute correlation part of the self-energy
|
||||||
|
|
||||||
|
if(regularize) call GW_regularization(nBas,nC,nO,nV,nR,nS,eGW,Om,rho)
|
||||||
|
|
||||||
|
call GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
|
||||||
|
|
||||||
|
! Solve the quasi-particle equation
|
||||||
|
|
||||||
|
if(linearize) then
|
||||||
|
|
||||||
|
write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
eGW(:) = eHF(:) + SigC(:)
|
||||||
|
|
||||||
|
else
|
||||||
|
|
||||||
|
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
call GGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eOld,eOld,eGW,Z)
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
! Convergence criteria
|
||||||
|
|
||||||
|
Conv = maxval(abs(eGW - eOld))
|
||||||
|
|
||||||
|
! Print results
|
||||||
|
|
||||||
|
call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM)
|
||||||
|
|
||||||
|
! Linear mixing or DIIS extrapolation
|
||||||
|
|
||||||
|
if(linear_mixing) then
|
||||||
|
|
||||||
|
eGW(:) = alpha*eGW(:) + (1d0 - alpha)*eOld(:)
|
||||||
|
|
||||||
|
else
|
||||||
|
|
||||||
|
n_diis = min(n_diis+1,max_diis)
|
||||||
|
if(abs(rcond) > 1d-7) then
|
||||||
|
call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGW-eOld,eGW)
|
||||||
|
else
|
||||||
|
n_diis = 0
|
||||||
|
end if
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
! Save quasiparticles energy for next cycle
|
||||||
|
|
||||||
|
eOld(:) = eGW(:)
|
||||||
|
|
||||||
|
! Increment
|
||||||
|
|
||||||
|
nSCF = nSCF + 1
|
||||||
|
|
||||||
|
end do
|
||||||
|
!------------------------------------------------------------------------
|
||||||
|
! End main loop
|
||||||
|
!------------------------------------------------------------------------
|
||||||
|
|
||||||
|
! Did it actually converge?
|
||||||
|
|
||||||
|
if(nSCF == maxSCF+1) then
|
||||||
|
|
||||||
|
write(*,*)
|
||||||
|
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||||
|
write(*,*)' Convergence failed '
|
||||||
|
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
stop
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
! Deallocate memory
|
||||||
|
|
||||||
|
deallocate(Aph,Bph,eOld,Z,SigC,Om,XpY,XmY,rho,error_diis,e_diis)
|
||||||
|
|
||||||
|
! Perform BSE calculation
|
||||||
|
|
||||||
|
if(dophBSE) then
|
||||||
|
|
||||||
|
call GGW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGW,eGW,EcBSE)
|
||||||
|
|
||||||
|
write(*,*)
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy =',EcBSE
|
||||||
|
write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW total energy =',ENuc + ERHF + EcBSE
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
! Compute the BSE correlation energy via the adiabatic connection
|
||||||
|
|
||||||
|
! 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 GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcBSE)
|
||||||
|
|
||||||
|
! write(*,*)
|
||||||
|
! write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (singlet) =',EcBSE(1)
|
||||||
|
! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (triplet) =',EcBSE(2)
|
||||||
|
! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy =',EcBSE(1) + EcBSE(2)
|
||||||
|
! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
|
||||||
|
! write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
! write(*,*)
|
||||||
|
|
||||||
|
! end if
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
! if(doppBSE) then
|
||||||
|
|
||||||
|
! call GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
|
||||||
|
|
||||||
|
! write(*,*)
|
||||||
|
! write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
! write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy (singlet) =',EcBSE(1)
|
||||||
|
! write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy (triplet) =',3d0*EcBSE(2)
|
||||||
|
! write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy =',EcBSE(1) + 3d0*EcBSE(2)
|
||||||
|
! write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2)
|
||||||
|
! write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
! write(*,*)
|
||||||
|
|
||||||
|
! end if
|
||||||
|
|
||||||
|
end subroutine
|
||||||
@ -83,8 +83,6 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole)
|
|||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
|
|
||||||
print*,Sx2,Sy2,Sz2
|
|
||||||
|
|
||||||
S2 = Sx2 + Sy2 + Sz2
|
S2 = Sx2 + Sy2 + Sz2
|
||||||
|
|
||||||
! Dump results
|
! Dump results
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user