diff --git a/input/dft b/input/dft
index f2a0704..3220f9d 100644
--- a/input/dft
+++ b/input/dft
@@ -6,14 +6,14 @@
 # GGA          = 2: B88,G96,PBE
 # MGGA         = 3: 
 # Hybrid       = 4: HF,B3,PBE
-  2 B88
+  1 S51
 # correlation rung: 
 # Hartree      = 0: H
 # LDA          = 1: PW92,VWN3,VWN5,eVWN5
 # GGA          = 2: LYP,PBE
 # MGGA         = 3:
 # Hybrid       = 4: HF,LYP,PBE
-  1 PW92
+ 0 H
 # quadrature grid SG-n
   1
 # Number of states in ensemble (nEns)
diff --git a/input/methods b/input/methods
index 504e70c..2452cdc 100644
--- a/input/methods
+++ b/input/methods
@@ -1,5 +1,5 @@
 # RHF UHF KS  MOM
-  T   F   F   F
+  F   T   F   F
 # MP2* MP3 MP2-F12
   F   F   F
 # CCD DCD CCSD CCSD(T) 
@@ -13,7 +13,7 @@
 # G0F2 evGF2 G0F3 evGF3
   F    F     F    F
 # G0W0* evGW* qsGW*
-  F    F    T
+  T    F    T
 # G0T0 evGT qsGT 
   F    F    F
 # MCMP2
diff --git a/src/MBPT/print_qsGW.f90 b/src/MBPT/print_qsGW.f90
index df854d6..d6be17b 100644
--- a/src/MBPT/print_qsGW.f90
+++ b/src/MBPT/print_qsGW.f90
@@ -1,4 +1,4 @@
-subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z,EcRPA,EqsGW,dipole)
+subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z,EcGM,EcRPA,EqsGW,dipole)
 
 ! Print one-electron energies and other stuff for qsGW
 
@@ -7,9 +7,18 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z
 
 ! Input variables
 
-  integer,intent(in)                 :: nBas,nO,nSCF
-  double precision,intent(in)        :: ENuc,EcRPA,Conv,thresh
-  double precision,intent(in)        :: eHF(nBas),eGW(nBas),c(nBas),P(nBas,nBas) 
+  integer,intent(in)                 :: nBas
+  integer,intent(in)                 :: nO
+  integer,intent(in)                 :: nSCF
+  double precision,intent(in)        :: ENuc
+  double precision,intent(in)        :: EcGM
+  double precision,intent(in)        :: EcRPA
+  double precision,intent(in)        :: Conv
+  double precision,intent(in)        :: thresh
+  double precision,intent(in)        :: eHF(nBas)
+  double precision,intent(in)        :: eGW(nBas)
+  double precision,intent(in)        :: c(nBas)
+  double precision,intent(in)        :: P(nBas,nBas) 
   double precision,intent(in)        :: T(nBas,nBas),V(nBas,nBas)
   double precision,intent(in)        :: J(nBas,nBas),K(nBas,nBas),F(nBas,nBas)
   double precision,intent(in)        :: Z(nBas),SigC(nBas,nBas)
@@ -37,9 +46,8 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z
   EV = trace_matrix(nBas,matmul(P,V))
   EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
   Ex = 0.25d0*trace_matrix(nBas,matmul(P,K))
-  Ec = 0d0
 ! Ec = -0.50d0*trace_matrix(nBas,matmul(P,SigC))
-  EqsGW = ET + EV + EJ + Ex + Ec
+  EqsGW = ET + EV + EJ + Ex + EcGM
 
 ! Dump results
 
@@ -69,7 +77,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z
   write(*,*)'-------------------------------------------'
   write(*,'(2X,A30,F15.6,A3)') '    qsGW total       energy:',EqsGW + ENuc,' au'
   write(*,'(2X,A30,F15.6,A3)') '    qsGW exchange    energy:',Ex,' au'
-! write(*,'(2X,A30,F15.6,A3)') '    qsGW correlation energy:',Ec,' au'
+  write(*,'(2X,A30,F15.6,A3)') '    qsGW correlation energy:',EcGM,' au'
   write(*,'(2X,A30,F15.6,A3)') 'RPA@qsGW correlation energy:',EcRPA,' au'
   write(*,*)'-------------------------------------------'
   write(*,*)
@@ -89,7 +97,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z
     write(*,'(A32,1X,F16.10,A3)') ' Two-electron energy: ',EJ + Ex,' au'
     write(*,'(A32,1X,F16.10,A3)') ' Hartree      energy: ',EJ,' au'
     write(*,'(A32,1X,F16.10,A3)') ' Exchange     energy: ',Ex,' au'
