diff --git a/src/RPA/sort_ppRPA.f90 b/src/RPA/sort_ppRPA.f90
index 4028c7f..f24743f 100644
--- a/src/RPA/sort_ppRPA.f90
+++ b/src/RPA/sort_ppRPA.f90
@@ -23,7 +23,9 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
   double precision,allocatable  :: S1(:,:)
   double precision,allocatable  :: S2(:,:)
   double precision,allocatable  :: O1(:,:)
-  double precision,allocatable  :: O2(:,:)
+  double precision,allocatable  :: O2(:,:)  
+  double precision,allocatable  :: tmp1(:,:)
+  double precision,allocatable  :: tmp2(:,:)
 
   integer,allocatable           :: order1(:)
   integer,allocatable           :: order2(:)
@@ -86,7 +88,7 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
 
     end if
 
-  end do
+ end do
 
   if(minval(Omega1(:)) < 0d0 .or. ab /= nVV) call print_warning('You may have instabilities in pp-RPA!!')
   if(maxval(Omega2(:)) > 0d0 .or. ij /= nOO) call print_warning('You may have instabilities in pp-RPA!!')
@@ -112,7 +114,7 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
     call quick_sort(Omega2(:),order2(:),nOO)
     call set_order(Z2(:,:),order2(:),nOO+nVV,nOO)
 
-  end if
+ end if
 
 ! Orthogonalize eigenvectors
 
@@ -196,16 +198,23 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
 !     ij_start = ij 
 
 !   end if
-! end do
-
+! end do 
+ 
   allocate(S1(nVV,nVV),S2(nOO,nOO),O1(nVV,nVV),O2(nOO,nOO))
-  S1 = + matmul(transpose(Z1),matmul(M,Z1))
+  allocate(tmp1(nOO+nVV,nVV),tmp2(nOO+nVV,nOO))
+
+  call dgemm ('N', 'N', nOO+nVV, nVV, nOO+nVV, 1d0,  M, nOO+nVV, Z1, nOO+nVV, 0d0, tmp1, nOO+nVV)
+  call dgemm ('T', 'N', nVV    , nVV, nOO+nVV, 1d0, Z1, nOO+nVV, tmp1, nOO+nVV, 0d0, S1, nVV)
+ 
   S2 = - matmul(transpose(Z2),matmul(M,Z2))
 
   if(nVV > 0) call orthogonalization_matrix(1,nVV,S1,O1)
   if(nOO > 0) call orthogonalization_matrix(1,nOO,S2,O2)
 
-  Z1 = matmul(Z1,O1)
+  
+  write (*,*) 'OK SO FAR'
+
+  call dgemm ('N', 'N', nVV, nVV, nVV, 1d0,  Z1, nVV, O1, nVV, 0d0, Z1, nVV)
   Z2 = matmul(Z2,O2)
 
 ! Define submatrices X1, Y1, X2, & Y2