diff --git a/src/LR/ppGLR.f90 b/src/LR/ppGLR.f90
index 536015b..362fcd7 100644
--- a/src/LR/ppGLR.f90
+++ b/src/LR/ppGLR.f90
@@ -1,4 +1,4 @@
-subroutine ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
+subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
 
   ! 
   ! Solve the pp-RPA linear eigenvalue problem
@@ -28,98 +28,103 @@ subroutine ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
   implicit none
   include 'parameters.h'
 
-  logical,          intent(in)  :: TDA
-  integer,          intent(in)  :: nOO, nVV
-  double precision, intent(in)  :: Bpp(nVV,nOO), Cpp(nVV,nVV), Dpp(nOO,nOO)
-  double precision, intent(out) :: Om1(nVV), X1(nVV,nVV), Y1(nOO,nVV)
-  double precision, intent(out) :: Om2(nOO), X2(nVV,nOO), Y2(nOO,nOO)
-  double precision, intent(out) :: EcRPA
+  logical,intent(in)           :: TDA
+  integer,intent(in)           :: nOO,nVV
+  double precision,intent(in)  :: Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)
+  double precision,intent(out) :: Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV)
+  double precision,intent(out) :: Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO)
+  double precision,intent(out) :: EcRPA
   
-  logical                       :: imp_bio, verbose
-  integer                       :: i, j, N
-  double precision              :: EcRPA1, EcRPA2
-  double precision              :: thr_d, thr_nd, thr_deg
-  double precision,allocatable  :: M(:,:), Z(:,:), Om(:)
+  logical                       :: imp_bio,verbose
+  integer                       :: i,j
+  integer                       :: Npp
+  double precision              :: EcRPA1,EcRPA2
+  double precision              :: thr_d,thr_nd,thr_deg
+  double precision,allocatable  :: M(:,:),Z(:,:),Om(:)
 
-  double precision, external    :: trace_matrix
+  double precision,external     :: trace_matrix
 
+  nPP = nOO + nVV
 
+  allocate(M(nPP,nPP),Z(nPP,nPP),Om(nPP))
 
-  N = nOO + nVV
-
-  allocate(M(N,N), Z(N,N), Om(N))
+  ! Hermitian case for TDA
 
   if(TDA) then
 
     X1(:,:) = +Cpp(:,:)
     Y1(:,:) = 0d0
-    if(nVV > 0) call diagonalize_matrix(nVV, X1, Om1)
+    if(nVV > 0) call diagonalize_matrix(nVV,X1,Om1)
 
     X2(:,:) = 0d0
     Y2(:,:) = -Dpp(:,:)
-    if(nOO > 0) call diagonalize_matrix(nOO, Y2, Om2)
+    if(nOO > 0) call diagonalize_matrix(nOO,Y2,Om2)
 
   else
 
     ! Diagonal blocks 
-    M(    1:nVV    ,    1:nVV)     = + Cpp(1:nVV,1:nVV)
-    M(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - Dpp(1:nOO,1:nOO)
+
+    M(    1:nVV,    1:nVV) = + Cpp(1:nVV,1:nVV)
+    M(nVV+1:nPP,nVV+1:nPP) = - Dpp(1:nOO,1:nOO)
 
     ! Off-diagonal blocks
-    M(    1:nVV    ,nVV+1:nOO+nVV) = -           Bpp(1:nVV,1:nOO)
-    M(nVV+1:nOO+nVV,    1:nVV)     = + transpose(Bpp(1:nVV,1:nOO))
 
-    if((nOO .eq. 0) .or. (nVV .eq. 0)) then
+    M(    1:nVV,nVV+1:nPP) = +           Bpp(1:nVV,1:nOO)
+    M(nVV+1:nPP,    1:nVV) = - transpose(Bpp(1:nVV,1:nOO))
 
-      ! Diagonalize the p-p matrix
-      if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV, M, Om, Z)
-      ! Split the various quantities in p-p and h-h parts
-      call sort_ppRPA(nOO, nVV, Om, Z, Om1, X1, Y1, Om2, X2, Y2)
+!   if((nOO .eq. 0) .or. (nVV .eq. 0)) then
 
-    else
+      ! Diagonalize the pp matrix
 
