From 5dd286aa43700ef4f63d84a1450fa40bbcad2cc6 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 11 Feb 2021 22:32:48 +0100 Subject: [PATCH] PBE exchange --- input/dft | 12 ++--- input/methods | 2 +- input/options | 2 +- mol/h2.xyz | 2 +- src/eDFT/UPBE_gga_exchange_energy.f90 | 50 +++++++++++++++++ src/eDFT/UPBE_gga_exchange_potential.f90 | 68 ++++++++++++++++++++++++ 6 files changed, 127 insertions(+), 9 deletions(-) create mode 100644 src/eDFT/UPBE_gga_exchange_energy.f90 create mode 100644 src/eDFT/UPBE_gga_exchange_potential.f90 diff --git a/input/dft b/input/dft index 69e0fd2..73b5fde 100644 --- a/input/dft +++ b/input/dft @@ -2,23 +2,23 @@ eDFT-UKS # exchange rung: # Hartree = 0 -# LDA = 1: RS51,RMFL20 -# GGA = 2: RB88 +# LDA = 1: S51,MFL20 +# GGA = 2: B88 # Hybrid = 4 # Hartree-Fock = 666 666 HF # correlation rung: # Hartree = 0 -# LDA = 1: RVWN5,RMFL20 +# LDA = 1: VWN5,MFL20 # GGA = 2: # Hybrid = 4: # Hartree-Fock = 666 - 0 H + 0 H # quadrature grid SG-n 1 # Number of states in ensemble (nEns) 3 -# occupation numbers of orbitals nO and nO+1 +# 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 @@ -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.25 0.0 + 0.5 0.0 # N-centered? F # Parameters for CC weight-dependent exchange functional diff --git a/input/methods b/input/methods index 94ad703..9dac2ea 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF KS MOM - F F T F + T F F F # MP2* MP3 MP2-F12 F F F # CCD DCD CCSD CCSD(T) diff --git a/input/options b/input/options index c310e86..9d7b8f9 100644 --- a/input/options +++ b/input/options @@ -1,5 +1,5 @@ # HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess - 128 0.000001 T 5 1 1 T + 128 0.0000001 T 5 1 1 T # MP: # CC: maxSCF thresh DIIS n_diis diff --git a/mol/h2.xyz b/mol/h2.xyz index dfa658c..fe38514 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 H 0.0 0.0 0.0 -H 0.0 0.0 2.645875 +H 0.0 0.0 1.05835 diff --git a/src/eDFT/UPBE_gga_exchange_energy.f90 b/src/eDFT/UPBE_gga_exchange_energy.f90 new file mode 100644 index 0000000..0297501 --- /dev/null +++ b/src/eDFT/UPBE_gga_exchange_energy.f90 @@ -0,0 +1,50 @@ +subroutine UPBE_gga_exchange_energy(nGrid,weight,rho,drho,Ex) + +! Compute PBE GGA exchange energy + + implicit none + + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rho(nGrid) + double precision,intent(in) :: drho(3,nGrid) + +! Local variables + + integer :: iG + double precision :: alpha,mu,kappa + double precision :: r,g,x + +! Output variables + + double precision :: 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 + kappa = 0.804d0 + +! Compute GGA exchange energy + + Ex = 0d0 + + do iG=1,nGrid + + r = max(0d0,rho(iG)) + + if(r > threshold) then + g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2 + x = sqrt(g)/r**(4d0/3d0) + + Ex = Ex + weight(iG)*alpha*r**(4d0/3d0)*(1d0 + kappa - kappa/(1d0 + mu*x**2/kappa)) + + end if + + end do + +end subroutine UPBE_gga_exchange_energy diff --git a/src/eDFT/UPBE_gga_exchange_potential.f90 b/src/eDFT/UPBE_gga_exchange_potential.f90 new file mode 100644 index 0000000..2ab7921 --- /dev/null +++ b/src/eDFT/UPBE_gga_exchange_potential.f90 @@ -0,0 +1,68 @@ +subroutine UPBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + +! Compute PBE GGA exchange potential + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + integer,intent(in) :: nBas + double precision,intent(in) :: AO(nBas,nGrid) + double precision,intent(in) :: dAO(3,nBas,nGrid) + double precision,intent(in) :: rho(nGrid) + double precision,intent(in) :: drho(3,nGrid) + +! Local variables + + integer :: mu,nu,iG + double precision :: alpha,mupbe,kappa + double precision :: r,g,x,vAO,gAO + +! Output variables + + double precision,intent(out) :: Fx(nBas,nBas) + +! Coefficients for PBE exchange functional + + alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) + mupbe = ((1d0/2d0**(1d0/3d0))/(2d0*(3d0*pi**2)**(1d0/3d0)))**2*0.21951d0 + kappa = 0.804d0 + +! Compute GGA exchange matrix in the AO basis + + Fx(:,:) = 0d0 + + do mu=1,nBas + do nu=1,nBas + do iG=1,nGrid + + r = max(0d0,rho(iG)) + + if(r > threshold) then + + g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2 + x = sqrt(g)/r**(4d0/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)) + + 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)) & + + drho(3,iG)*(dAO(3,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(3,nu,iG)) + + gAO = weight(iG)*gAO + + Fx(mu,nu) = Fx(mu,nu) + 2d0*gAO*alpha*r**(4d0/3d0)*mupbe/(1d0 + mupbe*x**2/kappa)**2 + + end if + + end do + end do + end do + +end subroutine UPBE_gga_exchange_potential