diff --git a/input/basis b/input/basis
index 7a2f348..b2b2293 100644
--- a/input/basis
+++ b/input/basis
@@ -1,32 +1,30 @@
-1  6
-S   9
-1         6.665000E+03           6.920000E-04
-2         1.000000E+03           5.329000E-03
-3         2.280000E+02           2.707700E-02
-4         6.471000E+01           1.017180E-01
-5         2.106000E+01           2.747400E-01
-6         7.495000E+00           4.485640E-01
-7         2.797000E+00           2.850740E-01
-8         5.215000E-01           1.520400E-02
-9         1.596000E-01          -3.191000E-03
-S   9
-1         6.665000E+03          -1.460000E-04
-2         1.000000E+03          -1.154000E-03
-3         2.280000E+02          -5.725000E-03
-4         6.471000E+01          -2.331200E-02
-5         2.106000E+01          -6.395500E-02
-6         7.495000E+00          -1.499810E-01
-7         2.797000E+00          -1.272620E-01
-8         5.215000E-01           5.445290E-01
-9         1.596000E-01           5.804960E-01
+1   6
+S   8
+  1   1469.0000000              0.0007660        
+  2    220.5000000              0.0058920        
+  3     50.2600000              0.0296710        
+  4     14.2400000              0.1091800        
+  5      4.5810000              0.2827890        
+  6      1.5800000              0.4531230        
+  7      0.5640000              0.2747740        
+  8      0.0734500              0.0097510        
+S   8
+  1   1469.0000000             -0.0001200        
+  2    220.5000000             -0.0009230        
+  3     50.2600000             -0.0046890        
+  4     14.2400000             -0.0176820        
+  5      4.5810000             -0.0489020        
+  6      1.5800000             -0.0960090        
+  7      0.5640000             -0.1363800        
+  8      0.0734500              0.5751020        
 S   1
-1         1.596000E-01           1.000000E+00
-P   4
-1         9.439000E+00           3.810900E-02
-2         2.002000E+00           2.094800E-01
-3         5.456000E-01           5.085570E-01
-4         1.517000E-01           4.688420E-01
+  1      0.0280500              1.0000000        
+P   3
+  1      1.5340000              0.0227840        
+  2      0.2749000              0.1391070        
+  3      0.0736200              0.5003750        
 P   1
-1         1.517000E-01           1.000000E+00
+  1      0.0240300              1.0000000        
 D   1
-1         5.500000E-01           1.0000000
+  1      0.1239000              1.0000000
+
diff --git a/input/dft b/input/dft
index 49eac86..3d1826c 100644
--- a/input/dft
+++ b/input/dft
@@ -19,18 +19,18 @@
 # Number of states in ensemble (nEns)
   3
 # occupation numbers of orbitals
-1 1 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  
 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 
 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 
+1 1 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.00 0.00 
+  0.00 1.00
 # Ncentered ? 0 for NO 
- 0 
+ 1 
 # Parameters for CC weight-dependent exchange functional
 0.420431 0.069097 -0.295049
 0.135075 -0.00770826 -0.028057
diff --git a/input/methods b/input/methods
index eb2074b..ba17592 100644
--- a/input/methods
+++ b/input/methods
@@ -1,7 +1,7 @@
 # RHF UHF MOM
   F   T   F
 # MP2 MP3 MP2-F12
-  T   F   F
+  F   F   F
 # CCD CCSD CCSD(T) 
   F   F    F
 # drCCD rCCD lCCD pCCD
@@ -13,7 +13,7 @@
 # G0F2 evGF2 G0F3 evGF3
   F    F     F    F
 # G0W0 evGW qsGW 
-  F    F    F
+  T    F    F
 # G0T0 evGT qsGT 
   F    F    F
 # MCMP2
diff --git a/input/molecule b/input/molecule
index 7f71f46..d8493e9 100644
--- a/input/molecule
+++ b/input/molecule
@@ -1,4 +1,4 @@
 # nAt nEla nElb nCore nRyd
