PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-03 20:53:53 +01:00
Code Releases Activity

debug GM correlation energy

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-11-04 14:47:45 +01:00
parent 6061f1b507
commit f44ab1eadb
12 changed files with 61 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
input/methods
View File

@ -13,7 +13,7 @@
# G0F2 evGF2 G0F3 evGF3
F F F F
# G0W0* evGW* qsGW*
T F F
F F T
# G0T0 evGT qsGT
F F F
# MCMP2

8
input/options
View File

@ -5,14 +5,14 @@
# CC: maxSCF thresh DIIS n_diis
64 0.0000001 T 5
# spin: TDA singlet triplet spin_conserved spin_flip
F T T T T
T T T T T
# 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
256 0.00001 T 5 T 0.00367493 F F F F F
# GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0
256 0.00001 T 5 T 0.0 F F F F F
# ACFDT: AC Kx XBS
F F T
# BSE: BSE dBSE dTDA evDyn
F F T F
F T T F
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T

2
src/HF/print_RHF.f90
View File

@ -44,7 +44,7 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole)
write(*,'(A32,1X,F16.10,A3)') ' Potential energy: ',EV,' au'
write(*,'(A50)') '-----------------------------------------'
write(*,'(A32,1X,F16.10,A3)') ' Two-electron energy: ',EJ + EK,' au'
write(*,'(A32,1X,F16.10,A3)') ' Coulomb energy: ',EJ,' au'
write(*,'(A32,1X,F16.10,A3)') ' Hartree energy: ',EJ,' au'
write(*,'(A32,1X,F16.10,A3)') ' Exchange energy: ',EK,' au'
write(*,'(A50)') '-----------------------------------------'
write(*,'(A32,1X,F16.10,A3)') ' Electronic energy: ',ERHF,' au'

8
src/HF/print_UHF.f90
View File

@ -70,10 +70,10 @@ subroutine print_UHF(nBas,nO,Ov,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
write(*,'(A40,1X,F16.10,A3)') ' Two-electron aa energy: ',EJ(1) + Ex(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Two-electron ab energy: ',EJ(2),' au'
write(*,'(A40,1X,F16.10,A3)') ' Two-electron bb energy: ',EJ(3) + Ex(2),' au'
write(*,'(A40,1X,F16.10,A3)') ' Coulomb energy: ',sum(EJ(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Coulomb aa energy: ',EJ(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Coulomb ab energy: ',EJ(2),' au'
write(*,'(A40,1X,F16.10,A3)') ' Coulomb bb energy: ',EJ(3),' au'
write(*,'(A40,1X,F16.10,A3)') ' Hartree energy: ',sum(EJ(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Hartree aa energy: ',EJ(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Hartree ab energy: ',EJ(2),' au'
write(*,'(A40,1X,F16.10,A3)') ' Hartree bb energy: ',EJ(3),' au'
write(*,'(A40,1X,F16.10,A3)') ' Exchange energy: ',sum(Ex(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Exchange a energy: ',Ex(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Exchange b energy: ',Ex(2),' au'

5
src/MBPT/UG0W0.f90
View File

@ -49,6 +49,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
integer :: is
integer :: ispin
double precision :: EcRPA
double precision :: EcGM(nspin)
double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
double precision,allocatable :: SigC(:,:)
@ -131,7 +132,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
! Compute self-energy !
!---------------------!
call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC)
call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,EcGM)
!--------------------------------!
! Compute renormalization factor !
@ -170,7 +171,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
! Dump results
call print_UG0W0(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA)
call print_UG0W0(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA,EcGM)
! Free memory

7
src/MBPT/evUGW.f90
View File

@ -57,6 +57,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE
double precision :: rcond(nspin)
double precision :: Conv
double precision :: EcRPA
double precision :: EcGM(nspin)
double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
double precision :: alpha
@ -167,12 +168,12 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE
if(G0W) then
call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC)
call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,EcGM)
call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z)
else
call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC)
call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,EcGM)
call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,Z)
endif
@ -189,7 +190,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE
! Print results
call print_evUGW(nBas,nO,nSCF,Conv,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA)
call print_evUGW(nBas,nO,nSCF,Conv,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA,EcGM)
! Linear mixing or DIIS extrapolation