-!   write(*,'(A32,1X,F16.10,A3)') ' Correlation  energy: ',Ec,' au'
+    write(*,'(A32,1X,F16.10,A3)') ' Correlation  energy: ',EcGM,' au'
     write(*,'(A50)')           '---------------------------------------'
     write(*,'(A32,1X,F16.10,A3)') ' Electronic   energy: ',EqsGW,' au'
     write(*,'(A32,1X,F16.10,A3)') ' Nuclear   repulsion: ',ENuc,' au'
diff --git a/src/MBPT/print_qsUGW.f90 b/src/MBPT/print_qsUGW.f90
index d753e0e..144b2d9 100644
--- a/src/MBPT/print_qsUGW.f90
+++ b/src/MBPT/print_qsUGW.f90
@@ -1,5 +1,5 @@
 subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, & 
-                       ENuc,ET,EV,EJ,Ex,Ec,EcRPA,EqsGW,SigC,Z,dipole)
+                       ENuc,ET,EV,EJ,Ex,Ec,EcGM,EcRPA,EqsGW,SigC,Z,dipole)
 
 ! Print one-electron energies and other stuff for qsUGW
 
@@ -17,6 +17,7 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, &
   double precision,intent(in)        :: EJ(nsp)
   double precision,intent(in)        :: Ex(nspin)
   double precision,intent(in)        :: Ec(nsp)
+  double precision,intent(in)        :: EcGM(nspin)
   double precision,intent(in)        :: EcRPA
   double precision,intent(in)        :: EqsGW
   double precision,intent(in)        :: Conv
@@ -103,6 +104,7 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, &
   write(*,'(2X,A30,F15.6,A3)') '    qsUGW total       energy:',EqsGW + ENuc,' au'
   write(*,'(2X,A30,F15.6,A3)') '    qsUGW exchange    energy:',sum(Ex(:)),' au'
 ! write(*,'(2X,A30,F15.6,A3)') '    qsUGW correlation energy:',sum(Ec(:)),' au'
+  write(*,'(2X,A30,F15.6,A3)') ' GM@qsUGW correlation energy:',sum(EcGM(:)),' au'
   write(*,'(2X,A30,F15.6,A3)') 'RPA@qsUGW correlation energy:',EcRPA,' au'
   write(*,*)'-------------------------------------------------------------------------------& 
               -------------------------------------------------'
@@ -141,7 +143,7 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, &
     write(*,'(A40,1X,F16.10,A3)') ' Exchange        energy: ',sum(Ex(:)),' au'
     write(*,'(A40,1X,F16.10,A3)') ' Exchange     a  energy: ',Ex(1),' au'
     write(*,'(A40,1X,F16.10,A3)') ' Exchange     b  energy: ',Ex(2),' au'
-!   write(*,*)
+    write(*,*)
 !   write(*,'(A40,1X,F16.10,A3)') ' Correlation     energy: ',sum(Ec(:)),' au'
 !   write(*,'(A40,1X,F16.10,A3)') ' Correlation  aa energy: ',Ec(1),' au'
 !   write(*,'(A40,1X,F16.10,A3)') ' Correlation  ab energy: ',Ec(2),' au'
diff --git a/src/MBPT/qsGW.f90 b/src/MBPT/qsGW.f90
index e906391..db7f56d 100644
--- a/src/MBPT/qsGW.f90
+++ b/src/MBPT/qsGW.f90
@@ -235,7 +235,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
     ! Print results
 
     call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
-    call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigCp,Z,EcRPA,EqsGW,dipole)
+    call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigCp,Z,EcGM,EcRPA,EqsGW,dipole)
 
   enddo
 !------------------------------------------------------------------------
diff --git a/src/MBPT/qsUGW.f90 b/src/MBPT/qsUGW.f90
index 09d8cef..d5b84d2 100644
--- a/src/MBPT/qsUGW.f90
+++ b/src/MBPT/qsUGW.f90
@@ -77,6 +77,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
   double precision              :: Ex(nspin)
   double precision              :: Ec(nsp)
   double precision              :: EcRPA
+  double precision              :: EcGM(nspin)
   double precision              :: EqsGW
   double precision              :: EcBSE(nspin)
   double precision              :: EcAC(nspin)
@@ -230,12 +231,12 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
 
     if(G0W) then
 
-      call unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC)
+      call unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,EcGM)
       call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z)
 
      else
 