- 1   5     1   0     0
+  1   1     1   0     0
 # Znuc   x            y           z
-  C     0.           0.         0.
+  Li     0.0          0.0         0.0
diff --git a/input/molecule.xyz b/input/molecule.xyz
index 8162169..c9a5a65 100644
--- a/input/molecule.xyz
+++ b/input/molecule.xyz
@@ -1,3 +1,3 @@
   1
 
- C      0.0000000000    0.0000000000    0.0000000000
+ Li     0.0000000000    0.0000000000    0.0000000000
diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90
index bd76850..5438259 100644
--- a/src/QuAcK/G0W0.f90
+++ b/src/QuAcK/G0W0.f90
@@ -112,7 +112,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
 
   if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rhox(:,:,:,ispin))
 
-  call self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & 
+  call self_energy_correlation_diag(.false.,COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & 
                                     Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC)
 
 ! Compute renormalization factor
diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90
index 263f819..d7f4fbc 100644
--- a/src/QuAcK/QuAcK.f90
+++ b/src/QuAcK/QuAcK.f90
@@ -20,7 +20,11 @@ program QuAcK
   logical                       :: doBas
 
   integer                       :: nNuc,nBas,nBasCABS
-  integer                       :: nEl(nspin),nC(nspin),nO(nspin),nV(nspin),nR(nspin)
+  integer                       :: nEl(nspin)
+  integer                       :: nC(nspin)
+  integer                       :: nO(nspin)
+  integer                       :: nV(nspin)
+  integer                       :: nR(nspin)
   integer                       :: nS(nspin)
   double precision              :: ENuc,ERHF,EUHF,Norm
   double precision              :: EcMP2(3),EcMP3,EcMP2F12(3),EcMCMP2(3),Err_EcMCMP2(3),Var_EcMCMP2(3)
@@ -401,7 +405,18 @@ program QuAcK
   if(doMP3) then
 
     call cpu_time(start_MP3)
-    call MP3(nBas,nEl,ERI_MO,eHF,ENuc,ERHF)
+    
+    if(unrestricted) then
+
+      write(*,*) 'MP3 NYI for UHF reference'
+      stop
+
+    else
+
+      call MP3(nBas,nEl,ERI_MO,eHF,ENuc,ERHF)
+
+    end if
+
     call cpu_time(end_MP3)
 
     t_MP3 = end_MP3 - start_MP3
@@ -732,9 +747,9 @@ program QuAcK
     call cpu_time(start_G0W0)
     if(unrestricted) then 
 
-      call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,       & 
-                 dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, & 
-                 nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,PHF,cHF,eHF,eG0W0)
+      call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn, &
+                 singlet_manifold,triplet_manifold,linGW,eta_GW,nBas,nC,nO,nV,nR,nS, & 
+                 ENuc,EUHF,Hc,ERI_MO_aa,ERI_MO_ab,ERI_MO_bb,PHF,cHF,eHF,eG0W0)
     else
 
       call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,       & 
diff --git a/src/QuAcK/evGW.f90 b/src/QuAcK/evGW.f90
index 625c1fa..08e8ec3 100644
--- a/src/QuAcK/evGW.f90
+++ b/src/QuAcK/evGW.f90
@@ -140,14 +140,14 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
 
     if(G0W) then
 
-      call self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & 
+      call self_energy_correlation_diag(.false.,COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & 
                                         Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC)
       call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & 
                                   Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:))
 
     else 
 
-      call self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, & 
+      call self_energy_correlation_diag(.false.,COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, & 
                                         Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC)
       call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, & 
                                   Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:))
diff --git a/src/QuAcK/self_energy_correlation_diag.f90 b/src/QuAcK/self_energy_correlation_diag.f90
index ce96ebd..c44d27f 100644
--- a/src/QuAcK/self_energy_correlation_diag.f90
+++ b/src/QuAcK/self_energy_correlation_diag.f90
@@ -1,4 +1,4 @@
-subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,EcGM,SigC)
+subroutine self_energy_correlation_diag(unrestricted,COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,EcGM,SigC)
 
 ! Compute diagonal of the correlation part of the self-energy
 