-        thr_d   = 1d-6   ! to determine if     diagonal elements of L.T x R are close enouph to 1
-        thr_nd  = 1d-6   ! to determine if non-diagonal elements of L.T x R are close enouph to 1
-        thr_deg = 1d-8   ! to determine if two eigenvectors are degenerate or not
-        imp_bio = .True. ! impose bi-orthogonality
-        verbose = .False.
-        call diagonalize_nonsym_matrix(N, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose)
-    
-        do i = 1, nOO
-          Om2(i) = Om(i)
-          do j = 1, nVV
-            X2(j,i) = Z(j,i)
-          enddo
-          do j = 1, nOO
-            Y2(j,i) = Z(nVV+j,i)
-          enddo
-        enddo
-    
-        do i = 1, nVV
-          Om1(i) = Om(nOO+i)
-          do j = 1, nVV
-            X1(j,i) = M(j,nOO+i)
-          enddo
-          do j = 1, nOO
-            Y1(j,i) = M(nVV+j,nOO+i)
-          enddo
-        enddo
+      if(nPP > 0) call diagonalize_general_matrix(nPP,M,Om,Z)
 
-    endif
+      ! Split the various quantities in ee and hh parts
+
+      call sort_ppRPA(nOO,nVV,nPP,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
+
+!   else
+
+!       thr_d   = 1d-6   ! to determine if     diagonal elements of L.T x R are close enouph to 1
+!       thr_nd  = 1d-6   ! to determine if non-diagonal elements of L.T x R are close enouph to 1
+!       thr_deg = 1d-8   ! to determine if two eigenvectors are degenerate or not
+!       imp_bio = .True. ! impose bi-orthogonality
+!       verbose = .False.
+!       call diagonalize_nonsym_matrix(N, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose)
+!   
+!       do i = 1, nOO
+!         Om2(i) = Om(i)
+!         do j = 1, nVV
+!           X2(j,i) = Z(j,i)
+!         enddo
+!         do j = 1, nOO
+!           Y2(j,i) = Z(nVV+j,i)
+!         enddo
+!       enddo
+!   
+!       do i = 1, nVV
+!         Om1(i) = Om(nOO+i)
+!         do j = 1, nVV
+!           X1(j,i) = M(j,nOO+i)
+!         enddo
+!         do j = 1, nOO
+!           Y1(j,i) = M(nVV+j,nOO+i)
+!         enddo
+!       enddo
+
+!   endif
 
   end if
 
   ! Compute the RPA correlation energy
-  EcRPA = 0.5d0 * (sum(Om1) - sum(Om2) - trace_matrix(nVV, Cpp) - trace_matrix(nOO, Dpp))
-  EcRPA1 = +sum(Om1) - trace_matrix(nVV, Cpp)
-  EcRPA2 = -sum(Om2) - trace_matrix(nOO, Dpp)
+
+  EcRPA = 0.5d0 * (sum(Om1) - sum(Om2) - trace_matrix(nVV,Cpp) - trace_matrix(nOO,Dpp))
+  EcRPA1 = +sum(Om1) - trace_matrix(nVV,Cpp)
+  EcRPA2 = -sum(Om2) - trace_matrix(nOO,Dpp)
 
   if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) then
     print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!'
   endif
 
-  deallocate(M, Z, Om)
+  deallocate(M,Z,Om)
 
 end subroutine 
-
-
diff --git a/src/LR/ppLR.f90 b/src/LR/ppLR.f90
index 6d62908..362fcd7 100644
--- a/src/LR/ppLR.f90
+++ b/src/LR/ppLR.f90
@@ -28,54 +28,59 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
   implicit none
   include 'parameters.h'
 
-  logical,          intent(in)  :: TDA
-  integer,          intent(in)  :: nOO, nVV
-  double precision, intent(in)  :: Bpp(nVV,nOO), Cpp(nVV,nVV), Dpp(nOO,nOO)
-  double precision, intent(out) :: Om1(nVV), X1(nVV,nVV), Y1(nOO,nVV)
-  double precision, intent(out) :: Om2(nOO), X2(nVV,nOO), Y2(nOO,nOO)
-  double precision, intent(out) :: EcRPA
+  logical,intent(in)           :: TDA
+  integer,intent(in)           :: nOO,nVV
+  double precision,intent(in)  :: Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)
+  double precision,intent(out) :: Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV)
+  double precision,intent(out) :: Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO)
+  double precision,intent(out) :: EcRPA
   