-      call unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC)
+      call unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,EcGM)
       call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,Z)
 
      endif
@@ -321,7 +322,8 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
     ! Coulomb energy
 
     EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1)))
-    EJ(2) = 1.0d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2)))
+    EJ(2) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) &
+          + 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,1)))
     EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2)))
 
     ! Exchange energy
@@ -348,7 +350,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
     !------------------------------------------------------------------------
 
     call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
-    call print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,P,S,T,V,J,K,ENuc,ET,EV,EJ,Ex,Ec,EcRPA,EqsGW,SigCp,Z,dipole)
+    call print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,P,S,T,V,J,K,ENuc,ET,EV,EJ,Ex,Ec,EcGM,EcRPA,EqsGW,SigCp,Z,dipole)
 
   enddo
 !------------------------------------------------------------------------
diff --git a/src/MBPT/self_energy_correlation.f90 b/src/MBPT/self_energy_correlation.f90
index e47a9b4..369ee32 100644
--- a/src/MBPT/self_energy_correlation.f90
+++ b/src/MBPT/self_energy_correlation.f90
@@ -16,7 +16,9 @@ subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Ec
 
 ! Local variables
 
-  integer                       :: i,j,a,b,p,x,y,jb
+  integer                       :: i,j,a,b
+  integer                       :: p,q,r
+  integer                       :: jb
   double precision              :: eps
 
 ! Output variables
@@ -26,7 +28,7 @@ subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Ec
 
 ! Initialize 
 
-  SigC = 0d0
+  SigC(:,:) = 0d0
 
 !-----------------------------!
 ! COHSEX static approximation !
@@ -36,32 +38,32 @@ subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Ec
 
    ! COHSEX: SEX of the COHSEX correlation self-energy
 
-    do x=nC+1,nBas-nR
-      do y=nC+1,nBas-nR
+    do p=nC+1,nBas-nR
+      do q=nC+1,nBas-nR
         do i=nC+1,nO
           do jb=1,nS
-            SigC(x,y) = SigC(x,y) + 4d0*rho(x,i,jb)*rho(y,i,jb)/Omega(jb)
-          enddo
-        enddo
-      enddo
-    enddo
+            SigC(p,q) = SigC(p,q) + 4d0*rho(p,i,jb)*rho(q,i,jb)/Omega(jb)
+          end do
+        end do
+      end do
+    end do
  
     ! COHSEX: COH part of the COHSEX correlation self-energy
  
-    do x=nC+1,nBas-nR
-      do y=nC+1,nBas-nR
-        do p=nC+1,nBas-nR
+    do p=nC+1,nBas-nR
+      do q=nC+1,nBas-nR
+        do r=nC+1,nBas-nR
           do jb=1,nS
-            SigC(x,y) = SigC(x,y) - 2d0*rho(x,p,jb)*rho(y,p,jb)/Omega(jb)
-          enddo
-        enddo
-      enddo
-    enddo
+            SigC(p,q) = SigC(p,q) - 2d0*rho(p,r,jb)*rho(q,r,jb)/Omega(jb)
+          end do
+        end do
+      end do
+    end do
 
     EcGM = 0d0
     do i=nC+1,nO
       EcGM = EcGM + 0.5d0*SigC(i,i)
-    enddo
+    end do
 
   else
 
@@ -71,30 +73,42 @@ subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Ec
 
   ! Occupied part of the correlation self-energy
 
-    do x=nC+1,nBas-nR
-      do y=nC+1,nBas-nR
+    do p=nC+1,nBas-nR
+      do q=nC+1,nBas-nR
         do i=nC+1,nO
           do jb=1,nS
-            eps = e(x) - e(i) + Omega(jb)
-            SigC(x,y) = SigC(x,y) + 2d0*rho(x,i,jb)*rho(y,i,jb)*eps/(eps**2 + eta**2)
-          enddo
-        enddo
-      enddo
-    enddo
+            eps = e(p) - e(i) + Omega(jb)
+            SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,jb)*rho(q,i,jb)*eps/(eps**2 + eta**2)
+          end do
+        end do
+      end do
+    end do
  
     ! Virtual part of the correlation self-energy
  
-    do x=nC+1,nBas-nR
-      do y=nC+1,nBas-nR
+    do p=nC+1,nBas-nR
+      do q=nC+1,nBas-nR
         do a=nO+1,nBas-nR
           do jb=1,nS