@@ -7,6 +7,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O
 
 ! Input variables
 
+  logical,intent(in)            :: unrestricted
   logical,intent(in)            :: COHSEX
   logical,intent(in)            :: SOSEX
   double precision,intent(in)   :: eta
@@ -23,7 +24,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O
 
 ! Local variables
 
-  integer                       :: i,j,a,b,p,x,jb
+  integer                       :: i,j,a,b,p,q,jb
   double precision              :: eps
   double precision,external     :: SigC_dcgw
 
@@ -34,82 +35,135 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O
 
 ! Initialize 
 
-  SigC = 0d0
+  SigC(:) = 0d0
 
-! COHSEX static approximation
+!-----------------------------
+! COHSEX static self-energy
+!-----------------------------
 
   if(COHSEX) then
 
     ! COHSEX: SEX part of the COHSEX correlation self-energy
 
-    do x=nC+1,nBas-nR
+    do p=nC+1,nBas-nR
       do i=nC+1,nO
         jb = 0
         do j=nC+1,nO
           do b=nO+1,nBas-nR
             jb = jb + 1
-            SigC(x) = SigC(x) + 4d0*rho(x,i,jb)**2/Omega(jb)
-          enddo
-        enddo
-      enddo
-    enddo
+            SigC(p) = SigC(p) + 4d0*rho(p,i,jb)**2/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 p=nC+1,nBas-nR
+    do p=nC+1,nBas-nR
+      do q=nC+1,nBas-nR
         jb = 0
         do j=nC+1,nO
           do b=nO+1,nBas-nR
             jb = jb + 1
-            SigC(x) = SigC(x) - 2d0*rho(x,p,jb)**2/Omega(jb)
-          enddo
-        enddo
-      enddo
-    enddo
+            SigC(p) = SigC(p) - 2d0*rho(p,q,jb)**2/Omega(jb)
+          end do
+        end do
+      end do
+    end do
 
     ! GM correlation energy
 
     EcGM = 0d0
     do i=nC+1,nO
       EcGM = EcGM + 0.5d0*SigC(i)
-    enddo
- 
-  else
- 
-    ! Occupied part of the correlation self-energy
- 
-    do x=nC+1,nBas-nR
+    end do
+
+!-----------------------------
+! SOSEX self-energy *BUG*
+!-----------------------------
+
+  elseif(SOSEX) then
+
+    ! SOSEX: occupied part of the correlation self-energy
+
+    do p=nC+1,nBas-nR
       do i=nC+1,nO
         jb = 0
         do j=nC+1,nO
           do b=nO+1,nBas-nR
             jb = jb + 1
-            eps = e(x) - e(i) + Omega(jb)
-!           SigC(x) = SigC(x) + 4d0*rho(x,i,jb)**2/(eps + eps*sqrt(1d0 + rho(x,i,jb)**2/eps**2))
-            SigC(x) = SigC(x) + 2d0*rho(x,i,jb)**2*eps/(eps**2 + eta**2)
-!           SigC(x) = SigC(x) + 2d0*SigC_dcgw(eps,rho(x,i,jb))
-          enddo
-        enddo
-      enddo
-    enddo
- 
-    ! Virtual part of the correlation self-energy
- 
-    do x=nC+1,nBas-nR
+            eps = e(p) - e(i) + Omega(jb)
+            SigC(p) = SigC(p) - rho(p,i,jb)*rhox(p,i,jb)*eps/(eps**2 + eta**2)
+          end do
+        end do
+      end do
+    end do
+
+    ! SOSEX: virtual part of the correlation self-energy
+
+    do p=nC+1,nBas-nR
       do a=nO+1,nBas-nR
         jb = 0
         do j=nC+1,nO
           do b=nO+1,nBas-nR
             jb = jb + 1
