LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-06-01 19:05:27 +02:00
Code Releases Activity

debug qsGF2/qsUGF2

Browse Source
This commit is contained in:
Pierre-Francois Loos 2021-03-23 19:50:46 +01:00
parent b5a28cca0d
commit d01e7c248c
9 changed files with 42 additions and 42 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
input/dft
View File

@ -1,12 +1,12 @@
# Restricted or unrestricted KS calculation
UKS
eDFT-UKS
# exchange rung:
# Hartree = 0: H
# LDA = 1: S51,CC-S51
# GGA = 2: B88,G96,PBE
# MGGA = 3:
# Hybrid = 4: HF,B3,PBE
1 S51
4 HF
# correlation rung:
# Hartree = 0: H
# LDA = 1: PW92,VWN3,VWN5,eVWN5
@ -17,18 +17,18 @@
# quadrature grid SG-n
1
# Number of states in ensemble (nEns)
1
2
# occupation numbers
1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
# Ensemble weights: wEns(1),...,wEns(nEns-1)
0.0 0.0
0.5 0.0
# N-centered?
F
# Parameters for CC weight-dependent exchange functional

2
input/methods
View File

@ -11,7 +11,7 @@
# RPA* RPAx* ppRPA
F F F
# G0F2* evGF2* qsGF2* G0F3 evGF3
F T F F F
F F T F F
# G0W0* evGW* qsGW*
F F F
# G0T0 evGT qsGT

8
src/GF/print_qsGF2.f90
View File

@ -69,7 +69,7 @@ subroutine print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,P,T,V,J,K,F,SigC,Z, &
write(*,'(2X,A30,F15.6,A3)') 'qsGF2 LUMO energy:',eGF2(LUMO)*HaToeV,' eV'
write(*,'(2X,A30,F15.6,A3)') 'qsGF2 HOMO-LUMO gap :',Gap*HaToeV,' eV'
write(*,*)'-------------------------------------------'
write(*,'(2X,A30,F15.6,A3)') ' qsGF2 total energy:',ENuc + EqsGF2 + Ec,' au'
write(*,'(2X,A30,F15.6,A3)') ' qsGF2 total energy:',ENuc + EqsGF2,' au'
write(*,'(2X,A30,F15.6,A3)') ' qsGF2 exchange energy:',Ex,' au'
write(*,'(2X,A30,F15.6,A3)') ' qsGF2 correlation energy:',Ec,' au'
write(*,*)'-------------------------------------------'
@ -87,14 +87,14 @@ subroutine print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,P,T,V,J,K,F,SigC,Z, &
write(*,'(A32,1X,F16.10,A3)') ' Kinetic energy: ',ET,' au'
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)') ' Two-electron energy: ',EJ + Ex + Ec,' 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)') '---------------------------------------'
write(*,'(A32,1X,F16.10,A3)') ' Electronic energy: ',EqsGF2 + Ec,' au'
write(*,'(A32,1X,F16.10,A3)') ' Electronic energy: ',EqsGF2,' au'
write(*,'(A32,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au'
write(*,'(A32,1X,F16.10,A3)') ' qsGF2 energy: ',ENuc + EqsGF2 + Ec,' au'
write(*,'(A32,1X,F16.10,A3)') ' qsGF2 energy: ',ENuc + EqsGF2,' au'
write(*,'(A50)') '---------------------------------------'
write(*,'(A35)') ' Dipole moment (Debye) '
write(*,'(10X,4A10)') 'X','Y','Z','Tot.'

10
src/GF/print_qsUGF2.f90
View File

@ -99,7 +99,7 @@ subroutine print_qsUGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,cGF2,PGF2,Ov,T,V,J,K,
write(*,'(2X,A30,F15.6,A3)') 'qsUGF2 HOMO-LUMO gap :',(minval(LUMO(:))-maxval(HOMO(:)))*HaToeV,' eV'
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,'(2X,A30,F15.6,A3)') ' qsUGF2 total energy:',ENuc + EqsGF2 + sum(Ec(:)),' au'
write(*,'(2X,A30,F15.6,A3)') ' qsUGF2 total energy:',ENuc + EqsGF2,' au'
write(*,'(2X,A30,F15.6,A3)') ' qsUGF2 exchange energy:',sum(Ex(:)),' au'
write(*,'(2X,A30,F15.6,A3)') ' qsUGF2 correlation energy:',sum(Ec(:)),' au'
write(*,*)'-------------------------------------------------------------------------------&
@ -126,10 +126,10 @@ subroutine print_qsUGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,cGF2,PGF2,Ov,T,V,J,K,
write(*,'(A40,1X,F16.10,A3)') ' Potential a energy: ',EV(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Potential b energy: ',EV(2),' au'
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,1X,F16.10,A3)') ' Two-electron energy: ',sum(EJ(:)) + sum(Ex(:)),' au'
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)') ' Two-electron energy: ',sum(EJ(:)) + sum(Ex(:)) + sum(Ec(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Two-electron aa energy: ',EJ(1) + Ex(1) + Ec(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Two-electron ab energy: ',EJ(2) + Ec(2),' au'
write(*,'(A40,1X,F16.10,A3)') ' Two-electron bb energy: ',EJ(3) + Ex(2) + Ec(3),' au'
write(*,*)
write(*,'(A40,1X,F16.10,A3)') ' Hartree energy: ',sum(EJ(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Hartree aa energy: ',EJ(1),' au'

15
src/GF/qsGF2.f90
View File

@ -197,16 +197,17 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,
! Exchange energy
Ex = -0.25d0*trace_matrix(nBas,matmul(P,K))
! Total energy
EqsGF2 = ET + EV + EJ + Ex
Ex = 0.25d0*trace_matrix(nBas,matmul(P,K))
! Correlation energy
call MP2(nBas,nC,nO,nV,nR,ERI_MO,ENuc,EqsGF2,eGF2,Ec)
! Total energy
EqsGF2 = ET + EV + EJ + Ex + Ec
!------------------------------------------------------------------------
! Print results
!------------------------------------------------------------------------
@ -248,8 +249,8 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGF2 correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGF2 correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGF2 correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGF2 total energy =',ENuc + EqsGF2 + EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGF2 correlation energy =',sum(EcBSE(:))
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGF2 total energy =',ENuc + EqsGF2 + sum(EcBSE(:))
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)

10
src/GF/qsUGF2.f90
View File

@ -64,7 +64,7 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,
double precision :: EV(nspin)
double precision :: EJ(nsp)
double precision :: Ex(nspin)
double precision :: Ec(nspin)
double precision :: Ec(nsp)
double precision :: EqsGF2
double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
@ -268,14 +268,14 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,
Ex(is) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,is),K(:,:,is)))
end do
! Total energy
EqsGF2 = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:))
! Correlation energy
call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EqsGF2,eGF2,Ec)
! Total energy
EqsGF2 = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:)) + sum(Ec(:))
!------------------------------------------------------------------------
! Print results
!------------------------------------------------------------------------

