diff --git a/src/GW/self_energy_correlation.f90 b/src/GW/self_energy_correlation.f90
index eae9a49..b965cdb 100644
--- a/src/GW/self_energy_correlation.f90
+++ b/src/GW/self_energy_correlation.f90
@@ -78,29 +78,43 @@ subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Ec
 
   ! Occupied part of the correlation self-energy
 
-    do jb=1,nS
-      do i=nC+1,nO
-        do q=nC+1,nBas-nR
-          do p=nC+1,nBas-nR
+!$OMP PARALLEL &
+!$OMP SHARED(SigC,rho,eta,nS,nC,nO,nBas,nR,e,Omega) &
+!$OMP PRIVATE(jb,i,q,p,eps) &
+!$OMP DEFAULT(NONE)
+!$OMP DO
+do q=nC+1,nBas-nR
+   do p=nC+1,nBas-nR
+      do jb=1,nS
+         do i=nC+1,nO
             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
-    end do
- 
+   end do
+end do
+!$OMP END DO
+!$OMP END PARALLEL
+
     ! Virtual part of the correlation self-energy
- 
-    do jb=1,nS
-      do a=nO+1,nBas-nR
-        do q=nC+1,nBas-nR
-          do p=nC+1,nBas-nR
+
+!$OMP PARALLEL &
+!$OMP SHARED(SigC,rho,eta,nS,nC,nO,nBas,nR,e,Omega) &
+!$OMP PRIVATE(jb,a,q,p,eps) &
+!$OMP DEFAULT(NONE)
+!$OMP DO  
+do q=nC+1,nBas-nR
+   do p=nC+1,nBas-nR
+      do jb=1,nS
+         do a=nO+1,nBas-nR
             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
-    end do
+   end do
+end do
+!$OMP END DO
+!$OMP END PARALLEL
 
     ! Galitskii-Migdal correlation energy
 
diff --git a/src/LR/linear_response.f90 b/src/LR/linear_response.f90
index e239f54..d966589 100644
--- a/src/LR/linear_response.f90
+++ b/src/LR/linear_response.f90
@@ -21,7 +21,9 @@ subroutine linear_response(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,E
   double precision,intent(in)   :: e(nBas)
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
 
-! Local variables
+  ! Local variables
+
+  integer :: i,j,k
 
   double precision              :: trace_matrix
   double precision,allocatable  :: A(:,:)
@@ -31,6 +33,7 @@ subroutine linear_response(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,E
   double precision,allocatable  :: AmBSq(:,:)
   double precision,allocatable  :: AmBIv(:,:)
   double precision,allocatable  :: Z(:,:)
+  double precision,allocatable  :: tmp(:,:)
 
 ! Output variables
 
@@ -41,7 +44,7 @@ subroutine linear_response(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,E
 
 ! Memory allocation
 
-  allocate(A(nS,nS),B(nS,nS),ApB(nS,nS),AmB(nS,nS),AmBSq(nS,nS),AmBIv(nS,nS),Z(nS,nS))
+  allocate(A(nS,nS),B(nS,nS),ApB(nS,nS),AmB(nS,nS),AmBSq(nS,nS),AmBIv(nS,nS),Z(nS,nS),tmp(nS,nS))
 
 ! Build A and B matrices 
 
@@ -77,10 +80,11 @@ subroutine linear_response(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,E
 !     if(Omega(ia) < 0d0) Omega(ia) = 0d0
 !   end do
 
-    call ADAt(nS,AmB,1d0*sqrt(Omega),AmBSq)
-    call ADAt(nS,AmB,1d0/sqrt(Omega),AmBIv)
+    call ADAt(nS,AmB,1d0*dsqrt(Omega),AmBSq)
+    call ADAt(nS,AmB,1d0/dsqrt(Omega),AmBIv)
 
-    Z = matmul(AmBSq,matmul(ApB,AmBSq))
+    call dgemm('N','N',nS,nS,nS,1d0,ApB,size(ApB,1),AmBSq,size(AmBSq,1),0d0,tmp,size(tmp,1))
+    call dgemm('N','N',nS,nS,nS,1d0,AmBSq,size(AmBSq,1),tmp,size(tmp,1),0d0,Z,size(Z,1))
 
    call diagonalize_matrix(nS,Z,Omega)
 
@@ -91,13 +95,13 @@ subroutine linear_response(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,E
   !   if(Omega(ia) < 0d0) Omega(ia) = 0d0
   ! end do
 
-    Omega = sqrt(Omega)
+    Omega = dsqrt(Omega)
 
-    XpY = matmul(transpose(Z),AmBSq)
-    call DA(nS,1d0/sqrt(Omega),XpY)
+    call dgemm('T','N',nS,nS,nS,1d0,Z,size(Z,1),AmBSq,size(AmBSq,1),0d0,XpY,size(XpY,1))
+    call DA(nS,1d0/dsqrt(Omega),XpY)
 
-    XmY = matmul(transpose(Z),AmBIv)
-    call DA(nS,1d0*sqrt(Omega),XmY)
+    call dgemm('T','N',nS,nS,nS,1d0,Z,size(Z,1),AmBIv,size(AmBIv,1),0d0,XmY,size(XmY,1))
+    call DA(nS,1d0*dsqrt(Omega),XmY)
  
   end if