From a97b1881b7885871d41bfb1c66b759f28cf32dc8 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 12 Feb 2021 10:41:01 +0100 Subject: [PATCH] OK for PBE and G96 exchange --- input/dft | 12 +++---- input/methods | 2 +- src/eDFT/UG96_gga_exchange_energy.f90 | 2 +- src/eDFT/UG96_gga_exchange_potential.f90 | 1 - src/eDFT/UPBE_gga_exchange_energy.f90 | 10 +++--- src/eDFT/UPBE_gga_exchange_potential.f90 | 9 +++--- src/eDFT/US51_lda_exchange_energy.f90 | 11 ------- src/eDFT/unrestricted_exchange_energy.f90 | 10 +++++- ..._gga_exchange_derivative_discontinuity.f90 | 10 +++++- src/eDFT/unrestricted_gga_exchange_energy.f90 | 4 +++ .../unrestricted_gga_exchange_potential.f90 | 4 +++ .../unrestricted_mgga_exchange_energy.f90 | 32 +++++++++++++++++++ 12 files changed, 76 insertions(+), 31 deletions(-) create mode 100644 src/eDFT/unrestricted_mgga_exchange_energy.f90 diff --git a/input/dft b/input/dft index 73b5fde..42c6dab 100644 --- a/input/dft +++ b/input/dft @@ -2,14 +2,14 @@ eDFT-UKS # exchange rung: # Hartree = 0 -# LDA = 1: S51,MFL20 -# GGA = 2: B88 +# LDA = 1: S51,CC-S51 +# GGA = 2: B88,G96,PBE # Hybrid = 4 # Hartree-Fock = 666 - 666 HF + 2 PBE # correlation rung: # Hartree = 0 -# LDA = 1: VWN5,MFL20 +# LDA = 1: VWN5,eVWN5 # GGA = 2: # Hybrid = 4: # Hartree-Fock = 666 @@ -20,7 +20,7 @@ 3 # occupation numbers 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 + 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 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 @@ -28,7 +28,7 @@ 1 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 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 0.5 0.0 + 0.0 0.0 # N-centered? F # Parameters for CC weight-dependent exchange functional diff --git a/input/methods b/input/methods index 9dac2ea..94ad703 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF KS MOM - T F F F + F F T F # MP2* MP3 MP2-F12 F F F # CCD DCD CCSD CCSD(T) diff --git a/src/eDFT/UG96_gga_exchange_energy.f90 b/src/eDFT/UG96_gga_exchange_energy.f90 index 48f0a8c..1876cc5 100644 --- a/src/eDFT/UG96_gga_exchange_energy.f90 +++ b/src/eDFT/UG96_gga_exchange_energy.f90 @@ -11,7 +11,7 @@ subroutine UG96_gga_exchange_energy(nGrid,weight,rho,drho,Ex) integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rho(nGrid) - double precision,intent(in) :: drho(3,nGrid) + double precision,intent(in) :: drho(ncart,nGrid) ! Local variables diff --git a/src/eDFT/UG96_gga_exchange_potential.f90 b/src/eDFT/UG96_gga_exchange_potential.f90 index 2dd097b..c8d9b40 100644 --- a/src/eDFT/UG96_gga_exchange_potential.f90 +++ b/src/eDFT/UG96_gga_exchange_potential.f90 @@ -29,7 +29,6 @@ subroutine UG96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) beta = +1d0/137d0 - beta = 0d0 ! Compute GGA exchange matrix in the AO basis diff --git a/src/eDFT/UPBE_gga_exchange_energy.f90 b/src/eDFT/UPBE_gga_exchange_energy.f90 index 0297501..d06f537 100644 --- a/src/eDFT/UPBE_gga_exchange_energy.f90 +++ b/src/eDFT/UPBE_gga_exchange_energy.f90 @@ -16,8 +16,8 @@ subroutine UPBE_gga_exchange_energy(nGrid,weight,rho,drho,Ex) ! Local variables integer :: iG - double precision :: alpha,mu,kappa - double precision :: r,g,x + double precision :: alpha,mupbe,kappa + double precision :: r,g,s2 ! Output variables @@ -26,7 +26,7 @@ subroutine UPBE_gga_exchange_energy(nGrid,weight,rho,drho,Ex) ! Coefficients for PBE exchange functional alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) - mu = ((1d0/2d0**(1d0/3d0))/(2d0*(3d0*pi**2)**(1d0/3d0)))**2*0.21951d0 + mupbe = ((1d0/2d0**(1d0/3d0))/(2d0*(3d0*pi**2)**(1d0/3d0)))**2*0.21951d0 kappa = 0.804d0 ! Compute GGA exchange energy @@ -39,9 +39,9 @@ subroutine UPBE_gga_exchange_energy(nGrid,weight,rho,drho,Ex) if(r > threshold) then g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2 - x = sqrt(g)/r**(4d0/3d0) + s2 = g/r**(8d0/3d0) - Ex = Ex + weight(iG)*alpha*r**(4d0/3d0)*(1d0 + kappa - kappa/(1d0 + mu*x**2/kappa)) + Ex = Ex + weight(iG)*alpha*r**(4d0/3d0)*(1d0 + kappa - kappa/(1d0 + mupbe*s2/kappa)) end if diff --git a/src/eDFT/UPBE_gga_exchange_potential.f90 b/src/eDFT/UPBE_gga_exchange_potential.f90 index 2ab7921..f05e0b6 100644 --- a/src/eDFT/UPBE_gga_exchange_potential.f90 +++ b/src/eDFT/UPBE_gga_exchange_potential.f90 @@ -19,7 +19,7 @@ subroutine UPBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) integer :: mu,nu,iG double precision :: alpha,mupbe,kappa - double precision :: r,g,x,vAO,gAO + double precision :: r,g,s2,vAO,gAO ! Output variables @@ -44,12 +44,13 @@ subroutine UPBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) if(r > threshold) then g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2 - x = sqrt(g)/r**(4d0/3d0) + s2 = g/r**(8d0/3d0) vAO = weight(iG)*AO(mu,iG)*AO(nu,iG) Fx(mu,nu) = Fx(mu,nu) & - + vAO*4d0/3d0*alpha*r**(1d0/3d0)*(1d0 + kappa - kappa/(1d0 + mupbe*x**2/kappa)) + + vAO*4d0/3d0*alpha*r**(1d0/3d0)*(1d0 + kappa - kappa/(1d0 + mupbe*s2/kappa)) & + - vAO*8d0/3d0*alpha*r**(1d0/3d0)*mupbe*s2/(1d0 + mupbe*s2/kappa)**2 gAO = drho(1,iG)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) & + drho(2,iG)*(dAO(2,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(2,nu,iG)) & @@ -57,7 +58,7 @@ subroutine UPBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) gAO = weight(iG)*gAO - Fx(mu,nu) = Fx(mu,nu) + 2d0*gAO*alpha*r**(4d0/3d0)*mupbe/(1d0 + mupbe*x**2/kappa)**2 + Fx(mu,nu) = Fx(mu,nu) + 2d0*gAO*alpha*r**(-4d0/3d0)*mupbe/(1d0 + mupbe*s2/kappa)**2 end if diff --git a/src/eDFT/US51_lda_exchange_energy.f90 b/src/eDFT/US51_lda_exchange_energy.f90 index b21d673..c0b8702 100644 --- a/src/eDFT/US51_lda_exchange_energy.f90 +++ b/src/eDFT/US51_lda_exchange_energy.f90 @@ -22,21 +22,10 @@ subroutine US51_lda_exchange_energy(nGrid,weight,rho,Ex) double precision :: Ex -! Cxw2 parameters for He N->N+1 -! a2 = 0.135068d0 -! b2 = -0.00774769d0 -! c2 = -0.0278205d0 - -! Cxw1 parameters for He N->N-1 -! a1 = 0.420243d0 -! b1 = 0.0700561d0 -! c1 = -0.288301d0 - ! Cx coefficient for Slater LDA exchange alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) -! alphaw = alpha*(1d0 - wEns(2)*(1d0 - wEns(2))*(a1 + b1*(wEns(2) - 0.5d0) + c1*(wEns(2) - 0.5d0)**2)) ! Compute LDA exchange energy Ex = 0d0 diff --git a/src/eDFT/unrestricted_exchange_energy.f90 b/src/eDFT/unrestricted_exchange_energy.f90 index 91c8e52..56cf438 100644 --- a/src/eDFT/unrestricted_exchange_energy.f90 +++ b/src/eDFT/unrestricted_exchange_energy.f90 @@ -26,7 +26,7 @@ subroutine unrestricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,a ! Local variables - double precision :: ExLDA,ExGGA,ExHF + double precision :: ExLDA,ExGGA,ExMGGA,ExHF double precision :: cX,aX,aC ! Output variables @@ -58,6 +58,14 @@ subroutine unrestricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,a Ex = ExGGA +! MGGA functionals + + case(3) + + call unrestricted_mgga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,ExMGGA) + + Ex = ExMGGA + ! Hybrid functionals case(4) diff --git a/src/eDFT/unrestricted_gga_exchange_derivative_discontinuity.f90 b/src/eDFT/unrestricted_gga_exchange_derivative_discontinuity.f90 index bb0c4a7..7f451c5 100644 --- a/src/eDFT/unrestricted_gga_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/unrestricted_gga_exchange_derivative_discontinuity.f90 @@ -22,14 +22,22 @@ subroutine unrestricted_gga_exchange_derivative_discontinuity(DFA,nEns,wEns,nGri double precision,intent(out) :: ExDD(nEns) -! Select correlation functional +! Select exchange functional select case (DFA) + case ('G96') + + ExDD(:) = 0d0 + case ('B88') ExDD(:) = 0d0 + case ('PBE') + + ExDD(:) = 0d0 + case default call print_warning('!!! GGA exchange derivative discontinuity not available !!!') diff --git a/src/eDFT/unrestricted_gga_exchange_energy.f90 b/src/eDFT/unrestricted_gga_exchange_energy.f90 index b7d2679..f184555 100644 --- a/src/eDFT/unrestricted_gga_exchange_energy.f90 +++ b/src/eDFT/unrestricted_gga_exchange_energy.f90 @@ -30,6 +30,10 @@ subroutine unrestricted_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho, call UB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) + case ('PBE') + + call UPBE_gga_exchange_energy(nGrid,weight,rho,drho,Ex) + case default call print_warning('!!! GGA exchange energy not available !!!') diff --git a/src/eDFT/unrestricted_gga_exchange_potential.f90 b/src/eDFT/unrestricted_gga_exchange_potential.f90 index 368a0e4..4c72163 100644 --- a/src/eDFT/unrestricted_gga_exchange_potential.f90 +++ b/src/eDFT/unrestricted_gga_exchange_potential.f90 @@ -34,6 +34,10 @@ subroutine unrestricted_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,A call UB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + case ('PBE') + + call UPBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + case default call print_warning('!!! GGA exchange potential not available !!!') diff --git a/src/eDFT/unrestricted_mgga_exchange_energy.f90 b/src/eDFT/unrestricted_mgga_exchange_energy.f90 new file mode 100644 index 0000000..18a9cf7 --- /dev/null +++ b/src/eDFT/unrestricted_mgga_exchange_energy.f90 @@ -0,0 +1,32 @@ +subroutine unrestricted_mgga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) + +! Select MGGA exchange functional for energy calculation + + implicit none + + include 'parameters.h' + +! Input variables + + character(len=12),intent(in) :: DFA + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rho(nGrid) + double precision,intent(in) :: drho(3,nGrid) + +! Output variables + + double precision :: Ex + + select case (DFA) + + case default + + call print_warning('!!! MGGA exchange energy not available !!!') + stop + + end select + +end subroutine unrestricted_mgga_exchange_energy