diff --git a/input/dft b/input/dft index 4bbd579..4e028f5 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 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 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 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 aa4c46e..bcf936d 100644 --- a/src/eDFT/eDFT_UKS.f90 +++ b/src/eDFT/eDFT_UKS.f90 @@ -50,7 +50,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max double precision :: rcond(nspin) double precision :: ET(nspin) double precision :: EV(nspin) - double precision :: EJ(nsp) + double precision :: EH(nsp) double precision :: Ex(nspin) double precision :: Ec(nsp) double precision :: dipole(ncart) @@ -250,7 +250,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max ! Build Coulomb repulsion do ispin=1,nspin - call hartree_coulomb(nBas,Pw(:,:,ispin),ERI(:,:,:,:),J(:,:,ispin)) + call unrestricted_hartree_potential(nBas,Pw(:,:,ispin),ERI(:,:,:,:),J(:,:,ispin)) end do ! Compute exchange potential @@ -327,10 +327,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max ! Coulomb energy - EJ(1) = 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1))) - EJ(2) = 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) & - + 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1))) - EJ(3) = 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2))) + call unrestricted_hartree_energy(nBas,Pw,J,EH) ! Exchange energy @@ -346,7 +343,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max ! Total energy - Ew = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:)) + sum(Ec(:)) + Ew = sum(ET(:)) + sum(EV(:)) + sum(EH(:)) + sum(Ex(:)) + sum(Ec(:)) ! Check the grid accuracy by computing the number of electrons @@ -382,7 +379,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max ! Compute final KS energy call dipole_moment(nBas,Pw(:,:,1)+Pw(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole) - call print_UKS(nBas,nEns,nO,S,wEns,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew,dipole) + call print_UKS(nBas,nEns,nO,S,wEns,eps,c,ENuc,ET,EV,EH,Ex,Ec,Ew,dipole) ! Compute Vxc for post-HF calculations diff --git a/src/eDFT/hartree_coulomb.f90 b/src/eDFT/hartree_coulomb.f90 deleted file mode 100644 index 5c92756..0000000 --- a/src/eDFT/hartree_coulomb.f90 +++ /dev/null @@ -1,33 +0,0 @@ -subroutine hartree_coulomb(nBas,P,ERI,J) - -! Compute Coulomb matrix - - implicit none - -! Input variables - - integer,intent(in) :: nBas - double precision,intent(in) :: P(nBas,nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - -! Local variables - - integer :: mu,nu,la,si - -! Output variables - - double precision,intent(out) :: J(nBas,nBas) - - J = 0d0 - do si=1,nBas - do la=1,nBas - do nu=1,nBas - do mu=1,nBas - J(mu,nu) = J(mu,nu) + P(la,si)*ERI(mu,la,nu,si) - enddo - enddo - enddo - enddo - - -end subroutine hartree_coulomb diff --git a/src/eDFT/print_UKS.f90 b/src/eDFT/print_UKS.f90 index be30958..2615291 100644 --- a/src/eDFT/print_UKS.f90 +++ b/src/eDFT/print_UKS.f90 @@ -1,4 +1,4 @@ -subroutine print_UKS(nBas,nEns,nO,Ov,wEns,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew,dipole) +subroutine print_UKS(nBas,nEns,nO,Ov,wEns,eps,c,ENuc,ET,EV,EH,Ex,Ec,Ew,dipole) ! Print one- and two-electron energies and other stuff for KS calculation @@ -17,7 +17,7 @@ subroutine print_UKS(nBas,nEns,nO,Ov,wEns,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew,dipole) double precision,intent(in) :: ENuc double precision,intent(in) :: ET(nspin) double precision,intent(in) :: EV(nspin) - double precision,intent(in) :: EJ(nsp) + double precision,intent(in) :: EH(nsp) double precision,intent(in) :: Ex(nspin) double precision,intent(in) :: Ec(nsp) double precision,intent(in) :: Ew @@ -69,14 +69,14 @@ subroutine print_UKS(nBas,nEns,nO,Ov,wEns,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew,dipole) write(*,'(A40,1X,F16.10,A3)') ' Potential a energy: ',EV(1),' au' write(*,'(A40,1X,F16.10,A3)') ' Potential b energy: ',EV(2),' au' write(*,'(A60)') '-------------------------------------------------' - write(*,'(A40,1X,F16.10,A3)') ' Two-electron a energy: ',sum(EJ(:)) + sum(Ex(:)) + sum(Ec(:)),' au' - write(*,'(A40,1X,F16.10,A3)') ' Two-electron aa energy: ',EJ(1) + Ex(1) + Ec(1),' au' - write(*,'(A40,1X,F16.10,A3)') ' Two-electron ab energy: ',EJ(2) + Ec(2),' au' - write(*,'(A40,1X,F16.10,A3)') ' Two-electron bb energy: ',EJ(3) + Ex(2) + Ec(3),' au' - write(*,'(A40,1X,F16.10,A3)') ' Coulomb energy: ',sum(EJ(:)),' au' - write(*,'(A40,1X,F16.10,A3)') ' Coulomb aa energy: ',EJ(1),' au' - write(*,'(A40,1X,F16.10,A3)') ' Coulomb ab energy: ',EJ(2),' au' - write(*,'(A40,1X,F16.10,A3)') ' Coulomb bb energy: ',EJ(3),' au' + write(*,'(A40,1X,F16.10,A3)') ' Two-electron a energy: ',sum(EH(:)) + sum(Ex(:)) + sum(Ec(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Two-electron aa energy: ',EH(1) + Ex(1) + Ec(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Two-electron ab energy: ',EH(2) + Ec(2),' au' + write(*,'(A40,1X,F16.10,A3)') ' Two-electron bb energy: ',EH(3) + Ex(2) + Ec(3),' au' + write(*,'(A40,1X,F16.10,A3)') ' Hartree energy: ',sum(EH(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Hartree aa energy: ',EH(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Hartree ab energy: ',EH(2),' au' + write(*,'(A40,1X,F16.10,A3)') ' Hartree bb energy: ',EH(3),' au' write(*,'(A40,1X,F16.10,A3)') ' Exchange energy: ',sum(Ex(:)),' au' write(*,'(A40,1X,F16.10,A3)') ' Exchange a energy: ',Ex(1),' au' write(*,'(A40,1X,F16.10,A3)') ' Exchange b energy: ',Ex(2),' au' diff --git a/src/eDFT/unrestricted_individual_energy.f90 b/src/eDFT/unrestricted_individual_energy.f90 index 82c705e..bdedfb8 100644 --- a/src/eDFT/unrestricted_individual_energy.f90 +++ b/src/eDFT/unrestricted_individual_energy.f90 @@ -117,11 +117,12 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered ! Individual Hartree energy !------------------------------------------------------------------------ + do ispin=1,nspin + call unrestricted_hartree_potential(nBas,Pw(:,:,ispin),ERI,J(:,:,ispin)) + end do + ! do iEns=1,nEns -! do ispin=1,nspin -! call hartree_coulomb(nBas,Pw(:,:,ispin),ERI,J(:,:,ispin)) -! end do ! if(doNcentered) then ! @@ -139,7 +140,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered ! EJ(1,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1))) & -! - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),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))) &