-  logical                       :: imp_bio, verbose
-  integer                       :: i, j, N
-  double precision              :: EcRPA1, EcRPA2
-  double precision              :: thr_d, thr_nd, thr_deg
-  double precision,allocatable  :: M(:,:), Z(:,:), Om(:)
+  logical                       :: imp_bio,verbose
+  integer                       :: i,j
+  integer                       :: Npp
+  double precision              :: EcRPA1,EcRPA2
+  double precision              :: thr_d,thr_nd,thr_deg
+  double precision,allocatable  :: M(:,:),Z(:,:),Om(:)
 
-  double precision, external    :: trace_matrix
+  double precision,external     :: trace_matrix
 
-  N = nOO + nVV
+  nPP = nOO + nVV
 
-  allocate(M(N,N), Z(N,N), Om(N))
+  allocate(M(nPP,nPP),Z(nPP,nPP),Om(nPP))
+
+  ! Hermitian case for TDA
 
   if(TDA) then
 
     X1(:,:) = +Cpp(:,:)
     Y1(:,:) = 0d0
-    if(nVV > 0) call diagonalize_matrix(nVV, X1, Om1)
+    if(nVV > 0) call diagonalize_matrix(nVV,X1,Om1)
 
     X2(:,:) = 0d0
     Y2(:,:) = -Dpp(:,:)
-    if(nOO > 0) call diagonalize_matrix(nOO, Y2, Om2)
+    if(nOO > 0) call diagonalize_matrix(nOO,Y2,Om2)
 
   else
 
     ! Diagonal blocks 
-    M(    1:nVV    ,    1:nVV)     = + Cpp(1:nVV,1:nVV)
-    M(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - Dpp(1:nOO,1:nOO)
+
+    M(    1:nVV,    1:nVV) = + Cpp(1:nVV,1:nVV)
+    M(nVV+1:nPP,nVV+1:nPP) = - Dpp(1:nOO,1:nOO)
 
     ! Off-diagonal blocks
-    M(    1:nVV    ,nVV+1:nOO+nVV) = -           Bpp(1:nVV,1:nOO)
-    M(nVV+1:nOO+nVV,    1:nVV)     = + transpose(Bpp(1:nVV,1:nOO))
+
+    M(    1:nVV,nVV+1:nPP) = +           Bpp(1:nVV,1:nOO)
+    M(nVV+1:nPP,    1:nVV) = - transpose(Bpp(1:nVV,1:nOO))
 
 !   if((nOO .eq. 0) .or. (nVV .eq. 0)) then
 
       ! Diagonalize the pp matrix
 
-      if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Om,Z)
+      if(nPP > 0) call diagonalize_general_matrix(nPP,M,Om,Z)
 
-      ! Split the various quantities in p-p and h-h parts
+      ! Split the various quantities in ee and hh parts
 
-      call sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
+      call sort_ppRPA(nOO,nVV,nPP,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
 
 !   else
 
@@ -112,9 +117,9 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
 
   ! Compute the RPA correlation energy
 
-  EcRPA = 0.5d0 * (sum(Om1) - sum(Om2) - trace_matrix(nVV, Cpp) - trace_matrix(nOO, Dpp))
-  EcRPA1 = +sum(Om1) - trace_matrix(nVV, Cpp)
-  EcRPA2 = -sum(Om2) - trace_matrix(nOO, Dpp)
+  EcRPA = 0.5d0 * (sum(Om1) - sum(Om2) - trace_matrix(nVV,Cpp) - trace_matrix(nOO,Dpp))
+  EcRPA1 = +sum(Om1) - trace_matrix(nVV,Cpp)
+  EcRPA2 = -sum(Om2) - trace_matrix(nOO,Dpp)
 
   if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) then
     print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!'
@@ -123,5 +128,3 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
   deallocate(M,Z,Om)
 
 end subroutine 
