UVWN3 for Gaussian B3LYP

input/dft

@@ -9,7 +9,7 @@
   4 BHHLYP
 # correlation rung: 
 # Hartree      = 0: H
-# LDA          = 1: VWN5,eVWN5
+# LDA          = 1: VWN3,VWN5,eVWN5
 # GGA          = 2: LYP,PBE
 # MGGA         = 3:
 # Hybrid       = 4: HF,LYP,PBE

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) 
@@ -15,7 +15,7 @@
 # G0W0* evGW* qsGW*
   F    F    F
 # G0T0 evGT qsGT 
-  T    F    F
+  F    F    F
 # MCMP2
   F
 # * unrestricted version available

input/options

@@ -9,7 +9,7 @@
 # GF: maxSCF thresh  DIIS n_diis lin eta renorm 
       256    0.00001 T    5      T   0.0  3      
 # GW/GT: maxSCF thresh  DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 
-         256   0.0000001   T    5     T   0.000 F      F     F     F   F  
+         256   0.0000001   T    5     T  0.00367493  F      F     F     F   F  
 # ACFDT: AC Kx  XBS
          F  F   T
 # BSE: BSE dBSE dTDA evDyn

src/MBPT/Bethe_Salpeter.f90

@@ -88,7 +88,7 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,
                                                            OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
       else
  
-        call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, &
+        call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,dipole_int,OmRPA,rho_RPA, &
                                                  OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
       end if
 
@@ -127,7 +127,7 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,
                                                            OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
       else
      
-        call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, &
+        call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,dipole_int,OmRPA,rho_RPA, &
                                                  OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
       end if
 

src/MBPT/Bethe_Salpeter_dynamic_perturbation.f90

@@ -1,4 +1,4 @@
-subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA,OmBSE,XpY,XmY)
+subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,dipole_int,OmRPA,rho_RPA,OmBSE,XpY,XmY)
 
 ! Compute dynamical effects via perturbation theory for BSE
 
@@ -16,6 +16,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,
   integer,intent(in)            :: nR
   integer,intent(in)            :: nS
 
+  double precision,intent(in)   :: eW(nBas)
   double precision,intent(in)   :: eGW(nBas)
   double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
   double precision,intent(in)   :: OmRPA(nS)

src/MBPT/qsUGW.f90

@@ -375,14 +375,14 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
 
 ! Deallocate memory
 
-  deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error,error_diis,F_diis)
+  deallocate(cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error,error_diis,F_diis)
 
 ! Perform BSE calculation
 
   if(BSE) then
 
     call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, &
-                                     S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eGW,eGW,EcBSE)
+                                     S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,c,eGW,eGW,EcBSE)
 
     if(exchange_kernel) then
 

src/MBPT/unrestricted_Bethe_Salpeter.f90

@@ -112,7 +112,7 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
 
     if(dBSE) &
       call unrestricted_Bethe_Salpeter_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,nS_sc, &
-                                                            eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, &
+                                                            eW,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, &
                                                             OmRPA,rho_RPA,OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
 
     deallocate(OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
@@ -148,7 +148,7 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
 
     if(dBSE) &
       call unrestricted_Bethe_Salpeter_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,nS_sc, &
-                                                            eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, &
+                                                            eW,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, &
                                                             OmRPA,rho_RPA,OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf)
 
     deallocate(OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf)

src/MBPT/unrestricted_Bethe_Salpeter_dynamic_perturbation.f90

@@ -1,4 +1,4 @@
-subroutine unrestricted_Bethe_Salpeter_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,nS_sc,eGW, &
+subroutine unrestricted_Bethe_Salpeter_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,nS_sc,eW,eGW, &
                                                             ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,      & 
                                                             OmRPA,rho_RPA,OmBSE,XpY_BSE,XmY_BSE)
 
