From 91a0120eee9469cf768a6a06fc4bab16da4e72a6 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 8 Jul 2020 10:55:03 +0200 Subject: [PATCH] fix problem with individual energies --- input/dft | 20 +++++++++--------- input/int | 6 +++--- input/options | 2 +- src/QuAcK/G0W0.f90 | 4 ++-- .../UCC_lda_exchange_individual_energy.f90 | 11 ++++++++-- ...VWN5_lda_correlation_individual_energy.f90 | 4 ++-- src/eDFT/print_UKS.f90 | 21 +++++++++++++++---- .../print_unrestricted_individual_energy.f90 | 2 +- 8 files changed, 45 insertions(+), 25 deletions(-) diff --git a/input/dft b/input/dft index 9d15b73..eefb2f9 100644 --- a/input/dft +++ b/input/dft @@ -1,25 +1,25 @@ # Restricted or unrestricted KS calculation eDFT-UKS # exchange rung: -# Hartree = 0 -# LDA = 1: RS51,RMFL20 -# GGA = 2: RB88 -# Hybrid = 4 -# Hartree-Fock = 666 - 666 HF +# Hartree = 0: H +# LDA = 1: S51,CC +# GGA = 2: B88 +# Hybrid = 4: +# Hartree-Fock = 666: HF + 1 S51 # correlation rung: -# Hartree = 0 -# LDA = 1: RVWN5,RMFL20 +# Hartree = 0: H +# LDA = 1: VWN5,MFL20 # GGA = 2: # Hybrid = 4: -# Hartree-Fock = 666 +# Hartree-Fock = 666: HF 0 H # quadrature grid SG-n 1 # Number of states in ensemble (nEns) 3 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 0.0 0.0 + 0.000000 0.0 # Parameters for CC weight-dependent exchange functional 0.000000 0.0000000 0.000000 0.000000 0.0000000 0.0000000 diff --git a/input/int b/input/int index 9341cab..86aa017 100644 --- a/input/int +++ b/input/int @@ -1,13 +1,13 @@ # Debuggin mode? F # Chemist notation for two-electron integral? - F + T # Exposant of the Slater geminal - 0.5 + 1.0 # One-electron integrals: Ov Kin Nuc T T T # Two-electron integrals: ERI F12 Yuk Erf - T F F T + T F F F # Three-electron integrals: Type1 Type2 Type3 F F F # Four-electron integrals: Type1 Type2 Type3 diff --git a/input/options b/input/options index a8ee178..348d4b6 100644 --- a/input/options +++ b/input/options @@ -9,7 +9,7 @@ # GF: maxSCF thresh DIIS n_diis lin eta renorm 256 0.00001 T 5 T 0.00367493 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 T F F + 256 0.00001 T 5 T 0.00367493 F F F F F # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index 96b12dc..bd76850 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -104,6 +104,8 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + if(print_W) call print_excitation('RPA@HF ',ispin,nS,Omega(:,ispin)) + ! Compute correlation part of the self-energy call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) @@ -143,8 +145,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, ! Dump results - if(print_W) call print_excitation('RPA@G0W0 ',ispin,nS,Omega(:,ispin)) - call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM) ! Compute the RPA correlation energy diff --git a/src/eDFT/UCC_lda_exchange_individual_energy.f90 b/src/eDFT/UCC_lda_exchange_individual_energy.f90 index 26e08de..76560f6 100644 --- a/src/eDFT/UCC_lda_exchange_individual_energy.f90 +++ b/src/eDFT/UCC_lda_exchange_individual_energy.f90 @@ -86,11 +86,18 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig r = max(0d0,rhow(iG)) rI = max(0d0,rho(iG)) - if(r > threshold .and. rI > threshold) then + if(r > threshold) then e_p = Cx*r**(1d0/3d0) dedr = 1d0/3d0*Cx*r**(-2d0/3d0) - Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI - dedr*r*r) + + Ex = Ex - weight(iG)*dedr*r*r + + if(rI > threshold) then + + Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI) + + endif endif diff --git a/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 b/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 index 7341477..80c6990 100644 --- a/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 +++ b/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 @@ -61,7 +61,7 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) ! spin-up contribution - if(ra > threshold .and. raI > threshold) then + if(ra > threshold .or. raI > threshold) then r = ra + rb rI = raI + rbI @@ -131,7 +131,7 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) ! spin-down contribution - if(rb > threshold .and. rbI > threshold) then + if(rb > threshold .or. rbI > threshold) then r = ra + rb rI = raI + rbI diff --git a/src/eDFT/print_UKS.f90 b/src/eDFT/print_UKS.f90 index 80be911..2c0865e 100644 --- a/src/eDFT/print_UKS.f90 +++ b/src/eDFT/print_UKS.f90 @@ -17,18 +17,31 @@ subroutine print_UKS(nBas,nO,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew) double precision,intent(in) :: Ec(nsp) double precision,intent(in) :: Ew + integer :: ispin integer :: HOMO(nspin) integer :: LUMO(nspin) double precision :: Gap(nspin) ! HOMO and LUMO - HOMO(:) = nO(:) - LUMO(:) = HOMO(:) + 1 + do ispin=1,nspin - Gap(1) = eps(LUMO(1),1) - eps(HOMO(1),1) - Gap(2) = eps(LUMO(2),2) - eps(HOMO(2),2) + if(nO(ispin) > 0) then + + HOMO(ispin) = nO(ispin) + LUMO(ispin) = HOMO(ispin) + 1 + Gap(ispin) = eps(LUMO(ispin),ispin) - eps(HOMO(ispin),ispin) + + else + + HOMO(ispin) = 0 + LUMO(ispin) = 0 + Gap(ispin) = 0d0 + + end if + + end do ! Dump results diff --git a/src/eDFT/print_unrestricted_individual_energy.f90 b/src/eDFT/print_unrestricted_individual_energy.f90 index c2abf33..fa4f3cc 100644 --- a/src/eDFT/print_unrestricted_individual_energy.f90 +++ b/src/eDFT/print_unrestricted_individual_energy.f90 @@ -137,7 +137,7 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc, !------------------------------------------------------------------------ write(*,'(A60)') '-------------------------------------------------' - write(*,'(A60)') ' IP and EA FROM AUXILIARY ENERGIES ' + write(*,'(A60)') ' IP AND EA FROM AUXILIARY ENERGIES ' write(*,'(A60)') '-------------------------------------------------' write(*,'(A43,F16.10,A4)') ' Ionization Potential 1 -> 2:',Omaux(2)+OmxcDD(2),' au'