print qsUGW

This commit is contained in:
Pierre-Francois Loos 2020-10-21 20:53:36 +02:00
parent d6723f332d
commit e26c052a17
3 changed files with 57 additions and 49 deletions

View File

@ -9,7 +9,7 @@
# GF: maxSCF thresh DIIS n_diis lin eta renorm
256 0.00001 T 5 T 0.0 3
# GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0
128 0.00001 T 5 T 0.00367493 F F F F F
128 0.00001 T 5 T 0.00367493 F F T F F
# ACFDT: AC Kx XBS
F F T
# BSE: BSE dBSE dTDA evDyn

View File

@ -1,4 +1,5 @@
subroutine print_qsUGW(nBas,nO,Ov,nSCF,Conv,thresh,eGW,cGW,PGW,T,V,J,K,ENuc,EHF,SigC,Z,EcRPA,dipole)
subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eGW,cGW,PGW,Ov,T,V,J,K, &
ENuc,ET,EV,EJ,Ex,Ec,EcRPA,EqsGW,SigC,Z,dipole)
! Print one-electron energies and other stuff for qsUGW
@ -9,16 +10,21 @@ subroutine print_qsUGW(nBas,nO,Ov,nSCF,Conv,thresh,eGW,cGW,PGW,T,V,J,K,ENuc,EHF,
integer,intent(in) :: nBas
integer,intent(in) :: nO(nspin)
double precision,intent(in) :: Ov(nBas,nBas)
integer,intent(in) :: nSCF
double precision,intent(in) :: ENuc
double precision,intent(in) :: EHF
double precision,intent(in) :: ET(nspin)
double precision,intent(in) :: EV(nspin)
double precision,intent(in) :: EJ(nsp)
double precision,intent(in) :: Ex(nspin)
double precision,intent(in) :: Ec(nsp)
double precision,intent(in) :: EcRPA
double precision,intent(in) :: EqsGW
double precision,intent(in) :: Conv
double precision,intent(in) :: thresh
double precision,intent(in) :: eGW(nBas,nspin)
double precision,intent(in) :: cGW(nBas,nBas,nspin)
double precision,intent(in) :: PGW(nBas,nBas,nspin)
double precision,intent(in) :: Ov(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: J(nBas,nBas,nspin)
@ -34,12 +40,6 @@ subroutine print_qsUGW(nBas,nO,Ov,nSCF,Conv,thresh,eGW,cGW,PGW,T,V,J,K,ENuc,EHF,
double precision :: HOMO(nspin)
double precision :: LUMO(nspin)
double precision :: Gap(nspin)
double precision :: ET(nspin)
double precision :: EV(nspin)
double precision :: EJ(nsp)
double precision :: Ex(nspin)
double precision :: Ec(nsp)
double precision :: EqsGW
double precision :: S_exact,S2_exact
double precision :: S,S2
double precision,external :: trace_matrix
@ -64,44 +64,6 @@ subroutine print_qsUGW(nBas,nO,Ov,nSCF,Conv,thresh,eGW,cGW,PGW,T,V,J,K,ENuc,EHF,
S_exact = 0.5d0*dble(nO(1) - nO(2))
S = -0.5d0 + 0.5d0*sqrt(1d0 + 4d0*S2)
!------------------------------------------------------------------------
! Compute total energy
!------------------------------------------------------------------------
! Kinetic energy
do ispin=1,nspin
ET(ispin) = trace_matrix(nBas,matmul(PGW(:,:,ispin),T(:,:)))
end do
! Potential energy
do ispin=1,nspin
EV(ispin) = trace_matrix(nBas,matmul(PGW(:,:,ispin),V(:,:)))
end do
! Coulomb energy
EJ(1) = 0.5d0*trace_matrix(nBas,matmul(PGW(:,:,1),J(:,:,1)))
EJ(2) = trace_matrix(nBas,matmul(PGW(:,:,1),J(:,:,2)))
EJ(3) = 0.5d0*trace_matrix(nBas,matmul(PGW(:,:,2),J(:,:,2)))
! Exchange energy
do ispin=1,nspin
Ex(ispin) = 0.5d0*trace_matrix(nBas,matmul(PGW(:,:,ispin),K(:,:,ispin)))
end do
! Correlation energy
Ec(1) = 0.5d0*trace_matrix(nBas,matmul(PGW(:,:,1),SigC(:,:,1)))
Ec(2) = trace_matrix(nBas,matmul(PGW(:,:,1),SigC(:,:,2)))
Ec(3) = 0.5d0*trace_matrix(nBas,matmul(PGW(:,:,2),SigC(:,:,2)))
! Total energy
EqsGW = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:)) + sum(Ec(:))
! Dump results
write(*,*)'-------------------------------------------------------------------------------&

View File

@ -71,7 +71,13 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
integer :: nS_aa,nS_bb,nS_sc
double precision :: dipole(ncart)
double precision :: ET(nspin)
double precision :: EV(nspin)
double precision :: EJ(nsp)
double precision :: Ex(nspin)
double precision :: Ec(nsp)
double precision :: EcRPA
double precision :: EqsGW
double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
double precision :: Conv
@ -296,10 +302,50 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
P(:,:,is) = matmul(c(:,1:nO(is),is),transpose(c(:,1:nO(is),is)))
end do
!------------------------------------------------------------------------
! Compute total energy
!------------------------------------------------------------------------
! Kinetic energy
do ispin=1,nspin
ET(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),T(:,:)))
end do
! Potential energy
do ispin=1,nspin
EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:)))
end do
! Coulomb energy
EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1)))
EJ(2) = trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2)))
EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2)))
! Exchange energy
do ispin=1,nspin
Ex(ispin) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,ispin),K(:,:,ispin)))
end do
! Correlation energy
Ec(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),SigC(:,:,1)))
Ec(2) = trace_matrix(nBas,matmul(P(:,:,1),SigC(:,:,2)))
Ec(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),SigC(:,:,2)))
! Total energy
EqsGW = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:)) + sum(Ec(:))
!------------------------------------------------------------------------
! Print results
!------------------------------------------------------------------------
call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
call print_qsUGW(nBas,nO,S,nSCF,Conv,thresh,eGW,c,P,T,V,J,K,ENuc,EHF,SigC,Z,EcRPA,dipole)
call print_qsUGW(nBas,nO,nSCF,Conv,thresh,eGW,c,P,S,T,V,J,K,ENuc,ET,EV,EJ,Ex,Ec,EcRPA,EqsGW,SigC,Z,dipole)
enddo
!------------------------------------------------------------------------