OK for PBE and G96 exchange

input/dft

@@ -2,14 +2,14 @@
   eDFT-UKS
 # exchange rung:
 # Hartree      = 0
-# LDA          = 1: S51,MFL20
+# LDA          = 1: S51,CC-S51
-# GGA          = 2: B88
+# 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

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) 

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
 

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
 

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              :: alpha,mupbe,kappa
-  double precision              :: r,g,x
+  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
 

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
 

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

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) 

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 !!!')

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 !!!')

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 !!!')

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