From c14f6c9aa1747f94dbecbbf44abe5e1457132237 Mon Sep 17 00:00:00 2001 From: Clotilde Marut Date: Mon, 6 Dec 2021 10:03:28 +0100 Subject: [PATCH] CC modifs --- input/dft | 18 ++--- ..._lda_exchange_derivative_discontinuity.f90 | 12 +++- .../UCC_lda_exchange_individual_energy.f90 | 68 +++++++++++++++---- .../US51_lda_exchange_individual_energy.f90 | 63 +++++++++++------ .../print_unrestricted_individual_energy.f90 | 60 ++++++++-------- ...nrestricted_exchange_individual_energy.f90 | 7 +- ...unrestricted_hartree_individual_energy.f90 | 44 +++++++++--- src/eDFT/unrestricted_individual_energy.f90 | 26 +++++-- ...tricted_lda_exchange_individual_energy.f90 | 10 +-- 9 files changed, 208 insertions(+), 100 deletions(-) diff --git a/input/dft b/input/dft index 14b762f..688fdad 100644 --- a/input/dft +++ b/input/dft @@ -6,21 +6,21 @@ # GGA = 2: B88,G96,PBE # MGGA = 3: # Hybrid = 4 HF,B3LYP,PBE -1 S51 +1 CC-S51 # correlation rung: # Hartree = 0: H # LDA = 1: PW92,VWN3,VWN5,eVWN5 # GGA = 2: LYP,PBE # MGGA = 3: # Hybrid = 4: HF,B3LYP,PBE -1 VWN5 +0 H # quadrature grid SG-n -0 + 1 # Number of states in ensemble (nEns) -2 +4 # occupation numbers -1 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 -1 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 +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 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 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 @@ -31,13 +31,13 @@ 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) - 0.25 0.0 0.0 + 1 0.0 0.0 # Ncentered ? T # Parameters for CC weight-dependent exchange functional 4 -0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 +0.642674 -0.07818 -0.0280307 0.00144198 +0.254939 -0.0893405 0.00765581 0. 0.0 0.0 0.0 0.0 # choice of UCC exchange coefficient : 1 for Cx1, 2 for Cx2, 3 for Cx1*Cx2 1 diff --git a/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 index 27fe75f..607643c 100644 --- a/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 @@ -154,11 +154,19 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,wei do iEns=1,nEns do jEns=2,nEns - ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*dExdw(jEns) +! 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(:) + ! 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 01e0dfb..4fc1b1c 100644 --- a/src/eDFT/UCC_lda_exchange_individual_energy.f90 +++ b/src/eDFT/UCC_lda_exchange_individual_energy.f90 @@ -1,4 +1,5 @@ -subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rhow,Cx_choice,doNcentered,Ex) +subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rhow,rho,Cx_choice,doNcentered,kappa,LZx,Ex) + ! Compute the unrestricted version of the curvature-corrected exchange functional @@ -13,15 +14,17 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho double precision,intent(in) :: aCC(nCC,nEns-1) integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) - double precision,intent(in) :: rhow(nGrid) + double precision,intent(in) :: rhow(nGrid,nspin) + 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 - integer :: iG - double precision :: r - double precision :: dedr + integer :: iG,iEns,ispin + double precision :: r,rI + double precision :: e,dedr double precision :: a1,b1,c1,d1,w1 double precision :: a2,b2,c2,d2,w2 @@ -29,7 +32,7 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho ! Output variables - double precision,intent(out) :: Ex + double precision,intent(out) :: LZx(nspin,nEns), Ex(nspin,nEns) ! Defining enhancements factor for weight-dependent functionals @@ -96,19 +99,56 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho ! Compute LDA exchange matrix in the AO basis - Ex = 0d0 + Ex(:,:) = 0d0 + LZx(:,:) = 0d0 - do iG=1,nGrid + do ispin=1,nspin - r = max(0d0,rhow(iG)) + do iG=1,nGrid - if(r > threshold) then + r = max(0d0,rhow(iG,ispin)) - dedr = 1d0/3d0*Cx*r**(-2d0/3d0) + if(doNcentered) then + + if(r > threshold) then + + e = Cx*r**(+1d0/3d0) + dedr = 1d0/3d0*Cx*r**(-2d0/3d0) - Ex = Ex - weight(iG)*dedr*r*r + do iEns=1,nEns - endif + 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 = 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 + + endif + + enddo enddo diff --git a/src/eDFT/US51_lda_exchange_individual_energy.f90 b/src/eDFT/US51_lda_exchange_individual_energy.f90 index 460715b..0016cee 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,LZx,Ex) +subroutine US51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,doNcentered,kappa,LZx,Ex) ! Compute the restricted version of Slater's LDA exchange individual energy @@ -12,6 +12,8 @@ subroutine US51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,LZx,Ex 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 @@ -25,10 +27,10 @@ subroutine US51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,LZx,Ex ! Output variables - double precision,intent(out) :: LZx(nspin) + double precision,intent(out) :: LZx(nspin,nEns) double precision,intent(out) :: Ex(nspin,nEns) - LZx(:) = 0d0 + LZx(:,:) = 0d0 Ex(:,:) = 0d0 do ispin=1,nspin @@ -37,22 +39,45 @@ subroutine US51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,LZx,Ex r = max(0d0,rhow(iG,ispin)) - 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 - - 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 +! 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 + + 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 + +! endif enddo diff --git a/src/eDFT/print_unrestricted_individual_energy.f90 b/src/eDFT/print_unrestricted_individual_energy.f90 index 1141746..7e679f8 100644 --- a/src/eDFT/print_unrestricted_individual_energy.f90 +++ b/src/eDFT/print_unrestricted_individual_energy.f90 @@ -46,14 +46,14 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux ! Individual energies !------------------------------------------------------------------------ - write(*,'(A60)') '-------------------------------------------------' - write(*,'(A60)') ' INDIVIDUAL TOTAL ENERGIES' - write(*,'(A60)') '-------------------------------------------------' - do iEns=1,nEns - write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns) + ENuc,' au' - end do - write(*,'(A60)') '-------------------------------------------------' - write(*,*) +! write(*,'(A60)') '-------------------------------------------------' +! write(*,'(A60)') ' INDIVIDUAL TOTAL ENERGIES' +! write(*,'(A60)') '-------------------------------------------------' +! do iEns=1,nEns +! write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns) + ENuc,' au' +! end do +! write(*,'(A60)') '-------------------------------------------------' +! write(*,*) !------------------------------------------------------------------------ ! Kinetic energy @@ -153,32 +153,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 ' diff --git a/src/eDFT/unrestricted_exchange_individual_energy.f90 b/src/eDFT/unrestricted_exchange_individual_energy.f90 index bb29a17..c9eea2d 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,LZx,Ex) + ERI,Pw,rhow,drhow,P,rho,drho,Cx_choice,doNcentered,kappa,LZx,Ex) ! Compute the exchange individual energy @@ -27,10 +27,11 @@ 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) + double precision,intent(out) :: LZx(nspin,nEns) double precision,intent(out) :: Ex(nspin,nEns) select case (rung) @@ -46,7 +47,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,LZx,Ex) + rhow,rho,Cx_choice,doNcentered,kappa,LZx,Ex) ! GGA functionals diff --git a/src/eDFT/unrestricted_hartree_individual_energy.f90 b/src/eDFT/unrestricted_hartree_individual_energy.f90 index fe6c3a2..46a1599 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,LZH,EH) +subroutine unrestricted_hartree_individual_energy(nBas,nEns,Pw,P,ERI,doNcentered,kappa,LZH,EH) ! Compute the hartree contribution to the individual energies @@ -12,6 +12,9 @@ subroutine unrestricted_hartree_individual_energy(nBas,nEns,Pw,P,ERI,LZH,EH) 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 @@ -23,28 +26,47 @@ subroutine unrestricted_hartree_individual_energy(nBas,nEns,Pw,P,ERI,LZH,EH) ! Output variables - double precision,intent(out) :: LZH(nsp) + double precision,intent(out) :: LZH(nsp,nEns) double precision,intent(out) :: EH(nsp,nEns) ! Compute HF exchange matrix allocate(J(nBas,nBas,nspin)) + 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 - 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))) +! 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 end do diff --git a/src/eDFT/unrestricted_individual_energy.f90 b/src/eDFT/unrestricted_individual_energy.f90 index eae516a..c2caa66 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) - double precision :: LZx(nspin) + double precision :: LZH(nsp,nEns) + double precision :: LZx(nspin,nEns) 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,LZH,EH) + call unrestricted_hartree_individual_energy(nBas,nEns,Pw,P,ERI,doNcentered,kappa,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,LZx,Ex) + Pw,rhow,drhow,P,rho,drho,Cx_choice,doNcentered,kappa,LZx,Ex) !------------------------------------------------------------------------ ! Individual correlation energy @@ -161,9 +161,21 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered do iEns=1,nEns E(iEns) = sum(Eaux(:,iEns)) & - + sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:)) & + + sum(LZH(:,iEns)) + sum(LZx(:,iEns)) + sum(LZc(:)) & + sum(ExDD(:,iEns)) + sum(EcDD(:,iEns)) end do + + !E(2) = (1.d0/2.d0)*(E(2)) + + +! 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 diff --git a/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 b/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 index ca71f85..ff766c7 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,LZx,Ex) + rho,Cx_choice,doNcentered,kappa,LZx,Ex) ! Compute LDA exchange energy for individual states @@ -20,11 +20,11 @@ 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) + double precision :: LZx(nspin,nEns) double precision :: Ex(nspin,nEns) ! Select correlation functional @@ -33,12 +33,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,LZx,Ex) + call US51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,doNcentered,kappa,LZx,Ex) case (2) call UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rhow,rho, & - Cx_choice,doNcentered,LZx,Ex) + Cx_choice,doNcentered,kappa,LZx,Ex) case default