-            eps = e(x) - e(a) - Omega(jb)
-!           SigC(x) = SigC(x) + 4d0*rho(x,a,jb)**2/(eps + eps*sqrt(1d0 + 4d0*rho(x,a,jb)**2/eps**2))
-            SigC(x) = SigC(x) + 2d0*rho(x,a,jb)**2*eps/(eps**2 + eta**2)
-!           SigC(x) = SigC(x) + 2d0*SigC_dcgw(eps,rho(x,a,jb))
-          enddo
-        enddo
-      enddo
-    enddo
+            eps = e(p) - e(a) - Omega(jb)
+            SigC(p) = SigC(p) - rho(p,a,jb)*rhox(p,a,jb)*eps/(eps**2 + eta**2)
+          end do
+        end do
+      end do
+    end do
+
+    ! GM correlation energy
+
+    do i=nC+1,nO
+      do a=nO+1,nBas-nR
+        jb = 0
+        do j=nC+1,nO
+          do b=nO+1,nBas-nR
+            jb = jb + 1
+            eps = e(a) - e(i) + Omega(jb)
+            EcGM = EcGM + rho(a,i,jb)*rhox(a,i,jb)*eps/(eps**2 + eta**2)
+          end do
+        end do
+      end do
+    end do
+
+!-----------------------------
+! GW self-energy
+!-----------------------------
+
+  else
+ 
+    ! Occupied part of the correlation self-energy
+ 
+    do p=nC+1,nBas-nR
+      do i=nC+1,nO
+        jb = 0
+        do j=nC+1,nO
+          do b=nO+1,nBas-nR
+            jb = jb + 1
+            eps = e(p) - e(i) + Omega(jb)
+            SigC(p) = SigC(p) + 2d0*rho(p,i,jb)**2*eps/(eps**2 + eta**2)
+          end do
+        end do
+      end do
+    end do
+ 
+    ! Virtual part of the correlation self-energy
+ 
+    do p=nC+1,nBas-nR
+      do a=nO+1,nBas-nR
+        jb = 0
+        do j=nC+1,nO
+          do b=nO+1,nBas-nR
+            jb = jb + 1
+            eps = e(p) - e(a) - Omega(jb)
+            SigC(p) = SigC(p) + 2d0*rho(p,a,jb)**2*eps/(eps**2 + eta**2)
+          end do
+        end do
+      end do
+    end do
 
     ! GM correlation energy
 
@@ -122,60 +176,20 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O
             jb = jb + 1
             eps = e(a) - e(i) + Omega(jb)
             EcGM = EcGM - 2d0*rho(a,i,jb)*rho(a,i,jb)*eps/(eps**2 + eta**2)
-          enddo
-        enddo
-      enddo
-    enddo
+          end do
+        end do
+      end do
+    end do
 
-    if(SOSEX) then
+  end if
 
-      ! SOSEX: occupied part of the correlation self-energy
+! Unrestricted reference
 
-      do x=nC+1,nBas-nR
-        do i=nC+1,nO
-          jb = 0
-          do j=nC+1,nO
-            do b=nO+1,nBas-nR
-              jb = jb + 1
-              eps = e(x) - e(i) + Omega(jb)
-              SigC(x) = SigC(x) - rho(x,i,jb)*rhox(x,i,jb)*eps/(eps**2 + eta**2)
-            enddo
-          enddo
-        enddo
-      enddo
+  if(unrestricted) then 
+ 
+    SigC(:) = 0.5d0*SigC(:)
+    EcGM = 0.5d0*EcGM
 