-            eps = e(x) - e(a) - Omega(jb)
-            SigC(x,y) = SigC(x,y) + 2d0*rho(x,a,jb)*rho(y,a,jb)*eps/(eps**2 + eta**2)
-          enddo
-        enddo
-      enddo
-    enddo
+            eps = e(p) - e(a) - Omega(jb)
+            SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,jb)*rho(q,a,jb)*eps/(eps**2 + eta**2)
+          end do
+        end do
+      end do
+    end do
 
-  endif
+    ! GM correlation energy
+
+    EcGM = 0d0
+    do i=nC+1,nO
+      do a=nO+1,nBas-nR
+        do jb=1,nS
+          eps = e(a) - e(i) + Omega(jb)
+          EcGM = EcGM - 4d0*rho(a,i,jb)*rho(a,i,jb)*eps/(eps**2 + eta**2)
+        end do
+      end do
+    end do
+
+  end if
 
 end subroutine self_energy_correlation
diff --git a/src/MBPT/unrestricted_self_energy_correlation.f90 b/src/MBPT/unrestricted_self_energy_correlation.f90
index 1fa2f5b..d6c5f98 100644
--- a/src/MBPT/unrestricted_self_energy_correlation.f90
+++ b/src/MBPT/unrestricted_self_energy_correlation.f90
@@ -1,4 +1,4 @@
-subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC)
+subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC,EcGM)
 
 ! Compute diagonal of the correlation part of the self-energy
 
@@ -26,10 +26,12 @@ subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega
 ! Output variables
 
   double precision,intent(out)  :: SigC(nBas,nBas,nspin)
+  double precision              :: EcGM(nspin)
 
 ! Initialize 
 
   SigC(:,:,:) = 0d0
+  EcGM(:)   = 0d0
 
 !--------------!
 ! Spin-up part !
@@ -61,6 +63,17 @@ subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega
     end do
   end do
 
+  ! GM correlation energy
+
+  do i=nC(1)+1,nO(1)
+    do a=nO(1)+1,nBas-nR(1)
+      do jb=1,nSt
+        eps = e(a,1) - e(i,1) + Omega(jb)
+        EcGM(1) = EcGM(1) - rho(a,i,jb,1)**2*eps/(eps**2 + eta**2)
+      end do
+    end do
+  end do
+
 !----------------!
 ! Spin-down part !
 !----------------!
@@ -91,4 +104,15 @@ subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega
     end do
   end do
 
+  ! GM correlation energy
+
+  do i=nC(2)+1,nO(2)
+    do a=nO(2)+1,nBas-nR(2)
+      do jb=1,nSt
+        eps = e(a,2) - e(i,2) + Omega(jb)
+        EcGM(2) = EcGM(2) - rho(a,i,jb,2)**2*eps/(eps**2 + eta**2)
+      end do
+    end do
+  end do
+
 end subroutine unrestricted_self_energy_correlation
diff --git a/src/eDFT/UPW92_lda_correlation_energy.f90 b/src/eDFT/UPW92_lda_correlation_energy.f90
index f661846..9a70abc 100644
--- a/src/eDFT/UPW92_lda_correlation_energy.f90
+++ b/src/eDFT/UPW92_lda_correlation_energy.f90
@@ -58,6 +58,8 @@ subroutine UPW92_lda_correlation_energy(nGrid,weight,rho,Ec)
 
     ra = max(0d0,rho(iG,1))
     rb = max(0d0,rho(iG,2))
+    r = ra + rb
+    z = (ra - rb)/r
 
 !   alpha-alpha contribution
 
@@ -74,11 +76,9 @@ subroutine UPW92_lda_correlation_energy(nGrid,weight,rho,Ec)
 
 !   alpha-beta contribution
 
-    if(ra > threshold .or. rb > threshold) then
+    if(r > threshold) then
 
-      r = ra + rb
       rs = (4d0*pi*r/3d0)**(-1d0/3d0)
-      z = (ra - rb)/r
 
       fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
       fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
diff --git a/src/eDFT/UPW92_lda_correlation_potential.f90 b/src/eDFT/UPW92_lda_correlation_potential.f90
index 9579b62..d63d9e4 100644
--- a/src/eDFT/UPW92_lda_correlation_potential.f90
+++ b/src/eDFT/UPW92_lda_correlation_potential.f90
@@ -173,6 +173,7 @@ subroutine UPW92_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
 
           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
