From 2ea01fd77eb37f140bfa28b3e5d87ecd967c2208 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 29 Sep 2020 11:47:18 +0200 Subject: [PATCH] scaling N-centered --- input/dft | 4 +- ...unrestricted_spatial_to_spin_MO_energy.f90 | 30 ------------- ..._lda_exchange_derivative_discontinuity.f90 | 37 ++++++++-------- .../UCC_lda_exchange_individual_energy.f90 | 29 +++++++------ .../US51_lda_exchange_individual_energy.f90 | 42 ++++++++----------- src/eDFT/eDFT.f90 | 3 +- src/eDFT/eDFT_UKS.f90 | 2 +- src/eDFT/read_options.f90 | 11 +++-- src/eDFT/unrestricted_auxiliary_energy.f90 | 2 +- ...cted_exchange_derivative_discontinuity.f90 | 7 ++-- ...nrestricted_exchange_individual_energy.f90 | 7 ++-- src/eDFT/unrestricted_individual_energy.f90 | 34 ++++++++------- ..._lda_exchange_derivative_discontinuity.f90 | 7 ++-- ...tricted_lda_exchange_individual_energy.f90 | 10 +++-- 14 files changed, 102 insertions(+), 123 deletions(-) delete mode 100644 src/QuAcK/unrestricted_spatial_to_spin_MO_energy.f90 diff --git a/input/dft b/input/dft index adfe387..636957a 100644 --- a/input/dft +++ b/input/dft @@ -29,8 +29,8 @@ 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 # Ensemble weights: wEns(1),...,wEns(nEns-1) 0.00 1.00 -# Ncentered ? 0 for NO - 0 +# N-centered? + T # Parameters for CC weight-dependent exchange functional 0.445525 0.0901503 -0.286898 0.191734 -0.0364788 -0.017035 diff --git a/src/QuAcK/unrestricted_spatial_to_spin_MO_energy.f90 b/src/QuAcK/unrestricted_spatial_to_spin_MO_energy.f90 deleted file mode 100644 index b48fbe4..0000000 --- a/src/QuAcK/unrestricted_spatial_to_spin_MO_energy.f90 +++ /dev/null @@ -1,30 +0,0 @@ -subroutine unrestricted_spatial_to_spin_MO_energy(nBas,e,nBas2,se) - -! Convert MO energies from unrestricted spatial to spin orbitals - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nBas - integer,intent(in) :: nBas2 - double precision,intent(in) :: e(nBas,nspin) - -! Local variables - - integer :: p - -! Output variables - - double precision,intent(out) :: se(nBas2) - - do p=1,nBas2,2 - se(p) = e(p,1) - enddo - - do p=2,nBas2,2 - se(p) = e(p,2) - enddo - -end subroutine unrestricted_spatial_to_spin_MO_energy diff --git a/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 index 6e164d8..660e98b 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,aCC_w1,aCC_w2,nGrid,weight,rhow,Cx_choice,doNcentered,ExDD) +subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,Cx_choice,doNcentered,kappa,ExDD) ! Compute the unrestricted version of the curvature-corrected exchange ensemble derivative @@ -15,8 +15,8 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rhow(nGrid) integer,intent(in) :: Cx_choice - integer,intent(in) :: doNcentered - + logical,intent(in) :: doNcentered + double precision,intent(in) :: kappa(nEns) ! Local variables @@ -60,7 +60,6 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr b2 = aCC_w2(2) c2 = aCC_w2(3) - nElw = electron_number(nGrid,weight,rhow) @@ -72,18 +71,21 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr w2 = wEns(3) select case (Cx_choice) - case(1) - dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0))) - dCxdw2 = 0.d0 - case(2) - dCxdw1 = 0.d0 - dCxdw2 =(0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0))) - case(3) - dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0))) & - * (1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)) - dCxdw2 = (1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)) & - * (0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0))) + case(1) + dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0))) + dCxdw2 = 0.d0 + + case(2) + dCxdw1 = 0.d0 + dCxdw2 =(0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0))) + + case(3) + dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0))) & + * (1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)) + + dCxdw2 = (1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)) & + * (0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0))) end select @@ -116,9 +118,6 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr end do end do - if (doNcentered .NE. 0) then - ExDD(2) = ((nElw-1)/nElw)*ExDD(2) - ExDD(3) = ((nElw+1)/nElw)*ExDD(3) - end if + 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 019564d..98072ca 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,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Cx_choice,doNcentered,Ex) +subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Cx_choice,doNcentered,kappa,Ex) ! Compute the unrestricted version of the curvature-corrected exchange functional @@ -16,7 +16,8 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig double precision,intent(in) :: rhow(nGrid) double precision,intent(in) :: rho(nGrid) integer,intent(in) :: Cx_choice - integer,intent(in) :: doNcentered + logical,intent(in) :: doNcentered + double precision,intent(in) :: kappa(nEns) ! Local variables @@ -61,12 +62,16 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2) select case (Cx_choice) - case(1) - Cx = alpha*Fx1 - case(2) - Cx = alpha*Fx2 - case(3) - Cx = alpha*Fx2*Fx1 + + case(1) + Cx = alpha*Fx1 + + case(2) + Cx = alpha*Fx2 + + case(3) + Cx = alpha*Fx2*Fx1 + end select nEli = electron_number(nGrid,weight,rho) @@ -87,15 +92,15 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig e_p = Cx*r**(1d0/3d0) dedr = 1d0/3d0*Cx*r**(-2d0/3d0) - if (doNcentered == 0) then - Ex = Ex - weight(iG)*dedr*r*r - else + if (doNcentered) then Ex = Ex - weight(iG)*dedr*r*r*(nEli/nElw) + else + Ex = Ex - weight(iG)*dedr*r*r end if if(rI > threshold) then - if (doNcentered == 0) then + if (doNcentered) then Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI) else Ex = Ex + weight(iG)*((nEli/nElw)*e_p*rI + dedr*r*rI) diff --git a/src/eDFT/US51_lda_exchange_individual_energy.f90 b/src/eDFT/US51_lda_exchange_individual_energy.f90 index 389af78..36a8acf 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(nGrid,weight,rhow,rho,doNcentered,Ex) +subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered,kappa,Ex) ! Compute the restricted version of Slater's LDA exchange individual energy @@ -11,34 +11,28 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rhow(nGrid) double precision,intent(in) :: rho(nGrid) - integer,intent(in) :: doNcentered - + logical,intent(in) :: doNcentered + double precision,intent(in) :: kappa ! Local variables integer :: iG double precision :: r,rI,alpha double precision :: e,dedr - double precision :: nEli,nElw + double precision :: Exrr,ExrI,ExrrI ! Output variables double precision,intent(out) :: Ex -! External variable - - double precision,external :: electron_number - - nEli = electron_number(nGrid,weight,rho) - - nElw = electron_number(nGrid,weight,rhow) - - ! Compute LDA exchange matrix in the AO basis - alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) + alpha = - (3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) + + Exrr = 0d0 + ExrI = 0d0 + ExrrI = 0d0 - Ex = 0d0 do iG=1,nGrid r = max(0d0,rhow(iG)) @@ -49,19 +43,12 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered e = alpha*r**(1d0/3d0) dedr = 1d0/3d0*alpha*r**(-2d0/3d0) - if (doNcentered == 0) then - Ex = Ex - weight(iG)*dedr*r*r - else - Ex = Ex - weight(iG)*dedr*r*r*(nEli/nElw) - end if + Exrr = Exrr - weight(iG)*dedr*r*r if(rI > threshold) then - if (doNcentered == 0) then - Ex = Ex + weight(iG)*(e*rI + dedr*r*rI) - else - Ex = Ex + weight(iG)*((nEli/nElw)*e*rI + dedr*r*rI) - end if + ExrI = ExrI + weight(iG)*e*rI + ExrrI = ExrrI + weight(iG)*dedr*r*rI endif @@ -69,4 +56,9 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered enddo + Exrr = kappa*Exrr + ExrI = kappa*ExrI + + Ex = Exrr + ExrI + ExrrI + end subroutine US51_lda_exchange_individual_energy diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index 092bdab..6fe4fc5 100644 --- a/src/eDFT/eDFT.f90 +++ b/src/eDFT/eDFT.f90 @@ -50,7 +50,8 @@ program eDFT double precision :: start_KS,end_KS,t_KS double precision :: start_int,end_int,t_int - integer :: nEns,doNcentered + integer :: nEns + logical :: doNcentered double precision,allocatable :: wEns(:) integer :: maxSCF,max_diis diff --git a/src/eDFT/eDFT_UKS.f90 b/src/eDFT/eDFT_UKS.f90 index d5e3d2a..62588d4 100644 --- a/src/eDFT/eDFT_UKS.f90 +++ b/src/eDFT/eDFT_UKS.f90 @@ -31,7 +31,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig double precision,intent(in) :: ENuc double precision,intent(in) :: occnum(nBas,nspin,nEns) integer,intent(in) :: Cx_choice - integer,intent(in) :: doNcentered + logical,intent(in) :: doNcentered ! Local variables diff --git a/src/eDFT/read_options.f90 b/src/eDFT/read_options.f90 index 8afe913..cb1d05b 100644 --- a/src/eDFT/read_options.f90 +++ b/src/eDFT/read_options.f90 @@ -25,7 +25,7 @@ subroutine read_options(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_ character(len=12),intent(out) :: x_DFA, c_DFA integer,intent(out) :: SGn integer,intent(out) :: nEns - integer,intent(out) :: doNcentered + logical,intent(out) :: doNcentered double precision,intent(out) :: wEns(maxEns) double precision,intent(out) :: aCC_w1(3) double precision,intent(out) :: aCC_w2(3) @@ -131,14 +131,17 @@ subroutine read_options(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_ end do print*,'nEl' print*,nEl - + + doNcentered = .false. read(1,*) read(1,*) (wEns(iEns),iEns=2,nEns) read(1,*) - read(1,*) doNcentered + read(1,*) answer - if (doNcentered == 0) then + if(answer == 'T') doNcentered = .true. + + if (doNcentered) then wEns(1) = 1d0 - wEns(2) - wEns(3) diff --git a/src/eDFT/unrestricted_auxiliary_energy.f90 b/src/eDFT/unrestricted_auxiliary_energy.f90 index 34a6e8e..a5f2b9a 100644 --- a/src/eDFT/unrestricted_auxiliary_energy.f90 +++ b/src/eDFT/unrestricted_auxiliary_energy.f90 @@ -12,7 +12,7 @@ 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) - integer,intent(in) :: doNcentered + logical,intent(in) :: doNcentered ! Local variables diff --git a/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 b/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 index 39547c7..a52f093 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,aCC_w1,aCC_w2,nGrid,weight,rhow,drhow,& - Cx_choice,doNcentered,ExDD) + Cx_choice,doNcentered,kappa,ExDD) ! Compute the exchange part of the derivative discontinuity @@ -19,7 +19,8 @@ subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC double precision,intent(in) :: rhow(nGrid) double precision,intent(in) :: drhow(ncart,nGrid) integer,intent(in) :: Cx_choice - integer,intent(in) :: doNcentered + logical,intent(in) :: doNcentered + double precision,intent(in) :: kappa(nEns) ! Local variables @@ -41,7 +42,7 @@ subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC case(1) call unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),& - rhow(:),Cx_choice,doNcentered,ExDD(:)) + rhow(:),Cx_choice,doNcentered,kappa,ExDD(:)) ! GGA functionals diff --git a/src/eDFT/unrestricted_exchange_individual_energy.f90 b/src/eDFT/unrestricted_exchange_individual_energy.f90 index dba469f..65263be 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,aCC_w1,aCC_w2,nGrid,weight,nBas, & - ERI,Pw,P,rhow,drhow,rho,drho,Cx_choice,doNcentered,Ex) + ERI,Pw,P,rhow,drhow,rho,drho,Cx_choice,doNcentered,kappa,Ex) ! Compute the exchange individual energy @@ -26,7 +26,8 @@ subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wE double precision,intent(in) :: rho(nGrid) double precision,intent(in) :: drho(ncart,nGrid) integer,intent(in) :: Cx_choice - integer,intent(in) :: doNcentered + logical,intent(in) :: doNcentered + double precision,intent(in) :: kappa ! Local variables @@ -51,7 +52,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,aCC_w1,aCC_w2,nGrid,weight,& - rhow,rho,Cx_choice,doNcentered,ExLDA) + rhow,rho,Cx_choice,doNcentered,kappa,ExLDA) Ex = ExLDA diff --git a/src/eDFT/unrestricted_individual_energy.f90 b/src/eDFT/unrestricted_individual_energy.f90 index 10efd57..d0fcd99 100644 --- a/src/eDFT/unrestricted_individual_energy.f90 +++ b/src/eDFT/unrestricted_individual_energy.f90 @@ -43,7 +43,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered double precision :: Ew double precision,intent(in) :: occnum(nBas,nspin,nEns) integer,intent(in) :: Cx_choice - integer,intent(in) :: doNcentered + logical,intent(in) :: doNcentered ! Local variables @@ -68,6 +68,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered integer :: ispin,iEns,iBas double precision,allocatable :: nEl(:) + double precision,allocatable :: kappa(:) double precision,external :: electron_number @@ -78,12 +79,14 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered double precision,intent(out) :: Om(nEns) - allocate(nEl(nEns)) + allocate(nEl(nEns),kappa(nEns)) + nEl(:) = 0d0 do iEns=1,nEns do iBas=1,nBas nEl(iEns) = nEl(iEns) + occnum(iBas,1,iEns) + occnum(iBas,2,iEns) end do + kappa(iEns) = nEl(iEns)/nEl(1) end do print*,'test1' @@ -94,10 +97,10 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered do ispin=1,nspin do iEns=1,nEns - if (doNcentered == 0) then - ET(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:))) + if (doNcentered) then + ET(ispin,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:))) else - ET(ispin,iEns) = (nEl(iEns)/nEl(1))*trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:))) + ET(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:))) end if end do end do @@ -109,10 +112,10 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered do iEns=1,nEns do ispin=1,nspin - if (doNcentered == 0) then - EV(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:))) + if (doNcentered) then + EV(ispin,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:))) else - EV(ispin,iEns) = (nEl(iEns)/nEl(1))*trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:))) + EV(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:))) end if end do end do @@ -138,9 +141,9 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered EJ(3,iEns) = trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) & - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2))) - if (doNcentered .NE. 0) then - EJ(:,iEns) = (nEl(iEns)/nEl(1))*EJ(:,iEns) - end if + + if(doNcentered) EJ(:,iEns) = kappa(iEns)*EJ(:,iEns) + end do print*,'test4' !------------------------------------------------------------------------ @@ -164,9 +167,10 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered do iEns=1,nEns do ispin=1,nspin - call unrestricted_exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,ERI, & - Pw(:,:,ispin),P(:,:,ispin,iEns),rhow(:,ispin),drhow(:,:,ispin), & - rho(:,ispin,iEns),drho(:,:,ispin,iEns),Cx_choice,doNcentered,Ex(ispin,iEns)) + call unrestricted_exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,ERI, & + Pw(:,:,ispin),P(:,:,ispin,iEns),rhow(:,ispin),drhow(:,:,ispin), & + rho(:,ispin,iEns),drho(:,:,ispin,iEns),Cx_choice,doNcentered,kappa(iEns), & + Ex(ispin,iEns)) end do end do @@ -222,7 +226,7 @@ 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,aCC_w1,aCC_w2,nGrid,weight, & - rhow(:,ispin),drhow(:,:,ispin),Cx_choice,doNcentered,ExDD(ispin,:)) + rhow(:,ispin),drhow(:,:,ispin),Cx_choice,doNcentered,kappa,ExDD(ispin,:)) end do call unrestricted_correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,EcDD) diff --git a/src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90 index 0e80f3a..9404a61 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,aCC_w1,aCC_w2,nGrid,weight,rhow,& - Cx_choice,doNcentered,ExDD) + Cx_choice,doNcentered,kappa,ExDD) ! Compute the exchange LDA part of the derivative discontinuity @@ -18,7 +18,8 @@ subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_ double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rhow(nGrid) integer,intent(in) :: Cx_choice - integer,intent(in) :: doNcentered + logical,intent(in) :: doNcentered + double precision,intent(in) :: kappa(nEns) ! Local variables @@ -38,7 +39,7 @@ subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_ case ('CC') call UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),& - Cx_choice,doNcentered,ExDD(:)) + Cx_choice,doNcentered,kappa,ExDD(:)) case default diff --git a/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 b/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 index 9ac05b9..fe26721 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,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,& - Cx_choice,doNcentered,Ex) + Cx_choice,doNcentered,kappa,Ex) ! Compute LDA exchange energy for individual states @@ -19,7 +19,9 @@ subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEn double precision,intent(in) :: rhow(nGrid) double precision,intent(in) :: rho(nGrid) integer,intent(in) :: Cx_choice - integer,intent(in) ::doNcentered + logical,intent(in) :: doNcentered + double precision,intent(in) :: kappa + ! Output variables @@ -31,12 +33,12 @@ subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEn case ('S51') - call US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered,Ex) + call US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered,kappa,Ex) case ('CC') call UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,& - Cx_choice,doNcentered,Ex) + Cx_choice,doNcentered,kappa,Ex) case default