mirror of
https://github.com/pfloos/quack
synced 2025-05-06 07:05:33 +02:00
complex qsGF2 with(out) reg
This commit is contained in:
Loris Burth
parent
cb647905cc
commit
291153e162
@ -40,7 +40,7 @@ subroutine RGF2_reg_self_energy(eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
|
||||
! Parameters for regularized calculations !
|
||||
!-----------------------------------------!
|
||||
|
||||
s = 100d0
|
||||
s = 500d0
|
||||
|
||||
!----------------------------------------------------!
|
||||
! Compute GF2 self-energy and renormalization factor !
|
||||
|
||||
@ -1,7 +1,8 @@
|
||||
subroutine complex_RGF(dotest,docG0F2,doevGF2,maxSCF, &
|
||||
subroutine complex_RGF(dotest,docG0F2,doevGF2,doqsGF2,maxSCF, &
|
||||
thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize, &
|
||||
eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF, &
|
||||
S,X,T,V,Hc,ERI_AO,ERI_MO,CAP,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
|
||||
S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF, &
|
||||
CAP_AO,CAP_MO)
|
||||
|
||||
! Perform a one-shot second-order Green function calculation
|
||||
|
||||
@ -11,7 +12,7 @@ subroutine complex_RGF(dotest,docG0F2,doevGF2,maxSCF,
|
||||
! Input variables
|
||||
|
||||
logical,intent(in) :: dotest
|
||||
logical,intent(in) :: docG0F2,doevGF2
|
||||
logical,intent(in) :: docG0F2,doevGF2,doqsGF2
|
||||
|
||||
integer,intent(in) :: maxSCF
|
||||
double precision,intent(in) :: thresh
|
||||
@ -46,7 +47,8 @@ subroutine complex_RGF(dotest,docG0F2,doevGF2,maxSCF,
|
||||
complex*16,intent(in) :: cHF(nBas,nOrb)
|
||||
complex*16,intent(in) :: PHF(nBas,nBas)
|
||||
complex*16,intent(in) :: S(nBas,nBas)
|
||||
complex*16,intent(in) :: CAP(nBas,nBas)
|
||||
complex*16,intent(in) :: CAP_AO(nBas,nBas)
|
||||
complex*16,intent(in) :: CAP_MO(nBas,nBas)
|
||||
double precision,intent(in) :: T(nBas,nBas)
|
||||
double precision,intent(in) :: V(nBas,nBas)
|
||||
double precision,intent(in) :: Hc(nBas,nBas)
|
||||
@ -70,7 +72,7 @@ subroutine complex_RGF(dotest,docG0F2,doevGF2,maxSCF,
|
||||
call wall_time(start_GF)
|
||||
call complex_cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet, &
|
||||
linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS, &
|
||||
ENuc,ERHF,ERI_MO,CAP,dipole_int_MO,eHF)
|
||||
ENuc,ERHF,ERI_MO,CAP_MO,dipole_int_MO,eHF)
|
||||
call wall_time(end_GF)
|
||||
|
||||
t_GF = end_GF - start_GF
|
||||
@ -90,4 +92,21 @@ subroutine complex_RGF(dotest,docG0F2,doevGF2,maxSCF,
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
if(doqsGF2) then
|
||||
|
||||
call wall_time(start_GF)
|
||||
call complex_qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA, &
|
||||
dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc, &
|
||||
rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, &
|
||||
ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF, &
|
||||
CAP_AO,CAP_MO)
|
||||
call wall_time(end_GF)
|
||||
|
||||
t_GF = end_GF - start_GF
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF2 = ',t_GF,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
end subroutine
|
||||
|
||||
138
src/GF/complex_cRGF2_reg_self_energy.f90
Normal file
138
src/GF/complex_cRGF2_reg_self_energy.f90
Normal file
@ -0,0 +1,138 @@
|
||||
subroutine complex_cRGF2_reg_self_energy(flow,eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
|
||||
|
||||
! Compute diagonal part of the GF2 self-energy and its renormalization factor
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
double precision,intent(in) :: eta
|
||||
double precision,intent(in) :: flow
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
complex*16,intent(in) :: e(nBas)
|
||||
complex*16,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,j,a,b
|
||||
integer :: p,q
|
||||
double precision :: eps_p
|
||||
double precision :: eps_q
|
||||
double precision :: s
|
||||
double precision :: eta_tilde_p
|
||||
double precision :: eta_tilde_q
|
||||
complex*16 :: num
|
||||
double precision,allocatable :: Re_DS(:)
|
||||
double precision,allocatable :: Im_DS(:)
|
||||
complex*16 :: z_dummy
|
||||
double precision,allocatable :: Re_SigC(:,:)
|
||||
double precision,allocatable :: Im_SigC(:,:)
|
||||
double precision,allocatable :: Re_Z(:)
|
||||
double precision,allocatable :: Im_Z(:)
|
||||
|
||||
! Output variables
|
||||
|
||||
complex*16,intent(out) :: SigC(nBas,nBas)
|
||||
complex*16,intent(out) :: Z(nBas)
|
||||
|
||||
! Initialize
|
||||
allocate(Re_DS(nBas),Im_DS(nBas),Re_SigC(nBas,nBas),Im_SigC(nBas,nBas),&
|
||||
Re_Z(nBas),Im_Z(nBas))
|
||||
Re_SigC(:,:) = 0d0
|
||||
Im_SigC(:,:) = 0d0
|
||||
Re_DS(:) = 0d0
|
||||
Im_DS(:) = 0d0
|
||||
s = flow
|
||||
|
||||
! Compute GF2 self-energy
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP SHARED(Re_DS,Im_DS,Re_SigC,Im_SigC,ERI,eta,nC,nO,nBas,nR,e,s) &
|
||||
!$OMP PRIVATE(p,q,i,j,a,eps_p,eps_q,num,eta_tilde_p,eta_tilde_q,z_dummy) &
|
||||
!$OMP DEFAULT(NONE)
|
||||
!$OMP DO
|
||||
do p=nC+1,nBas-nR
|
||||
do q=nC+1,nBas-nR
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=nO+1,nBas-nR
|
||||
|
||||
eps_p = real(e(p)) + real(e(a)) - real(e(i)) - real(e(j))
|
||||
eps_q = real(e(q)) + real(e(a)) - real(e(i)) - real(e(j))
|
||||
eta_tilde_p = eta - aimag(e(p)) + aimag(e(i)) - (aimag(e(a)) - aimag(e(j)))
|
||||
eta_tilde_q = eta - aimag(e(q)) + aimag(e(i)) - (aimag(e(a)) - aimag(e(j)))
|
||||
num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(q,a,i,j)&
|
||||
*(1d0 - exp(-s*(eps_p**2+eta_tilde_p**2 + eps_q**2 + eta_tilde_q**2)))
|
||||
|
||||
z_dummy = num*cmplx((eps_p + eps_q)/(eps_p**2 + eta_tilde_p**2 +eps_q**2 + eta_tilde_q**2 ),&
|
||||
(eta_tilde_p + eta_tilde_q)/(eps_p**2 + eta_tilde_p**2 +eps_q**2 + eta_tilde_q**2 ),kind=8)
|
||||
Re_SigC(p,q) = Re_SigC(p,q) + real(z_dummy)
|
||||
Im_SigC(p,q) = Im_SigC(p,q) + aimag(z_dummy)
|
||||
if(p==q) then
|
||||
z_dummy = num*cmplx(-(eps_p**2 - eta_tilde_p**2)/(eps_p**2 + eta_tilde_p**2)**2,&
|
||||
-2*eta_tilde_p*eps_p/(eps_p**2 + eta_tilde_p**2)**2,kind=8)
|
||||
Re_DS(p) = Re_DS(p) + real(z_dummy)
|
||||
Im_DS(p) = Im_DS(p) + aimag(z_dummy)
|
||||
end if
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP SHARED(Re_DS,Im_DS,Re_SigC,Im_SigC,ERI,eta,nC,nO,nBas,nR,e,s) &
|
||||
!$OMP PRIVATE(p,q,i,a,b,eps_p,eps_q,num,eta_tilde_p,eta_tilde_q,z_dummy) &
|
||||
!$OMP DEFAULT(NONE)
|
||||
!$OMP DO
|
||||
do p=nC+1,nBas-nR
|
||||
do q=nC+1,nBas-nR
|
||||
do i=nC+1,nO
|
||||
do a=nO+1,nBas-nR
|
||||
do b=nO+1,nBas-nR
|
||||
|
||||
eps_p = real(e(p)) + real(e(i)) - real(e(a)) - real(e(b))
|
||||
eps_q = real(e(q)) + real(e(i)) - real(e(a)) - real(e(b))
|
||||
eta_tilde_p = eta + aimag(e(p)) - aimag(e(a)) - aimag(e(b)) + aimag(e(i))
|
||||
eta_tilde_q = eta + aimag(e(q)) - aimag(e(a)) - aimag(e(b)) + aimag(e(i))
|
||||
num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(q,i,a,b)&
|
||||
*(1d0 - exp(-s*(eps_p**2+eta_tilde_p**2 + eps_q**2 + eta_tilde_q**2)))
|
||||
|
||||
|
||||
z_dummy = num*cmplx((eps_p + eps_q)/(eps_p**2+eta_tilde_p**2 + eps_q**2 + eta_tilde_q**2),&
|
||||
-(eta_tilde_p + eta_tilde_q)/(eps_p**2+eta_tilde_p**2 + eps_q**2 + eta_tilde_q**2),kind=8)
|
||||
Re_SigC(p,q) = Re_SigC(p,q) + real(z_dummy)
|
||||
Im_SigC(p,q) = Im_SigC(p,q) + aimag(z_dummy)
|
||||
if(p==q) then
|
||||
z_dummy = num*cmplx(-(eps_p**2 - eta_tilde_p**2)/(eps_p**2 + eta_tilde_p**2)**2,&
|
||||
2*eta_tilde_p*eps_p/(eps_p**2 + eta_tilde_p**2)**2,kind=8)
|
||||
Re_DS(p) = Re_DS(p) + real(z_dummy)
|
||||
Im_DS(p) = Im_DS(p) + aimag(z_dummy)
|
||||
end if
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
Re_Z(:) = (1d0-Re_DS(:))/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
|
||||
Im_Z(:) = Im_DS(:)/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
|
||||
|
||||
Z = cmplx(Re_Z,Im_Z,kind=8)
|
||||
SigC = cmplx(Re_SigC,Im_SigC,kind=8)
|
||||
|
||||
|
||||
deallocate(Re_DS,Im_DS,Re_Z,Im_Z,Re_SigC,Im_SigC)
|
||||
end subroutine
|
||||
124
src/GF/complex_cRGF2_self_energy.f90
Normal file
124
src/GF/complex_cRGF2_self_energy.f90
Normal file
@ -0,0 +1,124 @@
|
||||
subroutine complex_cRGF2_self_energy(eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
|
||||
|
||||
! Compute diagonal part of the GF2 self-energy and its renormalization factor
|
||||
|
||||
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
|
||||
complex*16,intent(in) :: e(nBas)
|
||||
complex*16,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,j,a,b
|
||||
integer :: p,q
|
||||
double precision :: eps
|
||||
double precision :: eta_tilde
|
||||
complex*16 :: num
|
||||
double precision,allocatable :: Re_DS(:)
|
||||
double precision,allocatable :: Im_DS(:)
|
||||
complex*16 :: z_dummy
|
||||
double precision,allocatable :: Re_SigC(:,:)
|
||||
double precision,allocatable :: Im_SigC(:,:)
|
||||
double precision,allocatable :: Re_Z(:)
|
||||
double precision,allocatable :: Im_Z(:)
|
||||
|
||||
! Output variables
|
||||
|
||||
complex*16,intent(out) :: SigC(nBas,nBas)
|
||||
complex*16,intent(out) :: Z(nBas)
|
||||
|
||||
! Initialize
|
||||
allocate(Re_DS(nBas),Im_DS(nBas),Re_SigC(nBas,nBas),Im_SigC(nBas,nBas),&
|
||||
Re_Z(nBas),Im_Z(nBas))
|
||||
Re_SigC(:,:) = 0d0
|
||||
Im_SigC(:,:) = 0d0
|
||||
Re_DS(:) = 0d0
|
||||
Im_DS(:) = 0d0
|
||||
|
||||
|
||||
! Compute GF2 self-energy
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP SHARED(Re_DS,Im_DS,Im_SigC,Re_SigC,ERI,eta,nC,nO,nBas,nR,e) &
|
||||
!$OMP PRIVATE(p,i,j,a,eps,num,eta_tilde,z_dummy) &
|
||||
!$OMP DEFAULT(NONE)
|
||||
!$OMP DO
|
||||
do p=nC+1,nBas-nR
|
||||
do q=nC+1,nBas-nR
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=nO+1,nBas-nR
|
||||
|
||||
eps = real(e(p)) + real(e(a)) - real(e(i)) - real(e(j))
|
||||
eta_tilde = eta - aimag(e(p)) + aimag(e(i)) - (aimag(e(a)) - aimag(e(j)))
|
||||
num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(q,a,i,j)
|
||||
z_dummy = num*cmplx(eps/(eps**2 + eta_tilde**2),eta_tilde/(eps**2 + eta_tilde**2),kind=8)
|
||||
Re_SigC(p,q) = Re_SigC(p,q) + real(z_dummy)
|
||||
Im_SigC(p,q) = Im_SigC(p,q) + aimag(z_dummy)
|
||||
if(p==q) then
|
||||
z_dummy = 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(z_dummy)
|
||||
Im_DS(p) = Im_DS(p) + aimag(z_dummy)
|
||||
end if
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP SHARED(Re_DS,Im_DS,Re_SigC,Im_SigC,ERI,eta,nC,nO,nBas,nR,e) &
|
||||
!$OMP PRIVATE(p,i,a,b,eps,num,eta_tilde,z_dummy) &
|
||||
!$OMP DEFAULT(NONE)
|
||||
!$OMP DO
|
||||
do p=nC+1,nBas-nR
|
||||
do q=nC+1,nBas-nR
|
||||
do i=nC+1,nO
|
||||
do a=nO+1,nBas-nR
|
||||
do b=nO+1,nBas-nR
|
||||
|
||||
eps = real(e(p)) + real(e(i)) - real(e(a)) - real(e(b))
|
||||
eta_tilde = eta + aimag(e(p)) - aimag(e(a)) - aimag(e(b)) + aimag(e(i))
|
||||
num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(q,i,a,b)
|
||||
|
||||
z_dummy = num*cmplx(eps/(eps**2 + eta_tilde**2),-eta_tilde/(eps**2 + eta_tilde**2),kind=8)
|
||||
Re_SigC(p,q) = Re_SigC(p,q) + real(z_dummy)
|
||||
Im_SigC(p,q) = Im_SigC(p,q) + aimag(z_dummy)
|
||||
if(p==q) then
|
||||
z_dummy = 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(z_dummy)
|
||||
Im_DS(p) = Im_DS(p) + aimag(z_dummy)
|
||||
end if
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
Re_Z(:) = (1d0-Re_DS(:))/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
|
||||
Im_Z(:) = Im_DS(:)/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
|
||||
|
||||
Z = cmplx(Re_Z,Im_Z,kind=8)
|
||||
SigC = cmplx(Re_SigC,Im_SigC,kind=8)
|
||||
|
||||
|
||||
deallocate(Re_DS,Im_DS,Re_Z,Im_Z,Re_SigC,Im_SigC)
|
||||
end subroutine
|
||||
346
src/GF/complex_qsRGF2.f90
Normal file
346
src/GF/complex_qsRGF2.f90
Normal file
@ -0,0 +1,346 @@
|
||||
subroutine complex_qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA, &
|
||||
dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc, &
|
||||
rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, &
|
||||
ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF, &
|
||||
CAP_AO,CAP_MO)
|
||||
|
||||
! Perform a quasiparticle self-consistent GF2 calculation
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
logical,intent(in) :: dotest
|
||||
|
||||
integer,intent(in) :: maxSCF
|
||||
integer,intent(in) :: max_diis
|
||||
double precision,intent(in) :: thresh
|
||||
logical,intent(in) :: dophBSE
|
||||
logical,intent(in) :: doppBSE
|
||||
logical,intent(in) :: TDA
|
||||
logical,intent(in) :: dBSE
|
||||
logical,intent(in) :: dTDA
|
||||
logical,intent(in) :: singlet
|
||||
logical,intent(in) :: triplet
|
||||
double precision,intent(in) :: eta
|
||||
logical,intent(in) :: regularize
|
||||
|
||||
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,nOrb,nC,nO,nV,nR,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) :: CAP_AO(nBas,nBas)
|
||||
complex*16,intent(inout) :: CAP_MO(nBas,nBas)
|
||||
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
|
||||
complex*16,intent(inout) :: ERI_MO(nOrb,nOrb,nOrb,nOrb)
|
||||
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
|
||||
complex*16,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: nSCF
|
||||
integer :: nBas_Sq
|
||||
integer :: ispin
|
||||
integer :: n_diis
|
||||
complex*16 :: EqsGF2
|
||||
double precision :: Conv
|
||||
double precision :: flow
|
||||
double precision :: rcond
|
||||
complex*16,external :: complex_trace_matrix
|
||||
complex*16 :: dipole(ncart)
|
||||
complex*16 :: ET
|
||||
complex*16 :: EV
|
||||
complex*16 :: EW
|
||||
complex*16 :: EJ
|
||||
complex*16 :: Ex
|
||||
complex*16 :: Ec
|
||||
complex*16 :: EcBSE(nspin)
|
||||
|
||||
complex*16,allocatable :: error_diis(:,:)
|
||||
complex*16,allocatable :: F_diis(:,:)
|
||||
complex*16,allocatable :: c(:,:)
|
||||
complex*16,allocatable :: cp(:,:)
|
||||
complex*16,allocatable :: eGF(:)
|
||||
complex*16,allocatable :: eOld(:)
|
||||
complex*16,allocatable :: P(:,:)
|
||||
complex*16,allocatable :: F(:,:)
|
||||
complex*16,allocatable :: Fp(:,:)
|
||||
complex*16,allocatable :: J(:,:)
|
||||
complex*16,allocatable :: K(:,:)
|
||||
complex*16,allocatable :: SigC(:,:)
|
||||
complex*16,allocatable :: SigCp(:,:)
|
||||
complex*16,allocatable :: Z(:)
|
||||
complex*16,allocatable :: error(:,:)
|
||||
|
||||
! Hello world
|
||||
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'********************************'
|
||||
write(*,*)'* Restricted qsGF2 Calculation *'
|
||||
write(*,*)'********************************'
|
||||
write(*,*)
|
||||
|
||||
! Warning
|
||||
|
||||
write(*,*) '!! ERIs in MO basis will be overwritten in qsGF2 !!'
|
||||
write(*,*)
|
||||
|
||||
! Stuff
|
||||
|
||||
nBas_Sq = nBas*nBas
|
||||
flow = 500d0
|
||||
! TDA
|
||||
|
||||
if(TDA) then
|
||||
write(*,*) 'Tamm-Dancoff approximation activated!'
|
||||
write(*,*)
|
||||
end if
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(eGF(nOrb))
|
||||
allocate(eOld(nOrb))
|
||||
|
||||
allocate(c(nBas,nOrb))
|
||||
|
||||
allocate(cp(nOrb,nOrb))
|
||||
allocate(Fp(nOrb,nOrb))
|
||||
|
||||
allocate(P(nBas,nBas))
|
||||
allocate(F(nBas,nBas))
|
||||
allocate(J(nBas,nBas))
|
||||
allocate(K(nBas,nBas))
|
||||
allocate(error(nBas,nBas))
|
||||
|
||||
allocate(Z(nOrb))
|
||||
allocate(SigC(nOrb,nOrb))
|
||||
|
||||
allocate(SigCp(nBas,nBas))
|
||||
|
||||
allocate(error_diis(nBas_Sq,max_diis))
|
||||
allocate(F_diis(nBas_Sq,max_diis))
|
||||
|
||||
! Initialization
|
||||
|
||||
nSCF = -1
|
||||
n_diis = 0
|
||||
ispin = 1
|
||||
Conv = 1d0
|
||||
P(:,:) = PHF(:,:)
|
||||
eOld(:) = eHF(:)
|
||||
eGF(:) = eHF(:)
|
||||
c(:,:) = cHF(:,:)
|
||||
F_diis(:,:) = 0d0
|
||||
error_diis(:,:) = 0d0
|
||||
rcond = 0d0
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Main loop
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
do while(Conv > thresh .and. nSCF <= maxSCF)
|
||||
|
||||
! Increment
|
||||
|
||||
nSCF = nSCF + 1
|
||||
|
||||
! Buid Hartree matrix
|
||||
|
||||
call complex_Hartree_matrix_AO_basis(nBas, P, ERI_AO, J)
|
||||
|
||||
! Compute exchange part of the self-energy
|
||||
|
||||
call complex_exchange_matrix_AO_basis(nBas, P, ERI_AO, K)
|
||||
|
||||
! AO to MO transformation of two-electron integrals
|
||||
|
||||
call complex_AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO)
|
||||
|
||||
! Compute self-energy and renormalization factor
|
||||
|
||||
if(regularize) then
|
||||
|
||||
call complex_cRGF2_reg_self_energy(flow,eta, nOrb, nC, nO, nV, nR, eGF, ERI_MO, SigC, Z)
|
||||
|
||||
else
|
||||
|
||||
call complex_cRGF2_self_energy(eta, nOrb, nC, nO, nV, nR, eGF, ERI_MO, SigC, Z)
|
||||
|
||||
end if
|
||||
|
||||
! Make correlation self-energy Hermitian and transform it back to AO basis
|
||||
|
||||
SigC = 0.5d0*(SigC + transpose(SigC))
|
||||
|
||||
call complex_MOtoAO(nBas, nOrb, S, c, SigC, SigCp)
|
||||
|
||||
! Solve the quasi-particle equation
|
||||
|
||||
F(:,:) = cmplx(Hc(:,:),CAP_AO(:,:),kind=8) + J(:,:) + 0.5d0*K(:,:) + SigCp(:,:)
|
||||
if(nBas .ne. nOrb) then
|
||||
call complex_complex_AOtoMO(nBas, nOrb, c(1,1), F(1,1), Fp(1,1))
|
||||
call complex_MOtoAO(nBas, nOrb, S(1,1), c(1,1), Fp(1,1), F(1,1))
|
||||
endif
|
||||
|
||||
! Compute commutator and convergence criteria
|
||||
|
||||
error = matmul(F, matmul(P, S)) - matmul(matmul(S, P), F)
|
||||
|
||||
! DIIS extrapolation
|
||||
|
||||
n_diis = min(n_diis+1, max_diis)
|
||||
if(abs(rcond) > 1d-7) then
|
||||
call complex_DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,error_diis,F_diis,error,F)
|
||||
else
|
||||
n_diis = 0
|
||||
end if
|
||||
|
||||
! Diagonalize Hamiltonian in AO basis
|
||||
|
||||
if(nBas .eq. nOrb) then
|
||||
Fp = matmul(transpose(X), matmul(F, X))
|
||||
cp(:,:) = Fp(:,:)
|
||||
call complex_diagonalize_matrix(nOrb, cp, eGF)
|
||||
c = matmul(X, cp)
|
||||
else
|
||||
Fp = matmul(transpose(c), matmul(F, c))
|
||||
cp(:,:) = Fp(:,:)
|
||||
call complex_diagonalize_matrix(nOrb, cp, eGF)
|
||||
c = matmul(c, cp)
|
||||
endif
|
||||
|
||||
|
||||
! Compute new density matrix in the AO basis
|
||||
|
||||
P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO)))
|
||||
|
||||
! Save quasiparticles energy for next cycle
|
||||
|
||||
Conv = maxval(abs(eGF - eOld))
|
||||
eOld(:) = eGF(:)
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute total energy
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
! Kinetic energy
|
||||
|
||||
ET = complex_trace_matrix(nBas, matmul(P, T))
|
||||
|
||||
! Potential energy
|
||||
|
||||
EV = complex_trace_matrix(nBas, matmul(P, V))
|
||||
|
||||
! CAP
|
||||
|
||||
EW = complex_trace_matrix(nBas,matmul(P,(0d0,1d0)*CAP_AO))
|
||||
|
||||
! Hartree energy
|
||||
|
||||
EJ = 0.5d0*complex_trace_matrix(nBas, matmul(P, J))
|
||||
|
||||
! Exchange energy
|
||||
|
||||
Ex = 0.25d0*complex_trace_matrix(nBas, matmul(P, K))
|
||||
|
||||
! Correlation energy
|
||||
|
||||
!call RMP2(.false., regularize, nOrb, nC, nO, nV, nR, ERI_MO, ENuc, EqsGF2, eGF, Ec)
|
||||
|
||||
! Total energy
|
||||
|
||||
EqsGF2 = ET + EV + EJ + Ex + Ec
|
||||
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Print results
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
!call dipole_moment(nBas, P, nNuc, ZNuc, rNuc, dipole_int_AO, dipole)
|
||||
call print_complex_qsRGF2(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGF, &
|
||||
c, SigC, Z, ENuc, ET, EV,EW, EJ, Ex, Ec, EqsGF2, dipole)
|
||||
end do
|
||||
!------------------------------------------------------------------------
|
||||
! End main loop
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
! Did it actually converge?
|
||||
|
||||
if(nSCF == maxSCF+1) then
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)' Convergence failed '
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)
|
||||
|
||||
deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, error, error_diis, F_diis)
|
||||
stop
|
||||
|
||||
end if
|
||||
|
||||
! Deallocate memory
|
||||
|
||||
deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, error, error_diis, F_diis)
|
||||
|
||||
!! Perform phBSE@GF2 calculation
|
||||
!
|
||||
! if(dophBSE) then
|
||||
!
|
||||
! call RGF2_phBSE(TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, nC, nO, &
|
||||
! nV, nR, nS, ERI_MO, dipole_int_MO, eGF, EcBSE)
|
||||
!
|
||||
! write(*,*)
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGF2 correlation energy (singlet) =',EcBSE(1)
|
||||
! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGF2 correlation energy (triplet) =',EcBSE(2)
|
||||
! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGF2 correlation energy =',sum(EcBSE(:))
|
||||
! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGF2 total energy =',ENuc + EqsGF2 + sum(EcBSE(:))
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,*)
|
||||
!
|
||||
! end if
|
||||
|
||||
|
||||
! Perform ppBSE@GF2 calculation
|
||||
|
||||
if(doppBSE) then
|
||||
|
||||
call RGF2_ppBSE(TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, &
|
||||
nC, nO, nV, nR, ERI_MO, dipole_int_MO, eGF, EcBSE)
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGF2 correlation energy (singlet) =',EcBSE(1),' au'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGF2 correlation energy (triplet) =',3d0*EcBSE(2),' au'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGF2 correlation energy =',EcBSE(1) + 3d0*EcBSE(2),' au'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGF2 total energy =',ENuc + EqsGF2 + EcBSE(1) + 3d0*EcBSE(2),' au'
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
! Testing zone
|
||||
|
||||
if(dotest) then
|
||||
|
||||
call dump_test_value('R','qsGF2 correlation energy',Ec)
|
||||
call dump_test_value('R','qsGF2 HOMO energy',eGF(nO))
|
||||
call dump_test_value('R','qsGF2 LUMO energy',eGF(nO+1))
|
||||
|
||||
end if
|
||||
|
||||
end subroutine
|
||||
140
src/GF/print_complex_qsRGF2.f90
Normal file
140
src/GF/print_complex_qsRGF2.f90
Normal file
@ -0,0 +1,140 @@
|
||||
|
||||
! ---
|
||||
|
||||
subroutine print_complex_qsRGF2(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGW, c, SigC, &
|
||||
Z, ENuc, ET, EV,EW, EJ, EK, EcGM, EcRPA, EqsGW, dipole)
|
||||
|
||||
! Print useful information about qsRGW calculation
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nBas, nOrb
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nSCF
|
||||
double precision,intent(in) :: ENuc
|
||||
complex*16,intent(in) :: ET
|
||||
complex*16,intent(in) :: EV
|
||||
complex*16,intent(in) :: EW
|
||||
complex*16,intent(in) :: EJ
|
||||
complex*16,intent(in) :: EK
|
||||
complex*16,intent(in) :: EcGM
|
||||
complex*16,intent(in) :: EcRPA
|
||||
double precision,intent(in) :: Conv
|
||||
double precision,intent(in) :: thresh
|
||||
complex*16,intent(in) :: eHF(nOrb)
|
||||
complex*16,intent(in) :: eGW(nOrb)
|
||||
complex*16,intent(in) :: c(nBas,nOrb)
|
||||
complex*16,intent(in) :: SigC(nOrb,nOrb)
|
||||
complex*16,intent(in) :: Z(nOrb)
|
||||
complex*16,intent(in) :: EqsGW
|
||||
complex*16,intent(in) :: dipole(ncart)
|
||||
|
||||
! Local variables
|
||||
|
||||
logical :: dump_orb = .false.
|
||||
integer :: p,ixyz,HOMO,LUMO
|
||||
complex*16 :: Gap
|
||||
double precision,external :: complex_trace_matrix
|
||||
|
||||
! Output variables
|
||||
|
||||
! HOMO and LUMO
|
||||
|
||||
HOMO = maxloc(real(eGW(1:nO)),1)
|
||||
LUMO = minloc(real(eGW(nO+1:nBas)),1) + nO
|
||||
Gap = eGW(LUMO)-eGW(HOMO)
|
||||
|
||||
! Compute energies
|
||||
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)' Self-consistent qsGF2 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,nOrb
|
||||
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,'|',real(eHF(p))*HaToeV,'|',real(SigC(p,p))*HaToeV,'|',real(Z(p)),'|',real(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,'|',aimag(eHF(p))*HaToeV,'|',aimag(SigC(p,p))*HaToeV,'|',aimag(Z(p)),'|',aimag(eGW(p))*HaToeV,'|'
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
if(p==nO) then
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A10,I3)') 'Iteration ',nSCF
|
||||
write(*,'(2X,A14,F15.5)')'Convergence = ',Conv
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A60,F15.6,A3)') 'qsGW@RHF HOMO real energy = ',real(eGW(HOMO))*HaToeV,' eV'
|
||||
write(*,'(2X,A60,F15.6,A3)') 'qsGW@RHF HOMO imag energy = ',aimag(eGW(HOMO))*HaToeV,' eV'
|
||||
write(*,'(2X,A60,F15.6,A3)') 'qsGW@RHF LUMO real energy = ',real(eGW(LUMO))*HaToeV,' eV'
|
||||
write(*,'(2X,A60,F15.6,A3)') 'qsGW@RHF LUMO imag energy = ',aimag(eGW(LUMO))*HaToeV,' eV'
|
||||
write(*,'(2X,A60,F15.6,A3)') 'qsGW@RHF HOMO-LUMO gap = ',real(Gap)*HaToeV,' eV'
|
||||
write(*,'(2X,A60,F15.6,A3)') 'qsGW@RHF HOMO-LUMO gap = ',aimag(Gap)*HaToeV,' eV'
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A60,F15.6,A3)') ' qsGW@RHF total real energy = ',ENuc + real(EqsGW),' au'
|
||||
write(*,'(2X,A60,F15.6,A3)') ' qsGW@RHF total imag energy = ',aimag(EqsGW),' au'
|
||||
write(*,'(2X,A60,F15.6,A3)') ' qsGW@RHF exchange energy = ',real(EK),' au'
|
||||
write(*,'(2X,A60,F15.6,A3)') ' qsGW@RHF exchange energy = ',aimag(EK),' au'
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
! Dump results for final iteration
|
||||
|
||||
if(Conv < thresh) then
|
||||
|
||||
write(*,*)
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
write(*,'(A33)') ' Summary '
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' One-electron energy = ',real(ET) + real(EV) + real(EW),' au'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' One-electron energy = ',aimag(ET) + aimag(EV) + aimag(EW),' au'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' Kinetic energy = ',real(ET),' au'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' Kinetic energy = ',aimag(ET),' au'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' Potential energy = ',real(EV),' au'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' Potential energy = ',aimag(EV),' au'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' CAP energy = ',real(EW),' au'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' CAP energy = ',aimag(EW),' au'
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' Two-electron energy = ',real(EJ + EK),' au'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' Two-electron energy = ',aimag(EJ + EK),' au'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' Hartree energy = ',real(EJ),' au'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' Hartree energy = ',aimag(EJ),' au'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' Exchange energy = ',real(EK),' au'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' Exchange energy = ',aimag(EK),' au'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' Correlation energy = ',real(EcGM),' au'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' Correlation energy = ',aimag(EcGM),' au'
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' Electronic energy = ',real(EqsGW),' au'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' Electronic energy = ',aimag(EqsGW),' au'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' Nuclear repulsion = ',ENuc,' au'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' qsRGW energy = ',ENuc + real(EqsGW),' au'
|
||||
write(*,'(A33,1X,F16.10,A3)') ' qsRGW energy = ',aimag(EqsGW),' au'
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
if(dump_orb) then
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
write(*,'(A50)') ' Restricted qsGW orbital coefficients'
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
call complex_matout(nBas, nOrb, c)
|
||||
write(*,*)
|
||||
end if
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
write(*,'(A50)') ' Restricted qsGW orbital energies (au) '
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
call complex_vecout(nOrb, eGW)
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
end subroutine
|
||||
@ -297,7 +297,7 @@ subroutine complex_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,d
|
||||
|
||||
! Print results
|
||||
|
||||
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
|
||||
!call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
|
||||
call print_complex_qsRGW(nBas,nOrb,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z, &
|
||||
ENuc,ET,EV,EW,EJ,EK,EcGM,EcRPA,EqsGW,dipole)
|
||||
|
||||
|
||||
@ -391,11 +391,11 @@ doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3
|
||||
|
||||
if(doGF .and. docRHF) then
|
||||
call wall_time(start_GF)
|
||||
call complex_RGF(dotest,docG0F2,doevGF2,maxSCF_GF, &
|
||||
call complex_RGF(dotest,docG0F2,doevGF2,doqsGF2,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,complex_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)
|
||||
S,X,T,V,Hc,ERI_AO,complex_ERI_MO,dipole_int_AO,complex_dipole_int_MO,&
|
||||
complex_PHF,complex_cHF,complex_eHF,CAP_AO, complex_CAP_MO)
|
||||
call wall_time(end_GF)
|
||||
|
||||
t_GF = end_GF - start_GF
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user