@@ -23,6 +23,7 @@ subroutine unrestricted_Bethe_Salpeter_dynamic_perturbation(ispin,dTDA,eta,nBas,
   integer,intent(in)            :: nSt
   integer,intent(in)            :: nS_sc
 
+  double precision,intent(in)   :: eW(nBas,nspin)
   double precision,intent(in)   :: eGW(nBas,nspin)
   double precision,intent(in)   :: ERI_aaaa(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ERI_aabb(nBas,nBas,nBas,nBas)

src/eDFT/UVWN3_lda_correlation_energy.f90

@@ -0,0 +1,137 @@
+subroutine UVWN3_lda_correlation_energy(nGrid,weight,rho,Ec)
+
+! Compute unrestricted VWN3 LDA correlation energy
+
+  implicit none
+
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nGrid
+  double precision,intent(in)   :: weight(nGrid)
+  double precision,intent(in)   :: rho(nGrid,nspin)
+
+! Local variables
+
+  integer                       :: iG
+  double precision              :: ra,rb,r,rs,x,z
+  double precision              :: a_p,x0_p,xx0_p,b_p,c_p,x_p,q_p
+  double precision              :: a_f,x0_f,xx0_f,b_f,c_f,x_f,q_f
+  double precision              :: a_a,x0_a,xx0_a,b_a,c_a,x_a,q_a
+
+  double precision              :: ec_z,ec_p,ec_f,ec_a
+  double precision              :: fz,d2fz
+
+! Output variables
+
+  double precision              :: Ec(nsp)
+
+! Parameters of the functional
+
+  a_p  = +0.0621814d0/2d0
+  x0_p = -0.409286d0
+  b_p  = +13.0720d0
+  c_p  = +42.7198d0
+
+  a_f  = +0.0621814d0/4d0
+  x0_f = -0.743294d0
+  b_f  = +20.1231d0
+  c_f  = +101.578d0
+
+  a_a  = -1d0/(6d0*pi**2)
+  x0_a = -0.0047584D0
+  b_a  = 1.13107d0
+  c_a  = 13.0045d0
+
+! Initialization
+
+  Ec(:) = 0d0
+
+  do iG=1,nGrid
+
+    ra = max(0d0,rho(iG,1))
+    rb = max(0d0,rho(iG,2))
+
+!   alpha-alpha contribution
+
+    if(ra > threshold) then
+
+      rs = (4d0*pi*ra/3d0)**(-1d0/3d0)
+      x = sqrt(rs)
+
+      x_f   = x*x + b_f*x + c_f
+      xx0_f = x0_f*x0_f + b_f*x0_f + c_f
+      q_f   = sqrt(4d0*c_f - b_f*b_f)
+
+      ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) &
+                 - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
+
+      Ec(1) = Ec(1) + weight(iG)*ec_f*ra
+
+    end if
+
+!   alpha-beta contribution
+
+    if(ra > threshold .or. rb > threshold) then
+
+      r = ra + rb
+      rs = (4d0*pi*r/3d0)**(-1d0/3d0)
+      z = (ra - rb)/r
+      x = sqrt(rs)
+
+      fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
+      fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
+
+      d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
+
+      x_p = x*x + b_p*x + c_p
+      x_f = x*x + b_f*x + c_f
+      x_a = x*x + b_a*x + c_a
+
+      xx0_p = x0_p*x0_p + b_p*x0_p + c_p
+      xx0_f = x0_f*x0_f + b_f*x0_f + c_f
+      xx0_a = x0_a*x0_a + b_a*x0_a + c_a
+
+      q_p = sqrt(4d0*c_p - b_p*b_p)
+      q_f = sqrt(4d0*c_f - b_f*b_f)
+      q_a = sqrt(4d0*c_a - b_a*b_a)
+
+      ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & 
+                 - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) )
+
+      ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) & 
+                 - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
+
+      ec_a = a_a*( log(x**2/x_a) + 2d0*b_a/q_a*atan(q_a/(2d0*x + b_a)) & 
+                 - b_a*x0_a/xx0_a*( log((x - x0_a)**2/x_a) + 2d0*(b_a + 2d0*x0_a)/q_a*atan(q_a/(2d0*x + b_a)) ) )
+
+      ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
+
+      Ec(2) = Ec(2) + weight(iG)*ec_z*r
+
+    end if
+
+!   beta-beta contribution
+
+    if(rb > threshold) then
+
+      rs = (4d0*pi*rb/3d0)**(-1d0/3d0)
+      x = sqrt(rs)
+
+      x_f   = x*x + b_f*x + c_f
+      xx0_f = x0_f*x0_f + b_f*x0_f + c_f
+      q_f   = sqrt(4d0*c_f - b_f*b_f)
+
+      ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) &
+                 - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
+
+      Ec(3) = Ec(3) + weight(iG)*ec_f*rb
+
+    end if
+
+  end do
+
+  Ec(2) = Ec(2) - Ec(1) - Ec(3)
+
+end subroutine UVWN3_lda_correlation_energy

