diff --git a/input/methods b/input/methods index 45d00ca..af2bd81 100644 --- a/input/methods +++ b/input/methods @@ -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 diff --git a/input/options b/input/options index 02bfdc4..4e7d31d 100644 --- a/input/options +++ b/input/options @@ -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 diff --git a/src/HF/print_RHF.f90 b/src/HF/print_RHF.f90 index 97c3437..9b88837 100644 --- a/src/HF/print_RHF.f90 +++ b/src/HF/print_RHF.f90 @@ -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' diff --git a/src/HF/print_UHF.f90 b/src/HF/print_UHF.f90 index 5ca6230..17834a6 100644 --- a/src/HF/print_UHF.f90 +++ b/src/HF/print_UHF.f90 @@ -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' diff --git a/src/MBPT/UG0W0.f90 b/src/MBPT/UG0W0.f90 index f252496..5f8b8b9 100644 --- a/src/MBPT/UG0W0.f90 +++ b/src/MBPT/UG0W0.f90 @@ -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 diff --git a/src/MBPT/evUGW.f90 b/src/MBPT/evUGW.f90 index a4f1529..5379763 100644 --- a/src/MBPT/evUGW.f90 +++ b/src/MBPT/evUGW.f90 @@ -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 diff --git a/src/MBPT/print_UG0W0.f90 b/src/MBPT/print_UG0W0.f90 index 1305461..4b45d70 100644 --- a/src/MBPT/print_UG0W0.f90 +++ b/src/MBPT/print_UG0W0.f90 @@ -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(*,*) diff --git a/src/MBPT/print_evUGW.f90 b/src/MBPT/print_evUGW.f90 index d52e43a..a4bf78a 100644 --- a/src/MBPT/print_evUGW.f90 +++ b/src/MBPT/print_evUGW.f90 @@ -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(*,*) diff --git a/src/MBPT/print_qsGW.f90 b/src/MBPT/print_qsGW.f90 index bbfd703..d8c6a3e 100644 --- a/src/MBPT/print_qsGW.f90 +++ b/src/MBPT/print_qsGW.f90 @@ -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)') '---------------------------------------' diff --git a/src/MBPT/print_qsUGW.f90 b/src/MBPT/print_qsUGW.f90 index d6413bd..14a6bdd 100644 --- a/src/MBPT/print_qsUGW.f90 +++ b/src/MBPT/print_qsUGW.f90 @@ -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' diff --git a/src/MBPT/qsUGW.f90 b/src/MBPT/qsUGW.f90 index 0da7d9c..a4a427c 100644 --- a/src/MBPT/qsUGW.f90 +++ b/src/MBPT/qsUGW.f90 @@ -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 diff --git a/src/MBPT/unrestricted_self_energy_correlation_diag.f90 b/src/MBPT/unrestricted_self_energy_correlation_diag.f90 index fa4084c..d683380 100644 --- a/src/MBPT/unrestricted_self_energy_correlation_diag.f90 +++ b/src/MBPT/unrestricted_self_energy_correlation_diag.f90 @@ -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