-
-
diff --git a/src/LR/ppULR.f90 b/src/LR/ppULR.f90
index 49afb79..f9e40af 100644
--- a/src/LR/ppULR.f90
+++ b/src/LR/ppULR.f90
@@ -81,8 +81,7 @@ subroutine ppULR(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,nHaa,nHab,nHbb,nH
 
 ! Split the various quantities in p-p and h-h parts
 
-  call sort_ppRPA(nHt,nPt,Om(:),Z(:,:),Om1(:),X1(:,:),Y1(:,:),Om2(:),X2(:,:),&
-                  Y2(:,:))
+  call sort_ppRPA(nHt,nPt,nHt+nPt,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
  
 ! Compute the RPA correlation energy
 
diff --git a/src/RPA/sort_ppRPA.f90 b/src/RPA/sort_ppRPA.f90
index 24c6f59..1c1c80a 100644
--- a/src/RPA/sort_ppRPA.f90
+++ b/src/RPA/sort_ppRPA.f90
@@ -1,4 +1,4 @@
-subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
+subroutine sort_ppRPA(nOO,nVV,nPP,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
 
 ! Compute the metric matrix for pp-RPA
 
@@ -9,14 +9,13 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
 
   integer,intent(in)            :: nOO
   integer,intent(in)            :: nVV
-  double precision,intent(in)   :: Om(nOO+nVV)
-  double precision,intent(in)   :: Z(nOO+nVV,nOO+nVV)
+  integer,intent(in)            :: nPP
+  double precision,intent(in)   :: Om(nPP)
+  double precision,intent(in)   :: Z(nPP,nPP)
   
 ! Local variables
 
   integer                       :: pq,ab,ij
-! integer                       :: deg1,ab_start,ab_end
-! integer                       :: deg2,ij_start,ij_end
   double precision,allocatable  :: M(:,:)
   double precision,allocatable  :: Z1(:,:)
   double precision,allocatable  :: Z2(:,:)
@@ -42,7 +41,7 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
 
 ! Memory allocation
 
-  allocate(M(nOO+nVV,nOO+nVV),Z1(nOO+nVV,nVV),Z2(nOO+nVV,nOO),order1(nVV),order2(nOO))
+  allocate(M(nPP,nPP),Z1(nPP,nVV),Z2(nPP,nOO),order1(nVV),order2(nOO))
 
 ! Initializatiom
 
@@ -71,19 +70,19 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
   ab = 0
   ij = 0
 
-  do pq=1,nOO+nVV
+  do pq=1,nPP
 
     if(Om(pq) > 0d0) then 
 
       ab = ab + 1
       Om1(ab) = Om(pq)
-      Z1(1:nOO+nVV,ab) = Z(1:nOO+nVV,pq)
+      Z1(1:nPP,ab) = Z(1:nPP,pq)
 
     else
 
       ij = ij + 1
       Om2(ij) = Om(pq)
-      Z2(1:nOO+nVV,ij) = Z(1:nOO+nVV,pq)
+      Z2(1:nPP,ij) = Z(1:nPP,pq)
 
     end if
 
@@ -99,7 +98,7 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
     end do
 
     call quick_sort(Om1,order1,nVV)
-    call set_order(Z1,order1,nOO+nVV,nVV)
+    call set_order(Z1,order1,nPP,nVV)
 
   end if
 
@@ -110,17 +109,17 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
     end do
 
     call quick_sort(Om2,order2,nOO)
-    call set_order(Z2,order2,nOO+nVV,nOO)
+    call set_order(Z2,order2,nPP,nOO)
 
   end if
  
   allocate(S1(nVV,nVV),S2(nOO,nOO),O1(nVV,nVV),O2(nOO,nOO))
-  allocate(tmp1(nOO+nVV,nVV),tmp2(nOO+nVV,nOO))
+  allocate(tmp1(nPP,nVV),tmp2(nPP,nOO))
 
-  if(nVV > 0) call dgemm ('N', 'N', nOO+nVV, nVV, nOO+nVV, 1d0,  M, nOO+nVV, Z1, nOO+nVV, 0d0, tmp1, nOO+nVV)
-  if(nVV > 0) call dgemm ('T', 'N', nVV    , nVV, nOO+nVV, 1d0, Z1, nOO+nVV, tmp1, nOO+nVV, 0d0, S1, nVV)
-  if(nOO > 0) call dgemm ('N', 'N', nOO+nVV, nOO, nOO+nVV, 1d0,  M, nOO+nVV, -1d0*Z2, nOO+nVV, 0d0, tmp2, nOO+nVV)
-  if(nOO > 0) call dgemm ('T', 'N', nOO    , nOO, nOO+nVV, 1d0, Z2, nOO+nVV, tmp2, nOO+nVV, 0d0, S2, nOO)
+  if(nVV > 0) call dgemm ('N','N',nPP,nVV,nPP,1d0, M,nPP,Z1,nPP,0d0,tmp1,nPP)
+  if(nVV > 0) call dgemm ('T','N',nVV,nVV,nPP,1d0,Z1,nPP,tmp1,nPP,0d0,S1,nVV)
+  if(nOO > 0) call dgemm ('N','N',nPP,nOO,nPP,1d0, M,nPP,-1d0*Z2,nPP,0d0,tmp2, nPP)
+  if(nOO > 0) call dgemm ('T','N',nOO,nOO,nPP,1d0,Z2,nPP,tmp2,nPP,0d0,S2,nOO)
 
 ! S1 = + matmul(transpose(Z1),matmul(M,Z1))
 ! S2 = - matmul(transpose(Z2),matmul(M,Z2))
@@ -128,9 +127,9 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
   if(nVV > 0) call orthogonalize_matrix(1,nVV,S1,O1)
   if(nOO > 0) call orthogonalize_matrix(1,nOO,S2,O2)
 
-  if(nVV > 0) call dgemm ('N', 'N', nOO+nVV,nVV,nVV, 1d0, Z1, nOO+nVV, O1, nVV,0d0, tmp1, nOO+nVV)
+  if(nVV > 0) call dgemm ('N','N',nPP,nVV,nVV,1d0,Z1,nPP,O1,nVV,0d0,tmp1,nPP)
   Z1 = tmp1
-  if(nOO > 0) call dgemm ('N', 'N', nOO+nVV,nOO,nOO, 1d0, Z2, nOO+nVV, O2, nOO,0d0, tmp2, nOO+nVV)
+  if(nOO > 0) call dgemm ('N','N',nPP,nOO,nOO,1d0,Z2,nPP,O2,nOO,0d0,tmp2,nPP)
   Z2 = tmp2
 
 ! Z1 = matmul(Z1,O1)
@@ -138,11 +137,11 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
 
 ! Define submatrices X1, Y1, X2, & Y2
 
-  X1(1:nVV,1:nVV) = + Z1(    1:    nVV,1:nVV)
-  Y1(1:nOO,1:nVV) = - Z1(nVV+1:nOO+nVV,1:nVV)
+  X1(1:nVV,1:nVV) = Z1(    1:    nVV,1:nVV)
+  Y1(1:nOO,1:nVV) = Z1(nVV+1:nPP,1:nVV)
 
-  X2(1:nVV,1:nOO) = + Z2(    1:    nVV,1:nOO)
-  Y2(1:nOO,1:nOO) = - Z2(nVV+1:nOO+nVV,1:nOO)
+  X2(1:nVV,1:nOO) = Z2(    1:    nVV,1:nOO)
+  Y2(1:nOO,1:nOO) = Z2(nVV+1:nPP,1:nOO)
 
 ! call matout(nVV,nVV,X1)
 ! call matout(nOO,nVV,Y1)
@@ -150,4 +149,11 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
 ! call matout(nVV,nOO,X2)
 ! call matout(nOO,nOO,Y2)
 
+! Check orthonormality
+
+! call matout(nVV,nVV,matmul(transpose(X1),X1) - matmul(transpose(Y1),Y1))
+! call matout(nOO,nOO,matmul(transpose(X2),X2) - matmul(transpose(Y2),Y2))
+! call matout(nVV,nOO,matmul(transpose(X1),X2) - matmul(transpose(Y1),Y2))
+! call matout(nOO,nVV,matmul(transpose(X2),X1) - matmul(transpose(Y2),Y1))
+
 end subroutine