PBE exchange

input/dft

@@ -2,23 +2,23 @@
   eDFT-UKS
 # exchange rung:
 # Hartree      = 0
-# LDA          = 1: RS51,RMFL20
+# LDA          = 1: S51,MFL20
-# GGA          = 2: RB88
+# 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

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) 

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

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

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

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