5
src/MBPT/print_UG0W0.f90
View File

@ -1,4 +1,4 @@
subroutine print_UG0W0(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA)
subroutine print_UG0W0(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA,EcGM)
! Print one-electron energies and other stuff for G0W0
@ -10,6 +10,7 @@ subroutine print_UG0W0(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA)
double precision,intent(in) :: ENuc
double precision,intent(in) :: EUHF
double precision,intent(in) :: EcRPA
double precision,intent(in) :: EcGM(nspin)
double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: SigC(nBas,nspin)
double precision,intent(in) :: Z(nBas,nspin)
@ -64,6 +65,8 @@ subroutine print_UG0W0(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA)
-------------------------------------------------'
write(*,'(2X,A30,F15.6,A3)') 'RPA@UG0W0 total energy :',ENuc + EUHF + EcRPA,' au'
write(*,'(2X,A30,F15.6,A3)') 'RPA@UG0W0 correlation energy:',EcRPA,' au'
write(*,'(2X,A30,F15.6,A3)') ' GM@UG0W0 total energy :',ENuc + EUHF + sum(EcGM(:)),' au'
write(*,'(2X,A30,F15.6,A3)') ' GM@UG0W0 correlation energy:',sum(EcGM(:)),' au'
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,*)

5
src/MBPT/print_evUGW.f90
View File

@ -1,4 +1,4 @@
subroutine print_evUGW(nBas,nO,nSCF,Conv,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA)
subroutine print_evUGW(nBas,nO,nSCF,Conv,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA,EcGM)
! Print one-electron energies and other stuff for evGW
@ -11,6 +11,7 @@ subroutine print_evUGW(nBas,nO,nSCF,Conv,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA)
double precision,intent(in) :: ENuc
double precision,intent(in) :: EUHF
double precision,intent(in) :: EcRPA
double precision,intent(in) :: EcGM(nspin)
double precision,intent(in) :: Conv
double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: SigC(nBas,nspin)
@ -74,6 +75,8 @@ subroutine print_evUGW(nBas,nO,nSCF,Conv,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA)
-------------------------------------------------'
write(*,'(2X,A30,F15.6,A3)') 'RPA@evGW total energy :',ENuc + EUHF + EcRPA,' au'
write(*,'(2X,A30,F15.6,A3)') 'RPA@evGW correlation energy:',EcRPA,' au'
write(*,'(2X,A30,F15.6,A3)') ' GM@evGW total energy :',ENuc + EUHF + sum(EcGM(:)),' au'
write(*,'(2X,A30,F15.6,A3)') ' GM@evGW correlation energy:',sum(EcGM(:)),' au'
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,*)

4
src/MBPT/print_qsGW.f90
View File

@ -37,7 +37,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z
EV = trace_matrix(nBas,matmul(P,V))
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
Ex = 0.25d0*trace_matrix(nBas,matmul(P,K))
Ec = 0.50d0*trace_matrix(nBas,matmul(P,SigC))
Ec = -0.50d0*trace_matrix(nBas,matmul(P,SigC))
EqsGW = ET + EV + EJ + Ex + Ec
! Dump results
@ -86,7 +86,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z
write(*,'(A32,1X,F16.10,A3)') ' Potential energy: ',EV,' au'
write(*,'(A50)') '---------------------------------------'
write(*,'(A32,1X,F16.10,A3)') ' Two-electron energy: ',EJ + Ex,' au'
write(*,'(A32,1X,F16.10,A3)') ' Coulomb energy: ',EJ,' au'
write(*,'(A32,1X,F16.10,A3)') ' Hartree energy: ',EJ,' au'
write(*,'(A32,1X,F16.10,A3)') ' Exchange energy: ',Ex,' au'
write(*,'(A32,1X,F16.10,A3)') ' Correlation energy: ',Ec,' au'
write(*,'(A50)') '---------------------------------------'