diff --git a/src/eDFT/UVWN3_lda_correlation_energy.f90 b/src/eDFT/UVWN3_lda_correlation_energy.f90
index 161f1ad..ff9d18e 100644
--- a/src/eDFT/UVWN3_lda_correlation_energy.f90
+++ b/src/eDFT/UVWN3_lda_correlation_energy.f90
@@ -52,6 +52,8 @@ subroutine UVWN3_lda_correlation_energy(nGrid,weight,rho,Ec)
 
     ra = max(0d0,rho(iG,1))
     rb = max(0d0,rho(iG,2))
+    r = ra + rb
+    z = (ra - rb)/r
 
 !   alpha-alpha contribution
 
@@ -73,11 +75,9 @@ subroutine UVWN3_lda_correlation_energy(nGrid,weight,rho,Ec)
 
 !   alpha-beta contribution
 
-    if(ra > threshold .or. rb > threshold) then
+    if(r > 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
diff --git a/src/eDFT/UVWN3_lda_correlation_potential.f90 b/src/eDFT/UVWN3_lda_correlation_potential.f90
index 47f8cd9..ce1c93a 100644
--- a/src/eDFT/UVWN3_lda_correlation_potential.f90
+++ b/src/eDFT/UVWN3_lda_correlation_potential.f90
@@ -59,21 +59,22 @@ subroutine UVWN3_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
 
         ra = max(0d0,rho(iG,1))
         rb = max(0d0,rho(iG,2))
+        r = ra + rb
+        z = (ra - rb)/r
+
+        fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
+        fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
+       
+        d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
+       
        
 !       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
@@ -132,16 +133,9 @@ subroutine UVWN3_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
        
         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
diff --git a/src/eDFT/UVWN5_lda_correlation_energy.f90 b/src/eDFT/UVWN5_lda_correlation_energy.f90
index 8dadbc4..6ec7eaf 100644
--- a/src/eDFT/UVWN5_lda_correlation_energy.f90
+++ b/src/eDFT/UVWN5_lda_correlation_energy.f90
@@ -52,6 +52,8 @@ subroutine UVWN5_lda_correlation_energy(nGrid,weight,rho,Ec)
 
     ra = max(0d0,rho(iG,1))
     rb = max(0d0,rho(iG,2))
+    r = ra + rb
+    z = (ra - rb)/r
 
 !   alpha-alpha contribution
 
@@ -73,11 +75,9 @@ subroutine UVWN5_lda_correlation_energy(nGrid,weight,rho,Ec)
 
 !   alpha-beta contribution
 
-    if(ra > threshold .or. rb > threshold) then
+    if(r > 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
diff --git a/src/eDFT/UVWN5_lda_correlation_potential.f90 b/src/eDFT/UVWN5_lda_correlation_potential.f90
index d36e8a1..c40ea59 100644
--- a/src/eDFT/UVWN5_lda_correlation_potential.f90
+++ b/src/eDFT/UVWN5_lda_correlation_potential.f90
@@ -59,21 +59,21 @@ subroutine UVWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
 
         ra = max(0d0,rho(iG,1))
         rb = max(0d0,rho(iG,2))
+        r  = ra + rb
+        z  = (ra - rb)/r
+
+        rs = (4d0*pi*r/3d0)**(-1d0/3d0)
+        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))
        
 !       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
@@ -132,16 +132,6 @@ subroutine UVWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
        
         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
diff --git a/src/eDFT/unrestricted_gga_correlation_energy.f90 b/src/eDFT/unrestricted_gga_correlation_energy.f90
index 34072e4..bcd3edd 100644
--- a/src/eDFT/unrestricted_gga_correlation_energy.f90
+++ b/src/eDFT/unrestricted_gga_correlation_energy.f90
@@ -30,6 +30,10 @@ subroutine unrestricted_gga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,dr
 
       call ULYP_gga_correlation_energy(nGrid,weight,rho,drho,Ec)
 
+    case ('PBE')
+
+      call UPBE_gga_correlation_energy(nGrid,weight,rho,drho,Ec)
+
     case default
 
       call print_warning('!!! GGA correlation energy not available !!!')
diff --git a/src/eDFT/unrestricted_gga_correlation_potential.f90 b/src/eDFT/unrestricted_gga_correlation_potential.f90
index fce504f..f52ebe9 100644
--- a/src/eDFT/unrestricted_gga_correlation_potential.f90
+++ b/src/eDFT/unrestricted_gga_correlation_potential.f90
@@ -34,7 +34,7 @@ subroutine unrestricted_gga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBa
 
     case ('PBE')
 
-!     call UPBE_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
+      call UPBE_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
 
     case default