-      ! SOSEX: virtual part of the correlation self-energy
-
-      do x=nC+1,nBas-nR
-        do a=nO+1,nBas-nR
-          jb = 0
-          do j=nC+1,nO
-            do b=nO+1,nBas-nR
-              jb = jb + 1
-              eps = e(x) - e(a) - Omega(jb)
-              SigC(x) = SigC(x) - rho(x,a,jb)*rhox(x,a,jb)*eps/(eps**2 + eta**2)
-            enddo
-          enddo
-        enddo
-      enddo
-
-      ! GM correlation energy
-
-      do i=nC+1,nO
-        do a=nO+1,nBas-nR
-          jb = 0
-          do j=nC+1,nO
-            do b=nO+1,nBas-nR
-              jb = jb + 1
-              eps = e(a) - e(i) + Omega(jb)
-              EcGM = EcGM + rho(a,i,jb)*rhox(a,i,jb)*eps/(eps**2 + eta**2)
-            enddo
-          enddo
-        enddo
-      enddo
-
-    endif
-
-  endif
+  end if
 
 end subroutine self_energy_correlation_diag
diff --git a/src/eDFT/print_unrestricted_individual_energy.f90 b/src/eDFT/print_unrestricted_individual_energy.f90
index 4c3db11..f5c57e0 100644
--- a/src/eDFT/print_unrestricted_individual_energy.f90
+++ b/src/eDFT/print_unrestricted_individual_energy.f90
@@ -147,12 +147,12 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,
     write(*,'(A44,      F16.10,A3)') '        c  ensemble derivative  : ',OmcDD(2), ' au'
     write(*,'(A44,      F16.10,A3)') '       xc  ensemble derivative  : ',OmxcDD(2),' au'
     write(*,*)
-    write(*,'(A43,F16.10,A4)') ' Electronic Affinity  1 -> 3:',-1.d0*(Omaux(3)+OmxcDD(3)),' au'
+    write(*,'(A43,F16.10,A4)') ' Electronic Affinity  1 -> 3:',-(Omaux(3)+OmxcDD(3)),' au'
     write(*,*)
-    write(*,'(A44,      F16.10,A3)') ' auxiliary energy contribution  : ',-1.d0*Omaux(3), ' au'
-    write(*,'(A44,      F16.10,A3)') '        x  ensemble derivative  : ',-1.d0*OmxDD(3), ' au'
-    write(*,'(A44,      F16.10,A3)') '        c  ensemble derivative  : ',-1.d0*OmcDD(3), ' au'
-    write(*,'(A44,      F16.10,A3)') '       xc  ensemble derivative  : ',-1.d0*OmxcDD(3),' au'
+    write(*,'(A44,      F16.10,A3)') ' auxiliary energy contribution  : ',-Omaux(3), ' au'
+    write(*,'(A44,      F16.10,A3)') '        x  ensemble derivative  : ',-OmxDD(3), ' au'
+    write(*,'(A44,      F16.10,A3)') '        c  ensemble derivative  : ',-OmcDD(3), ' au'
+    write(*,'(A44,      F16.10,A3)') '       xc  ensemble derivative  : ',-OmxcDD(3),' au'
     write(*,*)
     write(*,'(A43,F16.10,A4)') ' Fundamental Gap :',Omaux(2)+OmxcDD(2)+(Omaux(3)+OmxcDD(3)),' au'
     write(*,*)
@@ -172,12 +172,12 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,
     write(*,'(A44,      F16.10,A3)') '        c  ensemble derivative  : ',OmcDD(2)*HaToeV, ' eV'
     write(*,'(A44,      F16.10,A3)') '       xc  ensemble derivative  : ',OmxcDD(2)*HaToeV,' eV'
     write(*,*)
-    write(*,'(A40,F16.10,A3)') ' Electronic Affinity  1 -> 3:',-1.d0*(Omaux(3)+OmxcDD(3))*HaToeV,' eV'
+    write(*,'(A40,F16.10,A3)') ' Electronic Affinity  1 -> 3:',-(Omaux(3)+OmxcDD(3))*HaToeV,' eV'
     write(*,*)
