diff --git a/input/dft b/input/dft index 4e028f5..37b6f7a 100644 --- a/input/dft +++ b/input/dft @@ -22,8 +22,8 @@ 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 +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 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diff --git a/src/eDFT/eDFT_UKS.f90 b/src/eDFT/eDFT_UKS.f90 index bcf936d..94322d8 100644 --- a/src/eDFT/eDFT_UKS.f90 +++ b/src/eDFT/eDFT_UKS.f90 @@ -247,7 +247,11 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max end if -! Build Coulomb repulsion +!------------------------------------------------------------------------ +! Compute Hxc potential and Fock operator +!------------------------------------------------------------------------ + +! Compute Hartree potential do ispin=1,nspin call unrestricted_hartree_potential(nBas,Pw(:,:,ispin),ERI(:,:,:,:),J(:,:,ispin)) @@ -266,6 +270,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max call unrestricted_correlation_potential(c_rung,c_DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rhow,drhow,Fc) ! Build Fock operator + do ispin=1,nspin F(:,:,ispin) = Hc(:,:) + J(:,:,ispin) + J(:,:,mod(ispin,2)+1) + Fx(:,:,ispin) + Fc(:,:,ispin) end do @@ -313,28 +318,28 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max ! Compute KS energy !------------------------------------------------------------------------ -! Kinetic energy +! Kinetic energy do ispin=1,nspin ET(ispin) = trace_matrix(nBas,matmul(Pw(:,:,ispin),T(:,:))) end do -! Potential energy +! Potential energy do ispin=1,nspin EV(ispin) = trace_matrix(nBas,matmul(Pw(:,:,ispin),V(:,:))) end do -! Coulomb energy +! Hartree energy call unrestricted_hartree_energy(nBas,Pw,J,EH) ! Exchange energy do ispin=1,nspin - call unrestricted_exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, & - Pw(:,:,ispin),FxHF(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Ex(ispin)& - ,Cx_choice,doNcentered) + call unrestricted_exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, & + Pw(:,:,ispin),FxHF(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Ex(ispin), & + Cx_choice,doNcentered) end do ! Correlation energy @@ -358,6 +363,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max end do write(*,*)'------------------------------------------------------------------------------------------' + !------------------------------------------------------------------------ ! End of SCF loop !------------------------------------------------------------------------ diff --git a/src/eDFT/print_unrestricted_individual_energy.f90 b/src/eDFT/print_unrestricted_individual_energy.f90 index 5ce1641..703d29b 100644 --- a/src/eDFT/print_unrestricted_individual_energy.f90 +++ b/src/eDFT/print_unrestricted_individual_energy.f90 @@ -1,5 +1,5 @@ -subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,Eaux,ExDD,EcDD,ExcDD,E, & - Om,Omx,Omc,Omxc,Omaux,OmxDD,OmcDD,OmxcDD) +subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux,ExDD,EcDD,E, & + Om,Omx,Omc,Omaux,OmxDD,OmcDD) ! Print individual energies for eDFT calculation @@ -13,14 +13,19 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc, double precision,intent(in) :: Ew double precision,intent(in) :: ET(nspin,nEns) double precision,intent(in) :: EV(nspin,nEns) - double precision,intent(in) :: EJ(nsp,nEns) - double precision,intent(in) :: Ex(nspin,nEns), Ec(nsp,nEns), Exc(nEns) + double precision,intent(in) :: EH(nsp,nEns) + double precision,intent(in) :: Ex(nspin,nEns) + double precision,intent(in) :: Ec(nsp,nEns) double precision,intent(in) :: Eaux(nspin,nEns) - double precision,intent(in) :: ExDD(nspin,nEns), EcDD(nsp,nEns), ExcDD(nsp,nEns) - double precision,intent(in) :: Omx(nEns), Omc(nEns), Omxc(nEns) - double precision,intent(in) :: Omaux(nEns) - double precision,intent(in) :: OmxDD(nEns),OmcDD(nEns),OmxcDD(nEns) + double precision,intent(in) :: ExDD(nspin,nEns) + double precision,intent(in) :: EcDD(nsp,nEns) double precision,intent(in) :: E(nEns) + + double precision,intent(in) :: Omx(nEns) + double precision,intent(in) :: Omc(nEns) + double precision,intent(in) :: Omaux(nEns) + double precision,intent(in) :: OmxDD(nEns) + double precision,intent(in) :: OmcDD(nEns) double precision,intent(in) :: Om(nEns) ! Local variables @@ -55,66 +60,66 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc, ! Kinetic energy !------------------------------------------------------------------------ -! write(*,'(A60)') '-------------------------------------------------' -! write(*,'(A60)') ' INDIVIDUAL KINETIC ENERGIES' -! write(*,'(A60)') '-------------------------------------------------' -! do iEns=1,nEns -! write(*,'(A40,I2,A2,F16.10,A3)') ' Kinetic energy state ',iEns,': ',sum(ET(:,iEns)),' au' -! end do -! write(*,'(A60)') '-------------------------------------------------' -! write(*,*) + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') ' INDIVIDUAL KINETIC ENERGIES' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Kinetic energy state ',iEns,': ',sum(ET(:,iEns)),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) !------------------------------------------------------------------------ ! Potential energy !------------------------------------------------------------------------ -! write(*,'(A60)') '-------------------------------------------------' -! write(*,'(A60)') ' INDIVIDUAL POTENTIAL ENERGIES' -! write(*,'(A60)') '-------------------------------------------------' -! do iEns=1,nEns -! write(*,'(A40,I2,A2,F16.10,A3)') ' Potential energy state ',iEns,': ',sum(EV(:,iEns)),' au' -! end do -! write(*,'(A60)') '-------------------------------------------------' -! write(*,*) + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') ' INDIVIDUAL POTENTIAL ENERGIES' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Potential energy state ',iEns,': ',sum(EV(:,iEns)),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) !------------------------------------------------------------------------ ! Hartree energy !------------------------------------------------------------------------ -! write(*,'(A60)') '-------------------------------------------------' -! write(*,'(A60)') ' INDIVIDUAL HARTREE ENERGIES' -! write(*,'(A60)') '-------------------------------------------------' -! do iEns=1,nEns -! write(*,'(A40,I2,A2,F16.10,A3)') ' Hartree energy state ',iEns,': ',sum(EJ(:,iEns)),' au' -! end do -! write(*,'(A60)') '-------------------------------------------------' -! write(*,*) + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') ' INDIVIDUAL HARTREE ENERGIES' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Hartree energy state ',iEns,': ',sum(EH(:,iEns)),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) !------------------------------------------------------------------------ ! Exchange energy !------------------------------------------------------------------------ -! write(*,'(A60)') '-------------------------------------------------' -! write(*,'(A60)') ' INDIVIDUAL EXCHANGE ENERGIES' -! write(*,'(A60)') '-------------------------------------------------' -! do iEns=1,nEns -! write(*,'(A40,I2,A2,F16.10,A3)') ' Exchange energy state ',iEns,': ',sum(Ex(:,iEns)),' au' -! end do -! write(*,'(A60)') '-------------------------------------------------' -! write(*,*) + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') ' INDIVIDUAL EXCHANGE ENERGIES' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Exchange energy state ',iEns,': ',sum(Ex(:,iEns)),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) !------------------------------------------------------------------------ ! Correlation energy !------------------------------------------------------------------------ -! write(*,'(A60)') '-------------------------------------------------' -! write(*,'(A60)') ' INDIVIDUAL CORRELATION ENERGIES' -! write(*,'(A60)') '-------------------------------------------------' -! do iEns=1,nEns -! write(*,'(A40,I2,A2,F16.10,A3)') ' Correlation energy state ',iEns,': ',sum(Ec(:,iEns)),' au' -! end do -! write(*,'(A60)') '-------------------------------------------------' -! write(*,*) + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') ' INDIVIDUAL CORRELATION ENERGIES' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Correlation energy state ',iEns,': ',sum(Ec(:,iEns)),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) !------------------------------------------------------------------------ ! Auxiliary energies @@ -140,7 +145,7 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc, write(*,*) write(*,'(A40,I2,A2,F16.10,A3)') ' x ensemble derivative state ',iEns,': ',sum(ExDD(:,iEns)), ' au' write(*,'(A40,I2,A2,F16.10,A3)') ' c ensemble derivative state ',iEns,': ',sum(EcDD(:,iEns)), ' au' - write(*,'(A40,I2,A2,F16.10,A3)') ' xc ensemble derivative state ',iEns,': ',sum(ExcDD(:,iEns)),' au' + write(*,'(A40,I2,A2,F16.10,A3)') ' xc ensemble derivative state ',iEns,': ',sum(ExDD(:,iEns))+sum(EcDD(:,iEns)),' au' end do write(*,'(A60)') '-------------------------------------------------' write(*,*) @@ -154,22 +159,22 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc, write(*,'(A60)') '-------------------------------------------------' do iEns=2,nEns - write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',Omaux(iEns)+OmxcDD(iEns),' au' + 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 : ',OmxcDD(iEns),' au' + write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxDD(iEns)+OmcDD(iEns), ' au' write(*,*) write(*,'(A60)') '-------------------------------------------------' - write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',(Omaux(iEns)+OmxcDD(iEns))*HaToeV,' eV' + 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 : ',OmxcDD(iEns)*HaToeV,' eV' + write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',(OmxDD(iEns)+OmcDD(iEns))*HaToeV,' eV' write(*,*) end do @@ -179,34 +184,34 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc, write(*,'(A60)') '-------------------------------------------------' write(*,'(A60)') ' ENERGY DIFFERENCES FROM INDIVIDUAL 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)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns) + ENuc,' au' + end do + write(*,'(A60)') '-------------------------------------------------' do iEns=2,nEns write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 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 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(*,*) 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' + write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxDD(iEns)+OmcDD(iEns),' au' write(*,*) write(*,'(A60)') '-------------------------------------------------' write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 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 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(*,*) 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' + write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',(OmxDD(iEns)+OmcDD(iEns))*HaToeV,' eV' write(*,*) end do write(*,'(A60)') '-------------------------------------------------' diff --git a/src/eDFT/unrestricted_correlation_derivative_discontinuity.f90 b/src/eDFT/unrestricted_correlation_derivative_discontinuity.f90 index 7b33ceb..b9dccaa 100644 --- a/src/eDFT/unrestricted_correlation_derivative_discontinuity.f90 +++ b/src/eDFT/unrestricted_correlation_derivative_discontinuity.f90 @@ -19,8 +19,6 @@ subroutine unrestricted_correlation_derivative_discontinuity(rung,DFA,nEns,wEns, ! Local variables - double precision :: aC - ! Output variables double precision,intent(out) :: Ec(nsp,nEns) diff --git a/src/eDFT/unrestricted_hybrid_correlation_derivative_discontinuity.f90 b/src/eDFT/unrestricted_hybrid_correlation_derivative_discontinuity.f90 index 89fecd0..09e608b 100644 --- a/src/eDFT/unrestricted_hybrid_correlation_derivative_discontinuity.f90 +++ b/src/eDFT/unrestricted_hybrid_correlation_derivative_discontinuity.f90 @@ -16,8 +16,6 @@ subroutine unrestricted_hybrid_correlation_derivative_discontinuity(DFA,nEns,wEn ! Local variables - double precision :: aC - ! Output variables double precision,intent(out) :: Ec(nsp,nEns) diff --git a/src/eDFT/unrestricted_individual_energy.f90 b/src/eDFT/unrestricted_individual_energy.f90 index bdedfb8..1932167 100644 --- a/src/eDFT/unrestricted_individual_energy.f90 +++ b/src/eDFT/unrestricted_individual_energy.f90 @@ -36,10 +36,10 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered double precision,intent(in) :: rho(nGrid,nspin,nEns) double precision,intent(in) :: drho(ncart,nGrid,nspin,nEns) - double precision,intent(in) :: J(nBas,nBas,nspin) - double precision,intent(in) :: Fx(nBas,nBas,nspin) - double precision,intent(in) :: FxHF(nBas,nBas,nspin) - double precision,intent(in) :: Fc(nBas,nBas,nspin) + double precision,intent(inout):: J(nBas,nBas,nspin) + double precision,intent(inout):: Fx(nBas,nBas,nspin) + double precision,intent(inout):: FxHF(nBas,nBas,nspin) + double precision,intent(inout):: Fc(nBas,nBas,nspin) double precision,intent(in) :: Ew double precision,intent(in) :: occnum(nBas,nspin,nEns) integer,intent(in) :: Cx_choice @@ -50,23 +50,23 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered double precision :: ET(nspin,nEns) double precision :: EV(nspin,nEns) - double precision :: EJ(nsp,nEns) + double precision :: EH(nsp,nEns) double precision :: Ex(nspin,nEns) double precision :: Ec(nsp,nEns) - double precision :: Exc(nEns) - double precision :: LZH(nspin) + double precision :: LZH(nsp) double precision :: LZx(nspin) - double precision :: LZc(nspin) - double precision :: LZHxc(nspin) + double precision :: LZc(nsp) double precision :: Eaux(nspin,nEns) double precision :: ExDD(nspin,nEns) double precision :: EcDD(nsp,nEns) - double precision :: ExcDD(nsp,nEns) - double precision :: Omx(nEns),Omc(nEns),Omxc(nEns) + double precision :: OmH(nEns) + double precision :: Omx(nEns) + double precision :: Omc(nEns) double precision :: Omaux(nEns) - double precision :: OmxDD(nEns),OmcDD(nEns),OmxcDD(nEns) + double precision :: OmxDD(nEns) + double precision :: OmcDD(nEns) double precision,external :: trace_matrix @@ -97,21 +97,21 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered ! Kinetic energy !------------------------------------------------------------------------ -! do ispin=1,nspin -! do iEns=1,nEns -! ET(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:))) -! end do -! end do + do ispin=1,nspin + do iEns=1,nEns + ET(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:))) + end do + end do !------------------------------------------------------------------------ ! Potential energy !------------------------------------------------------------------------ -! do iEns=1,nEns -! do ispin=1,nspin -! EV(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:))) -! end do -! end do + do iEns=1,nEns + do ispin=1,nspin + EV(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:))) + end do + end do !------------------------------------------------------------------------ ! Individual Hartree energy @@ -120,51 +120,68 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered do ispin=1,nspin call unrestricted_hartree_potential(nBas,Pw(:,:,ispin),ERI,J(:,:,ispin)) end do + + do iEns=1,nEns -! do iEns=1,nEns + if(doNcentered) then + + 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 -! 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 + 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 if -! EJ(1,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1))) & -! LZH(ispin) - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1))) - -! 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 - -!------------------------------------------------------------------------ -! Individual Hartree energy -!------------------------------------------------------------------------ - - do ispin=1,nspin - LZH(ispin) = -0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1)+Pw(:,:,2),J(:,:,ispin))) end do + ! Levy-Zahariev shif for Hartree + + if(doNcentered) then + + ! Do I need to scale this term? +! LZH(:) = 0d0 + + else + +! 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))) + + LZH(1) = - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1))) & + - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1))) + LZH(2) = 0d0 +! print*,- 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1)+Pw(:,:,2),J(:,:,2))) + print*,- 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) + print*,- 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2))) + LZH(3) = - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2))) + + print*,LZH(3) + + + end if + !------------------------------------------------------------------------ ! Individual exchange energy !------------------------------------------------------------------------ + do iEns=1,nEns + do ispin=1,nspin + call unrestricted_exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, & + P(:,:,ispin,iEns),FxHF(:,:,ispin),rho(:,ispin,iEns),drho(:,:,ispin,iEns), & + Ex(ispin,iEns),Cx_choice,doNcentered) + end do + end do + + ! Levy-Zahariev shif for exchange + do ispin=1,nspin call unrestricted_exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,ERI, & Pw(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Cx_choice,doNcentered, & @@ -175,14 +192,14 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered ! Individual correlation energy !------------------------------------------------------------------------ + do iEns=1,nEns + call unrestricted_correlation_energy(c_rung,c_DFA,nEns,wEns,nGrid,weight,rho(:,:,iEns),drho(:,:,:,iEns),Ec(:,iEns)) + end do + + ! Levy-Zahariev shif for correlation + call unrestricted_correlation_individual_energy(c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,drhow,doNcentered,LZc) -!------------------------------------------------------------------------ -! Individual exchange-correlation energy -!------------------------------------------------------------------------ - - LZHxc(:) = LZH(:) + LZx(:) + LZc(:) - !------------------------------------------------------------------------ ! Compute auxiliary energies !------------------------------------------------------------------------ @@ -194,31 +211,33 @@ 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,:)) end do call unrestricted_correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,kappa,EcDD) - ExcDD(1,:) = ExDD(1,:) + EcDD(1,:) - ExcDD(2,:) = EcDD(2,:) - ExcDD(3,:) = ExDD(2,:) + EcDD(3,:) - !------------------------------------------------------------------------ ! Total energy !------------------------------------------------------------------------ -! do iEns=1,nEns -! Exc(iEns) = sum(Ex(:,iEns)) + sum(Ec(:,iEns)) -! E(iEns) = sum(ET(:,iEns)) + sum(EV(:,iEns)) + sum(EJ(:,iEns)) & -! + sum(Ex(:,iEns)) + sum(Ec(:,iEns)) + sum(ExcDD(:,iEns)) -! end do + 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 do iEns=1,nEns - E(iEns) = sum(Eaux(:,iEns)) + sum(LZHxc(:)) + sum(ExcDD(:,iEns)) + E(iEns) = sum(Eaux(:,iEns)) & + + sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:)) & + + sum(ExDD(:,iEns)) + sum(EcDD(:,iEns)) end do + print*,E + !------------------------------------------------------------------------ ! Excitation energies !------------------------------------------------------------------------ @@ -226,16 +245,14 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered do iEns=1,nEns Om(iEns) = E(iEns) - E(1) -! Omx(iEns) = sum(Ex(:,iEns)) - sum(Ex(:,1)) -! Omc(iEns) = sum(Ec(:,iEns)) - sum(Ec(:,1)) -! Omxc(iEns) = Exc(iEns) - Exc(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)) OmxDD(iEns) = sum(ExDD(:,iEns)) - sum(ExDD(:,1)) OmcDD(iEns) = sum(EcDD(:,iEns)) - sum(EcDD(:,1)) - OmxcDD(iEns) = sum(ExcDD(:,iEns)) - sum(ExcDD(:,1)) - end do @@ -243,7 +260,6 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered ! Dump results !------------------------------------------------------------------------ - call print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,Eaux,ExDD,EcDD,ExcDD,E, & - Om,Omx,Omc,Omxc,Omaux,OmxDD,OmcDD,OmxcDD) + call print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux,ExDD,EcDD,E,Om,Omx,Omc,Omaux,OmxDD,OmcDD) end subroutine unrestricted_individual_energy