8
src/MBPT/print_qsUGW.f90
View File

@ -133,10 +133,10 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, &
write(*,'(A40,1X,F16.10,A3)') ' Two-electron ab energy: ',EJ(2),' au'
write(*,'(A40,1X,F16.10,A3)') ' Two-electron bb energy: ',EJ(3) + Ex(2),' au'
write(*,*)
write(*,'(A40,1X,F16.10,A3)') ' Coulomb energy: ',sum(EJ(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Coulomb aa energy: ',EJ(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Coulomb ab energy: ',EJ(2),' au'
write(*,'(A40,1X,F16.10,A3)') ' Coulomb bb energy: ',EJ(3),' au'
write(*,'(A40,1X,F16.10,A3)') ' Hartree energy: ',sum(EJ(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Hartree aa energy: ',EJ(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Hartree ab energy: ',EJ(2),' au'
write(*,'(A40,1X,F16.10,A3)') ' Hartree bb energy: ',EJ(3),' au'
write(*,*)
write(*,'(A40,1X,F16.10,A3)') ' Exchange energy: ',sum(Ex(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Exchange a energy: ',Ex(1),' au'

8
src/MBPT/qsUGW.f90
View File

@ -332,10 +332,10 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
! Correlation energy
Ec(1) = 0.25d0*trace_matrix(nBas,matmul(P(:,:,1),SigCp(:,:,1)))
Ec(2) = 0.25d0*trace_matrix(nBas,matmul(P(:,:,1),SigCp(:,:,2))) &
+ 0.25d0*trace_matrix(nBas,matmul(P(:,:,2),SigCp(:,:,1)))
Ec(3) = 0.25d0*trace_matrix(nBas,matmul(P(:,:,2),SigCp(:,:,2)))
Ec(1) = - 0.25d0*trace_matrix(nBas,matmul(P(:,:,1),SigCp(:,:,1)))
Ec(2) = - 0.25d0*trace_matrix(nBas,matmul(P(:,:,1),SigCp(:,:,2))) &
- 0.25d0*trace_matrix(nBas,matmul(P(:,:,2),SigCp(:,:,1)))
Ec(3) = - 0.25d0*trace_matrix(nBas,matmul(P(:,:,2),SigCp(:,:,2)))
! Total energy

27
src/MBPT/unrestricted_self_energy_correlation_diag.f90
View File

@ -1,4 +1,4 @@
subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC)
subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC,EcGM)
! Compute diagonal of the correlation part of the self-energy
@ -26,10 +26,12 @@ subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSt,e,
! Output variables
double precision,intent(out) :: SigC(nBas,nspin)
double precision :: EcGM(nspin)
! Initialize
SigC(:,:) = 0d0
EcGM(:) = 0d0
!--------------!
! Spin-up part !
@ -57,6 +59,17 @@ subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSt,e,
end do
end do
! GM correlation energy
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
do jb=1,nSt
eps = e(a,1) - e(i,1) + Omega(jb)
EcGM(1) = EcGM(1) - rho(a,i,jb,1)**2*eps/(eps**2 + eta**2)
end do
end do
end do
!----------------!
! Spin-down part !
!----------------!
@ -83,4 +96,16 @@ subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSt,e,
end do
end do
! GM correlation energy
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
do jb=1,nSt
eps = e(a,2) - e(i,2) + Omega(jb)
EcGM(2) = EcGM(2) - rho(a,i,jb,2)**2*eps/(eps**2 + eta**2)
end do
end do
end do
end subroutine unrestricted_self_energy_correlation_diag