src/eDFT/UVWN3_lda_correlation_individual_energy.f90

@@ -0,0 +1,202 @@
+subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec)
+
+! Compute VWN3 LDA correlation potential
+
+  implicit none
+
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nGrid
+  double precision,intent(in)   :: weight(nGrid)
+  double precision,intent(in)   :: rhow(nGrid,nspin)
+  double precision,intent(in)   :: rho(nGrid,nspin)
+
+! Local variables
+
+  integer                       :: iG
+  double precision              :: ra,rb,r,raI,rbI,rI,rs,x,z
+  double precision              :: a_p,x0_p,xx0_p,b_p,c_p,x_p,q_p
+  double precision              :: a_f,x0_f,xx0_f,b_f,c_f,x_f,q_f
+  double precision              :: a_a,x0_a,xx0_a,b_a,c_a,x_a,q_a
+  double precision              :: dfzdz,dxdrs,dxdx_p,dxdx_f,dxdx_a,decdx_p,decdx_f,decdx_a
+  double precision              :: dzdra,dfzdra,drsdra,decdra_p,decdra_f,decdra_a,decdra
+  double precision              :: dzdrb,dfzdrb,drsdrb,decdrb_p,decdrb_f,decdrb_a,decdrb
+  double precision              :: ec_z,ec_p,ec_f,ec_a
+  double precision              :: fz,d2fz
+
+! Output variables
+
+  double precision              :: Ec(nsp)
+
+! Parameters of the functional
+
+  a_p  = +0.0621814d0/2d0
+  x0_p = -0.409286d0
+  b_p  = +13.0720d0
+  c_p  = +42.7198d0
+
+  a_f  = +0.0621814d0/4d0
+  x0_f = -0.743294d0
+  b_f  = +20.1231d0
+  c_f  = +101.578d0
+
+  a_a  = -1d0/(6d0*pi**2)
+  x0_a = -0.0047584D0
+  b_a  = +1.13107d0
+  c_a  = +13.0045d0
+
+! Initialization
+
+  Ec(:) = 0d0
+
+  do iG=1,nGrid
+
+    ra = max(0d0,rhow(iG,1))
+    rb = max(0d0,rhow(iG,2))
+
+    raI = max(0d0,rho(iG,1))
+    rbI = max(0d0,rho(iG,2))
+   
+!   spin-up contribution
+   
+    if(ra > threshold .or. raI > threshold) then
+   
+      r  = ra
+      rI = raI
+
+      rs = (4d0*pi*r/3d0)**(-1d0/3d0)
+      x = sqrt(rs)
+   
+      x_f = x*x + b_f*x + c_f
+   
+      xx0_f = x0_f*x0_f + b_f*x0_f + c_f
+   
+      q_f = sqrt(4d0*c_f - b_f*b_f)
+   
+      ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) & 
+                 - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
+   
+      drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
+      dxdrs  = 0.5d0/sqrt(rs)
+
+      dxdx_f = 2d0*x + b_f
+
+      decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f &
+                      - b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
+
+      decdra_f = drsdra*dxdrs*decdx_f
+
+      Ec(1) = Ec(1) + weight(iG)*(ec_f + decdra_f*r)*rI
+   
+    end if
+
+!   up-down contribution
+   
+    if(ra > threshold .or. raI > threshold) then
+   
+      r = ra + rb
+      rI = raI + rbI
+
+      rs = (4d0*pi*r/3d0)**(-1d0/3d0)
+      z = (ra - rb)/r
+      x = sqrt(rs)
+   
+      fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
+      fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
+   
+      d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
+   
+      x_p = x*x + b_p*x + c_p
+      x_f = x*x + b_f*x + c_f
+      x_a = x*x + b_a*x + c_a
+   
+      xx0_p = x0_p*x0_p + b_p*x0_p + c_p
+      xx0_f = x0_f*x0_f + b_f*x0_f + c_f
+      xx0_a = x0_a*x0_a + b_a*x0_a + c_a
+   
+      q_p = sqrt(4d0*c_p - b_p*b_p)
+      q_f = sqrt(4d0*c_f - b_f*b_f)
+      q_a = sqrt(4d0*c_a - b_a*b_a)
+   
+      ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & 
+                 - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) )
+   
+      ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) & 
+                 - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
+   
+      ec_a = a_a*( log(x**2/x_a) + 2d0*b_a/q_a*atan(q_a/(2d0*x + b_a)) & 
+                 - b_a*x0_a/xx0_a*( log((x - x0_a)**2/x_a) + 2d0*(b_a + 2d0*x0_a)/q_a*atan(q_a/(2d0*x + b_a)) ) )
+   
+      ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
+
+      dzdra = (1d0 - z)/r
+      dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
+      dfzdra = dzdra*dfzdz
+
+      drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
+      dxdrs  = 0.5d0/sqrt(rs)
+
+      dxdx_p = 2d0*x + b_p
+      dxdx_f = 2d0*x + b_f
+      dxdx_a = 2d0*x + b_a
+
+      decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p &
+                      - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) )
+
+      decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f &
+                      - b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
+
+      decdx_a = a_a*( 2d0/x - 4d0*b_a/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a &
+                      - b_a*x0_a/xx0_a*( 2/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
+
+      decdra_p = drsdra*dxdrs*decdx_p
+      decdra_f = drsdra*dxdrs*decdx_f
+      decdra_a = drsdra*dxdrs*decdx_a
+
+      decdra = decdra_p + decdra_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdra/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdra*z**3 &
+             + (decdra_f - decdra_p)*fz*z**4 + (ec_f - ec_p)*dfzdra*z**4 + 4d0*(ec_f - ec_p)*fz*dzdra*z**3
+
+      Ec(2) = Ec(2) + weight(iG)*(ec_z + decdra*r)*rI
+   
+    end if
+
+!   spin-down contribution
+   
+    if(rb > threshold .or. rbI > threshold) then
+ 
+      r = rb
+      rI = rbI
+
+      rs = (4d0*pi*r/3d0)**(-1d0/3d0)
+      x = sqrt(rs)
+
+      x_f = x*x + b_f*x + c_f
+
+      xx0_f = x0_f*x0_f + b_f*x0_f + c_f
+
+      q_f = sqrt(4d0*c_f - b_f*b_f)
+
+      ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) &
+                 - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
+
+      drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
+      dxdrs  = 0.5d0/sqrt(rs)
+
+      dxdx_f = 2d0*x + b_f
+
+      decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f &
+                      - b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
+
+      decdra_f = drsdra*dxdrs*decdx_f
+
+      Ec(3) = Ec(3) + weight(iG)*(ec_f + decdra_f*r)*rI  
+
+    end if
+
+  end do
+
+  Ec(2) = Ec(2) - Ec(1) - Ec(3)
+
+end subroutine UVWN3_lda_correlation_individual_energy

src/eDFT/UVWN3_lda_correlation_potential.f90

@@ -0,0 +1,202 @@
+subroutine UVWN3_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
+
+! Compute unrestricted VWN3 LDA correlation 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)   :: rho(nGrid,nspin)
+
+! Local variables
+
+  integer                       :: mu,nu,iG
+  double precision              :: ra,rb,r,rs,x,z
+  double precision              :: a_p,x0_p,xx0_p,b_p,c_p,x_p,q_p
+  double precision              :: a_f,x0_f,xx0_f,b_f,c_f,x_f,q_f
+  double precision              :: a_a,x0_a,xx0_a,b_a,c_a,x_a,q_a
+  double precision              :: dfzdz,dxdrs,dxdx_p,dxdx_f,dxdx_a,decdx_p,decdx_f,decdx_a
+  double precision              :: dzdra,dfzdra,drsdra,decdra_p,decdra_f,decdra_a,decdra
+  double precision              :: dzdrb,dfzdrb,drsdrb,decdrb_p,decdrb_f,decdrb_a,decdrb
+
+  double precision              :: ec_z,ec_p,ec_f,ec_a
+  double precision              :: fz,d2fz
+
+! Output variables
+
+  double precision              :: Fc(nBas,nBas,nspin)
+
+! Parameters of the functional
+
+  a_p  = +0.0621814d0/2d0
+  x0_p = -0.409286d0
+  b_p  = +13.0720d0
+  c_p  = +42.7198d0
+
+  a_f  = +0.0621814d0/4d0
+  x0_f = -0.743294d0
+  b_f  = +20.1231d0
+  c_f  = +101.578d0
+
+  a_a  = -1d0/(6d0*pi**2)
+  x0_a = -0.0047584D0
+  b_a  = +1.13107d0
+  c_a  = +13.0045d0
+
+! Initialization
+
+  Fc(:,:,:) = 0d0
+
+  do mu=1,nBas
+    do nu=1,nBas
+      do iG=1,nGrid
+
+        ra = max(0d0,rho(iG,1))
+        rb = max(0d0,rho(iG,2))
+       
+!       spin-up contribution
+       
+        if(ra > threshold) then
+       
+          r = ra + rb
+          rs = (4d0*pi*r/3d0)**(-1d0/3d0)
+          z = (ra - rb)/r
+          x = sqrt(rs)
+       
+          fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
+          fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
+       
+          d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
+       
+          x_p = x*x + b_p*x + c_p
+          x_f = x*x + b_f*x + c_f
+          x_a = x*x + b_a*x + c_a
+       
+          xx0_p = x0_p*x0_p + b_p*x0_p + c_p
+          xx0_f = x0_f*x0_f + b_f*x0_f + c_f
+          xx0_a = x0_a*x0_a + b_a*x0_a + c_a
+       
+          q_p = sqrt(4d0*c_p - b_p*b_p)
+          q_f = sqrt(4d0*c_f - b_f*b_f)
+          q_a = sqrt(4d0*c_a - b_a*b_a)
+       
+          ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & 
+                     - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) )
+       
+          ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) & 
+                     - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
+       
+          ec_a = a_a*( log(x**2/x_a) + 2d0*b_a/q_a*atan(q_a/(2d0*x + b_a)) & 
+                     - b_a*x0_a/xx0_a*( log((x - x0_a)**2/x_a) + 2d0*(b_a + 2d0*x0_a)/q_a*atan(q_a/(2d0*x + b_a)) ) )
+       
+          ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
+
+          dzdra = (1d0 - z)/r
+          dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
+          dfzdra = dzdra*dfzdz
+
+          drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
+          dxdrs  = 0.5d0/sqrt(rs)
+
+          dxdx_p = 2d0*x + b_p
+          dxdx_f = 2d0*x + b_f
+          dxdx_a = 2d0*x + b_a
+
+          decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p &
+                          - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) )
+
+          decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f &
+                          - b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
+
+          decdx_a = a_a*( 2d0/x - 4d0*b_a/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a &
+                          - b_a*x0_a/xx0_a*( 2/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
+
+          decdra_p = drsdra*dxdrs*decdx_p
+          decdra_f = drsdra*dxdrs*decdx_f
+          decdra_a = drsdra*dxdrs*decdx_a
+
+          decdra = decdra_p + decdra_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdra/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdra*z**3 &
+                 + (decdra_f - decdra_p)*fz*z**4 + (ec_f - ec_p)*dfzdra*z**4 + 4d0*(ec_f - ec_p)*fz*dzdra*z**3
+
+          Fc(mu,nu,1) = Fc(mu,nu,1) + weight(iG)*AO(mu,iG)*AO(nu,iG)*(ec_z + decdra*r)
+       
+        end if
+
+!       spin-down contribution
+       
+        if(rb > threshold) then
+       
+          r = ra + rb
+          rs = (4d0*pi*r/3d0)**(-1d0/3d0)
+          z = (ra - rb)/r
+          x = sqrt(rs)
+       
+          fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
+          fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
+       
+          d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
+       
+          x_p = x*x + b_p*x + c_p
+          x_f = x*x + b_f*x + c_f
+          x_a = x*x + b_a*x + c_a
+
+          xx0_p = x0_p*x0_p + b_p*x0_p + c_p
+          xx0_f = x0_f*x0_f + b_f*x0_f + c_f
+          xx0_a = x0_a*x0_a + b_a*x0_a + c_a
+       
+          q_p = sqrt(4d0*c_p - b_p*b_p)
+          q_f = sqrt(4d0*c_f - b_f*b_f)
+          q_a = sqrt(4d0*c_a - b_a*b_a)
+       
+          ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & 
+                     - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) )
+       
+          ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) & 
+                     - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
+       
+          ec_a = a_a*( log(x**2/x_a) + 2d0*b_a/q_a*atan(q_a/(2d0*x + b_a)) & 
+                     - b_a*x0_a/xx0_a*( log((x - x0_a)**2/x_a) + 2d0*(b_a + 2d0*x0_a)/q_a*atan(q_a/(2d0*x + b_a)) ) )
+
+          ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
+       
+          dzdrb = - (1d0 + z)/r
+          dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
+          dfzdrb = dzdrb*dfzdz
+
+          drsdrb = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
+          dxdrs  = 0.5d0/sqrt(rs)
+
+          dxdx_p = 2d0*x + b_p
+          dxdx_f = 2d0*x + b_f
+          dxdx_a = 2d0*x + b_a
+
+          decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p &
+                          - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) )
+
+          decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f &
+                          - b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
+
+          decdx_a = a_a*( 2d0/x - 4d0*b_a/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a &
+                          - b_a*x0_a/xx0_a*( 2/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
+
+          decdrb_p = drsdrb*dxdrs*decdx_p
+          decdrb_f = drsdrb*dxdrs*decdx_f
+          decdrb_a = drsdrb*dxdrs*decdx_a
+
+          decdrb = decdrb_p + decdrb_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdrb/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdrb*z**3 &
+                 + (decdrb_f - decdrb_p)*fz*z**4 + (ec_f - ec_p)*dfzdrb*z**4 + 4d0*(ec_f - ec_p)*fz*dzdrb*z**3
+          Fc(mu,nu,2) = Fc(mu,nu,2) + weight(iG)*AO(mu,iG)*AO(nu,iG)*(ec_z + decdrb*r)
+       
+        end if
+
+      end do
+    end do
+  end do
+
+end subroutine UVWN3_lda_correlation_potential

src/eDFT/unrestricted_hybrid_correlation_energy.f90

@@ -35,7 +35,7 @@ subroutine unrestricted_hybrid_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho
 
       aC = 0.81d0
 
-      call unrestricted_lda_correlation_energy('VWN5        ',nEns,wEns,nGrid,weight,rho,EcLDA)
+      call unrestricted_lda_correlation_energy('VWN3        ',nEns,wEns,nGrid,weight,rho,EcLDA)
       call unrestricted_gga_correlation_energy('LYP         ',nEns,wEns,nGrid,weight,rho,drho,EcGGA)
 
       Ec(:) = EcLDA(:) + aC*(EcGGA(:) - EcLDA(:)) 

src/eDFT/unrestricted_hybrid_correlation_potential.f90

@@ -42,7 +42,7 @@ subroutine unrestricted_hybrid_correlation_potential(DFA,nEns,wEns,nGrid,weight,
 
       aC = 0.81d0
 
-      call unrestricted_lda_correlation_potential('VWN5        ',nEns,wEns,nGrid,weight,nBas,AO,rho,FcLDA)
+      call unrestricted_lda_correlation_potential('VWN3        ',nEns,wEns,nGrid,weight,nBas,AO,rho,FcLDA)
       call unrestricted_gga_correlation_potential('LYP         ',nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FcGGA)
 
       Fc(:,:,:) = FcLDA(:,:,:) + aC*(FcGGA(:,:,:) - FcLDA(:,:,:))

src/eDFT/unrestricted_lda_correlation_energy.f90

@@ -34,6 +34,10 @@ subroutine unrestricted_lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec
 
 !   Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200
 
+    case ('VWN3')
+
+      call UVWN3_lda_correlation_energy(nGrid,weight,rho,Ec)
+
     case ('VWN5')
 
       call UVWN5_lda_correlation_energy(nGrid,weight,rho,Ec)

src/eDFT/unrestricted_lda_correlation_individual_energy.f90

@@ -26,6 +26,10 @@ subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,
 
 !   Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200
 
+    case ('VWN3')
+
+      call UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec)
+
     case ('VWN5')
 
       call UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec)

src/eDFT/unrestricted_lda_correlation_potential.f90

@@ -38,6 +38,10 @@ include 'parameters.h'
 
 !   Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200
 
+    case ('VWN3')
+
+      call UVWN3_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:))
+
     case ('VWN5')
 
       call UVWN5_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:))