-    write(*,'(A44,      F16.10,A3)') ' auxiliary energy contribution  : ',-1.d0*Omaux(3)*HaToeV, ' eV'
-    write(*,'(A44,      F16.10,A3)') '        x  ensemble derivative  : ',-1.d0*OmxDD(3)*HaToeV, ' eV'
-    write(*,'(A44,      F16.10,A3)') '        c  ensemble derivative  : ',-1.d0*OmcDD(3)*HaToeV, ' eV'
-    write(*,'(A44,      F16.10,A3)') '       xc  ensemble derivative  : ',-1.d0*OmxcDD(3)*HaToeV,' eV'
+    write(*,'(A44,      F16.10,A3)') ' auxiliary energy contribution  : ',-Omaux(3)*HaToeV, ' eV'
+    write(*,'(A44,      F16.10,A3)') '        x  ensemble derivative  : ',-OmxDD(3)*HaToeV, ' eV'
+    write(*,'(A44,      F16.10,A3)') '        c  ensemble derivative  : ',-OmcDD(3)*HaToeV, ' eV'
+    write(*,'(A44,      F16.10,A3)') '       xc  ensemble derivative  : ',-OmxcDD(3)*HaToeV,' eV'
     write(*,*)
     write(*,'(A43,F16.10,A4)') ' Fundamental Gap :',(Omaux(2)+OmxcDD(2)+(Omaux(3)+OmxcDD(3)))*HaToeV,' eV'
     write(*,*)
@@ -208,15 +208,15 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,
     write(*,'(A44,      F16.10,A3)') '     c  ensemble derivative     : ',OmcDD(2), ' au'
     write(*,'(A44,      F16.10,A3)') '    xc  ensemble derivative     : ',OmxcDD(2),' au'
     write(*,*)
-    write(*,'(A43,F16.10,A4)') ' Electronic Affinity  1 -> 3:',-1.d0*Om(3),    ' au'
+    write(*,'(A43,F16.10,A4)') ' Electronic Affinity  1 -> 3:',-Om(3),    ' au'
     write(*,*)
-    write(*,'(A44,      F16.10,A3)') '     x  energy contribution     : ',-1.d0*Omx(3),   ' au'
-    write(*,'(A44,      F16.10,A3)') '     c  energy contribution     : ',-1.d0*Omc(3),   ' au'
-    write(*,'(A44,      F16.10,A3)') '    xc  energy contribution     : ',-1.d0*Omxc(3),  ' au'
+    write(*,'(A44,      F16.10,A3)') '     x  energy contribution     : ',-Omx(3),   ' au'
+    write(*,'(A44,      F16.10,A3)') '     c  energy contribution     : ',-Omc(3),   ' au'
+    write(*,'(A44,      F16.10,A3)') '    xc  energy contribution     : ',-Omxc(3),  ' au'
     write(*,*)
-    write(*,'(A44,      F16.10,A3)') '     x  ensemble derivative     : ',-1.d0*OmxDD(3), ' au'
-    write(*,'(A44,      F16.10,A3)') '     c  ensemble derivative     : ',-1.d0*OmcDD(3), ' au'
-    write(*,'(A44,      F16.10,A3)') '    xc  ensemble derivative     : ',-1.d0*OmxcDD(3),' au'
+    write(*,'(A44,      F16.10,A3)') '     x  ensemble derivative     : ',-OmxDD(3), ' au'
+    write(*,'(A44,      F16.10,A3)') '     c  ensemble derivative     : ',-OmcDD(3), ' au'
+    write(*,'(A44,      F16.10,A3)') '    xc  ensemble derivative     : ',-OmxcDD(3),' au'
     write(*,*)
     write(*,'(A43,F16.10,A4)') ' Fundamental Gap :',Om(2)+Om(3),    ' au'
     write(*,*)
@@ -241,15 +241,15 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,
     write(*,'(A44,      F16.10,A3)') '     c  ensemble derivative     : ',OmcDD(2)*HaToeV, ' eV'
     write(*,'(A44,      F16.10,A3)') '    xc  ensemble derivative     : ',OmxcDD(2)*HaToeV,' eV'
     write(*,*)
-    write(*,'(A43,F16.10,A4)') ' Electronic Affinity 1 -> 3:',-1.d0*Om(3)*HaToeV,    ' eV'
+    write(*,'(A43,F16.10,A4)') ' Electronic Affinity 1 -> 3:',-Om(3)*HaToeV,    ' eV'
     write(*,*)
