diff --git a/input/dft b/input/dft index 688fdad..20323d4 100644 --- a/input/dft +++ b/input/dft @@ -6,18 +6,18 @@ # GGA = 2: B88,G96,PBE # MGGA = 3: # Hybrid = 4 HF,B3LYP,PBE -1 CC-S51 +1 S51 # correlation rung: # Hartree = 0: H # LDA = 1: PW92,VWN3,VWN5,eVWN5 # GGA = 2: LYP,PBE # MGGA = 3: # Hybrid = 4: HF,B3LYP,PBE -0 H +1 VWN5 # quadrature grid SG-n 1 # Number of states in ensemble (nEns) -4 +2 # occupation numbers 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 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 @@ -31,7 +31,7 @@ 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 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 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 1 0.0 0.0 + 0.55 0.0 0.0 # Ncentered ? T # Parameters for CC weight-dependent exchange functional diff --git a/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 index 607643c..8da35b4 100644 --- a/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 @@ -1,4 +1,4 @@ -subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weight,rhow,Cx_choice,doNcentered,kappa,ExDD) +subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weight,rhow,Cx_choice,doNcentered,ExDD) ! Compute the unrestricted version of the curvature-corrected exchange ensemble derivative @@ -16,7 +16,6 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,wei double precision,intent(in) :: rhow(nGrid) integer,intent(in) :: Cx_choice logical,intent(in) :: doNcentered - double precision,intent(in) :: kappa(nEns) ! Local variables @@ -154,19 +153,9 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,wei do iEns=1,nEns do jEns=2,nEns -! if (doNcentered) then - -! ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - kappa(iEns)*wEns(jEns))*dExdw(jEns) - -! else - ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*dExdw(jEns) -! endif - end do end do - ! if(doNcentered) ExDD(:) = kappa(:)*ExDD(:) - end subroutine UCC_lda_exchange_derivative_discontinuity diff --git a/src/eDFT/UCC_lda_exchange_individual_energy.f90 b/src/eDFT/UCC_lda_exchange_individual_energy.f90 index 4fc1b1c..06835df 100644 --- a/src/eDFT/UCC_lda_exchange_individual_energy.f90 +++ b/src/eDFT/UCC_lda_exchange_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rhow,rho,Cx_choice,doNcentered,kappa,LZx,Ex) +subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rhow,rho,Cx_choice,doNcentered,LZx,Ex) ! Compute the unrestricted version of the curvature-corrected exchange functional @@ -18,7 +18,6 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho double precision,intent(in) :: rho(nGrid,nspin,nEns) integer,intent(in) :: Cx_choice logical,intent(in) :: doNcentered - double precision,intent(in) :: kappa(nEns) ! Local variables @@ -32,7 +31,8 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho ! Output variables - double precision,intent(out) :: LZx(nspin,nEns), Ex(nspin,nEns) + double precision,intent(out) :: LZx(nspin) + double precision,intent(out) :: Ex(nspin,nEns) ! Defining enhancements factor for weight-dependent functionals @@ -52,7 +52,6 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho c2 = aCC(3,2) d2 = aCC(4,2) - w1 = wEns(2) Fx1 = 1d0 + a1*w1 + b1*w1**2 + c1*w1**3 + d1*w1**4 @@ -100,7 +99,7 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho ! Compute LDA exchange matrix in the AO basis Ex(:,:) = 0d0 - LZx(:,:) = 0d0 + LZx(:) = 0d0 do ispin=1,nspin @@ -108,43 +107,20 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho r = max(0d0,rhow(iG,ispin)) - if(doNcentered) then + if(r > threshold) then - if(r > threshold) then + e = Cx*r**(+1d0/3d0) + dedr = 1d0/3d0*Cx*r**(-2d0/3d0) - e = Cx*r**(+1d0/3d0) - dedr = 1d0/3d0*Cx*r**(-2d0/3d0) - - do iEns=1,nEns + LZx(ispin) = LZx(ispin) - weight(iG)*dedr*r*r - rI = max(0d0,rho(iG,ispin,iEns)) + do iEns=1,nEns - LZx(ispin,iEns) = LZx(ispin,iEns) - weight(iG)*kappa(iEns)*dedr*r*r + rI = max(0d0,rho(iG,ispin,iEns)) - if(rI > threshold) Ex(ispin,iEns) = Ex(ispin,iEns) + weight(iG)*(kappa(iEns)*e+dedr*r)*rI - - end do + if(rI > threshold) Ex(ispin,iEns) = Ex(ispin,iEns) + weight(iG)*(e+dedr*r)*rI - endif - - else - - if(r > threshold) then - - e = Cx*r**(+1d0/3d0) - dedr = 1d0/3d0*Cx*r**(-2d0/3d0) - - LZx(ispin,:) = LZx(ispin,:) - weight(iG)*dedr*r*r - - do iEns=1,nEns - - rI = max(0d0,rho(iG,ispin,iEns)) - - if(rI > threshold) Ex(ispin,iEns) = Ex(ispin,iEns) + weight(iG)*(e+dedr*r)*rI - - end do - - endif + end do endif diff --git a/src/eDFT/US51_lda_exchange_individual_energy.f90 b/src/eDFT/US51_lda_exchange_individual_energy.f90 index 0016cee..7eb6221 100644 --- a/src/eDFT/US51_lda_exchange_individual_energy.f90 +++ b/src/eDFT/US51_lda_exchange_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine US51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,doNcentered,kappa,LZx,Ex) +subroutine US51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,LZx,Ex) ! Compute the restricted version of Slater's LDA exchange individual energy @@ -12,8 +12,6 @@ subroutine US51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,doNcen double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rhow(nGrid,nspin) double precision,intent(in) :: rho(nGrid,nspin,nEns) - logical,intent(in) :: doNcentered - double precision,intent(in) :: kappa(nEns) ! Local variables @@ -27,10 +25,10 @@ subroutine US51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,doNcen ! Output variables - double precision,intent(out) :: LZx(nspin,nEns) + double precision,intent(out) :: LZx(nspin) double precision,intent(out) :: Ex(nspin,nEns) - LZx(:,:) = 0d0 + LZx(:) = 0d0 Ex(:,:) = 0d0 do ispin=1,nspin @@ -39,33 +37,12 @@ subroutine US51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,doNcen r = max(0d0,rhow(iG,ispin)) -! if(doNcentered) then - -! if(r > threshold) then - -! e = CxLSDA*r**(+1d0/3d0) -! dedr = 1d0/3d0*CxLSDA*r**(-2d0/3d0) - -! do iEns=1,nEns - -! rI = max(0d0,rho(iG,ispin,iEns)) - -! LZx(ispin,iEns) = LZx(ispin,iEns) - weight(iG)*kappa(iEns)*dedr*r*r - -! if(rI > threshold) Ex(ispin,iEns) = Ex(ispin,iEns) + weight(iG)*(kappa(iEns)*e+dedr*r)*rI - -! end do - -! endif - -! else - if(r > threshold) then e = CxLSDA*r**(+1d0/3d0) dedr = 1d0/3d0*CxLSDA*r**(-2d0/3d0) - LZx(ispin,:) = LZx(ispin,:) - weight(iG)*dedr*r*r + LZx(ispin) = LZx(ispin) - weight(iG)*dedr*r*r do iEns=1,nEns @@ -77,8 +54,6 @@ subroutine US51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,doNcen endif -! endif - enddo enddo diff --git a/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 b/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 index 0161573..05447ea 100644 --- a/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 +++ b/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine UVWN5_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,doNcentered,LZc,Ec) +subroutine UVWN5_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,LZc,Ec) ! Compute VWN5 LDA correlation potential @@ -13,7 +13,6 @@ subroutine UVWN5_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,do double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rhow(nGrid,nspin) double precision,intent(in) :: rho(nGrid,nspin,nEns) - logical,intent(in) :: doNcentered ! Local variables diff --git a/src/eDFT/print_unrestricted_individual_energy.f90 b/src/eDFT/print_unrestricted_individual_energy.f90 index 7e679f8..f27e907 100644 --- a/src/eDFT/print_unrestricted_individual_energy.f90 +++ b/src/eDFT/print_unrestricted_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux,ExDD,EcDD,E,Om,Omx,Omc,Omaux,OmxDD,OmcDD) +subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux,ExDD,EcDD,E,Om,OmH,Omx,Omc,Omaux,OmxDD,OmcDD) ! Print individual energies for eDFT calculation @@ -20,6 +20,7 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux double precision,intent(in) :: EcDD(nsp,nEns) double precision,intent(in) :: E(nEns) + double precision,intent(in) :: OmH(nEns) double precision,intent(in) :: Omx(nEns) double precision,intent(in) :: Omc(nEns) double precision,intent(in) :: Omaux(nEns) @@ -153,32 +154,32 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux ! Total Energy and IP and EA !------------------------------------------------------------------------ -! write(*,'(A60)') '-------------------------------------------------' -! write(*,'(A60)') ' ENERGY DIFFERENCES FROM AUXILIARY ENERGIES ' -! write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') ' ENERGY DIFFERENCES FROM AUXILIARY ENERGIES ' + write(*,'(A60)') '-------------------------------------------------' -! do iEns=2,nEns -! write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',Omaux(iEns)+OmxDD(iEns)+OmcDD(iEns),' au' -! write(*,*) -! write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(iEns), ' au' -! 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 : ',OmxDD(iEns)+OmcDD(iEns), ' au' -! write(*,*) + do iEns=2,nEns + write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',Omaux(iEns)+OmxDD(iEns)+OmcDD(iEns),' au' + write(*,*) + write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(iEns), ' au' + 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 : ',OmxDD(iEns)+OmcDD(iEns), ' au' + write(*,*) -! write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') '-------------------------------------------------' -! write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',(Omaux(iEns)+OmxDD(iEns)+OmcDD(iEns))*HaToeV,' eV' -! write(*,*) -! write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(iEns)*HaToeV, ' eV' -! 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 : ',(OmxDD(iEns)+OmcDD(iEns))*HaToeV,' eV' -! write(*,*) -! end do + write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',(Omaux(iEns)+OmxDD(iEns)+OmcDD(iEns))*HaToeV,' eV' + write(*,*) + write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(iEns)*HaToeV, ' eV' + 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 : ',(OmxDD(iEns)+OmcDD(iEns))*HaToeV,' eV' + write(*,*) + end do -! write(*,'(A60)') '-------------------------------------------------' -! write(*,*) + write(*,'(A60)') '-------------------------------------------------' + write(*,*) write(*,'(A60)') '-------------------------------------------------' write(*,'(A60)') ' ENERGY DIFFERENCES FROM INDIVIDUAL ENERGIES ' @@ -191,9 +192,10 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux do iEns=2,nEns write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',Om(iEns), ' au' write(*,*) + write(*,'(A44, F16.10,A3)') ' H energy contribution : ',OmH(iEns), ' au' 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 : ',Omx(iEns)+Omc(iEns), ' au' + write(*,'(A44, F16.10,A3)') ' Hxc energy contribution : ',OmH(iEns)+Omx(iEns)+Omc(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' @@ -204,9 +206,10 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',Om(iEns)*HaToeV, ' eV' write(*,*) + write(*,'(A44, F16.10,A3)') ' H energy contribution : ',OmH(iEns)*HaToeV, ' eV' 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 : ',(Omx(iEns)+Omc(iEns))*HaToeV, ' eV' + write(*,'(A44, F16.10,A3)') ' Hxc energy contribution : ',(OmH(iEns)+Omx(iEns)+Omc(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' diff --git a/src/eDFT/unrestricted_auxiliary_energy.f90 b/src/eDFT/unrestricted_auxiliary_energy.f90 index c609caf..00046e8 100644 --- a/src/eDFT/unrestricted_auxiliary_energy.f90 +++ b/src/eDFT/unrestricted_auxiliary_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_auxiliary_energy(nBas,nEns,eps,occnum,doNcentered,Eaux) +subroutine unrestricted_auxiliary_energy(nBas,nEns,eps,occnum,Eaux) ! Compute the auxiliary KS energies @@ -12,7 +12,6 @@ subroutine unrestricted_auxiliary_energy(nBas,nEns,eps,occnum,doNcentered,Eaux) integer,intent(in) :: nEns double precision,intent(in) :: eps(nBas,nspin) double precision,intent(in) :: occnum(nBas,nspin,nEns) - logical,intent(in) :: doNcentered ! Local variables diff --git a/src/eDFT/unrestricted_correlation_derivative_discontinuity.f90 b/src/eDFT/unrestricted_correlation_derivative_discontinuity.f90 index b9dccaa..38e8d1c 100644 --- a/src/eDFT/unrestricted_correlation_derivative_discontinuity.f90 +++ b/src/eDFT/unrestricted_correlation_derivative_discontinuity.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_correlation_derivative_discontinuity(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow,kappa,Ec) +subroutine unrestricted_correlation_derivative_discontinuity(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ec) ! Compute the correlation part of the derivative discontinuity @@ -15,7 +15,6 @@ subroutine unrestricted_correlation_derivative_discontinuity(rung,DFA,nEns,wEns, double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rhow(nGrid,nspin) double precision,intent(in) :: drhow(ncart,nGrid,nspin) - double precision,intent(in) :: kappa(nEns) ! Local variables diff --git a/src/eDFT/unrestricted_correlation_individual_energy.f90 b/src/eDFT/unrestricted_correlation_individual_energy.f90 index 5183a82..a7f96be 100644 --- a/src/eDFT/unrestricted_correlation_individual_energy.f90 +++ b/src/eDFT/unrestricted_correlation_individual_energy.f90 @@ -1,5 +1,5 @@ subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight, & - rhow,drhow,rho,drho,doNcentered,LZc,Ec) + rhow,drhow,rho,drho,LZc,Ec) ! Compute the correlation energy of individual states @@ -19,7 +19,6 @@ subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns double precision,intent(in) :: drhow(ncart,nGrid,nspin) double precision,intent(in) :: rho(nGrid,nspin,nEns) double precision,intent(in) :: drho(ncart,nGrid,nspin,nEns) - logical,intent(in) :: doNcentered ! Output variables @@ -38,7 +37,7 @@ subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns case(1) - call unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,doNcentered,LZc,Ec) + call unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,LZc,Ec) ! GGA functionals @@ -56,7 +55,7 @@ subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns case(4) - call unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,doNcentered,LZc,Ec) + call unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,LZc,Ec) end select diff --git a/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 b/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 index 5bc7eec..c465f8b 100644 --- a/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 @@ -1,5 +1,5 @@ subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,nCC,aCC,nGrid,weight,rhow,drhow,& - Cx_choice,doNcentered,kappa,ExDD) + Cx_choice,doNcentered,ExDD) ! Compute the exchange part of the derivative discontinuity @@ -20,7 +20,6 @@ subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,nCC double precision,intent(in) :: drhow(ncart,nGrid) integer,intent(in) :: Cx_choice logical,intent(in) :: doNcentered - double precision,intent(in) :: kappa(nEns) ! Local variables @@ -42,7 +41,7 @@ subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,nCC case(1) call unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nCC,aCC,nGrid,weight(:),& - rhow(:),Cx_choice,doNcentered,kappa,ExDD(:)) + rhow(:),Cx_choice,doNcentered,ExDD(:)) ! GGA functionals case(2) @@ -60,7 +59,7 @@ subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,nCC case(4) call unrestricted_hybrid_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nCC,aCC,nGrid,weight(:),& - rhow(:),Cx_choice,doNcentered,kappa,ExDD(:)) + rhow(:),Cx_choice,doNcentered,ExDD(:)) end select diff --git a/src/eDFT/unrestricted_exchange_individual_energy.f90 b/src/eDFT/unrestricted_exchange_individual_energy.f90 index c9eea2d..bb29a17 100644 --- a/src/eDFT/unrestricted_exchange_individual_energy.f90 +++ b/src/eDFT/unrestricted_exchange_individual_energy.f90 @@ -1,5 +1,5 @@ subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, & - ERI,Pw,rhow,drhow,P,rho,drho,Cx_choice,doNcentered,kappa,LZx,Ex) + ERI,Pw,rhow,drhow,P,rho,drho,Cx_choice,doNcentered,LZx,Ex) ! Compute the exchange individual energy @@ -27,11 +27,10 @@ subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wE double precision,intent(in) :: drho(ncart,nGrid,nspin,nEns) integer,intent(in) :: Cx_choice logical,intent(in) :: doNcentered - double precision,intent(in) :: kappa(nEns) ! Output variables - double precision,intent(out) :: LZx(nspin,nEns) + double precision,intent(out) :: LZx(nspin) double precision,intent(out) :: Ex(nspin,nEns) select case (rung) @@ -47,7 +46,7 @@ subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wE case(1) call unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,& - rhow,rho,Cx_choice,doNcentered,kappa,LZx,Ex) + rhow,rho,Cx_choice,doNcentered,LZx,Ex) ! GGA functionals diff --git a/src/eDFT/unrestricted_hartree_individual_energy.f90 b/src/eDFT/unrestricted_hartree_individual_energy.f90 index 46a1599..04b7bcf 100644 --- a/src/eDFT/unrestricted_hartree_individual_energy.f90 +++ b/src/eDFT/unrestricted_hartree_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_hartree_individual_energy(nBas,nEns,Pw,P,ERI,doNcentered,kappa,LZH,EH) +subroutine unrestricted_hartree_individual_energy(nBas,nEns,Pw,P,ERI,LZH,EH) ! Compute the hartree contribution to the individual energies @@ -12,8 +12,6 @@ subroutine unrestricted_hartree_individual_energy(nBas,nEns,Pw,P,ERI,doNcentered double precision,intent(in) :: Pw(nBas,nBas,nspin) double precision,intent(in) :: P(nBas,nBas,nspin,nEns) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: kappa(nEns) - logical,intent(in) :: doNcentered ! Local variables @@ -26,47 +24,31 @@ subroutine unrestricted_hartree_individual_energy(nBas,nEns,Pw,P,ERI,doNcentered ! Output variables - double precision,intent(out) :: LZH(nsp,nEns) + double precision,intent(out) :: LZH(nsp) double precision,intent(out) :: EH(nsp,nEns) ! Compute HF exchange matrix allocate(J(nBas,nBas,nspin)) - LZH(:,:) = 0.d0 - EH(:,:) = 0.d0 + LZH(:) = 0.d0 + EH(:,:) = 0.d0 do ispin=1,nspin call unrestricted_hartree_potential(nBas,Pw(:,:,ispin),ERI,J(:,:,ispin)) end do + LZH(1) = - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1))) + LZH(2) = - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) & + - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1))) + LZH(3) = - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2))) + do iEns=1,nEns -! if(doNcentered) then - -! LZH(1,iEns) = - 0.5d0*kappa(iEns)*kappa(iEns)*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1))) -! LZH(2,iEns) = - 0.5d0*kappa(iEns)*kappa(iEns)*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) & -! - 0.5d0*kappa(iEns)*kappa(iEns)*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1))) -! LZH(3,iEns) = - 0.5d0*kappa(iEns)*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2))) - -! EH(1,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1))) -! EH(2,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) & -! + kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,1))) -! EH(3,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) - -! else - - LZH(1,iEns) = - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1))) - LZH(2,iEns) = - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) & - - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1))) - LZH(3,iEns) = - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2))) - - EH(1,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1))) - EH(2,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) & - + trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,1))) - EH(3,iEns) = trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) - -! endif + EH(1,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1))) + EH(2,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) & + + trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,1))) + EH(3,iEns) = trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) end do diff --git a/src/eDFT/unrestricted_hybrid_correlation_individual_energy.f90 b/src/eDFT/unrestricted_hybrid_correlation_individual_energy.f90 index 4473ded..54b73f2 100644 --- a/src/eDFT/unrestricted_hybrid_correlation_individual_energy.f90 +++ b/src/eDFT/unrestricted_hybrid_correlation_individual_energy.f90 @@ -1,5 +1,5 @@ subroutine unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight, & - rhow,drhow,rho,drho,doNcentered,LZc,Ec) + rhow,drhow,rho,drho,LZc,Ec) ! Compute the hybrid correlation energy for individual states @@ -17,7 +17,6 @@ subroutine unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid double precision,intent(in) :: drhow(ncart,nGrid) double precision,intent(in) :: rho(nGrid,nEns) double precision,intent(in) :: drho(ncart,nGrid,nEns) - logical,intent(in) :: doNcentered ! Output variables diff --git a/src/eDFT/unrestricted_hybrid_exchange_derivative_discontinuity.f90 b/src/eDFT/unrestricted_hybrid_exchange_derivative_discontinuity.f90 index 5006cb6..8fe8e55 100644 --- a/src/eDFT/unrestricted_hybrid_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/unrestricted_hybrid_exchange_derivative_discontinuity.f90 @@ -1,5 +1,5 @@ subroutine unrestricted_hybrid_exchange_derivative_discontinuity(DFA,nEns,wEns,nCC,aCC,nGrid,weight,rhow,& - Cx_choice,doNcentered,kappa,ExDD) + Cx_choice,doNcentered,ExDD) ! Compute the exchange part of the derivative discontinuity for hybrid functionals @@ -19,7 +19,6 @@ subroutine unrestricted_hybrid_exchange_derivative_discontinuity(DFA,nEns,wEns,n double precision,intent(in) :: rhow(nGrid) integer,intent(in) :: Cx_choice logical,intent(in) :: doNcentered - double precision,intent(in) :: kappa(nEns) ! Local variables diff --git a/src/eDFT/unrestricted_individual_energy.f90 b/src/eDFT/unrestricted_individual_energy.f90 index c2caa66..71bf04e 100644 --- a/src/eDFT/unrestricted_individual_energy.f90 +++ b/src/eDFT/unrestricted_individual_energy.f90 @@ -52,8 +52,8 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered double precision :: EH(nsp,nEns) double precision :: Ex(nspin,nEns) double precision :: Ec(nsp,nEns) - double precision :: LZH(nsp,nEns) - double precision :: LZx(nspin,nEns) + double precision :: LZH(nsp) + double precision :: LZx(nspin) double precision :: LZc(nsp) double precision :: Eaux(nspin,nEns) @@ -116,18 +116,18 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered ! Individual Hartree energy !------------------------------------------------------------------------ - LZH(:,:) = 0d0 + LZH(:) = 0d0 EH(:,:) = 0d0 - call unrestricted_hartree_individual_energy(nBas,nEns,Pw,P,ERI,doNcentered,kappa,LZH,EH) + call unrestricted_hartree_individual_energy(nBas,nEns,Pw,P,ERI,LZH,EH) !------------------------------------------------------------------------ ! Individual exchange energy !------------------------------------------------------------------------ - LZx(:,:) = 0d0 + LZx(:) = 0d0 Ex(:,:) = 0d0 call unrestricted_exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,ERI, & - Pw,rhow,drhow,P,rho,drho,Cx_choice,doNcentered,kappa,LZx,Ex) + Pw,rhow,drhow,P,rho,drho,Cx_choice,doNcentered,LZx,Ex) !------------------------------------------------------------------------ ! Individual correlation energy @@ -136,13 +136,13 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered LZc(:) = 0d0 Ec(:,:) = 0d0 call unrestricted_correlation_individual_energy(c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight, & - rhow,drhow,rho,drho,doNcentered,LZc,Ec) + rhow,drhow,rho,drho,LZc,Ec) !------------------------------------------------------------------------ ! Compute auxiliary energies !------------------------------------------------------------------------ - call unrestricted_auxiliary_energy(nBas,nEns,eKS,occnum,doNcentered,Eaux) + call unrestricted_auxiliary_energy(nBas,nEns,eKS,occnum,Eaux) !------------------------------------------------------------------------ ! Compute derivative discontinuities @@ -150,64 +150,84 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered do ispin=1,nspin call unrestricted_exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns,nCC,aCC,nGrid,weight, & - rhow(:,ispin),drhow(:,:,ispin),Cx_choice,doNcentered,kappa,ExDD(ispin,:)) + rhow(:,ispin),drhow(:,:,ispin),Cx_choice,doNcentered,ExDD(ispin,:)) end do - call unrestricted_correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,kappa,EcDD) + call unrestricted_correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,EcDD) !------------------------------------------------------------------------ ! Total energy !------------------------------------------------------------------------ - do iEns=1,nEns - E(iEns) = sum(Eaux(:,iEns)) & - + sum(LZH(:,iEns)) + sum(LZx(:,iEns)) + sum(LZc(:)) & - + sum(ExDD(:,iEns)) + sum(EcDD(:,iEns)) - end do + if(doNcentered) then + + do iEns=1,nEns + E(iEns) = sum(Eaux(:,iEns)) & + + kappa(iEns)*(sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:))) & + + sum(ExDD(:,iEns)) + sum(EcDD(:,iEns)) + end do + + print*,E + + else + + do iEns=1,nEns + E(iEns) = sum(Eaux(:,iEns)) & + + sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:)) & + + sum(ExDD(:,iEns)) + sum(EcDD(:,iEns)) + end do + + end if - !E(2) = (1.d0/2.d0)*(E(2)) + do iEns=1,nEns + E(iEns) = sum(ET(:,iEns)) + sum(EV(:,iEns)) & + + sum(EH(:,iEns)) + sum(Ex(:,iEns)) + sum(Ec(:,iEns)) & + + kappa(iEns)*(sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:))) & + + sum(ExDD(:,iEns)) + sum(EcDD(:,iEns)) + end do - -! print*,'test shift =',(1.d0+wEns(2)/2.d0)*(sum(Eaux(:,2)) + sum(LZH(:,2)) & -! + sum(LZx(:,2)) + sum(LZc(:)))+(1.d0-kappa(2)*wEns(2))*(sum(Eaux(:,1))& -! +sum(LZH(:,1)) + sum(LZx(:,1)) + sum(LZc(:))) - -! print*,'test=',(1.d0+wEns(2)/2.d0)*(sum(Eaux(:,2)))+(1.d0-kappa(2)*wEns(2))*(sum(Eaux(:,1))& -! ) - -! print*, 'ensemble energy=',Ew - -! Alternative way of calculating individual energies - -! do iEns=1,nEns -! E(iEns) = sum(ET(:,iEns)) + sum(EV(:,iEns)) & -! + sum(EH(:,iEns)) + sum(Ex(:,iEns)) + sum(Ec(:,iEns)) & -! + sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:)) & -! + sum(ExDD(:,iEns)) + sum(EcDD(:,iEns)) -! end do + print*,E !------------------------------------------------------------------------ ! Excitation energies !------------------------------------------------------------------------ do iEns=1,nEns + Om(iEns) = E(iEns) - E(1) - OmH(iEns) = sum(EH(:,iEns)) - sum(EH(:,1)) + OmH(iEns) = sum(EH(:,iEns)) - sum(EH(:,1)) Omx(iEns) = sum(Ex(:,iEns)) - sum(Ex(:,1)) Omc(iEns) = sum(Ec(:,iEns)) - sum(Ec(:,1)) - Omaux(iEns) = sum(Eaux(:,iEns)) - sum(Eaux(:,1)) + Omaux(iEns) = sum(Eaux(:,iEns)) - sum(Eaux(:,1)) - OmxDD(iEns) = sum(ExDD(:,iEns)) - sum(ExDD(:,1)) - OmcDD(iEns) = sum(EcDD(:,iEns)) - sum(EcDD(:,1)) + OmxDD(iEns) = sum(ExDD(:,iEns)) - sum(ExDD(:,1)) + OmcDD(iEns) = sum(EcDD(:,iEns)) - sum(EcDD(:,1)) end do + if(doNcentered) then + + do iEns=1,nEns + OmH(iEns) = OmH(iEns) + (kappa(iEns) - kappa(1))*sum(LZH(:)) + Omx(iEns) = Omx(iEns) + (kappa(iEns) - kappa(1))*sum(LZx(:)) + Omc(iEns) = Omc(iEns) + (kappa(iEns) - kappa(1))*sum(LZc(:)) + + Omaux(iEns) = Omaux(iEns) & + + (kappa(iEns) - kappa(1))*(sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:))) + + OmxDD(iEns) = kappa(iEns)*sum(ExDD(:,iEns)) - kappa(1)*sum(ExDD(:,1)) + OmcDD(iEns) = kappa(iEns)*sum(EcDD(:,iEns)) - kappa(1)*sum(EcDD(:,1)) + + end do + + end if + !------------------------------------------------------------------------ ! Dump results !------------------------------------------------------------------------ - call print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux,ExDD,EcDD,E,Om,Omx,Omc,Omaux,OmxDD,OmcDD) + call print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux,ExDD,EcDD,E,Om,OmH,Omx,Omc,Omaux,OmxDD,OmcDD) end subroutine unrestricted_individual_energy diff --git a/src/eDFT/unrestricted_lda_correlation_individual_energy.f90 b/src/eDFT/unrestricted_lda_correlation_individual_energy.f90 index b2a5645..68bec93 100644 --- a/src/eDFT/unrestricted_lda_correlation_individual_energy.f90 +++ b/src/eDFT/unrestricted_lda_correlation_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,doNcentered,LZc,Ec) +subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,LZc,Ec) ! Compute LDA correlation energy for individual states @@ -15,7 +15,6 @@ subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns, double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rhow(nGrid,nspin) double precision,intent(in) :: rho(nGrid,nspin,nEns) - logical,intent(in) :: doNcentered ! Output variables @@ -28,19 +27,19 @@ subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns, case (1) -! call UW38_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcentered,kappa,Ec) +! call UW38_lda_correlation_individual_energy(nGrid,weight,rhow,rho,LZc,Ec) case (2) -! call UPW92_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcentered,kappa,Ec) +! call UPW92_lda_correlation_individual_energy(nGrid,weight,rhow,rho,LZc,Ec) case (3) -! call UVWN3_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,doNcentered,LZc,Ec) +! call UVWN3_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,LZc,Ec) case (4) - call UVWN5_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,doNcentered,LZc,Ec) + call UVWN5_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,LZc,Ec) case default diff --git a/src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90 index 09b0b9d..0617ab7 100644 --- a/src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90 @@ -1,5 +1,5 @@ subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nCC,aCC,nGrid,weight,rhow,& - Cx_choice,doNcentered,kappa,ExDD) + Cx_choice,doNcentered,ExDD) ! Compute the exchange LDA part of the derivative discontinuity @@ -19,7 +19,6 @@ subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nCC, double precision,intent(in) :: rhow(nGrid) integer,intent(in) :: Cx_choice logical,intent(in) :: doNcentered - double precision,intent(in) :: kappa(nEns) ! Local variables @@ -38,8 +37,8 @@ subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nCC, case (2) - call UCC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weight(:),rhow(:),& - Cx_choice,doNcentered,kappa,ExDD(:)) + call UCC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weight,rhow,& + Cx_choice,doNcentered,ExDD) case default diff --git a/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 b/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 index ff766c7..67600b8 100644 --- a/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 +++ b/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 @@ -1,5 +1,5 @@ subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,rhow,& - rho,Cx_choice,doNcentered,kappa,LZx,Ex) + rho,Cx_choice,doNcentered,LZx,Ex) ! Compute LDA exchange energy for individual states @@ -20,11 +20,10 @@ subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEn double precision,intent(in) :: rho(nGrid,nspin,nEns) integer,intent(in) :: Cx_choice logical,intent(in) :: doNcentered - double precision,intent(in) :: kappa(nEns) ! Output variables - double precision :: LZx(nspin,nEns) + double precision :: LZx(nspin) double precision :: Ex(nspin,nEns) ! Select correlation functional @@ -33,12 +32,12 @@ subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEn case (1) - call US51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,doNcentered,kappa,LZx,Ex) + call US51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,LZx,Ex) case (2) call UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rhow,rho, & - Cx_choice,doNcentered,kappa,LZx,Ex) + Cx_choice,doNcentered,LZx,Ex) case default