diff --git a/input/dft b/input/dft index 9537434..a3810e3 100644 --- a/input/dft +++ b/input/dft @@ -1,14 +1,24 @@ # Restricted or unrestricted KS calculation GOK-RKS -# exchange rung: Hartree = 0, LDA = 1 (RS51,S51), GGA = 2(G96,B88), hybrid = 4, Hartree-Fock = 666 - 1 RS51 -# correlation rung: Hartree = 0, LDA = 1 (W38,VWN5,C16,LF19), GGA = 2(LYP), hybrid = 4(B3LYP), Hartree-Fock = 666 +# exchange rung: +# Hartree = 0 +# LDA = 1: RS51,S51,RMFL20 +# GGA = 2: G96,B88 +# Hybrid = 4 +# Hartree-Fock = 666 + 1 RMFL20 +# correlation rung: +# Hartree = 0 +# LDA = 1: W38,VWN5,C16,RMFL20 +# GGA = 2: LYP +# Hybrid = 4: B3LYP +# Hartree-Fock = 666 1 RVWN5 # quadrature grid SG-n 1 # Number of states in ensemble (nEns) 2 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 0.00000 0.00000 -# eKS: maxSCF thresh DIIS n_diis guess_type ortho_type - 64 0.0000001 T 5 1 1 + 0.50000 0.00000 +# GOK-DFT: maxSCF thresh DIIS n_diis guess_type ortho_type + 64 0.0000001 T 5 1 1 diff --git a/src/eDFT/RMFL20_lda_exchange_energy.f90 b/src/eDFT/RMFL20_lda_exchange_energy.f90 index 2fd6a33..84ed417 100644 --- a/src/eDFT/RMFL20_lda_exchange_energy.f90 +++ b/src/eDFT/RMFL20_lda_exchange_energy.f90 @@ -29,11 +29,11 @@ subroutine RMFL20_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) ! Cx coefficient for Slater LDA exchange - Cx0 = -(4d0/3d0)*(1d0/pi)**(1d0/3d0) - Cx1 = -(176d0/105d0)*(1d0/pi)**(1d0/3d0) - CxLDA = -(3d0/4d0)*(3d0/pi)**(1d0/3d0) + Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0) + Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0) + CxLDA = - (3d0/4d0)*(3d0/pi)**(1d0/3d0) - Cxw = CxLDA + wEns(1)*(Cx1 - Cx0) + Cxw = CxLDA + wEns(2)*(Cx1 - Cx0) ! Compute LDA exchange energy diff --git a/src/eDFT/RMFL20_lda_exchange_potential.f90 b/src/eDFT/RMFL20_lda_exchange_potential.f90 index eebb5d4..dee66f1 100644 --- a/src/eDFT/RMFL20_lda_exchange_potential.f90 +++ b/src/eDFT/RMFL20_lda_exchange_potential.f90 @@ -28,13 +28,13 @@ subroutine RMFL20_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) double precision,intent(out) :: Fx(nBas,nBas) -! Cx coefficient for Slater LDA exchange +! Weight-dependent Cx coefficient for RMFL20 exchange functional Cx0 = -(4d0/3d0)*(1d0/pi)**(1d0/3d0) Cx1 = -(176d0/105d0)*(1d0/pi)**(1d0/3d0) CxLDA = -(3d0/4d0)*(3d0/pi)**(1d0/3d0) - Cxw = CxLDA + wEns(1)*(Cx1 - Cx0) + Cxw = CxLDA + wEns(2)*(Cx1 - Cx0) ! Compute LDA exchange matrix in the AO basis diff --git a/src/eDFT/RS51_lda_exchange_energy.f90 b/src/eDFT/RS51_lda_exchange_energy.f90 index 1c6c965..07c37c0 100644 --- a/src/eDFT/RS51_lda_exchange_energy.f90 +++ b/src/eDFT/RS51_lda_exchange_energy.f90 @@ -34,7 +34,7 @@ subroutine RS51_lda_exchange_energy(nGrid,weight,rho,Ex) if(r > threshold) then - Ex = Ex + weight(iG)*Cx*r**(4d0/3d0) + Ex = Ex + weight(iG)*Cx*r**(4d0/3d0) endif diff --git a/src/eDFT/lda_exchange_derivative_discontinuity.f90 b/src/eDFT/lda_exchange_derivative_discontinuity.f90 index c607139..cb4527c 100644 --- a/src/eDFT/lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/lda_exchange_derivative_discontinuity.f90 @@ -35,8 +35,7 @@ subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow case ('RMFL20') -! call MFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),ExDD(:)) - ExDD(:) = 0d0 + call RMFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),ExDD(:)) case default diff --git a/src/eDFT/lda_exchange_individual_energy.f90 b/src/eDFT/lda_exchange_individual_energy.f90 index 6fbb190..8b80244 100644 --- a/src/eDFT/lda_exchange_individual_energy.f90 +++ b/src/eDFT/lda_exchange_individual_energy.f90 @@ -29,7 +29,7 @@ subroutine lda_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,rho,Ex case ('RMFL20') -! call RMFL20_lda_exchange_individual_energy(nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ex) + call RMFL20_lda_exchange_individual_energy(nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ex) case default diff --git a/src/eDFT/print_restricted_individual_energy.f90 b/src/eDFT/print_restricted_individual_energy.f90 index 45b8a3e..671b436 100644 --- a/src/eDFT/print_restricted_individual_energy.f90 +++ b/src/eDFT/print_restricted_individual_energy.f90 @@ -1,4 +1,5 @@ -subroutine print_restricted_individual_energy(nEns,ET,EV,EJ,Ex,Ec,Exc,ExDD,EcDD,ExcDD,E,Om) +subroutine print_restricted_individual_energy(nEns,ET,EV,EJ,Ex,Ec,Exc,ExDD,EcDD,ExcDD,E, & + Om,Omx,Omc,Omxc,OmxDD,OmcDD,OmxcDD) ! Print individual energies for eDFT calculation @@ -11,8 +12,10 @@ subroutine print_restricted_individual_energy(nEns,ET,EV,EJ,Ex,Ec,Exc,ExDD,EcDD, double precision,intent(in) :: ET(nEns) double precision,intent(in) :: EV(nEns) double precision,intent(in) :: EJ(nEns) - double precision,intent(in) :: Ex(nEns),Ec(nEns),Exc(nEns) - double precision,intent(in) :: ExDD(nEns),EcDD(nEns),ExcDD(nEns) + double precision,intent(in) :: Ex(nEns), Ec(nEns), Exc(nEns) + double precision,intent(in) :: ExDD(nEns), EcDD(nEns), ExcDD(nEns) + double precision,intent(in) :: Omx(nEns), Omc(nEns), Omxc(nEns) + double precision,intent(in) :: OmxDD(nEns),OmcDD(nEns),OmxcDD(nEns) double precision,intent(in) :: E(nEns) double precision,intent(in) :: Om(nEns) @@ -25,7 +28,7 @@ subroutine print_restricted_individual_energy(nEns,ET,EV,EJ,Ex,Ec,Exc,ExDD,EcDD, !------------------------------------------------------------------------ write(*,'(A60)') '-------------------------------------------------' - write(*,'(A50)') ' Individual Kinetic energies' + write(*,'(A50)') ' Individual kinetic energies' write(*,'(A60)') '-------------------------------------------------' do iEns=1,nEns write(*,'(A40,I2,A2,F16.10,A3)') ' Kinetic energy state ',iEns,': ',ET(iEns),' au' @@ -38,7 +41,7 @@ subroutine print_restricted_individual_energy(nEns,ET,EV,EJ,Ex,Ec,Exc,ExDD,EcDD, !------------------------------------------------------------------------ write(*,'(A60)') '-------------------------------------------------' - write(*,'(A50)') ' Individual Potential energies' + write(*,'(A50)') ' Individual potential energies' write(*,'(A60)') '-------------------------------------------------' do iEns=1,nEns write(*,'(A40,I2,A2,F16.10,A3)') ' Potential energy state ',iEns,': ',EV(iEns),' au' @@ -93,12 +96,12 @@ subroutine print_restricted_individual_energy(nEns,ET,EV,EJ,Ex,Ec,Exc,ExDD,EcDD, write(*,'(A50)') ' Derivative discontinuities (DD) ' write(*,'(A60)') '-------------------------------------------------' do iEns=1,nEns - write(*,'(A40,I2,A2,F16.10,A3)') ' x ensemble derivative ',iEns,': ',ExDD(iEns), ' au' - write(*,'(A40,I2,A2,F16.10,A3)') ' c ensemble derivative ',iEns,': ',EcDD(iEns), ' au' - write(*,'(A40,I2,A2,F16.10,A3)') ' xc ensemble derivative ',iEns,': ',ExcDD(iEns),' au' + write(*,'(A40,I2,A2,F16.10,A3)') ' x ensemble derivative state ',iEns,': ',ExDD(iEns), ' au' + write(*,'(A40,I2,A2,F16.10,A3)') ' c ensemble derivative state ',iEns,': ',EcDD(iEns), ' au' + write(*,'(A40,I2,A2,F16.10,A3)') ' xc ensemble derivative state ',iEns,': ',ExcDD(iEns),' au' + write(*,*) end do write(*,'(A60)') '-------------------------------------------------' - write(*,*) !------------------------------------------------------------------------ ! Total and Excitation energies @@ -111,12 +114,30 @@ subroutine print_restricted_individual_energy(nEns,ET,EV,EJ,Ex,Ec,Exc,ExDD,EcDD, write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns),' au' end do write(*,'(A60)') '-------------------------------------------------' + do iEns=2,nEns - write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns),' au' + write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns), ' au' + write(*,*) + write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(iEns), ' au' + write(*,'(A44, F16.10,A3)') ' c energy contribution : ',Omc(iEns), ' au' + write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',Omxc(iEns), ' au' + write(*,*) + write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns), ' au' + write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns), ' au' + write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iEns),' au' end do write(*,'(A60)') '-------------------------------------------------' + do iEns=2,nEns - write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns)*HaToeV,' eV' + write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns)*HaToeV, ' eV' + write(*,*) + write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(iEns)*HaToeV, ' eV' + write(*,'(A44, F16.10,A3)') ' c energy contribution : ',Omc(iEns)*HaToeV, ' eV' + write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',Omxc(iEns)*HaToeV, ' eV' + write(*,*) + write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns)*HaToeV, ' eV' + write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns)*HaToeV, ' eV' + write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iEns)*HaToeV,' eV' end do write(*,'(A60)') '-------------------------------------------------' write(*,*) diff --git a/src/eDFT/read_options.f90 b/src/eDFT/read_options.f90 index 19a0a52..080485e 100644 --- a/src/eDFT/read_options.f90 +++ b/src/eDFT/read_options.f90 @@ -51,11 +51,21 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,maxSCF,th ! EXCHANGE: read rung of Jacob's ladder + read(1,*) + read(1,*) + read(1,*) + read(1,*) + read(1,*) read(1,*) read(1,*) x_rung,x_DFA ! CORRELATION: read rung of Jacob's ladder + read(1,*) + read(1,*) + read(1,*) + read(1,*) + read(1,*) read(1,*) read(1,*) c_rung,c_DFA diff --git a/src/eDFT/restricted_individual_energy.f90 b/src/eDFT/restricted_individual_energy.f90 index d67a2bd..0b15e40 100644 --- a/src/eDFT/restricted_individual_energy.f90 +++ b/src/eDFT/restricted_individual_energy.f90 @@ -42,8 +42,10 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri double precision :: ET(nEns) double precision :: EV(nEns) double precision :: EJ(nEns) - double precision :: Ex(nEns),Ec(nEns),Exc(nEns) - double precision :: ExDD(nEns),EcDD(nEns),ExcDD(nEns) + double precision :: Ex(nEns), Ec(nEns), Exc(nEns) + double precision :: ExDD(nEns), EcDD(nEns), ExcDD(nEns) + double precision :: Omx(nEns), Omc(nEns), Omxc(nEns) + double precision :: OmxDD(nEns),OmcDD(nEns),OmxcDD(nEns) double precision,external :: trace_matrix @@ -126,7 +128,17 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri !------------------------------------------------------------------------ do iEns=1,nEns - Om(iEns) = E(iEns) - E(1) + + Om(iEns) = E(iEns) - E(1) + + Omx(iEns) = Ex(iEns) - Ex(1) + Omc(iEns) = Ec(iEns) - Ec(1) + Omxc(iEns) = Exc(iEns) - Exc(1) + + OmxDD(iEns) = ExDD(iEns) - ExDD(1) + OmcDD(iEns) = EcDD(iEns) - EcDD(1) + OmxcDD(iEns) = ExcDD(iEns) - ExcDD(1) + end do !------------------------------------------------------------------------ @@ -134,6 +146,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri !------------------------------------------------------------------------ call print_restricted_individual_energy(nEns,ET(:),EV(:),EJ(:),Ex(:),Ec(:),Exc(:), & - ExDD(:),EcDD(:),ExcDD(:),E(:),Om(:)) + ExDD(:),EcDD(:),ExcDD(:),E(:), & + Om(:),Omx(:),Omc(:),Omxc(:),OmxDD(:),OmcDD(:),OmxcDD(:)) end subroutine restricted_individual_energy diff --git a/src/eDFT/restricted_lda_correlation_derivative_discontinuity.f90 b/src/eDFT/restricted_lda_correlation_derivative_discontinuity.f90 index 4b94b53..2c91968 100644 --- a/src/eDFT/restricted_lda_correlation_derivative_discontinuity.f90 +++ b/src/eDFT/restricted_lda_correlation_derivative_discontinuity.f90 @@ -30,7 +30,7 @@ subroutine restricted_lda_correlation_derivative_discontinuity(DFA,nEns,wEns,nGr Ec(:) = 0d0 - case ('VWN5') + case ('RVWN5') Ec(:) = 0d0