-    write(*,'(A44,      F16.10,A3)') '     x  energy contribution     : ',-1.d0*Omx(3)*HaToeV,   ' eV'
-    write(*,'(A44,      F16.10,A3)') '     c  energy contribution     : ',-1.d0*Omc(3)*HaToeV,   ' eV'
-    write(*,'(A44,      F16.10,A3)') '    xc  energy contribution     : ',-1.d0*Omxc(3)*HaToeV,  ' eV'
+    write(*,'(A44,      F16.10,A3)') '     x  energy contribution     : ',-Omx(3)*HaToeV,   ' eV'
+    write(*,'(A44,      F16.10,A3)') '     c  energy contribution     : ',-Omc(3)*HaToeV,   ' eV'
+    write(*,'(A44,      F16.10,A3)') '    xc  energy contribution     : ',-Omxc(3)*HaToeV,  ' eV'
     write(*,*)
-    write(*,'(A44,      F16.10,A3)') '     x  ensemble derivative     : ',-1.d0*OmxDD(3)*HaToeV, ' eV'
-    write(*,'(A44,      F16.10,A3)') '     c  ensemble derivative     : ',-1.d0*OmcDD(3)*HaToeV, ' eV'
-    write(*,'(A44,      F16.10,A3)') '    xc  ensemble derivative     : ',-1.d0*OmxcDD(3)*HaToeV,' eV'
+    write(*,'(A44,      F16.10,A3)') '     x  ensemble derivative     : ',-OmxDD(3)*HaToeV, ' eV'
+    write(*,'(A44,      F16.10,A3)') '     c  ensemble derivative     : ',-OmcDD(3)*HaToeV, ' eV'
+    write(*,'(A44,      F16.10,A3)') '    xc  ensemble derivative     : ',-OmxcDD(3)*HaToeV,' eV'
     write(*,*)
     write(*,'(A43,F16.10,A4)') ' Fundamental Gap :',(Om(2)+Om(3))*HaToeV,    ' eV'
     write(*,*)
diff --git a/src/eDFT/read_options.f90 b/src/eDFT/read_options.f90
index 4eb7f62..8afe913 100644
--- a/src/eDFT/read_options.f90
+++ b/src/eDFT/read_options.f90
@@ -138,7 +138,7 @@ subroutine read_options(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_
   read(1,*)
   read(1,*) doNcentered
 
-  if (doNcentered==0) then
+  if (doNcentered == 0) then
 
     wEns(1) = 1d0 - wEns(2) - wEns(3) 
 
diff --git a/src/utils/read_basis.f90 b/src/utils/read_basis.f90
index 35add6e..75a9f49 100644
--- a/src/utils/read_basis.f90
+++ b/src/utils/read_basis.f90
@@ -8,7 +8,7 @@ subroutine read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KSh
 
 ! Input variables
 
-  integer,intent(in)            :: nNuc,nO
+  integer,intent(in)            :: nNuc,nO(nspin)
   double precision,intent(in)   :: rNuc(nNuc,ncart)
 
 ! Local variables
@@ -23,7 +23,7 @@ subroutine read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KSh
   double precision,intent(out)  :: CenterShell(maxShell,ncart)
   integer,intent(out)           :: TotAngMomShell(maxShell),KShell(maxShell)
   double precision,intent(out)  :: DShell(maxShell,maxK),ExpShell(maxShell,maxK)
-  integer,intent(out)           :: nV
+  integer,intent(out)           :: nV(nspin)
 
   integer,intent(out)           :: max_ang_mom(nNuc)
   double precision,intent(out)  :: min_exponent(nNuc,maxL+1)
@@ -180,6 +180,6 @@ subroutine read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KSh
 
 ! Number of virtual orbitals
 
-  nV = nBas - nO
+  nV(:) = nBas - nO(:)
 
 end subroutine read_basis