mirror of
https://github.com/pfloos/quack
synced 2025-05-06 15:14:55 +02:00
first qsGW version + cleaner data allocation to avoid warnings
This commit is contained in:
Loris Burth
parent
f7df08e82d
commit
cdc1c25351
39
src/AOtoMO/complex_MOtoAO.f90
Normal file
39
src/AOtoMO/complex_MOtoAO.f90
Normal file
@ -0,0 +1,39 @@
|
||||
subroutine complex_MOtoAO(nBas, nOrb, S, C, M_MOs, M_AOs)
|
||||
|
||||
! Perform MO to AO transformation of a matrix M_AOs for a given metric S
|
||||
! and coefficients c
|
||||
!
|
||||
! M_AOs = S C M_MOs (S C).T
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: nBas, nOrb
|
||||
double precision, intent(in) :: S(nBas,nBas)
|
||||
complex*16, intent(in) :: C(nBas,nOrb)
|
||||
complex*16, intent(in) :: M_MOs(nOrb,nOrb)
|
||||
complex*16, intent(out) :: M_AOs(nBas,nBas)
|
||||
|
||||
complex*16, allocatable :: SC(:,:),BSC(:,:),cS(:,:)
|
||||
|
||||
|
||||
allocate(SC(nBas,nOrb), BSC(nOrb,nBas),cS(nBas,nBas))
|
||||
cS(:,:) = (0d0,1d0)*S(:,:)
|
||||
!SC = matmul(S, C)
|
||||
!BSC = matmul(M_MOs, transpose(SC))
|
||||
!M_AOs = matmul(SC, BSC)
|
||||
|
||||
call zgemm("N", "N", nBas, nOrb, nBas, 1.d0, &
|
||||
cS(1,1), nBas, C(1,1), nBas, &
|
||||
0.d0, SC(1,1), nBas)
|
||||
|
||||
call zgemm("N", "T", nOrb, nBas, nOrb, 1.d0, &
|
||||
M_MOs(1,1), nOrb, SC(1,1), nBas, &
|
||||
0.d0, BSC(1,1), nOrb)
|
||||
|
||||
call zgemm("N", "N", nBas, nBas, nOrb, 1.d0, &
|
||||
SC(1,1), nBas, BSC(1,1), nOrb, &
|
||||
0.d0, M_AOs(1,1), nBas)
|
||||
|
||||
deallocate(SC, BSC)
|
||||
|
||||
end subroutine
|
||||
24
src/AOtoMO/complex_complex_AOtoMO.f90
Normal file
24
src/AOtoMO/complex_complex_AOtoMO.f90
Normal file
@ -0,0 +1,24 @@
|
||||
subroutine complex_complex_AOtoMO(nBas, nOrb, C, M_AOs, M_MOs)
|
||||
|
||||
! Perform AO to MO transformation of a matrix M_AOs for given coefficients c
|
||||
! M_MOs = C.T M_AOs C
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: nBas, nOrb
|
||||
complex*16, intent(in) :: C(nBas,nOrb)
|
||||
complex*16, intent(in) :: M_AOs(nBas,nBas)
|
||||
|
||||
complex*16, intent(out) :: M_MOs(nOrb,nOrb)
|
||||
|
||||
complex*16, allocatable :: AC(:,:)
|
||||
complex*16, allocatable :: complex_C(:,:)
|
||||
|
||||
allocate(AC(nBas,nOrb))
|
||||
|
||||
AC = matmul(M_AOs, C)
|
||||
M_MOs = matmul(transpose(C), AC)
|
||||
|
||||
deallocate(AC)
|
||||
|
||||
end subroutine
|
||||
@ -1,7 +1,7 @@
|
||||
subroutine complex_RGW(dotest,docG0W0,doevGW,maxSCF,thresh,max_diis,doACFDT, &
|
||||
subroutine complex_RGW(dotest,docG0W0,doevGW,doqsGW,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)
|
||||
S,X,T,V,Hc,ERI_AO,ERI_MO,CAP_AO,CAP_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
|
||||
|
||||
! Restricted GW module
|
||||
|
||||
@ -12,7 +12,7 @@ subroutine complex_RGW(dotest,docG0W0,doevGW,maxSCF,thresh,max_diis,doACFDT,
|
||||
|
||||
logical,intent(in) :: dotest
|
||||
|
||||
logical,intent(in) :: docG0W0,doevGW
|
||||
logical,intent(in) :: docG0W0,doevGW,doqsGW
|
||||
|
||||
integer,intent(in) :: maxSCF
|
||||
integer,intent(in) :: max_diis
|
||||
@ -57,6 +57,7 @@ subroutine complex_RGW(dotest,docG0W0,doevGW,maxSCF,thresh,max_diis,doACFDT,
|
||||
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)
|
||||
double precision,intent(in) :: CAP_AO(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)
|
||||
@ -100,4 +101,24 @@ subroutine complex_RGW(dotest,docG0W0,doevGW,maxSCF,thresh,max_diis,doACFDT,
|
||||
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Perform qsGW calculation
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(doqsGW) then
|
||||
|
||||
call wall_time(start_GW)
|
||||
call complex_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, &
|
||||
TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet,eta,doSRG,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_GW)
|
||||
|
||||
t_GW = end_GW - start_GW
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_GW,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
end subroutine
|
||||
|
||||
99
src/GW/complex_RGW_SRG_self_energy_diag.f90
Normal file
99
src/GW/complex_RGW_SRG_self_energy_diag.f90
Normal file
@ -0,0 +1,99 @@
|
||||
subroutine complex_RGW_SRG_self_energy_diag(flow,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
|
||||
double precision,intent(in) :: flow
|
||||
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),&
|
||||
-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/(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
|
||||
178
src/GW/complex_RGW_self_energy.f90
Normal file
178
src/GW/complex_RGW_self_energy.f90
Normal file
@ -0,0 +1,178 @@
|
||||
subroutine complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho,EcGM,Re_Sig,Im_Sig,Re_Z,Im_Z)
|
||||
|
||||
! Compute 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(nOrb)
|
||||
double precision,intent(in) :: Im_e(nOrb)
|
||||
complex*16,intent(in) :: Om(nS)
|
||||
complex*16,intent(in) :: rho(nOrb,nOrb,nS)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,j,a,b
|
||||
integer :: p,q,m
|
||||
double precision :: eps,eta_tilde
|
||||
complex*16 :: num,tmp
|
||||
double precision, allocatable :: Re_DS(:)
|
||||
double precision, allocatable :: Im_DS(:)
|
||||
|
||||
! Output variables
|
||||
|
||||
complex*16,intent(out) :: EcGM
|
||||
double precision,intent(out) :: Re_Sig(nOrb,nOrb)
|
||||
double precision,intent(out) :: Im_Sig(nOrb,nOrb)
|
||||
double precision,intent(out) :: Re_Z(nOrb)
|
||||
double precision,intent(out) :: Im_Z(nOrb)
|
||||
|
||||
!----------------!
|
||||
! GW self-energy !
|
||||
!----------------!
|
||||
allocate(Re_DS(nBas),Im_DS(nBas))
|
||||
|
||||
Re_Sig(:,:) = 0d0
|
||||
Im_Sig(:,:) = 0d0
|
||||
|
||||
! Occupied part of the correlation self-energy
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP SHARED(Re_Sig,Im_Sig,rho,eta,nS,nC,nO,nOrb,nR,Re_e,Im_e,Om) &
|
||||
!$OMP PRIVATE(m,i,q,p,eps,num,eta_tilde,tmp) &
|
||||
!$OMP DEFAULT(NONE)
|
||||
!$OMP DO
|
||||
do q=nC+1,nOrb-nR
|
||||
do p=nC+1,nOrb-nR
|
||||
do m=1,nS
|
||||
do i=nC+1,nO
|
||||
|
||||
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)*rho(q,i,m)
|
||||
tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),&
|
||||
eta_tilde/(eps**2+eta_tilde**2),kind=8)
|
||||
Re_Sig(p,q) = Re_Sig(p,q) + real(tmp)
|
||||
Im_Sig(p,q) = Im_Sig(p,q) + aimag(tmp)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
! Virtual part of the correlation self-energy
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP SHARED(Re_Sig,Im_Sig,rho,eta,nS,nC,nO,nOrb,nR,Re_e,Im_e,Om) &
|
||||
!$OMP PRIVATE(m,a,q,p,eps,num,eta_tilde,tmp) &
|
||||
!$OMP DEFAULT(NONE)
|
||||
!$OMP DO
|
||||
do q=nC+1,nOrb-nR
|
||||
do p=nC+1,nOrb-nR
|
||||
do m=1,nS
|
||||
do a=nO+1,nOrb-nR
|
||||
|
||||
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)*rho(q,a,m)
|
||||
tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),&
|
||||
-eta_tilde/(eps**2 + eta_tilde**2),kind=8)
|
||||
Re_Sig(p,q) = Re_Sig(p,q) + real(tmp)
|
||||
Im_Sig(p,q) = Im_Sig(p,q) + aimag(tmp)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
!------------------------!
|
||||
! Renormalization factor !
|
||||
!------------------------!
|
||||
|
||||
Re_DS(:) = 0d0
|
||||
Im_DS(:) = 0d0
|
||||
|
||||
! Occupied part of the renormalization factor
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP SHARED(Re_DS,Im_DS,rho,eta,nS,nC,nO,nOrb,nR,Re_e,Im_e,Om) &
|
||||
!$OMP PRIVATE(m,i,p,eps,num,eta_tilde,tmp) &
|
||||
!$OMP DEFAULT(NONE)
|
||||
!$OMP DO
|
||||
do p=nC+1,nOrb-nR
|
||||
do m=1,nS
|
||||
do i=nC+1,nO
|
||||
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)*rho(p,i,m)
|
||||
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
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
! Virtual part of the renormalization factor
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP SHARED(Re_DS,Im_DS,rho,eta,nS,nC,nO,nOrb,nR,Re_e,Im_e,Om) &
|
||||
!$OMP PRIVATE(m,a,p,eps,num,eta_tilde,tmp) &
|
||||
!$OMP DEFAULT(NONE)
|
||||
!$OMP DO
|
||||
do p=nC+1,nOrb-nR
|
||||
do m=1,nS
|
||||
do a=nO+1,nOrb-nR
|
||||
|
||||
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)*rho(p,a,m)
|
||||
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
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
! 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)
|
||||
|
||||
!!-------------------------------------!
|
||||
!! Galitskii-Migdal correlation energy !
|
||||
!!-------------------------------------!
|
||||
!
|
||||
! EcGM = 0d0
|
||||
! do m=1,nS
|
||||
! do a=nO+1,nOrb-nR
|
||||
! do i=nC+1,nO
|
||||
!
|
||||
! eps = e(a) - e(i) + Om(m)
|
||||
! num = 4d0*rho(a,i,m)*rho(a,i,m)
|
||||
! EcGM = EcGM - num*eps/(eps**2 + eta**2)
|
||||
!
|
||||
! end do
|
||||
! end do
|
||||
! end do
|
||||
!
|
||||
end subroutine
|
||||
@ -14,7 +14,7 @@ subroutine complex_evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,d
|
||||
integer,intent(in) :: max_diis
|
||||
double precision,intent(in) :: thresh
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: ERHF
|
||||
complex*16,intent(in) :: ERHF
|
||||
logical,intent(in) :: doACFDT
|
||||
logical,intent(in) :: exchange_kernel
|
||||
logical,intent(in) :: doXBS
|
||||
|
||||
388
src/GW/complex_qsRGW.f90
Normal file
388
src/GW/complex_qsRGW.f90
Normal file
@ -0,0 +1,388 @@
|
||||
subroutine complex_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, &
|
||||
TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet,eta,doSRG,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 GW 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) :: 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
|
||||
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,nBas)
|
||||
double precision,intent(in) :: CAP_AO(nBas,nBas)
|
||||
double precision,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(inout) :: dipole_int_MO(nOrb,nOrb,ncart)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: nSCF
|
||||
integer :: nBas_Sq
|
||||
integer :: ispin
|
||||
integer :: ixyz
|
||||
integer :: n_diis
|
||||
complex*16 :: ET
|
||||
complex*16 :: EV
|
||||
complex*16 :: EJ
|
||||
complex*16 :: EK
|
||||
complex*16 :: EqsGW
|
||||
complex*16 :: EW
|
||||
complex*16 :: EcRPA
|
||||
complex*16 :: EcBSE(nspin)
|
||||
complex*16 :: EcGM
|
||||
double precision :: Conv
|
||||
double precision :: rcond
|
||||
double precision,external :: complex_trace_matrix
|
||||
complex*16 :: dipole(ncart)
|
||||
|
||||
double precision :: flow
|
||||
|
||||
logical :: dRPA_W = .true.
|
||||
logical :: print_W = .false.
|
||||
complex*16,allocatable :: err_diis(:,:)
|
||||
complex*16,allocatable :: F_diis(:,:)
|
||||
complex*16,allocatable :: Aph(:,:)
|
||||
complex*16,allocatable :: Bph(:,:)
|
||||
complex*16,allocatable :: Om(:)
|
||||
complex*16,allocatable :: XpY(:,:)
|
||||
complex*16,allocatable :: XmY(:,:)
|
||||
complex*16,allocatable :: rho(:,:,:)
|
||||
complex*16,allocatable :: c(:,:)
|
||||
complex*16,allocatable :: cp(:,:)
|
||||
complex*16,allocatable :: eGW(:)
|
||||
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 :: err(:,:)
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'*******************************'
|
||||
write(*,*)'* Restricted qsGW Calculation *'
|
||||
write(*,*)'*******************************'
|
||||
write(*,*)
|
||||
|
||||
! Warning
|
||||
|
||||
write(*,*) '!! ERIs and CAP in MO basis will be overwritten in qsGW !!'
|
||||
write(*,*)
|
||||
|
||||
! Stuff
|
||||
|
||||
nBas_Sq = nBas*nBas
|
||||
|
||||
! TDA for W
|
||||
|
||||
if(TDA_W) then
|
||||
write(*,*) 'Tamm-Dancoff approximation for dynamical screening!'
|
||||
write(*,*)
|
||||
end if
|
||||
|
||||
! SRG regularization
|
||||
|
||||
flow = 500d0
|
||||
|
||||
if(doSRG) then
|
||||
|
||||
write(*,*) '*** SRG regularized qsGW scheme ***'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(eGW(nOrb))
|
||||
allocate(Z(nOrb))
|
||||
|
||||
allocate(c(nBas,nOrb))
|
||||
|
||||
allocate(cp(nOrb,nOrb))
|
||||
allocate(Fp(nOrb,nOrb))
|
||||
allocate(SigC(nOrb,nOrb))
|
||||
|
||||
allocate(P(nBas,nBas))
|
||||
allocate(F(nBas,nBas))
|
||||
allocate(J(nBas,nBas))
|
||||
allocate(K(nBas,nBas))
|
||||
allocate(err(nBas,nBas))
|
||||
allocate(SigCp(nBas,nBas))
|
||||
|
||||
allocate(Aph(nS,nS))
|
||||
allocate(Bph(nS,nS))
|
||||
allocate(Om(nS))
|
||||
allocate(XpY(nS,nS))
|
||||
allocate(XmY(nS,nS))
|
||||
allocate(rho(nOrb,nOrb,nS))
|
||||
|
||||
allocate(err_diis(nBas_Sq,max_diis))
|
||||
allocate(F_diis(nBas_Sq,max_diis))
|
||||
|
||||
! Initialization
|
||||
|
||||
nSCF = -1
|
||||
n_diis = 0
|
||||
ispin = 1
|
||||
Conv = 1d0
|
||||
P(:,:) = PHF(:,:)
|
||||
eGW(:) = eHF(:)
|
||||
c(:,:) = cHF(:,:)
|
||||
F_diis(:,:) = 0d0
|
||||
err_diis(:,:) = 0d0
|
||||
rcond = 0d0
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Main loop
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
do while(Conv > thresh .and. nSCF <= maxSCF)
|
||||
|
||||
! Increment
|
||||
|
||||
nSCF = nSCF + 1
|
||||
|
||||
! Build Hartree-exchange matrix
|
||||
|
||||
call complex_Hartree_matrix_AO_basis(nBas,P,ERI_AO,J)
|
||||
call complex_exchange_matrix_AO_basis(nBas,P,ERI_AO,K)
|
||||
|
||||
! AO to MO transformation of two-electron integrals
|
||||
|
||||
do ixyz=1,ncart
|
||||
call complex_AOtoMO(nBas,nOrb,c,dipole_int_AO(1,1,ixyz),dipole_int_MO(1,1,ixyz))
|
||||
end do
|
||||
|
||||
call complex_AOtoMO_ERI_RHF(nBas,nOrb,c,ERI_AO,ERI_MO)
|
||||
|
||||
! Compute linear response
|
||||
|
||||
call complex_phRLR_A(ispin,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph)
|
||||
if(.not.TDA_W) call complex_phRLR_B(ispin,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
|
||||
|
||||
call complex_phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
|
||||
if(print_W) call print_excitation_energies('phRPA@GW@RHF','singlet',nS,Om)
|
||||
|
||||
call complex_RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI_MO,XpY,rho)
|
||||
|
||||
if(doSRG) then
|
||||
write(*,*) "SRG not implemented"
|
||||
!call complex_RGW_SRG_self_energy(flow,nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
|
||||
else
|
||||
call complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,real(eGW),aimag(eGW),Om,rho,&
|
||||
EcGM,real(SigC),aimag(SigC),real(Z),aimag(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
|
||||
|
||||
err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
|
||||
|
||||
if(nSCF > 1) Conv = maxval(abs(err))
|
||||
|
||||
! Kinetic energy
|
||||
|
||||
ET = complex_trace_matrix(nBas,matmul(P,T))
|
||||
|
||||
! Potential energy
|
||||
|
||||
EV = complex_trace_matrix(nBas,matmul(P,V))
|
||||
|
||||
! Hartree energy
|
||||
|
||||
EJ = 0.5d0*complex_trace_matrix(nBas,matmul(P,J))
|
||||
|
||||
! Exchange energy
|
||||
|
||||
EK = 0.25d0*complex_trace_matrix(nBas,matmul(P,K))
|
||||
|
||||
! CAP energy
|
||||
|
||||
EW = complex_trace_matrix(nBas,matmul(P,(0d0,1d0)*CAP_AO))
|
||||
|
||||
! Total energy
|
||||
|
||||
EqsGW = ET + EV + EJ + EK + EW
|
||||
|
||||
! DIIS extrapolation
|
||||
|
||||
if(max_diis > 1) then
|
||||
|
||||
n_diis = min(n_diis+1,max_diis)
|
||||
call complex_DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,F)
|
||||
|
||||
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,eGW)
|
||||
c = matmul(X,cp)
|
||||
else
|
||||
Fp = matmul(transpose(c),matmul(F,c))
|
||||
cp(:,:) = Fp(:,:)
|
||||
call complex_diagonalize_matrix(nOrb,cp,eGW)
|
||||
c = matmul(c,cp)
|
||||
endif
|
||||
|
||||
call complex_complex_AOtoMO(nBas,nOrb,c,SigCp,SigC)
|
||||
|
||||
! Density matrix
|
||||
|
||||
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
|
||||
|
||||
! Print results
|
||||
|
||||
!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)
|
||||
|
||||
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,Om,XpY,XmY,rho,err,err_diis,F_diis)
|
||||
! stop
|
||||
!
|
||||
! end if
|
||||
!
|
||||
!! Deallocate memory
|
||||
!
|
||||
! deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,Z,Om,XpY,XmY,rho,err,err_diis,F_diis)
|
||||
!
|
||||
!! Perform BSE calculation
|
||||
!
|
||||
! if(dophBSE) then
|
||||
!
|
||||
! call RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta, &
|
||||
! nOrb,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE)
|
||||
!
|
||||
! write(*,*)
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@qsGW@RHF correlation energy (singlet) = ',EcBSE(1),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@qsGW@RHF correlation energy (triplet) = ',EcBSE(2),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@qsGW@RHF correlation energy = ',sum(EcBSE),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@qsGW@RHF total energy = ',ENuc + EqsGW + sum(EcBSE),' au'
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,*)
|
||||
!
|
||||
!! Compute the BSE correlation energy via the adiabatic connection
|
||||
!
|
||||
! if(doACFDT) then
|
||||
!
|
||||
! call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcBSE)
|
||||
!
|
||||
! write(*,*)
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@qsGW@RHF correlation energy (singlet) = ',EcBSE(1),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@qsGW@RHF correlation energy (triplet) = ',EcBSE(2),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@qsGW@RHF correlation energy = ',sum(EcBSE),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@qsGW@RHF total energy = ',ENuc + EqsGW + sum(EcBSE),' au'
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,*)
|
||||
!
|
||||
! end if
|
||||
!
|
||||
! end if
|
||||
!
|
||||
! if(doppBSE) then
|
||||
!
|
||||
! call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eHF,eGW,EcBSE)
|
||||
!
|
||||
! write(*,*)
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGW@RHF correlation energy (singlet) = ',EcBSE(1),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGW@RHF correlation energy (triplet) = ',EcBSE(2),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGW@RHF correlation energy = ',sum(EcBSE),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGW@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au'
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,*)
|
||||
!
|
||||
! end if
|
||||
!
|
||||
!! Testing zone
|
||||
!
|
||||
! if(dotest) then
|
||||
!
|
||||
! call dump_test_value('R','qsGW correlation energy',EcRPA)
|
||||
! call dump_test_value('R','qsGW HOMO energy',eGW(nO))
|
||||
! call dump_test_value('R','qsGW LUMO energy',eGW(nO+1))
|
||||
!
|
||||
! end if
|
||||
!
|
||||
end subroutine
|
||||
148
src/GW/print_complex_qsRGW.f90
Normal file
148
src/GW/print_complex_qsRGW.f90
Normal file
@ -0,0 +1,148 @@
|
||||
|
||||
! ---
|
||||
|
||||
subroutine print_complex_qsRGW(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
|
||||
|
||||
! Dump results
|
||||
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
if(nSCF < 10) then
|
||||
write(*,'(1X,A20,I1,A1,I1,A17)')' Self-consistent qsG',nSCF,'W',nSCF,'@cRHF calculation'
|
||||
elseif(nSCF < 100) then
|
||||
write(*,'(1X,A20,I2,A1,I2,A17)')' Self-consistent qsG',nSCF,'W',nSCF,'@cRHF calculation'
|
||||
else
|
||||
write(*,'(1X,A20,I3,A1,I3,A17)')' Self-consistent qsG',nSCF,'W',nSCF,'@cRHF calculation'
|
||||
end if
|
||||
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
|
||||
@ -137,16 +137,28 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
|
||||
allocate(complex_cHF(nBas,nOrb))
|
||||
allocate(complex_FHF(nBas,nBas))
|
||||
allocate(complex_dipole_int_MO(nOrb,nOrb,ncart))
|
||||
allocate(dipole_int_MO(0,0,0))
|
||||
allocate(complex_ERI_MO(nOrb,nOrb,nOrb,nOrb))
|
||||
if (doCAP) allocate(complex_CAP_MO(nOrb,nOrb))
|
||||
allocate(CAP_MO(0,0))
|
||||
if (doCAP) then
|
||||
allocate(complex_CAP_MO(nOrb,nOrb))
|
||||
else
|
||||
allocate(complex_CAP_MO(0,0))
|
||||
end if
|
||||
else
|
||||
allocate(PHF(nBas,nBas))
|
||||
allocate(eHF(nOrb))
|
||||
allocate(cHF(nBas,nOrb))
|
||||
allocate(FHF(nBas,nBas))
|
||||
allocate(dipole_int_MO(nOrb,nOrb,ncart))
|
||||
allocate(complex_dipole_int_MO(0,0,0))
|
||||
allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb))
|
||||
if (doCAP) allocate(CAP_MO(nOrb,nOrb))
|
||||
allocate(complex_CAP_MO(0,0))
|
||||
if (doCAP) then
|
||||
allocate(CAP_MO(nOrb,nOrb))
|
||||
else
|
||||
allocate(CAP_MO(0,0))
|
||||
end if
|
||||
end if
|
||||
|
||||
allocate(ERI_AO(nBas,nBas,nBas,nBas))
|
||||
@ -418,10 +430,11 @@ doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3
|
||||
|
||||
if(doGW .and. docRHF) then
|
||||
call wall_time(start_GW)
|
||||
call complex_RGW(dotest,docG0W0,doevGW,maxSCF_GW,thresh_GW,max_diis_GW, &
|
||||
call complex_RGW(dotest,docG0W0,doevGW,doqsGW,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)
|
||||
V,Hc,ERI_AO,complex_ERI_MO,CAP_AO,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
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user