diff --git a/src/eDFT/US51_lda_exchange_individual_energy.f90 b/src/eDFT/US51_lda_exchange_individual_energy.f90 index f1541a5..d7da92e 100644 --- a/src/eDFT/US51_lda_exchange_individual_energy.f90 +++ b/src/eDFT/US51_lda_exchange_individual_energy.f90 @@ -54,8 +54,13 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered enddo - Exrr = kappa*Exrr - ExrI = kappa*ExrI + if(doNcentered) then + + Exrr = kappa*Exrr + ExrI = kappa*ExrI + + endif + Ex = Exrr + ExrI + ExrrI diff --git a/src/eDFT/read_options_dft.f90 b/src/eDFT/read_options_dft.f90 index 8c0fbc0..a188ae1 100644 --- a/src/eDFT/read_options_dft.f90 +++ b/src/eDFT/read_options_dft.f90 @@ -136,15 +136,16 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns, if (doNcentered) then - wEns(1) = 1d0 - wEns(2) - wEns(3) + wEns(1) = 1d0 - wEns(2) - wEns(3) + wEns(2) = (nEl(1)/nEl(2))*wEns(2) + wEns(3) = (nEl(1)/nEl(3))*wEns(3) + else ! wEns(1) = 1d0 - nEl(2)/nEl(1)*wEns(2) - nEl(3)/nEl(1)*wEns(3) wEns(1) = 1d0 - wEns(2) - wEns(3) - wEns(2) = nEl(1)/nEl(2)*wEns(2) - wEns(3) = nEl(1)/nEl(3)*wEns(3) end if diff --git a/src/eDFT/unrestricted_individual_energy.f90 b/src/eDFT/unrestricted_individual_energy.f90 index f397140..4cb7d11 100644 --- a/src/eDFT/unrestricted_individual_energy.f90 +++ b/src/eDFT/unrestricted_individual_energy.f90 @@ -97,11 +97,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered do ispin=1,nspin do iEns=1,nEns - if (doNcentered) then - ET(ispin,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:))) - else ET(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:))) - end if end do end do @@ -111,11 +107,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered do iEns=1,nEns do ispin=1,nspin - if (doNcentered) then - EV(ispin,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:))) - else EV(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:))) - end if end do end do @@ -129,18 +121,32 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered call hartree_coulomb(nBas,Pw(:,:,ispin),ERI,J(:,:,ispin)) end do - EJ(1,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1))) & - - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1))) + if(doNcentered) then + + EJ(1,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1))) & + -0.5d0*kappa(iEns)*kappa(iEns)*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1))) + + EJ(2,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) & + + kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,1))) & + -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))) + + EJ(3,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) & + -0.5d0*kappa(iEns)*kappa(iEns)*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2))) + else - EJ(2,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) & - + trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,1))) & - - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) & - - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1))) - EJ(3,iEns) = trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) & - - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2))) + EJ(1,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1))) & + - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1))) - if(doNcentered) EJ(:,iEns) = kappa(iEns)*EJ(:,iEns) + EJ(2,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) & + + trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,1))) & + - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) & + - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1))) + + EJ(3,iEns) = trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) & + - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2))) + end if end do