2
src/GF/self_energy_GF2.f90
View File

@ -28,7 +28,7 @@ subroutine self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z)
! Initialize
SigC(:,:) = 0d0
Z(:) = 0d0
Z(:) = 0d0
! Compute GF2 self-energy and renormalization factor

1
src/GF/unrestricted_self_energy_GF2.f90
View File

@ -36,6 +36,7 @@ subroutine unrestricted_self_energy_GF2(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_b
!---------------------!
SigC(:,:,:) = 0d0
Z(:,:) = 0d0
!----------------!
! Spin-up sector

16
src/MP/MP2.f90
View File

@ -23,7 +23,7 @@ subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
! Output variables
double precision,intent(out) :: EcMP2(3)
double precision,intent(out) :: EcMP2
! Hello world
@ -57,20 +57,18 @@ subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
enddo
enddo
EcMP2(2) = 2d0*E2a
EcMP2(3) = -E2b
EcMP2(1) = EcMP2(2) + EcMP2(3)
EcMP2 = 2d0*E2a - E2b
write(*,*)
write(*,'(A32)') '--------------------------'
write(*,'(A32)') ' MP2 calculation '
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',EcMP2(1)
write(*,'(A32,1X,F16.10)') ' Direct part = ',EcMP2(2)
write(*,'(A32,1X,F16.10)') ' Exchange part = ',EcMP2(3)
write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',EcMP2
write(*,'(A32,1X,F16.10)') ' Direct part = ',2d0*E2a
write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2b
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ',EHF + EcMP2(1)
write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + EcMP2(1)
write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ',EHF + EcMP2
write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + EcMP2
write(*,'(A32)') '--------------------------'
write(*,*)